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 10799 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90 – NEMO

Ignore:
Timestamp:
2019-03-25T09:38:58+01:00 (5 years ago)
Author:
davestorkey
Message:

Temporarily revert all changes to do with the time-level swapping so the
branch just does variable renaming of the 3D state variables and changes
the interfaces to the DYN and TRA routines from stp.

File:
1 edited

Legend:

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

    r10795 r10799  
    6464CONTAINS 
    6565 
    66    SUBROUTINE dyn_nxt ( kt, kNm1, kNp1, puu, pvv, pe3t, pe3u, pe3v ) 
     66   SUBROUTINE dyn_nxt ( kt ) 
    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 
    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 
     94      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    9895      ! 
    9996      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    10097      INTEGER  ::   ikt          ! local integers 
    101       REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars 
    102       REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
     98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars 
     99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    103100      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    104       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
     101      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
    105102      !!---------------------------------------------------------------------- 
    106103      ! 
     
    118115         ! Ensure below that barotropic velocities match time splitting estimate 
    119116         ! Compute actual transport and replace it with ts estimate at "after" time step 
    120          zue(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1) 
    121          zve(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1) 
     117         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
     118         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
    122119         DO jk = 2, jpkm1 
    123             zue(:,:) = zue(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk) 
    124             zve(:,:) = zve(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk) 
     120            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     121            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
    125122         END DO 
    126123         DO jk = 1, jpkm1 
    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) 
     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) 
    129126         END DO 
    130127         ! 
     
    135132            ! so that asselin contribution is removed at the same time  
    136133            DO jk = 1, jpkm1 
    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) 
     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) 
    139136            END DO   
    140137         ENDIF 
     
    147144# endif 
    148145      ! 
    149       CALL lbc_lnk_multi( 'dynnxt', uu(:,:,:,kNp1), 'U', -1., vv(:,:,:,kNp1), 'V', -1. )     !* local domain boundaries 
     146      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries 
    150147      ! 
    151148      !                                !* BDY open boundaries 
     
    160157         ! 
    161158         !                                  ! Kinetic energy and Conversion 
    162          IF( ln_KE_trd  )   CALL trd_dyn( uu(:,:,:,kNp1), vv(:,:,:,kNp1), jpdyn_ken, kt ) 
     159         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt ) 
    163160         ! 
    164161         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
    165             zua(:,:,:) = ( uu(:,:,:,kNp1) - uu(:,:,:,kNm1) ) * z1_2dt 
    166             zva(:,:,:) = ( vv(:,:,:,kNp1) - vv(:,:,:,kNm1) ) * z1_2dt 
     162            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
     163            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
    167164            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
    168165            CALL iom_put( "vtrd_tot", zva ) 
    169166         ENDIF 
    170167         ! 
    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 
     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 
    173170         !                                  !  computation of the asselin filter trends) 
    174171      ENDIF 
     
    176173      ! Time filter and swap of dynamics arrays 
    177174      ! ------------------------------------------ 
     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 
     181!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
     182            DO jk = 1, jpkm1 
     183!               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
     184!               e3u_b(:,:,jk) = e3u_n(:,:,jk) 
     185!               e3v_b(:,:,jk) = e3v_n(:,:,jk) 
     186               ! 
     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!!gm BUG end 
     192         ENDIF 
     193                                                            !  
    178194          
    179       IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin filter and swap 
     195      ELSE                                             !* Leap-Frog : Asselin filter and swap 
    180196         !                                ! =============! 
    181197         IF( ln_linssh ) THEN             ! Fixed volume ! 
     
    184200               DO jj = 1, jpj 
    185201                  DO ji = 1, jpi     
    186                      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) ) 
    187                      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) ) 
     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) ) 
     204                     ! 
     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) 
    188209                  END DO 
    189210               END DO 
     
    192213         ELSE                             ! Variable volume ! 
    193214            !                             ! ================! 
    194             ! Time-filtered scale factor at t-points 
     215            ! Before scale factor at t-points 
     216            ! (used as a now filtered scale factor until the swap) 
    195217            ! ---------------------------------------------------- 
    196             ALLOCATE( ze3t_f(jpi,jpj,jpk) ) 
    197218            DO jk = 1, jpkm1 
    198                ze3t_f(:,:,jk) = pe3t(:,:,jk) + atfp * ( e3t(:,:,jk,kNm1) - 2._wp * pe3t(:,:,jk) + e3t(:,:,jk,kNp1) ) 
     219               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
    199220            END DO 
    200221            ! Add volume filter correction: compatibility with tracer advection scheme 
     
    202223            zcoef = atfp * rdt * r1_rau0 
    203224 
    204             ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     225            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    205226 
    206227            IF ( ln_rnf ) THEN 
     
    210231                        DO ji = 1, jpi 
    211232                           IF( jk <=  nk_rnf(ji,jj)  ) THEN 
    212                                ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
    213                                       &          * ( pe3t(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
     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) 
    214235                           ENDIF 
    215236                        ENDDO 
     
    217238                  ENDDO 
    218239               ELSE 
    219                   ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
     240                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
    220241               ENDIF 
    221242            END IF 
     
    226247                     DO ji = 1, jpi 
    227248                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN 
    228                            ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    229                                 &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
     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) 
    230251                        ELSEIF ( jk==misfkb(ji,jj) ) THEN 
    231                            ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    232                                 &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
     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) 
    233254                        ENDIF 
    234255                     END DO 
     
    237258            END IF 
    238259            ! 
    239             pe3t(:,:,1:jpkm1) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points 
    240             ! 
    241260            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    242261               ! Before filtered scale factor at (u/v)-points 
    243                CALL dom_vvl_interpol( pe3t(:,:,:), pe3u(:,:,:), 'U' ) 
    244                CALL dom_vvl_interpol( pe3t(:,:,:), pe3v(:,:,:), 'V' ) 
     262               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
     263               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
    245264               DO jk = 1, jpkm1 
    246265                  DO jj = 1, jpj 
    247266                     DO ji = 1, jpi 
    248                         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) ) 
    249                         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) ) 
     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) ) 
     269                        ! 
     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) 
    250274                     END DO 
    251275                  END DO 
     
    255279               ! 
    256280               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    257                ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    258                CALL dom_vvl_interpol( pe3t(:,:,:), ze3u_f, 'U' ) 
    259                CALL dom_vvl_interpol( pe3t(:,:,:), ze3v_f, 'V' ) 
     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' ) 
    260284               DO jk = 1, jpkm1 
    261285                  DO jj = 1, jpj 
    262286                     DO ji = 1, jpi                   
    263                         zue3a = e3u(ji,jj,jk,kNp1) * uu(ji,jj,jk,kNp1) 
    264                         zve3a = e3v(ji,jj,jk,kNp1) * vv(ji,jj,jk,kNp1) 
    265                         zue3n = pe3u(ji,jj,jk)     * puu(ji,jj,jk) 
    266                         zve3n = pe3v(ji,jj,jk)     * pvv(ji,jj,jk) 
    267                         zue3b = e3u(ji,jj,jk,kNm1) * uu(ji,jj,jk,kNm1) 
    268                         zve3b = e3v(ji,jj,jk,kNm1) * vv(ji,jj,jk,kNm1) 
     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) 
    269293                        ! 
    270                         puu(ji,jj,jk) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    271                         pvv(ji,jj,jk) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     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) 
     296                        ! 
     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) 
    272301                     END DO 
    273302                  END DO 
    274303               END DO 
    275                pe3u(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)   
    276                pe3v(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     304               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
     305               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
    277306               ! 
    278307               DEALLOCATE( ze3u_f , ze3v_f ) 
    279308            ENDIF 
    280309            ! 
    281             DEALLOCATE( ze3t_f ) 
    282310         ENDIF 
    283311         ! 
    284312         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    285             ! Revert filtered "now" velocities to time split estimate 
     313            ! Revert "before" velocities to time split estimate 
    286314            ! Doing it here also means that asselin filter contribution is removed   
    287             zue(:,:) = pe3u(:,:,1) * puu(:,:,1) * umask(:,:,1) 
    288             zve(:,:) = pe3v(:,:,1) * pvv(:,:,1) * vmask(:,:,1)     
     315            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
     316            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)     
    289317            DO jk = 2, jpkm1 
    290                zue(:,:) = zue(:,:) + pe3u(:,:,jk) * puu(:,:,jk) * umask(:,:,jk) 
    291                zve(:,:) = zve(:,:) + pe3v(:,:,jk) * pvv(:,:,jk) * vmask(:,:,jk)     
     318               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
     319               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
    292320            END DO 
    293321            DO jk = 1, jpkm1 
    294                puu(:,:,jk) = puu(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 
    295                pvv(:,:,jk) = pvv(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
    296             END DO 
    297          ENDIF 
    298          ! 
    299          ikt = Nnn 
    300       ELSE 
    301          ikt = kNm1 
    302          puu(:,:,:)  =  uu(:,:,:,kNp1) 
    303          pvv(:,:,:)  =  vv(:,:,:,kNp1) 
    304          pe3t(:,:,:) = e3t(:,:,:,kNp1) 
    305          pe3u(:,:,:) = e3u(:,:,:,kNp1) 
    306          pe3v(:,:,:) = e3v(:,:,:,kNp1) 
     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) 
     324            END DO 
     325         ENDIF 
     326         ! 
    307327      ENDIF ! neuler =/0 
    308328      ! 
     
    311331      ! integration 
    312332      ! 
    313       ! DS IMMERSE :: This is very confusing as it stands. But should  
    314       ! hopefully be simpler when we do the time-level swapping for the 
    315       ! external mode solution.  
    316333      ! 
    317334      IF(.NOT.ln_linssh ) THEN 
    318          hu(:,:,ikt) = e3u(:,:,1,ikt ) * umask(:,:,1) 
    319          hv(:,:,ikt) = e3v(:,:,1,ikt ) * vmask(:,:,1) 
     335         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
     336         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
    320337         DO jk = 2, jpkm1 
    321             hu(:,:,ikt) = hu(:,:,ikt) + e3u(:,:,jk,ikt ) * umask(:,:,jk) 
    322             hv(:,:,ikt) = hv(:,:,ikt) + e3v(:,:,jk,ikt ) * vmask(:,:,jk) 
     338            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
     339            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
    323340         END DO 
    324          r1_hu(:,:,ikt) = ssumask(:,:) / ( hu(:,:,ikt) + 1._wp - ssumask(:,:) ) 
    325          r1_hv(:,:,ikt) = ssvmask(:,:) / ( hv(:,:,ikt) + 1._wp - ssvmask(:,:) ) 
    326       ENDIF 
    327       ! 
    328       un_b(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1) 
    329       ub_b(:,:) = e3u(:,:,1,ikt ) * uu(:,:,1,ikt ) * umask(:,:,1) 
    330       vn_b(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1) 
    331       vb_b(:,:) = e3v(:,:,1,ikt ) * vv(:,:,1,ikt ) * vmask(:,:,1) 
     341         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     342         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
     343      ENDIF 
     344      ! 
     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) 
    332349      DO jk = 2, jpkm1 
    333          un_b(:,:) = un_b(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk) 
    334          ub_b(:,:) = ub_b(:,:) + e3u(:,:,jk,ikt ) * uu(:,:,jk,ikt ) * umask(:,:,jk) 
    335          vn_b(:,:) = vn_b(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk) 
    336          vb_b(:,:) = vb_b(:,:) + e3v(:,:,jk,ikt ) * vv(:,:,jk,ikt ) * vmask(:,:,jk) 
     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) 
    337354      END DO 
    338       un_b(:,:) = un_b(:,:) * r1_hu(:,:,kNp1) 
    339       vn_b(:,:) = vn_b(:,:) * r1_hv(:,:,kNp1) 
    340       ub_b(:,:) = ub_b(:,:) * r1_hu(:,:,ikt) 
    341       vb_b(:,:) = vb_b(:,:) * r1_hv(:,:,ikt) 
     355      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 
     356      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 
     357      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 
     358      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 
    342359      ! 
    343360      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents 
     
    346363      ENDIF 
    347364      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
    348          zua(:,:,:) = ( puu(:,:,:) - zua(:,:,:) ) * z1_2dt 
    349          zva(:,:,:) = ( pvv(:,:,:) - zva(:,:,:) ) * z1_2dt 
     365         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
     366         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
    350367         CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 
    351368      ENDIF 
    352369      ! 
    353       IF(ln_ctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,kNp1), clinfo1=' nxt  - uu(:,:,:,kNp1): ', mask1=umask,   & 
    354          &                       tab3d_2=vv(:,:,:,kNp1), clinfo2=' vv(:,:,:,kNp1): '       , mask2=vmask ) 
     370      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
     371         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    355372      !  
    356373      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
Note: See TracChangeset for help on using the changeset viewer.