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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/traadv_qck.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • 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_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/traadv_qck.F90

    r10425 r13463  
    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" 
     43#  include "domzgr_substitute.h90" 
    4244   !!---------------------------------------------------------------------- 
    4345   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4749CONTAINS 
    4850 
    49    SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    50       &                                       ptb, ptn, pta, kjpt ) 
     51   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs ) 
    5152      !!---------------------------------------------------------------------- 
    5253      !!                  ***  ROUTINE tra_adv_qck  *** 
    … …  
    7273      !!         dt = 2*rdtra and the scalar values are tb and sb 
    7374      !! 
    74       !!       On the vertical, the simple centered scheme used ptn 
     75      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm) 
    7576      !! 
    7677      !!               The fluxes are bounded by the ULTIMATE limiter to 
    … …  
    7879      !!            prevent the appearance of spurious numerical oscillations 
    7980      !! 
    80       !! ** Action : - update pta  with the now advective tracer trends 
     81      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    8182      !!             - 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) 
     83      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
    8384      !! 
    8485      !! ** Reference : Leonard (1979, 1991) 
    8586      !!---------------------------------------------------------------------- 
    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  
     87      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     88      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     89      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     90      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     91      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     92      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
     94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    9495      !!---------------------------------------------------------------------- 
    9596      ! 
    … …  
    103104      l_trd = .FALSE. 
    104105      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.  
     106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
     107      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.  
    107108      ! 
    108109      ! 
    109110      !        ! 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 )  
     111      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )  
     112      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )  
    112113 
    113114      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    114       CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt ) 
     115      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
    115116      ! 
    116117   END SUBROUTINE tra_adv_qck 
    117118 
    118119 
    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  
     120   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
     121      !!---------------------------------------------------------------------- 
     122      !! 
     123      !!---------------------------------------------------------------------- 
     124      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     125      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     126      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     127      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     128      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components 
     130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    131131      !! 
    132132      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    142142         ! 
    143143!!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 
    152          CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
     144         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     145            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer 
     146            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer 
     147         END_3D 
     148         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions  
    153149          
    154150         ! 
    155151         ! Horizontal advective fluxes 
    156152         ! --------------------------- 
    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  
     153         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     154            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     155            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T  
     156         END_3D 
     157         ! 
     158         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     159            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     160            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     161            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     162            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 
     163            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 
     164         END_3D 
    177165         !--- Lateral boundary conditions  
    178          CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. ) 
     166         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp ) 
    179167 
    180168         !--- QUICKEST scheme 
    … …  
    182170         ! 
    183171         ! 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 
    191          CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions  
     172         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     173            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 
     174         END_3D 
     175         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions  
    192176 
    193177         ! 
    … …  
    195179         DO jk = 1, jpkm1   
    196180            ! 
    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 
     181            DO_2D( 0, 0, 0, 0 ) 
     182               zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     183               !--- If the second ustream point is a land point 
     184               !--- the flux is computed by the 1st order UPWIND scheme 
     185               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 
     186               zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
     187               zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 
     188            END_2D 
    207189         END DO 
    208190         ! 
    209          CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions 
     191         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 
    210192         ! 
    211193         ! 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 
     194         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     195            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     196            ! horizontal advective trends 
     197            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 
     198            !--- add it to the general tracer trends 
     199            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
     200         END_3D 
    223201         !                                 ! trend diagnostics 
    224          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     202         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) ) 
    225203         ! 
    226204      END DO 
    … …  
    229207 
    230208 
    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  
     209   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
     210      !!---------------------------------------------------------------------- 
     211      !! 
     212      !!---------------------------------------------------------------------- 
     213      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     214      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     215      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     216      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     217      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     218      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components 
     219      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    243220      !! 
    244221      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices 
    … …  
    256233            !                                              
    257234            !--- 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 
     235            DO_2D( 0, 0, 0, 0 ) 
     236               ! Upstream in the x-direction for the tracer 
     237               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 
     238               ! Downstream in the x-direction for the tracer 
     239               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 
     240            END_2D 
    266241         END DO 
    267          CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions  
     242         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions  
    268243 
    269244          
    … …  
    272247         ! --------------------------- 
    273248         ! 
    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 
     249         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     250            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     251            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T  
     252         END_3D 
     253         ! 
     254         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     255            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     256            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     257            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     258            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 
     259            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 
     260         END_3D 
    294261 
    295262         !--- Lateral boundary conditions  
    296          CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. ) 
     263         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 
    297264 
    298265         !--- QUICKEST scheme 
    … …  
    300267         ! 
    301268         ! 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 
    309          CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions  
     269         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     270            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
     271         END_3D 
     272         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions  
    310273         ! 
    311274         ! Tracer flux on the x-direction 
    312275         DO jk = 1, jpkm1   
    313276            ! 
    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 
     277            DO_2D( 0, 0, 0, 0 ) 
     278               zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     279               !--- If the second ustream point is a land point 
     280               !--- the flux is computed by the 1st order UPWIND scheme 
     281               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 
     282               zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 
     283               zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 
     284            END_2D 
    324285         END DO 
    325286         ! 
    326          CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions 
     287         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 
    327288         ! 
    328289         ! 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 
     290         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     291            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     292            ! horizontal advective trends 
     293            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 
     294            !--- add it to the general tracer trends 
     295            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
     296         END_3D 
    340297         !                                 ! trend diagnostics 
    341          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     298         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) 
    342299         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    343300         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    … …  
    348305 
    349306 
    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  
     307   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
     308      !!---------------------------------------------------------------------- 
     309      !! 
     310      !!---------------------------------------------------------------------- 
     311      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index 
     312      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices 
     313      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     314      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
     315      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity  
     316      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    361317      ! 
    362318      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    371327         !                                                       ! =========== 
    372328         ! 
    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 
     329         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     330            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) 
     331         END_3D 
    380332         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask) 
    381333            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    
     334               DO_2D( 1, 1, 1, 1 ) 
     335                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface  
     336               END_2D 
    387337            ELSE                                   ! no ocean cavities (only ocean surface) 
    388                zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
     338               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 
    389339            ENDIF 
    390340         ENDIF 
    391341         ! 
    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 
     342         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     343            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     344               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     345         END_3D 
    400346         !                                 ! Send trends for diagnostic 
    401          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     347         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) 
    402348         ! 
    403349      END DO 
    … …  
    423369      !---------------------------------------------------------------------- 
    424370      ! 
    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 
     371      DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
     372         zc     = puc(ji,jj,jk)                         ! Courant number 
     373         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 
     374         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) ) 
     375         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 
     376         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 
     377         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST  
     378         ! 
     379         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 
     380         zcoef2 = ABS( zcoef1 ) 
     381         zcoef3 = ABS( zcurv ) 
     382         IF( zcoef3 >= zcoef2 ) THEN 
     383            zfho = pfc(ji,jj,jk)  
     384         ELSE 
     385            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF 
     386            IF( zcoef1 >= 0. ) THEN 
     387               zfho = MAX( pfc(ji,jj,jk), zfho )  
     388               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )  
     389            ELSE 
     390               zfho = MIN( pfc(ji,jj,jk), zfho )  
     391               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )  
     392            ENDIF 
     393         ENDIF 
     394         puc(ji,jj,jk) = zfho 
     395      END_3D 
    454396      ! 
    455397   END SUBROUTINE quickest 
Note: See TracChangeset for help on using the changeset viewer.