Ignore:
Timestamp:
2019-05-24T10:27:58+02:00 (17 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Rewrite dynnxt.F90 and rename dynatf.F90 with time level swapping in stp. This doesn't work yet but committing for the record before I merge in Andrew's changes to the main branch.

File:
1 moved

Legend:

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

    r11049 r11050  
    1 MODULE dynnxt 
     1MODULE dynatf 
    22   !!========================================================================= 
    3    !!                       ***  MODULE  dynnxt  *** 
    4    !! Ocean dynamics: time stepping 
     3   !!                       ***  MODULE  dynatf  *** 
     4   !! Ocean dynamics: time filtering 
    55   !!========================================================================= 
    66   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code 
     
    2020   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
    2121   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification 
     22   !!            4.1  !  2019-05  (D. Storkey) Rename dynnxt -> dynatf. Now just does time filtering. 
    2223   !!------------------------------------------------------------------------- 
    2324   
    24    !!------------------------------------------------------------------------- 
    25    !!   dyn_nxt       : obtain the next (after) horizontal velocity 
    26    !!------------------------------------------------------------------------- 
     25   !!---------------------------------------------------------------------------------------------- 
     26   !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors 
     27   !!---------------------------------------------------------------------------------------------- 
    2728   USE oce            ! ocean dynamics and tracers 
    2829   USE dom_oce        ! ocean space and time domain 
     
    5051#if defined key_agrif 
    5152   USE agrif_oce_interp 
    52 #endif    
     53#endif 
    5354 
    5455   IMPLICIT NONE 
    5556   PRIVATE 
    5657 
    57    PUBLIC    dyn_nxt   ! routine called by step.F90 
     58   PUBLIC    dyn_atf   ! routine called by step.F90 
    5859 
    5960   !!---------------------------------------------------------------------- 
     
    6465CONTAINS 
    6566 
    66    SUBROUTINE dyn_nxt ( kt, Kbb, Kmm, Kaa ) 
     67   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 
    6768      !!---------------------------------------------------------------------- 
    68       !!                  ***  ROUTINE dyn_nxt  *** 
     69      !!                  ***  ROUTINE dyn_atf  *** 
    6970      !!                    
     71      !!  !!!!!!!!!!!!!!!!!!  REWRITE HEADER COMMENTS !!!!!!!!!!!!!!!!! 
     72      !! 
    7073      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
    7174      !!             condition on the after velocity, achieve the time stepping  
     
    8386      !!              * Apply the time filter applied and swap of the dynamics 
    8487      !!             arrays to start the next time step: 
    85       !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 
    86       !!                (un,vn) = (ua,va). 
     88      !!                (puu(:,:,:,Kbb),pvv(:,:,:,Kbb)) = (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) + atfp [ (puu(:,:,:,Kbb),pvv(:,:,:,Kbb)) + (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) - 2 (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) ] 
     89      !!                (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) = (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)). 
    8790      !!             Note that with flux form advection and non linear free surface, 
    8891      !!             the time filter is applied on thickness weighted velocity. 
    89       !!             As a result, dyn_nxt MUST be called after tra_nxt. 
    90       !! 
    91       !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step 
    92       !!               un,vn   now horizontal velocity of next time-step 
     92      !!             As a result, dyn_atf MUST be called after tra_nxt. 
     93      !! 
     94      !! ** Action :   uu(Kmm),vv(Kmm)   filtered now horizontal velocity  
    9395      !!---------------------------------------------------------------------- 
    94       INTEGER, INTENT( in ) ::   kt             ! ocean time-step index 
    95       INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa  ! time level indices 
     96      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
     97      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices 
     98      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered 
     99      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 
    96100      ! 
    97101      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    98       INTEGER  ::   ikt          ! local integers 
    99       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars 
    100       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     102      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars 
     103      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    101104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    102       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
     105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
    103106      !!---------------------------------------------------------------------- 
    104107      ! 
    105       IF( ln_timing    )   CALL timing_start('dyn_nxt') 
     108      IF( ln_timing    )   CALL timing_start('dyn_atf') 
    106109      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     ) 
    107110      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) ) 
     
    109112      IF( kt == nit000 ) THEN 
    110113         IF(lwp) WRITE(numout,*) 
    111          IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping' 
     114         IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering' 
    112115         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    113116      ENDIF 
     
    116119         ! Ensure below that barotropic velocities match time splitting estimate 
    117120         ! Compute actual transport and replace it with ts estimate at "after" time step 
    118          zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
    119          zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
     121         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
     122         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
    120123         DO jk = 2, jpkm1 
    121             zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    122             zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     124            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
     125            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
    123126         END DO 
    124127         DO jk = 1, jpkm1 
    125             ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
    126             va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     128            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 
     129            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
    127130         END DO 
    128131         ! 
     
    133136            ! so that asselin contribution is removed at the same time  
    134137            DO jk = 1, jpkm1 
    135                un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 
    136                vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     138               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
     139               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    137140            END DO   
    138141         ENDIF 
     
    145148# endif 
    146149      ! 
    147       CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries 
     150      CALL lbc_lnk_multi( 'dynnxt', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries 
    148151      ! 
    149152      !                                !* BDY open boundaries 
    150       !! IMMERSE development : Following the general pattern for the new code we want to pass in the  
    151       !!                       velocities to bdy_dyn as arguments so here we use "uu" and "vv" even  
    152       !!                       though we haven't converted the velocity names in the rest of dynnxt.F90 
    153       !!                       because it will be completely rewritten. DS. 
    154       IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa ) 
    155       IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa, dyn3d_only=.true. ) 
     153      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa ) 
     154      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. ) 
    156155 
    157156!!$   Do we need a call to bdy_vol here?? 
     
    162161         ! 
    163162         !                                  ! Kinetic energy and Conversion 
    164          IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt, Kmm ) 
     163         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm ) 
    165164         ! 
    166165         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
    167             zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
    168             zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
     166            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt 
     167            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt 
    169168            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
    170169            CALL iom_put( "vtrd_tot", zva ) 
    171170         ENDIF 
    172171         ! 
    173          zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter 
    174          zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
     172         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter 
     173         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the 
    175174         !                                  !  computation of the asselin filter trends) 
    176175      ENDIF 
     
    178177      ! Time filter and swap of dynamics arrays 
    179178      ! ------------------------------------------ 
    180       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap 
    181          DO jk = 1, jpkm1 
    182             un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua 
    183             vn(:,:,jk) = va(:,:,jk) 
    184          END DO 
    185          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    186 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
    187             DO jk = 1, jpkm1 
    188 !               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    189 !               e3u_b(:,:,jk) = e3u_n(:,:,jk) 
    190 !               e3v_b(:,:,jk) = e3v_n(:,:,jk) 
    191                ! 
    192                e3t_n(:,:,jk) = e3t_a(:,:,jk) 
    193                e3u_n(:,:,jk) = e3u_a(:,:,jk) 
    194                e3v_n(:,:,jk) = e3v_a(:,:,jk) 
    195             END DO 
    196 !!gm BUG end 
    197          ENDIF 
    198                                                             !  
    199179          
    200       ELSE                                             !* Leap-Frog : Asselin filter and swap 
     180      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter  
    201181         !                                ! =============! 
    202182         IF( ln_linssh ) THEN             ! Fixed volume ! 
     
    205185               DO jj = 1, jpj 
    206186                  DO ji = 1, jpi     
    207                      zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    208                      zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
    209                      ! 
    210                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    211                      vb(ji,jj,jk) = zvf 
    212                      un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    213                      vn(ji,jj,jk) = va(ji,jj,jk) 
     187                     puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 
     188                     pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 
    214189                  END DO 
    215190               END DO 
     
    218193         ELSE                             ! Variable volume ! 
    219194            !                             ! ================! 
    220             ! Before scale factor at t-points 
    221             ! (used as a now filtered scale factor until the swap) 
     195            ! Time-filtered scale factor at t-points 
    222196            ! ---------------------------------------------------- 
     197            ALLOCATE( ze3t_f(jpi,jpj,jpk) ) 
    223198            DO jk = 1, jpkm1 
    224                e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
     199               ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) ) 
    225200            END DO 
    226201            ! Add volume filter correction: compatibility with tracer advection scheme 
     
    228203            zcoef = atfp * rdt * r1_rau0 
    229204 
    230             e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
     205            ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    231206 
    232207            IF ( ln_rnf ) THEN 
    233208               IF( ln_rnf_depth ) THEN 
    234                   DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too 
     209                  DO jk = 1, jpkm1 ! Deal with Rivers separately, as can be through depth too 
    235210                     DO jj = 1, jpj 
    236211                        DO ji = 1, jpi 
    237212                           IF( jk <=  nk_rnf(ji,jj)  ) THEN 
    238                                e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
    239                                       &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
     213                               ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
     214                                      &          * ( pe3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
    240215                           ENDIF 
    241216                        ENDDO 
     
    243218                  ENDDO 
    244219               ELSE 
    245                   e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
     220                  ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
    246221               ENDIF 
    247222            END IF 
     
    252227                     DO ji = 1, jpi 
    253228                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN 
    254                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    255                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
     229                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
     230                                &          * ( pe3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
    256231                        ELSEIF ( jk==misfkb(ji,jj) ) THEN 
    257                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    258                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
     232                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
     233                                &          * ( pe3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
    259234                        ENDIF 
    260235                     END DO 
     
    263238            END IF 
    264239            ! 
     240            pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points 
     241            ! 
    265242            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    266243               ! Before filtered scale factor at (u/v)-points 
    267                CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
    268                CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
     244               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 
     245               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 
    269246               DO jk = 1, jpkm1 
    270247                  DO jj = 1, jpj 
    271248                     DO ji = 1, jpi 
    272                         zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    273                         zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 
    274                         ! 
    275                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
    276                         vb(ji,jj,jk) = zvf 
    277                         un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
    278                         vn(ji,jj,jk) = va(ji,jj,jk) 
     249                        puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 
     250                        pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 
    279251                     END DO 
    280252                  END DO 
     
    284256               ! 
    285257               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    286                ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    287                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 
    288                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 
     258               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
     259               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 
     260               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 
    289261               DO jk = 1, jpkm1 
    290262                  DO jj = 1, jpj 
    291263                     DO ji = 1, jpi                   
    292                         zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 
    293                         zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 
    294                         zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 
    295                         zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    296                         zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 
    297                         zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 
     264                        zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 
     265                        zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa) 
     266                        zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm) 
     267                        zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm) 
     268                        zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb) 
     269                        zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) 
    298270                        ! 
    299                         zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    300                         zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    301                         ! 
    302                         ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    303                         vb(ji,jj,jk) = zvf 
    304                         un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    305                         vn(ji,jj,jk) = va(ji,jj,jk) 
     271                        puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
     272                        pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    306273                     END DO 
    307274                  END DO 
    308275               END DO 
    309                e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
    310                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     276               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)   
     277               pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 
    311278               ! 
    312279               DEALLOCATE( ze3u_f , ze3v_f ) 
    313280            ENDIF 
    314281            ! 
     282            DEALLOCATE( ze3t_f ) 
    315283         ENDIF 
    316284         ! 
    317285         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    318             ! Revert "before" velocities to time split estimate 
     286            ! Revert filtered "now" velocities to time split estimate 
    319287            ! Doing it here also means that asselin filter contribution is removed   
    320             zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
    321             zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)     
     288            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
     289            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)     
    322290            DO jk = 2, jpkm1 
    323                zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
    324                zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)     
     291               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
     292               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)     
    325293            END DO 
    326294            DO jk = 1, jpkm1 
    327                ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 
    328                vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 
     295               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk) 
     296               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk) 
    329297            END DO 
    330298         ENDIF 
    331299         ! 
    332       ENDIF ! neuler =/0 
     300      ENDIF ! neuler /= 0 
    333301      ! 
    334302      ! Set "now" and "before" barotropic velocities for next time step: 
     
    336304      ! integration 
    337305      ! 
    338       ! 
    339306      IF(.NOT.ln_linssh ) THEN 
    340          hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
    341          hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
     307         hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1) 
     308         hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1) 
    342309         DO jk = 2, jpkm1 
    343             hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    344             hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
     310            hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk) 
     311            hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk) 
    345312         END DO 
    346          r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
    347          r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    348       ENDIF 
    349       ! 
    350       un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 
    351       ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 
    352       vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 
    353       vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 
     313         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 
     314         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 
     315      ENDIF 
     316      ! 
     317      uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
     318      uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
     319      vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
     320      vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    354321      DO jk = 2, jpkm1 
    355          un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 
    356          ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 
    357          vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 
    358          vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 
     322         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
     323         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
     324         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
     325         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    359326      END DO 
    360       un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 
    361       vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 
    362       ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 
    363       vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 
     327      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa) 
     328      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa) 
     329      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
     330      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 
    364331      ! 
    365332      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents 
    366          CALL iom_put(  "ubar", un_b(:,:) ) 
    367          CALL iom_put(  "vbar", vn_b(:,:) ) 
     333         CALL iom_put(  "ubar", uu_b(:,:,Kmm) ) 
     334         CALL iom_put(  "vbar", vv_b(:,:,Kmm) ) 
    368335      ENDIF 
    369336      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
    370          zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
    371          zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
     337         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt 
     338         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt 
    372339         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm ) 
    373340      ENDIF 
    374341      ! 
    375       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    376          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
     342      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   & 
     343         &                       tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask ) 
    377344      !  
    378345      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
    379346      IF( l_trddyn     )   DEALLOCATE( zua, zva ) 
    380       IF( ln_timing    )   CALL timing_stop('dyn_nxt') 
    381       ! 
    382    END SUBROUTINE dyn_nxt 
     347      IF( ln_timing    )   CALL timing_stop('dyn_atf') 
     348      ! 
     349   END SUBROUTINE dyn_atf 
    383350 
    384351   !!========================================================================= 
    385 END MODULE dynnxt 
     352END MODULE dynatf 
Note: See TracChangeset for help on using the changeset viewer.