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.
dynspg_ts.F90 in NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynspg_ts.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 4 years ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp -> rn_atfp (use namelist parameter everywhere)
rdtbt -> rDt_e
nn_baro -> nn_e
rn_scal_load -> rn_load
rau0 -> rho0

  • Property svn:keywords set to Id
File size: 78.2 KB
RevLine 
[358]1MODULE dynspg_ts
[9023]2
[12377]3   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
[9023]4
[358]5   !!======================================================================
[6140]6   !!                   ***  MODULE  dynspg_ts  ***
7   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme
8   !!======================================================================
[1502]9   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code
10   !!              -   ! 2005-11  (V. Garnier, G. Madec)  optimization
11   !!              -   ! 2006-08  (S. Masson)  distributed restart using iom
12   !!             2.0  ! 2007-07  (D. Storkey) calls to BDY routines
13   !!              -   ! 2008-01  (R. Benshila)  change averaging method
14   !!             3.2  ! 2009-07  (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]15   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]16   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
[4292]17   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
18   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
[5930]19   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
[7646]20   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
[9019]21   !!             4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
[2724]22   !!---------------------------------------------------------------------
[6140]23
[358]24   !!----------------------------------------------------------------------
[6140]25   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme
26   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme
27   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not)
28   !!   ts_rst         : read/write time-splitting fields in restart file
[358]29   !!----------------------------------------------------------------------
30   USE oce             ! ocean dynamics and tracers
31   USE dom_oce         ! ocean space and time domain
[888]32   USE sbc_oce         ! surface boundary condition: ocean
[12377]33   USE isf_oce         ! ice shelf variable (fwfisf)
[9019]34   USE zdf_oce         ! vertical physics: variables
35   USE zdfdrg          ! vertical physics: top/bottom drag coef.
[6140]36   USE sbcapr          ! surface boundary condition: atmospheric pressure
37   USE dynadv    , ONLY: ln_dynadv_vec
[9528]38   USE dynvor          ! vortivity scheme indicators
[358]39   USE phycst          ! physical constants
40   USE dynvor          ! vorticity term
[6152]41   USE wet_dry         ! wetting/drying flux limter
[7646]42   USE bdy_oce         ! open boundary
[10481]43   USE bdyvol          ! open boundary volume conservation
[5930]44   USE bdytides        ! open boundary condition data
[3294]45   USE bdydyn2d        ! open boundary conditions on barotropic variables
[12377]46   USE tide_mod        !
[7646]47   USE sbcwave         ! surface wave
[9019]48#if defined key_agrif
[9570]49   USE agrif_oce_interp ! agrif
[9124]50   USE agrif_oce
[9019]51#endif
52#if defined key_asminc   
53   USE asminc          ! Assimilation increment
54#endif
[6140]55   !
56   USE in_out_manager  ! I/O manager
[358]57   USE lib_mpp         ! distributed memory computing library
58   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
59   USE prtctl          ! Print control
[2715]60   USE iom             ! IOM library
[4292]61   USE restart         ! only for lrst_oce
[358]62
[11536]63   USE iom   ! to remove
64
[358]65   IMPLICIT NONE
66   PRIVATE
67
[9124]68   PUBLIC dyn_spg_ts        ! called by dyn_spg
69   PUBLIC dyn_spg_ts_init   !    -    - dyn_spg_init
[358]70
[9019]71   !! Time filtered arrays at baroclinic time step:
72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
[9124]73   !
[12443]74   INTEGER, SAVE :: icycle      ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
75   REAL(wp),SAVE :: rDt_e       ! Barotropic time step
[9019]76   !
77   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
[9124]78   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwz                ! ff_f/h at F points
79   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftnw, ftne         ! triad of coriolis parameter
80   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ftsw, ftse         ! (only used with een vorticity scheme)
[4292]81
[9043]82   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
83   REAL(wp) ::   r1_8  = 0.125_wp         !
84   REAL(wp) ::   r1_4  = 0.25_wp          !
85   REAL(wp) ::   r1_2  = 0.5_wp           !
[508]86
[358]87   !! * Substitutions
[12377]88#  include "do_loop_substitute.h90"
[2715]89   !!----------------------------------------------------------------------
[9598]90   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]91   !! $Id$
[10068]92   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]93   !!----------------------------------------------------------------------
[358]94CONTAINS
95
[2715]96   INTEGER FUNCTION dyn_spg_ts_alloc()
97      !!----------------------------------------------------------------------
98      !!                  ***  routine dyn_spg_ts_alloc  ***
99      !!----------------------------------------------------------------------
[6140]100      INTEGER :: ierr(3)
[4292]101      !!----------------------------------------------------------------------
102      ierr(:) = 0
[6140]103      !
[12443]104      ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) )
[9528]105      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
[11536]106         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2)   )
[6140]107         !
108      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
109      !
110      dyn_spg_ts_alloc = MAXVAL( ierr(:) )
111      !
[10425]112      CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
113      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
[2715]114      !
115   END FUNCTION dyn_spg_ts_alloc
116
[5836]117
[12377]118   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
[358]119      !!----------------------------------------------------------------------
120      !!
[6140]121      !! ** Purpose : - Compute the now trend due to the explicit time stepping
122      !!              of the quasi-linear barotropic system, and add it to the
123      !!              general momentum trend.
[358]124      !!
[6140]125      !! ** Method  : - split-explicit schem (time splitting) :
[4374]126      !!      Barotropic variables are advanced from internal time steps
127      !!      "n"   to "n+1" if ln_bt_fw=T
128      !!      or from
129      !!      "n-1" to "n+1" if ln_bt_fw=F
130      !!      thanks to a generalized forward-backward time stepping (see ref. below).
[358]131      !!
[4374]132      !! ** Action :
[12377]133      !!      -Update the filtered free surface at step "n+1"      : pssh(:,:,Kaa)
134      !!      -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa)
[9023]135      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv
[4374]136      !!      These are used to advect tracers and are compliant with discrete
137      !!      continuity equation taken at the baroclinic time steps. This
138      !!      ensures tracers conservation.
[12377]139      !!      - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
[358]140      !!
[6140]141      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
[358]142      !!---------------------------------------------------------------------
[12377]143      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
144      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
145      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
146      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
[2715]147      !
[9554]148      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
[9019]149      LOGICAL  ::   ll_fw_start           ! =T : forward integration
[9554]150      LOGICAL  ::   ll_init               ! =T : special startup of 2d equations
[11536]151      INTEGER  ::   noffset               ! local integers  : time offset for bdy update
[12424]152      REAL(wp) ::   r1_Dt_b, z1_hu, z1_hv          ! local scalars
[9528]153      REAL(wp) ::   za0, za1, za2, za3              !   -      -
[12206]154      REAL(wp) ::   zztmp, zldg               !   -      -
[11536]155      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      -
156      REAL(wp) ::   zun_save, zvn_save              !   -      -
157      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
158      REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
159      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
160      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
[9019]161      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
[11536]162      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
[3294]163      !
[9023]164      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.
165
166      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point
167
168      REAL(wp) ::   zepsilon, zgamma            !   -      -
[9019]169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef.
[9023]170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask
[12377]172      REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep
[358]173      !!----------------------------------------------------------------------
[3294]174      !
[9023]175      IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
176      !                                         !* Allocate temporary arrays
177      IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
[3294]178      !
[9023]179      zwdramp = r_rn_wdmin1               ! simplest ramp
180!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
[11536]181      !                                         ! inverse of baroclinic time step
[12424]182      r1_Dt_b = 1._wp / rDt 
[4292]183      !
[9023]184      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart
[4292]185      ll_fw_start = .FALSE.
[9023]186      !                                         ! time offset in steps for bdy data update
[12443]187      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_e
[6140]188      ELSE                       ;   noffset =   0 
189      ENDIF
[4292]190      !
[9023]191      IF( kt == nit000 ) THEN                   !* initialisation
[508]192         !
[358]193         IF(lwp) WRITE(numout,*)
194         IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
195         IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
[4354]196         IF(lwp) WRITE(numout,*)
[1502]197         !
[12397]198         IF( l_1st_euler )   ll_init=.TRUE.
[1502]199         !
[12397]200         IF( ln_bt_fw .OR. l_1st_euler ) THEN
[6140]201            ll_fw_start =.TRUE.
202            noffset     = 0
[4292]203         ELSE
[6140]204            ll_fw_start =.FALSE.
[4292]205         ENDIF
[11536]206         !                    ! Set averaging weights and cycle length:
[6140]207         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]208         !
[12397]209      ELSEIF( kt == nit000 + 1 ) THEN           !* initialisation 2nd time-step
210         !
211         IF( .NOT.ln_bt_fw ) THEN
212            ! If we did an Euler timestep on the first timestep we need to reset ll_fw_start
213            ! and the averaging weights. We don't have an easy way of telling whether we did
214            ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false.
215            ! at the end of the first timestep) so just do this in all cases.
216            ll_fw_start = .FALSE.
217            CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
218         ENDIF
219         !
[4292]220      ENDIF
221      !
[358]222      ! -----------------------------------------------------------------------------
223      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
224      ! -----------------------------------------------------------------------------
[1502]225      !     
[4292]226      !
[11536]227      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
228      !                                   !  ---------------------------  !
[12377]229      zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
230      zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
[1502]231      !
[4292]232      !
[12377]233      !                                   !=  U(Krhs) => baroclinic trend  =!   (remove its vertical mean)
234      DO jk = 1, jpkm1                    !  -----------------------------  !
235         uu(:,:,jk,Krhs) = ( uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk)
236         vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)
[1502]237      END DO
[7646]238     
239!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
240!!gm  Is it correct to do so ?   I think so...
241     
[11536]242      !                                   !=  remove 2D Coriolis and pressure gradient trends  =!
243      !                                   !  -------------------------------------------------  !
[9528]244      !
[12377]245      IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init( Kmm )   ! Set zwz, the barotropic Coriolis force coefficient
[11536]246      !       ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
[1502]247      !
[11536]248      !                                         !* 2D Coriolis trends
[12377]249      zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:)        ! now fluxes
250      zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
[11536]251      !
[12377]252      CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in
253         &                                                                     zu_trd, zv_trd   )   ! ==>> out
[11536]254      !
255      IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only)
[508]256         !
[11536]257         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
[12377]258            CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
259            DO_2D_00_00
260               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj ,Kmm) )   &
261                  &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
262               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj ,Kmm) )   &
263                  &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
264            END_2D
[11536]265         ELSE                                      ! now suface pressure gradient
[12377]266            DO_2D_00_00
267               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e1u(ji,jj)
268               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj) 
269            END_2D
[9019]270         ENDIF
271         !
[1502]272      ENDIF
[9019]273      !
[12377]274      DO_2D_00_00
275          zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
276          zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
277      END_2D
[4292]278      !
[11536]279      !                                   !=  Add bottom stress contribution from baroclinic velocities  =!
280      !                                   !  -----------------------------------------------------------  !
[12377]281      CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients
[11536]282      !                                   !=  Add atmospheric pressure forcing  =!
283      !                                   !  ----------------------------------  !
284      IF( ln_apr_dyn ) THEN
285         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
[12377]286            DO_2D_00_00
287               zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
288               zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
289            END_2D
[11536]290         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
291            zztmp = grav * r1_2
[12377]292            DO_2D_00_00
293               zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
294                    &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
295               zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  &
296                    &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
297            END_2D
298         ENDIF
[9019]299      ENDIF
[11536]300      !
301      !                                   !=  Add atmospheric pressure forcing  =!
302      !                                   !  ----------------------------------  !
[9019]303      IF( ln_bt_fw ) THEN                        ! Add wind forcing
[12377]304         DO_2D_00_00
[12443]305            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
306            zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
[12377]307         END_2D
[2724]308      ELSE
[12443]309         zztmp = r1_rho0 * r1_2
[12377]310         DO_2D_00_00
311            zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
312            zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
313         END_2D
[4292]314      ENDIF 
315      !
[11536]316      !              !----------------!
317      !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain)
318      !              !----------------!
319      !                                   !=  Net water flux forcing applied to a water column  =!
320      !                                   ! ---------------------------------------------------  !
321      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
[12443]322         zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
[11536]323      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
[12443]324         zztmp = r1_rho0 * r1_2
[12377]325         zssh_frc(:,:) = zztmp * (  emp(:,:)        + emp_b(:,:)                    &
326                &                 - rnf(:,:)        - rnf_b(:,:)                    &
327                &                 + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)             &
328                &                 + fwfisf_par(:,:) + fwfisf_par_b(:,:)             )
[4292]329      ENDIF
[11536]330      !                                   !=  Add Stokes drift divergence  =!   (if exist)
331      IF( ln_sdw ) THEN                   !  -----------------------------  !
[7753]332         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
[7646]333      ENDIF
334      !
[12377]335      !                                         ! ice sheet coupling
336      IF ( ln_isf .AND. ln_isfcpl ) THEN
337         !
338         ! ice sheet coupling
339         IF( ln_rstart .AND. kt == nit000 ) THEN
340            zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:)
341         END IF
342         !
343         ! conservation option
344         IF( ln_isfcpl_cons ) THEN
345            zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:)
346         END IF
347         !
348      END IF
349      !
[4292]350#if defined key_asminc
[11536]351      !                                   !=  Add the IAU weighted SSH increment  =!
352      !                                   !  ------------------------------------  !
[4292]353      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[7753]354         zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
[4292]355      ENDIF
356#endif
[11536]357      !                                   != Fill boundary data arrays for AGRIF
[5656]358      !                                   ! ------------------------------------
[4486]359#if defined key_agrif
360         IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
361#endif
[4292]362      !
[358]363      ! -----------------------------------------------------------------------
[4292]364      !  Phase 2 : Integration of the barotropic equations
[358]365      ! -----------------------------------------------------------------------
[1502]366      !
367      !                                             ! ==================== !
368      !                                             !    Initialisations   !
[4292]369      !                                             ! ==================== ! 
[4370]370      ! Initialize barotropic variables:     
[4770]371      IF( ll_init )THEN
[7753]372         sshbb_e(:,:) = 0._wp
373         ubb_e  (:,:) = 0._wp
374         vbb_e  (:,:) = 0._wp
375         sshb_e (:,:) = 0._wp
376         ub_e   (:,:) = 0._wp
377         vb_e   (:,:) = 0._wp
[4700]378      ENDIF
379      !
[11536]380      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
[12377]381         zhup2_e(:,:) = hu(:,:,Kmm)
382         zhvp2_e(:,:) = hv(:,:,Kmm)
383         zhtp2_e(:,:) = ht(:,:)
[11536]384      ENDIF
385      !
[4370]386      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                   
[12377]387         sshn_e(:,:) =    pssh(:,:,Kmm)           
388         un_e  (:,:) =    puu_b(:,:,Kmm)           
389         vn_e  (:,:) =    pvv_b(:,:,Kmm)
[7753]390         !
[12377]391         hu_e  (:,:) =    hu(:,:,Kmm)       
392         hv_e  (:,:) =    hv(:,:,Kmm) 
393         hur_e (:,:) = r1_hu(:,:,Kmm)   
394         hvr_e (:,:) = r1_hv(:,:,Kmm)
[4370]395      ELSE                                ! CENTRED integration: start from BEFORE fields
[12377]396         sshn_e(:,:) =    pssh(:,:,Kbb)
397         un_e  (:,:) =    puu_b(:,:,Kbb)         
398         vn_e  (:,:) =    pvv_b(:,:,Kbb)
[7753]399         !
[12377]400         hu_e  (:,:) =    hu(:,:,Kbb)       
401         hv_e  (:,:) =    hv(:,:,Kbb) 
402         hur_e (:,:) = r1_hu(:,:,Kbb)   
403         hvr_e (:,:) = r1_hv(:,:,Kbb)
[4292]404      ENDIF
405      !
406      ! Initialize sums:
[12377]407      puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)         
408      pvv_b  (:,:,Kaa) = 0._wp
409      pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level
[7753]410      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
411      vn_adv(:,:) = 0._wp
[9528]412      !
413      IF( ln_wd_dl ) THEN
[9023]414         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)
415         zvwdmask(:,:) = 0._wp  !
[9528]416         zuwdav2 (:,:) = 0._wp 
417         zvwdav2 (:,:) = 0._wp   
[9023]418      END IF 
419
[9528]420      !                                             ! ==================== !
[4292]421      DO jn = 1, icycle                             !  sub-time-step loop  !
[1502]422         !                                          ! ==================== !
[10425]423         !
424         l_full_nf_update = jn == icycle   ! false: disable full North fold update (performances) for jn = 1 to icycle-1
[4292]425         !
[11536]426         !                    !==  Update the forcing ==! (BDY and tides)
427         !
[12377]428         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) )
429         ! Update tide potential at the beginning of current time substep
430         IF( ln_tide_pot .AND. ln_tide ) THEN
[12443]431            zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_e, wp)
[12377]432            CALL upd_tide(zt0substep, Kmm)
433         END IF
[11536]434         !
435         !                    !==  extrapolation at mid-step  ==!   (jn+1/2)
436         !
437         !                       !* Set extrapolation coefficients for predictor step:
[4292]438         IF ((jn<3).AND.ll_init) THEN      ! Forward           
439           za1 = 1._wp                                         
440           za2 = 0._wp                       
441           za3 = 0._wp                       
442         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
443           za1 =  1.781105_wp              ! za1 =   3/2 +   bet
444           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
445           za3 =  0.281105_wp              ! za3 = bet
446         ENDIF
[11536]447         !
448         !                       !* Extrapolate barotropic velocities at mid-step (jn+1/2)
449         !--        m+1/2               m                m-1           m-2       --!
450         !--       u      = (3/2+beta) u   -(1/2+2beta) u      + beta u          --!
451         !-------------------------------------------------------------------------!
[7753]452         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
453         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
[4292]454
[6140]455         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only)
[4292]456            !                                             !  ------------------
457            ! Extrapolate Sea Level at step jit+0.5:
[11536]458            !--         m+1/2                 m                  m-1             m-2       --!
459            !--      ssh      = (3/2+beta) ssh   -(1/2+2beta) ssh      + beta ssh          --!
460            !--------------------------------------------------------------------------------!
[7753]461            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[9023]462           
[11536]463            ! set wetting & drying mask at tracer points for this barotropic mid-step
464            IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask )
[9528]465            !
[11536]466            !                          ! ocean t-depth at mid-step
467            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
468            !
469            !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
[12377]470            DO_2D_11_10
471               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
472                    &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
473                    &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
474            END_2D
475            DO_2D_10_11
476               zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
477                    &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
478                    &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
479            END_2D
[4292]480            !
481         ENDIF
482         !
[11536]483         !                    !==  after SSH  ==!   (jn+1)
484         !
485         !                             ! update (ua_e,va_e) to enforce volume conservation at open boundaries
486         !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
[10481]487         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
[12377]488         !     
[11536]489         !                             ! resulting flux at mid-step (not over the full domain)
490         zhU(1:jpim1,1:jpj  ) = e2u(1:jpim1,1:jpj  ) * ua_e(1:jpim1,1:jpj  ) * zhup2_e(1:jpim1,1:jpj  )   ! not jpi-column
491         zhV(1:jpi  ,1:jpjm1) = e1v(1:jpi  ,1:jpjm1) * va_e(1:jpi  ,1:jpjm1) * zhvp2_e(1:jpi  ,1:jpjm1)   ! not jpj-row
[4486]492         !
493#if defined key_agrif
[6140]494         ! Set fluxes during predictor step to ensure volume conservation
[12377]495         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
[4486]496#endif
[12443]497         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e)    !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
[9023]498
[11536]499         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where
500            !                          ! the direction of the flow is from dry cells
501            CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )   ! not jpi colomn for U, not jpj row for V
[9528]502            !
503         ENDIF   
[11536]504         !
505         !
506         !     Compute Sea Level at step jit+1
507         !--           m+1        m                               m+1/2          --!
508         !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
509         !-------------------------------------------------------------------------!
[12377]510         DO_2D_00_00
511            zhdiv = (   zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
[12443]512            ssha_e(ji,jj) = (  sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv )  ) * ssmask(ji,jj)
[12377]513         END_2D
[11536]514         !
515         CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp,  zhU, 'U', -1._wp,  zhV, 'V', -1._wp )
516         !
517         !                             ! Sum over sub-time-steps to compute advective velocities
518         za2 = wgtbtp2(jn)             ! zhU, zhV hold fluxes extrapolated at jn+0.5
519         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
520         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
521         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
522         IF ( ln_wd_dl_bc ) THEN
523            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
524            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
525         END IF
526         !
[6140]527         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
[7646]528         IF( ln_bdy )   CALL bdy_ssh( ssha_e )
[4292]529#if defined key_agrif
[6140]530         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn )
[4292]531#endif
532         
533         ! Sea Surface Height at u-,v-points (vvl case only)
[6140]534         IF( .NOT.ln_linssh ) THEN                               
[12377]535            DO_2D_00_00
536               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
537                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
538                  &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
539               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
540                  &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
541                  &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
542            END_2D
[4292]543         ENDIF   
[11536]544         !         
545         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
546         !--            m+1/2           m+1              m               m-1              m-2     --!
547         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --!
548         !------------------------------------------------------------------------------------------!
549         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation
[9528]550         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
551            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
[1502]552         !
[11536]553         !                             ! Surface pressure gradient
[12443]554         zldg = ( 1._wp - rn_load ) * grav    ! local factor
[12377]555         DO_2D_00_00
556            zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
557            zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
558         END_2D
[11536]559         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient
560            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters
561            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
562            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
[4292]563         ENDIF
564         !
565         ! Add Coriolis trend:
[6140]566         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]567         ! at each time step. We however keep them constant here for optimization.
[11536]568         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
569         CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
[4292]570         !
571         ! Add tidal astronomical forcing if defined
[7646]572         IF ( ln_tide .AND. ln_tide_pot ) THEN
[12377]573            DO_2D_00_00
574               zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
575               zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
576            END_2D
[4292]577         ENDIF
578         !
[9023]579         ! Add bottom stresses:
580!jth do implicitly instead
581         IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[12377]582            DO_2D_00_00
583               zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
584               zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
585            END_2D
[11536]586         ENDIF
[4292]587         !
588         ! Set next velocities:
[11536]589         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
590         !--                              VECTOR FORM
591         !--   m+1                 m               /                                                       m+1/2           \    --!
592         !--  u     =             u   + delta_t' * \         (1-r)*g * grad_x( ssh') -         f * k vect u      +     frc /    --!
593         !--                                                                                                                    --!
594         !--                             FLUX FORM                                                                              --!
595         !--  m+1   __1__  /  m    m               /  m+1/2                             m+1/2              m+1/2    n      \ \  --!
596         !-- u    =   m+1 |  h  * u   + delta_t' * \ h     * (1-r)*g * grad_x( ssh') - h     * f * k vect u      + h * frc /  | --!
597         !--         h     \                                                                                                 /  --!
598         !------------------------------------------------------------------------------------------------------------------------!
[9023]599         IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
[12377]600            DO_2D_00_00
601               ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
[12443]602                         &     + rDt_e * (                   zu_spg(ji,jj)   &
[12377]603                         &                                 + zu_trd(ji,jj)   &
604                         &                                 + zu_frc(ji,jj) ) & 
605                         &   ) * ssumask(ji,jj)
[358]606
[12377]607               va_e(ji,jj) = (                                 vn_e(ji,jj)   &
[12443]608                         &     + rDt_e * (                   zv_spg(ji,jj)   &
[12377]609                         &                                 + zv_trd(ji,jj)   &
610                         &                                 + zv_frc(ji,jj) ) &
611                         &   ) * ssvmask(ji,jj)
612            END_2D
[6140]613            !
[9023]614         ELSE                           !* Flux form
[12377]615            DO_2D_00_00
616               !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
617               !                    ! backward interpolated depth used in spg terms at jn+1/2
618               zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
619                    &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
620               zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
621                    &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
622               !                    ! inverse depth at jn+1
623               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
624               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
625               !
626               ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
[12443]627                    &            + rDt_e * (  zhu_bck        * zu_spg (ji,jj)  &   !
[12377]628                    &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
629                    &                       +  hu(ji,jj,Kmm) * zu_frc (ji,jj)  )   ) * z1_hu
630               !
631               va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      &
[12443]632                    &            + rDt_e * (  zhv_bck        * zv_spg (ji,jj)  &   !
[12377]633                    &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
634                    &                       +  hv(ji,jj,Kmm) * zv_frc (ji,jj)  )   ) * z1_hv
635            END_2D
[4292]636         ENDIF
[10272]637!jth implicit bottom friction:
638         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
[12377]639            DO_2D_00_00
[12443]640                  ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
641                  va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
[12377]642            END_2D
[10272]643         ENDIF
[11536]644       
645         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only)
646            hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
647            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
648            hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
649            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
650            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  &
[12072]651                 &                         , hu_e , 'U',  1._wp, hv_e , 'V',  1._wp  &
652                 &                         , hur_e, 'U',  1._wp, hvr_e, 'V',  1._wp  )
[11536]653         ELSE
654            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  )
[1438]655         ENDIF
[4292]656         !
[11536]657         !
[6140]658         !                                                 ! open boundaries
[7646]659         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
[4486]660#if defined key_agrif                                                           
661         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
[4292]662#endif
663         !                                             !* Swap
664         !                                             !  ----
[7753]665         ubb_e  (:,:) = ub_e  (:,:)
666         ub_e   (:,:) = un_e  (:,:)
667         un_e   (:,:) = ua_e  (:,:)
668         !
669         vbb_e  (:,:) = vb_e  (:,:)
670         vb_e   (:,:) = vn_e  (:,:)
671         vn_e   (:,:) = va_e  (:,:)
672         !
673         sshbb_e(:,:) = sshb_e(:,:)
674         sshb_e (:,:) = sshn_e(:,:)
675         sshn_e (:,:) = ssha_e(:,:)
[4292]676
677         !                                             !* Sum over whole bt loop
678         !                                             !  ----------------------
679         za1 = wgtbtp1(jn)                                   
[6140]680         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities
[12377]681            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) 
682            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) 
[9023]683         ELSE                                       ! Sum transports
684            IF ( .NOT.ln_wd_dl ) THEN 
[12377]685               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:)
686               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:)
[9023]687            ELSE
[12377]688               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:)
689               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:)
[9023]690            END IF
[4292]691         ENDIF
[9023]692         !                                          ! Sum sea level
[12377]693         pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
[9023]694
[358]695         !                                                 ! ==================== !
696      END DO                                               !        end loop      !
697      !                                                    ! ==================== !
[1438]698      ! -----------------------------------------------------------------------------
[1502]699      ! Phase 3. update the general trend with the barotropic trend
[1438]700      ! -----------------------------------------------------------------------------
[1502]701      !
[4292]702      ! Set advection velocity correction:
[9023]703      IF (ln_bt_fw) THEN
[12397]704         IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
[12377]705            DO_2D_11_11
706               zun_save = un_adv(ji,jj)
707               zvn_save = vn_adv(ji,jj)
708               !                          ! apply the previously computed correction
[12443]709               un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
710               vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
[12377]711               !                          ! Update corrective fluxes for next time step
[12443]712               un_bf(ji,jj)  = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
713               vn_bf(ji,jj)  = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
[12377]714               !                          ! Save integrated transport for next computation
715               ub2_b(ji,jj) = zun_save
716               vb2_b(ji,jj) = zvn_save
717            END_2D
[9023]718         ELSE
[11536]719            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
720            vn_bf(:,:) = 0._wp
721            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
722            vb2_b(:,:) = vn_adv(:,:)
723         END IF
[4292]724      ENDIF
[9023]725
726
[4292]727      !
728      ! Update barotropic trend:
[6140]729      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]730         DO jk=1,jpkm1
[12424]731            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
732            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
[4292]733         END DO
734      ELSE
[12377]735         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
736         DO_2D_10_10
737            zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) &
738               &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      &
739               &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
740            zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) &
741               &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      &
742               &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
743         END_2D
[10425]744         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
[5930]745         !
[4292]746         DO jk=1,jpkm1
[12424]747            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
748            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
[4292]749         END DO
750         ! Save barotropic velocities not transport:
[12377]751         puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
752         pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]753      ENDIF
[9023]754
755
756      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 
[4292]757      DO jk = 1, jpkm1
[12377]758         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk)
759         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
[358]760      END DO
[9023]761
762      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
763         DO jk = 1, jpkm1
[12377]764            puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
765                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 
766            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 
767                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 
[9023]768         END DO
769      END IF
770
771     
[12377]772      CALL iom_put(  "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) )    ! barotropic i-current
773      CALL iom_put(  "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) )    ! barotropic i-current
[1502]774      !
[4486]775#if defined key_agrif
776      ! Save time integrated fluxes during child grid integration
[5656]777      ! (used to update coarse grid transports at next time step)
[4486]778      !
[6140]779      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
780         IF( Agrif_NbStepint() == 0 ) THEN
[7753]781            ub2_i_b(:,:) = 0._wp
782            vb2_i_b(:,:) = 0._wp
[4486]783         END IF
784         !
785         za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]786         ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
787         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]788      ENDIF
789#endif     
[1502]790      !                                   !* write time-spliting arrays in the restart
[6140]791      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
[508]792      !
[9023]793      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
794      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]795      !
[12377]796      CALL iom_put( "baro_u" , puu_b(:,:,Kmm) )  ! Barotropic  U Velocity
797      CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) )  ! Barotropic  V Velocity
[2715]798      !
[508]799   END SUBROUTINE dyn_spg_ts
800
[6140]801
[4292]802   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
803      !!---------------------------------------------------------------------
804      !!                   ***  ROUTINE ts_wgt  ***
805      !!
806      !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
807      !!----------------------------------------------------------------------
808      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
809      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
810      INTEGER, INTENT(inout) :: jpit      ! cycle length   
[12443]811      REAL(wp), DIMENSION(3*nn_e), INTENT(inout) ::   zwgt1, & ! Primary weights
[4292]812                                                         zwgt2    ! Secondary weights
813     
814      INTEGER ::  jic, jn, ji                      ! temporary integers
815      REAL(wp) :: za1, za2
816      !!----------------------------------------------------------------------
[508]817
[4292]818      zwgt1(:) = 0._wp
819      zwgt2(:) = 0._wp
820
821      ! Set time index when averaged value is requested
822      IF (ll_fw) THEN
[12443]823         jic = nn_e
[4292]824      ELSE
[12443]825         jic = 2 * nn_e
[4292]826      ENDIF
827
828      ! Set primary weights:
829      IF (ll_av) THEN
830           ! Define simple boxcar window for primary weights
[12443]831           ! (width = nn_e, centered around jic)     
[4292]832         SELECT CASE ( nn_bt_flt )
833              CASE( 0 )  ! No averaging
834                 zwgt1(jic) = 1._wp
835                 jpit = jic
836
[12443]837              CASE( 1 )  ! Boxcar, width = nn_e
838                 DO jn = 1, 3*nn_e
839                    za1 = ABS(float(jn-jic))/float(nn_e) 
[4292]840                    IF (za1 < 0.5_wp) THEN
841                      zwgt1(jn) = 1._wp
842                      jpit = jn
843                    ENDIF
844                 ENDDO
845
[12443]846              CASE( 2 )  ! Boxcar, width = 2 * nn_e
847                 DO jn = 1, 3*nn_e
848                    za1 = ABS(float(jn-jic))/float(nn_e) 
[4292]849                    IF (za1 < 1._wp) THEN
850                      zwgt1(jn) = 1._wp
851                      jpit = jn
852                    ENDIF
853                 ENDDO
854              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
855         END SELECT
856
857      ELSE ! No time averaging
858         zwgt1(jic) = 1._wp
859         jpit = jic
860      ENDIF
861   
862      ! Set secondary weights
863      DO jn = 1, jpit
864        DO ji = jn, jpit
865             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
866        END DO
867      END DO
868
869      ! Normalize weigths:
870      za1 = 1._wp / SUM(zwgt1(1:jpit))
871      za2 = 1._wp / SUM(zwgt2(1:jpit))
872      DO jn = 1, jpit
873        zwgt1(jn) = zwgt1(jn) * za1
874        zwgt2(jn) = zwgt2(jn) * za2
875      END DO
876      !
877   END SUBROUTINE ts_wgt
878
[6140]879
[508]880   SUBROUTINE ts_rst( kt, cdrw )
881      !!---------------------------------------------------------------------
882      !!                   ***  ROUTINE ts_rst  ***
883      !!
884      !! ** Purpose : Read or write time-splitting arrays in restart file
885      !!----------------------------------------------------------------------
[9528]886      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
887      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
[508]888      !!----------------------------------------------------------------------
889      !
[9506]890      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
891         !                                   ! ---------------
[12397]892         IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN    !* Read the restart file
[9506]893            CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:), ldxios = lrxios )   
894            CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:), ldxios = lrxios ) 
895            CALL iom_get( numror, jpdom_autoglo, 'un_bf'  , un_bf  (:,:), ldxios = lrxios )   
896            CALL iom_get( numror, jpdom_autoglo, 'vn_bf'  , vn_bf  (:,:), ldxios = lrxios ) 
897            IF( .NOT.ln_bt_av ) THEN
898               CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:), ldxios = lrxios )   
899               CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:), ldxios = lrxios )   
900               CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:), ldxios = lrxios )
901               CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:), ldxios = lrxios ) 
902               CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:), ldxios = lrxios )   
903               CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:), ldxios = lrxios )
904            ENDIF
[4486]905#if defined key_agrif
[9506]906            ! Read time integrated fluxes
907            IF ( .NOT.Agrif_Root() ) THEN
908               CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lrxios )   
909               CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lrxios )
910            ENDIF
911#endif
912         ELSE                                   !* Start from rest
913            IF(lwp) WRITE(numout,*)
914            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
915            ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif
916            un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif
917            un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif
918#if defined key_agrif
919            IF ( .NOT.Agrif_Root() ) THEN
920               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
921            ENDIF
922#endif
[4486]923         ENDIF
[9506]924         !
925      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
926         !                                   ! -------------------
927         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]928         IF( lwxios ) CALL iom_swap(      cwxios_context          )
929         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:), ldxios = lwxios )
930         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:), ldxios = lwxios )
931         CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:), ldxios = lwxios )
932         CALL iom_rstput( kt, nitrst, numrow, 'vn_bf'   , vn_bf  (:,:), ldxios = lwxios )
[4292]933         !
934         IF (.NOT.ln_bt_av) THEN
[9367]935            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:), ldxios = lwxios ) 
936            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:), ldxios = lwxios )
937            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:), ldxios = lwxios )
938            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:), ldxios = lwxios )
939            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:), ldxios = lwxios )
940            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:), ldxios = lwxios )
[4292]941         ENDIF
[4486]942#if defined key_agrif
943         ! Save time integrated fluxes
944         IF ( .NOT.Agrif_Root() ) THEN
[9367]945            CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b'  , ub2_i_b(:,:), ldxios = lwxios )
946            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:), ldxios = lwxios )
[4486]947         ENDIF
948#endif
[9367]949         IF( lwxios ) CALL iom_swap(      cxios_context          )
[4292]950      ENDIF
951      !
952   END SUBROUTINE ts_rst
[2528]953
[6140]954
955   SUBROUTINE dyn_spg_ts_init
[4292]956      !!---------------------------------------------------------------------
957      !!                   ***  ROUTINE dyn_spg_ts_init  ***
958      !!
959      !! ** Purpose : Set time splitting options
960      !!----------------------------------------------------------------------
[6140]961      INTEGER  ::   ji ,jj              ! dummy loop indices
962      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
[9019]963      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
[4292]964      !!----------------------------------------------------------------------
[4370]965      !
[5930]966      ! Max courant number for ext. grav. waves
[4370]967      !
[12377]968      DO_2D_11_11
969         zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
970         zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
971         zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
972      END_2D
[5930]973      !
[5836]974      zcmax = MAXVAL( zcu(:,:) )
[10425]975      CALL mpp_max( 'dynspg_ts', zcmax )
[2528]976
[4370]977      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[12443]978      IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
[4292]979     
[12443]980      rDt_e = rn_Dt / REAL( nn_e , wp )
981      zcmax = zcmax * rDt_e
[9023]982      ! Print results
[4292]983      IF(lwp) WRITE(numout,*)
[9169]984      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
985      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]986      IF( ln_bt_auto ) THEN
[12443]987         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_e '
[4370]988         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
[4292]989      ELSE
[12443]990         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
[358]991      ENDIF
[4292]992
993      IF(ln_bt_av) THEN
[12443]994         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
[4292]995      ELSE
[9169]996         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]997      ENDIF
[508]998      !
[4292]999      !
1000      IF(ln_bt_fw) THEN
[4370]1001         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
[4292]1002      ELSE
[4370]1003         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]1004      ENDIF
1005      !
[4486]1006#if defined key_agrif
1007      ! Restrict the use of Agrif to the forward case only
[9023]1008!!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]1009#endif
1010      !
[4370]1011      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]1012      SELECT CASE ( nn_bt_flt )
[6140]1013         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
[12443]1014         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
1015         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
[9169]1016         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]1017      END SELECT
1018      !
[4370]1019      IF(lwp) WRITE(numout,*) ' '
[12443]1020      IF(lwp) WRITE(numout,*) '     nn_e = ', nn_e
1021      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rDt_e
[4370]1022      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
1023      !
[9023]1024      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1025      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1026         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1027      ENDIF
1028      !
[6140]1029      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]1030         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1031      ENDIF
[6140]1032      IF( zcmax>0.9_wp ) THEN
[12443]1033         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )         
[4292]1034      ENDIF
1035      !
[9124]1036      !                             ! Allocate time-splitting arrays
1037      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
1038      !
1039      !                             ! read restart when needed
[9506]1040      CALL ts_rst( nit000, 'READ' )
[9124]1041      !
[9367]1042      IF( lwxios ) THEN
1043! define variables in restart file when writing with XIOS
1044         CALL iom_set_rstw_var_active('ub2_b')
1045         CALL iom_set_rstw_var_active('vb2_b')
1046         CALL iom_set_rstw_var_active('un_bf')
1047         CALL iom_set_rstw_var_active('vn_bf')
1048         !
1049         IF (.NOT.ln_bt_av) THEN
1050            CALL iom_set_rstw_var_active('sshbb_e')
1051            CALL iom_set_rstw_var_active('ubb_e')
1052            CALL iom_set_rstw_var_active('vbb_e')
1053            CALL iom_set_rstw_var_active('sshb_e')
1054            CALL iom_set_rstw_var_active('ub_e')
1055            CALL iom_set_rstw_var_active('vb_e')
1056         ENDIF
1057#if defined key_agrif
1058         ! Save time integrated fluxes
1059         IF ( .NOT.Agrif_Root() ) THEN
1060            CALL iom_set_rstw_var_active('ub2_i_b')
1061            CALL iom_set_rstw_var_active('vb2_i_b')
1062         ENDIF
1063#endif
1064      ENDIF
1065      !
[4292]1066   END SUBROUTINE dyn_spg_ts_init
[508]1067
[11536]1068   
[12377]1069   SUBROUTINE dyn_cor_2D_init( Kmm )
[11536]1070      !!---------------------------------------------------------------------
[12377]1071      !!                   ***  ROUTINE dyn_cor_2D_init  ***
[11536]1072      !!
1073      !! ** Purpose : Set time splitting options
1074      !! Set arrays to remove/compute coriolis trend.
1075      !! Do it once during initialization if volume is fixed, else at each long time step.
1076      !! Note that these arrays are also used during barotropic loop. These are however frozen
1077      !! although they should be updated in the variable volume case. Not a big approximation.
1078      !! To remove this approximation, copy lines below inside barotropic loop
1079      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1080      !!
1081      !! Compute zwz = f / ( height of the water colomn )
1082      !!----------------------------------------------------------------------
[12377]1083      INTEGER,  INTENT(in)         ::  Kmm  ! Time index
[11536]1084      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
1085      REAL(wp) ::   z1_ht
1086      REAL(wp), DIMENSION(jpi,jpj) :: zhf
1087      !!----------------------------------------------------------------------
1088      !
1089      SELECT CASE( nvor_scheme )
1090      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
1091         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
1092         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
[12377]1093            DO_2D_10_10
1094               zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
1095                    &           ht(ji  ,jj  ) + ht(ji+1,jj  )   ) * 0.25_wp 
1096               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1097            END_2D
[11536]1098         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
[12377]1099            DO_2D_10_10
1100               zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      &
1101                    &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   &
1102                    &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
1103                    &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
1104               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1105            END_2D
[11536]1106         END SELECT
1107         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
1108         !
1109         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[12377]1110         DO_2D_01_01
1111            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
1112            ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
1113            ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
1114            ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
1115         END_2D
[11536]1116         !
1117      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
1118         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[12377]1119         DO_2D_01_01
1120            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
1121            ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
1122            ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
1123            ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1124            ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
1125         END_2D
[11536]1126         !
1127      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
1128         !
1129         zwz(:,:) = 0._wp
1130         zhf(:,:) = 0._wp
1131         
1132         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1133!!gm    A priori a better value should be something like :
1134!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1135!!gm                     divided by the sum of the corresponding mask
1136!!gm
1137!!           
1138         IF( .NOT.ln_sco ) THEN
1139 
1140   !!gm  agree the JC comment  : this should be done in a much clear way
1141 
1142   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1143   !     Set it to zero for the time being
1144   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
1145   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
1146   !              ENDIF
1147   !              zhf(:,:) = gdepw_0(:,:,jk+1)
1148            !
1149         ELSE
1150            !
1151            !zhf(:,:) = hbatf(:,:)
[12377]1152            DO_2D_10_10
1153               zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
1154                    &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
1155                    &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
1156                    &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
1157            END_2D
[11536]1158         ENDIF
1159         !
1160         DO jj = 1, jpjm1
1161            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1162         END DO
1163         !
1164         DO jk = 1, jpkm1
1165            DO jj = 1, jpjm1
[12377]1166               zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[11536]1167            END DO
1168         END DO
1169         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
1170         ! JC: TBC. hf should be greater than 0
[12377]1171         DO_2D_11_11
1172            IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
1173         END_2D
[11536]1174         zwz(:,:) = ff_f(:,:) * zwz(:,:)
1175      END SELECT
1176     
1177   END SUBROUTINE dyn_cor_2d_init
1178
1179
1180
[12377]1181   SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
[11536]1182      !!---------------------------------------------------------------------
1183      !!                   ***  ROUTINE dyn_cor_2d  ***
1184      !!
1185      !! ** Purpose : Compute u and v coriolis trends
1186      !!----------------------------------------------------------------------
1187      INTEGER  ::   ji ,jj                             ! dummy loop indices
1188      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
[12377]1189      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV
[11536]1190      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
1191      !!----------------------------------------------------------------------
1192      SELECT CASE( nvor_scheme )
1193      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
[12377]1194         DO_2D_00_00
1195            z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
1196            z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
1197            zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
1198               &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
1199               &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
1200               !
1201            zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
1202               &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
1203               &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
1204         END_2D
[11536]1205         !         
1206      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
[12377]1207         DO_2D_00_00
1208            zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
1209            zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1210            zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
1211            zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1212            ! energy conserving formulation for planetary vorticity term
1213            zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1214            zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
1215         END_2D
[11536]1216         !
1217      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
[12377]1218         DO_2D_00_00
1219            zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
1220              &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
1221            zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
1222              &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
1223            zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
1224            zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
1225         END_2D
[11536]1226         !
1227      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
[12377]1228         DO_2D_00_00
1229            zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
1230             &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
1231             &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
1232             &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
1233            zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
1234             &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
1235             &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
1236             &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
1237         END_2D
[11536]1238         !
1239      END SELECT
1240      !
1241   END SUBROUTINE dyn_cor_2D
1242
1243
1244   SUBROUTINE wad_tmsk( pssh, ptmsk )
1245      !!----------------------------------------------------------------------
1246      !!                  ***  ROUTINE wad_lmt  ***
1247      !!                   
1248      !! ** Purpose :   set wetting & drying mask at tracer points
1249      !!              for the current barotropic sub-step
1250      !!
1251      !! ** Method  :   ???
1252      !!
1253      !! ** Action  :  ptmsk : wetting & drying t-mask
1254      !!----------------------------------------------------------------------
1255      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
1256      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
1257      !
1258      INTEGER  ::   ji, jj   ! dummy loop indices
1259      !!----------------------------------------------------------------------
1260      !
1261      IF( ln_wd_dl_rmp ) THEN     
[12377]1262         DO_2D_11_11
1263            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
1264               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN
1265               ptmsk(ji,jj) = 1._wp
1266            ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
1267               ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
1268            ELSE
1269               ptmsk(ji,jj) = 0._wp
1270            ENDIF
1271         END_2D
[11536]1272      ELSE 
[12377]1273         DO_2D_11_11
1274            IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
1275            ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
1276            ENDIF
1277         END_2D
[11536]1278      ENDIF
1279      !
1280   END SUBROUTINE wad_tmsk
1281
1282
1283   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1284      !!----------------------------------------------------------------------
1285      !!                  ***  ROUTINE wad_lmt  ***
1286      !!                   
1287      !! ** Purpose :   set wetting & drying mask at tracer points
1288      !!              for the current barotropic sub-step
1289      !!
1290      !! ** Method  :   ???
1291      !!
1292      !! ** Action  :  ptmsk : wetting & drying t-mask
1293      !!----------------------------------------------------------------------
1294      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
1295      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
1296      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
1297      !
1298      INTEGER  ::   ji, jj   ! dummy loop indices
1299      !!----------------------------------------------------------------------
1300      !
[12377]1301      DO_2D_11_10
1302         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
1303         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj) 
1304         ENDIF
1305         phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1306         pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1307      END_2D
[11536]1308      !
[12377]1309      DO_2D_10_11
1310         IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
1311         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1) 
1312         ENDIF
1313         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
1314         pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1315      END_2D
[11536]1316      !
1317   END SUBROUTINE wad_Umsk
1318
1319
[12377]1320   SUBROUTINE wad_spg( pshn, zcpx, zcpy )
[11536]1321      !!---------------------------------------------------------------------
1322      !!                   ***  ROUTINE  wad_sp  ***
1323      !!
1324      !! ** Purpose :
1325      !!----------------------------------------------------------------------
1326      INTEGER  ::   ji ,jj               ! dummy loop indices
1327      LOGICAL  ::   ll_tmp1, ll_tmp2
[12377]1328      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
[11536]1329      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1330      !!----------------------------------------------------------------------
[12377]1331      DO_2D_00_00
1332         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
1333              &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
1334              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji+1,jj) + ht_0(ji+1,jj) )  &
1335              &                                                         > rn_wdmin1 + rn_wdmin2
1336         ll_tmp2 = ( ABS( pshn(ji+1,jj)            -  pshn(ji  ,jj))  > 1.E-12 ).AND.( &
1337              &      MAX(   pshn(ji,jj)              ,  pshn(ji+1,jj) ) >                &
1338              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1339         IF(ll_tmp1) THEN
1340            zcpx(ji,jj) = 1.0_wp
1341         ELSEIF(ll_tmp2) THEN
1342            ! no worries about  pshn(ji+1,jj) -  pshn(ji  ,jj) = 0, it won't happen ! here
1343            zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
1344                 &           / (pshn(ji+1,jj) - pshn(ji  ,jj)) )
1345            zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1346         ELSE
1347            zcpx(ji,jj) = 0._wp
1348         ENDIF
1349         !
1350         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji,jj+1) ) >                &
1351              &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
1352              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji,jj+1) + ht_0(ji,jj+1) )  &
1353              &                                                       > rn_wdmin1 + rn_wdmin2
1354         ll_tmp2 = ( ABS( pshn(ji,jj)              -  pshn(ji,jj+1))  > 1.E-12 ).AND.( &
1355              &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
1356              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1357         
1358         IF(ll_tmp1) THEN
1359            zcpy(ji,jj) = 1.0_wp
1360         ELSE IF(ll_tmp2) THEN
1361            ! no worries about  pshn(ji,jj+1) -  pshn(ji,jj  ) = 0, it won't happen ! here
1362            zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
1363                 &           / (pshn(ji,jj+1) - pshn(ji,jj  )) )
1364            zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
1365         ELSE
1366            zcpy(ji,jj) = 0._wp
1367         ENDIF
1368      END_2D
[11536]1369           
1370   END SUBROUTINE wad_spg
1371     
1372
1373
[12377]1374   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
[11536]1375      !!----------------------------------------------------------------------
1376      !!                  ***  ROUTINE dyn_drg_init  ***
1377      !!                   
1378      !! ** Purpose : - add the baroclinic top/bottom drag contribution to
1379      !!              the baroclinic part of the barotropic RHS
1380      !!              - compute the barotropic drag coefficients
1381      !!
1382      !! ** Method  :   computation done over the INNER domain only
1383      !!----------------------------------------------------------------------
[12377]1384      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
1385      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
1386      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
1387      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
1388      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
[11536]1389      !
1390      INTEGER  ::   ji, jj   ! dummy loop indices
1391      INTEGER  ::   ikbu, ikbv, iktu, iktv
1392      REAL(wp) ::   zztmp
1393      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
1394      !!----------------------------------------------------------------------
1395      !
1396      !                    !==  Set the barotropic drag coef.  ==!
1397      !
1398      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
1399         
[12377]1400         DO_2D_00_00
1401            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1402            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1403         END_2D
[11536]1404      ELSE                          ! bottom friction only
[12377]1405         DO_2D_00_00
1406            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1407            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1408         END_2D
[11536]1409      ENDIF
1410      !
1411      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
1412      !
1413      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
1414         
[12377]1415         DO_2D_00_00
1416            ikbu = mbku(ji,jj)       
1417            ikbv = mbkv(ji,jj)   
1418            zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
1419            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
1420         END_2D
[11536]1421      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
1422         
[12377]1423         DO_2D_00_00
1424            ikbu = mbku(ji,jj)       
1425            ikbv = mbkv(ji,jj)   
1426            zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
1427            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
1428         END_2D
[11536]1429      ENDIF
1430      !
1431      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
[12443]1432         zztmp = -1._wp / rDt_e
[12377]1433         DO_2D_00_00
1434            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
1435                 &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
1436            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
1437                 &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
1438         END_2D
[11536]1439      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1440         
[12377]1441         DO_2D_00_00
1442            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
1443            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
1444         END_2D
[11536]1445      END IF
1446      !
1447      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
1448      !
1449      IF( ln_isfcav ) THEN
1450         !
1451         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
1452           
[12377]1453            DO_2D_00_00
1454               iktu = miku(ji,jj)
1455               iktv = mikv(ji,jj)
1456               zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
1457               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
1458            END_2D
[11536]1459         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
1460           
[12377]1461            DO_2D_00_00
1462               iktu = miku(ji,jj)
1463               iktv = mikv(ji,jj)
1464               zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
1465               zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
1466            END_2D
[11536]1467         ENDIF
1468         !
1469         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1470         
[12377]1471         DO_2D_00_00
1472            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
1473            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
1474         END_2D
[11536]1475         !
1476      ENDIF
1477      !
1478   END SUBROUTINE dyn_drg_init
1479
1480   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
1481      &                      za0, za1, za2, za3 )   ! ==> out
1482      !!----------------------------------------------------------------------
1483      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
1484      LOGICAL ,INTENT(in   ) ::   ll_init              !
1485      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
1486      !
1487      REAL(wp) ::   zepsilon, zgamma                   !   -      -
1488      !!----------------------------------------------------------------------
1489      !                             ! set Half-step back interpolation coefficient
1490      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
1491         za0 = 1._wp                       
1492         za1 = 0._wp                           
1493         za2 = 0._wp
1494         za3 = 0._wp
1495      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1496         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
1497         za1 =-0.1666666666666_wp                 ! za1 = gam
1498         za2 = 0.0833333333333_wp                 ! za2 = eps
1499         za3 = 0._wp             
1500      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1501         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion 
1502            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
1503            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
1504            za2 = 0.088_wp                        ! za2 = gam
1505            za3 = 0.013_wp                        ! za3 = eps
1506         ELSE                                 ! no time diffusion
1507            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1508            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1509            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1510            za1 = 1._wp - za0 - zgamma - zepsilon
1511            za2 = zgamma
1512            za3 = zepsilon
1513         ENDIF
1514      ENDIF
1515   END SUBROUTINE ts_bck_interp
1516
1517
[358]1518   !!======================================================================
1519END MODULE dynspg_ts
Note: See TracBrowser for help on using the repository browser.