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 10825 – NEMO

Changeset 10825


Ignore:
Timestamp:
2019-04-02T16:40:54+02:00 (5 years ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : trazdf.F90 and dynzdf.F90

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90

    r10364 r10825  
    4545CONTAINS 
    4646    
    47    SUBROUTINE dyn_zdf( kt ) 
     47   SUBROUTINE dyn_zdf( kt, ktlev1, ktlev2, ktlev3, kt2lev, pu_rhs, pv_rhs ) 
    4848      !!---------------------------------------------------------------------- 
    4949      !!                  ***  ROUTINE dyn_zdf  *** 
    … …  
    5454      !! 
    5555      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing 
    56       !!         ua =         ub + 2*dt *       ua             vector form or linear free surf. 
    57       !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise 
     56      !!         pu_rhs =         uu(:,:,:,ktlev1) + 2*dt *       pu_rhs             vector form or linear free surf. 
     57      !!         pu_rhs = ( e3u_b*uu(:,:,:,ktlev1) + 2*dt * e3u_n*pu_rhs ) / e3u(:,:,:,ktlev3)   otherwise 
    5858      !!               - update the after velocity with the implicit vertical mixing. 
    5959      !!      This requires to solver the following system:  
    60       !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 
     60      !!         pu_rhs = pu_rhs + 1/e3u(:,:,:,ktlev3) dk+1[ mi(avm) / e3uw_a dk[ua] ] 
    6161      !!      with the following surface/top/bottom boundary condition: 
    6262      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
    6363      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
    6464      !! 
    65       !! ** Action :   (ua,va)   after velocity  
     65      !! ** Action :   (pu_rhs,pv_rhs)   after velocity  
    6666      !!--------------------------------------------------------------------- 
    67       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     67      INTEGER, INTENT(in) ::   kt                       ! ocean time-step index 
     68      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms 
     69      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms 
     70      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends -> momentum after fields 
    6871      ! 
    6972      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    … …  
    9699      ! 
    97100      !                             !* explicit top/bottom drag case 
    98       IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va) 
     101      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs )  ! add top/bottom friction trend to (pu_rhs,pv_rhs) 
    99102      ! 
    100103      ! 
    101104      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    102105         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
    103          ztrdu(:,:,:) = ua(:,:,:) 
    104          ztrdv(:,:,:) = va(:,:,:) 
    105       ENDIF 
    106       ! 
    107       !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va) 
     106         ztrdu(:,:,:) = pu_rhs(:,:,:) 
     107         ztrdv(:,:,:) = pv_rhs(:,:,:) 
     108      ENDIF 
     109      ! 
     110      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in pu_rhs,pv_rhs) 
    108111      ! 
    109112      !                    ! time stepping except vertical diffusion 
    110113      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
    111114         DO jk = 1, jpkm1 
    112             ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    113             va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     115            pu_rhs(:,:,jk) = ( uu(:,:,jk,ktlev1) + r2dt * pu_rhs(:,:,jk) ) * umask(:,:,jk) 
     116            pv_rhs(:,:,jk) = ( vv(:,:,jk,ktlev1) + r2dt * pv_rhs(:,:,jk) ) * vmask(:,:,jk) 
    114117         END DO 
    115118      ELSE                                      ! applied on thickness weighted velocity 
    116119         DO jk = 1, jpkm1 
    117             ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
    118                &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
    119             va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
    120                &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
     120            pu_rhs(:,:,jk) = (         e3u(:,:,jk,ktlev1) * uu(:,:,jk,ktlev1)  & 
     121               &          + r2dt * e3u(:,:,jk,ktlev2) * pu_rhs(:,:,jk)  ) / e3u(:,:,jk,ktlev3) * umask(:,:,jk) 
     122            pv_rhs(:,:,jk) = (         e3v(:,:,jk,ktlev1) * vv(:,:,jk,ktlev1)  & 
     123               &          + r2dt * e3v(:,:,jk,ktlev2) * pv_rhs(:,:,jk)  ) / e3v(:,:,jk,ktlev3) * vmask(:,:,jk) 
    121124         END DO 
    122125      ENDIF 
    … …  
    125128      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
    126129      !                not lead to the effective stress seen over the whole barotropic loop.  
    127       !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
     130      !     G. Madec : in linear free surface, e3u(:,:,:,ktlev3) = e3u(:,:,:,ktlev2) = e3u_0, so systematic use of e3u(:,:,:,ktlev3) 
    128131      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    129132         DO jk = 1, jpkm1        ! remove barotropic velocities 
    130             ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 
    131             va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
     133            pu_rhs(:,:,jk) = ( pu_rhs(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 
     134            pv_rhs(:,:,jk) = ( pv_rhs(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
    132135         END DO 
    133136         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
    … …  
    135138               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    136139               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    137                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    138                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    139                ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    140                va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     140               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) 
     141               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) 
     142               pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     143               pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
    141144            END DO 
    142145         END DO 
    … …  
    146149                  iku = miku(ji,jj)         ! top ocean level at u- and v-points  
    147150                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    148                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    149                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    150                   ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    151                   va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     151                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) 
     152                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) 
     153                  pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     154                  pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    152155               END DO 
    153156            END DO 
    … …  
    165168               DO jj = 2, jpjm1  
    166169                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     170                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3)   ! after scale factor at U-point 
    168171                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    169                         &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     172                        &         / ( ze3ua * e3uw(ji,jj,jk  ,kt2lev) ) * wumask(ji,jj,jk  ) 
    170173                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    171                         &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     174                        &         / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 
    172175                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    173176                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
    … …  
    182185               DO jj = 2, jpjm1  
    183186                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    184                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    185                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    186                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     187                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3)   ! after scale factor at U-point 
     188                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,kt2lev) ) * wumask(ji,jj,jk  ) 
     189                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 
    187190                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    188191                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
    … …  
    197200            DO ji = fs_2, fs_jpim1   ! vector opt. 
    198201               zwi(ji,jj,1) = 0._wp 
    199                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
    200                zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 
     202               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3) 
     203               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,kt2lev) ) * wumask(ji,jj,2) 
    201204               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) 
    202205               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
    … …  
    210213               DO jj = 2, jpjm1  
    211214                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    212                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
     215                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3)   ! after scale factor at U-point 
    213216                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
    214                         &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     217                        &         / ( ze3ua * e3uw(ji,jj,jk  ,kt2lev) ) * wumask(ji,jj,jk  ) 
    215218                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    216                         &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     219                        &         / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 
    217220                     zwi(ji,jj,jk) = zzwi 
    218221                     zws(ji,jj,jk) = zzws 
    … …  
    225228               DO jj = 2, jpjm1  
    226229                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    227                      ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point 
    228                      zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    229                      zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     230                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3)   ! after scale factor at U-point 
     231                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,kt2lev) ) * wumask(ji,jj,jk  ) 
     232                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 
    230233                     zwi(ji,jj,jk) = zzwi 
    231234                     zws(ji,jj,jk) = zzws 
    … …  
    254257            DO ji = 2, jpim1 
    255258               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
    256                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
     259               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3)   ! after scale factor at T-point 
    257260               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    258261            END DO 
    … …  
    263266                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
    264267                  iku = miku(ji,jj)       ! ocean top level at u- and v-points  
    265                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
     268                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3)   ! after scale factor at T-point 
    266269                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
    267270               END DO 
    … …  
    282285      !   m is decomposed in the product of an upper and a lower triangular matrix 
    283286      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    284       !   The solution (the after velocity) is in ua 
     287      !   The solution (the after velocity) is in pu_rhs 
    285288      !----------------------------------------------------------------------- 
    286289      ! 
    … …  
    295298      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    296299         DO ji = fs_2, fs_jpim1   ! vector opt. 
    297             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
    298             ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     300            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3)  
     301            pu_rhs(ji,jj,1) = pu_rhs(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    299302               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    300303         END DO 
    … …  
    303306         DO jj = 2, jpjm1 
    304307            DO ji = fs_2, fs_jpim1 
    305                ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
     308               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pu_rhs(ji,jj,jk-1) 
    306309            END DO 
    307310         END DO 
    … …  
    310313      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
    311314         DO ji = fs_2, fs_jpim1   ! vector opt. 
    312             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     315            pu_rhs(ji,jj,jpkm1) = pu_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    313316         END DO 
    314317      END DO 
    … …  
    316319         DO jj = 2, jpjm1 
    317320            DO ji = fs_2, fs_jpim1 
    318                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     321               pu_rhs(ji,jj,jk) = ( pu_rhs(ji,jj,jk) - zws(ji,jj,jk) * pu_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    319322            END DO 
    320323         END DO 
    … …  
    331334               DO jj = 2, jpjm1  
    332335                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    333                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     336                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3)   ! after scale factor at V-point 
    334337                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    335                         &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     338                        &         / ( ze3va * e3vw(ji,jj,jk  ,kt2lev) ) * wvmask(ji,jj,jk  ) 
    336339                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    337                         &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     340                        &         / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 
    338341                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
    339342                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
    … …  
    348351               DO jj = 2, jpjm1  
    349352                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    350                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    351                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    352                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     353                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3)   ! after scale factor at V-point 
     354                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,kt2lev) ) * wvmask(ji,jj,jk  ) 
     355                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 
    353356                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
    354357                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
    … …  
    363366            DO ji = fs_2, fs_jpim1   ! vector opt. 
    364367               zwi(ji,jj,1) = 0._wp 
    365                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
    366                zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 
     368               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3) 
     369               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,kt2lev) ) * wvmask(ji,jj,2) 
    367370               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) 
    368371               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
    … …  
    376379               DO jj = 2, jpjm1    
    377380                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    378                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
     381                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3)   ! after scale factor at V-point 
    379382                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
    380                         &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     383                        &         / ( ze3va * e3vw(ji,jj,jk  ,kt2lev) ) * wvmask(ji,jj,jk  ) 
    381384                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    382                         &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     385                        &         / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 
    383386                     zwi(ji,jj,jk) = zzwi 
    384387                     zws(ji,jj,jk) = zzws 
    … …  
    391394               DO jj = 2, jpjm1    
    392395                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    393                      ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point 
    394                      zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    395                      zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     396                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3)   ! after scale factor at V-point 
     397                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,kt2lev) ) * wvmask(ji,jj,jk  ) 
     398                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 
    396399                     zwi(ji,jj,jk) = zzwi 
    397400                     zws(ji,jj,jk) = zzws 
    … …  
    419422            DO ji = 2, jpim1 
    420423               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    421                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
     424               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3)   ! after scale factor at T-point 
    422425               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    423426            END DO 
    … …  
    427430               DO ji = 2, jpim1 
    428431                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    429                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
     432                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3)   ! after scale factor at T-point 
    430433                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
    431434               END DO 
    … …  
    446449      !   m is decomposed in the product of an upper and lower triangular matrix 
    447450      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    448       !   The solution (after velocity) is in 2d array va 
     451      !   The solution (after velocity) is in 2d array pv_rhs 
    449452      !----------------------------------------------------------------------- 
    450453      ! 
    … …  
    459462      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    460463         DO ji = fs_2, fs_jpim1   ! vector opt.           
    461             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
    462             va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     464            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3)  
     465            pv_rhs(ji,jj,1) = pv_rhs(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    463466               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
    464467         END DO 
    … …  
    467470         DO jj = 2, jpjm1 
    468471            DO ji = fs_2, fs_jpim1   ! vector opt. 
    469                va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
     472               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pv_rhs(ji,jj,jk-1) 
    470473            END DO 
    471474         END DO 
    … …  
    474477      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
    475478         DO ji = fs_2, fs_jpim1   ! vector opt. 
    476             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     479            pv_rhs(ji,jj,jpkm1) = pv_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    477480         END DO 
    478481      END DO 
    … …  
    480483         DO jj = 2, jpjm1 
    481484            DO ji = fs_2, fs_jpim1 
    482                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     485               pv_rhs(ji,jj,jk) = ( pv_rhs(ji,jj,jk) - zws(ji,jj,jk) * pv_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    483486            END DO 
    484487         END DO 
    … …  
    486489      ! 
    487490      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    488          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    489          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
     491         ztrdu(:,:,:) = ( pu_rhs(:,:,:) - uu(:,:,:,ktlev1) ) / r2dt - ztrdu(:,:,:) 
     492         ztrdv(:,:,:) = ( pv_rhs(:,:,:) - vv(:,:,:,ktlev1) ) / r2dt - ztrdv(:,:,:) 
    490493         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    491494         DEALLOCATE( ztrdu, ztrdv )  

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90

    r10425 r10825  
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_zdf( kt ) 
     46   SUBROUTINE tra_zdf( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs ) 
    4747      !!---------------------------------------------------------------------- 
    4848      !!                  ***  ROUTINE tra_zdf  *** 
    … …  
    5151      !!--------------------------------------------------------------------- 
    5252      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     53      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms 
     54      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms 
     55      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    5356      ! 
    5457      INTEGER  ::   jk   ! Dummy loop indices 
    … …  
    7073      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
    7174         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    72          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    73          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     75         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
     76         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
    7477      ENDIF 
    7578      ! 
    7679      !                                      !* compute lateral mixing trend and add it to the general trend 
    77       CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )  
     80      CALL tra_zdf_imp( kt, nit000, ktlev1, ktlev2, ktlev3, kt2lev, 'TRA', r2dt, ts(:,:,:,:,ktlev1), pts_rhs, jpts )  
    7881 
    7982!!gm WHY here !   and I don't like that ! 
    … …  
    8184      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    8285      ! JMM : restore negative salinities to small salinities: 
    83       WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
     86      WHERE( pts_rhs(:,:,:,jp_sal) < 0._wp )   pts_rhs(:,:,:,jp_sal) = 0.1_wp 
    8487!!gm 
    8588 
    8689      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    8790         DO jk = 1, jpkm1 
    88             ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 
    89                &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 
    90             ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 
    91               &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 
     91            ztrdt(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_tem)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_tem,ktlev1)*e3t(:,:,jk,ktlev1) ) & 
     92               &          / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrdt(:,:,jk) 
     93            ztrds(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_sal)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_sal,ktlev1)*e3t(:,:,jk,ktlev1) ) & 
     94              &           / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrds(:,:,jk) 
    9295         END DO 
    9396!!gm this should be moved in trdtra.F90 and done on all trends 
    … …  
    107110 
    108111  
    109    SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )  
     112   SUBROUTINE tra_zdf_imp( kt, kit000, ktlev1, ktlev2, ktlev3, kt2lev, cdtype, p2dt, pt, pt_rhs, kjpt )  
    110113      !!---------------------------------------------------------------------- 
    111114      !!                  ***  ROUTINE tra_zdf_imp  *** 
    … …  
    125128      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
    126129      !! 
    127       !! ** Action  : - pta  becomes the after tracer 
     130      !! ** Action  : - pt_rhs  becomes the after tracer 
    128131      !!--------------------------------------------------------------------- 
    129132      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    130133      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index 
     134      INTEGER                              , INTENT(in   ) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms 
     135      INTEGER                              , INTENT(in   ) ::   kt2lev                   ! time level index for 2-time-level source terms 
    131136      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    132137      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    133138      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step 
    134       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields 
    135       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field 
     139      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt      ! before and now tracer fields 
     140      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! in: tracer trend ; out: after tracer field 
    136141      ! 
    137142      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    … …  
    181186                  DO jj = 2, jpjm1 
    182187                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
    183                         zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
    184                         zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
    185                         zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   & 
     188                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,kt2lev) 
     189                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev) 
     190                        zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zzwi - zzws   & 
    186191                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
    187192                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
    … …  
    194199                  DO jj = 2, jpjm1 
    195200                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    196                         zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk) 
    197                         zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
    198                         zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     201                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,kt2lev) 
     202                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev) 
     203                        zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    199204                    END DO 
    200205                  END DO 
    … …  
    216221            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal 
    217222            !   and "superior" (above diagonal) components of the tridiagonal system. 
    218             !   The solution will be in the 4d array pta. 
     223            !   The solution will be in the 4d array pt_rhs. 
    219224            !   The 3d array zwt is used as a work space array. 
    220             !   En route to the solution pta is used a to evaluate the rhs and then  
     225            !   En route to the solution pt_rhs is used a to evaluate the rhs and then  
    221226            !   used as a work space array: its value is modified. 
    222227            ! 
    … …  
    238243         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    239244            DO ji = fs_2, fs_jpim1 
    240                pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 
     245               pt_rhs(ji,jj,1,jn) = e3t(ji,jj,1,ktlev1) * pt(ji,jj,1,jn) + p2dt * e3t(ji,jj,1,ktlev2) * pt_rhs(ji,jj,1,jn) 
    241246            END DO 
    242247         END DO 
    … …  
    244249            DO jj = 2, jpjm1 
    245250               DO ji = fs_2, fs_jpim1 
    246                   zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side 
    247                   pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
     251                  zrhs = e3t(ji,jj,jk,ktlev1) * pt(ji,jj,jk,jn) + p2dt * e3t(ji,jj,jk,ktlev2) * pt_rhs(ji,jj,jk,jn)   ! zrhs=right hand side 
     252                  pt_rhs(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt_rhs(ji,jj,jk-1,jn) 
    248253               END DO 
    249254            END DO 
    … …  
    252257         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    253258            DO ji = fs_2, fs_jpim1 
    254                pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     259               pt_rhs(ji,jj,jpkm1,jn) = pt_rhs(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    255260            END DO 
    256261         END DO 
    … …  
    258263            DO jj = 2, jpjm1 
    259264               DO ji = fs_2, fs_jpim1 
    260                   pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
     265                  pt_rhs(ji,jj,jk,jn) = ( pt_rhs(ji,jj,jk,jn) - zws(ji,jj,jk) * pt_rhs(ji,jj,jk+1,jn) )   & 
    261266                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    262267               END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10806 r10825  
    203203      ENDIF 
    204204       
    205                          CALL dyn_zdf       ( kstp )  ! vertical diffusion 
     205                         CALL dyn_zdf       ( kstp, Nm1, Nnn, Np1, Nnn_2lev, uu(:,:,:,Np1), vv(:,:,:,Np1) )  ! vertical diffusion 
    206206 
    207207      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    256256      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics 
    257257!!gm 
    258                          CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields 
     258                         CALL tra_zdf       ( kstp, Nm1, Nnn, Np1, Nnn_2lev, ts(:,:,:,:,Np1) )  ! vertical mixing and after tracer fields 
    259259      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection 
    260260 
Note: See TracChangeset for help on using the changeset viewer.