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/releases/r4.0/r4.0-HEAD/src/OCE/DYN – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynspg_ts.F90 @ 12737

Last change on this file since 12737 was 12737, checked in by jchanut, 4 years ago

Fixes AGRIF reproductibility with land processors removal, i.e. #2240. Trunk is not concerned by this problem since nbondi/nbondj variables are not used anymore.

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