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 3865 for branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 – NEMO

Ignore:
Timestamp:
2013-04-09T18:34:38+02:00 (11 years ago)
Author:
acc
Message:

Branch 2013/dev_r3858_NOC_ZTC, #863. Nearly complete port of 2011/dev_r2739_LOCEAN8_ZTC development branch into v3.5aplha base. Compiles and runs but currently unstable after 8 timesteps with ORCA2_LIM reference configuration.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r3764 r3865  
    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    !!---------------------------------------------------------------------- 
    11  
    12    !!---------------------------------------------------------------------- 
    13    !!   ssh_wzv        : after ssh & now vertical velocity 
    14    !!   ssh_nxt        : filter ans swap the ssh arrays 
     10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
     11   !!---------------------------------------------------------------------- 
     12 
     13   !!---------------------------------------------------------------------- 
     14   !!   ssh_nxt        : after ssh 
     15   !!   ssh_swp        : filter ans swap the ssh arrays 
     16   !!   wzv            : compute now vertical velocity 
    1517   !!---------------------------------------------------------------------- 
    1618   USE oce             ! ocean dynamics and tracers variables 
     
    2022   USE divcur          ! hor. divergence and curl      (div & cur routines) 
    2123   USE iom             ! I/O library 
     24   USE restart         ! only for lrst_oce 
    2225   USE in_out_manager  ! I/O manager 
    2326   USE prtctl          ! Print control 
     
    4447   PRIVATE 
    4548 
    46    PUBLIC   ssh_wzv    ! called by step.F90 
    4749   PUBLIC   ssh_nxt    ! called by step.F90 
     50   PUBLIC   wzv        ! called by step.F90 
     51   PUBLIC   ssh_swp    ! called by step.F90 
    4852 
    4953   !! * Substitutions 
     
    5761CONTAINS 
    5862 
    59    SUBROUTINE ssh_wzv( kt )  
    60       !!---------------------------------------------------------------------- 
    61       !!                ***  ROUTINE ssh_wzv  *** 
     63   SUBROUTINE ssh_nxt( kt ) 
     64      !!---------------------------------------------------------------------- 
     65      !!                ***  ROUTINE ssh_nxt  *** 
    6266      !!                    
    63       !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
    64       !!              and update the now vertical coordinate (lk_vvl=T). 
    65       !! 
    66       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    67       !!      velocity is computed by integrating the horizontal divergence   
    68       !!      from the bottom to the surface minus the scale factor evolution. 
    69       !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     67      !! ** Purpose :   compute the after ssh (ssha) 
     68      !! 
     69      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment 
     70      !!      is computed by integrating the horizontal divergence and multiply by 
     71      !!      by the time step. 
    7072      !! 
    7173      !! ** action  :   ssha    : after sea surface height 
    72       !!                wn      : now vertical velocity 
    73       !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T) 
    74       !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points 
    7574      !! 
    7675      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7776      !!---------------------------------------------------------------------- 
    78       INTEGER, INTENT(in) ::   kt   ! time step 
    79       ! 
    80       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    81       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
    82       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv 
    83       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d 
    84       !!---------------------------------------------------------------------- 
    85       ! 
    86       IF( nn_timing == 1 )  CALL timing_start('ssh_wzv') 
    87       ! 
    88       CALL wrk_alloc( jpi, jpj, z2d, zhdiv )  
     77      ! 
     78      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv 
     79      INTEGER, INTENT(in) ::   kt                      ! time step 
     80      !  
     81      INTEGER             ::   jk                      ! dummy loop indice 
     82      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
     83      !!---------------------------------------------------------------------- 
     84      ! 
     85      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     86      ! 
     87      CALL wrk_alloc( jpi, jpj, zhdiv )  
    8988      ! 
    9089      IF( kt == nit000 ) THEN 
    9190         ! 
    9291         IF(lwp) WRITE(numout,*) 
    93          IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 
     92         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 
    9493         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    9594         ! 
    96          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    97          ! 
    98          IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only) 
    99             DO jj = 1, jpjm1 
    100                DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
    101                   zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
    102                   zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
    103                   zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
    104                   sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    105                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    106                   sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    107                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    108                   sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    109                      &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    110                   sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    111                      &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    112                END DO 
    113             END DO 
    114             CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
    115             CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
    116             DO jj = 1, jpjm1 
    117                DO ji = 1, jpim1      ! NO Vector Opt. 
    118                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
    119                        &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    120                        &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    121                        &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    122                END DO 
    123             END DO 
    124             CALL lbc_lnk( sshf_n, 'F', 1. ) 
    125          ENDIF 
    126          ! 
    127       ENDIF 
    128  
    129       !                                           !------------------------------------------! 
    130       IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case) 
    131          !                                        !------------------------------------------! 
    132          DO jk = 1, jpkm1 
    133             fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays 
    134             fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 
    135             fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 
    136             ! 
    137             fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays 
    138             fse3u (:,:,jk) = fse3u_n (:,:,jk) 
    139             fse3v (:,:,jk) = fse3v_n (:,:,jk) 
    140             fse3f (:,:,jk) = fse3f_n (:,:,jk) 
    141             fse3w (:,:,jk) = fse3w_n (:,:,jk) 
    142             fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 
    143             fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 
    144          END DO 
    145          ! 
    146          hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points) 
    147          hv(:,:) = hv_0(:,:) + sshv_n(:,:) 
    148          !                                            ! now masked inverse of the ocean depth (at u- and v-points) 
    149          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) ) 
    150          hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) ) 
    151          !  
    15295      ENDIF 
    15396      ! 
     
    162105      zhdiv(:,:) = 0._wp 
    163106      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    164         zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 
     107        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 
    165108      END DO 
    166109      !                                                ! Sea surface elevation time stepping 
     
    181124#if defined key_bdy 
    182125      ssha(:,:) = ssha(:,:) * bdytmask(:,:) 
    183       CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm) 
     126      CALL lbc_lnk( ssha, 'T', 1. )  
    184127#endif 
    185128#if defined key_asminc 
     
    190133      ENDIF 
    191134#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) ) 
     135 
     136      !                                           !------------------------------! 
     137      !                                           !           outputs            ! 
     138      !                                           !------------------------------! 
     139      CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
     140      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     141      ! 
     142      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
     143      ! 
     144      CALL wrk_dealloc( jpi, jpj, zhdiv )  
     145      ! 
     146      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
     147      ! 
     148   END SUBROUTINE ssh_nxt 
     149 
     150    
     151   SUBROUTINE wzv( kt ) 
     152      !!---------------------------------------------------------------------- 
     153      !!                ***  ROUTINE wzv  *** 
     154      !!                    
     155      !! ** Purpose :   compute the now vertical velocity 
     156      !! 
     157      !! ** Method  : - Using the incompressibility hypothesis, the vertical  
     158      !!      velocity is computed by integrating the horizontal divergence   
     159      !!      from the bottom to the surface minus the scale factor evolution. 
     160      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     161      !! 
     162      !! ** action  :   wn      : now vertical velocity 
     163      !! 
     164      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     165      !!---------------------------------------------------------------------- 
     166      ! 
     167      INTEGER, INTENT(in) ::   kt           ! time step 
     168      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
     169      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
     170      ! 
     171      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     172      REAL(wp)            ::   z1_2dt       ! local scalars 
     173      !!---------------------------------------------------------------------- 
     174       
     175      IF( nn_timing == 1 )  CALL timing_start('wzv') 
     176      ! 
     177      CALL wrk_alloc( jpi, jpj, z2d )  
     178      CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv )  
     179      ! 
     180      IF( kt == nit000 ) THEN 
     181         ! 
     182         IF(lwp) WRITE(numout,*) 
     183         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
     184         IF(lwp) WRITE(numout,*) '~~~ ' 
     185         ! 
     186         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     187         ! 
     188      ENDIF 
     189      !                                           !------------------------------! 
     190      !                                           !     Now Vertical Velocity    ! 
     191      !                                           !------------------------------! 
     192      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog) 
     193      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
     194      ! 
     195      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
     196         DO jk = 1, jpkm1 
     197            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     198            ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     199            DO jj = 2, jpjm1 
     200               DO ji = fs_2, fs_jpim1   ! vector opt. 
     201                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     202               END DO 
    202203            END DO 
    203204         END DO 
    204          CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions 
    205       ENDIF 
    206  
    207       !                                           !------------------------------! 
    208       !                                           !     Now Vertical Velocity    ! 
    209       !                                           !------------------------------! 
    210       z1_2dt = 1.e0 / z2dt 
    211       DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    212          ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    213          wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    214             &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    215             &                         * tmask(:,:,jk) * z1_2dt 
     205         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     206         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
     207         !                             ! Same question holds for hdivn. Perhaps just for security 
     208         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     209            ! computation of w 
     210            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    & 
     211               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     212         END DO 
     213         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     214      ELSE   ! z_star and linear free surface cases 
     215         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     216            ! computation of w 
     217            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   & 
     218               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     219         END DO 
     220      ENDIF 
     221 
    216222#if defined key_bdy 
    217223         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    218224#endif 
    219       END DO 
    220225 
    221226      !                                           !------------------------------! 
    222227      !                                           !           outputs            ! 
    223228      !                                           !------------------------------! 
    224       CALL iom_put( "woce", wn                    )   ! vertical velocity 
    225       CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
    226       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     229      CALL iom_put( "woce", wn )   ! vertical velocity 
    227230      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    228231         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    229          CALL wrk_alloc( jpi,jpj,jpk, z3d ) 
    230          z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 
     232         z2d(:,:) = rau0 * e12t(:,:) 
    231233         DO jk = 1, jpk 
    232234            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
     
    234236         CALL iom_put( "w_masstr" , z3d                     )   
    235237         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    236          CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 
    237       ENDIF 
    238       ! 
    239       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    240       ! 
    241       CALL wrk_dealloc( jpi, jpj, z2d, zhdiv )  
    242       ! 
    243       IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv') 
    244       ! 
    245    END SUBROUTINE ssh_wzv 
    246  
    247  
    248    SUBROUTINE ssh_nxt( kt ) 
     238      ENDIF 
     239      ! 
     240      CALL wrk_dealloc( jpi, jpj, z2d  )  
     241      CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv )  
     242      ! 
     243      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
     244 
     245 
     246   END SUBROUTINE wzv 
     247 
     248 
     249   SUBROUTINE ssh_swp( kt ) 
    249250      !!---------------------------------------------------------------------- 
    250251      !!                    ***  ROUTINE ssh_nxt  *** 
     
    252253      !! ** Purpose :   achieve the sea surface  height time stepping by  
    253254      !!              applying Asselin time filter and swapping the arrays 
    254       !!              ssha  already computed in ssh_wzv   
     255      !!              ssha  already computed in ssh_nxt   
    255256      !! 
    256257      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
     
    266267      !!---------------------------------------------------------------------- 
    267268      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    268       !! 
    269       INTEGER  ::   ji, jj   ! dummy loop indices 
    270       REAL(wp) ::   zec      ! temporary scalar 
    271       !!---------------------------------------------------------------------- 
    272       ! 
    273       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     269      !!---------------------------------------------------------------------- 
     270      ! 
     271      IF( nn_timing == 1 )  CALL timing_start('ssh_swp') 
    274272      ! 
    275273      IF( kt == nit000 ) THEN 
    276274         IF(lwp) WRITE(numout,*) 
    277          IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 
     275         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
    278276         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    279277      ENDIF 
    280278 
    281       !                       !--------------------------! 
    282       IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points) 
    283          !                    !--------------------------! 
    284          ! 
    285          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    286             sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now) 
    287             sshu_n(:,:) = sshu_a(:,:) 
    288             sshv_n(:,:) = sshv_a(:,:) 
    289             DO jj = 1, jpjm1                                ! ssh now at f-point 
    290                DO ji = 1, jpim1      ! NO Vector Opt. 
    291                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    292                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    293                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    294                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    295                END DO 
    296             END DO 
    297             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    298             ! 
    299          ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
    300             zec = atfp * rdt / rau0 
    301             DO jj = 1, jpj 
    302                DO ji = 1, jpi                               ! before <-- now filtered 
    303                   sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   & 
    304                      &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
    305                   sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after 
    306                   sshu_n(ji,jj) = sshu_a(ji,jj) 
    307                   sshv_n(ji,jj) = sshv_a(ji,jj) 
    308                END DO 
    309             END DO 
    310             DO jj = 1, jpjm1                                ! ssh now at f-point 
    311                DO ji = 1, jpim1      ! NO Vector Opt. 
    312                   sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 & 
    313                      &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
    314                      &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
    315                      &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
    316                END DO 
    317             END DO 
    318             CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions 
    319             ! 
    320             DO jj = 1, jpjm1                                ! ssh before at u- & v-points 
    321                DO ji = 1, jpim1      ! NO Vector Opt. 
    322                   sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
    323                      &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
    324                      &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
    325                   sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
    326                      &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    327                      &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    328                END DO 
    329             END DO 
    330             CALL lbc_lnk( sshu_b, 'U', 1. ) 
    331             CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions 
    332             ! 
    333          ENDIF 
    334          !                    !--------------------------! 
    335       ELSE                    !        fixed levels      !     (ssh at t-point only) 
    336          !                    !--------------------------! 
    337          ! 
    338          IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
    339             sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
    340             ! 
    341          ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap 
    342             DO jj = 1, jpj 
    343                DO ji = 1, jpi                               ! before <-- now filtered 
    344                   sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 
    345                   sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after 
    346                END DO 
    347             END DO 
    348          ENDIF 
    349          ! 
     279      IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter 
     280         sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now) 
     281      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
     282         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered 
     283         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     284         sshn(:,:) = ssha(:,:)                           ! now <-- after 
    350285      ENDIF 
    351286      ! 
     
    357292      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    358293      ! 
    359       IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
    360       ! 
    361    END SUBROUTINE ssh_nxt 
     294      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp') 
     295      ! 
     296   END SUBROUTINE ssh_swp 
    362297 
    363298   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.