source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/step.F90 @ 10455

Last change on this file since 10455 was 10364, checked in by acc, 2 years ago

Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.

  • Property svn:keywords set to Id
File size: 19.9 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   !!            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   !!----------------------------------------------------------------------
34
35   !!----------------------------------------------------------------------
36   !!   stp             : OPA system time-stepping
37   !!----------------------------------------------------------------------
38   USE step_oce         ! time stepping definition modules
39   !
40   USE iom              ! xIOs server
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   stp   ! called by nemogcm.F90
46
47   !!----------------------------------------------------------------------
48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54#if defined key_agrif
55   RECURSIVE SUBROUTINE stp( )
56      INTEGER             ::   kstp   ! ocean time-step index
57#else
58   SUBROUTINE stp( kstp )
59      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
60#endif
61      !!----------------------------------------------------------------------
62      !!                     ***  ROUTINE stp  ***
63      !!
64      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
65      !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
66      !!              - Time stepping of TRC  (passive tracer eqs.)
67      !!
68      !! ** Method  : -1- Update forcings and data
69      !!              -2- Update ocean physics
70      !!              -3- Compute the t and s trends
71      !!              -4- Update t and s
72      !!              -5- Compute the momentum trends
73      !!              -6- Update the horizontal velocity
74      !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w)
75      !!              -8- Outputs and diagnostics
76      !!----------------------------------------------------------------------
77      INTEGER ::   ji, jj, jk   ! dummy loop indice
78      INTEGER ::   indic        ! error indicator if < 0
79!!gm kcall can be removed, I guess
80      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
81      !! ---------------------------------------------------------------------
82#if defined key_agrif
83      kstp = nit000 + Agrif_Nb_Step()
84      IF( lk_agrif_debug ) THEN
85         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
86         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
87      ENDIF
88      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
89# if defined key_iomput
90      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
91# endif
92#endif
93      !
94      IF( ln_timing )   CALL timing_start('stp')
95      !
96      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
97      ! update I/O and calendar
98      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
99                             indic = 0                ! reset to no error condition
100                             
101      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
102                             CALL iom_init(      cxios_context          )  ! for model grid (including passible AGRIF zoom)
103         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
104      ENDIF
105      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
106                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
107      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
108
109      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
110      ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
111      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
112      IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
113      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
114      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
115                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
116
117      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
118      ! Update stochastic parameters and random T/S fluctuations
119      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
120      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
121      IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
122
123      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124      ! Ocean physics update
125      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
126      !  THERMODYNAMICS
127                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
128                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
129                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
130                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
131
132      !  VERTICAL PHYSICS
133                         CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
134
135      !  LATERAL  PHYSICS
136      !
137      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
138                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
139
140         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
141            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
142            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
143
144         IF( ln_zps .AND.       ln_isfcav)                               &
145            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
146            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
147         IF( ln_traldf_triad ) THEN
148                         CALL ldf_slp_triad( kstp )                       ! before slope for triad operator
149         ELSE     
150                         CALL ldf_slp     ( kstp, rhd, rn2b )             ! before slope for standard operator
151         ENDIF
152      ENDIF
153      !                                                                   ! eddy diffusivity coeff.
154      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp )       !       and/or eiv coeff.
155      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff.
156
157      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
158      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
159      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
160
161                            CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
162      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
163                            CALL wzv           ( kstp )  ! now cross-level velocity
164      IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning
165                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
166                           
167!!jc: fs simplification
168!!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)
169!!                                         but ensures reproductible results
170!!                                         with previous versions using split-explicit free surface         
171            IF( ln_zps .AND. .NOT. ln_isfcav )                               &
172               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
173               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level
174            IF( ln_zps .AND.       ln_isfcav )                                          &
175               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
176               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level
177!!jc: fs simplification
178                           
179                         ua(:,:,:) = 0._wp            ! set dynamics trends to zero
180                         va(:,:,:) = 0._wp
181
182      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
183               &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
184      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends
185#if defined key_agrif
186      IF(.NOT. Agrif_Root())  & 
187               &         CALL Agrif_Sponge_dyn        ! momentum sponge
188#endif
189                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form)
190                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
191                         CALL dyn_ldf       ( kstp )  ! lateral mixing
192      IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
193                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
194                         CALL dyn_spg       ( kstp )  ! surface pressure gradient
195
196                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well
197      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
198                            CALL div_hor    ( kstp )              ! Horizontal divergence  (2nd call in time-split case)
199         IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
200                            CALL wzv        ( kstp )              ! now cross-level velocity
201         IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning
202      ENDIF
203     
204                         CALL dyn_zdf       ( kstp )  ! vertical diffusion
205
206      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
207      ! cool skin
208      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
209      IF ( ln_diurnal )  CALL stp_diurnal( kstp )
210     
211      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
212      ! diagnostics and outputs
213      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
214      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats
215      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
216      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
217      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports
218                         CALL dia_ar5 ( kstp )        ! ar5 diag
219      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
220                         CALL dia_wri ( kstp )        ! ocean model: outputs
221      !
222      IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output
223     
224#if defined key_top
225      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
226      ! Passive Tracer Model
227      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
228                         CALL trc_stp       ( kstp )  ! time-stepping
229#endif
230
231      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
232      ! Active tracers                             
233      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
234                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
235
236      IF(  lk_asminc .AND. ln_asmiau .AND. &
237         & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment
238                         CALL tra_sbc       ( kstp )  ! surface boundary condition
239      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr
240      IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux
241      IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme
242      IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends
243      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
244#if defined key_agrif
245      IF(.NOT. Agrif_Root())  & 
246               &         CALL Agrif_Sponge_tra        ! tracers sponge
247#endif
248                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection
249      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
250      IF( lrst_oce .AND. ln_zdfosm ) &
251           &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
252                         CALL tra_ldf       ( kstp )  ! lateral mixing
253
254!!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
255      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
256!!gm
257                         CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields
258      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection
259
260      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
261      ! Set boundary conditions and Swap
262      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
263!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
264!!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
265!!    If so:
266!!    (i) no need to call agrif update at initialization time
267!!    (ii) no need to update "before" fields
268!!
269!!    Apart from creating new tra_swp/dyn_swp routines, this however:
270!!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
271!!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
272!!    e.g. a shift of the feedback interface inside child domain.
273!!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
274!!    place.
275!!
276!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
277                         CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
278                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
279                         CALL ssh_swp       ( kstp )  ! swap of sea surface height
280      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
281      !
282      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics
283
284!!gm : This does not only concern the dynamics ==>>> add a new title
285!!gm2: why ouput restart before AGRIF update?
286!!
287!!jc: That would be better, but see comment above
288!!
289      IF( lrst_oce   )   CALL rst_write    ( kstp )   ! write output ocean restart file
290      IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
291
292#if defined key_agrif
293      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
294      ! AGRIF
295      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
296                         CALL Agrif_Integrate_ChildGrids( stp )  ! allows to finish all the Child Grids before updating
297
298                         IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components
299#endif
300      IF( ln_diaobs  )   CALL dia_obs      ( kstp )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
301
302      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
303      ! Control
304      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
305                         CALL stp_ctl      ( kstp, indic )
306                         
307      IF( kstp == nit000 ) THEN                          ! 1st time step only
308                                        CALL iom_close( numror )   ! close input  ocean restart file
309         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
310         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
311      ENDIF
312
313      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
314      ! Coupled mode
315      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
316!!gm why lk_oasis and not lk_cpl ????
317      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
318      !
319#if defined key_iomput
320      IF( kstp == nitend .OR. indic < 0 ) THEN
321                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
322                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
323         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
324      ENDIF
325#endif
326      !
327      IF( ln_timing )   CALL timing_stop('stp')
328      !
329   END SUBROUTINE stp
330   
331   !!======================================================================
332END MODULE step
Note: See TracBrowser for help on using the repository browser.