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.
step.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 24.8 KB
Line 
1MODULE step
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   !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs
27   !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state
28   !!----------------------------------------------------------------------
29
30   !!----------------------------------------------------------------------
31   !!   stp             : OPA system time-stepping
32   !!----------------------------------------------------------------------
33   USE step_oce         ! time stepping definition modules
34   USE iom
35   USE lbclnk
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   stp   ! called by opa.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44!!gm   #  include "zdfddm_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52#if defined key_agrif
53   RECURSIVE SUBROUTINE stp( )
54      INTEGER             ::   kstp   ! ocean time-step index
55#else
56   SUBROUTINE stp( kstp )
57      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
58#endif
59      !!----------------------------------------------------------------------
60      !!                     ***  ROUTINE stp  ***
61      !!
62      !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
63      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
64      !!              - Tme stepping  of TRC (passive tracer eqs.)
65      !!
66      !! ** Method  : -1- Update forcings and data
67      !!              -2- Update ocean physics
68      !!              -3- Compute the t and s trends
69      !!              -4- Update t and s
70      !!              -5- Compute the momentum trends
71      !!              -6- Update the horizontal velocity
72      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
73      !!              -8- Outputs and diagnostics
74      !!----------------------------------------------------------------------
75      INTEGER ::   jk       ! dummy loop indice
76      INTEGER ::   tind     ! tracer loop index
77      INTEGER ::   indic    ! error indicator if < 0
78      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
79      !! ---------------------------------------------------------------------
80
81#if defined key_agrif
82      kstp = nit000 + Agrif_Nb_Step()
83      IF ( lk_agrif_debug ) THEN
84         IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
85         IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp, 'int tstep',Agrif_NbStepint()
86      ENDIF
87
88      IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
89
90# if defined key_iomput
91      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
92# endif
93#endif
94                             indic = 0           ! reset to no error condition
95      IF( kstp == nit000 ) THEN
96         ! must be done after nemo_init for AGRIF+XIOS+OASIS
97                      CALL iom_init(      cxios_context          )  ! iom_put initialization
98         IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" )  ! initialize context for coarse grid
99      ENDIF
100
101      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
102                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell iom we are at time step kstp
103      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp
104
105      IF( ln_bias )          CALL bias_opn( kstp )
106
107      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
108      ! Update data, open boundaries, surface boundary condition (including sea-ice)
109      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
110      IF( lk_tide    )   CALL sbc_tide( kstp )
111      IF( lk_bdy     )  THEN
112         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib
113                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
114      ENDIF
115
116      ! We must ensure that tsb halos are up to date on EVERY timestep.
117      DO tind = 1, jpts
118         CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. )
119      END DO
120
121                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
122                                                      ! clem: moved here for bdy ice purpose
123      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124      ! Update stochastic parameters and random T/S fluctuations
125      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
126       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
127       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
128
129      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
130      ! Ocean physics update                (ua, va, tsa used as workspace)
131      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
132      !  THERMODYNAMICS
133                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
134                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
135                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
136                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
137      !
138      !  VERTICAL PHYSICS
139                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
140      !                                               ! Vertical eddy viscosity and diffusivity coefficients
141      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
142      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
143      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
144      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
145      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
146         avt (:,:,:) = rn_avt0 * wmask (:,:,:)
147         avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
148         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
149      ENDIF
150      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
151         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
152      ENDIF
153      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
154
155      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
156
157      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
158         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
159
160                         CALL zdf_mxl( kstp )         ! mixed layer depth
161
162                                                      ! write TKE or GLS information in the restart file
163      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
164      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
165      !
166      !  LATERAL  PHYSICS
167      !
168      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
169                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
170         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
171            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
172            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
173         IF( ln_zps .AND.       ln_isfcav)                               &
174            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
175            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
176            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
177         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
178                         CALL ldf_slp_grif( kstp )
179         ELSE
180                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
181         ENDIF
182      ENDIF
183#if defined key_traldf_c2d
184      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
185#endif
186#if defined key_traldf_c3d && defined key_traldf_smag
187                          CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
188#  endif
189#if defined key_dynldf_c3d && defined key_dynldf_smag
190                          CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
191#  endif
192
193      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
194      !  Ocean dynamics : hdiv, rot, ssh, e3, wn
195      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
196                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
197      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
198                         CALL wzv           ( kstp )  ! now cross-level velocity
199
200      IF( lk_dynspg_ts ) THEN 
201          ! In case the time splitting case, update almost all momentum trends here:
202          ! Note that the computation of vertical velocity above, hence "after" sea level
203          ! is necessary to compute momentum advection for the rhs of barotropic loop:
204                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
205            IF( ln_zps .AND. .NOT. ln_isfcav)                               &
206               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
207               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
208            IF( ln_zps .AND.       ln_isfcav)                               &
209               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
210               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
211               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
212
213                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero
214                                  va(:,:,:) = 0.e0
215          IF(  lk_asminc .AND. ln_asmiau .AND. &
216             & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment
217          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified)
218          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends
219                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form)
220                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis
221                                  CALL dyn_ldf      ( kstp )   ! lateral mixing
222          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
223#if defined key_agrif
224          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
225#endif
226                                  CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
227                                  CALL dyn_spg( kstp, indic )  ! surface pressure gradient
228
229                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated at this stage)
230                                  va_sv(:,:,:) = va(:,:,:)
231
232                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
233          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
234                                  CALL wzv           ( kstp )  ! now cross-level velocity
235      ENDIF
236
237      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
238      ! diagnostics and outputs             (ua, va, tsa used as workspace)
239      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240      IF( lk_floats  )      CALL flo_stp( kstp )         ! drifting Floats
241      IF( lk_diahth  )      CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
242      IF( .NOT. ln_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
243      IF( lk_diadct  )      CALL dia_dct( kstp )         ! Transports
244      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
245      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis
246                            CALL dia_prod( kstp )        ! ocean model: product diagnostics
247                            CALL dia_wri( kstp )         ! ocean model: outputs
248      !
249      IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
250
251#if defined key_top
252      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
253      ! Passive Tracer Model
254      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255      IF( lk_asminc .AND. ln_asmiau .AND. lk_bgcinc ) &
256         &               CALL bgc_asm_inc( kstp )     ! biogeochemistry assimilation
257                         CALL trc_stp( kstp )         ! time-stepping
258#endif
259
260
261      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
262      ! Active tracers                              (ua, va used as workspace)
263      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
264                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
265
266      IF(  lk_asminc .AND. ln_asmiau .AND. &
267         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
268                             CALL tra_sbc    ( kstp )       ! surface boundary condition
269      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
270      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
271      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
272      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
273      IF( ln_bias        )   CALL tra_bias   ( kstp )
274      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
275                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
276      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
277                             CALL tra_ldf    ( kstp )       ! lateral mixing
278
279      IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics
280
281#if defined key_agrif
282      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
283#endif
284                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
285
286      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
287         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
288                             CALL tra_nxt( kstp )                ! tracer fields at next time step
289                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
290            IF( ln_zps .AND. .NOT. ln_isfcav)                                &
291               &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
292               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
293            IF( ln_zps .AND.       ln_isfcav)                                &
294               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
295               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
296               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
297            IF( ln_bias )    CALL dyn_bias( kstp )
298      ELSE                                                  ! centered hpg  (eos then time stepping)
299         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
300                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
301         IF( ln_zps .AND. .NOT. ln_isfcav)                                   &
302               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
303               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
304         IF( ln_zps .AND.       ln_isfcav)                                   & 
305               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
306               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
307               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
308         ENDIF
309         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
310                             CALL tra_nxt( kstp )                ! tracer fields at next time step
311         IF( ln_bias )       CALL dyn_bias( kstp )
312      ENDIF
313
314      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
315      ! Dynamics                                    (tsa used as workspace)
316      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
317      IF( lk_dynspg_ts   )  THEN
318                                                             ! revert to previously computed momentum tendencies
319                                                             ! (not using ua, va as temporary arrays during tracers' update could avoid that)
320                               ua(:,:,:) = ua_sv(:,:,:)
321                               va(:,:,:) = va_sv(:,:,:)
322                                                             ! Revert now divergence and rotational to previously computed ones
323                                                             !(needed because of the time swap in div_cur, at the beginning of each time step)
324                               hdivn(:,:,:) = hdivb(:,:,:)
325                               rotn(:,:,:)  = rotb(:,:,:) 
326
327                               CALL dyn_bfr( kstp )         ! bottom friction
328                               CALL dyn_zdf( kstp )         ! vertical diffusion
329      ELSE
330                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero
331                               va(:,:,:) = 0.e0
332
333        IF(  lk_asminc .AND. ln_asmiau .AND. &
334           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
335        IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
336        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
337        IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
338                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
339                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
340                               CALL dyn_ldf( kstp )         ! lateral mixing
341        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
342#if defined key_agrif
343        IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
344#endif
345                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
346                               CALL dyn_bfr( kstp )         ! bottom friction
347                               CALL dyn_zdf( kstp )         ! vertical diffusion
348                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
349      ENDIF
350                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
351
352                               CALL ssh_swp( kstp )         ! swap of sea surface height
353      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
354      !
355      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
356      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
357
358#if defined key_agrif
359      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
360      ! AGRIF
361      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
362                               CALL Agrif_Integrate_ChildGrids( stp ) 
363
364      IF ( Agrif_NbStepint().EQ.0 ) THEN
365                               CALL Agrif_Update_Tra()      ! Update active tracers
366                               CALL Agrif_Update_Dyn()      ! Update momentum
367      ENDIF
368#endif
369      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
370      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
371
372      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
373      ! Control
374      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
375                               CALL stp_ctl( kstp, indic )
376      IF( indic < 0        )   THEN
377                               CALL ctl_stop( 'step: indic < 0' )
378                               CALL dia_wri_state( 'output.abort', kstp )
379      ENDIF
380      IF( kstp == nit000   )   THEN
381                 CALL iom_close( numror )     ! close input  ocean restart file
382         IF(lwm .AND. numond /= -1) CALL FLUSH    ( numond )     ! flush output namelist oce
383         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
384      ENDIF
385
386
387      IF( lrst_bias )          CALL bias_wrt     ( kstp )
388
389      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
390      ! Coupled mode
391      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
392      !IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
393      !
394#if defined key_iomput
395      IF( kstp == nitend .OR. indic < 0 ) THEN
396                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
397         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
398      ENDIF
399#endif
400      !
401      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
402      !     
403      !
404   END SUBROUTINE stp
405
406   !!======================================================================
407END MODULE step
Note: See TracBrowser for help on using the repository browser.