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

Changeset 10883


Ignore:
Timestamp:
2019-04-18T14:29:58+02:00 (4 years ago)
Author:
davestorkey
Message:

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: some of the ocean physics modules.

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

Legend:

Unmodified
Added
Removed

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

    r10425 r10883  
    124124 
    125125 
    126    SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 
     126   SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    127127      !!---------------------------------------------------------------------- 
    128128      !!                   ***  ROUTINE zdf_gls  *** 
    … …  
    134134      !! 
    135135      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     136      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    136137      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    137138      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    … …  
    176177          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    177178          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0) 
    178           ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
    179              &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
     179          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  & 
     180             &                                         + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
    180181         END DO 
    181182      END DO 
    … …  
    185186               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    186187               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0) 
    187                ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    188                   &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
     188               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  & 
     189                  &                                         + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
    189190            END DO 
    190191         END DO 
    … …  
    222223            DO jj = 2, jpjm1  
    223224               DO ji = fs_2, fs_jpim1   ! vector opt. 
    224                   zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    225                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 
     225                  zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
     226                  zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) 
    226227                  zcoef = ( zup / MAX( zdown, rsmall ) ) 
    227228                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 
    … …  
    276277               zcof = rfact_tke * tmask(ji,jj,jk) 
    277278               !                                        ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 
    278                zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
     279               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    279280               !                                        ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 
    280                zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
     281               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    281282               !                                        ! diagonal 
    282283               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
    … …  
    306307      !  
    307308      ! One level below 
    308       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   & 
     309      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))   & 
    309310         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   ) 
    310311      zd_lw(:,:,2) = 0._wp  
    … …  
    325326      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag 
    326327      zd_lw(:,:,2) = 0._wp 
    327       zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
     328      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    328329      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    329           &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
    330 !!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    331       en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     330          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
     331!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 
     332      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
    332333      ! 
    333334      ! 
    … …  
    526527               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
    527528               !                                               ! lower diagonal 
    528                zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
     529               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    529530               !                                               ! upper diagonal 
    530                zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
     531               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    531532               !                                               ! diagonal 
    532533               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
    … …  
    554555         ! 
    555556         ! One level below 
    556          zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
    557          zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 
     557         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 
     558         zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 
    558559         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    559560         zd_lw(:,:,2) = 0._wp 
    … …  
    575576         ! 
    576577         ! Set psi vertical flux at the surface: 
    577          zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    578          zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     578         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 
     579         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 
    579580         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    580581         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    581             &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 
     582            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 
    582583         zflxs(:,:)   = zdep(:,:) * zflxs(:,:) 
    583          psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     584         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 
    584585         ! 
    585586      END SELECT 
    … …  
    607608               ! 
    608609               ! Just above last level, Dirichlet condition again (GOTM like) 
    609                zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 
     610               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) ) 
    610611               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
    611612               zd_lw(ji,jj,ibotm1) = 0._wp 
    … …  
    635636               ! 
    636637               ! Set psi vertical flux at the bottom: 
    637                zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 
     638               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm) 
    638639               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   & 
    639640                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
    640                psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 
     641               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm) 
    641642            END DO 
    642643         END DO 

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

    r10425 r10883  
    4242   !!---------------------------------------------------------------------- 
    4343   USE oce            ! ocean dynamics and active tracers 
    44                       ! uses wn from previous time step (which is now wb) to calculate hbl 
     44                      ! uses ww from previous time step (which is now wb) to calculate hbl 
    4545   USE dom_oce        ! ocean space and time domain 
    4646   USE zdf_oce        ! ocean vertical physics 
    … …  
    122122 
    123123 
    124    SUBROUTINE zdf_osm( kt, p_avm, p_avt ) 
     124   SUBROUTINE zdf_osm( kt, Kbb, Kmm, p_avm, p_avt ) 
    125125      !!---------------------------------------------------------------------- 
    126126      !!                   ***  ROUTINE zdf_osm  *** 
    … …  
    157157      !!         the equation number. (LMD94, here after) 
    158158      !!---------------------------------------------------------------------- 
    159       INTEGER                   , INTENT(in   ) ::   kt            ! ocean time step 
     159      INTEGER                   , INTENT(in   ) ::  kt             ! ocean time step 
     160      INTEGER                   , INTENT(in   ) ::  Kbb, Kmm       ! ocean time level indices 
    160161      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points) 
    161162      !! 
    … …  
    314315           zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 
    315316           ! Upwards surface salinity flux for non-local term 
    316            zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
     317           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
    317318           ! Non radiative upwards surface buoyancy flux 
    318319           zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
    … …  
    407408     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    408409     ! BL must be always 2 levels deep. 
    409       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,3) ) 
     410      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 
    410411      ibld(:,:) = 3 
    411412      DO jk = 4, jpkm1 
    412413         DO jj = 2, jpjm1 
    413414            DO ji = 2, jpim1 
    414                IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
     415               IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    415416                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 
    416417               ENDIF 
    … …  
    430431               zthick=0._wp 
    431432               DO jm = 2, ibld(ji,jj) 
    432                   zthick=zthick+e3t_n(ji,jj,jm) 
    433                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    434                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    435                   zu   = zu  + e3t_n(ji,jj,jm) & 
    436                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
     433                  zthick=zthick+e3t(ji,jj,jm,Kmm) 
     434                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     435                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     436                  zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     437                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    437438                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    438                   zv   = zv  + e3t_n(ji,jj,jm) & 
    439                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
     439                  zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     440                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    440441                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    441442               END DO 
    … …  
    444445               zu_bl(ji,jj) = zu / zthick 
    445446               zv_bl(ji,jj) = zv / zthick 
    446                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    447                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    448                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
     447               zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     448               zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     449               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    449450                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    450                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
     451               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    451452                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    452453               zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    … …  
    487488      ibld(:,:) = 3 
    488489 
    489       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
     490      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
    490491      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 
    491       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     492      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
    492493 
    493494      DO jk = 4, jpkm1 
    494495         DO jj = 2, jpjm1 
    495496            DO ji = 2, jpim1 
    496                IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN 
     497               IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 
    497498                  ibld(ji,jj) =  MIN(mbkt(ji,jj), jk) 
    498499               ENDIF 
    … …  
    520521 
    521522                  DO jk = imld(ji,jj), ibld(ji,jj) 
    522                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) & 
    523                           & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max 
    524  
    525                      zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) ) 
     523                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 
     524                          & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 
     525 
     526                     zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 
    526527                     zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    527528 
    528                      IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1 
     529                     IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
    529530                  END DO 
    530531                  hbl(ji,jj) = zhbl_s 
    … …  
    534535! stable 
    535536                  DO jk = imld(ji,jj), ibld(ji,jj) 
    536                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          & 
    537                           &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) & 
     537                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          & 
     538                          &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 
    538539                          & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 
    539540 
    … …  
    543544                          &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   & 
    544545                          &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) & 
    545                           &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here 
     546                          &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here 
    546547 
    547548                     zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 
    548                      IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1 
     549                     IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
    549550                  END DO 
    550                   hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) ) 
     551                  hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 
    551552                  ibld(ji,jj) = MAX(jm, 3 ) 
    552553                  IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    … …  
    558559                  hbli(ji,jj) = hbl(ji,jj) 
    559560               ELSE 
    560                   hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) ) 
     561                  hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 
    561562                  IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 
    562563               ENDIF 
    563564            ENDIF 
    564             zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     565            zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    565566         END DO 
    566567      END DO 
    … …  
    581582               zthick=0._wp 
    582583               DO jm = 2, ibld(ji,jj) 
    583                   zthick=zthick+e3t_n(ji,jj,jm) 
    584                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    585                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    586                   zu   = zu  + e3t_n(ji,jj,jm) & 
    587                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
     584                  zthick=zthick+e3t(ji,jj,jm,Kmm) 
     585                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     586                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     587                  zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     588                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    588589                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    589                   zv   = zv  + e3t_n(ji,jj,jm) & 
    590                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
     590                  zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     591                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    591592                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    592593               END DO 
    … …  
    595596               zu_bl(ji,jj) = zu / zthick 
    596597               zv_bl(ji,jj) = zv / zthick 
    597                zdt_bl(ji,jj) = zt_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    598                zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    599                zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
     598               zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     599               zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     600               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    600601                      &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    601                zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
     602               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    602603                      &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    603604               zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 
    604                zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) 
     605               zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    605606               IF ( lconv(ji,jj) ) THEN 
    606607                  IF ( zdb_bl(ji,jj) > 0._wp )THEN 
    … …  
    616617                        zwb_ent(ji,jj) = 0._wp 
    617618                     ENDIF 
    618                      inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     619                     inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 
    619620                     imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    620                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     621                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    621622                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    622623                  ELSE  ! IF (zdb_bl) 
    623624                     imld(ji,jj) = ibld(ji,jj) - 1 
    624                      zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     625                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    625626                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    626627                  ENDIF 
    … …  
    632633                        zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 
    633634                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 ) 
    634                         inhml = MAX( INT( zari * zhbl(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 ) 
     635                        inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 
    635636                        imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 
    636                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     637                        zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    637638                        zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    638639                     ELSE 
    639640                        imld(ji,jj) = ibld(ji,jj) - 1 
    640                         zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) 
     641                        zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 
    641642                        zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 
    642643                     ENDIF ! IF (zdb_bl > 0.0) 
    … …  
    665666               zthick=0._wp 
    666667               DO jm = 2, imld(ji,jj) 
    667                   zthick=zthick+e3t_n(ji,jj,jm) 
    668                   zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    669                   zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    670                   zu   = zu  + e3t_n(ji,jj,jm) & 
    671                      &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
     668                  zthick=zthick+e3t(ji,jj,jm,Kmm) 
     669                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     670                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     671                  zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     672                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    672673                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    673                   zv   = zv  + e3t_n(ji,jj,jm) & 
    674                      &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
     674                  zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     675                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    675676                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    676677               END DO 
    … …  
    679680               zu_ml(ji,jj) = zu / zthick 
    680681               zv_ml(ji,jj) = zv / zthick 
    681                zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    682                zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    683                zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
     682               zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     683               zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     684               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    684685                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    685                zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
     686               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    686687                     &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    687688               zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    … …  
    696697                  zthick=0._wp 
    697698                  DO jm = 2, imld(ji,jj) 
    698                      zthick=zthick+e3t_n(ji,jj,jm) 
    699                      zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem) 
    700                      zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal) 
    701                      zu   = zu  + e3t_n(ji,jj,jm) & 
    702                         &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) & 
     699                     zthick=zthick+e3t(ji,jj,jm,Kmm) 
     700                     zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 
     701                     zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 
     702                     zu   = zu  + e3t(ji,jj,jm,Kmm) & 
     703                        &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 
    703704                        &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 
    704                      zv   = zv  + e3t_n(ji,jj,jm) & 
    705                         &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) & 
     705                     zv   = zv  + e3t(ji,jj,jm,Kmm) & 
     706                        &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 
    706707                        &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 
    707708                  END DO 
    … …  
    710711                  zu_ml(ji,jj) = zu / zthick 
    711712                  zv_ml(ji,jj) = zv / zthick 
    712                   zdt_ml(ji,jj) = zt_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_tem) 
    713                   zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal) 
    714                   zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) & 
     713                  zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 
     714                  zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 
     715                  zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 
    715716                        &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 
    716                   zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) & 
     717                  zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 
    717718                        &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 
    718719                  zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) 
    … …  
    775776                   zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 
    776777                   DO jk = 2 , ibld(ji,jj) 
    777                       znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     778                      znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    778779                      zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    779780                      zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    … …  
    791792                         zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 
    792793                         DO jk = 2, ibld(ji,jj) 
    793                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     794                            znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    794795                            zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    795796                            zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 
    … …  
    801802                         zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 
    802803                         DO jk = 2, ibld(ji,jj) 
    803                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     804                            znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    804805                            zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    805806                            zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    … …  
    824825              & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    825826                DO jk = 2 , ibld(ji,jj)-1 
    826                    znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     827                   znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    827828                   zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    828829                   zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 
    … …  
    833834                zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 
    834835                DO jk = 2, ibld(ji,jj) 
    835                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     836                   znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    836837                   IF ( znd < 1.0 ) THEN 
    837838                      zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) 
    … …  
    880881             IF ( lconv(ji,jj) ) THEN 
    881882                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity 
    882                     zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
     883                    zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    883884                    ! 
    884885                    zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5 
    … …  
    890891                IF ( zdh(ji,jj) > 0._wp ) THEN 
    891892                   DO jk = imld(ji,jj)+1 , ibld(ji,jj) 
    892                        zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     893                       zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    893894                       ! 
    894895                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) 
    … …  
    897898                   END DO 
    898899                ENDIF 
    899                 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out 
    900                 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj)) 
     900                ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 
     901                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 
    901902                ! could be taken out, take account of entrainment represents as a diffusivity 
    902903                ! should remove w from here, represents entrainment 
    … …  
    904905             ! stable conditions 
    905906                DO jk = 2, ibld(ji,jj) 
    906                    zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     907                   zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    907908                   zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 
    908909                   zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) 
    … …  
    932933            IF ( lconv(ji,jj) ) THEN 
    933934              DO jk = 2, imld(ji,jj) 
    934                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     935                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    935936                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 
    936937                 ! 
    … …  
    967968             IF ( lconv(ji,jj) ) THEN 
    968969                DO jk = 2, imld(ji,jj) 
    969                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     970                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    970971                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   & 
    971972                        &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) & 
    … …  
    978979! Stable conditions 
    979980                DO jk = 2, ibld(ji,jj) ! corrected to ibld 
    980                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     981                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    981982                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       & 
    982983                        &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) 
    … …  
    10011002             IF (lconv(ji,jj) ) THEN 
    10021003                DO jk = 2, imld(ji,jj) 
    1003                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
     1004                   zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
    10041005                   ! calculate turbulent length scale 
    10051006                   zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           & 
    … …  
    10351036             IF ( lconv(ji,jj) ) THEN 
    10361037                DO jk = 2 , imld(ji,jj) 
    1037                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     1038                   zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    10381039                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     & 
    10391040                        &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   & 
    … …  
    10781079            ELSE 
    10791080               DO jk = 2, ibld(ji,jj) 
    1080                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
    1081                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1081                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
     1082                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    10821083                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 
    10831084               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) 
    … …  
    11041105             IF ( lconv(ji,jj) ) THEN 
    11051106               DO jk = 2, imld(ji,jj) 
    1106                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj) 
    1107                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     1107                  zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 
     1108                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    11081109                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 
    11091110                       & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) 
    … …  
    11141115             ELSE 
    11151116               DO jk = 2, ibld(ji,jj) 
    1116                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    1117                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj) 
     1117                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
     1118                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 
    11181119                  IF ( zznd_d <= 2.0 ) THEN 
    11191120                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & 
    … …  
    11411142            IF ( lconv(ji,jj) ) THEN 
    11421143               DO jk = 2, ibld(ji,jj) 
    1143                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1144                  znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    11441145                  IF ( znd >= 0.0 ) THEN 
    11451146                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) 
    … …  
    11521153            ELSE 
    11531154               DO jk = 2, ibld(ji,jj) 
    1154                   znd = ( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
     1155                  znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 
    11551156                  IF ( znd >= 0.0 ) THEN 
    11561157                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
    … …  
    11711172          DO ji = 2, jpim1 
    11721173             DO jk= 2, ibld(ji,jj) 
    1173                 znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
     1174                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 
    11741175                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    11751176                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
    … …  
    11941195               END DO 
    11951196               DO jk = imld(ji,jj), ibld(ji,jj) 
    1196                   znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj) 
     1197                  znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 
    11971198                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 
    11981199                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) 
    … …  
    12391240             DO jj = 1, jpjm1 
    12401241                DO ji = 1, jpim1   ! vector opt. 
    1241                    z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   & 
    1242                         &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) & 
    1243                         &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    1244                    z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   & 
    1245                         &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) & 
    1246                         &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
     1242                   z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   & 
     1243                        &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 
     1244                        &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 
     1245                   z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   & 
     1246                        &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) & 
     1247                        &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 
    12471248                END DO 
    12481249             END DO 
    … …  
    13641365 
    13651366 
    1366    SUBROUTINE zdf_osm_init 
     1367   SUBROUTINE zdf_osm_init( Kmm )  
    13671368     !!---------------------------------------------------------------------- 
    13681369     !!                  ***  ROUTINE zdf_osm_init  *** 
    … …  
    13761377     !! ** input   :   Namlist namosm 
    13771378     !!---------------------------------------------------------------------- 
     1379     INTEGER, INTENT(in)    :: Kmm ! time level index (middle) 
     1380     ! 
    13781381     INTEGER  ::   ios            ! local integer 
    13791382     INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    … …  
    14231426     IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 
    14241427 
    1425      call osm_rst( nit000, 'READ' ) !* read or initialize hbl 
     1428     call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 
    14261429 
    14271430     IF( ln_zdfddm) THEN 
    … …  
    15171520 
    15181521 
    1519    SUBROUTINE osm_rst( kt, cdrw ) 
     1522   SUBROUTINE osm_rst( kt, Kmm, cdrw ) 
    15201523     !!--------------------------------------------------------------------- 
    15211524     !!                   ***  ROUTINE osm_rst  *** 
    … …  
    15271530     !!---------------------------------------------------------------------- 
    15281531 
    1529      INTEGER, INTENT(in) :: kt 
     1532     INTEGER         , INTENT(in) ::   kt     ! ocean time step index 
     1533     INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index (middle) 
    15301534     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    15311535 
    … …  
    15451549        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. ) 
    15461550        IF( id1 > 0 ) THEN                       ! 'wn' exists; read 
    1547            CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios ) 
    1548            WRITE(numout,*) ' ===>>>> :  wn read from restart file' 
     1551           CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios ) 
     1552           WRITE(numout,*) ' ===>>>> :  ww read from restart file' 
    15491553        ELSE 
    1550            wn(:,:,:) = 0._wp 
    1551            WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially' 
     1554           ww(:,:,:) = 0._wp 
     1555           WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially' 
    15521556        END IF 
    15531557        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. ) 
    … …  
    15681572     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return 
    15691573        IF(lwp) WRITE(numout,*) '---- osm-rst ----' 
    1570          CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios ) 
     1574         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  , ldxios = lwxios ) 
    15711575         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios ) 
    15721576         CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios ) 
    … …  
    15801584     ALLOCATE( imld_rst(jpi,jpj) ) 
    15811585     ! w-level of the mixing and mixed layers 
    1582      CALL eos_rab( tsn, rab_n ) 
    1583      CALL bn2(tsn, rab_n, rn2) 
     1586     CALL eos_rab( ts(:,:,:,:,Kmm), rab_n ) 
     1587     CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2) 
    15841588     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point 
    15851589     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    … …  
    15911595           DO ji = 1, jpi 
    15921596              ikt = mbkt(ji,jj) 
    1593               hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk) 
     1597              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 
    15941598              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level 
    15951599           END DO 
    … …  
    16001604        DO ji = 1, jpi 
    16011605           iiki = imld_rst(ji,jj) 
    1602            hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
     1606           hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth 
    16031607        END DO 
    16041608     END DO 
    … …  
    16101614 
    16111615 
    1612    SUBROUTINE tra_osm( kt ) 
     1616   SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) 
    16131617      !!---------------------------------------------------------------------- 
    16141618      !!                  ***  ROUTINE tra_osm  *** 
    … …  
    16201624      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace 
    16211625      !!---------------------------------------------------------------------- 
    1622       INTEGER, INTENT(in) :: kt 
     1626      INTEGER                                  , INTENT(in)    :: kt        ! time step index 
     1627      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices 
     1628      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation 
     1629      ! 
    16231630      INTEGER :: ji, jj, jk 
    16241631      ! 
    … …  
    16301637 
    16311638      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    1632          ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    1633          ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     1639         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
     1640         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 
    16341641      ENDIF 
    16351642 
    … …  
    16381645         DO jj = 2, jpjm1 
    16391646            DO ji = 2, jpim1 
    1640                tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      & 
     1647               pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      & 
    16411648                  &                 - (  ghamt(ji,jj,jk  )  & 
    1642                   &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk) 
    1643                tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      & 
     1649                  &                    - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) 
     1650               pts(ji,jj,jk,jp_sal,Krhs) =  pts(ji,jj,jk,jp_sal,Krhs)                      & 
    16441651                  &                 - (  ghams(ji,jj,jk  )  & 
    1645                   &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
     1652                  &                    - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 
    16461653            END DO 
    16471654         END DO 
    … …  
    16511658      ! save the non-local tracer flux trends for diagnostic 
    16521659      IF( l_trdtra )   THEN 
    1653          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    1654          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     1660         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     1661         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
    16551662!!bug gm jpttdzdf ==> jpttosm 
    16561663         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    … …  
    16601667 
    16611668      IF(ln_ctl) THEN 
    1662          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   & 
    1663          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     1669         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm  - Ta: ', mask1=tmask,   & 
     1670         &             tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    16641671      ENDIF 
    16651672      ! 
    … …  
    16841691 
    16851692 
    1686    SUBROUTINE dyn_osm( kt ) 
     1693   SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) 
    16871694      !!---------------------------------------------------------------------- 
    16881695      !!                  ***  ROUTINE dyn_osm  *** 
    … …  
    16931700      !! ** Method  :   ??? 
    16941701      !!---------------------------------------------------------------------- 
    1695       INTEGER, INTENT(in) ::   kt   ! 
     1702      INTEGER                             , INTENT( in )  ::  kt          ! ocean time step index 
     1703      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     1704      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    16961705      ! 
    16971706      INTEGER :: ji, jj, jk   ! dummy loop indices 
    … …  
    17081717         DO jj = 2, jpjm1 
    17091718            DO ji = 2, jpim1 
    1710                ua(ji,jj,jk) =  ua(ji,jj,jk)                      & 
     1719               puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      & 
    17111720                  &                 - (  ghamu(ji,jj,jk  )  & 
    1712                   &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) 
    1713                va(ji,jj,jk) =  va(ji,jj,jk)                      & 
     1721                  &                    - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 
     1722               pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs)                      & 
    17141723                  &                 - (  ghamv(ji,jj,jk  )  & 
    1715                   &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) 
     1724                  &                    - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 
    17161725            END DO 
    17171726         END DO 

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

    r10874 r10883  
    6161CONTAINS 
    6262 
    63    SUBROUTINE zdf_phy_init 
     63   SUBROUTINE zdf_phy_init( Kmm ) 
    6464      !!---------------------------------------------------------------------- 
    6565      !!                  ***  ROUTINE zdf_phy_init  *** 
    … …  
    7070      !!                set horizontal shape and vertical profile of background mixing coef. 
    7171      !!---------------------------------------------------------------------- 
     72      INTEGER, INTENT(in)    :: Kmm ! time level index (middle) 
     73      ! 
    7274      INTEGER ::   jk            ! dummy loop indices 
    7375      INTEGER ::   ioptio, ios   ! local integers 
    … …  
    193195      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF 
    194196      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF 
    195       IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init   ;   ENDIF 
     197      IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF 
    196198      ! 
    197199      IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' ) 
    … …  
    218220 
    219221 
    220    SUBROUTINE zdf_phy( kt ) 
     222   SUBROUTINE zdf_phy( kt, Kbb, Kmm ) 
    221223      !!---------------------------------------------------------------------- 
    222224      !!                     ***  ROUTINE zdf_phy  *** 
    … …  
    230232      !!                bottom stress.....                               <<<<====verifier ! 
    231233      !!---------------------------------------------------------------------- 
    232       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     234      INTEGER, INTENT(in) ::   kt         ! ocean time-step index 
     235      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
    233236      ! 
    234237      INTEGER ::   ji, jj, jk   ! dummy loop indice 
    … …  
    254257      ! 
    255258      IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form) 
    256          CALL zdf_sh2( ub, vb, un, vn, avm_k,   &     ! <<== in 
    257             &                           zsh2    )     ! ==>> out : shear production 
     259         CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in 
     260            &                     zsh2    )     ! ==>> out : shear production 
    258261      ! 
    259262      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points 
    260       CASE( np_RIC )   ;   CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz 
    261       CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz 
    262       CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz 
    263       CASE( np_OSM )   ;   CALL zdf_osm( kt               , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz 
     263      CASE( np_RIC )   ;   CALL zdf_ric( kt,      Kmm, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz 
     264      CASE( np_TKE )   ;   CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz 
     265      CASE( np_GLS )   ;   CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz 
     266      CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm      , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz 
    264267!     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value) 
    265268!         ! avt_k and avm_k set one for all at initialisation phase 
    … …  
    318321         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
    319322         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' )  
    320          ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after wn has been updated 
     323         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 
    321324      ENDIF 
    322325      ! 

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

    r10068 r10883  
    112112 
    113113 
    114    SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt ) 
     114   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt ) 
    115115      !!---------------------------------------------------------------------- 
    116116      !!                 ***  ROUTINE zdfric  *** 
    … …  
    125125      !!                    avt = avm0 / (1 + rn_alp*ri) 
    126126      !!                with ri  = N^2 / dz(u)**2 
    127       !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 
     127      !!                         = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ] 
    128128      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 
    129129      !!                where ri is the before local Richardson number, 
    … …  
    152152      !!---------------------------------------------------------------------- 
    153153      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step 
    154       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdept          ! depth of t-point  [m] 
     154      INTEGER                   , INTENT(in   ) ::   Kmm            ! ocean time level index 
    155155      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    156156      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points) 
    … …  
    189189            DO jj = 2, jpjm1 
    190190               DO ji = 2, jpim1 
    191                   IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 
     191                  IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN 
    192192                     p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk) 
    193193                     p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk) 

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

    r10069 r10883  
    1111   !!   zdf_sh2       : compute mixing the shear production term of TKE 
    1212   !!---------------------------------------------------------------------- 
     13   USE oce 
    1314   USE dom_oce        ! domain: ocean 
    1415   ! 
    … …  
    2829CONTAINS 
    2930 
    30    SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2  )  
     31   SUBROUTINE zdf_sh2( Kbb, Kmm, p_avm, p_sh2  )  
    3132      !!---------------------------------------------------------------------- 
    3233      !!                   ***  ROUTINE zdf_sh2  *** 
    … …  
    4748      !! References :   Bruchard, OM 2002 
    4849      !! --------------------------------------------------------------------- 
    49       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities 
     50      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm             ! ocean time level indices 
    5051      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points) 
    5152      REAL(wp), DIMENSION(:,:,:) , INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points) 
    … …  
    5960            DO ji = 1, jpim1 
    6061               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 
    61                   &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) & 
    62                   &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk) 
     62                  &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) & 
     63                  &         * (   uu(ji,jj,jk-1,Kbb) -   uu(ji,jj,jk,Kbb) ) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 
    6364               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 
    64                   &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) & 
    65                   &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk) 
     65                  &         * (   vv(ji,jj,jk-1,Kmm) -   vv(ji,jj,jk,Kmm) ) & 
     66                  &         * (   vv(ji,jj,jk-1,Kbb) -   vv(ji,jj,jk,Kbb) ) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 
    6667            END DO 
    6768         END DO 

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

    r10874 r10883  
    109109 
    110110 
    111    SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 
     111   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    112112      !!---------------------------------------------------------------------- 
    113113      !!                   ***  ROUTINE zdf_tke  *** 
    … …  
    155155      !!---------------------------------------------------------------------- 
    156156      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     157      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    157158      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
    158159      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    159160      !!---------------------------------------------------------------------- 
    160161      ! 
    161       CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )   ! now tke (en) 
    162       ! 
    163       CALL tke_avn( gdepw_n, e3t_n, e3w_n,        p_avm, p_avt )   ! now avt, avm, dissl 
     162      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en) 
     163      ! 
     164      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl 
    164165      ! 
    165166  END SUBROUTINE zdf_tke 
    166167 
    167168 
    168    SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt ) 
     169   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) 
    169170      !!---------------------------------------------------------------------- 
    170171      !!                   ***  ROUTINE tke_tke  *** 
    … …  
    186187      USE zdf_oce , ONLY : en   ! ocean vertical physics 
    187188      !! 
    188       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! depth of w-points 
    189       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     189      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    190190      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term 
    191191      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    … …  
    243243               zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
    244244               !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
    245                zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  & 
    246                   &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
     245               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  & 
     246                  &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
    247247               en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
    248248            END DO 
    … …  
    254254                  zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
    255255                  !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
    256                   zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    257                      &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
     256                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  & 
     257                     &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
    258258                  en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface 
    259259               END DO 
    … …  
    268268         ! 
    269269         !                        !* total energy produce by LC : cumulative sum over jk 
    270          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 
     270         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 
    271271         DO jk = 2, jpk 
    272             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 
     272            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 
    273273         END DO 
    274274         !                        !* finite Langmuir Circulation depth 
    … …  
    286286         DO jj = 1, jpj  
    287287            DO ji = 1, jpi 
    288                zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 
     288               zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 
    289289            END DO 
    290290         END DO 
    … …  
    302302                  IF ( zfr_i(ji,jj) /= 0. ) THEN                
    303303                     ! vertical velocity due to LC    
    304                      IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
     304                     IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    305305                        !                                           ! vertical velocity due to LC 
    306                         zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     306                        zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
    307307                        !                                           ! TKE Langmuir circulation source term 
    308308                        en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    … …  
    342342               !                                   ! eddy coefficient (ensure numerical stability) 
    343343               zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
    344                   &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  ) 
     344                  &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    345345               zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
    346                   &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  ) 
     346                  &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    347347               ! 
    348348               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    … …  
    402402      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    403403!!gm BUG : in the exp  remove the depth of ssh !!! 
    404 !!gm       i.e. use gde3w in argument (pdepw) 
     404!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 
    405405       
    406406       
    … …  
    409409            DO jj = 2, jpjm1 
    410410               DO ji = fs_2, fs_jpim1   ! vector opt. 
    411                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
     411                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    412412                     &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    413413               END DO 
    … …  
    418418            DO ji = fs_2, fs_jpim1   ! vector opt. 
    419419               jk = nmln(ji,jj) 
    420                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
     420               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    421421                  &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    422422            END DO 
    … …  
    431431                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    432432                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    433                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
     433                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    434434                     &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    435435               END DO 
    … …  
    441441 
    442442 
    443    SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 
     443   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 
    444444      !!---------------------------------------------------------------------- 
    445445      !!                   ***  ROUTINE tke_avn  *** 
    … …  
    477477      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics 
    478478      !! 
    479       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! depth (w-points) 
    480       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     479      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices 
    481480      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    482481      ! 
    … …  
    526525      ! 
    527526 !!gm Not sure of that coding for ISF.... 
    528       ! where wmask = 0 set zmxlm == p_e3w 
     527      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 
    529528      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    530529         DO jk = 2, jpkm1 
    531530            DO jj = 2, jpjm1 
    532531               DO ji = fs_2, fs_jpim1   ! vector opt. 
    533                   zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    534                   &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 
     532                  zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   & 
     533                  &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 
    535534                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
    536                   zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
    537                   zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
     535                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
     536                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 
    538537               END DO 
    539538            END DO 
    … …  
    544543            DO jj = 2, jpjm1 
    545544               DO ji = fs_2, fs_jpim1   ! vector opt. 
    546                   zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
     545                  zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 
    547546                  zmxlm(ji,jj,jk) = zemxl 
    548547                  zmxld(ji,jj,jk) = zemxl 
    … …  
    555554            DO jj = 2, jpjm1 
    556555               DO ji = fs_2, fs_jpim1   ! vector opt. 
    557                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
     556                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
    558557               END DO 
    559558            END DO 
    … …  
    562561            DO jj = 2, jpjm1 
    563562               DO ji = fs_2, fs_jpim1   ! vector opt. 
    564                   zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
     563                  zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
    565564                  zmxlm(ji,jj,jk) = zemxl 
    566565                  zmxld(ji,jj,jk) = zemxl 
    … …  
    573572            DO jj = 2, jpjm1 
    574573               DO ji = fs_2, fs_jpim1   ! vector opt. 
    575                   zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
     574                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 
    576575               END DO 
    577576            END DO 
    … …  
    580579            DO jj = 2, jpjm1 
    581580               DO ji = fs_2, fs_jpim1   ! vector opt. 
    582                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
     581                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 
    583582               END DO 
    584583            END DO 

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

    r10876 r10883  
    448448 
    449449      !                                      ! Ocean physics 
    450                            CALL zdf_phy_init    ! Vertical physics 
     450                           CALL zdf_phy_init( Nnn )    ! Vertical physics 
    451451                                      
    452452      !                                         ! Lateral physics 

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

    r10880 r10883  
    137137 
    138138      !  VERTICAL PHYSICS 
    139                          CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
     139                         CALL zdf_phy( kstp, Nbb, Nnn )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
    140140 
    141141      !  LATERAL  PHYSICS 
    … …  
    196196                         CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    197197                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    198       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes 
     198      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
    199199                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    200200                         CALL dyn_spg       ( kstp )  ! surface pressure gradient 
    … …  
    253253#endif 
    254254                         CALL tra_adv( kstp, Nbb, Nnn      , ts, Nrhs )  ! hor. + vert. advection  ==> RHS 
    255       IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes 
     255      IF( ln_zdfosm  )   CALL tra_osm( kstp,      Nnn      , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS 
    256256      IF( lrst_oce .AND. ln_zdfosm ) & 
    257            &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts 
     257           &             CALL osm_rst( kstp, Nnn, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts 
    258258                         CALL tra_ldf       ( kstp )  ! lateral mixing 
    259259 
Note: See TracChangeset for help on using the changeset viewer.