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.
stprk3_stg.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90 @ 15510

Last change on this file since 15510 was 15510, checked in by techene, 8 months ago

#2715 RK3 & TOP: musl adv scheme is only called at stage 3

File size: 21.4 KB
Line 
1MODULE stprk3_stg
2   !!======================================================================
3   !!                       ***  MODULE stprk3_stg  ***
4   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
5   !!                   using a 3rd order Runge-Kutta  with fixed or quasi-eulerian coordinate
6   !!======================================================================
7   !! History :  4.5  !  2021-01  (S. Techene, G. Madec, N. Ducousso, F. Lemarie)  Original code
8   !!   NEMO     
9   !!----------------------------------------------------------------------
10#if defined key_qco   ||   defined key_linssh
11   !!----------------------------------------------------------------------
12   !!   'key_qco'                        Quasi-Eulerian vertical coordinate
13   !!                          OR
14   !!   'key_linssh                       Fixed in time vertical coordinate
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!    stp_RK3_stg  : NEMO 3rd order Runge-Kutta stage with qco or linssh
19   !!----------------------------------------------------------------------
20   USE step_oce       ! time stepping used modules
21   USE domqco         ! quasi-eulerian coordinate      (dom_qco_r3c routine)
22   USE dynspg_ts, ONLY:   un_adv , vn_adv   ! advective transport from N to N+1
23   USE bdydyn         ! ocean open boundary conditions (define bdy_dyn)
24# if defined key_top
25   USE trc            ! ocean passive tracers variables
26   USE trcadv         ! passive tracers advection      (trc_adv routine)
27   USE trcsms         ! passive tracers source and sink
28   USE trctrp         ! passive tracers transport
29   USE trcsbc         ! passive tracers surface boundary condition !!st WARNING USELESS TO BE REMOVED
30   USE trcbdy         ! passive tracers transport open boundary
31   USE trcstp_rk3
32# endif
33# if defined key_agrif
34   USE agrif_oce_interp
35# endif
36   
37   !
38   USE prtctl         ! print control
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   stp_RK3_stg   ! called by nemogcm.F90
44
45   REAL(wp) ::   r1_2 = 1._wp / 2._wp
46   REAL(wp) ::   r1_3 = 1._wp / 3._wp
47   REAL(wp) ::   r2_3 = 2._wp / 3._wp
48
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssha         ! sea-surface height  at N+1
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b   ! barotropic velocity at N+1
51
52   !! * Substitutions
53#  include "do_loop_substitute.h90"
54#  include "domzgr_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
57   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
58   !! Software governed by the CeCILL license (see ./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE stp_RK3_stg( kstg, kstp, Kbb, Kmm, Krhs, Kaa )
63      !!----------------------------------------------------------------------
64      !!                     ***  ROUTINE stp_RK3_stg  ***
65      !!
66      !! ** Purpose : - stage of RK3 time stepping of OCE and TOP
67      !!
68      !! ** Method  :   input: computed in dynspg_ts
69      !!              ssh             shea surface height at N+1           (oce.F90)
70      !!              (uu_b,vv_b)     barotropic velocity at N, N+1        (oce.F90)
71      !!              (un_adv,vn_adv) barotropic transport from N to N+1   (dynspg_ts.F90)
72      !!              ,
73      !!              -1- set ssh(Naa) (Naa=N+1/3, N+1/2, or N)
74      !!              -2- set the advective velocity (zadU,zaV)
75      !!              -4- Compute the after (Naa) T-S
76      !!              -5- Update now
77      !!              -6- Update the horizontal velocity
78      !!----------------------------------------------------------------------
79      INTEGER, INTENT(in) ::   kstg                        ! RK3 stage
80      INTEGER, INTENT(in) ::   kstp, Kbb, Kmm, Krhs, Kaa   ! ocean time-step and time-level indices
81      !
82      INTEGER  ::   ji, jj, jk, jn, jtile                  ! dummy loop indices
83      REAL(wp) ::   ze3Tb, ze3Sb, z1_e3t     ! local scalars
84      REAL(wp) ::   ze3Tr, ze3Sr             !   -      -
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zaU, zaV       ! advective horizontal velocity
86      REAL(wp), DIMENSION(jpi,jpj)     ::   zub, zvb       ! advective transport
87      !! ---------------------------------------------------------------------
88      !
89      IF( ln_timing )   CALL timing_start('stp_RK3_stg')
90      !
91      IF( kstp == nit000 ) THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'stp_RK3_stg : Runge Kutta 3rd order at stage ', kstg
94         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
95      ENDIF
96      !
97      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
98      !  ssh, uu_b, vv_b, and  ssh/h0 at Kaa
99      !  3D advective velocity at Kmm
100      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
101      !     
102      SELECT CASE( kstg )
103      !                    !---------------!
104      CASE ( 1 )           !==  Stage 1  ==!   Kbb = Kmm = N  ;  Kaa = N+1/3
105         !                 !---------------!
106         !
107         ALLOCATE( ssha(jpi,jpj) , ua_b(jpi,jpj) , va_b(jpi,jpj) )
108         !
109         rDt = r1_3 * rn_Dt            ! set time-step : rn_Dt/3
110         r1_Dt = 1._wp / rDt
111         !
112         ssha(:,:) = ssh (:,:,Kaa)     ! save ssh, uu_b, vv_b at N+1  (computed in dynspg_ts)
113         ua_b(:,:) = uu_b(:,:,Kaa)
114         va_b(:,:) = vv_b(:,:,Kaa)
115         !                             ! interpolated ssh and (uu_b,vv_b) at Kaa (N+1/3)
116         ssh (:,:,Kaa) = r2_3 * ssh (:,:,Kbb) + r1_3 * ssha(:,:)
117         uu_b(:,:,Kaa) = r2_3 * uu_b(:,:,Kbb) + r1_3 * ua_b(:,:)
118         vv_b(:,:,Kaa) = r2_3 * vv_b(:,:,Kbb) + r1_3 * va_b(:,:)
119         !
120         !                 !---------------!
121      CASE ( 2 )           !==  Stage 2  ==!   Kbb = N   ;   Kmm = N+1/3   ;   Kaa = N+1/2
122         !                 !---------------!
123         !
124         rDt = r1_2 * rn_Dt            ! set time-step : rn_Dt/2
125         r1_Dt = 1._wp / rDt
126         !
127         !                             ! set ssh and (uu_b,vv_b) at N+1/2  (Kaa)
128         ssh (:,:,Kaa) = r1_2 * ( ssh (:,:,Kbb) + ssha(:,:) )
129         uu_b(:,:,Kaa) = r1_2 * ( uu_b(:,:,Kbb) + ua_b(:,:) )
130         vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
131         !
132         !                 !---------------!
133      CASE ( 3 )           !==  Stage 3  ==!   Kbb = N   ;   Kmm = N+1/2   ;   Kaa = N+1
134         !                 !---------------!
135         !
136         rDt = rn_Dt                   ! set time-step : rn_Dt
137         r1_Dt = 1._wp / rDt
138         !
139         ssh (:,:,Kaa) = ssha(:,:)     ! recover ssh and (uu_b,vv_b) at N + 1
140         uu_b(:,:,Kaa) = ua_b(:,:)                     
141         vv_b(:,:,Kaa) = va_b(:,:)
142         !
143         DEALLOCATE( ssha , ua_b , va_b )
144         !
145      END SELECT
146      !
147      !                     !==  ssh/h0 ratio at Kaa  ==!
148      !
149      IF( .NOT.lk_linssh )   CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! "after" ssh/h_0 ratio at t,u,v-column
150      !
151      !
152      !                     !==  advective velocity at Kmm  ==!
153      !
154      !                                            !- horizontal components -!   (zaU,zaV)
155      zub(:,:) = un_adv(:,:)*r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)    ! barotropic velocity correction
156      zvb(:,:) = vn_adv(:,:)*r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)
157      DO jk = 1, jpkm1                                         ! horizontal advective velocity
158         zaU(:,:,jk) = uu(:,:,jk,Kmm) + zub(:,:)*umask(:,:,jk)
159         zaV(:,:,jk) = vv(:,:,jk,Kmm) + zvb(:,:)*vmask(:,:,jk)
160      END DO
161      !                                            !- vertical components -!   ww
162                         CALL wzv  ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )     ! ww cross-level velocity
163
164!!st      IF( ln_zad_Aimp )  CALL wAimp( kstp, ww, wi       )     ! Adaptive-implicit vertical advection partitioning
165      !                                                            !==>>>  implicite for stages 1 & 2 ????
166      !
167
168      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
169      !  RHS of ocean dynamics : ADV + VOR/COR + HPG (+ ASM )      <<<===  Question:  Stokes drift ?
170      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
171      !
172      uu(:,:,:,Krhs) = 0._wp           ! set dynamics trends to zero
173      vv(:,:,:,Krhs) = 0._wp
174      !
175!===>>>>>> Modify dyn_adv_... dyn_keg routines so that Krhs to zero useless
176      !                                         ! advection (VF or FF)  ==> RHS
177      CALL dyn_adv( kstp, Kbb, Kmm      , uu, vv, Krhs, zaU, zaV, ww )
178      !                                         ! Coriolis / vorticity  ==> RHS
179      CALL dyn_vor( kstp,      Kmm      , uu, vv, Krhs )
180      !
181!!gm à appeler que pour ln_zad_Aimp=T   et en ne faisant que wi  par zdf
182!!                                              ! ZAD (implicit part)   ==> RHS
183!!    CALL dyn_zdf    ( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa  )               
184
185!===>>>>>> Modify dyn_hpg & dyn_hpg_...  routines : rhd computed in dyn_hpg and pass in argument to dyn_hpg_...
186
187      CALL eos    ( ts, Kmm, rhd )              ! Kmm in situ density anomaly for hpg computation
188
189!!gm end
190      CALL dyn_hpg( kstp,      Kmm      , uu, vv, Krhs )
191      !
192!!gm ===>>>>>> Probably useless since uu_b(Kaa) will be imposed at the end of stage 1 and 2
193!                   but may be necessary in stage 3 due to implicite in dynzdf.
194!                   except if my idea for the matrice construction is OK !
195!      !                                         ! grad_h of ps          ==> RHS       
196!      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
197!         uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) - grav * ( ssh(ji+1,jj  ,Kmm) - ssh(ji,jj,Kmm) )
198!         vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) - grav * ( ssh(ji  ,jj+1,Kmm) - ssh(ji,jj,Kmm) )
199!      END_3D
200!!gm     
201     
202      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
203      ! RHS of tracers : ADV only using (zaU,zaV,ww)
204      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
205
206# if defined key_top
207      !
208      !                       !==  Passive Tracer  ==!
209      !
210      SELECT CASE( kstg )
211      !                    !-------------------!
212      CASE ( 1 , 2 )       !==  Stage 1 & 2  ==!   stg1:  Kbb = N  ;  Kaa = N+1/3
213         !                 !-------------------!   stg2:  Kbb = N  ;  Kmm = N+1/3  ;  Kaa = N+1/2
214         !
215         IF( kstg == 1 ) THEN
216            CALL trc_stp_start( kstp, Kbb, Kmm, Krhs, Kaa )
217         ENDIF
218         !
219         IF(.NOT. ln_trcadv_mus ) THEN
220            !
221            DO jn = 1, jptra
222               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
223               tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero
224               END_3D
225            END DO
226            !                                      !==  advection of passive tracers  ==!
227            rDt_trc = rDt
228            !
229            CALL trc_sbc_RK3( kstp,      Kmm, tr, Krhs, kstg )              ! surface boundary condition
230            !
231            CALL trc_adv    ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww )      ! horizontal & vertical advection
232            !
233            !                                      !==  time integration  ==!   ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
234            DO jn = 1, jptra
235               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
236               ze3Tb = e3t(ji,jj,jk,Kbb) * tr(ji,jj,jk,jn,Kbb )
237               ze3Tr = e3t(ji,jj,jk,Kmm) * tr(ji,jj,jk,jn,Krhs)
238               z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
239               tr(ji,jj,jk,jn,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
240               END_3D
241            END DO
242            !
243         ENDIF
244         !                 !---------------!
245      CASE ( 3 )           !==  Stage 3  ==!   add all RHS terms but advection (=> Kbb only)
246         !                 !---------------!
247         !
248         DO jn = 1, jptra
249            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
250            tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero
251            END_3D
252         END DO
253         !                                         !==  advection of passive tracers  ==!
254         rDt_trc = rDt
255         !
256         CALL trc_sbc_RK3( kstp,      Kmm, tr, Krhs, kstg )              ! surface boundary condition
257         !
258         CALL trc_adv    ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww )      ! horizontal & vertical advection
259         !
260         CALL trc_sms    ( kstp, Kbb, Kbb, Krhs      )       ! tracers: sinks and sources
261         CALL trc_trp    ( kstp, Kbb, Kmm, Krhs, Kaa )       ! transport of passive tracers (without advection)
262         !
263         !
264         CALL trc_stp_end( kstp, Kbb, Kmm,       Kaa )
265         !
266      END SELECT
267# endif
268
269      !                       !==  T-S Tracers  ==!
270
271!===>>>>>> Modify tra_adv_...  routines so that Krhs to zero useless
272      DO jn = 1, jpts
273         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
274            ts(ji,jj,jk,jn,Krhs) = 0._wp                                    ! set tracer trends to zero
275         END_3D
276      END DO
277
278!===>>> CAUTION here may be without GM velocity but stokes drift required ! 0 barotropic divergence for GM  != 0 barotropic divergence for SD
279!!st consistence 2D / 3D - flux de masse
280      CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww, kstg )       ! hor. + vert. advection  ==> RHS
281
282!===>>>>>> stg1&2:  Verify the necessity of these trends (we may need it as there are in the RHS of dynspg_ts ?)
283!!gm ====>>>>   needed for heat and salt fluxes associated with mass/volume flux
284                        CALL tra_sbc_RK3( kstp,      Kmm, ts, Krhs, kstg )   ! surface boundary condition
285
286      IF( ln_isf )      CALL tra_isf    ( kstp,      Kmm, ts, Krhs )   ! ice shelf heat flux
287      IF( ln_traqsr )   CALL tra_qsr    ( kstp,      Kmm, ts, Krhs )   ! penetrative solar radiation qsr
288!!gm
289
290      !                               
291!!gm ===>>>>>>  Verify the necessity of these trends  at stages 1 and 2
292!           (we may need it as they are in the RHS of dynspg_ts ?)
293!      IF(  lk_asminc .AND. ln_asmiau ) THEN               ! apply assimilation increment
294!         IF( ln_dyninc )   CALL dyn_asm_inc( kstp, Kbb, Kmm, uu, vv, Krhs )   ! dynamics   ==> RHS
295!         IF( ln_trainc )   CALL tra_asm_inc( kstp, Kbb, Kmm, ts    , Krhs )   ! tracers    ==> RHS
296!      ENDIF
297!!gm  end Verif
298
299      !
300      SELECT CASE( kstg )
301      !                    !-------------------!
302      CASE ( 1 , 2 )       !==  Stage 1 & 2  ==!   stg1:  Kbb = N  ;  Kaa = N+1/3
303         !                 !-------------------!   stg2:  Kbb = N  ;  Kmm = N+1/3  ;  Kaa = N+1/2
304         !
305         !                                      !==  time integration  ==!   ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
306         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
307            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
308               uu(ji,jj,jk,Kaa) = ( uu(ji,jj,jk,Kbb) + rDt * uu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
309               vv(ji,jj,jk,Kaa) = ( vv(ji,jj,jk,Kbb) + rDt * vv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
310            END_3D
311         ELSE                                      ! applied on thickness weighted velocity
312            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
313               uu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb )  &
314                  &                 + rDt * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Krhs)  ) &
315                  &                       / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
316               vv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb )  &
317                  &                 + rDt * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Krhs)  ) &
318                  &                       / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
319            END_3D
320         ENDIF
321         !
322         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
323            ze3Tb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_tem,Kbb )
324            ze3Sb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_sal,Kbb )
325            ze3Tr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Krhs)
326            ze3Sr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Krhs)
327            z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
328            ts(ji,jj,jk,jp_tem,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
329            ts(ji,jj,jk,jp_sal,Kaa) = ( ze3Sb + rDt * ze3Sr*tmask(ji,jj,jk) ) * z1_e3t
330         END_3D
331         !
332         !
333         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,Kaa), clinfo1='stp stg   - Ua: ', mask1=umask,   &
334            &                                  tab3d_2=vv(:,:,:,Kaa), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
335         !
336         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kaa), clinfo1='stp stg   - Ta: ', mask1=tmask,   &
337            &                                  tab3d_2=ts(:,:,:,jp_sal,Kaa), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
338         !
339         !                 !---------------!
340      CASE ( 3 )           !==  Stage 3  ==!   add all remaining RHS terms
341         !                 !---------------!
342         !
343         !                                      !==  complete the momentum RHS ==!   except ZDF (implicit)
344         !                                                   ! lateral mixing                    ==> RHS
345                            CALL dyn_ldf( kstp, Kbb, Kmm, uu, vv, Krhs )
346         !                                                   ! OSMOSIS non-local velocity fluxes ==> RHS
347         IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Kmm, uu, vv, Krhs )
348         !                   
349         IF( ln_bdy     ) THEN                               ! bdy damping trends     ==> RHS
350                            CALL bdy_dyn3d_dmp ( kstp, Kbb, uu, vv, Krhs )
351                            CALL bdy_tra_dmp   ( kstp, Kbb, ts    , Krhs )
352         ENDIF
353
354# if defined key_agrif
355         IF(.NOT. Agrif_Root() ) THEN                        ! AGRIF:   sponge ==> momentum and tracer RHS
356            CALL Agrif_Sponge_dyn
357            CALL Agrif_Sponge_tra
358         ENDIF
359# endif
360         !                                      !==  complete the tracers RHS  ==!   except ZDF (implicit)
361         !                                            !*  T-S Tracer  *!
362         !
363                            CALL tra_ldf( kstp, Kbb, Kmm, ts, Krhs )  ! lateral mixing
364         IF( ln_trabbc  )   CALL tra_bbc( kstp,      Kmm, ts, Krhs )  ! bottom heat flux
365         IF( ln_trabbl  )   CALL tra_bbl( kstp, Kbb, Kmm, ts, Krhs )  ! advective (and/or diffusive) bottom boundary layer scheme
366         IF( ln_tradmp  )   CALL tra_dmp( kstp, Kbb, Kmm, ts, Krhs )  ! internal damping trends
367
368         IF( ln_zdfmfc  )   CALL tra_mfc( kstp, Kbb,      ts, Krhs )  ! Mass Flux Convection
369         IF( ln_zdfosm  ) THEN
370                            CALL tra_osm( kstp,      Kmm, ts, Krhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
371            IF( lrst_oce )  CALL osm_rst( kstp,      Kmm, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
372         ENDIF
373         !
374         !                                      !==  DYN & TRA time integration + ZDF  ==!   ∆t = rDt
375         !
376                            CALL dyn_zdf( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa  )  ! vertical diffusion and time integration
377                            CALL tra_zdf( kstp, Kbb, Kmm, Krhs, ts    , Kaa  )  ! vertical mixing and after tracer fields
378         IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Kmm, Krhs, ts    , Kaa  )  ! update after fields by non-penetrative convection
379         !
380      END SELECT     
381      !                                         !==  correction of the barotropic (all stages)  ==!    at Kaa = N+1/3, N+1/2 or N+1
382      !                                                           ! barotropic velocity correction
383      zub(A2D(0)) = uu_b(A2D(0),Kaa) - SUM( e3u_0(A2D(0),:)*uu(A2D(0),:,Kaa), 3 ) * r1_hu_0(A2D(0))
384      zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0))
385      !
386      DO jk = 1, jpkm1                                            ! corrected horizontal velocity
387         uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk)
388         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk)
389      END DO
390         
391
392      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
393      ! Set boundary conditions
394      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
395      !
396# if defined key_agrif
397            CALL Agrif_tra( kstp, kstg )             !* AGRIF zoom boundaries
398            CALL Agrif_dyn( kstp, kstg )
399# endif
400      !                                              !* local domain boundaries  (T-point, unchanged sign)
401      CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &
402         &                             , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. )
403      !
404      !                                              !* BDY open boundaries
405      IF( ln_bdy )   THEN
406                               CALL bdy_tra( kstp, Kbb, ts,     Kaa )
407         IF( ln_dynspg_exp )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa )
408         IF( ln_dynspg_ts  )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
409      ENDIF
410      !
411      IF( ln_timing )   CALL timing_stop('stp_RK3_stg')
412      !
413   END SUBROUTINE stp_RK3_stg
414
415#else
416   !!----------------------------------------------------------------------
417   !!   default option             EMPTY MODULE           qco not activated
418   !!----------------------------------------------------------------------
419#endif
420   
421   !!======================================================================
422END MODULE stprk3_stg
Note: See TracBrowser for help on using the repository browser.