Changeset 3970


Ignore:
Timestamp:
2013-07-11T15:59:14+02:00 (7 years ago)
Author:
cbricaud
Message:

Time splitting update, see ticket #1079

Location:
branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC
Files:
18 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3851 r3970  
    2929   USE iom             ! IOM library 
    3030   USE in_out_manager  ! I/O logical units 
     31   ! bg jchanut tschanges 
     32   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 
     33   ! end jchanut tschanges 
     34 
    3135#if defined key_lim2 
    3236   USE ice_2 
     
    314318                     END DO 
    315319                  ENDIF 
    316                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    317                      CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  & 
    318                                         & td=tides(ib_bdy), time_offset=time_offset ) 
    319                   ENDIF 
     320                  ! bg jchanut tschanges 
     321                  !IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     322                  !   CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  & 
     323                  !                      & td=tides(ib_bdy), time_offset=time_offset ) 
     324                  !ENDIF 
     325                  ! end jchanut tschanges 
    320326               ENDIF 
    321327            ENDIF 
     
    323329         END IF ! nn_dta(ib_bdy) = 1 
    324330      END DO  ! ib_bdy 
     331 
     332      ! bg jchanut tschanges 
     333#if defined key_tide 
     334      ! Add tides if not split-explicit free surface else this is done in ts loop 
     335      IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
     336#endif 
     337      ! end jchanut tschanges 
    325338 
    326339      IF ( ln_apr_obc ) THEN 
     
    476489            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
    477490 
    478                IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
     491               IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
    479492                  jfld = jfld + 1 
    480493                  blf_i(jfld) = bn_ssh 
     
    572585            ! Recalculate field counts 
    573586            !------------------------- 
    574             nb_bdy_fld_sum = 0 
    575587            IF( ib_bdy .eq. 1 ) THEN  
     588               nb_bdy_fld_sum = 0 
    576589               nb_bdy_fld(ib_bdy) = jfld 
    577590               nb_bdy_fld_sum     = jfld               
     
    616629               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
    617630               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    618                IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN 
     631               IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN 
    619632                  jfld = jfld + 1 
    620633                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r3294 r3970  
    3030   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3131   USE in_out_manager  ! 
     32   USE domvvl          ! variable volume 
    3233 
    3334   IMPLICIT NONE 
     
    8485      pu2d(:,:) = 0.e0 
    8586      pv2d(:,:) = 0.e0 
    86       DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    87           pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    88           pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    89       END DO 
    90       pu2d(:,:) = pu2d(:,:) * phur(:,:) 
    91       pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     87      ! bg jchanut tschanges (not specifically related to ts; this is a bug) 
     88      IF (lk_vvl) THEN 
     89         DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     90            pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     91            pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     92         END DO 
     93         pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) 
     94         pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     95      ! end jchanut tschanges 
     96      ELSE 
     97         DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     98            pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     99            pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     100         END DO 
     101         pu2d(:,:) = pu2d(:,:) * phur(:,:) 
     102         pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     103      ENDIF 
     104 
    92105      DO jk = 1 , jpkm1 
    93          ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) 
    94          va(:,:,jk) = va(:,:,jk) - pv2d(:,:) 
     106         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk) 
     107         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk) 
    95108      END DO 
    96109 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r3970  
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite 
    77   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
     8   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_bdy  
     
    2728 
    2829   PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn 
     30   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv 
    2931 
    3032   !!---------------------------------------------------------------------- 
     
    135137      REAL(wp) ::   zcorr                            ! Flather correction 
    136138      REAL(wp) ::   zforc                            ! temporary scalar 
     139      REAL(wp) ::   zflag, z1_2                      !    "        " 
    137140      !!---------------------------------------------------------------------- 
    138141 
    139142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 
     143 
     144      z1_2 = 0.5_wp 
    140145 
    141146      ! ---------------------------------! 
     
    164169         ! 
    165170         zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    166          zforc = dta%u2d(jb) 
    167          pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
     171         ! bg jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 
     172!!         zforc = dta%u2d(jb) 
     173         zflag = ABS(idx%flagu(jb)) 
     174         iim1 = ii + idx%flagu(jb) 
     175         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pu2d(iim1,ij) 
     176         pu2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
     177         ! end jchanut tschanges 
    168178      END DO 
    169179      ! 
     
    177187         ! 
    178188         zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    179          zforc = dta%v2d(jb) 
    180          pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
     189         ! bg jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 
     190!!         zforc = dta%v2d(jb) 
     191         zflag = ABS(idx%flagv(jb)) 
     192         ijm1 = ij + idx%flagv(jb) 
     193         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pv2d(ii,ijm1) 
     194         pv2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
     195         ! end jchanut tschanges 
    181196      END DO 
    182197      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     
    186201      ! 
    187202   END SUBROUTINE bdy_dyn2d_fla 
     203 
     204   SUBROUTINE bdy_ssh( zssh ) 
     205      !!---------------------------------------------------------------------- 
     206      !!                  ***  SUBROUTINE bdy_ssh  *** 
     207      !! 
     208      !! ** Purpose : Duplicate sea level across open boundaries 
     209      !! 
     210      !!---------------------------------------------------------------------- 
     211      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level 
     212      !! 
     213      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers 
     214      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       " 
     215 
     216      igrd = 1                       ! Everything is at T-points here 
     217 
     218      DO ib_bdy = 1, nb_bdy 
     219         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     220            ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     221            ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     222            ! Set gradient direction: 
     223            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
     224            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
     225            IF ( zcoef1+zcoef2 == 0 ) THEN 
     226               ! corner 
     227!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1) 
     228!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + & 
     229!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + & 
     230!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + & 
     231!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1) 
     232               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1) 
     233               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + & 
     234                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + & 
     235                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + & 
     236                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1) 
     237               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 
     238            ELSE 
     239               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     240               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     241               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 
     242            ENDIF 
     243         END DO 
     244 
     245         ! Boundary points should be updated 
     246         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy ) 
     247      END DO 
     248 
     249   END SUBROUTINE bdy_ssh 
    188250#else 
    189251   !!---------------------------------------------------------------------- 
     
    192254CONTAINS 
    193255   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine 
    194       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
     256      INTEGER, intent(in) :: kt 
     257      WRITE(*,*) 'bdy_dyn2: You should not have seen this print! error?', kt 
    195258   END SUBROUTINE bdy_dyn2d 
     259 
    196260#endif 
    197261 
    198262   !!====================================================================== 
    199263END MODULE bdydyn2d 
     264 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3703 r3970  
    161161 
    162162      DO ib_bdy = 1,nb_bdy 
     163        icount = 0 ! flag to set max rimwidth to 1 if no relaxation 
    163164        IF(lwp) WRITE(numout,*) ' '  
    164165        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
     
    175176          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    176177          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     178          icount = icount + 1 
    177179          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    178180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
     
    196198          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    197199          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     200          icount = icount + 1 
    198201          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    199202          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     
    215218              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    216219           ELSE 
     220              icount = icount + 1 
    217221              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    218222              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     
    228232          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    229233          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     234          icount = icount + 1 
    230235          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    231236          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     
    248253              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    249254           ELSE 
     255              icount = icount + 1 
    250256              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    251257              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     
    262268          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    263269          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     270          icount = icount + 1 
    264271          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    265272        END SELECT 
     
    273280        IF(lwp) WRITE(numout,*) 
    274281#endif 
    275  
    276         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    277         IF(lwp) WRITE(numout,*) 
     282        IF ( icount>0 ) THEN 
     283           IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     284           IF(lwp) WRITE(numout,*) 
     285        ELSE 
     286           nn_rimwidth(ib_bdy) = 1 ! no relaxation 
     287        ENDIF 
    278288 
    279289      ENDDO 
     
    391401            ENDDO 
    392402            CALL iom_close( inum ) 
    393  
    394403         ENDIF  
    395404 
     
    398407      IF (nb_bdy>0) THEN 
    399408         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    400  
    401409         ! Allocate arrays 
    402410         !--------------- 
     
    446454         ENDIF  
    447455 
    448       ENDDO       
     456      ENDDO      
    449457     
    450458      ! 2. Now fill indices corresponding to straight open boundary arrays: 
     
    752760               ! check if point is in local domain 
    753761               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
    754                   & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN 
     762                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   & 
     763                  & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy)     ) THEN       
    755764                  ! 
    756765                  icount = icount  + 1 
     
    765774         ! Allocate index arrays for this boundary set 
    766775         !-------------------------------------------- 
    767          ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 
     776 
     777         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd)) 
     778         ilen1 = MAX(1,ilen1) 
    768779         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 
    769780         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
     
    773784         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    774785         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
    775          ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     786         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )       
    776787 
    777788         ! Dispatch mapping indices and discrete distances on each processor 
     
    10081019      ! bdytmask = 1  on the computational domain AND on open boundaries 
    10091020      !          = 0  elsewhere    
    1010   
     1021      bdytmask(:,:) = 1.e0 
     1022      bdyumask(:,:) = 1.e0 
     1023      bdyvmask(:,:) = 1.e0 
     1024 
    10111025      IF( ln_mask_file ) THEN 
    10121026         CALL iom_open( cn_mask_file, inum ) 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r3651 r3970  
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    1010   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
     11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_bdy 
     
    3233!   USE tide_mod       ! Useless ?? 
    3334   USE fldread, ONLY: fld_map 
     35   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3436 
    3537   IMPLICIT NONE 
     
    3840   PUBLIC   bdytide_init     ! routine called in bdy_init 
    3941   PUBLIC   bdytide_update   ! routine called in bdy_dta 
     42   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    4043 
    4144   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     
    4952 
    5053   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
     54   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    5155 
    5256   !!---------------------------------------------------------------------- 
     
    252256            ENDIF 
    253257            ! 
     258            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 
     259                                     ! time splitting integration 
     260               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
     261               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
     262               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
     263               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 
     264               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 
     265               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 
     266            ENDIF 
     267            ! 
    254268         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
    255269         ! 
     
    318332          
    319333      IF( PRESENT(jit) ) THEN   
    320          z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 
     334         z_arg = ((kt-kt_tide) * rdt + (jit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    321335      ELSE                               
    322336         z_arg = ((kt-kt_tide)+time_add) * rdt 
     
    324338 
    325339      ! Linear ramp on tidal component at open boundaries  
    326       zramp = 1. 
    327       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     340      zramp = 1._wp 
     341      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
    328342 
    329343      DO itide = 1, nb_harmo 
     
    351365      ! 
    352366   END SUBROUTINE bdytide_update 
     367 
     368   SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 
     369      !!---------------------------------------------------------------------- 
     370      !!                 ***  SUBROUTINE bdy_dta_tides  *** 
     371      !!                 
     372      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     373      !!                 
     374      !!---------------------------------------------------------------------- 
     375      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter 
     376      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option) 
     377      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit 
     378                                                        ! is present then units = subcycle timesteps. 
     379                                                        ! time_offset = 0  => get data at "now"    time level 
     380                                                        ! time_offset = -1 => get data at "before" time level 
     381                                                        ! time_offset = +1 => get data at "after"  time level 
     382                                                        ! etc. 
     383      !! 
     384      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step 
     385      INTEGER,          DIMENSION(jpbgrd) :: ilen0  
     386      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts 
     387      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices 
     388      INTEGER  :: time_add                                    ! time offset in units of timesteps 
     389      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist       
     390      !!---------------------------------------------------------------------- 
     391 
     392      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides') 
     393 
     394      lk_first_btstp=.TRUE. 
     395      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF 
     396 
     397      time_add = 0 
     398      IF( PRESENT(time_offset) ) THEN 
     399         time_add = time_offset 
     400      ENDIF 
     401       
     402      ! Absolute time from model initialization:    
     403      IF( PRESENT(kit) ) THEN   
     404         z_arg = ( kt + (kit+0.5*(time_add-1)) / REAL(nn_baro,wp) ) * rdt 
     405      ELSE                               
     406         z_arg = ( kt + time_add ) * rdt 
     407      ENDIF 
     408 
     409      ! Linear ramp on tidal component at open boundaries  
     410      zramp = 1. 
     411      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     412 
     413      DO ib_bdy = 1,nb_bdy 
     414 
     415         ! line below should be simplified (runoff case) 
     416         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN 
     417 
     418            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd) 
     419            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd) 
     420 
     421            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     422               ilen0(:)=nblen(:) 
     423            ELSE 
     424               ilen0(:)=nblenrim(:) 
     425            ENDIF      
     426 
     427            ! We refresh nodal factors every day below 
     428            ! This should be done somewhere else 
     429            IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. lk_first_btstp ) THEN 
     430               ! 
     431               kt_tide = kt                
     432               ! 
     433               IF(lwp) THEN 
     434               WRITE(numout,*) 
     435               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt 
     436               WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     437               ENDIF 
     438               ! 
     439               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     440               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     441               ! 
     442            ENDIF 
     443            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
     444            ! 
     445            ! If time splitting, save data at first barotropic iteration 
     446            IF ( PRESENT(kit) ) THEN 
     447               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 
     448                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 
     449                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 
     450                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 
     451 
     452               ELSE ! Initialize arrays from slow varying open boundary data:             
     453                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
     454                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
     455                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
     456               ENDIF 
     457            ENDIF 
     458            ! 
     459            ! Update open boundary data arrays: 
     460            DO itide = 1, nb_harmo 
     461               ! 
     462               z_sarg = (z_arg + zoff) * omega_tide(itide) 
     463               z_cost = zramp * COS( z_sarg ) 
     464               z_sist = zramp * SIN( z_sarg ) 
     465               ! 
     466               igrd=1                              ! SSH on tracer grid 
     467               DO ib = 1, ilen0(igrd) 
     468                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 
     469                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & 
     470                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist ) 
     471               END DO 
     472               ! 
     473               igrd=2                              ! U grid 
     474               DO ib = 1, ilen0(igrd) 
     475                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 
     476                     &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 
     477                     &                        tides(ib_bdy)%u(ib,itide,2)*z_sist ) 
     478               END DO 
     479               ! 
     480               igrd=3                              ! V grid 
     481               DO ib = 1, ilen0(igrd)  
     482                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 
     483                     &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & 
     484                     &                        tides(ib_bdy)%v(ib,itide,2)*z_sist ) 
     485               END DO 
     486            END DO 
     487         END IF 
     488      END DO 
     489      ! 
     490      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides') 
     491      ! 
     492   END SUBROUTINE bdy_dta_tides 
    353493 
    354494   SUBROUTINE tide_init_elevation( idx, td ) 
     
    457597      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
    458598   END SUBROUTINE bdytide_update 
     599   SUBROUTINE bdy_dta_tides( kt, jit )   ! Empty routine 
     600      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit 
     601   END SUBROUTINE bdy_dta_tides 
    459602#endif 
    460603 
    461604   !!====================================================================== 
    462605END MODULE bdytides 
     606 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r3851 r3970  
    158158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m) 
    159159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu   , hv     !: depth at u- and v-points (meters) 
    160    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters) 
    161  
     160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0   !: reference depth at u- and v-points (meters) 
     161   ! bg jchanut tschanges  
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ht   , hf     !: depth at t- and f-points (meters) 
     163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ht_0 , hf_0   !: reference depth at t- and f-points (meters) 
     164   ! end jchanut tschanges  
    162165   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    163166   INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1)  
     
    264267   INTEGER FUNCTION dom_oce_alloc() 
    265268      !!---------------------------------------------------------------------- 
    266       INTEGER, DIMENSION(11) :: ierr 
     269      ! bg jchanut tschanges 
     270      INTEGER, DIMENSION(13) :: ierr 
     271      ! end jchanut tschanges 
    267272      !!---------------------------------------------------------------------- 
    268273      ierr(:) = 0 
     
    314319      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) ) 
    315320#endif 
     321      ! bg jchanut tschanges 
     322#if defined key_dynspg_ts 
     323      ALLOCATE( ht(jpi,jpj)  , hf(jpi,jpj)  , STAT=ierr(12) ) 
     324#if defined key_vvl 
     325      ALLOCATE( ht_0(jpi,jpj), hf_0(jpi,jpj), STAT=ierr(13) ) 
     326#endif 
     327#endif 
     328      ! end jchanut tschanges 
    316329      ! 
    317330      dom_oce_alloc = MAXVAL(ierr) 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r3764 r3970  
    6868      !!              - 1D configuration, move Coriolis, u and v at T-point 
    6969      !!---------------------------------------------------------------------- 
    70       INTEGER ::   jk          ! dummy loop argument 
     70      INTEGER ::   jj, jk      ! dummy loop arguments 
    7171      INTEGER ::   iconf = 0   ! local integers 
    7272      !!---------------------------------------------------------------------- 
     
    100100      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
    101101 
     102      ! bg jchanut tschanges 
     103#if defined key_dynspg_ts 
     104      ! 
     105      ht(:,:) = 0._wp                         ! Ocean depth at T-point 
     106      DO jk = 1, jpk 
     107         ht(:,:) = ht(:,:) + fse3t(:,:,jk) * tmask(:,:,jk) 
     108      END DO  
     109 
     110                                              ! Ocean depth at F-point  
     111                                              ! Ensure that depth is non-zero over land 
     112      IF ( .not. ln_sco ) THEN 
     113         IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     114         ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     115         ENDIF 
     116         hf(:,:) = gdepw_0(jk+1) 
     117      ELSE 
     118         hf(:,:) = hbatf(:,:) 
     119      END IF 
     120 
     121      DO jj = 1, jpjm1 
     122         hf(:,jj) = hf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     123      END DO  
     124             
     125      DO jk = 1, jpkm1 
     126         DO jj = 1, jpjm1 
     127            hf(:,jj) = hf(:,jj) + fse3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     128         END DO 
     129      END DO 
     130      CALL lbc_lnk( hf, 'F', 1._wp )   
     131#endif  
     132      ! end jchanut tschanges 
    102133                             CALL dom_stp      ! time step 
    103134      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file 
     
    218249         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt 
    219250         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp 
    220          WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro 
     251         ! bg jchanut tschanges 
     252!         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro 
     253         ! end jchanut tschanges 
    221254         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc 
    222255         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r3294 r3970  
    141141      END DO 
    142142       
     143#if defined key_dynspg_ts 
     144      ! bg jchanut tschanges 
     145      ht_0(:,:) = 0._wp                         ! Reference ocean depth at T-points 
     146      DO jk = 1, jpk 
     147         ht_0(:,:) = ht_0(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 
     148      END DO  
     149                                                ! Reference ocean depth at F-points  
     150                                                ! Ensure that depth is non-zero over land 
     151      IF ( .not. ln_sco ) THEN 
     152         IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     153         ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     154         ENDIF 
     155         hf_0(:,:) = gdepw_0(jk+1) 
     156      ELSE 
     157         hf_0(:,:) = hbatf(:,:) 
     158      END IF 
     159 
     160      DO jj = 1, jpjm1 
     161         hf_0(:,jj) = hf_0(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     162      END DO   
     163            
     164      DO jk = 1, jpkm1 
     165         DO jj = 1, jpjm1 
     166            hf_0(:,jj) = hf_0(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     167         END DO 
     168      END DO 
     169      CALL lbc_lnk( hf_0, 'F', 1._wp ) 
     170      ! end jchanut tschanges 
     171#endif   
     172 
    143173      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points 
    144174         DO ji = 1, jpim1   ! NO vector opt. 
     
    192222      INTEGER  ::   iku, ikv     ! local integers     
    193223      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
    194       REAL(wp) ::   zvt          ! local scalars 
     224      REAL(wp) ::   zvt, zvtip1, zvtjp1  ! local scalars 
    195225      !!---------------------------------------------------------------------- 
    196226      ! 
     
    202232         WRITE(numout,*) '~~~~~~~~~ ' 
    203233         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 
    204          pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 
     234         pe3v_b(:,:,jpk) = fse3v_0(:,:,jpk) 
    205235      ENDIF 
    206236       
     
    208238         DO jj = 2, jpjm1 
    209239            DO ji = fs_2, fs_jpim1 
    210                zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 
    211                pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    212                pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     240               zvt    = ( fse3t_b(ji  ,jj  ,jk) - fse3t_0(ji  ,jj  ,jk) ) * e1e2t(ji  ,jj  ) 
     241               zvtip1 = ( fse3t_b(ji+1,jj  ,jk) - fse3t_0(ji+1,jj  ,jk) ) * e1e2t(ji+1,jj  ) 
     242               zvtjp1 = ( fse3t_b(ji  ,jj+1,jk) - fse3t_0(ji  ,jj+1,jk) ) * e1e2t(ji  ,jj+1) 
     243               pe3u_b(ji,jj,jk) = fse3u_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtip1 ) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     244               pe3v_b(ji,jj,jk) = fse3v_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtjp1 ) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    213245            END DO 
    214246         END DO 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

    r3294 r3970  
    2929 
    3030   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS 
    31    REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8 4th order centred 
     31   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred 
    3232 
    3333   PUBLIC   dyn_adv_ubs   ! routine called by step.F90 
     
    5757      !!                       = 1/3  3rd order Upstream biased scheme 
    5858      !!                gamma2 = 0    2nd order finite differencing  
    59       !!                       = 1/8 4th order finite differencing 
     59      !!                       = 1/32 4th order finite differencing 
    6060      !!      For stability reasons, the first term of the fluxes which cor- 
    6161      !!      responds to a second order centered scheme is evaluated using   
     
    6464      !!      before velocity (forward in time).  
    6565      !!      Default value (hard coded in the begining of the module) are  
    66       !!      gamma1=1/3 and gamma2=1/8. 
     66      !!      gamma1=1/3 and gamma2=1/32. 
    6767      !! 
    6868      !! ** Action : - (ua,va) updated with the 3D advective momentum trends 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r3764 r3970  
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
    1818   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 
     19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1920   !!------------------------------------------------------------------------- 
    2021   
     
    4243   USE wrk_nemo        ! Memory Allocation 
    4344   USE prtctl          ! Print control 
     45   USE dynspg_ts       ! Barotropic velocities 
     46 
    4447#if defined key_agrif 
    4548   USE agrif_opa_interp 
     
    103106      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
    104107      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
     108      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva 
    105109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
    106110      !!---------------------------------------------------------------------- 
     
    109113      ! 
    110114      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     115      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva ) 
    111116      ! 
    112117      IF( kt == nit000 ) THEN 
     
    127132      ! 
    128133#else 
     134 
     135#if defined key_dynspg_exp 
    129136      ! Next velocity :   Leap-frog time stepping 
    130137      ! ------------- 
     
    147154         END DO 
    148155      ENDIF 
    149  
     156#endif 
     157 
     158#if defined key_dynspg_ts 
     159      ! Ensure below that barotropic velocities match time splitting estimate 
     160      ! Compute actual transport and replace it with ts estimate at "after" time step 
     161      zua(:,:) = 0._wp 
     162      zva(:,:) = 0._wp 
     163      IF (lk_vvl) THEN 
     164         DO jk = 1, jpkm1 
     165            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     166            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)     
     167         END DO 
     168         DO jk = 1, jpkm1 
     169            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) & 
     170            & / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) + ua_b(:,:) ) * umask(:,:,jk) 
     171            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) &  
     172            & / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) + va_b(:,:) ) * vmask(:,:,jk) 
     173         END DO 
     174      ELSE 
     175         DO jk = 1, jpkm1 
     176            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     177            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)     
     178         END DO 
     179         DO jk = 1, jpkm1 
     180            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk) 
     181            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk) 
     182         END DO 
     183      ENDIF 
     184 
     185      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 
     186         ! Remove advective velocity from "now velocities"  
     187         ! prior to asselin filtering      
     188         ! In the forward case, this is done below after asselin filtering     
     189         DO jk = 1, jpkm1 
     190            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
     191            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     192         END DO   
     193      ENDIF 
     194#endif 
    150195 
    151196      ! Update after velocity on domain lateral boundaries 
     
    194239            vn(:,:,jk) = va(:,:,jk) 
    195240         END DO 
     241         IF (lk_vvl) THEN 
     242            DO jk = 1, jpkm1 
     243               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     244               fse3u_b(:,:,jk) = fse3u_n(:,:,jk) 
     245               fse3v_b(:,:,jk) = fse3v_n(:,:,jk) 
     246            ENDDO 
     247         ENDIF 
    196248      ELSE                                             !* Leap-Frog : Asselin filter and swap 
    197249         !                                ! =============! 
     
    215267            !                             ! ================! 
    216268            ! 
    217             DO jk = 1, jpkm1                 ! Before scale factor at t-points 
    218                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    219                   &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
    220                   &                         - 2._wp * fse3t_n(:,:,jk)            ) 
    221             END DO 
    222             zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors 
    223             fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     269            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     270               ! Remove asselin filtering on thicknesses if forward time splitting 
     271               DO jk = 1, jpkm1    
     272                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     273               ENDDO 
     274            ELSE 
     275 
     276               DO jk = 1, jpkm1                 ! Before scale factor at t-points 
     277                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
     278                     &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
     279                     &                         - 2._wp * fse3t_n(:,:,jk)            ) 
     280               END DO 
     281               zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors 
     282               fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     283            ENDIF 
    224284            ! 
    225285            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation) 
     
    272332         ENDIF 
    273333         ! 
    274       ENDIF 
     334#if defined key_dynspg_ts 
     335         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     336         ! Remove asselin filtering of barotropic velocities if forward time splitting 
     337         ! note that we replace barotropic velocities by advective velocities        
     338            zua(:,:) = 0._wp 
     339            zva(:,:) = 0._wp 
     340            IF (lk_vvl) THEN 
     341               DO jk = 1, jpkm1 
     342                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     343                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     344               END DO 
     345               DO jk = 1, jpkm1 
     346                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) & 
     347                  & / ( hu_0(:,:) + sshu_n(:,:) + 1._wp - umask(:,:,1) ) - un_b(:,:)) * umask(:,:,jk) 
     348                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) &  
     349                  & / ( hv_0(:,:) + sshv_n(:,:) + 1._wp - vmask(:,:,1) ) - vn_b(:,:)) * vmask(:,:,jk) 
     350               END DO 
     351            ELSE 
     352               DO jk = 1, jpkm1 
     353                  zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     354                  zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     355               END DO 
     356               DO jk = 1, jpkm1 
     357                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 
     358                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
     359               END DO 
     360            ENDIF 
     361         ENDIF 
     362#endif 
     363         ! 
     364      ENDIF ! neuler =/0 
    275365 
    276366      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
     
    278368      !  
    279369      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     370      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva ) 
    280371      ! 
    281372      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r3625 r3970  
    2727   USE in_out_manager ! I/O manager 
    2828   USE lib_mpp        ! MPP library 
    29    USE solver          ! solver initialization 
    30    USE wrk_nemo        ! Memory Allocation 
    31    USE timing          ! Timing 
    32  
     29   USE solver         ! solver initialization 
     30   USE wrk_nemo       ! Memory Allocation 
     31   USE timing         ! Timing 
     32   USE dynhpg, ONLY: ln_dynhpg_imp ! semi-implicit hpg flag 
    3333 
    3434   IMPLICIT NONE 
     
    9999         ztrdv(:,:,:) = va(:,:,:) 
    100100      ENDIF 
    101  
    102       IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==! 
     101      ! bg jchanut tschanges 
     102      IF( ( ln_apr_dyn ).AND.(.NOT.lk_dynspg_ts) ) THEN  !==  Atmospheric pressure gradient 
     103      ! (According to time splitting option, lines below are done latter) 
     104      ! end jchanut tschanges 
    103105         zg_2 = grav * 0.5 
    104106         DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     
    207209         WRITE(numout,*) '     Filtered free surface cst volume       lk_dynspg_flt = ', lk_dynspg_flt 
    208210      ENDIF 
    209  
     211      ! 
     212      ! bg jchanut tschanges 
     213      IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 
     214      ! (do it now, to set nn_baro, used to allocate some arrays later on) 
     215      ! end jchanut tschanges 
     216      !  
    210217      !                        ! allocate dyn_spg arrays 
    211218      IF( lk_dynspg_ts ) THEN 
     
    247254      ENDIF 
    248255 
    249       !                        ! Control of momentum formulation 
    250       IF( lk_dynspg_ts .AND. lk_vvl ) THEN 
    251          IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' ) 
    252       ENDIF 
     256      ! bg jchanut tschanges 
     257      !               ! Control of hydrostatic pressure choice 
     258      IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 
     259         CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
     260      ENDIF 
     261!      !                        ! Control of momentum formulation 
     262!      IF( lk_dynspg_ts .AND. lk_vvl ) THEN 
     263!         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' ) 
     264!      ENDIF 
     265      ! end jchanut tschanges 
    253266      ! 
    254267      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init') 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3680 r3970  
    99   !!             3.3  ! 2010-09  (D. Storkey, E. O'Dea) update for BDY for Shelf configurations 
    1010   !!             3.3  ! 2011-03  (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b 
     11   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
    1112   !!--------------------------------------------------------------------- 
    1213#if defined key_dynspg_ts   ||   defined key_esopa 
     
    1617   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    1718   !!                 splitting scheme and add to the general trend  
    18    !!   ts_rst      : read/write the time-splitting restart fields in the ocean restart file 
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce             ! ocean dynamics and tracers 
     
    2424   USE phycst          ! physical constants 
    2525   USE domvvl          ! variable volume 
    26    USE zdfbfr          ! bottom friction 
    2726   USE dynvor          ! vorticity term 
    28    USE obc_oce         ! Lateral open boundary condition 
    29    USE obc_par         ! open boundary condition parameters 
    30    USE obcdta          ! open boundary condition data      
    31    USE obcfla          ! Flather open boundary condition   
    3227   USE bdy_par         ! for lk_bdy 
    3328   USE bdy_oce         ! Lateral open boundary condition 
    34    USE bdydta          ! open boundary condition data      
     29   USE bdytides        ! open boundary condition data      
    3530   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    36    USE sbctide 
    37    USE updtide 
     31   USE sbctide         ! tides 
     32   USE updtide         ! tide potential 
    3833   USE lib_mpp         ! distributed memory computing library 
    3934   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    4136   USE in_out_manager  ! I/O manager 
    4237   USE iom             ! IOM library 
     38   USE restart         ! only for lrst_oce 
    4339   USE zdf_oce         ! Vertical diffusion 
    4440   USE wrk_nemo        ! Memory Allocation 
    45    USE timing          ! Timing 
     41   USE timing          ! Timing     
     42   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     43   USE dynadv, ONLY: ln_dynadv_vec 
    4644 
    4745 
     
    4947   PRIVATE 
    5048 
    51    PUBLIC dyn_spg_ts        ! routine called by step.F90 
    52    PUBLIC ts_rst            ! routine called by istate.F90 
    53    PUBLIC dyn_spg_ts_alloc  ! routine called by dynspg.F90 
    54  
    55  
     49   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90  
     50   PUBLIC dyn_spg_ts_alloc  !    "      "     "    " 
     51   PUBLIC dyn_spg_ts_init   !    "      "     "    " 
     52 
     53   ! Potential namelist parameters below to be read in dyn_spg_ts_init 
     54   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping 
     55   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables 
     56   LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically 
     57   INTEGER,  PRIVATE, PARAMETER :: nn_bt_flt=1            !: Filter choice 
     58   REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp      !: Max. courant number (used if ln_bt_nn_auto=T) 
     59   !!INTEGER, PUBLIC :: nn_baro=30                          !: Number of barotropic time steps / baroclinic step (rdt) 
     60   ! End namelist parameters 
     61 
     62   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 
     63   REAL(wp),SAVE :: rdtbt   ! Barotropic time step 
     64 
     65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 
     66                    wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables 
     67                    wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables 
     68 
     69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points 
    5670   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter 
    5771   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
    5872 
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b   ! now    averaged velocity 
    60    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b   ! before averaged velocity 
     73   ! Would be convenient to have arrays below defined whatever the free surface option ? 
     74   ! These could be computed once for all at the beginning of the each baroclinic time step 
     75   ! and eventually swapped at the end 
     76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b     ! after  averaged velocities 
     77   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b     ! now    averaged velocities 
     78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b     ! before averaged velocities 
     79 
     80   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step 
     81   REAL(wp), PRIVATE,ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step 
     82 
     83   ! Arrays below are saved to allow testing of the "no time averaging" option 
     84   ! If this option is not retained, these could be replaced by temporary arrays 
     85   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 
     86                                                   ubb_e, ub_e,     & 
     87                                                   vbb_e, vb_e 
    6188 
    6289   !! * Substitutions 
     
    6491#  include "vectopt_loop_substitute.h90" 
    6592   !!---------------------------------------------------------------------- 
    66    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    67    !! $Id$ 
     93   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
     94   !! $Id: dynspg_ts.F90 
    6895   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6996   !!---------------------------------------------------------------------- 
     
    74101      !!                  ***  routine dyn_spg_ts_alloc  *** 
    75102      !!---------------------------------------------------------------------- 
    76       ALLOCATE( ftnw  (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) ,     & 
    77          &      ftsw  (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc ) 
    78          ! 
     103      INTEGER :: ierr(3) 
     104      !!---------------------------------------------------------------------- 
     105      ierr(:) = 0 
     106 
     107      ALLOCATE( ub_b(jpi,jpj)  , vb_b(jpi,jpj)   , & 
     108         &      un_b(jpi,jpj)  , vn_b(jpi,jpj)   , & 
     109         &      ua_b(jpi,jpj)  , va_b(jpi,jpj)   , & 
     110         &      un_adv(jpi,jpj), vn_adv(jpi,jpj) , & 
     111         &      sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 
     112         &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , & 
     113         &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , & 
     114         &      wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), &  
     115         &      zwz(jpi,jpj), STAT= ierr(1) ) 
     116 
     117      IF( ln_bt_fw )      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), STAT=ierr(2) ) 
     118 
     119      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     120                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
     121 
     122      dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
     123 
    79124      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc ) 
    80125      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') 
     
    82127   END FUNCTION dyn_spg_ts_alloc 
    83128 
    84  
    85129   SUBROUTINE dyn_spg_ts( kt ) 
    86130      !!---------------------------------------------------------------------- 
    87       !!                  ***  routine dyn_spg_ts  *** 
    88131      !! 
    89       !! ** Purpose :   Compute the now trend due to the surface pressure 
    90       !!      gradient in case of free surface formulation with time-splitting. 
    91       !!      Add it to the general trend of momentum equation. 
     132      !! ** Purpose :    
     133      !!      -Compute the now trend due to the explicit time stepping 
     134      !!      of the quasi-linear barotropic system. Barotropic variables are 
     135      !!      advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F) 
     136      !!      or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a 
     137      !!      generalized forward-backward (see ref. below) time stepping. 
     138      !!      -Update the free surface at step "n+1" (ssha, sshu_a, sshv_a). 
     139      !!      -Compute barotropic advective velocities at step "n" to be used  
     140      !!      to advect tracers latter on. These are compliant with discrete 
     141      !!      continuity equation taken at the baroclinic time steps, thus  
     142      !!      ensuring tracers conservation. 
    92143      !! 
    93       !! ** Method  :   Free surface formulation with time-splitting 
    94       !!      -1- Save the vertically integrated trend. This general trend is 
    95       !!          held constant over the barotropic integration. 
    96       !!          The Coriolis force is removed from the general trend as the 
    97       !!          surface gradient and the Coriolis force are updated within 
    98       !!          the barotropic integration. 
    99       !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and  
    100       !!          barotropic velocity (ua_e and va_e) through barotropic  
    101       !!          momentum and continuity integration. Barotropic former  
    102       !!          variables are time averaging over the full barotropic cycle 
    103       !!          (= 2 * baroclinic time step) and saved in uX_b  
    104       !!          and vX_b (X specifying after, now or before). 
    105       !!      -3- The new general trend becomes : 
    106       !!          ua = ua - sum_k(ua)/H + ( un_b - ub_b ) 
     144      !! ** Method  :   
    107145      !! 
    108       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
     146      !! ** Action : - Update barotropic velocities: ua_b, va_b 
     147      !!             - Update trend (ua,va) with barotropic component 
     148      !!             - Update ssha, sshu_a, sshv_a 
     149      !!             - Update barotropic advective velocity at kt=now 
    109150      !! 
    110       !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 
     151      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:  
     152      !!              The regional oceanic modeling system (ROMS):  
     153      !!              a split-explicit, free-surface, 
     154      !!              topography-following-coordinate oceanic model.  
     155      !!              Ocean Modelling, 9, 347-404.  
    111156      !!--------------------------------------------------------------------- 
    112157      ! 
    113158      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    114159      ! 
    115       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    116       INTEGER  ::   icycle           ! local scalar 
    117       INTEGER  ::   ikbu, ikbv       ! local scalar 
    118       REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars 
    119       REAL(wp) ::   z1_8, zx1, zy1                            !   -      - 
    120       REAL(wp) ::   z1_4, zx2, zy2                            !   -      - 
    121       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
    122       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
    123       REAL(wp) ::   ua_btm, va_btm                            !   -      - 
    124       ! 
    125       REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv  
    126       REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e  
    127       REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 
     160      LOGICAL  ::   lk_fw_start        ! if true, forward integration  
     161      LOGICAL  ::   lk_init         ! if true, special startup of 2d equations 
     162      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
     163      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     164      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     165      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     166      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      - 
     167      REAL(wp) ::   zu_spg, zv_spg     !   -      - 
     168      REAL(wp) ::   zhura, zhvra          !   -      - 
     169      REAL(wp) ::   za0, za1, za2, za3    !   -      - 
     170      ! 
     171      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     172      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
     173      REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 
     174      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
    128175      !!---------------------------------------------------------------------- 
    129176      ! 
    130177      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
    131178      ! 
    132       CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    133       CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    134       CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    135       ! 
    136       IF( kt == nit000 ) THEN             !* initialisation 
    137          ! 
    138          IF(lwp) WRITE(numout,*) 
    139          IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend' 
    140          IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting' 
    141          IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro 
    142          ! 
    143          CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b   
    144          ! 
    145          ua_e  (:,:) = un_b (:,:) 
    146          va_e  (:,:) = vn_b (:,:) 
    147          hu_e  (:,:) = hu   (:,:) 
    148          hv_e  (:,:) = hv   (:,:) 
    149          hur_e (:,:) = hur  (:,:) 
    150          hvr_e (:,:) = hvr  (:,:) 
    151          IF( ln_dynvor_een ) THEN 
    152             ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
     179      !                                         !* Allocate temporay arrays 
     180      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 
     181      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  ) 
     182      CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 
     183      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
     184      ! 
     185      !                                         !* Local constant initialization 
     186      z1_12 = 1._wp / 12._wp  
     187      z1_8  = 0.125_wp                                    
     188      z1_4  = 0.25_wp 
     189      z1_2  = 0.5_wp      
     190      zraur = 1._wp / rau0 
     191      ! 
     192      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     193        z2dt_bf = rdt 
     194      ELSE 
     195        z2dt_bf = 2.0_wp * rdt 
     196      ENDIF 
     197      z1_2dt_b = 1.0_wp / z2dt_bf  
     198      ! 
     199      lk_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     200      lk_fw_start = .FALSE. 
     201      ! 
     202                                                       ! time offset in steps for bdy data update 
     203      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
     204      ! 
     205      IF( kt == nit000 ) THEN                !* initialisation 
     206         ! 
     207         IF (neuler==0) lk_init=.TRUE. 
     208         ! 
     209         IF (ln_bt_fw.OR.(neuler==0)) THEN 
     210           lk_fw_start=.TRUE. 
     211           noffset = 0 
     212         ELSE 
     213           lk_fw_start=.FALSE. 
     214         ENDIF 
     215         ! 
     216         ! Set averaging weights and cycle length: 
     217         CALL ts_wgt(ln_bt_av, lk_fw_start, icycle, wgtbtp1, wgtbtp2) 
     218         ! 
     219         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )  
     220         ! 
     221      ENDIF 
     222      ! 
     223      ! Set arrays to remove/compute coriolis trend. 
     224      ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 
     225      ! Note that these arrays are also used during barotropic loop. These are however frozen 
     226      ! although they should be updated in variable volume case. Not a big approximation. 
     227      ! To remove this approximation, copy lines below inside barotropic loop 
     228      ! and update depths at T-F points (ht and hf resp.) at each barotropic time step 
     229      ! 
     230      IF ( kt == nit000 .OR. lk_vvl ) THEN 
     231         IF ( ln_dynvor_een ) THEN 
     232            DO jj = 1, jpjm1 
     233               DO ji = 1, jpim1 
     234                  zwz(ji,jj) = ( ht(ji  ,jj+1)*tmask(ji  ,jj+1,1) + &  
     235                        &        ht(ji+1,jj+1)*tmask(ji+1,jj+1,1) + & 
     236                        &        ht(ji  ,jj  )*tmask(ji  ,jj  ,1) + &  
     237                        &        ht(ji+1,jj  )*tmask(ji+1,jj  ,1)   & 
     238                        &      ) * z1_4 
     239                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
     240               END DO 
     241            END DO 
     242            CALL lbc_lnk( zwz, 'F', 1._wp ) 
     243            zwz(:,:) = ff(:,:) * zwz(:,:) 
     244 
     245            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    153246            DO jj = 2, jpj 
    154247               DO ji = fs_2, jpi   ! vector opt. 
    155                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
    156                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
    157                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
    158                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    159                END DO 
    160             END DO 
    161          ENDIF 
    162          ! 
    163       ENDIF 
    164  
    165       !                                                     !* Local constant initialization 
    166       z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step 
    167       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic  
    168                                                                         ! time step (euler timestep) 
    169       z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates 
    170       z1_4     = 0.25_wp         
    171       zraur    = 1._wp / rau0                               ! 1 / volumic mass 
    172       ! 
    173       zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
    174       zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
    175       zv_sld = 0._wp   ;   zv_asp = 0._wp 
    176  
    177       IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
    178         z2dt_bf = rdt 
    179       ELSE 
    180         z2dt_bf = 2.0_wp * rdt 
    181       ENDIF 
    182  
     248                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     249                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     250                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     251                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     252               END DO 
     253            END DO 
     254         ELSE 
     255            zwz(:,:) = 0._wp 
     256            ! JC: TBC. hf should be greater than 0  
     257            DO jj = 1, jpj 
     258               DO ji = 1, jpi 
     259                  IF( hf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / hf(ji,jj) 
     260               END DO 
     261            END DO 
     262            zwz(:,:) = ff(:,:) * zwz(:,:) 
     263         ENDIF 
     264      ENDIF 
     265      ! 
     266      ! If forward start at previous time step, and centered integration,  
     267      ! then update averaging weights: 
     268      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 
     269         lk_fw_start=.FALSE. 
     270         CALL ts_wgt(ln_bt_av, lk_fw_start, icycle, wgtbtp1, wgtbtp2) 
     271      ENDIF 
     272                           
    183273      ! ----------------------------------------------------------------------------- 
    184274      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    185275      ! ----------------------------------------------------------------------------- 
    186276      !       
     277      ! Some vertical sums (at now and before time steps) below could be suppressed  
     278      ! if one swap barotropic arrays somewhere 
     279      ! 
    187280      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    188       !                                   ! -------------------------- 
    189       zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
    190       zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     281      !                                   ! -------------------------------------------------- 
     282      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp 
     283      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp 
    191284      ! 
    192285      DO jk = 1, jpkm1 
     
    198291            DO ji = 1, jpi 
    199292#endif 
    200                !                                                                              ! now trend 
    201                zua(ji,jj) = zua(ji,jj) + fse3u  (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    202                zva(ji,jj) = zva(ji,jj) + fse3v  (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    203                !                                                                              ! now velocity  
    204                zun(ji,jj) = zun(ji,jj) + fse3u  (ji,jj,jk) * un(ji,jj,jk) 
    205                zvn(ji,jj) = zvn(ji,jj) + fse3v  (ji,jj,jk) * vn(ji,jj,jk)                
    206                ! 
    207 #if defined key_vvl 
    208                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk)  
    209                vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk) 
    210 #else 
    211                ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    212                vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    213 #endif 
     293               !        ! now trend:                                                                    
     294               zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     295               zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)          
     296               !        ! now bt transp:                    
     297               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)         
     298               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     299               !  ! before bt transp: 
     300               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
     301               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk) 
    214302            END DO 
    215303         END DO 
    216304      END DO 
    217  
     305      ! 
     306      zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 
     307      zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 
     308      ! 
     309      IF( lk_vvl ) THEN 
     310          ub_b(:,:) = ub_b(:,:) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     311          vb_b(:,:) = vb_b(:,:) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
     312      ELSE 
     313          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
     314          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
     315      ENDIF 
     316      ! 
    218317      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    219       DO jk = 1, jpkm1                    ! -------------------------- 
     318      DO jk = 1, jpkm1                    ! ----------------------------------------------------------- 
    220319         DO jj = 2, jpjm1 
    221320            DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    223                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
     321               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 
     322               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
    224323            END DO 
    225324         END DO 
    226325      END DO 
    227  
    228       !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
    229       !                                   ! ---------------------------==== 
    230       zwx(:,:) = zun(:,:) * e2u(:,:)                   ! now transport  
    231       zwy(:,:) = zvn(:,:) * e1v(:,:) 
     326      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
     327      !                                   ! -------------------------------------------------------- 
     328      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport  
     329      zwy(:,:) = vn_b(:,:) * e1v(:,:) 
    232330      ! 
    233331      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     
    239337               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    240338               ! energy conserving formulation for planetary vorticity term 
    241                zcu(ji,jj) = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 
    242                zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) 
    243             END DO 
    244          END DO 
    245          ! 
    246       ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme 
     339               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     340               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     341            END DO 
     342         END DO 
     343         ! 
     344      ELSEIF ( ln_dynvor_ens ) THEN             ! enstrophy conserving scheme 
    247345         DO jj = 2, jpjm1 
    248346            DO ji = fs_2, fs_jpim1   ! vector opt. 
    249                zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    250                zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    251                zcu(ji,jj)  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) 
    252                zcv(ji,jj)  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) 
    253             END DO 
    254          END DO 
    255          ! 
    256       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme 
     347               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     348                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     349               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     350                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     351               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     352               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     353            END DO 
     354         END DO 
     355         ! 
     356      ELSEIF ( ln_dynvor_een ) THEN             ! enstrophy and energy conserving scheme 
    257357         DO jj = 2, jpjm1 
    258358            DO ji = fs_2, fs_jpim1   ! vector opt. 
    259                zcu(ji,jj) = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    260                   &                                + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    261                zcv(ji,jj) = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    262                   &                                + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    263             END DO 
    264          END DO 
    265          ! 
    266       ENDIF 
    267  
     359               zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     360                &                                      + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     361                &                                      + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
     362                &                                      + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     363               zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
     364                &                                      + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     365                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
     366                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     367            END DO 
     368         END DO 
     369         ! 
     370      ENDIF  
     371      ! 
     372      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities 
     373      vn_b (:,:) = vn_b(:,:) * hvr(:,:)   
    268374      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    269375      !                                   ! ---------------------------------------------------- 
    270       IF( lk_vvl ) THEN                         ! Variable volume : remove both Coriolis and Surface pressure gradient 
     376      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    271377         DO jj = 2, jpjm1  
    272378            DO ji = fs_2, fs_jpim1   ! vector opt. 
    273                zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    & 
    274                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
    275                zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    & 
    276                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
    277             END DO 
    278          END DO 
    279       ENDIF 
    280  
    281       DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
     379               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     380               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     381            END DO 
     382         END DO 
     383      ENDIF 
     384 
     385      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    282386         DO ji = fs_2, fs_jpim1 
    283              zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    284              zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     387             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 
     388             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    285389          END DO 
    286       END DO 
    287  
    288                      
    289       !                                             ! Remove barotropic contribution of bottom friction  
    290       !                                             ! from the barotropic transport trend 
    291       zcoef = -1._wp * z1_2dt_b 
    292  
    293       IF(ln_bfrimp) THEN 
    294       !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
    295       !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
    296         DO jj = 2, jpjm1          
    297            DO ji = fs_2, fs_jpim1 
    298               ikbu = mbku(ji,jj) 
    299               ikbv = mbkv(ji,jj) 
    300               ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
    301               va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
    302  
    303               zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
    304               zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
    305            END DO 
    306         END DO 
    307  
    308       ELSE 
    309  
    310 # if defined key_vectopt_loop 
    311         DO jj = 1, 1 
    312            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    313 # else 
    314         DO jj = 2, jpjm1 
    315            DO ji = 2, jpim1 
    316 # endif 
    317             ! Apply stability criteria for bottom friction 
    318             !RBbug for vvl and external mode we may need to use varying fse3 
    319             !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    320               zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    321               zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    322            END DO 
    323         END DO 
    324  
    325         IF( lk_vvl ) THEN 
    326            DO jj = 2, jpjm1 
    327               DO ji = fs_2, fs_jpim1   ! vector opt. 
    328                  zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    329                     &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 
    330                  zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    331                     &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 
    332               END DO 
    333            END DO 
    334         ELSE 
    335            DO jj = 2, jpjm1 
    336               DO ji = fs_2, fs_jpim1   ! vector opt. 
    337                  zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    338                  zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    339               END DO 
    340            END DO 
    341         ENDIF 
    342       END IF    ! end (ln_bfrimp) 
    343  
    344                      
    345       !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    346       !                                   ! --------------------------  
    347       zua(:,:) = zua(:,:) * hur(:,:) 
    348       zva(:,:) = zva(:,:) * hvr(:,:) 
    349       ! 
    350       IF( lk_vvl ) THEN 
    351          ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    352          vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    353       ELSE 
    354          ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    355          vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    356       ENDIF 
    357  
     390      END DO  
     391      ! 
     392      !                 ! Add bottom stress contribution from baroclinic velocities:       
     393      IF (ln_bt_fw) THEN 
     394         DO jj = 2, jpjm1                           
     395            DO ji = fs_2, fs_jpim1   ! vector opt. 
     396               ikbu = mbku(ji,jj)        
     397               ikbv = mbkv(ji,jj)     
     398               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
     399               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     400            END DO 
     401         END DO 
     402      ELSE 
     403         DO jj = 2, jpjm1 
     404            DO ji = fs_2, fs_jpim1   ! vector opt. 
     405               ikbu = mbku(ji,jj)        
     406               ikbv = mbkv(ji,jj)     
     407               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
     408               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     409            END DO 
     410         END DO 
     411      ENDIF 
     412      ! 
     413      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
     414      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     415      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     416      !        
     417      IF (ln_bt_fw) THEN                        ! Add wind forcing 
     418         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
     419         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 
     420      ELSE 
     421         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 
     422         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
     423      ENDIF   
     424      ! 
     425      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing 
     426         IF (ln_bt_fw) THEN 
     427            DO jj = 2, jpjm1               
     428               DO ji = fs_2, fs_jpim1   ! vector opt. 
     429                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 
     430                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj) 
     431                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     432                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     433               END DO 
     434            END DO 
     435         ELSE 
     436            DO jj = 2, jpjm1               
     437               DO ji = fs_2, fs_jpim1   ! vector opt. 
     438                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    & 
     439                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     440                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    & 
     441                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     442                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 
     443                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 
     444               END DO 
     445            END DO 
     446         ENDIF  
     447      ENDIF 
     448      !                                   !* Right-Hand-Side of the barotropic ssh equation 
     449      !                                   ! ----------------------------------------------- 
     450      !                                         ! Surface net water flux and rivers 
     451      IF (ln_bt_fw) THEN 
     452         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 
     453      ELSE 
     454         zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 
     455      ENDIF 
     456#if defined key_asminc 
     457      !                                         ! Include the IAU weighted SSH increment 
     458      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     459         zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:) 
     460      ENDIF 
     461#endif 
     462      ! 
    358463      ! ----------------------------------------------------------------------- 
    359       !  Phase 2 : Integration of the barotropic equations with time splitting 
     464      !  Phase 2 : Integration of the barotropic equations  
    360465      ! ----------------------------------------------------------------------- 
    361466      ! 
    362467      !                                             ! ==================== ! 
    363468      !                                             !    Initialisations   ! 
     469      !                                             ! ==================== !   
     470      ! Initialize barotropic variables:     
     471      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                              
     472         sshn_e (:,:) = sshn (:,:)             
     473         zun_e  (:,:) = un_b (:,:)             
     474         zvn_e  (:,:) = vn_b (:,:) 
     475      ELSE                                ! CENTERED integration: start from BEFORE fields 
     476         sshn_e (:,:) = sshb (:,:) 
     477         zun_e  (:,:) = ub_b (:,:)          
     478         zvn_e  (:,:) = vb_b (:,:) 
     479      ENDIF 
     480      ! 
     481      ! Initialize depths: 
     482      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 
     483         hu_e   (:,:) = hu_0(:,:) + sshu_b(:,:)  
     484         hv_e   (:,:) = hv_0(:,:) + sshv_b(:,:) 
     485         hur_e  (:,:) = umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     486         hvr_e  (:,:) = vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
     487      ELSE 
     488         hu_e   (:,:) = hu   (:,:)        
     489         hv_e   (:,:) = hv   (:,:)  
     490         hur_e  (:,:) = hur  (:,:)     
     491         hvr_e  (:,:) = hvr  (:,:) 
     492      ENDIF 
     493      ! 
     494      IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5: 
     495         zhup2_e (:,:) = hu(:,:) 
     496         zhvp2_e (:,:) = hv(:,:) 
     497      ENDIF 
     498      ! 
     499      ! Initialize sums: 
     500      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     501      va_b  (:,:) = 0._wp 
     502      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level 
     503      zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
     504      zv_sum(:,:) = 0._wp 
    364505      !                                             ! ==================== ! 
    365       icycle = 2  * nn_baro            ! Number of barotropic sub time-step 
    366        
    367       !                                ! Start from NOW field 
    368       hu_e   (:,:) = hu   (:,:)            ! ocean depth at u- and v-points 
    369       hv_e   (:,:) = hv   (:,:)  
    370       hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points 
    371       hvr_e  (:,:) = hvr  (:,:) 
    372 !RBbug     zsshb_e(:,:) = sshn (:,:)   
    373       zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now) 
    374       sshn_e (:,:) = sshn (:,:) 
    375        
    376       zun_e  (:,:) = un_b (:,:)            ! barotropic velocity (external) 
    377       zvn_e  (:,:) = vn_b (:,:) 
    378       zub_e  (:,:) = un_b (:,:) 
    379       zvb_e  (:,:) = vn_b (:,:) 
    380  
    381       zu_sum  (:,:) = un_b (:,:)           ! summation 
    382       zv_sum  (:,:) = vn_b (:,:) 
    383       zssh_sum(:,:) = sshn (:,:) 
    384  
    385 #if defined key_obc 
    386       ! set ssh corrections to 0 
    387       ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    388       IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp 
    389       IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp 
    390       IF( lp_obc_south )   sshfos_b(:,:) = 0._wp 
    391       IF( lp_obc_north )   sshfon_b(:,:) = 0._wp 
    392 #endif 
    393  
    394       !                                             ! ==================== ! 
    395       DO jn = 1, icycle                             !  sub-time-step loop  ! (from NOW to AFTER+1) 
     506      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    396507         !                                          ! ==================== ! 
    397          z2dt_e = 2. * ( rdt / nn_baro ) 
    398          IF( jn == 1 )   z2dt_e = rdt / nn_baro 
    399  
    400508         !                                                !* Update the forcing (BDY and tides) 
    401509         !                                                !  ------------------ 
    402          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    403          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    404          IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 
    405  
    406          !                                                !* after ssh_e 
     510         ! Update only tidal forcing at open boundaries 
     511#if defined key_tide 
     512         IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
     513         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 
     514#endif 
     515         ! 
     516         ! Set extrapolation coefficients for predictor step: 
     517         IF ((jn<3).AND.lk_init) THEN      ! Forward            
     518           za1 = 1._wp                                           
     519           za2 = 0._wp                         
     520           za3 = 0._wp                         
     521         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105  
     522           za1 =  1.781105_wp              ! za1 =   3/2 +   bet 
     523           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet) 
     524           za3 =  0.281105_wp              ! za3 = bet 
     525         ENDIF 
     526 
     527         ! Extrapolate barotropic velocities at step jit+0.5: 
     528         ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
     529         va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     530 
     531         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
     532            !                                             !  ------------------ 
     533            ! Extrapolate Sea Level at step jit+0.5: 
     534            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
     535            ! 
     536            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     537               DO ji = 2, fs_jpim1   ! Vector opt. 
     538                  zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     539                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     540                     &              +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 
     541                  zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     542                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     543                     &              +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     544               END DO 
     545            END DO 
     546            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp ) 
     547            ! 
     548            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)               ! Ocean depth at U- and V-points 
     549            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     550         ENDIF 
     551         !                                                !* after ssh 
    407552         !                                                !  ----------- 
    408          DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
     553         ! One should enforce volume conservation at open boundaries here 
     554         ! considering fluxes below: 
     555         ! 
     556         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
     557         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
     558         DO jj = 2, jpjm1                                  
    409559            DO ji = fs_2, fs_jpim1   ! vector opt. 
    410                zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    411                   &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
    412                   &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
    413                   &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
    414             END DO 
    415          END DO 
    416          ! 
    417 #if defined key_obc 
    418          !                                                     ! OBC : zhdiv must be zero behind the open boundary 
    419 !!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    420          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east 
    421          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west 
    422          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north 
    423          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south 
     560               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   & 
     561                  &             + zwy(ji,jj) - zwy(ji,jj-1)   & 
     562                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     563            END DO 
     564         END DO 
     565         ! 
     566         ! Sum over sub-time-steps to compute advective velocities 
     567         za2 = wgtbtp2(jn) 
     568         zu_sum  (:,:) = zu_sum  (:,:) + za2 * ua_e  (:,:) * zhup2_e  (:,:) 
     569         zv_sum  (:,:) = zv_sum  (:,:) + za2 * va_e  (:,:) * zhvp2_e  (:,:) 
     570         ! 
     571         ! Set next sea level: 
     572         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     573         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
     574 
     575#if defined key_bdy 
     576         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 
     577         IF (lk_bdy) CALL bdy_ssh( ssha_e ) 
    424578#endif 
    425 #if defined key_bdy 
    426          zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)               ! BDY mask 
    427 #endif 
    428          ! 
    429          DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
    430             DO ji = fs_2, fs_jpim1   ! vector opt. 
    431                ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
    432             END DO 
    433          END DO 
    434  
    435          !                                                !* after barotropic velocities (vorticity scheme dependent) 
    436          !                                                !  ---------------------------   
    437          zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport 
    438          zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
     579         !   
     580         ! Sea Surface Height at u-,v-points (vvl case only) 
     581         IF ( lk_vvl ) THEN                                 
     582            DO jj = 2, jpjm1 
     583               DO ji = 2, jpim1      ! NO Vector Opt. 
     584                  sshu_a(ji,jj) = z1_2  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                 & 
     585                     &                                   * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha_e(ji  ,jj) & 
     586                     &                                     + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) ) 
     587                  sshv_a(ji,jj) = z1_2  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                 & 
     588                     &                                   * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha_e(ji,jj  ) & 
     589                     &                                     + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) ) 
     590               END DO 
     591            END DO 
     592            CALL lbc_lnk( sshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( sshv_a, 'V', 1._wp ) 
     593         ENDIF    
     594         !                                  
     595         ! Half-step back interpolation of SSH for surface pressure computation: 
     596         !---------------------------------------------------------------------- 
     597         IF ((jn==1).AND.lk_init) THEN 
     598           za0=1._wp                        ! Forward-backward 
     599           za1=0._wp                            
     600           za2=0._wp 
     601           za3=0._wp 
     602         ELSEIF ((jn==2).AND.lk_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 
     603           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps 
     604           za1=-0.1666666666666_wp          ! za1 = gam 
     605           za2= 0.0833333333333_wp          ! za2 = eps 
     606           za3= 0._wp               
     607         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
     608           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps     
     609           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps  
     610           za2=0.088_wp                     ! za2 = gam 
     611           za3=0.013_wp                     ! za3 = eps 
     612         ENDIF 
     613 
     614         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
     615          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     616 
     617         ! 
     618         ! Compute associated depths at U and V points: 
     619         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN       
     620            !                                         
     621            DO jj = 2, jpjm1                             
     622               DO ji = 2, jpim1 
     623                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     624                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     625                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                  
     626                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     627                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     628                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     629                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1  
     630                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 
     631               END DO 
     632            END DO 
     633         ENDIF 
     634         ! 
     635         ! Add Coriolis trend: 
     636         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 
     637         ! at each time step. We however keep them constant here for optimization. 
     638         ! Recall that zwx and zwy arrays hold fluxes at this stage: 
     639         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5 
     640         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    439641         ! 
    440642         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==! 
    441643            DO jj = 2, jpjm1 
    442644               DO ji = fs_2, fs_jpim1   ! vector opt. 
    443                   ! surface pressure gradient 
    444                   IF( lk_vvl) THEN 
    445                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    446                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    447                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    448                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    449                   ELSE 
    450                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    451                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    452                   ENDIF 
    453                   ! add tidal astronomical forcing 
    454                   IF ( ln_tide_pot .AND. lk_tide ) THEN  
    455                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    456                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    457                   ENDIF 
    458                   ! energy conserving formulation for planetary vorticity term 
    459645                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    460646                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    461647                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    462648                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    463                   zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
    464                   zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    465                   ! after velocities with implicit bottom friction 
    466  
    467                   IF( ln_bfrimp ) THEN      ! implicit bottom friction 
    468                      !   A new method to implement the implicit bottom friction.  
    469                      !   H. Liu 
    470                      !   Sept 2011 
    471                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    472                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    473                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    474                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    475                      ! 
    476                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    477                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    478                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    479                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    480                      ! 
    481                   ELSE 
    482                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    483                       &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    484                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    485                       &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    486                   ENDIF 
     649                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     650                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    487651               END DO 
    488652            END DO 
     
    491655            DO jj = 2, jpjm1 
    492656               DO ji = fs_2, fs_jpim1   ! vector opt. 
    493                    ! surface pressure gradient 
    494                   IF( lk_vvl) THEN 
    495                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    496                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    497                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    498                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    499                   ELSE 
    500                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    501                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    502                   ENDIF 
    503                   ! add tidal astronomical forcing 
    504                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    505                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    506                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    507                   ENDIF 
    508                   ! enstrophy conserving formulation for planetary vorticity term 
    509                   zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    510                   zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    511                   zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
    512                   zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    513                   ! after velocities with implicit bottom friction 
    514                   IF( ln_bfrimp ) THEN 
    515                      !   A new method to implement the implicit bottom friction.  
    516                      !   H. Liu 
    517                      !   Sept 2011 
    518                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    519                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    520                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    521                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    522                      ! 
    523                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    524                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    525                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    526                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    527                      ! 
    528                   ELSE 
    529                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    530                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    531                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    532                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    533                   ENDIF 
     657                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) & 
     658                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     659                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) & 
     660                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
     661                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     662                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    534663               END DO 
    535664            END DO 
     
    538667            DO jj = 2, jpjm1 
    539668               DO ji = fs_2, fs_jpim1   ! vector opt. 
    540                   ! surface pressure gradient 
    541                   IF( lk_vvl) THEN 
    542                      zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    543                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    544                      zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    545                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    546                   ELSE 
    547                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    548                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    549                   ENDIF 
    550                   ! add tidal astronomical forcing 
    551                   IF ( ln_tide_pot .AND. lk_tide ) THEN 
    552                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    553                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
    554                   ENDIF 
    555                   ! energy/enstrophy conserving formulation for planetary vorticity term 
    556                   zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    557                      &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 
    558                   zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    559                      &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    560                   ! after velocities with implicit bottom friction 
    561                   IF( ln_bfrimp ) THEN 
    562                      !   A new method to implement the implicit bottom friction.  
    563                      !   H. Liu 
    564                      !   Sept 2011 
    565                      ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
    566                       &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
    567                       &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
    568                      ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    569                      ! 
    570                      va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
    571                       &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
    572                       &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
    573                      va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
    574                      ! 
    575                   ELSE 
    576                      ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1)   & 
    577                      &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    578                      va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1)   & 
    579                      &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    580                   ENDIF 
     669                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
     670                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
     671                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
     672                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     673                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
     674                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
     675                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
     676                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    581677               END DO 
    582678            END DO 
    583679            !  
    584680         ENDIF 
    585          !                                                     !* domain lateral boundary 
    586          !                                                     !  ----------------------- 
    587  
    588                                                                ! OBC open boundaries 
    589          IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 
    590  
    591                                                                ! BDY open boundaries 
    592 #if defined key_bdy 
    593          pssh => sshn_e 
     681         ! 
     682         ! Add tidal astronomical forcing if defined 
     683         IF ( lk_tide.AND.ln_tide_pot ) THEN 
     684            DO jj = 2, jpjm1 
     685               DO ji = fs_2, fs_jpim1   ! vector opt. 
     686                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     687                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     688                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 
     689                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 
     690               END DO 
     691            END DO 
     692         ENDIF 
     693         ! 
     694         ! Add bottom stresses: 
     695         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
     696         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     697         ! 
     698         ! Surface pressure trend: 
     699         DO jj = 2, jpjm1 
     700            DO ji = fs_2, fs_jpim1   ! vector opt. 
     701               ! Add surface pressure gradient 
     702               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     703               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     704               zwx(ji,jj) = zu_spg 
     705               zwy(ji,jj) = zv_spg 
     706            END DO 
     707         END DO 
     708         ! 
     709         ! Set next velocities: 
     710         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form 
     711            DO jj = 2, jpjm1 
     712               DO ji = fs_2, fs_jpim1   ! vector opt. 
     713                  ua_e(ji,jj) = (                                zun_e(ji,jj)   &  
     714                            &     + rdtbt * (                      zwx(ji,jj)   & 
     715                            &                                 + zu_trd(ji,jj)   & 
     716                            &                                 + zu_frc(ji,jj) ) &  
     717                            &   ) * umask(ji,jj,1) 
     718 
     719                  va_e(ji,jj) = (                                zvn_e(ji,jj)   & 
     720                            &     + rdtbt * (                      zwy(ji,jj)   & 
     721                            &                                 + zv_trd(ji,jj)   & 
     722                            &                                 + zv_frc(ji,jj) ) & 
     723                            &   ) * vmask(ji,jj,1) 
     724               END DO 
     725            END DO 
     726 
     727         ELSE                 ! Flux form 
     728            DO jj = 2, jpjm1 
     729               DO ji = fs_2, fs_jpim1   ! vector opt. 
     730 
     731                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + sshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
     732                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + sshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
     733 
     734                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
     735                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     736                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
     737                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     738                            &   ) * zhura 
     739 
     740                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   & 
     741                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
     742                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
     743                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) & 
     744                            &   ) * zhvra 
     745               END DO 
     746            END DO 
     747         ENDIF 
     748         ! 
     749         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
     750            !                                          !  ----------------------------------------------         
     751            hu_e (:,:) = hu_0(:,:) + sshu_a(:,:) 
     752            hv_e (:,:) = hv_0(:,:) + sshv_a(:,:) 
     753            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
     754            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     755            ! 
     756         ENDIF 
     757         !                                                 !* domain lateral boundary 
     758         !                                                 !  ----------------------- 
     759         ! 
     760         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries  
     761         CALL lbc_lnk( va_e  , 'V', -1._wp ) 
     762 
     763#if defined key_bdy   
     764 
     765         pssh => ssha_e 
    594766         phur => hur_e 
    595767         phvr => hvr_e 
    596768         pu2d => ua_e 
    597769         pv2d => va_e 
    598  
    599          IF( lk_bdy )   CALL bdy_dyn2d( kt )  
     770                                       
     771         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries 
    600772#endif 
    601  
    602          ! 
    603          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
    604          CALL lbc_lnk( va_e  , 'V', -1. ) 
    605          CALL lbc_lnk( ssha_e, 'T',  1. ) 
    606  
    607          zu_sum  (:,:) = zu_sum  (:,:) + ua_e  (:,:)           ! Sum over sub-time-steps 
    608          zv_sum  (:,:) = zv_sum  (:,:) + va_e  (:,:)  
    609          zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)  
    610  
    611          !                                                !* Time filter and swap 
    612          !                                                !  -------------------- 
    613          IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step) 
    614             zsshb_e(:,:) = sshn_e(:,:) 
    615             zub_e  (:,:) = zun_e (:,:) 
    616             zvb_e  (:,:) = zvn_e (:,:) 
    617             sshn_e (:,:) = ssha_e(:,:) 
    618             zun_e  (:,:) = ua_e  (:,:) 
    619             zvn_e  (:,:) = va_e  (:,:) 
    620          ELSE                                                   ! Swap + Filter 
    621             zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 
    622             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e (:,:) 
    623             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e (:,:) 
    624             sshn_e (:,:) = ssha_e(:,:) 
    625             zun_e  (:,:) = ua_e  (:,:) 
    626             zvn_e  (:,:) = va_e  (:,:) 
    627          ENDIF 
    628  
    629          IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
    630             !                                             !  ------------------ 
    631             DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
    632                DO ji = 1, fs_jpim1   ! Vector opt. 
    633                   zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
    634                      &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    635                      &                     +   e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    636                   zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
    637                      &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    638                      &                     +   e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    639                END DO 
    640             END DO 
    641             CALL lbc_lnk( zsshun_e, 'U', 1. )                   ! lateral boundaries conditions 
    642             CALL lbc_lnk( zsshvn_e, 'V', 1. )  
    643             ! 
    644             hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
    645             hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
    646             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    647             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
    648             ! 
    649          ENDIF 
     773         !                                             !* Swap 
     774         !                                             !  ---- 
     775         ubb_e  (:,:) = ub_e  (:,:) 
     776         ub_e   (:,:) = zun_e (:,:) 
     777         zun_e  (:,:) = ua_e  (:,:) 
     778         ! 
     779         vbb_e  (:,:) = vb_e  (:,:) 
     780         vb_e   (:,:) = zvn_e (:,:) 
     781         zvn_e  (:,:) = va_e  (:,:) 
     782         ! 
     783         sshbb_e(:,:) = sshb_e(:,:) 
     784         sshb_e (:,:) = sshn_e(:,:) 
     785         sshn_e (:,:) = ssha_e(:,:) 
     786 
     787         !                                             !* Sum over whole bt loop 
     788         !                                             !  ---------------------- 
     789         za1 = wgtbtp1(jn)                                     
     790         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities 
     791            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
     792            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
     793         ELSE                                ! Sum transports 
     794            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
     795            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     796         ENDIF 
     797         !                                   ! Sum sea level 
     798         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
    650799         !                                                 ! ==================== ! 
    651800      END DO                                               !        end loop      ! 
    652801      !                                                    ! ==================== ! 
    653  
    654 #if defined key_obc 
    655       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)     !!gm totally useless ????? 
    656       IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:) 
    657       IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:) 
    658       IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:) 
    659 #endif 
    660  
    661802      ! ----------------------------------------------------------------------------- 
    662803      ! Phase 3. update the general trend with the barotropic trend 
    663804      ! ----------------------------------------------------------------------------- 
    664805      ! 
    665       !                                   !* Time average ==> after barotropic u, v, ssh 
    666       zcoef =  1._wp / ( 2 * nn_baro  + 1 )  
    667       zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    668       zv_sum(:,:) = zcoef * zv_sum  (:,:)  
    669       !  
    670       !                                   !* update the general momentum trend 
    671       DO jk=1,jpkm1 
    672          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
    673          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
     806      ! At this stage ssha holds a time averaged value 
     807      !                                                ! Sea Surface Height at u-,v- and f-points 
     808      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     809         DO jj = 1, jpjm1 
     810            DO ji = 1, jpim1      ! NO Vector Opt. 
     811               sshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     812                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
     813                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     814               sshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     815                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
     816                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     817            END DO 
     818         END DO 
     819         CALL lbc_lnk( sshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( sshv_a, 'V', 1._wp ) ! Boundary conditions 
     820      ENDIF 
     821      ! 
     822      ! Set advection velocity correction: 
     823      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
     824         un_adv(:,:) = zu_sum(:,:)*hur(:,:) 
     825         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 
     826      ELSE 
     827         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 
     828         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 
     829      END IF 
     830 
     831      IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
     832         ub2_b(:,:) = zu_sum(:,:) 
     833         vb2_b(:,:) = zv_sum(:,:) 
     834      ENDIF 
     835      ! 
     836      ! Update barotropic trend: 
     837      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
     838         DO jk=1,jpkm1 
     839            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     840            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
     841         END DO 
     842      ELSE 
     843         hu_e(:,:) = hu_0(:,:) + sshu_b(:,:) 
     844         hv_e(:,:) = hv_0(:,:) + sshv_b(:,:) 
     845         DO jk=1,jpkm1 
     846            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b 
     847            va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b 
     848         END DO 
     849         ! Save barotropic velocities not transport: 
     850         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) 
     851         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     852      ENDIF 
     853      ! 
     854      !                                   !* write time-spliting arrays in the restart 
     855      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' ) 
     856      ! 
     857      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 
     858      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 
     859      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 
     860      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
     861      ! 
     862      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
     863      ! 
     864   END SUBROUTINE dyn_spg_ts 
     865 
     866   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 
     867      !!--------------------------------------------------------------------- 
     868      !!                   ***  ROUTINE ts_wgt  *** 
     869      !! 
     870      !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 
     871      !!---------------------------------------------------------------------- 
     872      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true. 
     873      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true. 
     874      INTEGER, INTENT(inout) :: jpit      ! cycle length     
     875      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights 
     876                                                         zwgt2    ! Secondary weights 
     877       
     878      INTEGER ::  jic, jn, ji                      ! temporary integers 
     879      REAL(wp) :: za1, za2 
     880      !!---------------------------------------------------------------------- 
     881 
     882      zwgt1(:) = 0._wp 
     883      zwgt2(:) = 0._wp 
     884 
     885      ! Set time index when averaged value is requested 
     886      IF (ll_fw) THEN  
     887         jic = nn_baro 
     888      ELSE 
     889         jic = 2 * nn_baro 
     890      ENDIF 
     891 
     892      ! Set primary weights: 
     893      IF (ll_av) THEN 
     894           ! Define simple boxcar window for primary weights  
     895           ! (width = nn_baro, centered around jic)      
     896         SELECT CASE ( nn_bt_flt ) 
     897              CASE( 0 )  ! No averaging 
     898                 zwgt1(jic) = 1._wp 
     899                 jpit = jic 
     900 
     901              CASE( 1 )  ! Boxcar, width = nn_baro 
     902                 DO jn = 1, 3*nn_baro 
     903                    za1 = ABS(float(jn-jic))/float(nn_baro)  
     904                    IF (za1 < 0.5_wp) THEN 
     905                      zwgt1(jn) = 1._wp 
     906                      jpit = jn 
     907                    ENDIF 
     908                 ENDDO 
     909 
     910              CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
     911                 DO jn = 1, 3*nn_baro 
     912                    za1 = ABS(float(jn-jic))/float(nn_baro)  
     913                    IF (za1 < 1._wp) THEN 
     914                      zwgt1(jn) = 1._wp 
     915                      jpit = jn 
     916                    ENDIF 
     917                 ENDDO 
     918              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
     919         END SELECT 
     920 
     921      ELSE ! No time averaging 
     922         zwgt1(jic) = 1._wp 
     923         jpit = jic 
     924      ENDIF 
     925     
     926      ! Set secondary weights 
     927      DO jn = 1, jpit 
     928        DO ji = jn, jpit 
     929             zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
     930        END DO 
    674931      END DO 
    675       un_b  (:,:) =  zu_sum(:,:)  
    676       vn_b  (:,:) =  zv_sum(:,:)  
    677       sshn_b(:,:) = zcoef * zssh_sum(:,:)  
    678       ! 
    679       !                                   !* write time-spliting arrays in the restart 
    680       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    681       ! 
    682       CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
    683       CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
    684       CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    685       ! 
    686       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    687       ! 
    688    END SUBROUTINE dyn_spg_ts 
    689  
     932 
     933      ! Normalize weigths: 
     934      za1 = 1._wp / SUM(zwgt1(1:jpit)) 
     935      za2 = 1._wp / SUM(zwgt2(1:jpit)) 
     936      DO jn = 1, jpit 
     937        zwgt1(jn) = zwgt1(jn) * za1 
     938        zwgt2(jn) = zwgt2(jn) * za2 
     939      END DO 
     940      ! 
     941   END SUBROUTINE ts_wgt 
    690942 
    691943   SUBROUTINE ts_rst( kt, cdrw ) 
     
    698950      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    699951      ! 
    700       INTEGER ::  ji, jk        ! dummy loop indices 
    701952      !!---------------------------------------------------------------------- 
    702953      ! 
    703954      IF( TRIM(cdrw) == 'READ' ) THEN 
    704          IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 
    705             CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! external velocity issued 
    706             CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    707          ELSE 
    708             un_b (:,:) = 0._wp 
    709             vn_b (:,:) = 0._wp 
    710             ! vertical sum 
    711             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    712                DO jk = 1, jpkm1 
    713                   DO ji = 1, jpij 
    714                      un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 
    715                      vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 
    716                   END DO 
    717                END DO 
    718             ELSE                             ! No  vector opt. 
    719                DO jk = 1, jpkm1 
    720                   un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    721                   vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
    722                END DO 
    723             ENDIF 
    724             un_b (:,:) = un_b(:,:) * hur(:,:) 
    725             vn_b (:,:) = vn_b(:,:) * hvr(:,:) 
    726          ENDIF 
    727  
    728          ! Vertically integrated velocity (before) 
    729          IF (neuler/=0) THEN 
    730             ub_b (:,:) = 0._wp 
    731             vb_b (:,:) = 0._wp 
    732  
    733             ! vertical sum 
    734             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    735                DO jk = 1, jpkm1 
    736                   DO ji = 1, jpij 
    737                      ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 
    738                      vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 
    739                   END DO 
    740                END DO 
    741             ELSE                             ! No  vector opt. 
    742                DO jk = 1, jpkm1 
    743                   ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 
    744                   vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 
    745                END DO 
    746             ENDIF 
    747  
    748             IF( lk_vvl ) THEN 
    749                ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    750                vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    751             ELSE 
    752                ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    753                vb_b(:,:) = vb_b(:,:) * hvr(:,:) 
    754             ENDIF 
    755          ELSE                                 ! neuler==0 
    756             ub_b (:,:) = un_b (:,:) 
    757             vb_b (:,:) = vn_b (:,:) 
    758          ENDIF 
    759  
    760          IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 
    761             CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) )   ! filtered ssh 
    762          ELSE 
    763             sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value    
    764          ENDIF  
     955         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )    
     956         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )  
     957         IF (.NOT.ln_bt_av) THEN 
     958            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
     959            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )    
     960            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) ) 
     961            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) )  
     962            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )    
     963            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) ) 
     964         ENDIF 
     965      ! 
    765966      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    766          CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh 
    767          CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop 
    768          CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !  
     967         CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) ) 
     968         CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) ) 
     969         ! 
     970         IF (.NOT.ln_bt_av) THEN 
     971            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )  
     972            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) ) 
     973            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) ) 
     974            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) ) 
     975            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) ) 
     976            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) ) 
     977         ENDIF 
    769978      ENDIF 
    770979      ! 
    771980   END SUBROUTINE ts_rst 
     981 
     982   SUBROUTINE dyn_spg_ts_init( kt ) 
     983      !!--------------------------------------------------------------------- 
     984      !!                   ***  ROUTINE dyn_spg_ts_init  *** 
     985      !! 
     986      !! ** Purpose : Set time splitting options 
     987      !!---------------------------------------------------------------------- 
     988      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     989      ! 
     990      INTEGER  :: ji ,jj 
     991      REAL(wp) :: zxr2, zyr2, zcmax 
     992      REAL(wp), POINTER, DIMENSION(:,:) :: zcu 
     993      !! 
     994!      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 
     995!      &                  nn_baro, rn_bt_cmax, nn_bt_flt 
     996      !!---------------------------------------------------------------------- 
     997!      REWIND( numnam )              !* Namelist namsplit: split-explicit free surface 
     998!      READ  ( numnam, namsplit ) 
     999      !         ! Max courant number for ext. grav. waves 
     1000      ! 
     1001      CALL wrk_alloc( jpi, jpj, zcu ) 
     1002      ! 
     1003      IF (lk_vvl) THEN 
     1004         DO jj = 1, jpj 
     1005            DO ji =1, jpi  
     1006               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
     1007               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
     1008               zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 
     1009            END DO 
     1010         END DO 
     1011      ELSE 
     1012         DO jj = 1, jpj 
     1013            DO ji =1, jpi  
     1014               zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 
     1015               zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 
     1016               zcu(ji,jj) = sqrt( grav*ht(ji,jj)*(zxr2 + zyr2) ) 
     1017            END DO 
     1018         END DO 
     1019      ENDIF 
     1020 
     1021      zcmax = MAXVAL(zcu(:,:)) 
     1022      IF( lk_mpp )   CALL mpp_max( zcmax ) 
     1023 
     1024      ! Estimate number of iterations to satisfy a max courant number=0.8  
     1025      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
     1026       
     1027      rdtbt = rdt / FLOAT(nn_baro) 
     1028      zcmax = zcmax * rdtbt 
     1029                     ! Print results 
     1030      IF(lwp) WRITE(numout,*) 
     1031      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 
     1032      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     1033      IF( ln_bt_nn_auto ) THEN 
     1034         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro ' 
     1035         IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 
     1036      ELSE 
     1037         IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist ' 
     1038      ENDIF 
     1039      IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro 
     1040      IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt 
     1041      IF(lwp) WRITE(numout,*) ' Maximum Courant number is   :', zcmax 
     1042 
     1043      IF(ln_bt_av) THEN 
     1044         IF(lwp) WRITE(numout,*) ' ln_bt_av=.true.  => Time averaging over nn_baro time steps is on ' 
     1045      ELSE 
     1046         IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables ' 
     1047      ENDIF 
     1048      ! 
     1049      ! 
     1050      IF(ln_bt_fw) THEN 
     1051         IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables ' 
     1052      ELSE 
     1053         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables ' 
     1054      ENDIF 
     1055      ! 
     1056      IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 
     1057      SELECT CASE ( nn_bt_flt ) 
     1058           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '      Dirac' 
     1059           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = nn_baro' 
     1060           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '      Boxcar: width = 2*nn_baro'  
     1061           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 
     1062      END SELECT 
     1063      ! 
     1064      IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN 
     1065         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 
     1066      ENDIF 
     1067      IF ( zcmax>0.9_wp ) THEN 
     1068         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )           
     1069      ENDIF 
     1070      ! 
     1071      CALL wrk_dealloc( jpi, jpj, zcu ) 
     1072      ! 
     1073   END SUBROUTINE dyn_spg_ts_init 
    7721074 
    7731075#else 
     
    7751077   !!   Default case :   Empty module   No standart free surface cst volume 
    7761078   !!---------------------------------------------------------------------- 
     1079 
     1080   LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping 
     1081 
    7771082CONTAINS 
     1083 
    7781084   INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function 
    7791085      dyn_spg_ts_alloc = 0 
     
    7821088      INTEGER, INTENT(in) :: kt 
    7831089      WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt 
    784    END SUBROUTINE dyn_spg_ts 
     1090   END SUBROUTINE dyn_spg_ts  
    7851091   SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine 
    7861092      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    7871093      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    7881094      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw 
    789    END SUBROUTINE ts_rst     
     1095   END SUBROUTINE ts_rst   
     1096   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine 
     1097      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     1098      WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt 
     1099   END SUBROUTINE dyn_spg_ts_init 
    7901100#endif 
    7911101    
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r3625 r3970  
    1616   USE oce             ! ocean dynamics and tracers 
    1717   USE dom_oce         ! ocean space and time domain 
     18   USE domvvl          ! variable volume 
    1819   USE sbc_oce         ! surface boundary condition: ocean 
    1920   USE zdf_oce         ! ocean vertical physics 
     
    2425   USE wrk_nemo        ! Memory Allocation 
    2526   USE timing          ! Timing 
     27   USE dynadv          ! dynamics: vector invariant versus flux form 
     28   USE dynspg_oce, ONLY: lk_dynspg_ts 
     29   USE dynspg_ts 
    2630 
    2731   IMPLICIT NONE 
     
    2933 
    3034   PUBLIC   dyn_zdf_imp   ! called by step.F90 
     35 
     36   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise  
    3137 
    3238   !! * Substitutions 
     
    6470      INTEGER  ::   ikbu, ikbv   ! local integers 
    6571      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     72      REAL(wp) ::   ze3ua, ze3va 
    6673      !!---------------------------------------------------------------------- 
    6774 
    6875      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    69       REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv 
    7076      !!---------------------------------------------------------------------- 
    7177      ! 
     
    7379      ! 
    7480      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )  
    75       CALL wrk_alloc( jpi,jpj, zavmu, zavmv )  
    7681      ! 
    7782      IF( kt == nit000 ) THEN 
     
    7984         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    8085         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     86         ! 
     87         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator 
     88         ELSE                ;    r_vvl = 0._wp        
     89         ENDIF 
    8190      ENDIF 
    8291 
     
    94103      IF( ln_bfrimp ) THEN 
    95104# if defined key_vectopt_loop 
    96       DO jj = 1, 1 
    97          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     105         DO jj = 1, 1 
     106            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    98107# else 
    99       DO jj = 2, jpjm1 
    100          DO ji = 2, jpim1 
     108         DO jj = 2, jpjm1 
     109            DO ji = 2, jpim1 
    101110# endif 
    102             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    103             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    104             zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 
    105             zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 
    106             avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
    107             avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    108          END DO 
    109       END DO 
    110       ENDIF 
     111               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
     112               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     113               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
     114               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     115            END DO 
     116         END DO 
     117      ENDIF 
     118 
     119#if defined key_dynspg_ts 
     120      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     121         DO jk = 1, jpkm1 
     122            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     123            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     124         END DO 
     125      ELSE                                            ! applied on thickness weighted velocity 
     126         DO jk = 1, jpkm1 
     127            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
     128               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
     129               &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
     130            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
     131               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
     132               &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     133         END DO 
     134      ENDIF 
     135 
     136      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 
     137         ! remove barotropic velocities: 
     138         DO jk = 1, jpkm1 
     139            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 
     140            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 
     141         ENDDO 
     142         ! Add bottom stress due to barotropic component only: 
     143         DO jj = 2, jpjm1         
     144            DO ji = fs_2, fs_jpim1   ! vector opt. 
     145               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     146               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     147               ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
     148               ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     149               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     150               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
     151            END DO 
     152         END DO 
     153      ENDIF 
     154#endif 
    111155 
    112156      ! 2. Vertical diffusion on u 
     
    119163         DO jj = 2, jpjm1  
    120164            DO ji = fs_2, fs_jpim1   ! vector opt. 
    121                zcoef = - p2dt / fse3u(ji,jj,jk) 
     165               ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
     166               zcoef = - p2dt / ze3ua       
    122167               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    123168               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    124169               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    125170               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    126                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    127             END DO 
    128          END DO 
    129       END DO 
    130       DO jj = 2, jpjm1        ! Surface boudary conditions 
     171               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 
     172            END DO 
     173         END DO 
     174      END DO 
     175      DO jj = 2, jpjm1        ! Surface boundary conditions 
    131176         DO ji = fs_2, fs_jpim1   ! vector opt. 
    132177            zwi(ji,jj,1) = 0._wp 
    133             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     178            zwd(ji,jj,1) = 1 - zws(ji,jj,1) 
    134179         END DO 
    135180      END DO 
     
    160205      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    161206         DO ji = fs_2, fs_jpim1   ! vector opt. 
    162             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    163                &                                                       * r1_rau0 / fse3u(ji,jj,1)       ) 
     207            ze3ua =  ( 1. - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     208#if defined key_dynspg_ts 
     209            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     210               &                                      / ( ze3ua * rau0 )  
     211#else 
     212            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     213               &                                                     / ( fse3u(ji,jj,1) * rau0     ) )  
     214#endif 
    164215         END DO 
    165216      END DO 
     
    167218         DO jj = 2, jpjm1    
    168219            DO ji = fs_2, fs_jpim1   ! vector opt. 
    169                zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
     220#if defined key_dynspg_ts 
     221               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
     222#else 
     223               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 
     224#endif 
    170225               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    171226            END DO 
     
    186241      END DO 
    187242 
     243#if ! defined key_dynspg_ts 
    188244      ! Normalization to obtain the general momentum trend ua 
    189245      DO jk = 1, jpkm1 
     
    194250         END DO 
    195251      END DO 
    196  
     252#endif 
    197253 
    198254      ! 3. Vertical diffusion on v 
     
    205261         DO jj = 2, jpjm1    
    206262            DO ji = fs_2, fs_jpim1   ! vector opt. 
    207                zcoef = -p2dt / fse3v(ji,jj,jk) 
     263               ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point 
     264               zcoef = - p2dt / ze3va 
    208265               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    209266               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    210267               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    211268               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    212                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    213             END DO 
    214          END DO 
    215       END DO 
    216       DO jj = 2, jpjm1        ! Surface boudary conditions 
     269               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 
     270            END DO 
     271         END DO 
     272      END DO 
     273      DO jj = 2, jpjm1        ! Surface boundary conditions 
    217274         DO ji = fs_2, fs_jpim1   ! vector opt. 
    218275            zwi(ji,jj,1) = 0._wp 
    219             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     276            zwd(ji,jj,1) = 1. - zws(ji,jj,1) 
    220277         END DO 
    221278      END DO 
     
    246303      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    247304         DO ji = fs_2, fs_jpim1   ! vector opt. 
    248             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    249                &                                                       * r1_rau0 / fse3v(ji,jj,1)       ) 
     305            ze3va =  ( 1. - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
     306#if defined key_dynspg_ts             
     307            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     308               &                                      / ( ze3va * rau0 )  
     309#else 
     310            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     311               &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
     312#endif 
    250313         END DO 
    251314      END DO 
     
    253316         DO jj = 2, jpjm1 
    254317            DO ji = fs_2, fs_jpim1   ! vector opt. 
    255                zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
     318#if defined key_dynspg_ts 
     319               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
     320#else 
     321               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 
     322#endif 
    256323               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    257324            END DO 
     
    273340 
    274341      ! Normalization to obtain the general momentum trend va 
     342#if ! defined key_dynspg_ts 
    275343      DO jk = 1, jpkm1 
    276344         DO jj = 2, jpjm1    
     
    280348         END DO 
    281349      END DO 
    282  
     350#endif 
     351 
     352      ! J. Chanut: Lines below are useless ? 
    283353      !! restore bottom layer avmu(v)  
    284354      IF( ln_bfrimp ) THEN 
     
    292362            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    293363            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    294             avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 
    295             avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 
     364            avmu(ji,jj,ikbu+1) = 0.e0 
     365            avmv(ji,jj,ikbv+1) = 0.e0 
    296366         END DO 
    297367      END DO 
     
    299369      ! 
    300370      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
    301       CALL wrk_dealloc( jpi,jpj, zavmu, zavmv)  
    302371      ! 
    303372      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r3764 r3970  
    88   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
     10   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2829   USE obc_oce 
    2930   USE bdy_oce 
     31   USE bdy_par          
     32   USE bdydyn2d        ! bdy_ssh routine 
    3033   USE diaar5, ONLY:   lk_diaar5 
    3134   USE iom 
    32    USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff  
     35   USE sbcrnf          ! River runoff  
    3336#if defined key_agrif 
    3437   USE agrif_opa_update 
     
    4043   USE wrk_nemo        ! Memory Allocation 
    4144   USE timing          ! Timing 
     45 
     46#if defined key_dynspg_ts 
     47   ! jchanut: Needed modules for dynamics update: 
     48   USE eosbn2           ! equation of state                (eos_bn2 routine) 
     49   USE zpshde           ! partial step: hor. derivative    (zps_hde routine) 
     50   USE dynadv           ! advection                        (dyn_adv routine) 
     51   USE dynvor           ! vorticity term                   (dyn_vor routine) 
     52   USE dynhpg           ! hydrostatic pressure grad.       (dyn_hpg routine) 
     53   USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine) 
     54   USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine) 
     55   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
     56   USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine) 
     57   USE asminc           ! assimilation increments      (dyn_asm_inc routine) 
     58   USE bdydyn3d         ! (bdy_dyn3d_dmp routine)   
     59   USE dynspg_ts        ! for advective velocities issued from time splitting 
     60   USE zdfbfr           ! Bottom friction 
     61#if defined key_agrif 
     62   USE agrif_opa_sponge ! Momemtum and tracers sponges 
     63#endif 
     64#endif 
    4265 
    4366   IMPLICIT NONE 
     
    79102      ! 
    80103      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     104      INTEGER  ::   indic 
    81105      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
    82106      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv 
     
    146170         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points) 
    147171         hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
     172         ! bg jchanut tschanges 
     173#if defined key_dynspg_ts 
     174         ht(:,:) = ht_0(:,:) + sshn(:,:)              ! now ocean depth (at t- and f-points) 
     175         hf(:,:) = hf_0(:,:) + sshf_n(:,:) 
     176#endif 
     177         ! end jchanut tschanges 
    148178         !                                            ! now masked inverse of the ocean depth (at u- and v-points) 
    149179         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) ) 
     
    180210#endif 
    181211#if defined key_bdy 
    182       ssha(:,:) = ssha(:,:) * bdytmask(:,:) 
    183       CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm) 
    184 #endif 
     212      ! bg jchanut tschanges 
     213      ! These lines are not necessary with time splitting since 
     214      ! boundary condition on sea level is set during ts loop 
     215      IF (lk_bdy) THEN 
     216         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 
     217         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 
     218      ENDIF 
     219#endif 
     220      ! end jchanut tschanges 
    185221#if defined key_asminc 
    186222      !                                                ! Include the IAU weighted SSH increment 
     
    190226      ENDIF 
    191227#endif 
    192       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
    193       IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
    194          DO jj = 1, jpjm1 
    195             DO ji = 1, jpim1      ! NO Vector Opt. 
    196                sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    197                   &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
    198                   &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
    199                sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    200                   &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
    201                   &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    202             END DO 
    203          END DO 
    204          CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions 
    205       ENDIF 
    206  
    207228      !                                           !------------------------------! 
    208229      !                                           !     Now Vertical Velocity    ! 
     
    219240      END DO 
    220241 
     242      ! bg jchanut tschanges 
     243#if defined key_dynspg_ts 
     244      ! In case the time splitting case, update most of all momentum trends here: 
     245      ! Note that the computation of vertical velocity above, hence "after" sea level is necessary  
     246      ! to compute momentum advection for the rhs of barotropic loop: 
     247                               CALL eos    ( tsn, rhd, rhop )                   ! now in situ density for hpg computation 
     248 
     249      IF( ln_zps      )        CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &      ! zps: now hor. derivative 
     250            &                                          rhd, gru , grv  )        ! of t, s, rd at the last ocean level 
     251                               CALL zdf_bfr( kt )         ! bottom friction (if quadratic) 
     252 
     253                               ua(:,:,:) = 0.e0           ! set dynamics trends to zero 
     254                               va(:,:,:) = 0.e0 
     255      IF(  ln_asmiau .AND. & 
     256         & ln_dyninc       )   CALL dyn_asm_inc( kt )     ! apply dynamics assimilation increment 
     257      IF( ln_neptsimp )        CALL dyn_nept_cor( kt )    ! subtract Neptune velocities (simplified) 
     258      IF( lk_bdy           )   CALL bdy_dyn3d_dmp( kt )   ! bdy damping trends 
     259                               CALL dyn_adv( kt )         ! advection (vector or flux form) 
     260                               CALL dyn_vor( kt )         ! vorticity term including Coriolis 
     261                               CALL dyn_ldf( kt )         ! lateral mixing 
     262      IF( ln_neptsimp )        CALL dyn_nept_cor( kt )    ! add Neptune velocities (simplified) 
     263#if defined key_agrif 
     264      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn      ! momentum sponge 
     265#endif 
     266                               CALL dyn_hpg( kt )         ! horizontal gradient of Hydrostatic pressure 
     267                               CALL dyn_spg( kt, indic )  ! surface pressure gradient 
     268 
     269                               ua_bak(:,:,:) = ua(:,:,:)  ! save next velocities (not trends !) 
     270                               va_bak(:,:,:) = va(:,:,:) 
     271 
     272      ! At this stage:  
     273      ! 1) ssha, sshu_a, sshv_a have been updated.  
     274      ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as  
     275      !   "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement  
     276      !   with continuity equation are available. 
     277      ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components. 
     278      ! Below => Update "now" velocities, divergence, then vertical velocity 
     279      ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal 
     280      !     momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied 
     281      !     some issues with UBS with the current method. Uncomment "divup" below to update the divergence. 
     282      ! 
     283      CALL wrk_alloc( jpi,jpj,jpk, z3d ) 
     284      ! 
     285      DO jk = 1, jpkm1  
     286         ! Correct velocities:                             
     287         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
     288         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
     289         ! 
     290         ! Compute divergence: 
     291         DO jj = 2, jpjm1 
     292            DO ji = fs_2, fs_jpim1   ! vector opt. 
     293               z3d(ji,jj,jk) =   & 
     294                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       & 
     295                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    & 
     296                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     297            END DO   
     298         END DO  
     299      END DO 
     300      ! 
     301      IF( ln_rnf      )   CALL sbc_rnf_div( z3d )      ! runoffs (update divergence) 
     302      ! 
     303      ! 
     304      ! Set new vertical velocities: 
     305      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
     306         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
     307         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * z3d(:,:,jk)        &     
     308            &                      -  (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )  & 
     309            &                          * tmask(:,:,jk) * z1_2dt 
     310#if defined key_bdy 
     311         ! JC: line below is purely cosmetic I guess: 
     312         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     313#endif 
     314      END DO 
     315      ! 
     316      CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 
     317 
     318!! bg divup 
     319!!      ! 
     320!!      DO jk = 1, jpkm1  
     321!!         ! Correct velocities:                             
     322!!         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
     323!!         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
     324!!         ! 
     325!!      END DO 
     326!!      ! 
     327!!      ! 
     328!!      CALL div_cur( kt )                               ! Horizontal divergence & Relative vorticity 
     329!!      ! 
     330!!      ! Set new vertical velocities: 
     331!!      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
     332!!         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
     333!!         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)      &   
     334!!            &                      -  (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )  & 
     335!!            &                          * tmask(:,:,jk) * z1_2dt 
     336!!#if defined key_bdy 
     337!!         ! JC: line below is purely cosmetic I guess: 
     338!!         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     339!!#endif 
     340!!      END DO 
     341!! end divup 
     342      ! 
     343      ! 
     344      ! End of time splitting case 
     345#else 
     346      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
     347      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     348         DO jj = 1, jpjm1 
     349            DO ji = 1, jpim1      ! NO Vector Opt. 
     350               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     351                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     & 
     352                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     353               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     354                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
     355                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     356            END DO 
     357         END DO 
     358         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions 
     359      ENDIF 
     360#endif 
    221361      !                                           !------------------------------! 
    222362      !                                           !           outputs            ! 
     
    283423         !                    !--------------------------! 
    284424         ! 
    285          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
     425#if defined key_dynspg_ts 
     426         IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN    !** Euler time-stepping: no filter 
     427#else 
     428         IF ( neuler == 0 .AND. kt == nit000 ) THEN 
     429#endif 
     430            sshb  (:,:) = sshn  (:,:) 
    286431            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now) 
     432            sshu_b(:,:) = sshu_n(:,:) 
     433            sshv_b(:,:) = sshv_n(:,:) 
    287434            sshu_n(:,:) = sshu_a(:,:) 
    288435            sshv_n(:,:) = sshv_a(:,:) 
     
    336483         !                    !--------------------------! 
    337484         ! 
    338          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
     485#if defined key_dynspg_ts 
     486         IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN    !** Euler time-stepping: no filter 
     487#else 
     488         IF( neuler == 0 .AND. kt == nit000 ) THEN 
     489#endif 
     490            sshb(:,:) = sshn(:,:) 
    339491            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
    340492            ! 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

    r3651 r3970  
    44  !! Initialization of tidal forcing 
    55  !! History :  9.0  !  07  (O. Le Galloudec)  Original code 
     6  !!            3.5  !  2013-07   (J. Chanut) Compliant with time splitting changes 
    67  !!================================================================================= 
    78#if defined key_tide 
     
    2627CONTAINS 
    2728 
    28   SUBROUTINE upd_tide (kt,kit) 
     29  SUBROUTINE upd_tide (kt, kit, time_offset) 
    2930    !!---------------------------------------------------------------------- 
    3031    !!                 ***  ROUTINE upd_tide  *** 
     
    3233    !! * Local declarations 
    3334 
    34     INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index 
    35     INTEGER  :: ji,jj,jk 
    36     REAL (wp) :: zramp 
    37     REAL (wp), DIMENSION(nb_harmo) :: zwt  
     35    INTEGER, INTENT( in ) ::   kt            ! ocean time-step index 
     36    INTEGER, INTENT( in ), OPTIONAL  :: kit  ! Barotropic timestep counter (ts option) 
     37    INTEGER, INTENT( in ), OPTIONAL  :: time_offset ! time offset in units of timesteps 
     38    INTEGER                          :: time_add    ! time offset in units of timesteps 
     39    INTEGER  :: ji, jj, jk                    
     40    REAL (wp) :: zramp, z_arg, z_t 
    3841    !............................................................................... 
    3942 
    4043    pot_astro(:,:)=0.e0 
    4144    zramp = 1.e0 
     45    time_add = 0 
    4246 
    43     IF (lk_dynspg_ts) THEN 
    44        zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp))) 
    45        IF (ln_tide_ramp) THEN 
    46           zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.) 
    47        ENDIF 
    48     ELSE 
    49        zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt 
    50        IF (ln_tide_ramp) THEN 
    51           zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.)  
    52        ENDIF   
     47    IF( PRESENT(time_offset) ) THEN 
     48       time_add = time_offset 
     49    ENDIF 
     50          
     51    IF( PRESENT(kit) ) THEN   
     52       z_arg = ((kt-kt_tide) * rdt + (kit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
     53    ELSE                               
     54       z_arg = ((kt-kt_tide)+time_add) * rdt 
    5355    ENDIF 
    5456 
    55     do jk=1,nb_harmo 
    56        do ji=1,jpi 
    57           do jj=1,jpj 
    58              pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
    59           enddo 
    60        enddo 
    61     enddo 
     57    IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     58 
     59    DO jk=1,nb_harmo 
     60       z_t = z_arg * omega_tide(jk) 
     61       DO ji=1,jpi 
     62          DO jj=1,jpj 
     63             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp * amp_pot(ji,jj,jk)*COS(z_t + phi_pot(ji,jj,jk))       
     64          END DO 
     65       END DO 
     66    END DO 
    6267 
    6368  END SUBROUTINE upd_tide 
     
    6873  !!---------------------------------------------------------------------- 
    6974CONTAINS 
    70   SUBROUTINE upd_tide( kt,kit )          ! Empty routine 
    71     INTEGER,INTENT (IN) :: kt, kit 
     75  SUBROUTINE upd_tide( kt)          ! Empty routine 
     76    INTEGER,INTENT (IN) :: kt 
    7277    WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt 
    7378  END SUBROUTINE upd_tide 
     
    7883 
    7984END MODULE updtide 
     85 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r3769 r3970  
    343343      IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation 
    344344      IF( lk_bdy        )   CALL     bdy_dta_init   ! Open boundaries initialisation of external data arrays 
    345       IF( lk_bdy        )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing 
     345      IF( lk_bdy.AND.lk_tide )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing 
    346346 
    347347                            CALL dyn_nept_init  ! simplified form of Neptune effect 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r3625 r3970  
    2222   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s] 
    2323   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
     24   ! bg jchanut tschanges 
     25   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ua_bak   ,  va_bak     !: Saved trends for mod. ts     [m/s2] 
     26   ! end jchanut tschanges 
    2427   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s] 
    2528   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity           [s-1] 
     
    7073      ! 
    7174      ALLOCATE( ub   (jpi,jpj,jpk)      , un   (jpi,jpj,jpk)      , ua(jpi,jpj,jpk)       ,     & 
    72          &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &       
     75         &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &  
     76      ! bg jchanut tschanges 
     77#if defined key_dynspg_ts   
     78      ! These temporary arrays are used to save tendencies computed before the time stepping of tracers. 
     79      ! These could be suppressed if ua and va would not have been used as temporary arrays  
     80      ! during tracers' update 
     81         &      ua_bak(jpi,jpj,jpk)     , va_bak(jpi,jpj,jpk)     ,                             & 
     82#endif 
     83      ! end jchanut tschanges 
    7384         &      wn   (jpi,jpj,jpk)      ,                                                       & 
    7485         &      rotb (jpi,jpj,jpk)      , rotn (jpi,jpj,jpk)      ,                             &    
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3769 r3970  
    108108      ! 
    109109      !  VERTICAL PHYSICS 
     110      ! bg jchanut tschanges 
     111      ! One need bottom friction parameter in ssh_wzv routine with time splitting. 
     112      ! The idea could be to move the call below before ssh_wzv. However, "now" scale factors 
     113      ! at U-V points (which are set thanks to sshu_n, sshv_n) are actually available in sshwzv. 
     114      ! These are needed for log bottom friction... 
     115#if ! defined key_dynspg_ts 
    110116                         CALL zdf_bfr( kstp )         ! bottom friction 
     117#endif 
     118      ! end jchanut tschanges 
    111119 
    112120      !                                               ! Vertical eddy viscosity and diffusivity coefficients 
     
    206214            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    207215 
    208       ELSE                                                  ! centered hpg  (eos then time stepping) 
     216      ELSE    
     217                                               ! centered hpg  (eos then time stepping) 
     218      ! bg jchanut tschanges 
     219#if ! defined key_dynspg_ts 
     220      ! eos already called 
    209221                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation 
    210222         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative 
    211223            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
     224#endif 
     225      ! end jchanut tschanges 
    212226         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    213227                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
     
    217231      ! Dynamics                                    (tsa used as workspace) 
    218232      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     233      ! bg jchanut tschanges 
     234#if defined key_dynspg_ts       
     235! revert to previously computed tendencies: 
     236! (not using ua, va as temporary arrays during tracers' update could avoid that) 
     237                               ua(:,:,:) = ua_bak(:,:,:)             
     238                               va(:,:,:) = va_bak(:,:,:) 
     239                               CALL dyn_bfr( kstp )         ! bottom friction 
     240                               CALL dyn_zdf( kstp )         ! vertical diffusion 
     241#else 
     242      ! end jchanut tschanges 
    219243                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    220244                               va(:,:,:) = 0.e0 
     
    236260                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    237261                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
     262      ! bg jchanut tschanges 
     263#endif 
     264      ! end jchanut tschanges 
    238265                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
    239266 
Note: See TracChangeset for help on using the changeset viewer.