New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6141 for branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2015-12-21T12:38:26+01:00 (9 years ago)
Author:
acc
Message:

Branch 2015/dev_r5803_NOC_WAD. Merge in dev_merge_2015 changes up to 6136. Conflicts all resolved

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r6138 r6141  
    11MODULE dynspg_ts 
     2   !!====================================================================== 
     3   !!                   ***  MODULE  dynspg_ts  *** 
     4   !! Ocean dynamics:  surface pressure gradient trend, split-explicit scheme 
    25   !!====================================================================== 
    36   !! History :   1.0  ! 2004-12  (L. Bessieres, G. Madec)  Original code 
     
    1316   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
    1417   !!--------------------------------------------------------------------- 
     18 
    1519   !!---------------------------------------------------------------------- 
    16    !!                       split explicit free surface 
    17    !!---------------------------------------------------------------------- 
    18    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    19    !!                 splitting scheme and add to the general trend  
     20   !!   dyn_spg_ts     : compute surface pressure gradient trend using a time-splitting scheme  
     21   !!   dyn_spg_ts_init: initialisation of the time-splitting scheme 
     22   !!   ts_wgt         : set time-splitting weights for temporal averaging (or not) 
     23   !!   ts_rst         : read/write time-splitting fields in restart file 
    2024   !!---------------------------------------------------------------------- 
    2125   USE oce             ! ocean dynamics and tracers 
    2226   USE dom_oce         ! ocean space and time domain 
    2327   USE sbc_oce         ! surface boundary condition: ocean 
     28   USE zdf_oce         ! Bottom friction coefts 
    2429   USE sbcisf          ! ice shelf variable (fwfisf) 
     30   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     31   USE dynadv    , ONLY: ln_dynadv_vec 
    2532   USE phycst          ! physical constants 
    2633   USE dynvor          ! vorticity term 
     34   USE wet_dry         ! wetting/drying flux limter 
    2735   USE bdy_par         ! for lk_bdy 
    2836   USE bdytides        ! open boundary condition data 
     
    3038   USE sbctide         ! tides 
    3139   USE updtide         ! tide potential 
     40   ! 
     41   USE in_out_manager  ! I/O manager 
    3242   USE lib_mpp         ! distributed memory computing library 
    3343   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3444   USE prtctl          ! Print control 
    35    USE in_out_manager  ! I/O manager 
    3645   USE iom             ! IOM library 
    3746   USE restart         ! only for lrst_oce 
    38    USE zdf_oce         ! Bottom friction coefts 
    3947   USE wrk_nemo        ! Memory Allocation 
    4048   USE timing          ! Timing     
    41    USE sbcapr          ! surface boundary condition: atmospheric pressure 
    42    USE wet_dry         ! wetting/drying flux limter 
    43    USE dynadv, ONLY: ln_dynadv_vec 
     49   USE diatmb          ! Top,middle,bottom output 
    4450#if defined key_agrif 
    4551   USE agrif_opa_interp ! agrif 
     
    4955#endif 
    5056 
     57 
    5158   IMPLICIT NONE 
    5259   PRIVATE 
     
    6067   REAL(wp),SAVE :: rdtbt   ! Barotropic time step 
    6168 
    62    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 
    63                     wgtbtp1, &              ! Primary weights used for time filtering of barotropic variables 
    64                     wgtbtp2                 ! Secondary weights used for time filtering of barotropic variables 
    65  
    66    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          ! ff/h at F points 
    67    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   ! triad of coriolis parameter 
    68    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
     69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields 
     70 
     71   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff/h at F points 
     72   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter 
     73   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme) 
     74 
     75   !! Time filtered arrays at baroclinic time step: 
     76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection vel. at "now" barocl. step 
    6977 
    7078   !! * Substitutions 
    71 #  include "domzgr_substitute.h90" 
    7279#  include "vectopt_loop_substitute.h90" 
    7380   !!---------------------------------------------------------------------- 
     
    8289      !!                  ***  routine dyn_spg_ts_alloc  *** 
    8390      !!---------------------------------------------------------------------- 
    84       INTEGER :: ierr(4) 
     91      INTEGER :: ierr(3) 
    8592      !!---------------------------------------------------------------------- 
    8693      ierr(:) = 0 
    87  
    88       ALLOCATE( ssha_e(jpi,jpj),  sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 
    89          &        ua_e(jpi,jpj),    un_e(jpi,jpj),   ub_e(jpi,jpj),   ubb_e(jpi,jpj), & 
    90          &        va_e(jpi,jpj),    vn_e(jpi,jpj),   vb_e(jpi,jpj),   vbb_e(jpi,jpj), & 
    91          &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT= ierr(1) ) 
    92  
    93       ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
    94  
    95       IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    96          &                          ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
    97  
    98       ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 
    99 #if defined key_agrif 
    100          &      ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                              , & 
    101 #endif 
    102          &      STAT= ierr(4)) 
    103  
    104       dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
    105  
     94      ! 
     95      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 
     96      ! 
     97      IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     98         &                            ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
     99         ! 
     100      ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) ) 
     101      ! 
     102      dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 
     103      ! 
    106104      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc ) 
    107105      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 
     
    113111      !!---------------------------------------------------------------------- 
    114112      !! 
    115       !! ** Purpose :    
    116       !!      -Compute the now trend due to the explicit time stepping 
    117       !!      of the quasi-linear barotropic system.  
     113      !! ** Purpose : - Compute the now trend due to the explicit time stepping 
     114      !!              of the quasi-linear barotropic system, and add it to the 
     115      !!              general momentum trend.  
    118116      !! 
    119       !! ** Method  :   
     117      !! ** Method  : - split-explicit schem (time splitting) : 
    120118      !!      Barotropic variables are advanced from internal time steps 
    121119      !!      "n"   to "n+1" if ln_bt_fw=T 
     
    131129      !!      continuity equation taken at the baroclinic time steps. This  
    132130      !!      ensures tracers conservation. 
    133       !!      -Update 3d trend (ua, va) with barotropic component. 
     131      !!      - (ua, va) momentum trend updated with barotropic component. 
    134132      !! 
    135       !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:  
    136       !!              The regional oceanic modeling system (ROMS):  
    137       !!              a split-explicit, free-surface, 
    138       !!              topography-following-coordinate oceanic model.  
    139       !!              Ocean Modelling, 9, 347-404.  
     133      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.  
    140134      !!--------------------------------------------------------------------- 
    141       ! 
    142135      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    143136      ! 
     
    147140      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    148141      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     142      INTEGER  ::   iktu, iktv               ! local integers 
     143      REAL(wp) ::   zmdi 
    149144      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    150145      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     
    167162      ! 
    168163      !                                         !* Allocate temporary arrays 
    169       CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 
    170       CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 
    171       CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 
    172       CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
    173       CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    174       CALL wrk_alloc( jpi, jpj, zhf ) 
     164      CALL wrk_alloc( jpi,jpj,  zsshp2_e, zhdiv ) 
     165      CALL wrk_alloc( jpi,jpj,  zu_trd, zv_trd) 
     166      CALL wrk_alloc( jpi,jpj,  zwx, zwy, zssh_frc, zu_frc, zv_frc) 
     167      CALL wrk_alloc( jpi,jpj,  zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
     168      CALL wrk_alloc( jpi,jpj,  zsshu_a, zsshv_a                                   ) 
     169      CALL wrk_alloc( jpi,jpj,  zhf ) 
    175170      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    176171      ! 
     172      zmdi=1.e+20                               !  missing data indicator for masking 
    177173      !                                         !* Local constant initialization 
    178174      z1_12 = 1._wp / 12._wp  
     
    181177      z1_2  = 0.5_wp      
    182178      zraur = 1._wp / rau0 
    183       ! 
    184       IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
    185         z2dt_bf = rdt 
    186       ELSE 
    187         z2dt_bf = 2.0_wp * rdt 
     179      !                                            ! reciprocal of baroclinic time step  
     180      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
     181      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
    188182      ENDIF 
    189183      z1_2dt_b = 1.0_wp / z2dt_bf  
    190184      ! 
    191       ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     185      ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart  
    192186      ll_fw_start = .FALSE. 
    193       ! 
    194                                                        ! time offset in steps for bdy data update 
    195       IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
     187      !                                            ! time offset in steps for bdy data update 
     188      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro 
     189      ELSE                       ;   noffset =   0  
     190      ENDIF 
    196191      ! 
    197192      IF( kt == nit000 ) THEN                !* initialisation 
     
    202197         IF(lwp) WRITE(numout,*) 
    203198         ! 
    204          IF (neuler==0) ll_init=.TRUE. 
    205          ! 
    206          IF (ln_bt_fw.OR.(neuler==0)) THEN 
    207            ll_fw_start=.TRUE. 
    208            noffset = 0 
     199         IF( neuler == 0 )  ll_init=.TRUE. 
     200         ! 
     201         IF( ln_bt_fw .OR. neuler == 0 ) THEN 
     202            ll_fw_start =.TRUE. 
     203            noffset    = 0 
    209204         ELSE 
    210            ll_fw_start=.FALSE. 
     205            ll_fw_start =.FALSE. 
    211206         ENDIF 
    212207         ! 
    213208         ! Set averaging weights and cycle length: 
    214          CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
    215          ! 
     209         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    216210         ! 
    217211      ENDIF 
     
    224218      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 
    225219      ! 
    226       IF ( kt == nit000 .OR. lk_vvl ) THEN 
    227          IF ( ln_dynvor_een ) THEN              !==  EEN scheme  ==! 
     220      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 
     221         IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==! 
    228222            SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point 
    229223            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    230224               DO jj = 1, jpjm1 
    231225                  DO ji = 1, jpim1 
    232                      zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
    233                         &             ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp   
     226                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    & 
     227                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp   
    234228                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 
    235229                  END DO 
     
    238232               DO jj = 1, jpjm1 
    239233                  DO ji = 1, jpim1 
    240                      zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     & 
    241                         &             ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   & 
     234                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     & 
     235                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   & 
    242236                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
    243237                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     
    281275            DO jk = 1, jpkm1 
    282276               DO jj = 1, jpjm1 
    283                   zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     277                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
    284278               END DO 
    285279            END DO 
     
    313307      ! 
    314308      DO jk = 1, jpkm1 
    315          zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    316          zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)          
     309         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     310         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)          
    317311      END DO 
    318312      ! 
    319       zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 
    320       zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 
     313      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 
     314      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 
    321315      ! 
    322316      ! 
     
    332326      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    333327      !                                   ! -------------------------------------------------------- 
    334       zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes  
    335       zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:) 
     328      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
     329      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 
    336330      ! 
    337331      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme 
     
    378372      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    379373      !                                   ! ---------------------------------------------------- 
    380       IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
     374      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    381375        IF(ln_wd) THEN                          ! Calculating and applying W/D gravity filters 
    382376          wduflt1(:,:) = 1.0_wp 
     
    444438      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    445439         DO ji = fs_2, fs_jpim1 
    446             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 
    447             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    448          END DO 
    449       END DO 
     440             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 
     441             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 
     442          END DO 
     443      END DO  
    450444      ! 
    451445      !                 ! Add bottom stress contribution from baroclinic velocities:       
     
    472466      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    473467      IF(ln_wd) THEN 
    474         zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
    475         zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     468        zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
     469        zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
    476470      ELSE 
    477         zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    478         zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     471        zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
     472        zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 
    479473      END IF 
    480474      ! 
     475      !                                         ! Add top stress contribution from baroclinic velocities:       
     476      IF (ln_bt_fw) THEN 
     477         DO jj = 2, jpjm1 
     478            DO ji = fs_2, fs_jpim1   ! vector opt. 
     479               iktu = miku(ji,jj) 
     480               iktv = mikv(ji,jj) 
     481               zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 
     482               zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 
     483            END DO 
     484         END DO 
     485      ELSE 
     486         DO jj = 2, jpjm1 
     487            DO ji = fs_2, fs_jpim1   ! vector opt. 
     488               iktu = miku(ji,jj) 
     489               iktv = mikv(ji,jj) 
     490               zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 
     491               zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 
     492            END DO 
     493         END DO 
     494      ENDIF 
     495      ! 
     496      ! Note that the "unclipped" top friction parameter is used even with explicit drag 
     497      zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 
     498      zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 
     499      !        
    481500      IF (ln_bt_fw) THEN                        ! Add wind forcing 
    482          zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
    483          zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 
     501         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 
     502         zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 
    484503      ELSE 
    485          zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 
    486          zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
     504         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 
     505         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 
    487506      ENDIF   
    488507      ! 
     
    555574      ! 
    556575      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    557          sshn_e(:,:) = sshn (:,:)             
    558          un_e  (:,:) = un_b (:,:)             
    559          vn_e  (:,:) = vn_b (:,:) 
    560          ! 
    561          hu_e  (:,:) = hu   (:,:)        
    562          hv_e  (:,:) = hv   (:,:)  
    563          hur_e (:,:) = hur  (:,:)     
    564          hvr_e (:,:) = hvr  (:,:) 
     576         sshn_e(:,:) =    sshn(:,:)             
     577         un_e  (:,:) =    un_b(:,:)             
     578         vn_e  (:,:) =    vn_b(:,:) 
     579         ! 
     580         hu_e  (:,:) =    hu_n(:,:)        
     581         hv_e  (:,:) =    hv_n(:,:)  
     582         hur_e (:,:) = r1_hu_n(:,:)     
     583         hvr_e (:,:) = r1_hv_n(:,:) 
    565584      ELSE                                ! CENTRED integration: start from BEFORE fields 
    566          sshn_e(:,:) = sshb (:,:) 
    567          un_e  (:,:) = ub_b (:,:)          
    568          vn_e  (:,:) = vb_b (:,:) 
    569          ! 
    570          hu_e  (:,:) = hu_b (:,:)        
    571          hv_e  (:,:) = hv_b (:,:)  
    572          hur_e (:,:) = hur_b(:,:)     
    573          hvr_e (:,:) = hvr_b(:,:) 
     585         sshn_e(:,:) =    sshb(:,:) 
     586         un_e  (:,:) =    ub_b(:,:)          
     587         vn_e  (:,:) =    vb_b(:,:) 
     588         ! 
     589         hu_e  (:,:) =    hu_b(:,:)        
     590         hv_e  (:,:) =    hv_b(:,:)  
     591         hur_e (:,:) = r1_hu_b(:,:)     
     592         hvr_e (:,:) = r1_hv_b(:,:) 
    574593      ENDIF 
    575594      ! 
     
    589608         ! Update only tidal forcing at open boundaries 
    590609#if defined key_tide 
    591          IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
    592          IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 
     610         IF( lk_bdy      .AND. lk_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 
     611         IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset  ) 
    593612#endif 
    594613         ! 
     
    608627         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
    609628 
    610          IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
     629         IF( .NOT.ln_linssh ) THEN                        !* Update ocean depth (variable volume case only) 
    611630            !                                             !  ------------------ 
    612631            ! Extrapolate Sea Level at step jit+0.5: 
     
    615634            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    616635               DO ji = 2, fs_jpim1   ! Vector opt. 
    617                   zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)     & 
     636                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     & 
    618637                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
    619638                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 
    620                   zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)     & 
     639                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     & 
    621640                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
    622641                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     
    632651            END IF 
    633652         ELSE 
    634             zhup2_e (:,:) = hu(:,:) 
    635             zhvp2_e (:,:) = hv(:,:) 
     653            zhup2_e (:,:) = hu_n(:,:) 
     654            zhvp2_e (:,:) = hv_n(:,:) 
    636655         ENDIF 
    637656         !                                                !* after ssh 
     
    644663         ! 
    645664#if defined key_agrif 
    646          ! Set fluxes during predictor step to ensure  
    647          ! volume conservation 
    648          IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN 
     665         ! Set fluxes during predictor step to ensure volume conservation 
     666         IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 
    649667            IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    650668               DO jj=1,jpj 
     
    683701            END DO 
    684702         END DO 
    685          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     703         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    686704         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
    687705         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    688706 
    689707#if defined key_bdy 
    690          ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 
    691          IF (lk_bdy) CALL bdy_ssh( ssha_e ) 
     708         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 
     709         IF( lk_bdy )  CALL bdy_ssh( ssha_e ) 
    692710#endif 
    693711#if defined key_agrif 
    694          IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 
     712         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn ) 
    695713#endif 
    696714         !   
    697715         ! Sea Surface Height at u-,v-points (vvl case only) 
    698          IF ( lk_vvl ) THEN                                 
     716         IF( .NOT.ln_linssh ) THEN                                 
    699717            DO jj = 2, jpjm1 
    700718               DO ji = 2, jpim1      ! NO Vector Opt. 
    701                   zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)  & 
    702                      &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) & 
    703                      &              +   e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) ) 
    704                   zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)  & 
    705                      &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) & 
    706                      &              +   e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) ) 
     719                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    & 
     720                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
     721                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) ) 
     722                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    & 
     723                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
     724                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) ) 
    707725               END DO 
    708726            END DO 
     
    728746           za3=0.013_wp                     ! za3 = eps 
    729747         ENDIF 
    730  
     748         ! 
    731749         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    732750          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    733  
    734751         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
    735752           wduflt1(:,:) = 1._wp 
     
    774791         ! 
    775792         ! Compute associated depths at U and V points: 
    776          IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN       
     793         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN   !* Vector form 
    777794            !                                         
    778795            DO jj = 2, jpjm1                             
    779796               DO ji = 2, jpim1 
    780                   zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e1e2u(ji  ,jj)    & 
     797                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    & 
    781798                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    & 
    782799                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) ) 
    783                   zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e1e2v(ji  ,jj  )  & 
     800                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  & 
    784801                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  & 
    785802                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) ) 
     
    798815         ! 
    799816         ! Add Coriolis trend: 
    800          ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 
     817         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 
    801818         ! at each time step. We however keep them constant here for optimization. 
    802819         ! Recall that zwx and zwy arrays hold fluxes at this stage: 
     
    804821         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    805822         ! 
    806          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==! 
     823         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==! 
    807824            DO jj = 2, jpjm1 
    808825               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    816833            END DO 
    817834            ! 
    818          ELSEIF ( ln_dynvor_ens ) THEN                    !==  enstrophy conserving scheme  ==! 
     835         ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==! 
    819836            DO jj = 2, jpjm1 
    820837               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    828845            END DO 
    829846            ! 
    830          ELSEIF ( ln_dynvor_een ) THEN !==  energy and enstrophy conserving scheme  ==! 
     847         ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==! 
    831848            DO jj = 2, jpjm1 
    832849               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    859876         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
    860877         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
     878         ! 
     879         ! Add top stresses: 
     880         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 
     881         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
    861882         ! 
    862883         ! Surface pressure trend: 
     
    886907         ! 
    887908         ! Set next velocities: 
    888          IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form 
     909         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   !* Vector form 
    889910            DO jj = 2, jpjm1 
    890911               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    893914                            &                                 + zu_trd(ji,jj)   & 
    894915                            &                                 + zu_frc(ji,jj) ) &  
    895                             &   ) * umask(ji,jj,1) 
     916                            &   ) * ssumask(ji,jj) 
    896917 
    897918                  va_e(ji,jj) = (                                 vn_e(ji,jj)   & 
     
    899920                            &                                 + zv_trd(ji,jj)   & 
    900921                            &                                 + zv_frc(ji,jj) ) & 
    901                             &   ) * vmask(ji,jj,1) 
    902                END DO 
    903             END DO 
    904  
    905          ELSE                 ! Flux form 
     922                            &   ) * ssvmask(ji,jj) 
     923               END DO 
     924            END DO 
     925            ! 
     926         ELSE                                      !* Flux form 
    906927            DO jj = 2, jpjm1 
    907928               DO ji = fs_2, fs_jpim1   ! vector opt. 
     929                  zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 
     930                  zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 
    908931 
    909932                  IF(ln_wd) THEN 
     
    914937                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    915938                  END IF 
    916  
    917939                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 
    918940                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 
     
    921943                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
    922944                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    923                             &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     945                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) & 
    924946                            &   ) * zhura 
    925947 
     
    927949                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
    928950                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
    929                             &               +      hv(ji,jj)  * zv_frc(ji,jj) ) & 
     951                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) & 
    930952                            &   ) * zhvra 
    931953               END DO 
     
    933955         ENDIF 
    934956         ! 
    935          IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
    936             !                                          !  ---------------------------------------------- 
     957         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    937958            IF(ln_wd) THEN 
    938959              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
     
    942963              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    943964            END IF 
    944  
    945             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    946             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     965            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
     966            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
    947967            ! 
    948968         ENDIF 
    949          !                                                 !* domain lateral boundary 
    950          !                                                 !  ----------------------- 
    951          ! 
     969         !                                             !* domain lateral boundary 
    952970         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 
    953  
     971         ! 
    954972#if defined key_bdy   
    955                                                            ! open boundaries 
    956          IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 
     973         !                                                 ! open boundaries 
     974         IF( lk_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 
    957975#endif 
    958976#if defined key_agrif                                                            
     
    976994         !                                             !  ---------------------- 
    977995         za1 = wgtbtp1(jn)                                     
    978          IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities 
     996         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities 
    979997            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
    980998            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
    981          ELSE                                ! Sum transports 
     999         ELSE                                              ! Sum transports 
    9821000            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    9831001            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     
    9951013      zwx(:,:) = un_adv(:,:) 
    9961014      zwy(:,:) = vn_adv(:,:) 
    997       IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
    998          un_adv(:,:) = zwx(:,:)*hur(:,:) 
    999          vn_adv(:,:) = zwy(:,:)*hvr(:,:) 
     1015      IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN      
     1016         un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 
     1017         vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 
    10001018      ELSE 
    1001          un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 
    1002          vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 
     1019         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 
     1020         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 
    10031021      END IF 
    10041022 
    1005       IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
     1023      IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 
    10061024         ub2_b(:,:) = zwx(:,:) 
    10071025         vb2_b(:,:) = zwy(:,:) 
     
    10091027      ! 
    10101028      ! Update barotropic trend: 
    1011       IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
     1029      IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
    10121030         DO jk=1,jpkm1 
    10131031            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     
    10291047         ! 
    10301048         DO jk=1,jpkm1 
    1031             ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 
    1032             va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 
     1049            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 
     1050            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 
    10331051         END DO 
    10341052         ! Save barotropic velocities not transport: 
    1035          ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 
    1036          va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     1053         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     1054         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    10371055      ENDIF 
    10381056      ! 
    10391057      DO jk = 1, jpkm1 
    10401058         ! Correct velocities: 
    1041          un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
    1042          vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
     1059         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 
     1060         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 
    10431061         ! 
    10441062      END DO 
     1063      ! 
     1064      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
     1065      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current 
    10451066      ! 
    10461067#if defined key_agrif 
     
    10481069      ! (used to update coarse grid transports at next time step) 
    10491070      ! 
    1050       IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 
    1051          IF ( Agrif_NbStepint().EQ.0 ) THEN 
    1052             ub2_i_b(:,:) = 0.e0 
    1053             vb2_i_b(:,:) = 0.e0 
     1071      IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 
     1072         IF( Agrif_NbStepint() == 0 ) THEN 
     1073            ub2_i_b(:,:) = 0._wp 
     1074            vb2_i_b(:,:) = 0._wp 
    10541075         END IF 
    10551076         ! 
     
    10581079         vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 
    10591080      ENDIF 
    1060       ! 
    1061       ! 
    10621081#endif       
    1063       ! 
    10641082      !                                   !* write time-spliting arrays in the restart 
    1065       IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' ) 
    1066       ! 
    1067       CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 
    1068       CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 
    1069       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 
    1070       CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
    1071       CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    1072       CALL wrk_dealloc( jpi, jpj, zhf ) 
     1083      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' ) 
     1084      ! 
     1085      CALL wrk_dealloc( jpi,jpj,  zsshp2_e, zhdiv ) 
     1086      CALL wrk_dealloc( jpi,jpj,  zu_trd, zv_trd ) 
     1087      CALL wrk_dealloc( jpi,jpj,  zwx, zwy, zssh_frc, zu_frc, zv_frc ) 
     1088      CALL wrk_dealloc( jpi,jpj,  zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
     1089      CALL wrk_dealloc( jpi,jpj,  zsshu_a, zsshv_a                                   ) 
     1090      CALL wrk_dealloc( jpi,jpj,  zhf ) 
    10731091      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    10741092      ! 
     1093      IF ( ln_diatmb ) THEN 
     1094         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity 
     1095         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity 
     1096      ENDIF 
    10751097      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    10761098      ! 
    10771099   END SUBROUTINE dyn_spg_ts 
     1100 
    10781101 
    10791102   SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 
     
    11541177   END SUBROUTINE ts_wgt 
    11551178 
     1179 
    11561180   SUBROUTINE ts_rst( kt, cdrw ) 
    11571181      !!--------------------------------------------------------------------- 
     
    12071231   END SUBROUTINE ts_rst 
    12081232 
    1209    SUBROUTINE dyn_spg_ts_init( kt ) 
     1233 
     1234   SUBROUTINE dyn_spg_ts_init 
    12101235      !!--------------------------------------------------------------------- 
    12111236      !!                   ***  ROUTINE dyn_spg_ts_init  *** 
     
    12131238      !! ** Purpose : Set time splitting options 
    12141239      !!---------------------------------------------------------------------- 
    1215       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1216       ! 
    1217       INTEGER  :: ji ,jj 
    1218       REAL(wp) :: zxr2, zyr2, zcmax 
    1219       REAL(wp), POINTER, DIMENSION(:,:) :: zcu 
    1220       !! 
     1240      INTEGER  ::   ji ,jj              ! dummy loop indices 
     1241      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar 
     1242      REAL(wp), POINTER, DIMENSION(:,:) ::   zcu 
    12211243      !!---------------------------------------------------------------------- 
    12221244      ! 
    12231245      ! Max courant number for ext. grav. waves 
    12241246      ! 
    1225       CALL wrk_alloc( jpi, jpj, zcu ) 
     1247      CALL wrk_alloc( jpi,jpj,  zcu ) 
    12261248      ! 
    12271249      DO jj = 1, jpj 
     
    12371259 
    12381260      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 
    1239       IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
     1261      IF( ln_bt_auto )  nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
    12401262       
    12411263      rdtbt = rdt / REAL( nn_baro , wp ) 
     
    12671289#if defined key_agrif 
    12681290      ! Restrict the use of Agrif to the forward case only 
    1269       IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
     1291      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )  CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
    12701292#endif 
    12711293      ! 
    12721294      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt 
    12731295      SELECT CASE ( nn_bt_flt ) 
    1274            CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac' 
    1275            CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro' 
    1276            CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro'  
    1277            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 
     1296         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac' 
     1297         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro' 
     1298         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro'  
     1299         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 
    12781300      END SELECT 
    12791301      ! 
     
    12831305      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax 
    12841306      ! 
    1285       IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN 
     1307      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 
    12861308         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 
    12871309      ENDIF 
    1288       IF ( zcmax>0.9_wp ) THEN 
     1310      IF( zcmax>0.9_wp ) THEN 
    12891311         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )           
    12901312      ENDIF 
    12911313      ! 
    1292       CALL wrk_dealloc( jpi, jpj, zcu ) 
     1314      CALL wrk_dealloc( jpi,jpj,  zcu ) 
    12931315      ! 
    12941316   END SUBROUTINE dyn_spg_ts_init 
Note: See TracChangeset for help on using the changeset viewer.