Ignore:
Timestamp:
2019-03-25T09:38:58+01:00 (2 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.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
Files:
3 edited

Legend:

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

    r10795 r10799  
    563563 
    564564 
    565    SUBROUTINE dom_vvl_sf_swp( kt, kNm1, kNnn, kNm1_2lev, kNnn_2lev ) 
     565   SUBROUTINE dom_vvl_sf_swp( kt ) 
    566566      !!---------------------------------------------------------------------- 
    567567      !!                ***  ROUTINE dom_vvl_sf_swp  *** 
     
    588588      !!---------------------------------------------------------------------- 
    589589      INTEGER, INTENT( in ) ::   kt   ! time step 
    590       INTEGER, INTENT( in ) ::   kNm1, kNnn, kNm1_2lev, kNnn_2lev ! time level indices 
    591590      ! 
    592591      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    606605      ! Time filter and swap of scale factors 
    607606      ! ===================================== 
    608       ! - ML - e3(t/u/v)_b are already computed in dynnxt. 
     607      ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 
    609608      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    610609         IF( neuler == 0 .AND. kt == nit000 ) THEN 
     
    616615         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    617616      ENDIF 
     617      gdept_b(:,:,:) = gdept_n(:,:,:) 
     618      gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     619 
     620      e3t_n(:,:,:) = e3t_a(:,:,:) 
     621      e3u_n(:,:,:) = e3u_a(:,:,:) 
     622      e3v_n(:,:,:) = e3v_a(:,:,:) 
    618623 
    619624      ! Compute all missing vertical scale factor and depths 
     
    624629      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    625630       
    626       CALL dom_vvl_interpol( e3u(:,:,:,kNnn), e3f(:,:,:), 'F'  ) 
     631      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
    627632       
    628633      ! Vertical scale factor interpolations 
    629       CALL dom_vvl_interpol( e3t(:,:,:,kNnn),  e3w(:,:,:,kNnn_2lev), 'W'  ) 
    630       CALL dom_vvl_interpol( e3u(:,:,:,kNnn), e3uw(:,:,:,kNnn_2lev), 'UW' ) 
    631       CALL dom_vvl_interpol( e3v(:,:,:,kNnn), e3vw(:,:,:,kNnn_2lev), 'VW' ) 
    632       CALL dom_vvl_interpol( e3t(:,:,:,kNm1),  e3w(:,:,:,kNm1_2lev), 'W'  ) 
    633       CALL dom_vvl_interpol( e3u(:,:,:,kNm1), e3uw(:,:,:,kNm1_2lev), 'UW' ) 
    634       CALL dom_vvl_interpol( e3v(:,:,:,kNm1), e3vw(:,:,:,kNm1_2lev), 'VW' ) 
     634      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  ) 
     635      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
     636      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     637      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  ) 
     638      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
     639      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    635640 
    636641      ! t- and w- points depth (set the isf depth as it is in the initial step) 
    637       gdept(:,:,1,kNnn_2lev) = 0.5_wp * e3w(:,:,1,kNnn_2lev) 
    638       gdepw(:,:,1,kNnn_2lev) = 0.0_wp 
    639       gde3w(:,:,1)           = gdept(:,:,1,kNnn_2lev) - sshn(:,:) 
     642      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
     643      gdepw_n(:,:,1) = 0.0_wp 
     644      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    640645      DO jk = 2, jpk 
    641646         DO jj = 1,jpj 
     
    644649                                                                 ! 1 for jk = mikt 
    645650               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    646                gdepw(ji,jj,jk,kNnn_2lev) = gdepw(ji,jj,jk-1,kNnn_2lev) + e3t(ji,jj,jk-1,kNnn) 
    647                gdept(ji,jj,jk,kNnn_2lev) =    zcoef  * ( gdepw(ji,jj,jk,kNnn_2lev  ) + 0.5 * e3w(ji,jj,jk,kNnn_2lev) )  & 
    648                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,kNnn_2lev) +       e3w(ji,jj,jk,kNnn_2lev) )  
    649                gde3w(ji,jj,jk) = gdept(ji,jj,jk,kNnn_2lev) - sshn(ji,jj) 
     651               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
     652               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
     653                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
     654               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    650655            END DO 
    651656         END DO 
     
    654659      ! Local depth and Inverse of the local depth of the water 
    655660      ! ------------------------------------------------------- 
    656       ! 
    657       ht(:,:) = e3t(:,:,1,kNnn) * tmask(:,:,1) 
     661      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:) 
     662      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
     663      ! 
     664      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
    658665      DO jk = 2, jpkm1 
    659          ht(:,:) = ht(:,:) + e3t(:,:,jk,kNnn) * tmask(:,:,jk) 
     666         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    660667      END DO 
    661668 
  • 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 ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10789 r10799  
    277277!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine 
    278278                         CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap 
    279                          CALL dyn_nxt       ( kstp, Nm1, Np1, uu(:,:,:,Nnn), vv(:,:,:,Nnn), e3t(:,:,:,Nnn), e3u(:,:,:,Nnn), e3v(:,:,:,Nnn) )  !  
     279                         CALL dyn_nxt       ( kstp )  !  
    280280                         CALL ssh_swp       ( kstp )  ! swap of sea surface height 
    281       ! 
    282       ! Swap time levels 
    283       IF( .NOT. (neuler == 0 .AND. kstp == nit000)  ) THEN 
    284          Nrhs = Nm1 
    285          Nm1 = Nnn 
    286          Nnn = Np1 
    287          Np1 = Nrhs 
    288       ENDIF 
    289       ! 
    290       ! Update temporary pointers 
    291       CALL update_pointers() 
    292       ! 
    293       ! Note that 2-time-level indices don't need to be swapped because both "before" and "now" fields are derived in dom_vvl_sf_swp 
    294       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nm1, Nnn, Nm1_2lev, Nnn_2lev )  ! interpolate vertical scale factors for Nnn time level 
     281      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! interpolate vertical scale factors for Nnn time level 
    295282      ! 
    296283      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics 
Note: See TracChangeset for help on using the changeset viewer.