Ignore:
Timestamp:
2019-03-12T12:06:29+01:00 (20 months ago)
Author:
davestorkey
Message:

First block of changes:

  1. New state variables uu, vv, e3X and corresponding temporary pointers.
  2. Updated code for time-filtering and time-level-swapping uu, vv, e3X.
  3. DYN routines not updated to use new variables at this stage.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90

    r10425 r10743  
    6464CONTAINS 
    6565 
    66    SUBROUTINE dyn_nxt ( kt ) 
     66   SUBROUTINE dyn_nxt ( kt, kNm1, kNp1, puu, pvv, pe3t, pe3u, pe3v ) 
    6767      !!---------------------------------------------------------------------- 
    6868      !!                  ***  ROUTINE dyn_nxt  *** 
     
    9292      !!               un,vn   now horizontal velocity of next time-step 
    9393      !!---------------------------------------------------------------------- 
    94       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     94      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
     95      INTEGER, INTENT( in ) ::   kNm1, kNp1 ! before and after time level indices 
     96      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: puu, pvv ! now velocities to be time filtered 
     97      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pe3t, pe3u, pe3v ! now scale factors to be time filtered 
    9598      ! 
    9699      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    97100      INTEGER  ::   ikt          ! local integers 
    98       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars 
    99       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     101      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars 
     102      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    100103      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    101       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
     104      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
    102105      !!---------------------------------------------------------------------- 
    103106      ! 
     
    115118         ! Ensure below that barotropic velocities match time splitting estimate 
    116119         ! Compute actual transport and replace it with ts estimate at "after" time step 
    117          zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
    118          zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
     120         zue(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1) 
     121         zve(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1) 
    119122         DO jk = 2, jpkm1 
    120             zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    121             zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     123            zue(:,:) = zue(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk) 
     124            zve(:,:) = zve(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk) 
    122125         END DO 
    123126         DO jk = 1, jpkm1 
    124             ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
    125             va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     127            uu(:,:,jk,kNp1) = ( uu(:,:,jk,kNp1) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
     128            vv(:,:,jk,kNp1) = ( vv(:,:,jk,kNp1) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
    126129         END DO 
    127130         ! 
     
    132135            ! so that asselin contribution is removed at the same time  
    133136            DO jk = 1, jpkm1 
    134                un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 
    135                vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     137               puu(:,:,jk) = ( puu(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 
     138               pvv(:,:,jk) = ( pvv(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
    136139            END DO   
    137140         ENDIF 
     
    144147# endif 
    145148      ! 
    146       CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries 
     149      CALL lbc_lnk_multi( 'dynnxt', uu(:,:,:,kNp1), 'U', -1., vv(:,:,:,kNp1), 'V', -1. )     !* local domain boundaries 
    147150      ! 
    148151      !                                !* BDY open boundaries 
     
    157160         ! 
    158161         !                                  ! Kinetic energy and Conversion 
    159          IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt ) 
     162         IF( ln_KE_trd  )   CALL trd_dyn( uu(:,:,:,kNp1), vv(:,:,:,kNp1), jpdyn_ken, kt ) 
    160163         ! 
    161164         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
    162             zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
    163             zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
     165            zua(:,:,:) = ( uu(:,:,:,kNp1) - uu(:,:,:,kNm1) ) * z1_2dt 
     166            zva(:,:,:) = ( vv(:,:,:,kNp1) - vv(:,:,:,kNm1) ) * z1_2dt 
    164167            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
    165168            CALL iom_put( "vtrd_tot", zva ) 
    166169         ENDIF 
    167170         ! 
    168          zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter 
    169          zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
     171         zua(:,:,:) = puu(:,:,:)             ! save the now velocity before the asselin filter 
     172         zva(:,:,:) = pvv(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
    170173         !                                  !  computation of the asselin filter trends) 
    171174      ENDIF 
     
    173176      ! Time filter and swap of dynamics arrays 
    174177      ! ------------------------------------------ 
    175       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap 
    176          DO jk = 1, jpkm1 
    177             un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua 
    178             vn(:,:,jk) = va(:,:,jk) 
    179          END DO 
    180          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
     178!! TO BE DELETED 
     179!!$      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap 
     180!!$         DO jk = 1, jpkm1 
     181!!$            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua 
     182!!$            vn(:,:,jk) = va(:,:,jk) 
     183!!$         END DO  
     184!!$         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    181185!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
    182             DO jk = 1, jpkm1 
     186!!$            DO jk = 1, jpkm1 
    183187!               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    184188!               e3u_b(:,:,jk) = e3u_n(:,:,jk) 
    185189!               e3v_b(:,:,jk) = e3v_n(:,:,jk) 
    186190               ! 
    187                e3t_n(:,:,jk) = e3t_a(:,:,jk) 
    188                e3u_n(:,:,jk) = e3u_a(:,:,jk) 
    189                e3v_n(:,:,jk) = e3v_a(:,:,jk) 
    190             END DO 
     191!!$               e3t_n(:,:,jk) = e3t_a(:,:,jk) 
     192!!$               e3u_n(:,:,jk) = e3u_a(:,:,jk) 
     193!!$               e3v_n(:,:,jk) = e3v_a(:,:,jk) 
     194!!$            END DO 
    191195!!gm BUG end 
    192          ENDIF 
     196!!$         ENDIF 
     197!! TO BE DELETED 
    193198                                                            !  
    194199          
    195       ELSE                                             !* Leap-Frog : Asselin filter and swap 
     200      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin filter and swap 
    196201         !                                ! =============! 
    197202         IF( ln_linssh ) THEN             ! Fixed volume ! 
     
    200205               DO jj = 1, jpj 
    201206                  DO ji = 1, jpi     
    202                      zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    203                      zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     207                     puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) ) 
     208                     pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) ) 
    204209                     ! 
    205                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    206                      vb(ji,jj,jk) = zvf 
    207                      un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    208                      vn(ji,jj,jk) = va(ji,jj,jk) 
     210!! TO BE DELETED 
     211!!$                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     212!!$                     vb(ji,jj,jk) = zvf 
     213!!$                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     214!!$                     vn(ji,jj,jk) = va(ji,jj,jk) 
     215!! TO BE DELETED 
    209216                  END DO 
    210217               END DO 
     
    213220         ELSE                             ! Variable volume ! 
    214221            !                             ! ================! 
    215             ! Before scale factor at t-points 
    216             ! (used as a now filtered scale factor until the swap) 
     222            ! Time-filtered scale factor at t-points 
    217223            ! ---------------------------------------------------- 
     224            ALLOCATE( ze3t_f(jpi,jpj,jpk) ) 
    218225            DO jk = 1, jpkm1 
    219                e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
     226               ze3t_f(:,:,jk) = pe3t(:,:,jk) + atfp * ( e3t(:,:,jk,kNm1) - 2._wp * pe3t(:,:,jk) + e3t(:,:,jk,kNp1) ) 
    220227            END DO 
    221228            ! Add volume filter correction: compatibility with tracer advection scheme 
     
    223230            zcoef = atfp * rdt * r1_rau0 
    224231 
    225             e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     232            ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    226233 
    227234            IF ( ln_rnf ) THEN 
     
    231238                        DO ji = 1, jpi 
    232239                           IF( jk <=  nk_rnf(ji,jj)  ) THEN 
    233                                e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
    234                                       &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
     240                               ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
     241                                      &          * ( pe3t(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
    235242                           ENDIF 
    236243                        ENDDO 
     
    238245                  ENDDO 
    239246               ELSE 
    240                   e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
     247                  ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
    241248               ENDIF 
    242249            END IF 
     
    247254                     DO ji = 1, jpi 
    248255                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN 
    249                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    250                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
     256                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
     257                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
    251258                        ELSEIF ( jk==misfkb(ji,jj) ) THEN 
    252                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    253                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
     259                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
     260                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
    254261                        ENDIF 
    255262                     END DO 
     
    258265            END IF 
    259266            ! 
     267            pe3t(:,:,1:jpkm1) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points 
     268            ! 
    260269            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    261270               ! Before filtered scale factor at (u/v)-points 
    262                CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
    263                CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
     271               CALL dom_vvl_interpol( pe3t(:,:,:), pe3u(:,:,:), 'U' ) 
     272               CALL dom_vvl_interpol( pe3t(:,:,:), pe3v(:,:,:), 'V' ) 
    264273               DO jk = 1, jpkm1 
    265274                  DO jj = 1, jpj 
    266275                     DO ji = 1, jpi 
    267                         zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    268                         zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
     276                        puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) ) 
     277                        pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) ) 
    269278                        ! 
    270                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    271                         vb(ji,jj,jk) = zvf 
    272                         un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    273                         vn(ji,jj,jk) = va(ji,jj,jk) 
     279!! TO BE DELETED 
     280!!$                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     281!!$                        vb(ji,jj,jk) = zvf 
     282!!$                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     283!!$                        vn(ji,jj,jk) = va(ji,jj,jk) 
     284!! TO BE DELETED 
    274285                     END DO 
    275286                  END DO 
     
    279290               ! 
    280291               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    281                ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    282                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 
    283                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 
     292               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
     293               CALL dom_vvl_interpol( pe3t(:,:,:), ze3u_f, 'U' ) 
     294               CALL dom_vvl_interpol( pe3t(:,:,:), ze3v_f, 'V' ) 
    284295               DO jk = 1, jpkm1 
    285296                  DO jj = 1, jpj 
    286297                     DO ji = 1, jpi                   
    287                         zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 
    288                         zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 
    289                         zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 
    290                         zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    291                         zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 
    292                         zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 
     298                        zue3a = e3u(ji,jj,jk,kNp1) * uu(ji,jj,jk,kNp1) 
     299                        zve3a = e3v(ji,jj,jk,kNp1) * vv(ji,jj,jk,kNp1) 
     300                        zue3n = pe3u(ji,jj,jk)     * puu(ji,jj,jk) 
     301                        zve3n = pe3v(ji,jj,jk)     * pvv(ji,jj,jk) 
     302                        zue3b = e3u(ji,jj,jk,kNm1) * uu(ji,jj,jk,kNm1) 
     303                        zve3b = e3v(ji,jj,jk,kNm1) * vv(ji,jj,jk,kNm1) 
    293304                        ! 
    294                         zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    295                         zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     305                        puu(ji,jj,jk) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
     306                        pvv(ji,jj,jk) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    296307                        ! 
    297                         ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    298                         vb(ji,jj,jk) = zvf 
    299                         un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    300                         vn(ji,jj,jk) = va(ji,jj,jk) 
     308!! TO BE DELETED 
     309!!$                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
     310!!$                        vb(ji,jj,jk) = zvf 
     311!!$                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
     312!!$                        vn(ji,jj,jk) = va(ji,jj,jk) 
     313!! TO BE DELETED 
    301314                     END DO 
    302315                  END DO 
    303316               END DO 
    304                e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
    305                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     317               pe3u(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)   
     318               pe3v(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
    306319               ! 
    307320               DEALLOCATE( ze3u_f , ze3v_f ) 
    308321            ENDIF 
    309322            ! 
     323            DEALLOCATE( ze3t_f ) 
    310324         ENDIF 
    311325         ! 
    312326         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    313             ! Revert "before" velocities to time split estimate 
     327            ! Revert filtered "now" velocities to time split estimate 
    314328            ! Doing it here also means that asselin filter contribution is removed   
    315             zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
    316             zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)     
     329            zue(:,:) = pe3u(:,:,1) * puu(:,:,1) * umask(:,:,1) 
     330            zve(:,:) = pe3v(:,:,1) * pvv(:,:,1) * vmask(:,:,1)     
    317331            DO jk = 2, jpkm1 
    318                zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
    319                zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     332               zue(:,:) = zue(:,:) + pe3u(:,:,jk) * puu(:,:,jk) * umask(:,:,jk) 
     333               zve(:,:) = zve(:,:) + pe3v(:,:,jk) * pvv(:,:,jk) * vmask(:,:,jk)     
    320334            END DO 
    321335            DO jk = 1, jpkm1 
    322                ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 
    323                vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
     336               puu(:,:,jk) = puu(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 
     337               pvv(:,:,jk) = pvv(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
    324338            END DO 
    325339         ENDIF 
     
    331345      ! integration 
    332346      ! 
     347      ! DS IMMERSE :: This is very confusing as it stands. But should  
     348      ! hopefully be simpler when we do the time-level swapping for the 
     349      ! external mode solution.  
    333350      ! 
    334351      IF(.NOT.ln_linssh ) THEN 
    335          hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
    336          hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
     352         hu_b(:,:) = pe3u(:,:,1) * umask(:,:,1) 
     353         hv_b(:,:) = pe3v(:,:,1) * vmask(:,:,1) 
    337354         DO jk = 2, jpkm1 
    338             hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    339             hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
     355            hu_b(:,:) = hu_b(:,:) + pe3u(:,:,jk) * umask(:,:,jk) 
     356            hv_b(:,:) = hv_b(:,:) + pe3v(:,:,jk) * vmask(:,:,jk) 
    340357         END DO 
    341358         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     
    343360      ENDIF 
    344361      ! 
    345       un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 
    346       ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
    347       vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 
    348       vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 
     362      un_b(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1) 
     363      ub_b(:,:) = pe3u(:,:,1)     * puu(:,:,1) * umask(:,:,1) 
     364      vn_b(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1) 
     365      vb_b(:,:) = pe3v(:,:,1)     * pvv(:,:,1) * vmask(:,:,1) 
    349366      DO jk = 2, jpkm1 
    350          un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 
    351          ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
    352          vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 
    353          vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 
     367         un_b(:,:) = un_b(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk) 
     368         ub_b(:,:) = ub_b(:,:) + pe3u(:,:,jk)     * puu(:,:,jk) * umask(:,:,jk) 
     369         vn_b(:,:) = vn_b(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk) 
     370         vb_b(:,:) = vb_b(:,:) + pe3v(:,:,jk)     * pvv(:,:,jk) * vmask(:,:,jk) 
    354371      END DO 
    355372      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 
     
    363380      ENDIF 
    364381      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
    365          zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
    366          zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
     382         zua(:,:,:) = ( puu(:,:,:) - zua(:,:,:) ) * z1_2dt 
     383         zva(:,:,:) = ( pvv(:,:,:) - zva(:,:,:) ) * z1_2dt 
    367384         CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 
    368385      ENDIF 
    369386      ! 
    370       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    371          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
     387      IF(ln_ctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,kNp1), clinfo1=' nxt  - uu(:,:,:,kNp1): ', mask1=umask,   & 
     388         &                       tab3d_2=vv(:,:,:,kNp1), clinfo2=' vv(:,:,:,kNp1): '       , mask2=vmask ) 
    372389      !  
    373390      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
Note: See TracChangeset for help on using the changeset viewer.