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.

stpmlf.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/stpmlf.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago
Merge trunk -r14642:HEAD
File size: 31.0 KB

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

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

Download in other formats: