source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/step.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 21 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File size: 27.3 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      IF( ln_stopack )   CALL stopack_pert( kstp )
122                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
123                                                      ! clem: moved here for bdy ice purpose
124
125      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
126      ! Update stochastic parameters and random T/S fluctuations
127      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
128
129      IF( ln_sto_eos )                 CALL sto_par( kstp )   ! Stochastic parameters
130      IF( ln_sto_eos )                 CALL sto_pts( tsn  )   ! Random T/S fluctuations
131      IF( ln_stopack .AND. ln_skeb  )  CALL skeb_comp( kstp ) ! Stochastic Energy Backscatter
132
133      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
134      ! Ocean physics update                (ua, va, tsa used as workspace)
135      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
136      !  THERMODYNAMICS
137                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
138                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
139                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
140                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
141      !
142      !  VERTICAL PHYSICS
143                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
144      !                                               ! Vertical eddy viscosity and diffusivity coefficients
145      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
146      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
147      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
148      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
149      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
150         avt (:,:,:) = rn_avt0 * wmask (:,:,:)
151         avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
152         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
153      ENDIF
154      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
155         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) THEN
156              rn_avt_rnf0 = rn_avt_rnf
157              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf )
158         ENDIF
159         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
160      ENDIF
161      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
162
163      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
164
165      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
166         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
167
168                         CALL zdf_mxl( kstp )         ! mixed layer depth
169
170                                                      ! write TKE or GLS information in the restart file
171      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
172      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
173      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' )
174      !
175      !  LATERAL  PHYSICS
176      !
177      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
178                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
179         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
180            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
181            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
182         IF( ln_zps .AND.       ln_isfcav)                               &
183            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
184            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
185            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
186         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
187                         CALL ldf_slp_grif( kstp )
188         ELSE
189                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
190         ENDIF
191      ENDIF
192#if defined key_traldf_c2d
193      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
194#endif
195#if defined key_traldf_c3d && defined key_traldf_smag
196                          CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
197#  endif
198#if defined key_dynldf_c3d && defined key_dynldf_smag
199                          CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
200#  endif
201
202      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
203      !  Ocean dynamics : hdiv, rot, ssh, e3, wn
204      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
205                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
206      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
207                         CALL wzv           ( kstp )  ! now cross-level velocity
208
209      IF( lk_dynspg_ts ) THEN 
210          ! In case the time splitting case, update almost all momentum trends here:
211          ! Note that the computation of vertical velocity above, hence "after" sea level
212          ! is necessary to compute momentum advection for the rhs of barotropic loop:
213            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
214                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
215            IF( ln_zps .AND. .NOT. ln_isfcav)                               &
216               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
217               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
218            IF( ln_zps .AND.       ln_isfcav)                               &
219               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
220               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
221               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
222
223                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero
224                                  va(:,:,:) = 0.e0
225          IF(  lk_asminc .AND. ln_asmiau .AND. &
226             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment
227          IF( ln_stopack .AND. ln_sppt_dyn ) &
228             &                    CALL dyn_sppt_apply( kstp, 0 )
229          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified)
230          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends
231                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form)
232                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis
233                                  CALL dyn_ldf       ( kstp )   ! lateral mixing
234          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified)
235#if defined key_agrif
236          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
237#endif
238                                  CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
239                                  CALL dyn_spg( kstp, indic )  ! surface pressure gradient
240
241                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated at this stage)
242                                  va_sv(:,:,:) = va(:,:,:)
243
244                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
245          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
246                                  CALL wzv           ( kstp )  ! now cross-level velocity
247      ENDIF
248
249      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
250      ! diagnostics and outputs             (ua, va, tsa used as workspace)
251      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
252      IF( lk_floats  )      CALL flo_stp( kstp )         ! drifting Floats
253      IF( lk_diahth  )      CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
254      IF( .NOT. ln_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
255      IF( lk_diadct  )      CALL dia_dct( kstp )         ! Transports
256      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
257      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis
258                            CALL dia_prod( kstp )        ! ocean model: product diagnostics
259                            CALL dia_wri( kstp )         ! ocean model: outputs
260      !
261      IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
262
263#if defined key_top
264      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
265      ! Passive Tracer Model
266      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
267      IF( lk_asminc .AND. ln_asmiau .AND. lk_bgcinc ) &
268         &               CALL bgc_asm_inc( kstp )     ! biogeochemistry assimilation
269                         CALL trc_stp( kstp )         ! time-stepping
270#endif
271
272
273      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
274      ! Active tracers                              (ua, va used as workspace)
275      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
276                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
277
278      IF( ln_stopack .AND. ln_sppt_tra ) &
279         &                   CALL tra_sppt_apply ( kstp, 0 )
280      IF( lk_asminc .AND. ln_asmiau .AND. &
281        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
282                             CALL tra_sbc    ( kstp )       ! surface boundary condition
283      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
284      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
285      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
286      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
287      IF( ln_bias        )   CALL tra_bias   ( kstp )
288      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
289                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
290      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
291                             CALL tra_ldf    ( kstp )       ! lateral mixing
292
293      IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics
294
295#if defined key_agrif
296      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
297#endif
298      IF( ln_stopack .AND. ln_sppt_tra ) &
299         &                   CALL tra_sppt_apply ( kstp, 1 )
300                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
301
302      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
303         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
304                             CALL tra_nxt( kstp )                ! tracer fields at next time step
305         IF( ln_stopack .AND. ln_sppt_tra ) &
306            &                CALL tra_sppt_apply ( kstp, 2 )
307                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
308            IF( ln_zps .AND. .NOT. ln_isfcav)                                &
309               &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
310               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
311            IF( ln_zps .AND.       ln_isfcav)                                &
312               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
313               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
314               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
315            IF( ln_bias )    CALL dyn_bias( kstp )
316      ELSE                                                  ! centered hpg  (eos then time stepping)
317         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
318                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
319         IF( ln_zps .AND. .NOT. ln_isfcav)                                   &
320               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
321               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
322         IF( ln_zps .AND.       ln_isfcav)                                   & 
323               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
324               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
325               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
326         ENDIF
327         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
328                             CALL tra_nxt( kstp )                ! tracer fields at next time step
329         IF( ln_bias )       CALL dyn_bias( kstp )
330         IF( ln_stopack .AND. ln_sppt_tra )   THEN
331                                CALL tra_sppt_apply ( kstp, 2 )
332                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
333                                IF( ln_zps .AND. .NOT. ln_isfcav)                &
334                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
335                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level
336                                IF( ln_zps .AND.       ln_isfcav)                                   & 
337                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
338                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
339                                &                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
340         ENDIF
341      ENDIF
342
343      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
344      ! Dynamics                                    (tsa used as workspace)
345      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
346      IF( lk_dynspg_ts   )  THEN
347                                                             ! revert to previously computed momentum tendencies
348                                                             ! (not using ua, va as temporary arrays during tracers' update could avoid that)
349                               ua(:,:,:) = ua_sv(:,:,:)
350                               va(:,:,:) = va_sv(:,:,:)
351                                                             ! Revert now divergence and rotational to previously computed ones
352                                                             !(needed because of the time swap in div_cur, at the beginning of each time step)
353                               hdivn(:,:,:) = hdivb(:,:,:)
354                               rotn(:,:,:)  = rotb(:,:,:) 
355
356                               CALL dyn_bfr( kstp )         ! bottom friction
357            IF( ln_stopack .AND. ln_sppt_dyn ) &
358               &               CALL dyn_sppt_apply ( kstp, 1 )
359                               CALL dyn_zdf( kstp )         ! vertical diffusion
360      ELSE
361                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero
362                               va(:,:,:) = 0.e0
363
364        IF(  lk_asminc .AND. ln_asmiau .AND. &
365           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
366        IF( ln_stopack .AND. ln_sppt_dyn ) &
367           &                   CALL dyn_sppt_apply ( kstp, 0 )
368        IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
369        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
370        IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
371                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
372                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
373                               CALL dyn_ldf( kstp )         ! lateral mixing
374        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
375#if defined key_agrif
376        IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
377#endif
378                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
379                               CALL dyn_bfr( kstp )         ! bottom friction
380            IF( ln_stopack .AND. ln_sppt_dyn ) &
381               &               CALL dyn_sppt_apply ( kstp, 1 )
382                               CALL dyn_zdf( kstp )         ! vertical diffusion
383                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
384      ENDIF
385                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
386            IF( ln_stopack .AND. ln_sppt_dyn ) &
387               &               CALL dyn_sppt_apply ( kstp, 2 )
388            IF( ln_stopack .AND. ln_skeb ) &
389               &               CALL skeb_apply ( kstp )
390
391                               CALL ssh_swp( kstp )         ! swap of sea surface height
392      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
393      !
394      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
395      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
396
397#if defined key_agrif
398      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
399      ! AGRIF
400      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
401                               CALL Agrif_Integrate_ChildGrids( stp ) 
402
403      IF ( Agrif_NbStepint().EQ.0 ) THEN
404                               CALL Agrif_Update_Tra()      ! Update active tracers
405                               CALL Agrif_Update_Dyn()      ! Update momentum
406      ENDIF
407#endif
408      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
409      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
410
411      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
412      ! Control
413      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
414                               CALL stp_ctl( kstp, indic )
415      IF( indic < 0        )   THEN
416                               CALL ctl_stop( 'step: indic < 0' )
417                               CALL dia_wri_state( 'output.abort', kstp )
418      ENDIF
419      IF( kstp == nit000   )   THEN
420                 CALL iom_close( numror )     ! close input  ocean restart file
421         IF(lwm) CALL FLUSH    ( numond )     ! flush output namelist oce
422         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
423      ENDIF
424
425      IF( lrst_bias )          CALL bias_wrt     ( kstp )
426
427      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
428      ! Coupled mode
429      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
430      !IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
431      !
432#if defined key_iomput
433      IF( kstp == nitend .OR. indic < 0 ) THEN
434                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
435         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
436      ENDIF
437#endif
438      !
439      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
440      !     
441      !
442   END SUBROUTINE stp
443
444   !!======================================================================
445END MODULE step
Note: See TracBrowser for help on using the repository browser.