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

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/stpMLF.F90 @ 13416

Last change on this file since 13416 was 13416, checked in by gm, 4 years ago
First dev of BVP, cond 3.4, gridF fixes
File size: 27.4 KB

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

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

Download in other formats: