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_r14273_HPC-02_Daley_Tiling/src/OCE – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/stpmlf.F90 @ 14787

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

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

Download in other formats: