Changeset 8055


Ignore:
Timestamp:
2017-05-20T13:49:38+02:00 (3 years ago)
Author:
gm
Message:

wrk_OMP: 2nd step: add OMP processor distribution in ZDF package

Location:
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r7990 r8055  
    3131 
    3232   !! * Substitutions 
    33 #  include "vectopt_loop_substitute.h90" 
     33#  include "domain_substitute.h90"    
    3434   !!---------------------------------------------------------------------- 
    3535   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     
    3939CONTAINS 
    4040 
    41    SUBROUTINE zdf_ddm( kt ) 
     41   SUBROUTINE zdf_ddm( ARG_2D, kt, p_avm, p_avt, p_avs ) 
    4242      !!---------------------------------------------------------------------- 
    4343      !!                  ***  ROUTINE zdf_ddm  *** 
     
    6969      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999. 
    7070      !!---------------------------------------------------------------------- 
    71       INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step 
     71      INTEGER, INTENT(in   ) ::   ARG_2D   ! inner domain start-end i-indices 
     72      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step indexocean time step 
     73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points) 
     74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points) 
     75      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points) 
    7276      ! 
    7377      INTEGER  ::   ji, jj , jk     ! dummy loop indices 
     
    7781      REAL(wp) ::   zavft, zavfs    !   -      - 
    7882      REAL(wp) ::   zavdt, zavds    !   -      - 
    79       REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
     83      REAL(wp), DIMENSION(WRK_2D) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
    8084      !!---------------------------------------------------------------------- 
    8185      ! 
     
    9195!!gm                            and many acces in memory 
    9296          
    93          DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
    94             DO ji = 1, jpi 
     97         DO jj = k_Jstr, k_Jend        !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==! 
     98            DO ji = k_Jstr, k_Iend 
    9599               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
    96100!!gm please, use e3w_n below  
     
    109113         END DO 
    110114 
    111          DO jj = 1, jpj                                     ! indicators: 
    112             DO ji = 1, jpi 
     115         DO jj = k_Jstr, k_Jend          !==  indicators  ==! 
     116            DO ji = k_Jstr, k_Iend 
    113117               ! stability indicator: msks=1 if rn2>0; 0 elsewhere 
    114118               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp 
     
    135139         END DO 
    136140         ! mask zmsk in order to have avt and avs masked 
    137          zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 
     141         zmsks(WRK_2D) = zmsks(WRK_2D) * wmask(WRK_2D,jk) 
    138142 
    139143 
     
    141145         ! ------------------ 
    142146         ! Constant eddy coefficient: reset to the background value 
    143          DO jj = 1, jpj 
    144             DO ji = 1, jpi 
     147         DO jj = k_Jstr, k_Jend 
     148            DO ji = k_Jstr, k_Iend 
    145149               zinr = 1._wp / zrau(ji,jj) 
    146150               ! salt fingering 
     
    154158                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    155159               ! add to the eddy viscosity coef. previously computed 
    156                avs(ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
    157                avt(ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt 
    158                avm(ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 
     160               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds 
     161               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt 
     162               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 
    159163            END DO 
    160164         END DO 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r7990 r8055  
    3131   PUBLIC   zdf_evd    ! called by step.F90 
    3232 
     33   !! * Substitutions 
     34#  include "domain_substitute.h90"    
    3335   !!---------------------------------------------------------------------- 
    3436   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     
    3840CONTAINS 
    3941 
    40    SUBROUTINE zdf_evd( kt ) 
     42   SUBROUTINE zdf_evd( ARG_2D, kt, p_avm, p_avt ) 
    4143      !!---------------------------------------------------------------------- 
    4244      !!                  ***  ROUTINE zdf_evd  *** 
     
    5557      !! ** Action  :   avt, avm   enhanced where static instability occurs 
    5658      !!---------------------------------------------------------------------- 
    57       INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step 
     59      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     60      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time-step indexocean time step 
     61      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    5862      ! 
    5963      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    60       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd 
     64!!OMP issue to be solved with outputs... 
     65!      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd 
     66!!OMP end 
    6167      !!---------------------------------------------------------------------- 
    6268      ! 
    6369      IF( nn_timing == 1 )  CALL timing_start('zdf_evd') 
    6470      ! 
     71!!OMP  this is to be suppressed 
    6572      IF( kt == nit000 ) THEN 
    6673         IF(lwp) WRITE(numout,*) 
     
    6976         IF(lwp) WRITE(numout,*) 
    7077      ENDIF 
     78!!OMP end 
    7179      ! 
    7280      ! 
    73       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
     81!!OMP issue with outputs 
     82!      zavt_evd(:,:,:) = p_avt(:,:,:)           ! set avt prior to evd application 
     83!!OMP end 
    7484      ! 
    7585      SELECT CASE ( nn_evdm ) 
     
    7787      CASE ( 1 )           !==  enhance tracer & momentum Kz  ==!   (if rn2<-1.e-12) 
    7888         ! 
    79          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application 
     89!!OMP issue with outputs 
     90!         zavm_evd(:,:,:) = p_avm(:,:,:)           ! set avm prior to evd application 
     91!!OMP end 
    8092         ! 
    81 !! change last digits results 
    82 !         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 ) THEN 
    83 !            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
    84 !            avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
     93!! change last digits results ???? 
     94!         WHERE( MAX( rn2(WRK_2D,2:jpkm1), rn2b(WRK_2D,2:jpkm1) )  <= -1.e-12 ) THEN 
     95!            p_avt(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1) 
     96!            p_avm(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1) 
    8597!         END WHERE 
    8698         ! 
    8799         DO jk = 1, jpkm1  
    88             DO jj = 2, jpjm1 
    89                DO ji = 2, jpim1 
     100            DO jj = k_Jstr, k_Jend 
     101               DO ji = k_Jstr, k_Iend 
    90102                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 
    91                      avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
    92                      avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
     103                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
     104                     p_avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
    93105                  ENDIF 
    94106               END DO 
     
    96108         END DO  
    97109         ! 
    98          zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd 
    99          CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
     110!!OMP issue with outputs 
     111!         zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd 
     112!         CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
     113!!gm end 
    100114         ! 
    101115      CASE DEFAULT         !==  enhance tracer Kz  ==!   (if rn2<-1.e-12)  
    102116!! change last digits results 
    103 !         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 ) 
    104 !            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
     117!         WHERE( MAX( rn2(WRK_2D,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 ) 
     118!            avt(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1) 
    105119!         END WHERE 
    106120 
     
    109123               DO ji = 2, jpim1 
    110124                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   & 
    111                      avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
     125                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
    112126               END DO 
    113127            END DO 
     
    116130      END SELECT  
    117131      ! 
    118       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
    119       CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
    120       IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 
     132!!OMP issue with outputs 
     133!      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
     134!      CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
     135!      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 
     136!!OMP end 
    121137      ! 
    122138      IF( nn_timing == 1 )  CALL timing_stop('zdf_evd') 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7991 r8055  
    1919   USE domvvl         ! ocean space and time domain : variable volume layer 
    2020   USE zdf_oce        ! ocean vertical physics 
    21    USE zdfsh2         ! vertical physics: shear production term of TKE 
    2221   USE zdfbfr         ! bottom friction (only for rn_bfrz0) 
    2322   USE sbc_oce        ! surface boundary condition: ocean 
     
    103102 
    104103   !! * Substitutions 
    105 #  include "vectopt_loop_substitute.h90" 
     104#  include "domain_substitute.h90" 
    106105   !!---------------------------------------------------------------------- 
    107106   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    123122 
    124123 
    125    SUBROUTINE zdf_gls( kt ) 
     124   SUBROUTINE zdf_gls( ARG_2D, kt, psh2, p_avm, p_avt ) 
    126125      !!---------------------------------------------------------------------- 
    127126      !!                   ***  ROUTINE zdf_gls  *** 
     
    130129      !!              coefficients using the GLS turbulent closure scheme. 
    131130      !!---------------------------------------------------------------------- 
    132       INTEGER, INTENT(in) ::   kt ! ocean time step 
     131      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     132      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     133      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   psh2           ! shear production term 
     134      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
     135      ! 
    133136      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments 
    134137      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars 
     
    137140      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      - 
    138141      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      - 
    139       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep 
    140       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar 
    141       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves  
    142       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves) 
    143       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before 
    144       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hmxl_b      ! mixing length at time before 
    145       REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear 
    146       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate 
    147       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi) 
    148       REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now 
    149       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal 
    150       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal 
    151       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal 
    152       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zstt, zstm  ! stability function on tracer and momentum 
     142      REAL(wp), DIMENSION(WRK_2D) ::   zdep 
     143      REAL(wp), DIMENSION(WRK_2D) ::   zkar 
     144      REAL(wp), DIMENSION(WRK_2D) ::   zflxs       ! Turbulence fluxed induced by internal waves  
     145      REAL(wp), DIMENSION(WRK_2D) ::   zhsro       ! Surface roughness (surface waves) 
     146      REAL(wp), DIMENSION(WRK_3D) ::   eb          ! tke at time before 
     147      REAL(wp), DIMENSION(WRK_3D) ::   hmxl_b      ! mixing length at time before 
     148      REAL(wp), DIMENSION(WRK_3D) ::   eps         ! dissipation rate 
     149      REAL(wp), DIMENSION(WRK_3D) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi) 
     150      REAL(wp), DIMENSION(WRK_3D) ::   psi         ! psi at time now 
     151      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_a    ! element of the first  matrix diagonal 
     152      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_b    ! element of the second matrix diagonal 
     153      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_c    ! element of the third  matrix diagonal 
     154      REAL(wp), DIMENSION(WRK_3D) ::   zstt, zstm  ! stability function on tracer and momentum 
    153155      !!-------------------------------------------------------------------- 
    154156      ! 
    155157      IF( nn_timing == 1 )   CALL timing_start('zdf_gls') 
    156158      ! 
    157       CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro ) 
    158       CALL wrk_alloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    159       CALL wrk_alloc( jpi,jpj,jpk,   zstt, zstm ) 
    160  
    161159      ! Preliminary computing 
    162160 
    163161      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    164162 
    165       ! restore before values computed GLS alone 
    166       avt(:,:,:) = avt_k(:,:,:) 
    167       avm(:,:,:) = avm_k(:,:,:) 
    168163 
    169164      ! Compute surface and bottom friction at T-points 
    170       DO jj = 2, jpjm1           
    171          DO ji = fs_2, fs_jpim1   ! vector opt.          
     165      DO jj = k_Jstr, k_Jend  
     166         DO ji = k_Istr, k_Iend 
    172167            ! 
    173168            ! surface friction 
     
    186181      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==! 
    187182      CASE ( 0 )                          ! Constant roughness           
    188          zhsro(:,:) = rn_hsro 
     183         zhsro(WRK_2D) = rn_hsro 
    189184      CASE ( 1 )             ! Standard Charnock formula 
    190          zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 
     185         zhsro(WRK_2D) = MAX(rsbc_zs1 * ustars2(WRK_2D), rn_hsro) 
    191186      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 
    192187!!gm         zcof = 2._wp * 0.6_wp / 28._wp 
    193 !!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(:,:),rsmall) )  )      ! Wave age (eq. 10) 
    194          zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    195          zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
     188!!gm         zdep(WRK_2D)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(WRK_2D),rsmall) )  )      ! Wave age (eq. 10) 
     189         zdep(WRK_2D)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(WRK_2D),rsmall))))             ! Wave age (eq. 10) 
     190         zhsro(WRK_2D) = MAX(rsbc_zs2 * ustars2(WRK_2D) * zdep(WRK_2D)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    196191      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
    197192!!gm  BUG missing a multiplicative coefficient.... 
    198          zhsro(:,:) = hsw(:,:) 
     193         zhsro(WRK_2D) = hsw(WRK_2D) 
    199194      END SELECT 
    200  
    201       !                             !==  Compute shear and dissipation rate  ==! 
    202       CALL zdf_sh2( shear ) 
    203  
    204  
     195      ! 
    205196      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==! 
    206          DO jj = 1, jpjm1 
    207             DO ji = 1, jpim1 
     197         DO jj = k_Jstr, k_Jend  
     198            DO ji = k_Istr, k_Iend 
    208199               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
    209200            END DO 
     
    212203 
    213204      ! Save tke at before time step 
    214       eb    (:,:,:) = en    (:,:,:) 
    215       hmxl_b(:,:,:) = hmxl_n(:,:,:) 
     205      eb    (WRK_2D,:) = en    (WRK_2D,:) 
     206      hmxl_b(WRK_2D,:) = hmxl_n(WRK_2D,:) 
    216207 
    217208      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    218209         DO jk = 2, jpkm1 
    219             DO jj = 2, jpjm1  
    220                DO ji = fs_2, fs_jpim1   ! vector opt. 
     210            DO jj = k_Jstr, k_Jend  
     211               DO ji = k_Istr, k_Iend 
    221212                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    222213                  zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 
     
    242233 
    243234      DO jk = 2, jpkm1 
    244          DO jj = 2, jpjm1 
    245             DO ji = 2, jpim1 
    246                ! 
    247                buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction 
     235         DO jj = k_Jstr, k_Jend  
     236            DO ji = k_Istr, k_Iend 
     237               ! 
     238               buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction 
    248239               ! 
    249240               diss = eps(ji,jj,jk)                         ! dissipation 
    250241               ! 
    251                dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    252                ! 
    253                zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term 
     242               dir = 0.5_wp + SIGN( 0.5_wp, psh2(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     243               ! 
     244               zesh2 = dir*(psh2(ji,jj,jk)+buoy)+(1._wp-dir)*psh2(ji,jj,jk)          ! production term 
    254245               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
    255246               ! 
     
    269260               zcof = rfact_tke * tmask(ji,jj,jk) 
    270261               !                                               ! lower diagonal 
    271                z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
     262               z_elem_a(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) ) 
    272263               !                                               ! upper diagonal 
    273                z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
     264               z_elem_c(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) ) 
    274265               !                                               ! diagonal 
    275266               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
     
    280271      END DO 
    281272      ! 
    282       z_elem_b(:,:,jpk) = 1._wp 
     273      z_elem_b(WRK_2D,jpk) = 1._wp 
    283274      ! 
    284275      ! Set surface condition on zwall_psi (1 at the bottom) 
    285       zwall_psi(:,:, 1 ) = zwall_psi(:,:,2) 
    286       zwall_psi(:,:,jpk) = 1._wp 
     276      zwall_psi(WRK_2D, 1 ) = zwall_psi(WRK_2D,2) 
     277      zwall_psi(WRK_2D,jpk) = 1._wp 
    287278      ! 
    288279      ! Surface boundary condition on tke 
     
    293284      CASE ( 0 )             ! Dirichlet case 
    294285      ! First level 
    295       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    296       en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    297       z_elem_a(:,:,1) = en(:,:,1) 
    298       z_elem_c(:,:,1) = 0._wp 
    299       z_elem_b(:,:,1) = 1._wp 
     286      en(WRK_2D,1) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     287      en(WRK_2D,1) = MAX(en(WRK_2D,1), rn_emin)  
     288      z_elem_a(WRK_2D,1) = en(WRK_2D,1) 
     289      z_elem_c(WRK_2D,1) = 0._wp 
     290      z_elem_b(WRK_2D,1) = 1._wp 
    300291      !  
    301292      ! One level below 
    302       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 
    303          &               / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    304       en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    305       z_elem_a(:,:,2) = 0._wp  
    306       z_elem_c(:,:,2) = 0._wp 
    307       z_elem_b(:,:,2) = 1._wp 
     293      en(WRK_2D,2) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1 * ((zhsro(WRK_2D)+gdepw_n(WRK_2D,2)) & 
     294         &               / zhsro(WRK_2D) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     295      en(WRK_2D,2) = MAX(en(WRK_2D,2), rn_emin ) 
     296      z_elem_a(WRK_2D,2) = 0._wp  
     297      z_elem_c(WRK_2D,2) = 0._wp 
     298      z_elem_b(WRK_2D,2) = 1._wp 
    308299      ! 
    309300      ! 
     
    311302      ! 
    312303      ! Dirichlet conditions at k=1 
    313       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    314       en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    315       z_elem_a(:,:,1) = en(:,:,1) 
    316       z_elem_c(:,:,1) = 0._wp 
    317       z_elem_b(:,:,1) = 1._wp 
     304      en(WRK_2D,1)       = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     305      en(WRK_2D,1)       = MAX(en(WRK_2D,1), rn_emin)       
     306      z_elem_a(WRK_2D,1) = en(WRK_2D,1) 
     307      z_elem_c(WRK_2D,1) = 0._wp 
     308      z_elem_b(WRK_2D,1) = 1._wp 
    318309      ! 
    319310      ! at k=2, set de/dz=Fw 
    320311      !cbr 
    321       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    322       z_elem_a(:,:,2) = 0._wp 
    323       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
    324       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    325           &                       * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    326  
    327       en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2) 
     312      z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) +  z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b 
     313      z_elem_a(WRK_2D,2) = 0._wp 
     314      zkar(WRK_2D)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(WRK_2D,1)/zhsro(WRK_2D)) )) 
     315      zflxs(WRK_2D)      = rsbc_tke2 * ustars2(WRK_2D)**1.5_wp * zkar(WRK_2D) & 
     316          &                       * ((zhsro(WRK_2D)+gdept_n(WRK_2D,1)) / zhsro(WRK_2D) )**(1.5_wp*ra_sf) 
     317 
     318      en(WRK_2D,2) = en(WRK_2D,2) + zflxs(WRK_2D)/e3w_n(WRK_2D,2) 
    328319      ! 
    329320      ! 
     
    338329         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 
    339330         !                      ! Balance between the production and the dissipation terms 
    340          DO jj = 2, jpjm1 
    341             DO ji = fs_2, fs_jpim1   ! vector opt. 
     331         DO jj = k_Jstr, k_Jend  
     332            DO ji = k_Istr, k_Iend 
    342333               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    343334               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     
    359350      CASE ( 1 )             ! Neumman boundary condition 
    360351         !                       
    361          DO jj = 2, jpjm1 
    362             DO ji = fs_2, fs_jpim1   ! vector opt. 
     352         DO jj = k_Jstr, k_Jend  
     353            DO ji = k_Istr, k_Iend 
    363354               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    364355               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     
    382373      ! 
    383374      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    384          DO jj = 2, jpjm1 
    385             DO ji = fs_2, fs_jpim1    ! vector opt. 
     375         DO jj = k_Jstr, k_Jend  
     376            DO ji = k_Istr, k_Iend 
    386377               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1) 
    387378            END DO 
     
    389380      END DO 
    390381      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    391          DO jj = 2, jpjm1 
    392             DO ji = fs_2, fs_jpim1    ! vector opt. 
     382         DO jj = k_Jstr, k_Jend  
     383            DO ji = k_Istr, k_Iend 
    393384               z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1) 
    394385            END DO 
     
    396387      END DO 
    397388      DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    398          DO jj = 2, jpjm1 
    399             DO ji = fs_2, fs_jpim1    ! vector opt. 
     389         DO jj = k_Jstr, k_Jend  
     390            DO ji = k_Istr, k_Iend 
    400391               en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk) 
    401392            END DO 
     
    403394      END DO 
    404395      !                                            ! set the minimum value of tke  
    405       en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
     396      en(WRK_3D) = MAX( en(WRK_3D), rn_emin ) 
    406397 
    407398      !!----------------------------------------!! 
     
    415406      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    416407         DO jk = 2, jpkm1 
    417             DO jj = 2, jpjm1 
    418                DO ji = fs_2, fs_jpim1   ! vector opt. 
     408            DO jj = k_Jstr, k_Jend  
     409               DO ji = k_Istr, k_Iend 
    419410                  psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 
    420411               END DO 
     
    424415      CASE( 1 )               ! k-eps 
    425416         DO jk = 2, jpkm1 
    426             DO jj = 2, jpjm1 
    427                DO ji = fs_2, fs_jpim1   ! vector opt. 
     417            DO jj = k_Jstr, k_Jend  
     418               DO ji = k_Istr, k_Iend 
    428419                  psi(ji,jj,jk)  = eps(ji,jj,jk) 
    429420               END DO 
     
    433424      CASE( 2 )               ! k-w 
    434425         DO jk = 2, jpkm1 
    435             DO jj = 2, jpjm1 
    436                DO ji = fs_2, fs_jpim1   ! vector opt. 
     426            DO jj = k_Jstr, k_Jend  
     427               DO ji = k_Istr, k_Iend 
    437428                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 
    438429               END DO 
     
    442433      CASE( 3 )               ! generic 
    443434         DO jk = 2, jpkm1 
    444             DO jj = 2, jpjm1 
    445                DO ji = fs_2, fs_jpim1   ! vector opt. 
     435            DO jj = k_Jstr, k_Jend  
     436               DO ji = k_Istr, k_Iend 
    446437                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn  
    447438               END DO 
     
    459450 
    460451      DO jk = 2, jpkm1 
    461          DO jj = 2, jpjm1 
    462             DO ji = fs_2, fs_jpim1   ! vector opt. 
     452         DO jj = k_Jstr, k_Jend  
     453            DO ji = k_Istr, k_Iend 
    463454               ! 
    464455               ! psi / k 
     
    471462               ! 
    472463               ! shear prod. - stratif. destruction 
    473                prod = rpsi1 * zratio * shear(ji,jj,jk) 
     464               prod = rpsi1 * zratio * psh2(ji,jj,jk) 
    474465               ! 
    475466               ! stratif. destruction 
    476                buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
     467               buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
    477468               ! 
    478469               ! shear prod. - stratif. destruction 
     
    487478               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
    488479               !                                               ! lower diagonal 
    489                z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 
     480               z_elem_a(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) ) 
    490481               !                                               ! upper diagonal 
    491                z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) ) 
     482               z_elem_c(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) ) 
    492483               !                                               ! diagonal 
    493484               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
     
    498489      END DO 
    499490      ! 
    500       z_elem_b(:,:,jpk) = 1._wp 
     491      z_elem_b(WRK_2D,jpk) = 1._wp 
    501492 
    502493      ! Surface boundary condition on psi 
     
    508499      ! 
    509500      ! Surface value 
    510       zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    511       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    512       z_elem_a(:,:,1) = psi(:,:,1) 
    513       z_elem_c(:,:,1) = 0._wp 
    514       z_elem_b(:,:,1) = 1._wp 
     501      zdep(WRK_2D)       = zhsro(WRK_2D) * rl_sf ! Cosmetic 
     502      psi (WRK_2D,1)     = rc0**rpp * en(WRK_2D,1)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1) 
     503      z_elem_a(WRK_2D,1) = psi(WRK_2D,1) 
     504      z_elem_c(WRK_2D,1) = 0._wp 
     505      z_elem_b(WRK_2D,1) = 1._wp 
    515506      ! 
    516507      ! One level below 
    517       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
    518       zdep(:,:)       = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 
    519       psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    520       z_elem_a(:,:,2) = 0._wp 
    521       z_elem_c(:,:,2) = 0._wp 
    522       z_elem_b(:,:,2) = 1._wp 
     508      zkar(WRK_2D)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(WRK_2D,2)/zhsro(WRK_2D) ))) 
     509      zdep(WRK_2D)       = (zhsro(WRK_2D) + gdepw_n(WRK_2D,2)) * zkar(WRK_2D) 
     510      psi (WRK_2D,2)     = rc0**rpp * en(WRK_2D,2)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1) 
     511      z_elem_a(WRK_2D,2) = 0._wp 
     512      z_elem_c(WRK_2D,2) = 0._wp 
     513      z_elem_b(WRK_2D,2) = 1._wp 
    523514      !  
    524515      ! 
     
    526517      ! 
    527518      ! Surface value: Dirichlet 
    528       zdep(:,:)       = zhsro(:,:) * rl_sf 
    529       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    530       z_elem_a(:,:,1) = psi(:,:,1) 
    531       z_elem_c(:,:,1) = 0._wp 
    532       z_elem_b(:,:,1) = 1._wp 
     519      zdep(WRK_2D)       = zhsro(WRK_2D) * rl_sf 
     520      psi (WRK_2D,1)     = rc0**rpp * en(WRK_2D,1)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1) 
     521      z_elem_a(WRK_2D,1) = psi(WRK_2D,1) 
     522      z_elem_c(WRK_2D,1) = 0._wp 
     523      z_elem_b(WRK_2D,1) = 1._wp 
    533524      ! 
    534525      ! Neumann condition at k=2 
    535       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    536       z_elem_a(:,:,2) = 0._wp 
     526      z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) +  z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b 
     527      z_elem_a(WRK_2D,2) = 0._wp 
    537528      ! 
    538529      ! Set psi vertical flux at the surface: 
    539       zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    540       zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    541       zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    542       zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
    543              & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 
    544       zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    545       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     530      zkar(WRK_2D) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(WRK_2D,1)/zhsro(WRK_2D) )) ! Lengh scale slope 
     531      zdep(WRK_2D) = ((zhsro(WRK_2D) + gdept_n(WRK_2D,1)) / zhsro(WRK_2D))**(rmm*ra_sf) 
     532      zflxs(WRK_2D) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(WRK_2D))*(1._wp + rsbc_tke1*zdep(WRK_2D))**(2._wp*rmm/3._wp-1_wp) 
     533      zdep(WRK_2D) =  rsbc_psi1 * (zwall_psi(WRK_2D,1)*p_avm(WRK_2D,1)+zwall_psi(WRK_2D,2)*p_avm(WRK_2D,2)) * & 
     534             & ustars2(WRK_2D)**rmm * zkar(WRK_2D)**rnn * (zhsro(WRK_2D) + gdept_n(WRK_2D,1))**(rnn-1.) 
     535      zflxs(WRK_2D) = zdep(WRK_2D) * zflxs(WRK_2D) 
     536      psi(WRK_2D,2) = psi(WRK_2D,2) + zflxs(WRK_2D) / e3w_n(WRK_2D,2) 
    546537      !    
    547538      ! 
     
    557548         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0 
    558549         !                      ! Balance between the production and the dissipation terms 
    559          DO jj = 2, jpjm1 
    560             DO ji = fs_2, fs_jpim1   ! vector opt. 
     550         DO jj = k_Jstr, k_Jend  
     551            DO ji = k_Istr, k_Iend 
    561552               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    562553               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     
    578569      CASE ( 1 )             ! Neumman boundary condition 
    579570         !                       
    580          DO jj = 2, jpjm1 
    581             DO ji = fs_2, fs_jpim1   ! vector opt. 
     571         DO jj = k_Jstr, k_Jend  
     572            DO ji = k_Istr, k_Iend 
    582573               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    583574               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
     
    597588               ! Set psi vertical flux at the bottom: 
    598589               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1) 
    599                zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   & 
     590               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   & 
    600591                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
    601592               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 
     
    609600      ! 
    610601      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    611          DO jj = 2, jpjm1 
    612             DO ji = fs_2, fs_jpim1    ! vector opt. 
     602         DO jj = k_Jstr, k_Jend  
     603            DO ji = k_Istr, k_Iend 
    613604               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1) 
    614605            END DO 
     
    616607      END DO 
    617608      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    618          DO jj = 2, jpjm1 
    619             DO ji = fs_2, fs_jpim1    ! vector opt. 
     609         DO jj = k_Jstr, k_Jend  
     610            DO ji = k_Istr, k_Iend 
    620611               z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1) 
    621612            END DO 
     
    623614      END DO 
    624615      DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    625          DO jj = 2, jpjm1 
    626             DO ji = fs_2, fs_jpim1    ! vector opt. 
     616         DO jj = k_Jstr, k_Jend  
     617            DO ji = k_Istr, k_Iend 
    627618               psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk) 
    628619            END DO 
     
    637628      CASE( 0 )               ! k-kl  (Mellor-Yamada) 
    638629         DO jk = 1, jpkm1 
    639             DO jj = 2, jpjm1 
    640                DO ji = fs_2, fs_jpim1   ! vector opt. 
     630            DO jj = k_Jstr, k_Jend  
     631               DO ji = k_Istr, k_Iend 
    641632                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 
    642633               END DO 
     
    646637      CASE( 1 )               ! k-eps 
    647638         DO jk = 1, jpkm1 
    648             DO jj = 2, jpjm1 
    649                DO ji = fs_2, fs_jpim1   ! vector opt. 
     639            DO jj = k_Jstr, k_Jend  
     640               DO ji = k_Istr, k_Iend 
    650641                  eps(ji,jj,jk) = psi(ji,jj,jk) 
    651642               END DO 
     
    655646      CASE( 2 )               ! k-w 
    656647         DO jk = 1, jpkm1 
    657             DO jj = 2, jpjm1 
    658                DO ji = fs_2, fs_jpim1   ! vector opt. 
     648            DO jj = k_Jstr, k_Jend  
     649               DO ji = k_Istr, k_Iend 
    659650                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
    660651               END DO 
     
    667658         zex2  = -1._wp / rnn 
    668659         DO jk = 1, jpkm1 
    669             DO jj = 2, jpjm1 
    670                DO ji = fs_2, fs_jpim1   ! vector opt. 
     660            DO jj = k_Jstr, k_Jend  
     661               DO ji = k_Istr, k_Iend 
    671662                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 
    672663               END DO 
     
    679670      ! -------------------------------------------------- 
    680671      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 
    681          DO jj = 2, jpjm1 
    682             DO ji = fs_2, fs_jpim1    ! vector opt. 
     672         DO jj = k_Jstr, k_Jend  
     673            DO ji = k_Istr, k_Iend 
    683674               ! limitation 
    684675               eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
     
    699690      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions 
    700691         DO jk = 2, jpkm1 
    701             DO jj = 2, jpjm1 
    702                DO ji = fs_2, fs_jpim1   ! vector opt. 
     692            DO jj = k_Jstr, k_Jend  
     693               DO ji = k_Istr, k_Iend 
    703694                  ! zcof =  l²/q² 
    704695                  zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
     
    720711      CASE ( 2, 3 )               ! Canuto stability functions 
    721712         DO jk = 2, jpkm1 
    722             DO jj = 2, jpjm1 
    723                DO ji = fs_2, fs_jpim1   ! vector opt. 
     713            DO jj = k_Jstr, k_Jend  
     714               DO ji = k_Istr, k_Iend 
    724715                  ! zcof =  l²/q² 
    725716                  zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
     
    730721                  gh = gh * rf6 
    731722                  ! Gm =  M²l²/q² Shear number 
    732                   shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall ) 
     723                  shr = psh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 
    733724                  gm = MAX( shr * zcof , 1.e-10 ) 
    734725                  gm = gm * rf6 
     
    751742      ! Lines below are useless if GOTM style Dirichlet conditions are used 
    752743 
    753       zstm(:,:,1) = zstm(:,:,2) 
     744      zstm(WRK_2D,1) = zstm(WRK_2D,2) 
    754745 
    755746!!gm should be done for ISF (top boundary cond.) 
    756       DO jj = 2, jpjm1 
    757          DO ji = fs_2, fs_jpim1   ! vector opt. 
     747      DO jj = k_Jstr, k_Jend  
     748         DO ji = k_Istr, k_Iend 
    758749            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 
    759750         END DO 
     
    763754      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used 
    764755      DO jk = 1, jpk 
    765          DO jj = 2, jpjm1 
    766             DO ji = fs_2, fs_jpim1   ! vector opt. 
     756         DO jj = k_Jstr, k_Jend  
     757            DO ji = k_Istr, k_Iend 
    767758               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
    768759               zav           = zsqen * zstt(ji,jj,jk) 
    769                avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
     760               p_avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
    770761               zav           = zsqen * zstm(ji,jj,jk) 
    771                avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
     762               p_avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
    772763            END DO 
    773764         END DO 
    774765      END DO 
    775       avt(:,:,1)  = 0._wp 
    776 !!gm I'm sure this is needed to compute the shear term at next time-step 
    777       CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. ) 
     766      p_avt(WRK_2D,1)  = 0._wp 
    778767      ! 
    779768      IF(ln_ctl) THEN 
     
    781770         CALL prt_ctl( tab3d_1=avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk ) 
    782771      ENDIF 
    783       ! 
    784       avt_k(:,:,:) = avt(:,:,:)      !==  store avt, avm values computed by GLS only  ==! 
    785       avm_k(:,:,:) = avm(:,:,:) 
    786       ! 
    787       IF( lrst_oce )   CALL gls_rst( kt, 'WRITE' )     ! write the GLS info in the restart file 
    788       ! 
    789       CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro ) 
    790       CALL wrk_dealloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    791       CALL wrk_dealloc( jpi,jpj,jpk,   zstt, zstm ) 
    792772      ! 
    793773      IF( nn_timing == 1 )   CALL timing_stop('zdf_gls') 
     
    837817      IF(lwp) THEN                     !* Control print 
    838818         WRITE(numout,*) 
    839          WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme' 
     819         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme' 
    840820         WRITE(numout,*) '~~~~~~~~~~~~' 
    841821         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters' 
     
    855835         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
    856836         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0 
     837         WRITE(numout,*) 
    857838      ENDIF 
    858839 
     
    861842 
    862843      !                                !* Check of some namelist values 
    863       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
    864       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
    865       IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 
    866       IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
    867       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    868       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
     844      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
     845      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 
     846      IF( nn_z0_met  < 0   .OR. nn_z0_met    > 3 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 
     847      IF( nn_z0_met == 3  .AND. .NOT.ln_wave     )  CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
     848      IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
     849      IF( nn_clos      < 0 .OR. nn_clos      > 3 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    869850 
    870851      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
     
    872853      CASE( 0 )                              ! k-kl  (Mellor-Yamada) 
    873854         ! 
    874          IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 
     855         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)' 
     856         IF(lwp) WRITE(numout,*) 
    875857         rpp     = 0._wp 
    876858         rmm     = 1._wp 
     
    890872      CASE( 1 )                              ! k-eps 
    891873         ! 
    892          IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 
     874         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen' 
     875         IF(lwp) WRITE(numout,*) 
    893876         rpp     =  3._wp 
    894877         rmm     =  1.5_wp 
     
    908891      CASE( 2 )                              ! k-omega 
    909892         ! 
    910          IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 
     893         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen' 
     894         IF(lwp) WRITE(numout,*) 
    911895         rpp     = -1._wp 
    912896         rmm     =  0.5_wp 
     
    926910      CASE( 3 )                              ! generic 
    927911         ! 
    928          IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 
     912         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen' 
     913         IF(lwp) WRITE(numout,*) 
    929914         rpp     = 2._wp 
    930915         rmm     = 1._wp 
     
    949934      CASE ( 0 )                             ! Galperin stability functions 
    950935         ! 
    951          IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 
     936         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin' 
    952937         rc2     =  0._wp 
    953938         rc3     =  0._wp 
     
    961946      CASE ( 1 )                             ! Kantha-Clayson stability functions 
    962947         ! 
    963          IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 
     948         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson' 
    964949         rc2     =  0.7_wp 
    965950         rc3     =  0.2_wp 
     
    973958      CASE ( 2 )                             ! Canuto A stability functions 
    974959         ! 
    975          IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 
     960         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A' 
    976961         rs0 = 1.5_wp * rl1 * rl5*rl5 
    977962         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 
     
    997982      CASE ( 3 )                             ! Canuto B stability functions 
    998983         ! 
    999          IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 
     984         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B' 
    1000985         rs0 = 1.5_wp * rm1 * rm5*rm5 
    1001986         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 
     
    10521037      IF(lwp) THEN                     !* Control print 
    10531038         WRITE(numout,*) 
    1054          WRITE(numout,*) 'Limit values' 
    1055          WRITE(numout,*) '~~~~~~~~~~~~' 
    1056          WRITE(numout,*) 'Parameter  m = ',rmm 
    1057          WRITE(numout,*) 'Parameter  n = ',rnn 
    1058          WRITE(numout,*) 'Parameter  p = ',rpp 
    1059          WRITE(numout,*) 'rpsi1   = ',rpsi1 
    1060          WRITE(numout,*) 'rpsi2   = ',rpsi2 
    1061          WRITE(numout,*) 'rpsi3m  = ',rpsi3m 
    1062          WRITE(numout,*) 'rpsi3p  = ',rpsi3p 
    1063          WRITE(numout,*) 'rsc_tke = ',rsc_tke 
    1064          WRITE(numout,*) 'rsc_psi = ',rsc_psi 
    1065          WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 
    1066          WRITE(numout,*) 'rc0     = ',rc0 
     1039         WRITE(numout,*) '   Limit values :' 
     1040         WRITE(numout,*) '      Parameter  m = ',rmm 
     1041         WRITE(numout,*) '      Parameter  n = ',rnn 
     1042         WRITE(numout,*) '      Parameter  p = ',rpp 
     1043         WRITE(numout,*) '      rpsi1   = ',rpsi1 
     1044         WRITE(numout,*) '      rpsi2   = ',rpsi2 
     1045         WRITE(numout,*) '      rpsi3m  = ',rpsi3m 
     1046         WRITE(numout,*) '      rpsi3p  = ',rpsi3p 
     1047         WRITE(numout,*) '      rsc_tke = ',rsc_tke 
     1048         WRITE(numout,*) '      rsc_psi = ',rsc_psi 
     1049         WRITE(numout,*) '      rsc_psi0 = ',rsc_psi0 
     1050         WRITE(numout,*) '      rc0     = ',rc0 
    10671051         WRITE(numout,*) 
    1068          WRITE(numout,*) 'Shear free turbulence parameters:' 
    1069          WRITE(numout,*) 'rcm_sf  = ',rcm_sf 
    1070          WRITE(numout,*) 'ra_sf   = ',ra_sf 
    1071          WRITE(numout,*) 'rl_sf   = ',rl_sf 
    1072          WRITE(numout,*) 
     1052         WRITE(numout,*) '   Shear free turbulence parameters:' 
     1053         WRITE(numout,*) '      rcm_sf  = ',rcm_sf 
     1054         WRITE(numout,*) '      ra_sf   = ',ra_sf 
     1055         WRITE(numout,*) '      rl_sf   = ',rl_sf 
    10731056      ENDIF 
    10741057 
     
    10901073      ! 
    10911074      !                                !* Wall proximity function 
    1092       zwall (:,:,:) = 1._wp * tmask(:,:,:) 
     1075      zwall(:,:,:) = 1._wp * tmask(:,:,:) 
    10931076 
    10941077      !                                !* read or initialize all required files   
     
    11331116               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 
    11341117            ELSE                         
    1135                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 
     1118               IF(lwp) WRITE(numout,*) 
     1119               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values' 
    11361120               en    (:,:,:) = rn_emin 
    11371121               hmxl_n(:,:,:) = 0.05_wp 
    1138                avt_k (:,:,:) = avt(:,:,:) 
    1139                avm_k (:,:,:) = avm(:,:,:) 
    1140                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO 
     1122               ! avt_k, avm_k already set to the background value in zdf_phy_init 
    11411123            ENDIF 
    11421124         ELSE                                   !* Start from rest 
    1143             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 
     1125            IF(lwp) WRITE(numout,*) 
     1126            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values' 
    11441127            en    (:,:,:) = rn_emin 
    11451128            hmxl_n(:,:,:) = 0.05_wp 
    1146             DO jk = 1, jpk 
    1147                avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 
    1148                avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk) 
    1149             END DO 
     1129            ! avt_k, avm_k already set to the background value in zdf_phy_init 
    11501130         ENDIF 
    11511131         ! 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

    r7990 r8055  
    3535   PUBLIC   zdf_iwm         ! called in step module  
    3636   PUBLIC   zdf_iwm_init    ! called in nemogcm module  
    37    PUBLIC   zdf_iwm_alloc   ! called in nemogcm module 
    3837 
    3938   !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing * 
     
    5655 
    5756   !! * Substitutions 
    58 #  include "vectopt_loop_substitute.h90" 
     57#  include "domain_substitute.h90"    
    5958   !!---------------------------------------------------------------------- 
    6059   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     
    7877 
    7978 
    80    SUBROUTINE zdf_iwm( kt ) 
     79   SUBROUTINE zdf_iwm( ARG_2D, kt, p_avm, p_avt, p_avs ) 
    8180      !!---------------------------------------------------------------------- 
    8281      !!                  ***  ROUTINE zdf_iwm  *** 
     
    127126      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
    128127      !!---------------------------------------------------------------------- 
    129       INTEGER, INTENT(in) ::   kt   ! ocean time-step  
     128      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     129      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     130      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
     131      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    130132      ! 
    131133      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    132       REAL(wp) ::   ztpc         ! scalar workspace 
    133       REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure 
    134       REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth 
    135       REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom 
    136       REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution 
    137       REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid) 
    138       REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid) 
    139       REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter 
     134      REAL(wp) ::   zztemp       ! scalar workspace 
     135      REAL(wp), DIMENSION(WRK_2D) ::  zfact     ! Used for vertical structure 
     136      REAL(wp), DIMENSION(WRK_2D) ::  zhdep     ! Ocean depth 
     137      REAL(wp), DIMENSION(WRK_3D) ::  zwkb      ! WKB-stretched height above bottom 
     138      REAL(wp), DIMENSION(WRK_3D) ::  zweight   ! Weight for high mode vertical distribution 
     139      REAL(wp), DIMENSION(WRK_3D) ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     140      REAL(wp), DIMENSION(WRK_3D) ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     141      REAL(wp), DIMENSION(WRK_3D) ::  zReb      ! Turbulence intensity parameter 
    140142      !!---------------------------------------------------------------------- 
    141143      ! 
    142144      IF( nn_timing == 1 )   CALL timing_start('zdf_iwm') 
    143145      ! 
    144       CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
    145       CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    146  
    147146      !                          ! ----------------------------- ! 
    148147      !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
     
    151150      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    152151      !                                                 using an exponential decay from the seafloor. 
    153       DO jj = 1, jpj                ! part independent of the level 
    154          DO ji = 1, jpi 
     152      DO jj = k_Jstr, k_Jend 
     153         DO ji = k_Istr, k_Iend     ! part independent of the level 
    155154            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    156155            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     
    158157         END DO 
    159158      END DO 
    160  
    161159      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    162          emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      & 
    163             &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   & 
     160         DO jj = k_Jstr, k_Jend 
     161            DO ji = k_Istr, k_Iend 
     162               emix_iwm(ji,jj,jk) = zfact(ji,jj) * wmask(ji,jj,jk)                              &  
     163                  &   * (  EXP( ( gde3w_n(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )      & 
     164                  &      - EXP( ( gde3w_n(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )  )   & 
     165                  &   / ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) 
     166 
    164167!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point 
    165             &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     168!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all 
     169 
     170            END DO 
     171         END DO 
    166172      END DO 
    167173 
    168174      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
    169175      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
    170       ! 
     176      !                                          ! (NB: N2 is masked, so no use of wmask here) 
    171177      SELECT CASE ( nn_zpyc ) 
    172178      ! 
    173179      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    174180         ! 
    175          zfact(:,:) = 0._wp 
    176          DO jk = 2, jpkm1              ! part independent of the level 
    177             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    178          END DO 
    179          ! 
    180          DO jj = 1, jpj 
    181             DO ji = 1, jpi 
    182                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    183             END DO 
    184          END DO 
     181         zfact(WRK_2D) = 0._wp         ! part independent of the level 
     182         DO jk = 2, jpkm1 
     183            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     184         END DO 
     185         ! 
     186         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     187!         DO jj = k_Jstr, k_Jend 
     188!            DO ji = k_Istr, k_Iend 
     189!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     190!            END DO 
     191!         END DO 
     192         !                             ! complete with the level-dependent part 
     193         DO jk = 2, jpkm1 
     194            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     195         END DO 
     196         ! 
     197      CASE ( 2 )               ! Dissipation scales as N^2 
     198         ! 
     199         zfact(WRK_2D) = 0._wp 
     200         DO jk = 2, jpkm1              ! part independent of the level  
     201            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * MAX( 0._wp, rn2(WRK_2D,jk) ) 
     202         END DO 
     203         ! 
     204         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     205!         DO jj = k_Jstr, k_Jend 
     206!            DO ji = k_Istr, k_Iend 
     207!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     208!            END DO 
     209!         END DO 
    185210         ! 
    186211         DO jk = 2, jpkm1              ! complete with the level-dependent part 
    187             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    188          END DO 
    189          ! 
    190       CASE ( 2 )               ! Dissipation scales as N^2 
    191          ! 
    192          zfact(:,:) = 0._wp 
    193          DO jk = 2, jpkm1              ! part independent of the level 
    194             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    195          END DO 
    196          ! 
    197          DO jj= 1, jpj 
    198             DO ji = 1, jpi 
    199                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    200             END DO 
    201          END DO 
    202          ! 
    203          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    204             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     212            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * MAX( 0._wp, rn2(WRK_2D,jk) ) 
    205213         END DO 
    206214         ! 
     
    210218      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    211219      ! 
    212       zwkb(:,:,:) = 0._wp 
    213       zfact(:,:) = 0._wp 
     220      zwkb (WRK_3D) = 0._wp 
     221      zfact(WRK_2D) = 0._wp 
    214222      DO jk = 2, jpkm1 
    215          zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    216          zwkb(:,:,jk) = zfact(:,:) 
     223!!gm  this is potentially dangerous  
     224!         zfact(WRK_2D)    = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     225!         zwkb (WRK_2D,jk) = zfact(WRK_2D) 
     226         DO jj = k_Jstr, k_Jend 
     227            DO ji = k_Istr, k_Iend 
     228               zztemp = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) 
     229               zfact(ji,jj)    = zztemp 
     230               zwkb (ji,jj,jk) = zztemp 
     231            END DO 
     232         END DO          
    217233      END DO 
    218234      ! 
    219235      DO jk = 2, jpkm1 
    220          DO jj = 1, jpj 
    221             DO ji = 1, jpi 
     236         DO jj = k_Jstr, k_Jend 
     237            DO ji = k_Istr, k_Iend 
    222238               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    223                   &                                     * tmask(ji,jj,jk) / zfact(ji,jj) 
     239                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    224240            END DO 
    225241         END DO 
    226242      END DO 
    227       zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
    228       ! 
    229       zweight(:,:,:) = 0._wp 
     243      zwkb(WRK_2D,1) = zhdep(WRK_2D) * wmask(WRK_2D,1) 
     244      ! 
     245      zweight(WRK_3D) = 0._wp 
    230246      DO jk = 2, jpkm1 
    231          zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    & 
    232             &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  ) 
    233       END DO 
    234       ! 
    235       zfact(:,:) = 0._wp 
     247         zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(:,:)     & 
     248            &               * (  EXP( -zwkb(WRK_2D,jk  ) / hbot_iwm(WRK_2D) )  & 
     249            &                  - EXP( -zwkb(WRK_2D,jk-1) / hbot_iwm(WRK_2D) )  ) 
     250      END DO 
     251      ! 
     252      zfact(WRK_2D) = 0._wp 
    236253      DO jk = 2, jpkm1              ! part independent of the level 
    237          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    238       END DO 
    239       ! 
    240       DO jj = 1, jpj 
    241          DO ji = 1, jpi 
    242             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    243          END DO 
    244       END DO 
     254         zfact(WRK_2D) = zfact(WRK_2D) + zweight(WRK_2D,jk) 
     255      END DO 
     256      ! 
     257      WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = ebot_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     258!      DO jj = k_Jstr, k_Jend 
     259!         DO ji = k_Istr, k_Iend 
     260!            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     261!         END DO 
     262!      END DO 
    245263      ! 
    246264      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    247          emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    248             &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    249       END DO 
    250       ! 
     265         emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zweight(WRK_2D,jk) * zfact(WRK_2D) * wmask(WRK_2D,jk)   & 
     266            &                                / ( gde3w_n(WRK_2D,jk) - gde3w_n(WRK_2D,jk-1) ) 
     267!!gm  use of e3t_n just above? 
     268      END DO 
     269      ! 
     270!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    251271      ! Calculate molecular kinematic viscosity 
    252       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
    253          &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     272      znu_t(WRK_3D) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(WRK_3D,jp_tem)             & 
     273         &                        + 0.00694_wp * tsn(WRK_3D,jp_tem) * tsn(WRK_3D,jp_tem)   & 
     274         &                        + 0.02305_wp * tsn(WRK_3D,jp_sal)  ) * tmask(WRK_3D) * r1_rau0 
    254275      DO jk = 2, jpkm1 
    255          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    256       END DO 
     276         znu_w(WRK_2D,jk) = 0.5_wp * ( znu_t(WRK_2D,jk-1) + znu_t(WRK_2D,jk) ) * wmask(WRK_2D,jk) 
     277      END DO 
     278!!gm end 
    257279      ! 
    258280      ! Calculate turbulence intensity parameter Reb 
    259281      DO jk = 2, jpkm1 
    260          zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     282         zReb(WRK_2D,jk) = emix_iwm(WRK_2D,jk) / MAX( 1.e-20_wp, znu_w(WRK_2D,jk) * rn2(WRK_2D,jk) ) 
    261283      END DO 
    262284      ! 
    263285      ! Define internal wave-induced diffusivity 
    264286      DO jk = 2, jpkm1 
    265          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     287         zav_wave(WRK_2D,jk) = znu_w(WRK_2D,jk) * zReb(WRK_2D,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    266288      END DO 
    267289      ! 
    268290      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    269291         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
    270             DO jj = 1, jpj 
    271                DO ji = 1, jpi 
     292            DO jj = k_Jstr, k_Jend 
     293               DO ji = k_Istr, k_Iend 
    272294                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    273295                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     
    281303      ! 
    282304      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    283          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
     305         zav_wave(WRK_2D,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(WRK_2D,jk) ), 1.e-2_wp  ) * wmask(WRK_2D,jk) 
    284306      END DO 
    285307      ! 
    286308      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    287          ztpc = 0._wp 
     309         zztemp = 0._wp 
    288310!!gm used of glosum 3D.... 
    289311         DO jk = 2, jpkm1 
    290             DO jj = 1, jpj 
    291                DO ji = 1, jpi 
    292                   ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
     312            DO jj = k_Jstr, k_Jend 
     313               DO ji = k_Istr, k_Iend 
     314                  zztemp = zztemp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
    293315                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    294316               END DO 
    295317            END DO 
    296318         END DO 
    297          IF( lk_mpp )   CALL mpp_sum( ztpc ) 
    298          ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     319         IF( lk_mpp )   CALL mpp_sum( zztemp ) 
     320         zztemp = rau0 * zztemp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    299321         ! 
    300322         IF(lwp) THEN 
     
    303325            WRITE(numout,*) '~~~~~~~ ' 
    304326            WRITE(numout,*) 
    305             WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     327            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztemp * 1.e-12_wp, 'TW' 
    306328         ENDIF 
    307329      ENDIF 
     
    313335      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
    314336         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
    315             DO jj = 1, jpj 
    316                DO ji = 1, jpi 
     337            DO jj = k_Jstr, k_Jend 
     338               DO ji = k_Istr, k_Iend 
    317339                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  & 
    318340                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   & 
     
    323345         CALL iom_put( "av_ratio", zav_ratio ) 
    324346         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    325             avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    326             avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
    327             avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     347            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk) * zav_ratio(WRK_2D,jk) 
     348            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     349            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
    328350         END DO 
    329351         ! 
    330352      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    331353         DO jk = 2, jpkm1 
    332             avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) 
    333             avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
    334             avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     354            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     355            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     356            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
    335357         END DO 
    336358      ENDIF 
     
    341363                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm) 
    342364      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    343          bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    344          pcmap_iwm(:,:) = 0._wp 
    345365         DO jk = 2, jpkm1 
    346             pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk) 
    347          END DO 
    348          pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:) 
    349          CALL iom_put( "bflx_iwm", bflx_iwm ) 
     366            bflx_iwm(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * zav_wave(WRK_2D,jk) 
     367         END DO 
     368         pcmap_iwm(WRK_2D) = 0._wp 
     369         DO jk = 2, jpkm1 
     370            pcmap_iwm(WRK_2D) = pcmap_iwm(WRK_2D) + e3w_n(WRK_2D,jk) * bflx_iwm(WRK_2D,jk) 
     371         END DO 
     372         pcmap_iwm(WRK_2D) = rau0 * pcmap_iwm(WRK_2D) 
     373         CALL iom_put( "bflx_iwm" , bflx_iwm  ) 
    350374         CALL iom_put( "pcmap_iwm", pcmap_iwm ) 
    351375      ENDIF 
     
    353377      CALL iom_put( "emix_iwm", emix_iwm ) 
    354378       
    355       CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
    356       CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    357  
    358379      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
    359380      ! 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

    r7990 r8055  
    1111   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff.  
    1212   !!---------------------------------------------------------------------- 
    13    USE par_oce        ! ocean parameters 
     13   USE oce            ! ocean dynamics and tracers variables 
    1414   USE zdf_oce        ! vertical physics: shared variables          
    1515   USE zdfbfr         ! vertical physics: bottom friction 
    16    USE zdfric         ! vertical physics: Richardson vertical mixing    
     16   USE zdfsh2         ! vertical physics: shear production term of TKE 
     17   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing    
    1718   USE zdftke         ! vertical physics: TKE vertical mixing 
    1819   USE zdfgls         ! vertical physics: GLS vertical mixing 
     
    4445   INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz 
    4546   INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz 
    46        
     47 
     48   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 
     49 
     50   !! * Substitutions 
     51#  include "domain_substitute.h90"    
    4752   !!---------------------------------------------------------------------- 
    4853   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
     
    6166      !!                set horizontal shape and vertical profile of background mixing coef. 
    6267      !!---------------------------------------------------------------------- 
     68      INTEGER ::   jk            ! dummy loop indices 
    6369      INTEGER ::   ioptio, ios   ! local integers 
    6470      !! 
     
    9096         WRITE(numout,*) '      vertical closure scheme' 
    9197         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst 
    92          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfric = ', ln_zdfric 
    93          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdftke = ', ln_zdftke 
    94          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfgls = ', ln_zdfgls 
     98         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric 
     99         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke 
     100         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls 
    95101         WRITE(numout,*) '      convection: ' 
    96102         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd 
     
    137143      ENDIF 
    138144      ! 
    139       !                                   ! set 1st/last level Av to zero one for all 
    140       avt(:,:,1) = 0._wp   ;   avs(:,:,1) = 0._wp   ;   avm(:,:,1) = 0._wp 
     145      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k) 
     146         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk) 
     147         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk) 
     148      END DO 
     149!!gm  to be tested only the 1st & last levels 
     150!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp 
     151!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp 
     152!!gm 
     153      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp 
    141154 
    142155      !                          !==  Convection  ==! 
     
    170183      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' ) 
    171184      ENDIF 
     185      !                                ! shear production term flag 
     186      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE. 
     187      ELSE                   ;   l_zdfsh2 = .TRUE. 
     188      ENDIF 
    172189 
    173190      !                          !== gravity wave-driven mixing  ==! 
     
    200217      ! 
    201218      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     219!!OMP      REAL(wp), DIMENSION(WRK_3D) ::   zsh2   ! shear production 
     220      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zsh2   ! shear production 
    202221      !! --------------------------------------------------------------------- 
     222      !       
     223!!OMP ===>>>>   Open-MP   to be defined elsewhere    
     224      ! 
     225      INTEGER ::   ARG_2D 
     226      ! 
     227#if defined key_vectopt_loop 
     228      k_Istr = 1   ;   k_Iend = jpi 
     229#else 
     230      k_Istr = 2   ;   k_Iend = jpim1 
     231#endif 
     232      k_Jstr = 2   ;   k_Jend = jpjm1 
     233      ! 
     234!!OMP <<<<===   end 
     235       
     236      ALLOCATE( zsh2(WRK_3D) ) 
     237       
     238       
    203239      ! 
    204240      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic) 
    205       !                        
     241      ! 
     242 
     243      !                          !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k) 
     244       
     245      IF( l_zdfsh2 )   &                        !* shear production at w-points (energy conserving form) 
     246         CALL zdf_sh2( ARG_2D, ub, vb, un, vn, avm_k,   &   ! <<== in  fields 
     247            &                                   zsh2    )   ! ==>> out field  : shear production 
     248 
     249      ! 
    206250      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points 
    207       CASE( np_RIC )   ;   CALL zdf_ric( kt )         ! Richardson number dependent Kz 
    208       CASE( np_TKE )   ;   CALL zdf_tke( kt )         ! TKE closure scheme for Kz 
    209       CASE( np_GLS )   ;   CALL zdf_gls( kt )         ! GLS closure scheme for Kz 
     251      CASE( np_RIC )   ;   CALL zdf_ric( ARG_2D, kt, gdept_n, zsh2, avm_k, avt_k )         ! Richardson number dependent Kz 
     252      CASE( np_TKE )   ;   CALL zdf_tke( ARG_2D, kt         , zsh2, avm_k, avt_k )         ! TKE closure scheme for Kz 
     253      CASE( np_GLS )   ;   CALL zdf_gls( ARG_2D, kt         , zsh2, avm_k, avt_k )         ! GLS closure scheme for Kz 
    210254      CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value) 
    211          avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
    212          avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
     255         ! avt_k and avm_k set one for all at initialisation phase 
     256!!gm         avt_k(WRK_3D) = rn_avt0 * wmask(WRK_3D) 
     257!!gm         avm_k(WRK_3D) = rn_avm0 * wmask(WRK_3D) 
    213258      END SELECT 
     259      !   
     260      !                          !==  ocean Kz  ==!   (avt, avs, avm) 
     261      ! 
     262      !                                         !* start from turbulent closure values 
     263      avt(WRK_2D,2:jpkm1) = avt_k(WRK_2D,2:jpkm1) 
     264      avm(WRK_2D,2:jpkm1) = avm_k(WRK_2D,2:jpkm1) 
    214265      ! 
    215266      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths 
    216267         DO jk = 2, nkrnf 
    217             avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk) 
     268            avt(WRK_2D,jk) = avt(WRK_2D,jk) + 2._wp * rn_avt_rnf * rnfmsk(WRK_2D) * wmask(WRK_2D,jk) 
    218269         END DO 
    219270      ENDIF 
    220271      ! 
    221       IF( ln_zdfevd )   CALL zdf_evd( kt )      !* convection: enhanced vertical eddy diffusivity 
     272      IF( ln_zdfevd )   CALL zdf_evd( ARG_2D, kt, avm, avt )      !* convection: enhanced vertical eddy diffusivity 
    222273      ! 
    223274      !                                         !* double diffusive mixing 
    224275      IF( ln_zdfddm ) THEN                            ! update avt and compute avs 
    225                         CALL zdf_ddm( kt ) 
     276                        CALL zdf_ddm( ARG_2D, kt, avm, avt, avs ) 
    226277      ELSE                                            ! same mixing on all tracers 
    227          avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 
     278         avs(WRK_3D) = avt(WRK_3D) 
    228279      ENDIF 
    229280      ! 
    230281      !                                         !* wave-induced mixing  
    231       IF( ln_zdfswm )   CALL zdf_swm( kt )            ! surface  wave (Qiao et al. 2004)  
    232       IF( ln_zdfiwm )   CALL zdf_iwm( kt )            ! internal wave (de Lavergne et al 2017) 
     282      IF( ln_zdfswm )   CALL zdf_swm( ARG_2D, kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)  
     283      IF( ln_zdfiwm )   CALL zdf_iwm( ARG_2D, kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017) 
    233284 
    234285 
    235286      !                                         !* Lateral boundary conditions (sign unchanged) 
    236       CALL lbc_lnk( avm, 'W', 1. ) 
    237       CALL lbc_lnk( avt, 'W', 1. ) 
    238       CALL lbc_lnk( avs, 'W', 1. ) 
    239       ! 
    240 !!gm  ===>>>   TO BE REMOVED  
    241       DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
    242          DO jj = 2, jpjm1 
    243             DO ji = 2, jpim1 
    244                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
    245                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    246             END DO 
    247          END DO 
    248       END DO 
    249       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    250 !!gm end 
    251  
     287      CALL lbc_lnk( avm_k, 'W', 1. ) 
     288      CALL lbc_lnk( avt_k, 'W', 1. ) 
     289      CALL lbc_lnk( avm  , 'W', 1. ) 
     290      CALL lbc_lnk( avt  , 'W', 1. )                  !!gm  a priori only avm_k and avm are required 
     291      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  To be tested 
     292      ! 
    252293 
    253294      CALL zdf_mxl( kt )                        !* mixed layer depth, and level 
    254295 
    255 !!gm  !==>> to be moved in zdftke & zdfgls     
    256                                                       ! write TKE or GLS information in the restart file 
    257       IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kt, 'WRITE' ) 
    258       IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
    259       !       
     296 
     297      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file 
     298         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' ) 
     299         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
     300         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' )  
     301      ENDIF 
     302      ! 
     303      DEALLOCATE( zsh2 ) 
     304      ! 
    260305   END SUBROUTINE zdf_phy 
    261306 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r7991 r8055  
    1616 
    1717   !!---------------------------------------------------------------------- 
     18   !!   zdf_ric_init  : initialization, namelist read, & parameters control 
    1819   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number 
    19    !!   zdf_ric_init  : initialization, namelist read, & parameters control 
     20   !!   ric_rst       : read/write RIC restart in ocean restart file 
    2021   !!---------------------------------------------------------------------- 
    2122   USE oce            ! ocean dynamics and tracers variables 
    2223   USE dom_oce        ! ocean space and time domain variables 
    2324   USE zdf_oce        ! vertical physics: variables 
    24    USE zdfsh2         ! vertical physics: shear production term of TKE 
    2525   USE phycst         ! physical constants 
    2626   USE sbc_oce,  ONLY :   taum 
     
    2828   ! 
    2929   USE in_out_manager ! I/O manager 
    30    USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    31    USE lib_mpp        ! MPP library 
     30   USE iom            ! I/O manager library 
    3231   USE timing         ! Timing 
    3332   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     
    3736   PRIVATE 
    3837 
    39    PUBLIC   zdf_ric         ! called by step.F90 
    40    PUBLIC   zdf_ric_init    ! called by opa.F90 
     38   PUBLIC   zdf_ric         ! called by zdfphy.F90 
     39   PUBLIC   ric_rst         ! called by zdfphy.F90 
     40   PUBLIC   zdf_ric_init    ! called by nemogcm.F90 
    4141 
    4242   !                        !!* Namelist namzdf_ric : Richardson number dependent Kz * 
     
    5252 
    5353   !! * Substitutions 
    54 #  include "vectopt_loop_substitute.h90" 
     54#  include "domain_substitute.h90" 
    5555   !!---------------------------------------------------------------------- 
    5656   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     
    5959   !!---------------------------------------------------------------------- 
    6060CONTAINS 
    61  
    62    SUBROUTINE zdf_ric( kt ) 
    63       !!---------------------------------------------------------------------- 
    64       !!                 ***  ROUTINE zdfric  *** 
    65       !!                     
    66       !! ** Purpose :   Compute the before eddy viscosity and diffusivity as 
    67       !!                a function of the local richardson number. 
    68       !! 
    69       !! ** Method  :   Local richardson number dependent formulation of the  
    70       !!                vertical eddy viscosity and diffusivity coefficients.  
    71       !!                The eddy coefficients are given by: 
    72       !!                    avm = avm0 + avmb 
    73       !!                    avt = avm0 / (1 + rn_alp*ri) 
    74       !!                with ri  = N^2 / dz(u)**2 
    75       !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 
    76       !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric 
    77       !!      Where ri is the before local Richardson number, 
    78       !!            rn_avmri is the maximum value reaches by avm and avt  
    79       !!            avmb and avtb are the background (or minimum) values 
    80       !!            and rn_alp, nn_ric are adjustable parameters. 
    81       !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s 
    82       !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2. 
    83       !!      a numerical threshold is impose on the vertical shear (1.e-20) 
    84       !!      As second step compute Ekman depth from wind stress forcing 
    85       !!      and apply namelist provided vertical coeff within this depth. 
    86       !!      The Ekman depth is: 
    87       !!              Ustar = SQRT(Taum/rho0) 
    88       !!              ekd= rn_ekmfc * Ustar / f0 
    89       !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation 
    90       !!      of the above equation indicates the value is somewhat arbitrary; therefore 
    91       !!      we allow the freedom to increase or decrease this value, if the 
    92       !!      Ekman depth estimate appears too shallow or too deep, respectively. 
    93       !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the 
    94       !!      namelist 
    95       !!        N.B. the mask are required for implicit scheme, and surface 
    96       !!      and bottom value already set in zdfphy.F90 
    97       !! 
    98       !! ** Action  :   avm, avt  mixing coeff (inner domain values only) 
    99       !! 
    100       !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 
    101       !!              PFJ Lermusiaux 2001. 
    102       !!---------------------------------------------------------------------- 
    103       INTEGER, INTENT(in) ::   kt   ! ocean time-step 
    104       !! 
    105       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    106       LOGICAL  ::   ll_av_weight = .TRUE.      ! local logical 
    107       REAL(wp) ::   zcfRi, zav, zustar   ! local scalars 
    108       REAL(wp), DIMENSION(jpi,jpj)     ::   zh_ekm   ! 2D workspace 
    109       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2     ! 3D     - 
    110       !!---------------------------------------------------------------------- 
    111       ! 
    112       IF( nn_timing == 1 )   CALL timing_start('zdf_ric') 
    113       ! 
    114       !                       !==  avm and avt = F(Richardson number)  ==! 
    115       ! 
    116       !                          !* Shear production at uw- and vw-points (energy conserving form) 
    117       CALL zdf_sh2( zsh2 ) 
    118       ! 
    119       DO jk = 2, jpkm1           !* Vertical eddy viscosity and diffusivity coefficients 
    120          DO jj = 1, jpjm1 
    121             DO ji = 1, jpim1 
    122                !                          ! coefficient = F(richardson number) (avm-weighted Ri) 
    123                zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) )  ) 
    124                zav   = rn_avmri * zcfRi**nn_ric 
    125                ! 
    126                !                          ! avm and avt coefficients 
    127                avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk) 
    128                avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk) 
    129             END DO 
    130          END DO 
    131       END DO 
    132       ! 
    133 !!gm BUG <<<<====  This param can't work at low latitude  
    134 !!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 
    135       ! 
    136       IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==! 
    137          ! 
    138          DO jj = 2, jpjm1        !* Ekman depth 
    139             DO ji = 2, jpim1 
    140                zustar = SQRT( taum(ji,jj) * r1_rau0 ) 
    141                zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth 
    142                zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax )  )   ! set allowed rang 
    143             END DO 
    144          END DO 
    145          DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer 
    146             DO jj = 2, jpjm1 
    147                DO ji = 2, jpim1 
    148                   IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 
    149                      avm(ji,jj,jk) = MAX(  avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk) 
    150                      avt(ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk) 
    151                   ENDIF 
    152                END DO 
    153             END DO 
    154          END DO 
    155       ENDIF 
    156       ! 
    157       IF( nn_timing == 1 )   CALL timing_stop('zdf_ric') 
    158       ! 
    159    END SUBROUTINE zdf_ric 
    160  
    16161 
    16262   SUBROUTINE zdf_ric_init 
     
    205105      ENDIF 
    206106      ! 
    207       DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value 
    208          avt(:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    209          avm(:,:,jk) = avmb(jk) * umask(:,:,jk) 
     107      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files 
     108      ! 
     109   END SUBROUTINE zdf_ric_init 
     110 
     111 
     112   SUBROUTINE zdf_ric( ARG_2D, kt, pdept, p_sh2, p_avm, p_avt ) 
     113      !!---------------------------------------------------------------------- 
     114      !!                 ***  ROUTINE zdfric  *** 
     115      !!                     
     116      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as 
     117      !!                a function of the local richardson number. 
     118      !! 
     119      !! ** Method  :   Local richardson number dependent formulation of the  
     120      !!                vertical eddy viscosity and diffusivity coefficients.  
     121      !!                The eddy coefficients are given by: 
     122      !!                    avm = avm0 + avmb 
     123      !!                    avt = avm0 / (1 + rn_alp*ri) 
     124      !!                with Ri  = N^2 / dz(u)**2 
     125      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 
     126      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 
     127      !!                where ri is the before local Richardson number, 
     128      !!                rn_avmri is the maximum value reaches by avm and avt  
     129      !!                and rn_alp, nn_ric are adjustable parameters. 
     130      !!                Typical values : rn_alp=5. and nn_ric=2. 
     131      !! 
     132      !!      As second step compute Ekman depth from wind stress forcing 
     133      !!      and apply namelist provided vertical coeff within this depth. 
     134      !!      The Ekman depth is: 
     135      !!              Ustar = SQRT(Taum/rho0) 
     136      !!              ekd= rn_ekmfc * Ustar / f0 
     137      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation 
     138      !!      of the above equation indicates the value is somewhat arbitrary; therefore 
     139      !!      we allow the freedom to increase or decrease this value, if the 
     140      !!      Ekman depth estimate appears too shallow or too deep, respectively. 
     141      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the 
     142      !!      namelist 
     143      !!        N.B. the mask are required for implicit scheme, and surface 
     144      !!      and bottom value already set in zdfphy.F90 
     145      !! 
     146      !! ** Action  :   avm, avt  mixing coeff (inner domain values only) 
     147      !! 
     148      !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 
     149      !!              PFJ Lermusiaux 2001. 
     150      !!---------------------------------------------------------------------- 
     151      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     152      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time-step 
     153      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdept          ! depth of t-point  [m] 
     154      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term 
     155      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points) 
     156      !! 
     157      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
     158      REAL(wp) ::   zcfRi, zav, zustar, zed     ! local scalars 
     159      REAL(wp), DIMENSION(WRK_2D) ::   zh_ekm   ! 2D workspace 
     160      !!---------------------------------------------------------------------- 
     161      ! 
     162      IF( nn_timing == 1 )   CALL timing_start('zdf_ric') 
     163      ! 
     164      !                       !==  avm and avt = F(Richardson number)  ==! 
     165      DO jk = 2, jpkm1 
     166         DO jj = k_Jstr, k_Jend 
     167            DO ji = k_Jstr, k_Iend        ! coefficient = F(richardson number) (avm-weighted Ri) 
     168               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  ) 
     169               zav   = rn_avmri * zcfRi**nn_ric 
     170               !                          ! avm and avt coefficients 
     171               p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk) 
     172               p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk) 
     173            END DO 
     174         END DO 
    210175      END DO 
    211176      ! 
    212    END SUBROUTINE zdf_ric_init 
     177!!gm BUG <<<<====  This param can't work at low latitude  
     178!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 
     179      ! 
     180      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==! 
     181         ! 
     182         DO jj = k_Jstr, k_Jend     !* Ekman depth 
     183            DO ji = k_Jstr, k_Iend 
     184               zustar = SQRT( taum(ji,jj) * r1_rau0 ) 
     185               zed    = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth 
     186               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zed , rn_mldmax )  )      ! set allowed range 
     187            END DO 
     188         END DO 
     189         DO jk = 2, jpkm1           !* minimum mixing coeff. within the Ekman layer 
     190            DO jj = k_Jstr, k_Jend 
     191               DO ji = k_Jstr, k_Iend 
     192                  IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 
     193                     p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk) 
     194                     p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk) 
     195                  ENDIF 
     196               END DO 
     197            END DO 
     198         END DO 
     199      ENDIF 
     200      ! 
     201      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric') 
     202      ! 
     203   END SUBROUTINE zdf_ric 
     204 
     205 
     206   SUBROUTINE ric_rst( kt, cdrw ) 
     207      !!--------------------------------------------------------------------- 
     208      !!                   ***  ROUTINE ric_rst  *** 
     209      !!                      
     210      !! ** Purpose :   Read or write TKE file (en) in restart file 
     211      !! 
     212      !! ** Method  :   use of IOM library 
     213      !!                if the restart does not contain TKE, en is either  
     214      !!                set to rn_emin or recomputed  
     215      !!---------------------------------------------------------------------- 
     216      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     217      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     218      ! 
     219      INTEGER ::   jit, jk    ! dummy loop indices 
     220      INTEGER ::   id1, id2   ! local integers 
     221      !!---------------------------------------------------------------------- 
     222      ! 
     223      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     224         !                                   ! --------------- 
     225         !           !* Read the restart file 
     226         IF( ln_rstart ) THEN 
     227            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
     228            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
     229            ! 
     230            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it 
     231               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
     232               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
     233            ENDIF 
     234         ENDIF 
     235         !           !* otherwise Kz already set to the background value in zdf_phy_init 
     236         ! 
     237      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     238         !                                   ! ------------------- 
     239         IF(lwp) WRITE(numout,*) '---- ric-rst ----' 
     240         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
     241         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
     242         ! 
     243      ENDIF 
     244      ! 
     245   END SUBROUTINE ric_rst 
    213246 
    214247   !!====================================================================== 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90

    r7997 r8055  
    55   !!===================================================================== 
    66   !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code 
    7    !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm 
     7   !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm + ready for coarse-grained Open-MP 
    88   !!---------------------------------------------------------------------- 
    99 
     
    1111   !!   zdf_sh2       : compute mixing the shear production term of TKE 
    1212   !!---------------------------------------------------------------------- 
    13    USE oce            ! ocean: shared variables 
    1413   USE dom_oce        ! domain: ocean 
    15    USE zdf_oce        ! vertical physics: variables 
    1614   ! 
    1715   USE in_out_manager ! I/O manager 
     
    2321 
    2422   PUBLIC   zdf_sh2        ! called by zdftke, zdfglf, and zdfric 
    25     
     23 
     24   !! * Substitutions 
     25#  include "domain_substitute.h90"    
    2626   !!---------------------------------------------------------------------- 
    2727   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
     
    3131CONTAINS 
    3232 
    33    SUBROUTINE zdf_sh2(  psh2  )  
     33   SUBROUTINE zdf_sh2( ARG_2D, pub, pvb, pun, pvn, p_avm, p_sh2  )  
    3434      !!---------------------------------------------------------------------- 
    3535      !!                   ***  ROUTINE zdf_sh2  *** 
     
    4646      !!                NB: wet-point only horizontal averaging of shear 
    4747      !! 
    48       !! ** Action  : - psh2 shear prod. term at w-point (interior ocean domain only) 
    49       !! 
     48      !! ** Action  : - p_sh2 shear prod. term at w-point (inner domain only) 
     49      !!                                                   ***** 
    5050      !! References :   Bruchard, OM 2002 
    5151      !! --------------------------------------------------------------------- 
    52       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   psh2   ! shear production of TKE (w-points) 
     52      INTEGER                    , INTENT(in   ) ::   ARG_2D               ! inner domain start-end i-indices 
     53      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities 
     54      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points) 
     55      REAL(wp), DIMENSION(WRK_3D), INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points) 
    5356      ! 
    54       INTEGER  ::   ji, jj, jk   ! dummy loop arguments 
    55       REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace 
     57      INTEGER  ::   ji, jj, jk                           ! dummy loop indices 
     58      REAL(wp), DIMENSION(WRK_2De(-1,)) ::   zsh2u, zsh2v   ! 2D workspace 
    5659      !!-------------------------------------------------------------------- 
     60      ! 
    5761      IF( nn_timing == 1 )  CALL timing_start('zdf_sh2') 
    5862      ! 
    5963      DO jk = 2, jpkm1 
    60          DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form) 
    61             DO ji = 1, jpim1 
    62                zsh2u(ji,jj) = ( avm(ji+1,jj,jk) + avm(ji,jj,jk) ) & 
    63                   &         * (  un(ji,jj,jk-1) -  un(ji,jj,jk) ) & 
    64                   &         * (  ub(ji,jj,jk-1) -  ub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)  
    65                zsh2v(ji,jj) = ( avm(ji,jj+1,jk) + avm(ji,jj,jk) ) & 
    66                   &         * (  vn(ji,jj,jk-1) -  vn(ji,jj,jk) ) & 
    67                   &         * (  vb(ji,jj,jk-1) -  vb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk) 
     64         DO jj = k_Jstr-1, k_Jend      !* 2 x shear production at uw- and vw-points (energy conserving form) 
     65            DO ji = k_Istr-1, k_Iend 
     66               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 
     67                  &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) & 
     68                  &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk) 
     69               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 
     70                  &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) & 
     71                  &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk) 
    6872            END DO 
    6973         END DO 
    70          DO jj = 2, jpjm1        !* shear production at w-point 
    71             DO ji = 2, jpim1           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 
    72                psh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   & 
    73                   &                      + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   ) 
     74         DO jj = k_Jstr, k_Jend        !* shear production at w-point 
     75            DO ji = k_Jstr, k_Iend           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 
     76               p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   & 
     77                  &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   ) 
    7478            END DO 
    7579         END DO 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90

    r7990 r8055  
    2727   PUBLIC zdf_swm_init    ! routine called in zdf_phy_init 
    2828 
     29   !! * Substitutions 
     30#  include "domain_substitute.h90"    
    2931   !!---------------------------------------------------------------------- 
    3032   !! NEMO/OPA 4.0 , NEMO Consortium (2017)  
     
    3436CONTAINS 
    3537 
    36    SUBROUTINE zdf_swm( kt ) 
     38   SUBROUTINE zdf_swm( ARG_2D, kt, p_avm, p_avt, p_avs ) 
    3739      !!--------------------------------------------------------------------- 
    3840      !!                     ***  ROUTINE zdf_swm *** 
     
    5153      !! reference : Qiao et al. GRL, 2004 
    5254      !!--------------------------------------------------------------------- 
    53       INTEGER, INTENT(in) ::   kt   ! ocean time step 
     55      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     56      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     57      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
     58      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    5459      ! 
    5560      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    5964      zcoef = 1._wp * 0.353553_wp 
    6065      DO jk = 2, jpkm1 
    61          DO jj = 2, jpjm1 
    62             DO ji = 2, jpim1 
     66         DO jj = k_Jstr, k_Jend 
     67            DO ji = k_Istr, k_Iend 
    6368               zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk) 
    6469               ! 
    65                avt(ji,jj,jk) = avt(ji,jj,jk) + zqb 
    66                avs(ji,jj,jk) = avs(ji,jj,jk) + zqb 
    67                avm(ji,jj,jk) = avm(ji,jj,jk) + zqb 
     70               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zqb 
     71               p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zqb 
     72               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zqb 
    6873            END DO 
    6974         END DO 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r7991 r8055  
    4343   USE sbc_oce        ! surface boundary condition: ocean 
    4444   USE zdf_oce        ! vertical physics: ocean variables 
    45    USE zdfsh2         ! vertical physics: shear production term of TKE 
    4645   USE zdfmxl         ! vertical physics: mixed layer 
    47    USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    48    USE prtctl         ! Print control 
    49    USE in_out_manager ! I/O manager 
    50    USE iom            ! I/O manager library 
    51    USE lib_mpp        ! MPP library 
    52    USE wrk_nemo       ! work arrays 
    53    USE timing         ! Timing 
    54    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    5546#if defined key_agrif 
    5647   USE agrif_opa_interp 
    5748   USE agrif_opa_update 
    5849#endif 
     50   ! 
     51   USE in_out_manager ! I/O manager 
     52   USE iom            ! I/O manager library 
     53   USE lib_mpp        ! MPP library 
     54   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     55   USE prtctl         ! Print control 
     56   USE wrk_nemo       ! work arrays 
     57   USE timing         ! Timing 
     58   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    5959 
    6060   IMPLICIT NONE 
     
    8787   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    8888 
    89    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    90    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    91    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation 
     89   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau) 
     90   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation 
     91   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation 
    9292 
    9393   !! * Substitutions 
    94 #  include "vectopt_loop_substitute.h90" 
     94#  include "domain_substitute.h90"    
    9595   !!---------------------------------------------------------------------- 
    9696   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     
    112112 
    113113 
    114    SUBROUTINE zdf_tke( kt ) 
     114   SUBROUTINE zdf_tke( ARG_2D, kt, p_sh2, p_avm, p_avt ) 
    115115      !!---------------------------------------------------------------------- 
    116116      !!                   ***  ROUTINE zdf_tke  *** 
     
    157157      !!              Bruchard OM 2002 
    158158      !!---------------------------------------------------------------------- 
    159       INTEGER, INTENT(in) ::   kt   ! ocean time step 
     159      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     160      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     161      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term 
     162      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    160163      !!---------------------------------------------------------------------- 
    161164      ! 
     
    165168#endif 
    166169      ! 
    167       avt(:,:,:) = avt_k(:,:,:)     ! restore before Kz computed with TKE alone 
    168       avm(:,:,:) = avm_k(:,:,:)  
    169       ! 
    170       CALL tke_tke                  ! now tke (en) 
    171       ! 
    172       CALL tke_avn                  ! now avt, avm, dissl 
    173       ! 
    174       avt_k(:,:,:) = avt(:,:,:)     ! store Kz computed with TKE alone 
    175       avm_k(:,:,:) = avm(:,:,:)  
     170      CALL tke_tke( ARG_2D, gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )        ! now tke (en) 
     171      ! 
     172 
     173!!gm not sure we need this lbc ....      ==>>>  certainly NOT !!! 
     174      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     175      ! 
     176      CALL tke_avn( ARG_2D, gdepw_n, e3t_n, e3w_n, p_avm, p_avt )        ! now avt, avm, dissl 
    176177      ! 
    177178#if defined key_agrif 
     
    179180      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    180181#endif       
    181       !  
    182       IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' ) 
    183182      ! 
    184183  END SUBROUTINE zdf_tke 
    185184 
    186185 
    187    SUBROUTINE tke_tke 
     186   SUBROUTINE tke_tke( ARG_2D, pdepw, p_e3t, p_e3w, p_sh2    & 
     187      &                                    , p_avm, p_avt ) 
    188188      !!---------------------------------------------------------------------- 
    189189      !!                   ***  ROUTINE tke_tke  *** 
     
    200200      !! ** Action  : - en : now turbulent kinetic energy) 
    201201      !! --------------------------------------------------------------------- 
    202       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
     202      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     203      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! vertical eddy viscosity (w-points) 
     204      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! vertical eddy viscosity (w-points) 
     205      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term 
     206      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   !  vertical eddy viscosity (w-points) 
     207      ! 
     208      INTEGER ::   ji, jj, jk              ! dummy loop arguments 
    203209!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! local integers 
    204210!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          !   -      - 
    205211!!bfr      REAL(wp) ::   zebot                           ! local scalars 
    206       REAL(wp) ::   zrhoa  = 1.22            ! Air density kg/m3 
    207       REAL(wp) ::   zcdrag = 1.5e-3          ! drag coefficient 
    208       REAL(wp) ::   zbbrau, zri              ! local scalars 
    209       REAL(wp) ::   zfact1, zfact2, zfact3   !   -         - 
    210       REAL(wp) ::   ztx2  , zty2  , zcof     !   -         - 
    211       REAL(wp) ::   ztau  , zdif             !   -         - 
    212       REAL(wp) ::   zus   , zwlc  , zind     !   -         - 
    213       REAL(wp) ::   zzd_up, zzd_lw           !   -         - 
    214       INTEGER , DIMENSION(jpi,jpj) ::   imlc 
    215       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc 
    216       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw 
    217       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2 
     212      REAL(wp)::   zrhoa  = 1.22            ! Air density kg/m3 
     213      REAL(wp)::   zcdrag = 1.5e-3          ! drag coefficient 
     214      REAL(wp)::   zbbrau, zri              ! local scalars 
     215      REAL(wp)::   zfact1, zfact2, zfact3   !   -         - 
     216      REAL(wp)::   ztx2  , zty2  , zcof     !   -         - 
     217      REAL(wp)::   ztau  , zdif             !   -         - 
     218      REAL(wp)::   zus   , zwlc  , zind     !   -         - 
     219      REAL(wp)::   zzd_up, zzd_lw           !   -         - 
     220      INTEGER , DIMENSION(WRK_2D) ::   imlc 
     221      REAL(wp), DIMENSION(WRK_2D) ::   zhlc 
     222      REAL(wp), DIMENSION(WRK_3D) ::   zpelc, zdiag, zd_up, zd_lw 
    218223      !!-------------------------------------------------------------------- 
    219224      ! 
    220225      IF( nn_timing == 1 )  CALL timing_start('tke_tke') 
    221       ! 
    222       CALL wrk_alloc( jpi,jpj,       zhlc )  
    223       CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw )  
    224226      ! 
    225227      zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     
    233235      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    234236      IF ( ln_isfcav ) THEN 
    235          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
    236             DO ji = fs_2, fs_jpim1   ! vector opt. 
     237         DO jj = k_Jstr, k_Jend            ! en(mikt(ji,jj))   = rn_emin 
     238            DO ji = k_Istr, k_Iend 
    237239               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 
    238240            END DO 
    239241         END DO 
    240242      END IF 
    241       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242          DO ji = fs_2, fs_jpim1   ! vector opt. 
     243      DO jj = k_Jstr, k_Jend            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
     244         DO ji = k_Istr, k_Iend 
    243245            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    244246         END DO 
     
    256258      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    257259      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
    258 !!    DO jj = 2, jpjm1 
    259 !!       DO ji = fs_2, fs_jpim1   ! vector opt. 
     260!!    DO jj = k_Jstr, k_Jend 
     261!!       DO ji = k_Jstr, k_Iend 
    260262!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 
    261263!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) ) 
     
    269271      ! 
    270272      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    271       IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002) 
     273      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002) 
    272274         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    273275         ! 
    274276         !                        !* total energy produce by LC : cumulative sum over jk 
    275          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1) 
     277         zpelc(WRK_2D,1) =  MAX( rn2b(WRK_2D,1), 0._wp ) * pdepw(WRK_2D,1) * p_e3w(WRK_2D,1) 
    276278         DO jk = 2, jpk 
    277             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk) 
     279            zpelc(WRK_2D,jk) = zpelc(WRK_2D,jk-1)   & 
     280               &             + MAX( rn2b(WRK_2D,jk), 0._wp )   & 
     281               &             * pdepw(WRK_2D,jk) * p_e3w(WRK_2D,jk) 
    278282         END DO 
    279283         !                        !* finite Langmuir Circulation depth 
    280284         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    281          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
     285         imlc(WRK_2D) = mbkt(WRK_2D) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    282286         DO jk = jpkm1, 2, -1 
    283             DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us  
    284                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1) 
     287            DO jj = k_Jstr, k_Jend               ! Last w-level at which zpelc>=0.5*us*us  
     288               DO ji = k_Istr, k_Iend            !      with us=0.016*wind(starting from jpk-1) 
    285289                  zus  = zcof * taum(ji,jj) 
    286290                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
     
    289293         END DO 
    290294         !                               ! finite LC depth 
    291          DO jj = 1, jpj  
    292             DO ji = 1, jpi 
    293                zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj)) 
     295         DO jj = k_Jstr, k_Jend  
     296            DO ji = k_Istr, k_Iend 
     297               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 
    294298            END DO 
    295299         END DO 
    296300         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    297301         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    298             DO jj = 2, jpjm1 
    299                DO ji = fs_2, fs_jpim1   ! vector opt. 
     302            DO jj = k_Jstr, k_Jend 
     303               DO ji = k_Istr, k_Iend 
    300304                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    301305                  !                                           ! vertical velocity due to LC 
    302                   zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) ) 
    303                   zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) ) 
     306                  zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 
     307                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    304308                  !                                           ! TKE Langmuir circulation source term 
    305309                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
     
    318322      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 
    319323      ! 
    320       !                           !* Shear production at uw- and vw-points (energy conserving form) 
    321       CALL zdf_sh2( zsh2 ) 
    322       ! 
    323324      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri ) 
    324325         DO jk = 2, jpkm1 
    325             DO jj = 2, jpjm1 
    326                DO ji = 2, jpim1 
     326            DO jj = k_Jstr, k_Jend 
     327               DO ji = k_Istr, k_Iend 
    327328                  !                             ! local Richardson number 
    328                   zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear ) 
     329                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 
    329330                  !                             ! inverse of Prandtl number 
    330331                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
     
    335336      !          
    336337      DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    337          DO jj = 2, jpjm1 
    338             DO ji = fs_2, fs_jpim1   ! vector opt. 
     338         DO jj = k_Jstr, k_Jend 
     339            DO ji = k_Istr, k_Iend 
    339340               zcof   = zfact1 * tmask(ji,jj,jk) 
    340341               !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
    341342               !                                   ! eddy coefficient (ensure numerical stability) 
    342                zzd_up = zcof * MAX(   avm(ji,jj,jk+1) +   avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
    343                   &          /    ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  ) 
    344                zzd_lw = zcof * MAX(   avm(ji,jj,jk  ) +   avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
    345                   &          /    ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  ) 
     343               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  )  ) 
     345               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  )  ) 
    346347               ! 
    347348               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    350351               ! 
    351352               !                                   ! right hand side in en 
    352                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zsh2(ji,jj,jk)                           &   ! shear 
    353                   &                                 - avt(ji,jj,jk) * rn2(ji,jj,jk)            &   ! stratification 
     353               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
     354                  &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
    354355                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
    355356                  &                                ) * wmask(ji,jj,jk) 
     
    359360      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    360361      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    361          DO jj = 2, jpjm1 
    362             DO ji = fs_2, fs_jpim1    ! vector opt. 
     362         DO jj = k_Jstr, k_Jend 
     363            DO ji = k_Istr, k_Iend 
    363364               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    364365            END DO 
    365366         END DO 
    366367      END DO 
    367       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    368          DO ji = fs_2, fs_jpim1   ! vector opt. 
     368      DO jj = k_Jstr, k_Jend                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     369         DO ji = k_Istr, k_Iend 
    369370            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    370371         END DO 
    371372      END DO 
    372373      DO jk = 3, jpkm1 
    373          DO jj = 2, jpjm1 
    374             DO ji = fs_2, fs_jpim1    ! vector opt. 
     374         DO jj = k_Jstr, k_Jend 
     375            DO ji = k_Istr, k_Iend 
    375376               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    376377            END DO 
    377378         END DO 
    378379      END DO 
    379       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    380          DO ji = fs_2, fs_jpim1   ! vector opt. 
     380      DO jj = k_Jstr, k_Jend                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     381         DO ji = k_Istr, k_Iend 
    381382            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    382383         END DO 
    383384      END DO 
    384385      DO jk = jpk-2, 2, -1 
    385          DO jj = 2, jpjm1 
    386             DO ji = fs_2, fs_jpim1    ! vector opt. 
     386         DO jj = k_Jstr, k_Jend 
     387            DO ji = k_Istr, k_Iend 
    387388               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    388389            END DO 
     
    390391      END DO 
    391392      DO jk = 2, jpkm1                             ! set the minimum value of tke 
    392          DO jj = 2, jpjm1 
    393             DO ji = fs_2, fs_jpim1   ! vector opt. 
     393         DO jj = k_Jstr, k_Jend 
     394            DO ji = k_Istr, k_Iend 
    394395               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    395396            END DO 
     
    405406      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    406407         DO jk = 2, jpkm1 
    407             DO jj = 2, jpjm1 
    408                DO ji = fs_2, fs_jpim1   ! vector opt. 
    409                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   & 
     408            DO jj = k_Jstr, k_Jend 
     409               DO ji = k_Istr, k_Iend 
     410                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    410411                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    411412               END DO 
     
    413414         END DO 
    414415      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction) 
    415          DO jj = 2, jpjm1 
    416             DO ji = fs_2, fs_jpim1   ! vector opt. 
     416         DO jj = k_Jstr, k_Jend 
     417            DO ji = k_Istr, k_Iend 
    417418               jk = nmln(ji,jj) 
    418                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   & 
     419               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    419420                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    420421            END DO 
     
    422423      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
    423424         DO jk = 2, jpkm1 
    424             DO jj = 2, jpjm1 
    425                DO ji = fs_2, fs_jpim1   ! vector opt. 
     425            DO jj = k_Jstr, k_Jend 
     426               DO ji = k_Istr, k_Iend 
    426427                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    427428                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
     
    429430                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    430431                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    431                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   & 
     432                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    432433                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    433434               END DO 
     
    435436         END DO 
    436437      ENDIF 
    437 !!gm not sure we need this lbc .... 
    438       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    439       ! 
    440       CALL wrk_dealloc( jpi,jpj,       zhlc )  
    441       CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw )  
    442438      ! 
    443439      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
     
    446442 
    447443 
    448    SUBROUTINE tke_avn 
     444   SUBROUTINE tke_avn( ARG_2D, pdepw, p_e3t, p_e3w, p_avm, p_avt ) 
    449445      !!---------------------------------------------------------------------- 
    450446      !!                   ***  ROUTINE tke_avn  *** 
     
    480476      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 
    481477      !!---------------------------------------------------------------------- 
     478      INTEGER                   , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     479      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! vertical eddy viscosity (w-points) 
     480      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w     ! vertical eddy viscosity (w-points) 
     481      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   !  vertical eddy viscosity (w-points) 
     482      ! 
    482483      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    483484      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
    484485      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
    485486      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
    486       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld 
     487      REAL(wp), DIMENSION(WRK_3D) ::  zmxlm, zmxld 
    487488      !!-------------------------------------------------------------------- 
    488489      ! 
    489490      IF( nn_timing == 1 )  CALL timing_start('tke_avn') 
    490491 
    491       CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
    492  
    493492      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    494493      !                     !  Mixing length 
     
    498497      ! 
    499498      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
    500       zmxlm(:,:,:)  = rmxl_min     
    501       zmxld(:,:,:)  = rmxl_min 
     499      zmxlm(WRK_3D)  = rmxl_min     
     500      zmxld(WRK_3D)  = rmxl_min 
    502501      ! 
    503502      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    504          DO jj = 2, jpjm1 
    505             DO ji = fs_2, fs_jpim1 
    506                zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     503         zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     504         DO jj = k_Jstr, k_Jend 
     505            DO ji = k_Istr, k_Iend 
    507506               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    508507            END DO 
    509508         END DO 
    510509      ELSE  
    511          zmxlm(:,:,1) = rn_mxl0 
     510         zmxlm(WRK_2D,1) = rn_mxl0 
    512511      ENDIF 
    513512      ! 
    514513      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    515          DO jj = 2, jpjm1 
    516             DO ji = fs_2, fs_jpim1   ! vector opt. 
     514         DO jj = k_Jstr, k_Jend 
     515            DO ji = k_Istr, k_Iend 
    517516               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    518517               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
     
    523522      !                     !* Physical limits for the mixing length 
    524523      ! 
    525       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value  
    526       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
     524      zmxld(WRK_2D, 1 ) = zmxlm(WRK_2D,1)   ! surface set to the minimum value  
     525      zmxld(WRK_2D,jpk) = rmxl_min          ! last level  set to the minimum value 
    527526      ! 
    528527      SELECT CASE ( nn_mxl ) 
    529528      ! 
    530529 !!gm Not sure of that coding for ISF.... 
    531       ! where wmask = 0 set zmxlm == e3w_n 
     530      ! where wmask = 0 set zmxlm == p_e3w 
    532531      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    533532         DO jk = 2, jpkm1 
    534             DO jj = 2, jpjm1 
    535                DO ji = fs_2, fs_jpim1   ! vector opt. 
    536                   zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    537                   &            gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) ) 
     533            DO jj = k_Jstr, k_Jend 
     534               DO ji = k_Istr, k_Iend 
     535                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
     536                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 
    538537                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
    539                   zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
    540                   zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
     538                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
     539                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
    541540               END DO 
    542541            END DO 
     
    545544      CASE ( 1 )           ! bounded by the vertical scale factor 
    546545         DO jk = 2, jpkm1 
    547             DO jj = 2, jpjm1 
    548                DO ji = fs_2, fs_jpim1   ! vector opt. 
    549                   zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) ) 
     546            DO jj = k_Jstr, k_Jend 
     547               DO ji = k_Istr, k_Iend 
     548                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    550549                  zmxlm(ji,jj,jk) = zemxl 
    551550                  zmxld(ji,jj,jk) = zemxl 
     
    556555      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    557556         DO jk = 2, jpkm1         ! from the surface to the bottom : 
    558             DO jj = 2, jpjm1 
    559                DO ji = fs_2, fs_jpim1   ! vector opt. 
    560                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
     557            DO jj = k_Jstr, k_Jend 
     558               DO ji = k_Istr, k_Iend 
     559                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    561560               END DO 
    562561            END DO 
    563562         END DO 
    564563         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
    565             DO jj = 2, jpjm1 
    566                DO ji = fs_2, fs_jpim1   ! vector opt. 
    567                   zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
     564            DO jj = k_Jstr, k_Jend 
     565               DO ji = k_Istr, k_Iend 
     566                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    568567                  zmxlm(ji,jj,jk) = zemxl 
    569568                  zmxld(ji,jj,jk) = zemxl 
     
    574573      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    575574         DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
    576             DO jj = 2, jpjm1 
    577                DO ji = fs_2, fs_jpim1   ! vector opt. 
    578                   zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
     575            DO jj = k_Jstr, k_Jend 
     576               DO ji = k_Istr, k_Iend 
     577                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    579578               END DO 
    580579            END DO 
    581580         END DO 
    582581         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
    583             DO jj = 2, jpjm1 
    584                DO ji = fs_2, fs_jpim1   ! vector opt. 
    585                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
     582            DO jj = k_Jstr, k_Jend 
     583               DO ji = k_Istr, k_Iend 
     584                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    586585               END DO 
    587586            END DO 
    588587         END DO 
    589588         DO jk = 2, jpkm1 
    590             DO jj = 2, jpjm1 
    591                DO ji = fs_2, fs_jpim1   ! vector opt. 
     589            DO jj = k_Jstr, k_Jend 
     590               DO ji = k_Istr, k_Iend 
    592591                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) ) 
    593592                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 
     
    605604      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    606605      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points 
    607          DO jj = 2, jpjm1 
    608             DO ji = fs_2, fs_jpim1   ! vector opt. 
     606         DO jj = k_Jstr, k_Jend 
     607            DO ji = k_Istr, k_Iend 
    609608               zsqen = SQRT( en(ji,jj,jk) ) 
    610609               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    611                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    612                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     610               p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     611               p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    613612               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    614613            END DO 
     
    619618      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    620619         DO jk = 2, jpkm1 
    621             DO jj = 2, jpjm1 
    622                DO ji = fs_2, fs_jpim1   ! vector opt. 
    623                   avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     620            DO jj = k_Jstr, k_Jend 
     621               DO ji = k_Istr, k_Iend 
     622                  p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    624623              END DO 
    625624            END DO 
    626625         END DO 
    627626      ENDIF 
    628 !!gm I'm sure this is needed to compute the shear term at next time-step 
    629       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    630       CALL lbc_lnk( avt, 'W', 1. )      ! Lateral boundary conditions on avt  (sign unchanged) 
    631627 
    632628      IF(ln_ctl) THEN 
     
    634630         CALL prt_ctl( tab3d_1=avm, clinfo1=' tke  - m: ', ovlap=1, kdim=jpk ) 
    635631      ENDIF 
    636       ! 
    637       CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
    638632      ! 
    639633      IF( nn_timing == 1 )  CALL timing_stop('tke_avn') 
     
    737731      !                               !* set vertical eddy coef. to the background value 
    738732      DO jk = 1, jpk 
    739          avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    740          avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     733         avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 
     734         avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 
    741735      END DO 
    742736      dissl(:,:,:) = 1.e-12_wp 
     
    748742 
    749743   SUBROUTINE tke_rst( kt, cdrw ) 
    750      !!--------------------------------------------------------------------- 
    751      !!                   ***  ROUTINE tke_rst  *** 
    752      !!                      
    753      !! ** Purpose :   Read or write TKE file (en) in restart file 
    754      !! 
    755      !! ** Method  :   use of IOM library 
    756      !!                if the restart does not contain TKE, en is either  
    757      !!                set to rn_emin or recomputed  
    758      !!---------------------------------------------------------------------- 
    759      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    760      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    761      ! 
    762      INTEGER ::   jit, jk              ! dummy loop indices 
    763      INTEGER ::   id1, id2, id3, id4   ! local integers 
    764      !!---------------------------------------------------------------------- 
    765      ! 
    766      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    767         !                                   ! --------------- 
    768         IF( ln_rstart ) THEN                   !* Read the restart file 
    769            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. ) 
    770            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
    771            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
    772            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 
    773            ! 
    774            IF( id1 > 0 ) THEN                       ! 'en' exists 
    775               CALL iom_get( numror, jpdom_autoglo, 'en', en ) 
    776               IF( MIN( id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist 
    777                  CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
    778                  CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
    779                  CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 
    780               ELSE                                                 ! one at least array is missing 
    781                  CALL tke_avn                                          ! compute avt, avm, and dissl (approximation) 
    782               ENDIF 
    783            ELSE                                     ! No TKE array found: initialisation 
    784               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 
    785               en (:,:,:) = rn_emin * wmask(:,:,:) 
    786               CALL tke_avn                               ! recompute avt, avm, and dissl (approximation) 
    787               avt_k(:,:,:) = avt(:,:,:) 
    788               avm_k(:,:,:) = avm(:,:,:) 
    789               ! 
    790               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    791               avt_k(:,:,:) = avt(:,:,:) 
    792               avm_k(:,:,:) = avm(:,:,:) 
    793            ENDIF 
    794         ELSE                                   !* Start from rest 
    795            en(:,:,:) = rn_emin * tmask(:,:,:) 
    796            DO jk = 1, jpk                           ! set the Kz to the background value 
    797               avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    798               avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    799            END DO 
    800         ENDIF 
    801         ! 
    802      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    803         !                                   ! ------------------- 
    804         IF(lwp) WRITE(numout,*) '---- tke-rst ----' 
    805         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
    806         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
    807         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
    808         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 
    809         ! 
    810      ENDIF 
    811      ! 
     744      !!--------------------------------------------------------------------- 
     745      !!                   ***  ROUTINE tke_rst  *** 
     746      !!                      
     747      !! ** Purpose :   Read or write TKE file (en) in restart file 
     748      !! 
     749      !! ** Method  :   use of IOM library 
     750      !!                if the restart does not contain TKE, en is either  
     751      !!                set to rn_emin or recomputed  
     752      !!---------------------------------------------------------------------- 
     753      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     754      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     755      ! 
     756      INTEGER ::   jit, jk              ! dummy loop indices 
     757      INTEGER ::   id1, id2, id3, id4   ! local integers 
     758      !!---------------------------------------------------------------------- 
     759      ! 
     760      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     761         !                                   ! --------------- 
     762         IF( ln_rstart ) THEN                   !* Read the restart file 
     763            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. ) 
     764            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
     765            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
     766            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 
     767            ! 
     768            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist 
     769               CALL iom_get( numror, jpdom_autoglo, 'en', en ) 
     770               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
     771               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
     772               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 
     773            ELSE                                          ! start TKE from rest 
     774               IF(lwp) WRITE(numout,*) '   ==>>   previous run without TKE scheme, set en to background values' 
     775               en(:,:,:) = rn_emin * wmask(:,:,:) 
     776               ! avt_k, avm_k already set to the background value in zdf_phy_init 
     777            ENDIF 
     778         ELSE                                   !* Start from rest 
     779            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set en to the background value' 
     780            en(:,:,:) = rn_emin * wmask(:,:,:) 
     781            ! avt_k, avm_k already set to the background value in zdf_phy_init 
     782         ENDIF 
     783         ! 
     784      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     785         !                                   ! ------------------- 
     786         IF(lwp) WRITE(numout,*) '---- tke-rst ----' 
     787         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
     788         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
     789         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
     790         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 
     791         ! 
     792      ENDIF 
     793      ! 
    812794   END SUBROUTINE tke_rst 
    813795 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7954 r8055  
    3636   !!---------------------------------------------------------------------- 
    3737   USE step_oce         ! time stepping definition modules 
     38!!gm to be removed when removing avmu, avmv 
     39   USE zdf_oce        ! ocean vertical physics variables 
     40!!gm 
    3841   ! 
    3942   USE iom              ! xIOs server 
     
    102105      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp 
    103106 
     107CALL FLUSH( numout) 
    104108      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    105109      ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice) 
     
    109113      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    110114                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice) 
     115CALL FLUSH( numout) 
    111116 
    112117      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    124129                         CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 
    125130                         CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency 
     131CALL FLUSH( numout) 
    126132 
    127133      !  VERTICAL PHYSICS 
    128134                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD) 
     135CALL FLUSH( numout) 
     136 
     137!!gm  ===>>>   TO BE REMOVED   (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri) 
     138      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
     139         DO jj = 2, jpjm1 
     140            DO ji = 2, jpim1 
     141               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     142               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
     143            END DO 
     144         END DO 
     145      END DO 
     146      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     147!!gm end 
    129148 
    130149      !  LATERAL  PHYSICS 
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/sette.sh

    r7931 r8055  
    10981098    export TEST_NAME="SHORT" 
    10991099    cd ${CONFIG_DIR0} 
    1100     . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top" 
     1100    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 
    11011101    cd ${SETTE_DIR} 
    11021102    . ./param.cfg 
     
    11451145    export TEST_NAME="SHORT_NOZOOM" 
    11461146    cd ${CONFIG_DIR0} 
    1147     . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top" 
     1147    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 
    11481148    cd ${SETTE_DIR} 
    11491149    . ./param.cfg 
     
    11821182    export TEST_NAME="SHORT_NOAGRIF" 
    11831183    cd ${CONFIG_DIR0} 
    1184     . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_zdftmx key_top" 
     1184    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_top" 
    11851185    cd ${SETTE_DIR} 
    11861186    . ./param.cfg 
     
    12201220    export TEST_NAME="LONG" 
    12211221    cd ${CONFIG_DIR0} 
    1222     . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top" 
     1222    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 
    12231223    cd ${SETTE_DIR} 
    12241224    . ./param.cfg 
     
    13191319    export TEST_NAME="REPRO_4_4" 
    13201320    cd ${CONFIG_DIR0} 
    1321     . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top" 
     1321    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 
    13221322    cd ${SETTE_DIR} 
    13231323    . ./param.cfg 
Note: See TracChangeset for help on using the changeset viewer.