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 10884 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps – NEMO

Ignore:
Timestamp:
2019-04-18T17:00:30+02:00 (5 years ago)
Author:
davestorkey
Message:

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : dynzdf, trazdf, trczdf

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

Legend:

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

    r10874 r10884  
    4545CONTAINS 
    4646    
    47    SUBROUTINE dyn_zdf( kt ) 
     47   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa ) 
    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      !!         puu(:,:,:,Kaa) =         puu(:,:,:,Kbb) + 2*dt *       puu(:,:,:,Krhs)             vector form or linear free surf. 
     57      !!         puu(:,:,:,Kaa) = ( e3u_b*puu(:,:,:,Kbb) + 2*dt * e3u_n*puu(:,:,:,Krhs) ) / e3u(:,:,:,Kaa)   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      !!         puu(:,:,:,Kaa) = puu(:,:,:,Kaa) + 1/e3u(:,:,:,Kaa) dk+1[ mi(avm) / e3uw(:,:,:,Kaa) 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 :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   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 )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
     69      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
    6870      ! 
    6971      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
     
    9698      ! 
    9799      !                             !* 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) 
     100      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 
    99101      ! 
    100102      ! 
    101103      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    102104         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) 
     105         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     106         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
     107      ENDIF 
     108      ! 
     109      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 
    108110      ! 
    109111      !                    ! time stepping except vertical diffusion 
    110112      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
    111113         DO jk = 1, jpkm1 
    112             ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    113             va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     114            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk) 
     115            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk) 
    114116         END DO 
    115117      ELSE                                      ! applied on thickness weighted velocity 
    116118         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) 
     119            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  & 
     120               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 
     121            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  & 
     122               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 
    121123         END DO 
    122124      ENDIF 
     
    125127      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
    126128      !                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 
     129      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 
    128130      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    129131         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) 
     132            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - ua_b(:,:) ) * umask(:,:,jk) 
     133            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - va_b(:,:) ) * vmask(:,:,jk) 
    132134         END DO 
    133135         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
     
    135137               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    136138               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 
     139               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
     140               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
     141               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     142               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
    141143            END DO 
    142144         END DO 
     
    146148                  iku = miku(ji,jj)         ! top ocean level at u- and v-points  
    147149                  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 
     150                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
     151                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
     152                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     153                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    152154               END DO 
    153155            END DO 
     
    165167               DO jj = 2, jpjm1  
    166168                  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 
     169                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
    168170                     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  ) 
     171                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    170172                     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) 
     173                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    172174                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    173175                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     
    182184               DO jj = 2, jpjm1  
    183185                  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) 
     186                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     187                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     188                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    187189                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) 
    188190                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 
     
    197199            DO ji = fs_2, fs_jpim1   ! vector opt. 
    198200               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) 
     201               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
     202               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 
    201203               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) 
    202204               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 
     
    210212               DO jj = 2, jpjm1  
    211213                  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 
     214                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
    213215                     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  ) 
     216                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
    215217                     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) 
     218                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    217219                     zwi(ji,jj,jk) = zzwi 
    218220                     zws(ji,jj,jk) = zzws 
     
    225227               DO jj = 2, jpjm1  
    226228                  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) 
     229                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point 
     230                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  ) 
     231                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 
    230232                     zwi(ji,jj,jk) = zzwi 
    231233                     zws(ji,jj,jk) = zzws 
     
    254256            DO ji = 2, jpim1 
    255257               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 
     258               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
    257259               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    258260            END DO 
     
    263265                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
    264266                  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 
     267                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point 
    266268                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
    267269               END DO 
     
    282284      !   m is decomposed in the product of an upper and a lower triangular matrix 
    283285      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    284       !   The solution (the after velocity) is in ua 
     286      !   The solution (the after velocity) is in puu(:,:,:,Kaa) 
    285287      !----------------------------------------------------------------------- 
    286288      ! 
     
    295297      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    296298         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) )   & 
     299            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)  
     300            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    299301               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    300302         END DO 
     
    303305         DO jj = 2, jpjm1 
    304306            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) 
     307               puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) 
    306308            END DO 
    307309         END DO 
     
    310312      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
    311313         DO ji = fs_2, fs_jpim1   ! vector opt. 
    312             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     314            puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
    313315         END DO 
    314316      END DO 
     
    316318         DO jj = 2, jpjm1 
    317319            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) 
     320               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
    319321            END DO 
    320322         END DO 
     
    331333               DO jj = 2, jpjm1  
    332334                  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 
     335                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    334336                     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  ) 
     337                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    336338                     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) 
     339                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    338340                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
    339341                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     
    348350               DO jj = 2, jpjm1  
    349351                  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) 
     352                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     353                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     354                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    353355                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  ) 
    354356                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) 
     
    363365            DO ji = fs_2, fs_jpim1   ! vector opt. 
    364366               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) 
     367               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
     368               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 
    367369               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) 
    368370               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 
     
    376378               DO jj = 2, jpjm1    
    377379                  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 
     380                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
    379381                     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  ) 
     382                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
    381383                     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) 
     384                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    383385                     zwi(ji,jj,jk) = zzwi 
    384386                     zws(ji,jj,jk) = zzws 
     
    391393               DO jj = 2, jpjm1    
    392394                  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) 
     395                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point 
     396                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  ) 
     397                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 
    396398                     zwi(ji,jj,jk) = zzwi 
    397399                     zws(ji,jj,jk) = zzws 
     
    419421            DO ji = 2, jpim1 
    420422               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 
     423               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
    422424               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    423425            END DO 
     
    427429               DO ji = 2, jpim1 
    428430                  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 
     431                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point 
    430432                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
    431433               END DO 
     
    446448      !   m is decomposed in the product of an upper and lower triangular matrix 
    447449      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    448       !   The solution (after velocity) is in 2d array va 
     450      !   The solution (after velocity) is in 2d array pvv(:,:,:,Kaa) 
    449451      !----------------------------------------------------------------------- 
    450452      ! 
     
    459461      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    460462         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) )   & 
     463            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)  
     464            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    463465               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
    464466         END DO 
     
    467469         DO jj = 2, jpjm1 
    468470            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) 
     471               pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) 
    470472            END DO 
    471473         END DO 
     
    474476      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
    475477         DO ji = fs_2, fs_jpim1   ! vector opt. 
    476             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     478            pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
    477479         END DO 
    478480      END DO 
     
    480482         DO jj = 2, jpjm1 
    481483            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) 
     484               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 
    483485            END DO 
    484486         END DO 
     
    486488      ! 
    487489      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    488          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    489          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
     490         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:) 
     491         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:) 
    490492         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    491493         DEALLOCATE( ztrdu, ztrdv )  
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90

    r10880 r10884  
    9292      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
    9393      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean transport components 
    94       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
     94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90

    r10874 r10884  
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_zdf( kt ) 
     46   SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 
    4747      !!---------------------------------------------------------------------- 
    4848      !!                  ***  ROUTINE tra_zdf  *** 
     
    5050      !! ** Purpose :   compute the vertical ocean tracer physics. 
    5151      !!--------------------------------------------------------------------- 
    52       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     52      INTEGER                                  , INTENT(in)    :: kt                  ! ocean time-step index 
     53      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs, Kaa ! time level indices 
     54      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts                 ! active tracers and RHS of tracer equation 
    5355      ! 
    5456      INTEGER  ::   jk   ! Dummy loop indices 
     
    7072      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
    7173         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    72          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    73          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     74         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
     75         ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 
    7476      ENDIF 
    7577      ! 
    7678      !                                      !* compute lateral mixing trend and add it to the general trend 
    77       CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )  
     79      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, pts, Kaa, jpts )  
    7880 
    7981!!gm WHY here !   and I don't like that ! 
     
    8183      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    8284      ! JMM : restore negative salinities to small salinities: 
    83       WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
     85      WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
    8486!!gm 
    8587 
    8688      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    8789         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) 
     90            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 
     91               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 
     92            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 
     93              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 
    9294         END DO 
    9395!!gm this should be moved in trdtra.F90 and done on all trends 
     
    107109 
    108110  
    109    SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )  
     111   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, pt, Kaa, kjpt )  
    110112      !!---------------------------------------------------------------------- 
    111113      !!                  ***  ROUTINE tra_zdf_imp  *** 
     
    125127      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
    126128      !! 
    127       !! ** Action  : - pta  becomes the after tracer 
    128       !!--------------------------------------------------------------------- 
    129       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    130       INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index 
    131       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    132       INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    133       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 
     129      !! ** Action  : - pt(:,:,:,:,Kaa)  becomes the after tracer 
     130      !!--------------------------------------------------------------------- 
     131      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index 
     132      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     133      INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index 
     134      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     135      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
     136      REAL(wp)                                 , INTENT(in   ) ::   p2dt     ! tracer time-step 
     137      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt       ! tracers and RHS of tracer equation 
    136138      ! 
    137139      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     
    181183                  DO jj = 2, jpjm1 
    182184                     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   & 
     185                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
     186                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     187                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
    186188                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
    187189                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
     
    194196                  DO jj = 2, jpjm1 
    195197                     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) 
     198                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
     199                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     200                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    199201                    END DO 
    200202                  END DO 
     
    218220            !   The solution will be in the 4d array pta. 
    219221            !   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  
     222            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then  
    221223            !   used as a work space array: its value is modified. 
    222224            ! 
     
    238240         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    239241            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) 
     242               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Kaa) 
    241243            END DO 
    242244         END DO 
     
    244246            DO jj = 2, jpjm1 
    245247               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) 
     248                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Kaa)   ! zrhs=right hand side 
     249                  pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
    248250               END DO 
    249251            END DO 
     
    252254         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    253255            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) 
     256               pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    255257            END DO 
    256258         END DO 
     
    258260            DO jj = 2, jpjm1 
    259261               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) )   & 
     262                  pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
    261263                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    262264               END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10883 r10884  
    208208      ENDIF 
    209209       
    210                          CALL dyn_zdf       ( kstp )  ! vertical diffusion 
     210                         CALL dyn_zdf( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion     ==> after 
    211211 
    212212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    261261      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics 
    262262!!gm 
    263                          CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields 
     263                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer   ==> after 
    264264      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection 
    265265 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90

    r10880 r10884  
    7777         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge 
    7878#endif 
    79                                 CALL trc_zdf    ( kt )      ! vertical mixing and after tracer fields 
     79                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer ==> after 
    8080                                CALL trc_nxt    ( kt )      ! tracer fields at next time step      
    8181         IF( ln_trcrad )        CALL trc_rad    ( kt )      ! Correct artificial negative concentrations 
     
    8686                                CALL trc_sbc( kt )            ! surface boundary condition 
    8787         IF( ln_trcdmp )        CALL trc_dmp( kt )            ! internal damping trends 
    88                                 CALL trc_zdf( kt )            ! vertical mixing and after tracer fields 
     88                                CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer ==> after 
    8989                                CALL trc_nxt( kt )            ! tracer fields at next time step      
    9090          IF( ln_trcrad )       CALL trc_rad( kt )            ! Correct artificial negative concentrations 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90

    r10068 r10884  
    3636CONTAINS 
    3737 
    38    SUBROUTINE trc_zdf( kt ) 
     38   SUBROUTINE trc_zdf( kt, Kbb, Kmm, Krhs, ptr, Kaa ) 
    3939      !!---------------------------------------------------------------------- 
    4040      !!                  ***  ROUTINE trc_zdf  *** 
     
    4343      !!              an implicit time-stepping scheme. 
    4444      !!--------------------------------------------------------------------- 
    45       INTEGER, INTENT( in ) ::  kt      ! ocean time-step index 
     45      INTEGER                                   , INTENT(in   ) ::   kt                   ! ocean time-step index 
     46      INTEGER                                   , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices 
     47      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr                  ! passive tracers and RHS of tracer equation 
    4648      ! 
    4749      INTEGER               ::  jk, jn 
     
    5254      IF( ln_timing )   CALL timing_start('trc_zdf') 
    5355      ! 
    54       IF( l_trdtrc )   ztrtrd(:,:,:,:)  = tra(:,:,:,:) 
     56      IF( l_trdtrc )   ztrtrd(:,:,:,:)  = ptr(:,:,:,:,Krhs) 
    5557      ! 
    56       CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, trb, tra, jptra )    !   implicit scheme           
     58      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, Kbb, Kmm, ptr, Kaa, jptra )    !   implicit scheme           
    5759      ! 
    5860      IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    5961         DO jn = 1, jptra 
    6062            DO jk = 1, jpkm1 
    61                ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dttrc ) - ztrtrd(:,:,jk,jn) 
     63               ztrtrd(:,:,jk,jn) = ( ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kbb) ) / r2dttrc ) - ztrtrd(:,:,jk,jn) 
    6264            END DO 
    6365            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) ) 
Note: See TracChangeset for help on using the changeset viewer.