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_mus.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_mus.F90

    r5930 r7351  
    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 
     
    7373      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    7474      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
    75       REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     75      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    7676      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    7777      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field 
     
    8282      REAL(wp) ::   zu, z0u, zzwx, zw    ! local scalars 
    8383      REAL(wp) ::   zv, z0v, zzwy, z0w   !   -      - 
    84       REAL(wp) ::   zdt, zalpha          !   -      - 
     84      REAL(wp) ::   zalpha               !   -      - 
    8585      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace 
    8686      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -  
     
    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 
    165             zdt  = p2dt(jk) 
     162         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes 
    166163            DO jj = 2, jpjm1 
    167164               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    169166                  z0u = SIGN( 0.5, pun(ji,jj,jk) ) 
    170167                  zalpha = 0.5 - z0u 
    171                   zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     168                  zu  = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    172169                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
    173170                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
     
    176173                  z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 
    177174                  zalpha = 0.5 - z0v 
    178                   zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 
     175                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    179176                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
    180177                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
     
    185182         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
    186183         ! 
    187          DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend 
     184         DO jk = 1, jpkm1                 !-- Tracer advective trend 
    188185            DO jj = 2, jpjm1       
    189186               DO ji = fs_2, fs_jpim1   ! vector opt. 
    190187                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
    191188                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
    192                   &                                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     189                  &                                   * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    193190               END DO 
    194191           END DO 
    195192         END DO         
    196          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     193         !                                ! trend diagnostics 
    197194         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
    198195            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN 
     
    200197            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    201198         END IF 
    202          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     199         !                                ! "Poleward" heat and salt transports 
    203200         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    204201            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    205202            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    206203         ENDIF 
    207  
    208          ! II. Vertical advective fluxes 
    209          ! ----------------------------- 
     204         ! 
     205         !                          !* Vertical advective fluxes 
     206         ! 
    210207         !                                !-- first guess of the slopes 
    211208         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions 
    212          zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions 
    213          DO jk = 2, jpkm1                                     ! interior values 
     209         zwx(:,:,jpk) = 0._wp 
     210         DO jk = 2, jpkm1                       ! interior values 
    214211            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 
    215212         END DO 
    216  
    217213         !                                !-- Slopes of tracer 
    218214         zslpx(:,:,1) = 0._wp                   ! surface values 
     
    220216            DO jj = 1, jpj 
    221217               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 
     218                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
     219                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
     220               END DO 
     221            END DO 
     222         END DO 
     223         DO jk = 2, jpkm1                 !-- Slopes limitation 
     224            DO jj = 1, jpj                      ! interior values 
    230225               DO ji = 1, jpi 
    231226                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     
    235230            END DO 
    236231         END DO 
    237          !                                !-- vertical advective flux 
    238          DO jk = 1, jpkm1                       ! interior values 
    239             zdt  = p2dt(jk) 
     232         DO jk = 1, jpk-2                 !-- vertical advective flux 
    240233            DO jj = 2, jpjm1       
    241234               DO ji = fs_2, fs_jpim1   ! vector opt. 
    242235                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    243236                  zalpha = 0.5 + z0w 
    244                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
     237                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
    245238                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    246239                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
     
    249242            END DO 
    250243         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) 
     244         IF( ln_linssh ) THEN                   ! top values, linear free surface only 
     245            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean) 
    256246               DO jj = 1, jpj 
    257247                  DO ji = 1, jpi 
     
    259249                  END DO 
    260250               END DO    
    261             ELSE                                             ! no cavities: only at the ocean surface 
     251            ELSE                                      ! no cavities: only at the ocean surface 
    262252               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    263253            ENDIF 
    264254         ENDIF 
    265255         ! 
    266          DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend 
     256         DO jk = 1, jpkm1                 !-- vertical advective trend 
    267257            DO jj = 2, jpjm1       
    268258               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 
     259                  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) 
     260               END DO 
     261            END DO 
     262         END DO 
     263         !                                ! send trends for diagnostic 
    274264         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
    275265            &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
    276266            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    277267         ! 
    278       END DO 
     268      END DO                     ! end of tracer loop 
    279269      ! 
    280270      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy ) 
Note: See TracChangeset for help on using the changeset viewer.