Context Navigation

step.F90 @ 11500

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: