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

Ignore:
Timestamp:
2019-11-22T15:29:17+01:00 (4 years ago)
Author:
acc
Message:

Merge in changes from 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. This just creates a fresh copy of this branch to use as the merge base. See ticket #2341

Location:
NEMO/branches/2019/dev_r11943_MERGE_2019/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src

    • Property svn:mergeinfo deleted

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

    r10425 r11949  
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    50       &                                       ptb, ptn, pta, kjpt ) 
     49   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs ) 
    5150      !!---------------------------------------------------------------------- 
    5251      !!                  ***  ROUTINE tra_adv_qck  *** 
    … …  
    7271      !!         dt = 2*rdtra and the scalar values are tb and sb 
    7372      !! 
    74       !!       On the vertical, the simple centered scheme used ptn 
     73      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm) 
    7574      !! 
    7675      !!               The fluxes are bounded by the ULTIMATE limiter to 
    … …  
    7877      !!            prevent the appearance of spurious numerical oscillations 
    7978      !! 
    80       !! ** Action : - update pta  with the now advective tracer trends 
     79      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    8180      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
    8281      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    … …  
    8483      !! ** Reference : Leonard (1979, 1991) 
    8584      !!---------------------------------------------------------------------- 
    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  
     85      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     86      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     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   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    9493      !!---------------------------------------------------------------------- 
    9594      ! 
    … …  
    108107      ! 
    109108      !        ! 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 )  
     109      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )  
     110      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )  
    112111 
    113112      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    114       CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt ) 
     113      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
    115114      ! 
    116115   END SUBROUTINE tra_adv_qck 
    117116 
    118117 
    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  
     118   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
     119      !!---------------------------------------------------------------------- 
     120      !! 
     121      !!---------------------------------------------------------------------- 
     122      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     123      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     124      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     125      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     126      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     127      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components 
     128      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    131129      !! 
    132130      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    145143            DO jj = 2, jpjm1 
    146144               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 
     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 
    149147               END DO 
    150148            END DO 
    … …  
    158156            DO jj = 2, jpjm1 
    159157               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  
     158                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    161159                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T  
    162160               END DO 
    … …  
    167165            DO jj = 2, jpjm1 
    168166               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 
     167                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     168                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     169                  zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     170                  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 
     171                  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 
    174172               END DO 
    175173            END DO 
    … …  
    197195            DO jj = 2, jpjm1 
    198196               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  
     197                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    200198                  !--- If the second ustream point is a land point 
    201199                  !--- the flux is computed by the 1st order UPWIND scheme 
    202200                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 
    203201                  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) 
     202                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 
    205203               END DO 
    206204            END DO 
    … …  
    213211            DO jj = 2, jpjm1 
    214212               DO ji = fs_2, fs_jpim1   ! vector opt.   
    215                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     213                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    216214                  ! horizontal advective trends 
    217215                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 
    218216                  !--- add it to the general tracer trends 
    219                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     217                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
    220218               END DO 
    221219            END DO 
    222220         END DO 
    223221         !                                 ! trend diagnostics 
    224          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     222         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) ) 
    225223         ! 
    226224      END DO 
    … …  
    229227 
    230228 
    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  
     229   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
     230      !!---------------------------------------------------------------------- 
     231      !! 
     232      !!---------------------------------------------------------------------- 
     233      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     234      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     235      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     236      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
     237      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
     238      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components 
     239      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    243240      !! 
    244241      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices 
    … …  
    259256               DO ji = fs_2, fs_jpim1   ! vector opt. 
    260257                  ! Upstream in the x-direction for the tracer 
    261                   zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn) 
     258                  zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 
    262259                  ! Downstream in the x-direction for the tracer 
    263                   zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn) 
     260                  zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 
    264261               END DO 
    265262            END DO 
    … …  
    275272            DO jj = 2, jpjm1 
    276273               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  
     274                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    278275                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T  
    279276               END DO 
    … …  
    284281            DO jj = 2, jpjm1 
    285282               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 
     283                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
     284                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     285                  zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction) 
     286                  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 
     287                  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 
    291288               END DO 
    292289            END DO 
    … …  
    314311            DO jj = 2, jpjm1 
    315312               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  
     313                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0  
    317314                  !--- If the second ustream point is a land point 
    318315                  !--- the flux is computed by the 1st order UPWIND scheme 
    319316                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 
    320317                  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) 
     318                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 
    322319               END DO 
    323320            END DO 
    … …  
    330327            DO jj = 2, jpjm1 
    331328               DO ji = fs_2, fs_jpim1   ! vector opt.   
    332                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     329                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    333330                  ! horizontal advective trends 
    334331                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 
    335332                  !--- add it to the general tracer trends 
    336                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     333                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra 
    337334               END DO 
    338335            END DO 
    339336         END DO 
    340337         !                                 ! trend diagnostics 
    341          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     338         IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) 
    342339         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    343340         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    … …  
    348345 
    349346 
    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  
     347   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
     348      !!---------------------------------------------------------------------- 
     349      !! 
     350      !!---------------------------------------------------------------------- 
     351      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index 
     352      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices 
     353      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     354      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
     355      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity  
     356      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    361357      ! 
    362358      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    374370            DO jj = 2, jpjm1 
    375371               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) 
     372                  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) 
    377373               END DO 
    378374            END DO 
    … …  
    382378               DO jj = 1, jpj 
    383379                  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  
     380                     zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface  
    385381                  END DO 
    386382               END DO    
    387383            ELSE                                   ! no ocean cavities (only ocean surface) 
    388                zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 
     384               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 
    389385            ENDIF 
    390386         ENDIF 
    … …  
    393389            DO jj = 2, jpjm1 
    394390               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) 
     391                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
     392                     &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    397393               END DO 
    398394            END DO 
    399395         END DO 
    400396         !                                 ! Send trends for diagnostic 
    401          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     397         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) 
    402398         ! 
    403399      END DO 
Note: See TracChangeset for help on using the changeset viewer.