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

    r11993 r12377  
    2121   USE trdtra          ! trends manager: tracers  
    2222   USE diaptr          ! poleward transport diagnostics 
     23   USE iom 
    2324   ! 
    2425   USE in_out_manager  ! I/O manager 
    … …  
    3940 
    4041   !! * Substitutions 
    41 #  include "vectopt_loop_substitute.h90" 
     42#  include "do_loop_substitute.h90" 
    4243   !!---------------------------------------------------------------------- 
    4344   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4748CONTAINS 
    4849 
    49    SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    50       &                                       ptb, ptn, pta, kjpt ) 
     50   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs ) 
    5151      !!---------------------------------------------------------------------- 
    5252      !!                  ***  ROUTINE tra_adv_qck  *** 
    … …  
    7272      !!         dt = 2*rdtra and the scalar values are tb and sb 
    7373      !! 
    74       !!       On the vertical, the simple centered scheme used ptn 
     74      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm) 
    7575      !! 
    7676      !!               The fluxes are bounded by the ULTIMATE limiter to 
    … …  
    7878      !!            prevent the appearance of spurious numerical oscillations 
    7979      !! 
    80       !! ** Action : - update pta  with the now advective tracer trends 
     80      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    8181      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
    82       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
     82      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
    8383      !! 
    8484      !! ** Reference : Leonard (1979, 1991) 
    8585      !!---------------------------------------------------------------------- 
    86       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    87       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    88       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    89       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    90       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    91       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    92       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
     86      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     87      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     88      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     89      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     90      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     91      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    9494      !!---------------------------------------------------------------------- 
    9595      ! 
    … …  
    103103      l_trd = .FALSE. 
    104104      l_ptr = .FALSE. 
    105       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    106       IF(   cdtype == 'TRA' .AND. ln_diaptr )                                              l_ptr = .TRUE.  
     105      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
     106      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.  
    107107      ! 
    108108      ! 
    109109      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
    110       CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )  
    111       CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )  
     110      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )  
     111      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )  
    112112 
    113113      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    114       CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt ) 
     114      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
    115115      ! 
    116116   END SUBROUTINE tra_adv_qck 
    117117 
    118118 
    119    SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  & 
    120       &                                        ptb, ptn, pta, kjpt   ) 
    121       !!---------------------------------------------------------------------- 
    122       !! 
    123       !!---------------------------------------------------------------------- 
    124       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    125       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    126       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    127       REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step 
    128       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components 
    129       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields 
    130       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     119   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
     120      !!---------------------------------------------------------------------- 
     121      !! 
     122      !!---------------------------------------------------------------------- 
     123      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     124      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     125      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     126      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     127      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     128      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    131130      !! 
    132131      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    142141         ! 
    143142!!gm why not using a SHIFT instruction... 
    144          DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask 
    145             DO jj = 2, jpjm1 
    146                DO ji = fs_2, fs_jpim1   ! vector opt. 
    147                   zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer 
    148                   zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer 
    149                END DO 
    150             END DO 
    151          END DO 
     143         DO_3D_00_00( 1, jpkm1 ) 
     144            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer 
     145            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer 
     146         END_3D 
    152147         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
    153148          
    … …  
    155150         ! Horizontal advective fluxes 
    156151         ! --------------------------- 
    157          DO jk = 1, jpkm1                              
    158             DO jj = 2, jpjm1 
    159                DO ji = fs_2, fs_jpim1   ! vector opt.          
    160                   zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    161                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T  
    162                END DO 
    163             END DO 
    164          END DO 
    165          ! 
    166          DO jk = 1, jpkm1   
    167             DO jj = 2, jpjm1 
    168                DO ji = fs_2, fs_jpim1   ! vector opt.    
    169                   zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    170                   zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk) 
    171                   zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
    172                   zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T 
    173                   zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T 
    174                END DO 
    175             END DO 
    176          END DO  
     152         DO_3D_00_00( 1, jpkm1 ) 
     153            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     154            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T  
     155         END_3D 
     156         ! 
     157         DO_3D_00_00( 1, jpkm1 ) 
     158            zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     159            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     160            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     161            zfc(ji,jj,jk)  = zdir * pt(ji  ,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji+1,jj,jk,jn,Kbb)  ! FC in the x-direction for T 
     162            zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T 
     163         END_3D 
    177164         !--- Lateral boundary conditions  
    178165         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. ) 
    … …  
    182169         ! 
    183170         ! Mask at the T-points in the x-direction (mask=0 or mask=1) 
    184          DO jk = 1, jpkm1   
    185             DO jj = 2, jpjm1 
    186                DO ji = fs_2, fs_jpim1   ! vector opt.                
    187                   zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 
    188                END DO 
    189             END DO 
    190          END DO 
     171         DO_3D_00_00( 1, jpkm1 ) 
     172            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 
     173         END_3D 
    191174         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions  
    192175 
    … …  
    195178         DO jk = 1, jpkm1   
    196179            ! 
    197             DO jj = 2, jpjm1 
    198                DO ji = fs_2, fs_jpim1   ! vector opt.                
    199                   zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    200                   !--- If the second ustream point is a land point 
    201                   !--- the flux is computed by the 1st order UPWIND scheme 
    202                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 
    203                   zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
    204                   zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk) 
    205                END DO 
    206             END DO 
     180            DO_2D_00_00 
     181               zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     182               !--- If the second ustream point is a land point 
     183               !--- the flux is computed by the 1st order UPWIND scheme 
     184               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 
     185               zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
     186               zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 
     187            END_2D 
    207188         END DO 
    208189         ! 
    … …  
    210191         ! 
    211192         ! Computation of the trend 
    212          DO jk = 1, jpkm1   
    213             DO jj = 2, jpjm1 
    214                DO ji = fs_2, fs_jpim1   ! vector opt.   
    215                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    216                   ! horizontal advective trends 
    217                   ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 
    218                   !--- add it to the general tracer trends 
    219                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    220                END DO 
    221             END DO 
    222          END DO 
     193         DO_3D_00_00( 1, jpkm1 ) 
     194            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     195            ! horizontal advective trends 
     196            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 
     197            !--- add it to the general tracer trends 
     198            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
     199         END_3D 
    223200         !                                 ! trend diagnostics 
    224          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     201         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) ) 
    225202         ! 
    226203      END DO 
    … …  
    229206 
    230207 
    231    SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                & 
    232       &                                        ptb, ptn, pta, kjpt ) 
    233       !!---------------------------------------------------------------------- 
    234       !! 
    235       !!---------------------------------------------------------------------- 
    236       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    237       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    238       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    239       REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step 
    240       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components 
    241       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields 
    242       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     208   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
     209      !!---------------------------------------------------------------------- 
     210      !! 
     211      !!---------------------------------------------------------------------- 
     212      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     213      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     214      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     215      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     216      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     217      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components 
     218      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    243219      !! 
    244220      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices 
    … …  
    256232            !                                              
    257233            !--- Computation of the ustream and downstream value of the tracer and the mask 
    258             DO jj = 2, jpjm1 
    259                DO ji = fs_2, fs_jpim1   ! vector opt. 
    260                   ! Upstream in the x-direction for the tracer 
    261                   zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn) 
    262                   ! Downstream in the x-direction for the tracer 
    263                   zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn) 
    264                END DO 
    265             END DO 
     234            DO_2D_00_00 
     235               ! Upstream in the x-direction for the tracer 
     236               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 
     237               ! Downstream in the x-direction for the tracer 
     238               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 
     239            END_2D 
    266240         END DO 
    267241         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
    … …  
    272246         ! --------------------------- 
    273247         ! 
    274          DO jk = 1, jpkm1                              
    275             DO jj = 2, jpjm1 
    276                DO ji = fs_2, fs_jpim1   ! vector opt.          
    277                   zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    278                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T  
    279                END DO 
    280             END DO 
    281          END DO 
    282          ! 
    283          DO jk = 1, jpkm1   
    284             DO jj = 2, jpjm1 
    285                DO ji = fs_2, fs_jpim1   ! vector opt.    
    286                   zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    287                   zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk) 
    288                   zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
    289                   zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T 
    290                   zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T 
    291                END DO 
    292             END DO 
    293          END DO 
     248         DO_3D_00_00( 1, jpkm1 ) 
     249            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     250            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T  
     251         END_3D 
     252         ! 
     253         DO_3D_00_00( 1, jpkm1 ) 
     254            zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     255            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     256            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     257            zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T 
     258            zfd(ji,jj,jk)  = zdir * pt(ji,jj+1,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj  ,jk,jn,Kbb)  ! FD in the x-direction for T 
     259         END_3D 
    294260 
    295261         !--- Lateral boundary conditions  
    … …  
    300266         ! 
    301267         ! Mask at the T-points in the x-direction (mask=0 or mask=1) 
    302          DO jk = 1, jpkm1   
    303             DO jj = 2, jpjm1 
    304                DO ji = fs_2, fs_jpim1   ! vector opt.                
    305                   zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
    306                END DO 
    307             END DO 
    308          END DO 
     268         DO_3D_00_00( 1, jpkm1 ) 
     269            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
     270         END_3D 
    309271         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions  
    310272         ! 
    … …  
    312274         DO jk = 1, jpkm1   
    313275            ! 
    314             DO jj = 2, jpjm1 
    315                DO ji = fs_2, fs_jpim1   ! vector opt.                
    316                   zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0  
    317                   !--- If the second ustream point is a land point 
    318                   !--- the flux is computed by the 1st order UPWIND scheme 
    319                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 
    320                   zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
    321                   zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk) 
    322                END DO 
    323             END DO 
     276            DO_2D_00_00 
     277               zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     278               !--- If the second ustream point is a land point 
     279               !--- the flux is computed by the 1st order UPWIND scheme 
     280               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 
     281               zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
     282               zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 
     283            END_2D 
    324284         END DO 
    325285         ! 
    … …  
    327287         ! 
    328288         ! Computation of the trend 
    329          DO jk = 1, jpkm1   
    330             DO jj = 2, jpjm1 
    331                DO ji = fs_2, fs_jpim1   ! vector opt.   
    332                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    333                   ! horizontal advective trends 
    334                   ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 
    335                   !--- add it to the general tracer trends 
    336                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    337                END DO 
    338             END DO 
    339          END DO 
     289         DO_3D_00_00( 1, jpkm1 ) 
     290            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     291            ! horizontal advective trends 
     292            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 
     293            !--- add it to the general tracer trends 
     294            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
     295         END_3D 
    340296         !                                 ! trend diagnostics 
    341          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     297         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) 
    342298         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    343299         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    … …  
    348304 
    349305 
    350    SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           & 
    351      &                                    ptn, pta, kjpt ) 
    352       !!---------------------------------------------------------------------- 
    353       !! 
    354       !!---------------------------------------------------------------------- 
    355       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    356       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    357       INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    358       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity  
    359       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields 
    360       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend  
     306   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
     307      !!---------------------------------------------------------------------- 
     308      !! 
     309      !!---------------------------------------------------------------------- 
     310      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index 
     311      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices 
     312      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     313      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
     314      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity  
     315      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    361316      ! 
    362317      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    371326         !                                                       ! =========== 
    372327         ! 
    373          DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux) 
    374             DO jj = 2, jpjm1 
    375                DO ji = fs_2, fs_jpim1   ! vector opt. 
    376                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
    377                END DO 
    378             END DO 
    379          END DO 
     328         DO_3D_00_00( 2, jpkm1 ) 
     329            zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk) 
     330         END_3D 
    380331         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask) 
    381332            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean) 
    382                DO jj = 1, jpj 
    383                   DO ji = 1, jpi 
    384                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    385                   END DO 
    386                END DO    
     333               DO_2D_11_11 
     334                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface  
     335               END_2D 
    387336            ELSE                                   ! no ocean cavities (only ocean surface) 
    388                zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
     337               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 
    389338            ENDIF 
    390339         ENDIF 
    391340         ! 
    392          DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    393             DO jj = 2, jpjm1 
    394                DO ji = fs_2, fs_jpim1   ! vector opt. 
    395                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
    396                      &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    397                END DO 
    398             END DO 
    399          END DO 
     341         DO_3D_00_00( 1, jpkm1 ) 
     342            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     343               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     344         END_3D 
    400345         !                                 ! Send trends for diagnostic 
    401          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     346         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) 
    402347         ! 
    403348      END DO 
    … …  
    423368      !---------------------------------------------------------------------- 
    424369      ! 
    425       DO jk = 1, jpkm1 
    426          DO jj = 1, jpj 
    427             DO ji = 1, jpi 
    428                zc     = puc(ji,jj,jk)                         ! Courant number 
    429                zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 
    430                zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) ) 
    431                zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 
    432                zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 
    433                zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST  
    434                ! 
    435                zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 
    436                zcoef2 = ABS( zcoef1 ) 
    437                zcoef3 = ABS( zcurv ) 
    438                IF( zcoef3 >= zcoef2 ) THEN 
    439                   zfho = pfc(ji,jj,jk)  
    440                ELSE 
    441                   zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF 
    442                   IF( zcoef1 >= 0. ) THEN 
    443                      zfho = MAX( pfc(ji,jj,jk), zfho )  
    444                      zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )  
    445                   ELSE 
    446                      zfho = MIN( pfc(ji,jj,jk), zfho )  
    447                      zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )  
    448                   ENDIF 
    449                ENDIF 
    450                puc(ji,jj,jk) = zfho 
    451             END DO 
    452          END DO 
    453       END DO 
     370      DO_3D_11_11( 1, jpkm1 ) 
     371         zc     = puc(ji,jj,jk)                         ! Courant number 
     372         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 
     373         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) ) 
     374         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 
     375         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 
     376         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST  
     377         ! 
     378         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 
     379         zcoef2 = ABS( zcoef1 ) 
     380         zcoef3 = ABS( zcurv ) 
     381         IF( zcoef3 >= zcoef2 ) THEN 
     382            zfho = pfc(ji,jj,jk)  
     383         ELSE 
     384            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF 
     385            IF( zcoef1 >= 0. ) THEN 
     386               zfho = MAX( pfc(ji,jj,jk), zfho )  
     387               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )  
     388            ELSE 
     389               zfho = MIN( pfc(ji,jj,jk), zfho )  
     390               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )  
     391            ENDIF 
     392         ENDIF 
     393         puc(ji,jj,jk) = zfho 
     394      END_3D 
    454395      ! 
    455396   END SUBROUTINE quickest 
Note: See TracChangeset for help on using the changeset viewer.