source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90 @ 10001

Last change on this file since 10001 was 10001, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step I.1 and I.2 (see wiki page)

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