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

Last change on this file since 15561 was 15561, checked in by techene, 3 years ago

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