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 @ 15373

Last change on this file since 15373 was 15373, checked in by techene, 9 months ago

#2715 RK3: adapt trc surface boundary management

File size: 20.6 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      IF( kstg == 1 ) THEN
211         CALL trc_stp_start( kstp, Kbb, Kmm, Krhs, Kaa )
212      ENDIF
213      !
214      DO jn = 1, jptra
215         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
216            tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero
217         END_3D
218      END DO
219      !                                         !==  advection of passive tracers  ==!
220      rDt_trc = rDt
221!!st
222      CALL trc_sbc_RK3( kstp,      Kmm, tr, Krhs, kstg )              ! surface boundary condition
223      !
224      CALL trc_adv    ( kstp, Kbb, Kmm, tr, Krhs, zaU, zaV, ww )      ! horizontal & vertical advection
225      !
226      !
227      SELECT CASE( kstg )
228      !                    !-------------------!
229      CASE ( 1 , 2 )       !==  Stage 1 & 2  ==!   stg1:  Kbb = N  ;  Kaa = N+1/3
230         !                 !-------------------!   stg2:  Kbb = N  ;  Kmm = N+1/3  ;  Kaa = N+1/2
231         !
232         !                                      !==  time integration  ==!   ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
233         DO jn = 1, jptra
234            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
235               ze3Tb = e3t(ji,jj,jk,Kbb) * tr(ji,jj,jk,jn,Kbb )
236               ze3Tr = e3t(ji,jj,jk,Kmm) * tr(ji,jj,jk,jn,Krhs)
237               z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
238               tr(ji,jj,jk,jn,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
239            END_3D
240         END DO
241         !                 !---------------!
242      CASE ( 3 )           !==  Stage 3  ==!   add all RHS terms but advection (=> Kbb only)
243         !                 !---------------!
244         CALL trc_sms    ( kstp, Kbb, Kbb, Krhs      )       ! tracers: sinks and sources
245         CALL trc_trp    ( kstp, Kbb, Kmm, Krhs, Kaa )       ! transport of passive tracers (without advection)
246         !
247         !
248         CALL trc_stp_end( kstp, Kbb, Kmm,       Kaa )
249         !
250      END SELECT
251# endif
252
253      !                       !==  T-S Tracers  ==!
254
255!===>>>>>> Modify tra_adv_...  routines so that Krhs to zero useless
256      DO jn = 1, jpts
257         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
258            ts(ji,jj,jk,jn,Krhs) = 0._wp                                    ! set tracer trends to zero
259         END_3D
260      END DO
261
262!===>>> CAUTION here may be without GM velocity but stokes drift required ! 0 barotropic divergence for GM  != 0 barotropic divergence for SD
263!!st consistence 2D / 3D - flux de masse
264      CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww )       ! hor. + vert. advection  ==> RHS
265
266!===>>>>>> stg1&2:  Verify the necessity of these trends (we may need it as there are in the RHS of dynspg_ts ?)
267!!gm ====>>>>   needed for heat and salt fluxes associated with mass/volume flux
268                        CALL tra_sbc_RK3( kstp,      Kmm, ts, Krhs, kstg )   ! surface boundary condition
269
270      IF( ln_isf )      CALL tra_isf    ( kstp,      Kmm, ts, Krhs )   ! ice shelf heat flux
271      IF( ln_traqsr )   CALL tra_qsr    ( kstp,      Kmm, ts, Krhs )   ! penetrative solar radiation qsr
272!!gm
273
274      !                               
275!!gm ===>>>>>>  Verify the necessity of these trends  at stages 1 and 2
276!           (we may need it as they are in the RHS of dynspg_ts ?)
277!      IF(  lk_asminc .AND. ln_asmiau ) THEN               ! apply assimilation increment
278!         IF( ln_dyninc )   CALL dyn_asm_inc( kstp, Kbb, Kmm, uu, vv, Krhs )   ! dynamics   ==> RHS
279!         IF( ln_trainc )   CALL tra_asm_inc( kstp, Kbb, Kmm, ts    , Krhs )   ! tracers    ==> RHS
280!      ENDIF
281!!gm  end Verif
282
283      !
284      SELECT CASE( kstg )
285      !                    !-------------------!
286      CASE ( 1 , 2 )       !==  Stage 1 & 2  ==!   stg1:  Kbb = N  ;  Kaa = N+1/3
287         !                 !-------------------!   stg2:  Kbb = N  ;  Kmm = N+1/3  ;  Kaa = N+1/2
288         !
289         !                                      !==  time integration  ==!   ∆t = rn_Dt/3 (stg1) or rn_Dt/2 (stg2)
290         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
291            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
292               uu(ji,jj,jk,Kaa) = ( uu(ji,jj,jk,Kbb) + rDt * uu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
293               vv(ji,jj,jk,Kaa) = ( vv(ji,jj,jk,Kbb) + rDt * vv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
294            END_3D
295         ELSE                                      ! applied on thickness weighted velocity
296            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
297               uu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb )  &
298                  &                 + rDt * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Krhs)  ) &
299                  &                       / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
300               vv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb )  &
301                  &                 + rDt * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Krhs)  ) &
302                  &                       / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
303            END_3D
304         ENDIF
305         !
306         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
307            ze3Tb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_tem,Kbb )
308            ze3Sb = e3t(ji,jj,jk,Kbb) * ts(ji,jj,jk,jp_sal,Kbb )
309            ze3Tr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Krhs)
310            ze3Sr = e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Krhs)
311            z1_e3t= 1._wp / e3t(ji,jj,jk, Kaa)
312            ts(ji,jj,jk,jp_tem,Kaa) = ( ze3Tb + rDt * ze3Tr*tmask(ji,jj,jk) ) * z1_e3t
313            ts(ji,jj,jk,jp_sal,Kaa) = ( ze3Sb + rDt * ze3Sr*tmask(ji,jj,jk) ) * z1_e3t
314         END_3D
315         !
316         !
317         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,Kaa), clinfo1='stp stg   - Ua: ', mask1=umask,   &
318            &                                  tab3d_2=vv(:,:,:,Kaa), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
319         !
320         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Kaa), clinfo1='stp stg   - Ta: ', mask1=tmask,   &
321            &                                  tab3d_2=ts(:,:,:,jp_sal,Kaa), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
322         !
323         !                 !---------------!
324      CASE ( 3 )           !==  Stage 3  ==!   add all remaining RHS terms
325         !                 !---------------!
326         !
327         !                                      !==  complete the momentum RHS ==!   except ZDF (implicit)
328         !                                                   ! lateral mixing                    ==> RHS
329                            CALL dyn_ldf( kstp, Kbb, Kmm, uu, vv, Krhs )
330         !                                                   ! OSMOSIS non-local velocity fluxes ==> RHS
331         IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Kmm, uu, vv, Krhs )
332         !                   
333         IF( ln_bdy     ) THEN                               ! bdy damping trends     ==> RHS
334                            CALL bdy_dyn3d_dmp ( kstp, Kbb, uu, vv, Krhs )
335                            CALL bdy_tra_dmp   ( kstp, Kbb, ts    , Krhs )
336         ENDIF
337
338# if defined key_agrif
339         IF(.NOT. Agrif_Root() ) THEN                        ! AGRIF:   sponge ==> momentum and tracer RHS
340            CALL Agrif_Sponge_dyn
341            CALL Agrif_Sponge_tra
342         ENDIF
343# endif
344         !                                      !==  complete the tracers RHS  ==!   except ZDF (implicit)
345         !                                            !*  T-S Tracer  *!
346         !
347                            CALL tra_ldf( kstp, Kbb, Kmm, ts, Krhs )  ! lateral mixing
348         IF( ln_trabbc  )   CALL tra_bbc( kstp,      Kmm, ts, Krhs )  ! bottom heat flux
349         IF( ln_trabbl  )   CALL tra_bbl( kstp, Kbb, Kmm, ts, Krhs )  ! advective (and/or diffusive) bottom boundary layer scheme
350         IF( ln_tradmp  )   CALL tra_dmp( kstp, Kbb, Kmm, ts, Krhs )  ! internal damping trends
351
352         IF( ln_zdfmfc  )   CALL tra_mfc( kstp, Kbb,      ts, Krhs )  ! Mass Flux Convection
353         IF( ln_zdfosm  ) THEN
354                            CALL tra_osm( kstp,      Kmm, ts, Krhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
355            IF( lrst_oce )  CALL osm_rst( kstp,      Kmm, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
356         ENDIF
357         !
358         !                                      !==  DYN & TRA time integration + ZDF  ==!   ∆t = rDt
359         !
360                            CALL dyn_zdf( kstp, Kbb, Kmm, Krhs, uu, vv, Kaa  )  ! vertical diffusion and time integration
361                            CALL tra_zdf( kstp, Kbb, Kmm, Krhs, ts    , Kaa  )  ! vertical mixing and after tracer fields
362         IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Kmm, Krhs, ts    , Kaa  )  ! update after fields by non-penetrative convection
363         !
364      END SELECT     
365      !                                         !==  correction of the barotropic (all stages)  ==!    at Kaa = N+1/3, N+1/2 or N+1
366      !                                                           ! barotropic velocity correction
367      zub(A2D(0)) = uu_b(A2D(0),Kaa) - SUM( e3u_0(A2D(0),:)*uu(A2D(0),:,Kaa), 3 ) * r1_hu_0(A2D(0))
368      zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0))
369      !
370      DO jk = 1, jpkm1                                            ! corrected horizontal velocity
371         uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk)
372         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk)
373      END DO
374         
375
376      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
377      ! Set boundary conditions
378      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
379      !
380# if defined key_agrif
381            CALL Agrif_tra( kstp, kstg )             !* AGRIF zoom boundaries
382            CALL Agrif_dyn( kstp, kstg )
383# endif
384      !                                              !* local domain boundaries  (T-point, unchanged sign)
385      CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &
386         &                             , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. )
387      !
388      !                                              !* BDY open boundaries
389      IF( ln_bdy )   THEN
390                               CALL bdy_tra( kstp, Kbb, ts,     Kaa )
391         IF( ln_dynspg_exp )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa )
392         IF( ln_dynspg_ts  )   CALL bdy_dyn( kstp, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
393      ENDIF
394      !
395      IF( ln_timing )   CALL timing_stop('stp_RK3_stg')
396      !
397   END SUBROUTINE stp_RK3_stg
398
399#else
400   !!----------------------------------------------------------------------
401   !!   default option             EMPTY MODULE           qco not activated
402   !!----------------------------------------------------------------------
403#endif
404   
405   !!======================================================================
406END MODULE stprk3_stg
Note: See TracBrowser for help on using the repository browser.