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 12377 for NEMO/trunk/src/OCE/TRA/traadv_mus.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/TRA/traadv_mus.F90

    r11993 r12377  
    4646 
    4747   !! * Substitutions 
    48 #  include "vectopt_loop_substitute.h90" 
     48#  include "do_loop_substitute.h90" 
    4949   !!---------------------------------------------------------------------- 
    5050   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5454CONTAINS 
    5555 
    56    SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             & 
    57       &                                              ptb, pta, kjpt, ld_msc_ups ) 
     56   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW,             & 
     57      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups ) 
    5858      !!---------------------------------------------------------------------- 
    5959      !!                    ***  ROUTINE tra_adv_mus  *** 
    … …  
    6666      !!              ld_msc_ups=T :  
    6767      !! 
    68       !! ** Action : - update pta  with the now advective tracer trends 
     68      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    6969      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
    70       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
     70      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
    7171      !! 
    7272      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 
    7373      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 
    7474      !!---------------------------------------------------------------------- 
    75       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    76       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    77       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    78       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    79       LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
    80       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    81       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    82       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field 
    83       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
     75      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     76      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     77      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     78      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     79      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     80      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
     81      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    8484      ! 
    8585      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    120120      l_ptr = .FALSE. 
    121121      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    122       IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     122      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.  
    123123      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    124124         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
    … …  
    131131         zwx(:,:,jpk) = 0._wp                   ! bottom values 
    132132         zwy(:,:,jpk) = 0._wp   
    133          DO jk = 1, jpkm1                       ! interior values 
    134             DO jj = 1, jpjm1       
    135                DO ji = 1, fs_jpim1   ! vector opt. 
    136                   zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) ) 
    137                   zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
    138                END DO 
    139            END DO 
    140          END DO 
     133         DO_3D_10_10( 1, jpkm1 ) 
     134            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     135            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     136         END_3D 
    141137         ! lateral boundary conditions   (changed sign) 
    142138         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. ) 
    … …  
    144140         zslpx(:,:,jpk) = 0._wp                 ! bottom values 
    145141         zslpy(:,:,jpk) = 0._wp 
    146          DO jk = 1, jpkm1                       ! interior values 
    147             DO jj = 2, jpj 
    148                DO ji = fs_2, jpi   ! vector opt. 
    149                   zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
    150                      &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
    151                   zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   & 
    152                      &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
    153                END DO 
    154             END DO 
    155          END DO 
    156          ! 
    157          DO jk = 1, jpkm1                 !-- Slopes limitation 
    158             DO jj = 2, jpj 
    159                DO ji = fs_2, jpi   ! vector opt. 
    160                   zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
    161                      &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
    162                      &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
    163                   zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
    164                      &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
    165                      &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
    166                END DO 
    167            END DO 
    168          END DO 
    169          ! 
    170          DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes 
    171             DO jj = 2, jpjm1 
    172                DO ji = fs_2, fs_jpim1   ! vector opt. 
    173                   ! MUSCL fluxes 
    174                   z0u = SIGN( 0.5, pun(ji,jj,jk) ) 
    175                   zalpha = 0.5 - z0u 
    176                   zu  = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    177                   zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
    178                   zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
    179                   zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    180                   ! 
    181                   z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 
    182                   zalpha = 0.5 - z0v 
    183                   zv  = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    184                   zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
    185                   zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
    186                   zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    187                END DO 
    188             END DO 
    189          END DO 
     142         DO_3D_01_01( 1, jpkm1 ) 
     143            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
     144               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
     145            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   & 
     146               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
     147         END_3D 
     148         ! 
     149         DO_3D_01_01( 1, jpkm1 ) 
     150            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
     151               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
     152               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
     153            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
     154               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
     155               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
     156         END_3D 
     157         ! 
     158         DO_3D_00_00( 1, jpkm1 ) 
     159            ! MUSCL fluxes 
     160            z0u = SIGN( 0.5, pU(ji,jj,jk) ) 
     161            zalpha = 0.5 - z0u 
     162            zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     163            zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
     164            zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
     165            zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     166            ! 
     167            z0v = SIGN( 0.5, pV(ji,jj,jk) ) 
     168            zalpha = 0.5 - z0v 
     169            zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     170            zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
     171            zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
     172            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     173         END_3D 
    190174         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
    191175         ! 
    192          DO jk = 1, jpkm1                 !-- Tracer advective trend 
    193             DO jj = 2, jpjm1       
    194                DO ji = fs_2, fs_jpim1   ! vector opt. 
    195                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
    196                   &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
    197                   &                                   * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    198                END DO 
    199            END DO 
    200          END DO         
     176         DO_3D_00_00( 1, jpkm1 ) 
     177            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
     178            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
     179            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     180         END_3D 
    201181         !                                ! trend diagnostics 
    202182         IF( l_trd )  THEN 
    203             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
    204             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
     183            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) ) 
     184            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) ) 
    205185         END IF 
    206186         !                                 ! "Poleward" heat and salt transports  
    … …  
    215195         zwx(:,:,jpk) = 0._wp 
    216196         DO jk = 2, jpkm1                       ! interior values 
    217             zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 
     197            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) ) 
    218198         END DO 
    219199         !                                !-- Slopes of tracer 
    220200         zslpx(:,:,1) = 0._wp                   ! surface values 
    221          DO jk = 2, jpkm1                       ! interior value 
    222             DO jj = 1, jpj 
    223                DO ji = 1, jpi 
    224                   zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
    225                      &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
    226                END DO 
    227             END DO 
    228          END DO 
    229          DO jk = 2, jpkm1                 !-- Slopes limitation 
    230             DO jj = 1, jpj                      ! interior values 
    231                DO ji = 1, jpi 
    232                   zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
    233                      &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
    234                      &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  ) 
    235                END DO 
    236             END DO 
    237          END DO 
    238          DO jk = 1, jpk-2                 !-- vertical advective flux 
    239             DO jj = 2, jpjm1       
    240                DO ji = fs_2, fs_jpim1   ! vector opt. 
    241                   z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    242                   zalpha = 0.5 + z0w 
    243                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
    244                   zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    245                   zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    246                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
    247                END DO  
    248             END DO 
    249          END DO 
     201         DO_3D_11_11( 2, jpkm1 ) 
     202            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
     203               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
     204         END_3D 
     205         DO_3D_11_11( 2, jpkm1 ) 
     206            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     207               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
     208               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  ) 
     209         END_3D 
     210         DO_3D_00_00( 1, jpk-2 ) 
     211            z0w = SIGN( 0.5, pW(ji,jj,jk+1) ) 
     212            zalpha = 0.5 + z0w 
     213            zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm) 
     214            zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
     215            zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
     216            zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
     217         END_3D 
    250218         IF( ln_linssh ) THEN                   ! top values, linear free surface only 
    251219            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean) 
    252                DO jj = 1, jpj 
    253                   DO ji = 1, jpi 
    254                      zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
    255                   END DO 
    256                END DO    
     220               DO_2D_11_11 
     221                  zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) 
     222               END_2D 
    257223            ELSE                                      ! no cavities: only at the ocean surface 
    258                zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     224               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 
    259225            ENDIF 
    260226         ENDIF 
    261227         ! 
    262          DO jk = 1, jpkm1                 !-- vertical advective trend 
    263             DO jj = 2, jpjm1       
    264                DO ji = fs_2, fs_jpim1   ! vector opt. 
    265                   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) 
    266                END DO 
    267             END DO 
    268          END DO 
     228         DO_3D_00_00( 1, jpkm1 ) 
     229            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     230         END_3D 
    269231         !                                ! send trends for diagnostic 
    270          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     232         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) ) 
    271233         ! 
    272234      END DO                     ! end of tracer loop 
Note: See TracChangeset for help on using the changeset viewer.