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/2021/dev_r14318_RK3_stage1/src/OCE – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90 @ 15193

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

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

Download in other formats: