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

Context Navigation

source: NEMO/trunk/src/OCE/stpmlf.F90 @ 14553

Last change on this file since 14553 was 14553, checked in by gsamson, 3 years ago
merge ticket2628_r14502_abl_restart_xios branch into trunk; sette identical between r14502 and r14544; ticket #2628
File size: 30.9 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 "do_loop_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, jtile ! dummy loop indice
96	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgdept
97	!! ---------------------------------------------------------------------
98	#if defined key_agrif
99	IF( nstop > 0 ) RETURN ! avoid to go further if an error was detected during previous time step (child grid)
100	kstp = nit000 + Agrif_Nb_Step()
101	Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
102	IF( lk_agrif_debug ) THEN
103	IF( Agrif_Root() .and. lwp) WRITE(,) '---'
104	IF(lwp) WRITE(,) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
105	ENDIF
106	IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.
107	# if defined key_xios
108	IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )
109	# endif
110	#endif
111	!
112	IF( ln_timing ) CALL timing_start('stp_MLF')
113	!
114	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
115	! model timestep
116	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
117	!
118	IF( l_1st_euler ) THEN ! start or restart with Euler 1st time-step
119	rDt = rn_Dt
120	r1_Dt = 1._wp / rDt
121	ENDIF
122	!
123	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124	! update I/O and calendar
125	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
126	!
127	IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
128	CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including possible AGRIF zoom)
129	IF( lk_diamlr ) CALL dia_mlr_iom_init ! with additional setup for multiple-linear-regression analysis
130	CALL iom_init_closedef
131	IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
132	ENDIF
133	IF( kstp == nitrst .AND. lwxios ) THEN
134	CALL iom_swap( cw_ocerst_cxt )
135	CALL iom_init_closedef( cw_ocerst_cxt )
136	CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
137	#if defined key_top
138	CALL iom_swap( cw_toprst_cxt )
139	CALL iom_init_closedef( cw_toprst_cxt )
140	CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
141	#endif
142	ENDIF
143	IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
144	#if defined key_si3
145	CALL iom_swap( cw_icerst_cxt )
146	CALL iom_init_closedef( cw_icerst_cxt )
147	CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
148	#endif
149	IF( ln_abl ) THEN
150	CALL iom_swap( cw_ablrst_cxt )
151	CALL iom_init_closedef( cw_ablrst_cxt )
152	CALL iom_setkt( kstp - nit000 + 1, cw_ablrst_cxt )
153	ENDIF
154	ENDIF
155	IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
156	CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
157	IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
158
159	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
160	! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
161	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
162	IF( ln_tide ) CALL tide_update( kstp ) ! update tide potential
163	IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
164	IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries
165	IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
166	CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice)
167
168	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
169	! Update stochastic parameters and random T/S fluctuations
170	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
171	IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
172	IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations
173
174	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
175	! Ocean physics update
176	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
177	! THERMODYNAMICS
178	CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points
179	CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points
180	CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
181	CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
182
183	! VERTICAL PHYSICS
184	CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
185
186	! LATERAL PHYSICS
187	!
188	IF( l_ldfslp ) THEN ! slope of lateral mixing
189	CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density
190
191	IF( ln_zps .AND. .NOT. ln_isfcav) &
192	& CALL zps_hde ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient
193	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
194
195	IF( ln_zps .AND. ln_isfcav) &
196	& CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
197	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
198	IF( ln_traldf_triad ) THEN
199	CALL ldf_slp_triad( kstp, Nbb, Nnn ) ! before slope for triad operator
200	ELSE
201	CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nnn ) ! before slope for standard operator
202	ENDIF
203	ENDIF
204	! ! eddy diffusivity coeff.
205	IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nnn ) ! and/or eiv coeff.
206	IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
207
208	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
209	! Ocean dynamics : hdiv, ssh, e3, u, v, w
210	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211
212	CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor)
213	IF( .NOT.lk_linssh ) THEN
214	CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) ) ! "after" ssh/h_0 ratio at t,u,v pts
215	IF( ln_dynspg_exp ) &
216	& 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
217	ENDIF
218	CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity
219	IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
220	ALLOCATE( zgdept(jpi,jpj,jpk) )
221	DO jk = 1, jpk
222	zgdept(:,:,jk) = gdept(:,:,jk,Nnn)
223	END DO
224	CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
225	DEALLOCATE( zgdept )
226
227	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
228	vv(:,:,:,Nrhs) = 0._wp
229
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 ) ! 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
255
256	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
257	! cool skin
258	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
259	IF ( ln_diurnal ) CALL diurnal_layers( kstp )
260
261	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
262	! diagnostics and outputs
263	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
264	IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats
265	IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics
266	CALL dia_hth ( kstp, Nnn ) ! Thermocline depth (20 degres isotherm depth)
267	IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports
268	CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag
269	CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics
270	CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
271	IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output
272	IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
273	IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis
274
275	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
276	! Now ssh filtering
277	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
278	CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height
279	IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh
280	#if defined key_top
281	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
282	! Passive Tracer Model
283	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
284	CALL trc_stp ( kstp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping
285	#endif
286
287	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
288	! Active tracers
289	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
290	! Loop over tile domains
291	DO jtile = 1, nijtile
292	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
293
294	DO_3D( 0, 0, 0, 0, 1, jpk )
295	ts(ji,jj,jk,:,Nrhs) = 0._wp ! set tracer trends to zero
296	END_3D
297
298	IF( lk_asminc .AND. ln_asmiau .AND. &
299	& ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment
300	CALL tra_sbc ( kstp, Nnn, ts, Nrhs ) ! surface boundary condition
301	IF( ln_traqsr ) CALL tra_qsr ( kstp, Nnn, ts, Nrhs ) ! penetrative solar radiation qsr
302	IF( ln_isf ) CALL tra_isf ( kstp, Nnn, ts, Nrhs ) ! ice shelf heat flux
303	IF( ln_trabbc ) CALL tra_bbc ( kstp, Nnn, ts, Nrhs ) ! bottom heat flux
304	IF( ln_trabbl ) CALL tra_bbl ( kstp, Nbb, Nnn, ts, Nrhs ) ! advective (and/or diffusive) bottom boundary layer scheme
305	IF( ln_tradmp ) CALL tra_dmp ( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends
306	IF( ln_bdy ) CALL bdy_tra_dmp( kstp, Nbb, ts, Nrhs ) ! bdy damping trends
307	END DO
308
309	#if defined key_agrif
310	IF(.NOT. Agrif_Root() ) THEN
311	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
312	CALL Agrif_Sponge_tra ! tracers sponge
313	ENDIF
314	#endif
315
316	! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
317	DO jtile = 1, nijtile
318	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
319
320	CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS
321	IF( ln_zdfmfc ) CALL tra_mfc ( kstp, Nbb, ts, Nrhs ) ! Mass Flux Convection
322	IF( ln_zdfosm ) THEN
323	CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS
324	IF( lrst_oce ) CALL osm_rst ( kstp, Nnn, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts
325	ENDIF
326	CALL tra_ldf ( kstp, Nbb, Nnn, ts, Nrhs ) ! lateral mixing
327
328	CALL tra_zdf ( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing and after tracer fields
329	IF( ln_zdfnpc ) CALL tra_npc ( kstp, Nnn, Nrhs, ts, Naa ) ! update after fields by non-penetrative convection
330	END DO
331
332	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Revert to tile over full domain
333	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
334	! Set boundary conditions, time filter and swap time levels
335	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
336	!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
337	!! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
338	!! If so:
339	!! (i) no need to call agrif update at initialization time
340	!! (ii) no need to update "before" fields
341	!!
342	!! Apart from creating new tra_swp/dyn_swp routines, this however:
343	!! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
344	!! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
345	!! e.g. a shift of the feedback interface inside child domain.
346	!! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
347	!! place.
348	!!
349	IF( ln_dynspg_ts ) CALL mlf_baro_corr ( Nnn, Naa, uu, vv ) ! barotrope adjustment
350	CALL finalize_lbc ( kstp, Nbb , Naa, uu, vv, ts ) ! boundary conditions
351	CALL tra_atf_qco ( kstp, Nbb, Nnn, Naa , ts ) ! time filtering of "now" tracer arrays
352	CALL dyn_atf_qco ( kstp, Nbb, Nnn, Naa, uu, vv ) ! time filtering of "now" velocities
353	IF(.NOT.lk_linssh) THEN
354	r3t(:,:,Nnn) = r3t_f(:,:) ! update now ssh/h_0 with time filtered values
355	r3u(:,:,Nnn) = r3u_f(:,:)
356	r3v(:,:,Nnn) = r3v_f(:,:)
357	ENDIF
358	!
359	! Swap time levels
360	Nrhs = Nbb
361	Nbb = Nnn
362	Nnn = Naa
363	Naa = Nrhs
364	!
365	!
366	IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics
367
368	!!gm : This does not only concern the dynamics ==>>> add a new title
369	!!gm2: why ouput restart before AGRIF update?
370	!!
371	!!jc: That would be better, but see comment above
372	!!
373	IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file
374	IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
375
376	#if defined key_agrif
377	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
378	! AGRIF recursive integration
379	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
380	Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
381	CALL Agrif_Integrate_ChildGrids( stp_MLF ) ! allows to finish all the Child Grids before updating
382
383	#endif
384	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
385	! Control
386	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
387	CALL stp_ctl ( kstp, Nnn )
388
389	#if defined key_agrif
390	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
391	! AGRIF update
392	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
393	IF( Agrif_NbStepint() == 0 .AND. nstop == 0 ) &
394	& CALL Agrif_update_all( ) ! Update all components
395
396	#endif
397	IF( ln_diaobs .AND. nstop == 0 ) &
398	& CALL dia_obs( kstp, Nnn ) ! obs-minus-model (assimilation) diags (after dynamics update)
399
400	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
401	! File manipulation at the end of the first time step
402	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
403	IF( kstp == nit000 ) THEN ! 1st time step only
404	CALL iom_close( numror ) ! close input ocean restart file
405	IF( lrxios ) CALL iom_context_finalize( cr_ocerst_cxt )
406	IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
407	IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
408	ENDIF
409
410	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
411	! Coupled mode
412	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
413	IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
414	!
415	#if defined key_xios
416	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
417	! Finalize contextes if end of simulation or error detected
418	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
419	IF( kstp == nitend .OR. nstop > 0 ) THEN
420	CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
421	IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
422	ENDIF
423	#endif
424	!
425	IF( l_1st_euler ) THEN ! recover Leap-frog timestep
426	rDt = 2._wp * rn_Dt
427	r1_Dt = 1._wp / rDt
428	l_1st_euler = .FALSE.
429	ENDIF
430	!
431	IF( ln_timing ) CALL timing_stop('stp_MLF')
432	!
433	END SUBROUTINE stp_MLF
434
435
436	SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv )
437	!!----------------------------------------------------------------------
438	!! * ROUTINE mlf_baro_corr *
439	!!
440	!! ** Purpose : Finalize after horizontal velocity.
441	!!
442	!! ** Method : * Ensure after velocities transport matches time splitting
443	!! estimate (ln_dynspg_ts=T)
444	!!
445	!! ** Action : puu(Kmm),pvv(Kmm) updated now horizontal velocity (ln_bt_fw=F)
446	!! puu(Kaa),pvv(Kaa) after horizontal velocity
447	!!----------------------------------------------------------------------
448	USE dynspg_ts, ONLY : un_adv, vn_adv ! updated Kmm barotropic transport
449	!!
450	INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices
451	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities
452	!
453	INTEGER :: jk ! dummy loop indices
454	REAL(wp), DIMENSION(jpi,jpj) :: zue, zve
455	!!----------------------------------------------------------------------
456
457	! Ensure below that barotropic velocities match time splitting estimate
458	! Compute actual transport and replace it with ts estimate at "after" time step
459	zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
460	zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
461	DO jk = 2, jpkm1
462	zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
463	zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
464	END DO
465	DO jk = 1, jpkm1
466	puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
467	pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
468	END DO
469	!
470	IF( .NOT.ln_bt_fw ) THEN
471	! Remove advective velocity from "now velocities"
472	! prior to asselin filtering
473	! In the forward case, this is done below after asselin filtering
474	! so that asselin contribution is removed at the same time
475	DO jk = 1, jpkm1
476	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )umask(:,:,jk)
477	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )vmask(:,:,jk)
478	END DO
479	ENDIF
480	!
481	END SUBROUTINE mlf_baro_corr
482
483
484	SUBROUTINE finalize_lbc( kt, Kbb, Kaa, puu, pvv, pts )
485	!!----------------------------------------------------------------------
486	!! * ROUTINE finalize_lbc *
487	!!
488	!! ** Purpose : Apply the boundary condition on the after velocity
489	!!
490	!! ** Method : * Apply lateral boundary conditions on after velocity
491	!! at the local domain boundaries through lbc_lnk call,
492	!! at the one-way open boundaries (ln_bdy=T),
493	!! at the AGRIF zoom boundaries (lk_agrif=T)
494	!!
495	!! ** Action : puu(Kaa),pvv(Kaa) after horizontal velocity and tracers
496	!!----------------------------------------------------------------------
497	#if defined key_agrif
498	USE agrif_oce_interp
499	#endif
500	USE bdydyn ! ocean open boundary conditions (define bdy_dyn)
501	!!
502	INTEGER , INTENT(in ) :: kt ! ocean time-step index
503	INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices
504	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt) , INTENT(inout) :: puu, pvv ! velocities to be time filtered
505	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers
506	!!----------------------------------------------------------------------
507	!
508	! Update after tracer and velocity on domain lateral boundaries
509	!
510	# if defined key_agrif
511	CALL Agrif_tra !* AGRIF zoom boundaries
512	CALL Agrif_dyn( kt )
513	# endif
514	! ! local domain boundaries (T-point, unchanged sign)
515	CALL lbc_lnk( 'finalize_lbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. &
516	& , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
517	!
518	! !* BDY open boundaries
519	IF( ln_bdy ) THEN
520	CALL bdy_tra( kt, Kbb, pts, Kaa )
521	IF( ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
522	IF( ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
523	ENDIF
524	!
525	END SUBROUTINE finalize_lbc
526
527	#else
528	!!----------------------------------------------------------------------
529	!! default option EMPTY MODULE qco not activated
530	!!----------------------------------------------------------------------
531	#endif
532
533	!!======================================================================
534	END MODULE stpmlf

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

Download in other formats: