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 6060 for branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90 – NEMO

Ignore:
Timestamp:
2015-12-16T10:25:22+01:00 (8 years ago)
Author:
timgraham
Message:

Merged dev_r5836_noc2_VVL_BY_DEFAULT into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

    r5930 r6060  
    4141    
    4242   !! * Substitutions 
    43 #  include "domzgr_substitute.h90" 
    4443#  include "vectopt_loop_substitute.h90" 
    4544   !!---------------------------------------------------------------------- 
     
    6261      !!              ld_msc_ups=T :  
    6362      !! 
    64       !! ** Action  : - update (ta,sa) with the now advective tracer trends 
    65       !!              - send trends to trdtra module for further diagnostcs 
     63      !! ** Action : - update pta  with the now advective tracer trends 
     64      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     65      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    6666      !! 
    6767      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 
     
    116116      ENDIF  
    117117      !       
    118       !                                                     ! =========== 
    119       DO jn = 1, kjpt                                       ! tracer loop 
    120          !                                                  ! =========== 
    121          ! I. Horizontal advective fluxes 
    122          ! ------------------------------ 
    123          ! first guess of the slopes 
    124          zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values 
    125          ! interior values 
    126          DO jk = 1, jpkm1 
     118      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
     119         ! 
     120         !                          !* Horizontal advective fluxes 
     121         ! 
     122         !                                !-- first guess of the slopes 
     123         zwx(:,:,jpk) = 0._wp                   ! bottom values 
     124         zwy(:,:,jpk) = 0._wp   
     125         DO jk = 1, jpkm1                       ! interior values 
    127126            DO jj = 1, jpjm1       
    128127               DO ji = 1, fs_jpim1   ! vector opt. 
     
    132131           END DO 
    133132         END DO 
    134          ! 
    135          CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign) 
     133         CALL lbc_lnk( zwx, 'U', -1. )          ! lateral boundary conditions   (changed sign) 
    136134         CALL lbc_lnk( zwy, 'V', -1. ) 
    137          !                                             !-- Slopes of tracer 
    138          zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values 
    139          DO jk = 1, jpkm1                                     ! interior values 
     135         !                                !-- Slopes of tracer 
     136         zslpx(:,:,jpk) = 0._wp                 ! bottom values 
     137         zslpy(:,:,jpk) = 0._wp 
     138         DO jk = 1, jpkm1                       ! interior values 
    140139            DO jj = 2, jpj 
    141140               DO ji = fs_2, jpi   ! vector opt. 
     
    148147         END DO 
    149148         ! 
    150          DO jk = 1, jpkm1                                     ! Slopes limitation 
     149         DO jk = 1, jpkm1                 !-- Slopes limitation 
    151150            DO jj = 2, jpj 
    152151               DO ji = fs_2, jpi   ! vector opt. 
     
    161160         END DO 
    162161         ! 
    163          !                                             !-- MUSCL horizontal advective fluxes 
    164          DO jk = 1, jpkm1                                     ! interior values 
     162         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes 
    165163            zdt  = p2dt(jk) 
    166164            DO jj = 2, jpjm1 
     
    169167                  z0u = SIGN( 0.5, pun(ji,jj,jk) ) 
    170168                  zalpha = 0.5 - z0u 
    171                   zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     169                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    172170                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
    173171                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
     
    176174                  z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 
    177175                  zalpha = 0.5 - z0v 
    178                   zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 
     176                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt * r1_e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
    179177                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
    180178                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
     
    185183         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
    186184         ! 
    187          DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend 
     185         DO jk = 1, jpkm1                 !-- Tracer advective trend 
    188186            DO jj = 2, jpjm1       
    189187               DO ji = fs_2, fs_jpim1   ! vector opt. 
    190188                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
    191189                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
    192                   &                                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     190                  &                                   * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    193191               END DO 
    194192           END DO 
    195193         END DO         
    196          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     194         !                                ! trend diagnostics 
    197195         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
    198196            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN 
     
    200198            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    201199         END IF 
    202          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     200         !                                ! "Poleward" heat and salt transports 
    203201         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    204202            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    205203            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    206204         ENDIF 
    207  
    208          ! II. Vertical advective fluxes 
    209          ! ----------------------------- 
     205         ! 
     206         !                          !* Vertical advective fluxes 
     207         ! 
    210208         !                                !-- first guess of the slopes 
    211209         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions 
    212          zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions 
    213          DO jk = 2, jpkm1                                     ! interior values 
     210         zwx(:,:,jpk) = 0._wp 
     211         DO jk = 2, jpkm1                       ! interior values 
    214212            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 
    215213         END DO 
    216  
    217214         !                                !-- Slopes of tracer 
    218215         zslpx(:,:,1) = 0._wp                   ! surface values 
     
    220217            DO jj = 1, jpj 
    221218               DO ji = 1, jpi 
    222                   zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   & 
    223                      &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 
    224                END DO 
    225             END DO 
    226          END DO 
    227          !                                !-- Slopes limitation 
    228          DO jk = 2, jpkm1                       ! interior values 
    229             DO jj = 1, jpj 
     219                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
     220                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
     221               END DO 
     222            END DO 
     223         END DO 
     224         DO jk = 2, jpkm1                 !-- Slopes limitation 
     225            DO jj = 1, jpj                      ! interior values 
    230226               DO ji = 1, jpi 
    231227                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     
    235231            END DO 
    236232         END DO 
    237          !                                !-- vertical advective flux 
    238          DO jk = 1, jpkm1                       ! interior values 
     233         DO jk = 1, jpk-2                 !-- vertical advective flux 
    239234            zdt  = p2dt(jk) 
    240235            DO jj = 2, jpjm1       
     
    242237                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    243238                  zalpha = 0.5 + z0w 
    244                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
     239                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
    245240                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    246241                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
     
    249244            END DO 
    250245         END DO 
    251          !                                      ! top values  (bottom already set to zero) 
    252          IF( lk_vvl ) THEN                            !  variable volume 
    253             zwx(:,:, 1 ) = 0._wp                            ! k=1 only as zwx has been multiplied by wmask 
    254          ELSE                                         ! linear free surface  
    255             IF( ln_isfcav ) THEN                            ! ice-shelf cavities (top of the ocean) 
     246         IF( ln_linssh ) THEN                   ! top values, linear free surface only 
     247            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean) 
    256248               DO jj = 1, jpj 
    257249                  DO ji = 1, jpi 
     
    259251                  END DO 
    260252               END DO    
    261             ELSE                                             ! no cavities: only at the ocean surface 
     253            ELSE                                      ! no cavities: only at the ocean surface 
    262254               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    263255            ENDIF 
    264256         ENDIF 
    265257         ! 
    266          DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend 
     258         DO jk = 1, jpkm1                 !-- vertical advective trend 
    267259            DO jj = 2, jpjm1       
    268260               DO ji = fs_2, fs_jpim1   ! vector opt. 
    269                   pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    270                END DO 
    271             END DO 
    272          END DO 
    273          !                                 ! Save the vertical advective trends for diagnostic 
     261                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     262               END DO 
     263            END DO 
     264         END DO 
     265         !                                ! send trends for diagnostic 
    274266         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
    275267            &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
    276268            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    277269         ! 
    278       END DO 
     270      END DO                     ! end of tracer loop 
    279271      ! 
    280272      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy ) 
Note: See TracChangeset for help on using the changeset viewer.