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 7351 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 – NEMO

Ignore:
Timestamp:
2016-11-28T17:04:10+01:00 (7 years ago)
Author:
emanuelaclementi
Message:

ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r5930 r7351  
    3939 
    4040   !! * Substitutions 
    41 #  include "domzgr_substitute.h90" 
    4241#  include "vectopt_loop_substitute.h90" 
    4342   !!---------------------------------------------------------------------- 
     
    5352      !!                  ***  ROUTINE tra_adv_fct  *** 
    5453      !!  
    55       !! **  Purpose :   Compute the now trend due to total advection of  
    56       !!       tracers and add it to the general trend of tracer equations 
     54      !! **  Purpose :   Compute the now trend due to total advection of tracers 
     55      !!               and add it to the general trend of tracer equations 
    5756      !! 
    5857      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction 
    5958      !!               (choice through the value of kn_fct) 
    60       !!               - 4th order compact scheme on the vertical  
     59      !!               - on the vertical the 4th order is a compact scheme  
    6160      !!               - corrected flux (monotonic correction)  
    6261      !! 
    63       !! ** Action : - update (pta) with the now advective tracer trends 
    64       !!             - send the trends for further diagnostics 
     62      !! ** Action : - update pta  with the now advective tracer trends 
     63      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     64      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    6565      !!---------------------------------------------------------------------- 
    6666      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
     
    7070      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4) 
    7171      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    72       REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     72      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    7373      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    7474      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
     
    7676      ! 
    7777      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices   
    78       REAL(wp) ::   z2dtt, ztra                              ! local scalar 
     78      REAL(wp) ::   ztra                                     ! local scalar 
    7979      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      - 
    8080      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      - 
     
    101101      ENDIF 
    102102      ! 
    103       !                                         ! surface & bottom value : flux set to zero one for all 
    104       IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp                ! except at the surface in linear free surface case 
     103      !                          ! surface & bottom value : flux set to zero one for all 
     104      zwz(:,:, 1 ) = 0._wp             
    105105      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp 
    106106      ! 
    107107      zwi(:,:,:) = 0._wp         
    108       !                                                          ! =========== 
    109       DO jn = 1, kjpt                                            ! tracer loop 
    110          !                                                       ! =========== 
     108      ! 
     109      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
    111110         ! 
    112111         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
     
    126125         END DO 
    127126         !                    !* upstream tracer flux in the k direction *! 
    128          DO jk = 2, jpkm1         ! Interior value ( multiplied by wmask) 
     127         DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    129128            DO jj = 1, jpj 
    130129               DO ji = 1, jpi 
     
    135134            END DO 
    136135         END DO 
    137          !                     
    138          IF(.NOT.lk_vvl ) THEN   ! top ocean value (only in linear free surface as zwz has been w-masked) 
     136         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked) 
    139137            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
    140138               DO jj = 1, jpj 
     
    149147         !                
    150148         DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
    151             z2dtt = p2dt(jk) 
    152149            DO jj = 2, jpjm1 
    153150               DO ji = fs_2, fs_jpim1   ! vector opt. 
    154                   ! total intermediate advective trends 
     151                  !                             ! total intermediate advective trends 
    155152                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    156153                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    157                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    158                   ! update and guess with monotonic sheme 
    159 !!gm why tmask added in the two following lines ???    the mask is done in tranxt ! 
    160                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk) 
    161                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     154                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     155                  !                             ! update and guess with monotonic sheme 
     156                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     157                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    162158               END DO 
    163159            END DO 
     
    166162         !                 
    167163         IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    168             ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
     164            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
    169165         END IF 
    170166         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     
    174170         ENDIF 
    175171         ! 
    176          ! 
    177172         !        !==  anti-diffusive flux : high order minus low order  ==! 
    178173         ! 
    179          SELECT CASE( kn_fct_h )         !* horizontal anti-diffusive fluxes 
    180          ! 
    181          CASE(  2  )                         ! 2nd order centered 
     174         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes 
     175         ! 
     176         CASE(  2  )                   !- 2nd order centered 
    182177            DO jk = 1, jpkm1 
    183178               DO jj = 1, jpjm1 
     
    189184            END DO 
    190185            ! 
    191          CASE(  4  )                         ! 4th order centered 
    192             zltu(:,:,jpk) = 0._wp                            ! Bottom value : flux set to zero 
     186         CASE(  4  )                   !- 4th order centered 
     187            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero 
    193188            zltv(:,:,jpk) = 0._wp 
    194             DO jk = 1, jpkm1                                 ! Laplacian 
    195                DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     189            DO jk = 1, jpkm1                 ! Laplacian 
     190               DO jj = 1, jpjm1                    ! 1st derivative (gradient) 
    196191                  DO ji = 1, fs_jpim1   ! vector opt. 
    197192                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     
    199194                  END DO 
    200195               END DO 
    201                DO jj = 2, jpjm1                                   !  
     196               DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6 
    202197                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    203198                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
     
    206201               END DO 
    207202            END DO 
    208             ! 
    209203            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    210204            ! 
    211             DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     205            DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    212206               DO jj = 1, jpjm1 
    213207                  DO ji = 1, fs_jpim1   ! vector opt. 
     
    221215            END DO          
    222216            ! 
    223          CASE(  41 )                         ! 4th order centered       ==>>   !!gm coding attempt   need to be tested 
    224             ztu(:,:,jpk) = 0._wp                             ! Bottom value : flux set to zero 
     217         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested 
     218            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
    225219            ztv(:,:,jpk) = 0._wp 
    226             DO jk = 1, jpkm1                                 ! gradient 
    227                DO jj = 1, jpjm1                                   ! First derivative (gradient) 
     220            DO jk = 1, jpkm1                 ! 1st derivative (gradient) 
     221               DO jj = 1, jpjm1 
    228222                  DO ji = 1, fs_jpim1   ! vector opt. 
    229223                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     
    234228            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn) 
    235229            ! 
    236             DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     230            DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    237231               DO jj = 2, jpjm1 
    238232                  DO ji = 2, fs_jpim1   ! vector opt. 
     
    250244            ! 
    251245         END SELECT 
    252          !                                !* vertical anti-diffusive fluxes 
    253          SELECT CASE( kn_fct_v )                ! Interior values (w-masked) 
    254          ! 
    255          CASE(  2  )                                  ! 2nd order centered 
     246         !                       
     247         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values) 
     248         ! 
     249         CASE(  2  )                   !- 2nd order centered 
    256250            DO jk = 2, jpkm1     
    257251               DO jj = 2, jpjm1 
    258252                  DO ji = fs_2, fs_jpim1 
    259                      zwz(ji,jj,jk) =  ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
    260                                        - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    261                   END DO 
    262                END DO 
    263             END DO 
    264             ! 
    265          CASE(  4  )                                  ! 4th order COMPACT 
    266             !     
    267             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! COMPACT interpolation of T at w-point 
    268             ! 
     253                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
     254                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
     255                  END DO 
     256               END DO 
     257            END DO 
     258            ! 
     259         CASE(  4  )                   !- 4th order COMPACT 
     260            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
    269261            DO jk = 2, jpkm1 
    270262               DO jj = 2, jpjm1 
     
    276268            ! 
    277269         END SELECT 
    278          !                                      ! top ocean value: high order = upstream  ==>>  zwz=0 
    279          zwz(:,:, 1 ) = 0._wp                   ! only ocean surface as interior zwz values have been w-masked 
     270         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0 
     271            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
     272         ENDIF 
    280273         ! 
    281274         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    282275         CALL lbc_lnk( zwz, 'W',  1. ) 
    283  
     276         ! 
    284277         !        !==  monotonicity algorithm  ==! 
    285278         ! 
    286279         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
    287  
    288  
     280         ! 
    289281         !        !==  final trend with corrected fluxes  ==! 
    290282         ! 
     
    295287                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    296288                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) & 
    297                      &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    298                END DO 
    299             END DO 
    300          END DO 
    301          ! 
    302          IF( l_trd ) THEN                 ! trend diagnostics (contribution of upstream fluxes) 
     289                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     290               END DO 
     291            END DO 
     292         END DO 
     293         ! 
     294         IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    303295            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    304296            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     
    311303            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    312304         END IF 
    313          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     305         !                    ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    314306         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    315            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 
    316            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 
     307           IF( jn == jp_tem )   htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 
     308           IF( jn == jp_sal )   str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 
    317309         ENDIF 
    318310         ! 
    319       END DO 
     311      END DO                     ! end of tracer loop 
    320312      ! 
    321313      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     
    348340      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    349341      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps 
    350       REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     342      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    351343      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    352344      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
     
    354346      ! 
    355347      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection 
    356       REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep 
     348      REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep 
    357349      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices   
    358350      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps 
    359351      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps 
    360352      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection 
    361       REAL(wp) ::   z2dtt, ztra              ! local scalar 
     353      REAL(wp) ::   ztra            ! local scalar 
    362354      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
    363355      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      - 
     
    371363      ! 
    372364      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
    373       CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
     365      CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
    374366      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
    375367      ! 
     
    390382      zwi(:,:,:) = 0._wp 
    391383      z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 
    392       zr_p2dt(:) = 1._wp / p2dt(:) 
     384      zr_p2dt = 1._wp / p2dt 
     385      ! 
     386      ! surface & Bottom value : flux set to zero for all tracers 
     387      zwz(:,:, 1 ) = 0._wp 
     388      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
     389      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
    393390      ! 
    394391      !                                                          ! =========== 
    395392      DO jn = 1, kjpt                                            ! tracer loop 
    396393         !                                                       ! =========== 
    397          ! 1. Bottom value : flux set to zero 
    398          ! ---------------------------------- 
    399          zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
    400          zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
    401  
    402          ! 2. upstream advection with initial mass fluxes & intermediate update 
    403          ! -------------------------------------------------------------------- 
    404          ! upstream tracer flux in the i and j direction 
    405          DO jk = 1, jpkm1 
     394         ! 
     395         ! Upstream advection with initial mass fluxes & intermediate update 
     396         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction 
    406397            DO jj = 1, jpjm1 
    407398               DO ji = 1, fs_jpim1   ! vector opt. 
     
    416407            END DO 
    417408         END DO 
    418  
    419          ! upstream tracer flux in the k direction 
    420          DO jk = 2, jpkm1         ! Interior value 
     409         !                       ! upstream tracer flux in the k direction 
     410         DO jk = 2, jpkm1              ! Interior value 
    421411            DO jj = 1, jpj 
    422412               DO ji = 1, jpi 
     
    427417            END DO 
    428418         END DO 
    429          !                       ! top value 
    430          IF( lk_vvl ) THEN             ! variable volume: only k=1 as zwz is multiplied by wmask 
    431             zwz(:,:, 1 ) = 0._wp 
    432          ELSE                          ! linear free surface 
    433             IF( ln_isfcav ) THEN             ! ice-shelf cavities 
     419         IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask) 
     420            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value 
    434421               DO jj = 1, jpj 
    435422                  DO ji = 1, jpi 
     
    437424                  END DO 
    438425               END DO    
    439             ELSE                             ! standard case 
     426            ELSE                             ! no cavities, surface value 
    440427               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    441428            ENDIF 
     
    443430         ! 
    444431         DO jk = 1, jpkm1         ! total advective trend 
    445             z2dtt = p2dt(jk) 
    446432            DO jj = 2, jpjm1 
    447433               DO ji = fs_2, fs_jpim1   ! vector opt. 
    448                   ! total intermediate advective trends 
     434                  !                             ! total intermediate advective trends 
    449435                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    450436                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    451                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    452                   ! update and guess with monotonic sheme 
    453                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
    454                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     437                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) 
     438                  !                             ! update and guess with monotonic sheme 
     439                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     440                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    455441               END DO 
    456442            END DO 
     
    497483            END DO 
    498484         END DO 
    499        
     485         ! 
    500486         !                                !* vertical anti-diffusive flux 
    501487         zwz_sav(:,:,:)   = zwz(:,:,:) 
    502488         ztrs   (:,:,:,1) = ptb(:,:,:,jn) 
     489         ztrs   (:,:,1,2) = ptb(:,:,1,jn) 
     490         ztrs   (:,:,1,3) = ptb(:,:,1,jn) 
    503491         zwzts  (:,:,:)   = 0._wp 
    504          IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp    ! surface value set to zero in vvl case 
    505492         ! 
    506493         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop 
     
    508495            IF( jl == 1 ) THEN                        ! Euler forward to kick things off 
    509496               jtb = 1   ;   jtn = 1   ;   jta = 2 
    510                zts(:) = p2dt(:) * z_rzts 
     497               zts(:) = p2dt * z_rzts 
    511498               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux 
    512499               !                                            ! starting at jl =1 if kn_fct_zts is odd;  
     
    514501            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step 
    515502               jtb = 1   ;   jtn = 2   ;   jta = 3 
    516                zts(:) = 2._wp * p2dt(:) * z_rzts 
     503               zts(:) = 2._wp * p2dt * z_rzts 
    517504            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps 
    518505               jtb = MOD(jtb,3) + 1 
     
    528515               END DO 
    529516            END DO 
    530             IF(.NOT.lk_vvl ) THEN                    ! top value (only in linear free surface case) 
     517            IF( ln_linssh ) THEN                    ! top value (only in linear free surface case) 
    531518               IF( ln_isfcav ) THEN                      ! ice-shelf cavities 
    532519                  DO jj = 1, jpj 
     
    535522                     END DO 
    536523                  END DO    
    537                ELSE                                      ! standard case 
     524               ELSE                                      ! no ocean cavities 
    538525                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    539526               ENDIF 
     
    547534                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 & 
    548535                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
    549                         &                         / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     536                        &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    550537                  END DO 
    551538               END DO 
     
    557544            DO jj = 2, jpjm1 
    558545               DO ji = fs_2, fs_jpim1 
    559                   zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 
     546                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 
    560547               END DO 
    561548            END DO 
     
    576563                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       & 
    577564                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   & 
    578                      &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     565                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    579566               END DO 
    580567            END DO 
     
    623610      !!       in-space based differencing for fluid 
    624611      !!---------------------------------------------------------------------- 
    625       REAL(wp), DIMENSION(jpk)         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     612      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step 
    626613      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
    627614      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions 
     
    629616      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    630617      INTEGER  ::   ikm1         ! local integer 
    631       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars 
     618      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars 
    632619      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    633620      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo 
     
    652639      DO jk = 1, jpkm1 
    653640         ikm1 = MAX(jk-1,1) 
    654          z2dtt = p2dt(jk) 
    655641         DO jj = 2, jpjm1 
    656642            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    679665 
    680666               ! up & down beta terms 
    681                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
     667               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 
    682668               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
    683669               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
Note: See TracChangeset for help on using the changeset viewer.