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.

step.F90 in NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/step.F90 @ 11355

Last change on this file since 11355 was 11355, checked in by dancopsey, 5 years ago
Add lots of print statements printing out values of files throughout the ocean_ice timestep.. New output will be in files crash_stats_419.out and crash_stats_68.out.
File size: 25.3 KB

Line
1	MODULE step
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	!!----------------------------------------------------------------------
34
35	!!----------------------------------------------------------------------
36	!! stp : OPA system time-stepping
37	!!----------------------------------------------------------------------
38	USE step_oce ! time stepping definition modules
39	!
40	USE iom ! xIOs server
41
42	IMPLICIT NONE
43	PRIVATE
44
45	PUBLIC stp ! called by nemogcm.F90
46
47	!!----------------------------------------------------------------------
48	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
49	!! $Id$
50	!! Software governed by the CeCILL license (see ./LICENSE)
51	!!----------------------------------------------------------------------
52	CONTAINS
53
54	#if defined key_agrif
55	RECURSIVE SUBROUTINE stp( )
56	INTEGER :: kstp ! ocean time-step index
57	#else
58	SUBROUTINE stp( kstp )
59	INTEGER, INTENT(in) :: kstp ! ocean time-step index
60	#endif
61	!!----------------------------------------------------------------------
62	!! * ROUTINE stp *
63	!!
64	!! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
65	!! - Time stepping of SI3 (dynamic and thermodynamic eqs.)
66	!! - Time stepping of TRC (passive tracer eqs.)
67	!!
68	!! ** Method : -1- Update forcings and data
69	!! -2- Update ocean physics
70	!! -3- Compute the t and s trends
71	!! -4- Update t and s
72	!! -5- Compute the momentum trends
73	!! -6- Update the horizontal velocity
74	!! -7- Compute the diagnostics variables (rd,N2, hdiv,w)
75	!! -8- Outputs and diagnostics
76	!!----------------------------------------------------------------------
77	INTEGER :: ji, jj, jk ! dummy loop indice
78	INTEGER :: indic ! error indicator if < 0
79	!!gm kcall can be removed, I guess
80	INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt)
81	!! ---------------------------------------------------------------------
82	#if defined key_agrif
83	kstp = nit000 + Agrif_Nb_Step()
84	IF( lk_agrif_debug ) THEN
85	IF( Agrif_Root() .and. lwp) WRITE(,) '---'
86	IF(lwp) WRITE(,) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
87	ENDIF
88	IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.
89	# if defined key_iomput
90	IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )
91	# endif
92	#endif
93	!
94	IF( ln_timing ) CALL timing_start('stp')
95
96	IF(narea == 419) THEN
97	WRITE(9419,*) 'max sshn step before update IO = ',MAXVAL( ABS( sshn(:,:) ) )
98	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
99	WRITE(9419,*) 'max emp_b = ',MAXVAL( ABS( emp_b(:,:) ) )
100	WRITE(9419,*) 'max emp = ',MAXVAL( ABS( emp(:,:) ) )
101	ENDIF
102
103	!
104	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
105	! update I/O and calendar
106	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
107	indic = 0 ! reset to no error condition
108
109	IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
110	CALL iom_init( cxios_context ) ! for model grid (including passible AGRIF zoom)
111	IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
112	ENDIF
113	IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
114	CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
115	IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
116
117	IF(narea == 419) THEN
118	WRITE(9419,*) 'max sshn step before external forcing = ',MAXVAL( ABS( sshn(:,:) ) )
119	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
120	WRITE(9419,*) 'max emp_b = ',MAXVAL( ABS( emp_b(:,:) ) )
121	WRITE(9419,*) 'max emp = ',MAXVAL( ABS( emp(:,:) ) )
122	ENDIF
123
124	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
125	! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
126	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
127	IF( ln_tide ) CALL sbc_tide( kstp ) ! update tide potential
128	IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
129	IF( ln_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries
130	CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice)
131
132
133
134
135	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
136	! Update stochastic parameters and random T/S fluctuations
137	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
138	IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
139	IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations
140
141	IF(narea == 419) THEN
142	WRITE(9419,*) 'max sshn step before ocean physics = ',MAXVAL( ABS( sshn(:,:) ) )
143	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
144	WRITE(9419,*) 'max emp_b = ',MAXVAL( ABS( emp_b(:,:) ) )
145	WRITE(9419,*) 'max emp = ',MAXVAL( ABS( emp(:,:) ) )
146	ENDIF
147
148	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
149	! Ocean physics update
150	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
151	! THERMODYNAMICS
152	CALL eos_rab( tsb, rab_b ) ! before local thermal/haline expension ratio at T-points
153	CALL eos_rab( tsn, rab_n ) ! now local thermal/haline expension ratio at T-points
154	CALL bn2 ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
155	CALL bn2 ( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency
156
157	! VERTICAL PHYSICS
158	CALL zdf_phy( kstp ) ! 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( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density
164
165	IF( ln_zps .AND. .NOT. ln_isfcav) &
166	& CALL zps_hde ( kstp, jpts, tsb, 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, jpts, tsb, 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 ) ! before slope for triad operator
174	ELSE
175	CALL ldf_slp ( kstp, rhd, rn2b ) ! 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 ) ! and/or eiv coeff.
180	IF( l_ldfdyn_time ) CALL ldf_dyn( kstp ) ! eddy viscosity coeff.
181
182	IF(narea == 419) THEN
183	WRITE(9419,*) 'max sshn step before ocean dynamics = ',MAXVAL( ABS( sshn(:,:) ) )
184	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
185	WRITE(9419,*) 'max emp_b = ',MAXVAL( ABS( emp_b(:,:) ) )
186	WRITE(9419,*) 'max emp = ',MAXVAL( ABS( emp(:,:) ) )
187	ENDIF
188
189	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
190	! Ocean dynamics : hdiv, ssh, e3, u, v, w
191	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
192
193	CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor)
194	IF(narea == 419) THEN
195	WRITE(9419,*) 'max ssha step before dom_vvl_sf_nxt = ',MAXVAL( ABS( ssha(:,:) ) )
196	CALL FLUSH(9419)
197	ENDIF
198
199	IF(narea == 68) THEN
200	WRITE(968,*) 'max ssha step before dom_vvl_sf_nxt = ',MAXVAL( ABS( ssha(:,:) ) )
201	CALL FLUSH(968)
202	ENDIF
203
204
205	IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors
206
207	IF(narea == 419) THEN
208	WRITE(9419,*) 'max ssha step before wzv = ',MAXVAL( ABS( ssha(:,:) ) )
209	ENDIF
210
211	CALL wzv ( kstp ) ! now cross-level velocity
212
213	IF(narea == 419) THEN
214	WRITE(9419,*) 'max ssha step before wAimp = ',MAXVAL( ABS( ssha(:,:) ) )
215	ENDIF
216
217	IF( ln_zad_Aimp ) CALL wAimp ( kstp ) ! Adaptive-implicit vertical advection partitioning
218	CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation
219
220	IF(narea == 419) THEN
221	WRITE(9419,*) 'max ssha step before simplification = ',MAXVAL( ABS( ssha(:,:) ) )
222	ENDIF
223
224
225	!!jc: fs simplification
226	!!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)
227	!! but ensures reproductible results
228	!! with previous versions using split-explicit free surface
229	IF( ln_zps .AND. .NOT. ln_isfcav ) &
230	& CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient
231	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
232	IF( ln_zps .AND. ln_isfcav ) &
233	& CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
234	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
235	!!jc: fs simplification
236
237	ua(:,:,:) = 0._wp ! set dynamics trends to zero
238	va(:,:,:) = 0._wp
239
240	IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &
241	& CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment
242
243	IF(narea == 419) THEN
244	WRITE(9419,*) 'max ssha step before bdy_dyn3d_dmp = ',MAXVAL( ABS( ssha(:,:) ) )
245	ENDIF
246
247	IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends
248	#if defined key_agrif
249	IF(.NOT. Agrif_Root()) &
250	& CALL Agrif_Sponge_dyn ! momentum sponge
251	#endif
252
253	IF(narea == 419) THEN
254	WRITE(9419,*) 'max ssha step before dom_adv = ',MAXVAL( ABS( ssha(:,:) ) )
255	ENDIF
256
257	CALL dyn_adv ( kstp ) ! advection (vector or flux form)
258
259	IF(narea == 419) THEN
260	WRITE(9419,*) 'max ssha step before dyn_vor = ',MAXVAL( ABS( ssha(:,:) ) )
261	ENDIF
262
263	CALL dyn_vor ( kstp ) ! vorticity term including Coriolis
264
265	IF(narea == 419) THEN
266	WRITE(9419,*) 'max ssha step before dyn_ldft = ',MAXVAL( ABS( ssha(:,:) ) )
267	ENDIF
268
269	CALL dyn_ldf ( kstp ) ! lateral mixing
270
271	IF(narea == 419) THEN
272	WRITE(9419,*) 'max ssha step before dyn_osm = ',MAXVAL( ABS( ssha(:,:) ) )
273	ENDIF
274
275	IF( ln_zdfosm ) CALL dyn_osm ( kstp ) ! OSMOSIS non-local velocity fluxes
276
277	IF(narea == 419) THEN
278	WRITE(9419,*) 'max ssha step before dyn_hpg = ',MAXVAL( ABS( ssha(:,:) ) )
279	ENDIF
280
281	CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure
282
283	IF(narea == 419) THEN
284	WRITE(9419,*) 'max ssha step before dyn_spg = ',MAXVAL( ABS( ssha(:,:) ) )
285	ENDIF
286
287	CALL dyn_spg ( kstp ) ! surface pressure gradient
288
289
290	IF(narea == 419) THEN
291	WRITE(9419,*) 'max ssha step before div_hor = ',MAXVAL( ABS( ssha(:,:) ) )
292	ENDIF
293
294	! With split-explicit free surface, since now transports have been updated and ssha as well
295	IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
296	CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case)
297	IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component)
298	CALL wzv ( kstp ) ! now cross-level velocity
299	IF( ln_zad_Aimp ) CALL wAimp ( kstp ) ! Adaptive-implicit vertical advection partitioning
300	ENDIF
301
302	IF(narea == 419) THEN
303	WRITE(9419,*) 'max ssha step before dyn_zdf = ',MAXVAL( ABS( ssha(:,:) ) )
304	ENDIF
305
306
307	CALL dyn_zdf ( kstp ) ! vertical diffusion
308
309	IF(narea == 419) THEN
310	WRITE(9419,*) 'max sshn step before cool skin = ',MAXVAL( ABS( sshn(:,:) ) )
311	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
312	ENDIF
313
314	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
315	! cool skin
316	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
317	IF ( ln_diurnal ) CALL stp_diurnal( kstp )
318
319	IF(narea == 419) THEN
320	WRITE(9419,*) 'max sshn step before diagnostics = ',MAXVAL( ABS( sshn(:,:) ) )
321	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
322	ENDIF
323
324	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
325	! diagnostics and outputs
326	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
327	IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats
328	IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics
329	IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth)
330	IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports
331	CALL dia_ar5 ( kstp ) ! ar5 diag
332	IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis
333	CALL dia_prod( kstp ) ! ocean model: product diagnostics
334	CALL dia_wri ( kstp ) ! ocean model: outputs
335	!
336	IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output
337
338	#if defined key_top
339	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
340	! Passive Tracer Model
341	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
342	CALL trc_stp ( kstp ) ! time-stepping
343	#endif
344
345	IF(narea == 419) THEN
346	WRITE(9419,*) 'max sshn step before active tracers = ',MAXVAL( ABS( sshn(:,:) ) )
347	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
348	ENDIF
349
350	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
351	! Active tracers
352	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
353	tsa(:,:,:,:) = 0._wp ! set tracer trends to zero
354
355	IF( lk_asminc .AND. ln_asmiau .AND. &
356	& ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment
357
358	IF(narea == 419) THEN
359	WRITE(9419,*) 'max sshn step before tra_sbc = ',MAXVAL( ABS( sshn(:,:) ) )
360	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
361	ENDIF
362
363	CALL tra_sbc ( kstp ) ! surface boundary condition
364
365	IF(narea == 419) THEN
366	WRITE(9419,*) 'max sshn step after tra_sbc = ',MAXVAL( ABS( sshn(:,:) ) )
367	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
368	ENDIF
369
370	IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr
371	IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux
372	IF( ln_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme
373	IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends
374	IF( ln_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends
375	#if defined key_agrif
376	IF(.NOT. Agrif_Root()) &
377	& CALL Agrif_Sponge_tra ! tracers sponge
378	#endif
379	CALL tra_adv ( kstp ) ! horizontal & vertical advection
380	IF( ln_zdfosm ) CALL tra_osm ( kstp ) ! OSMOSIS non-local tracer fluxes
381	IF( lrst_oce .AND. ln_zdfosm ) &
382	& CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
383	CALL tra_ldf ( kstp ) ! lateral mixing
384
385	!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
386	IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics
387	!!gm
388	CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields
389	IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection
390
391	IF(narea == 419) THEN
392	WRITE(9419,*) 'max sshn step before boundary conditions = ',MAXVAL( ABS( sshn(:,:) ) )
393	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
394	ENDIF
395
396	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
397	! Set boundary conditions and Swap
398	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
399	!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
400	!! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
401	!! If so:
402	!! (i) no need to call agrif update at initialization time
403	!! (ii) no need to update "before" fields
404	!!
405	!! Apart from creating new tra_swp/dyn_swp routines, this however:
406	!! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
407	!! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
408	!! e.g. a shift of the feedback interface inside child domain.
409	!! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
410	!! place.
411	!!
412	!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
413	CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap
414
415	IF(narea == 419) THEN
416	WRITE(9419,*) 'max sshn step before dyn_nxt = ',MAXVAL( ABS( sshn(:,:) ) )
417	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
418	ENDIF
419
420	CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
421
422	IF(narea == 419) THEN
423	WRITE(9419,*) 'max sshn step before ssh_swp = ',MAXVAL( ABS( sshn(:,:) ) )
424	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
425	ENDIF
426
427	CALL ssh_swp ( kstp ) ! swap of sea surface height
428
429	IF(narea == 419) THEN
430	WRITE(9419,*) 'max sshn step before dom_vvl_sf_swp = ',MAXVAL( ABS( sshn(:,:) ) )
431	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
432	ENDIF
433
434	IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors
435	!
436
437	IF(narea == 419) THEN
438	WRITE(9419,*) 'max sshn step before dia_hsb = ',MAXVAL( ABS( sshn(:,:) ) )
439	WRITE(9419,*) 'max ssha = ',MAXVAL( ABS( ssha(:,:) ) )
440	ENDIF
441
442	IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics
443
444	!!gm : This does not only concern the dynamics ==>>> add a new title
445	!!gm2: why ouput restart before AGRIF update?
446	!!
447	!!jc: That would be better, but see comment above
448	!!
449	IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file
450	IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
451
452	IF(narea == 419) THEN
453	WRITE(9419,*) 'max sshn step before agrif = ',MAXVAL( ABS( sshn(:,:) ) )
454	ENDIF
455
456	#if defined key_agrif
457	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
458	! AGRIF
459	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
460	CALL Agrif_Integrate_ChildGrids( stp ) ! allows to finish all the Child Grids before updating
461
462	IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components
463	#endif
464	IF( ln_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
465
466	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
467	! Control
468	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
469
470	IF(narea == 419) THEN
471	WRITE(9419,*) 'max sshn step before stp_ctl = ',MAXVAL( ABS( sshn(:,:) ) )
472	ENDIF
473
474	CALL stp_ctl ( kstp, indic )
475
476	IF( kstp == nit000 ) THEN ! 1st time step only
477	CALL iom_close( numror ) ! close input ocean restart file
478	IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
479	IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
480	ENDIF
481
482	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
483	! Coupled mode
484	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
485	!!gm why lk_oasis and not lk_cpl ????
486	IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges
487	!
488	#if defined key_iomput
489	IF( kstp == nitend .OR. indic < 0 ) THEN
490	CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
491	IF(lrxios) CALL iom_context_finalize( crxios_context )
492	IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
493	ENDIF
494	#endif
495	!
496	IF( ln_timing ) CALL timing_stop('stp')
497	!
498	END SUBROUTINE stp
499
500	!!======================================================================
501	END MODULE step

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

Download in other formats: