Changeset 7991


Ignore:
Timestamp:
2017-05-01T16:21:42+02:00 (3 years ago)
Author:
gm
Message:

#1880 (HPC-09) - step-5: OPA/ZDF shear production term computed outside closure schemes

Location:
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r7990 r7991  
    99   USE dom_oce         ! ocean space and time domain 
    1010   USE zdf_oce         ! ocean vertical physics     
    11    USE zdfgls   , ONLY : hmxn 
     11   USE zdfgls   , ONLY : hmxl_n 
    1212   USE in_out_manager  ! I/O units 
    1313   USE iom             ! I/0 library 
     
    2323 
    2424   ! variables for calculating 25-hourly means 
    25    INTEGER , SAVE ::   cnt_25h     ! Counter for 25 hour means 
     25   INTEGER , SAVE ::   cnt_25h           ! Counter for 25 hour means 
     26   REAL(wp), SAVE ::   r1_25 = 0.04_wp   ! =1/25  
    2627   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h 
    2728   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h  
     
    9495      ! ------------------------- ! 
    9596      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)  
    96       tn_25h(:,:,:) = tsb(:,:,:,jp_tem) 
    97       sn_25h(:,:,:) = tsb(:,:,:,jp_sal) 
    98       sshn_25h(:,:) = sshb(:,:) 
    99       un_25h(:,:,:) = ub(:,:,:) 
    100       vn_25h(:,:,:) = vb(:,:,:) 
    101       wn_25h(:,:,:) = wn(:,:,:) 
    102       avt_25h(:,:,:) = avt(:,:,:) 
    103       avm_25h(:,:,:) = avm(:,:,:) 
     97      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem) 
     98      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal) 
     99      sshn_25h(:,:)   = sshb(:,:) 
     100      un_25h  (:,:,:) = ub  (:,:,:) 
     101      vn_25h  (:,:,:) = vb  (:,:,:) 
     102      wn_25h  (:,:,:) = wn  (:,:,:) 
     103      avt_25h (:,:,:) = avt (:,:,:) 
     104      avm_25h (:,:,:) = avm (:,:,:) 
    104105      IF( ln_zdftke ) THEN 
    105106         en_25h(:,:,:) = en(:,:,:) 
    106107      ENDIF 
    107108      IF( ln_zdfgls ) THEN 
    108          en_25h(:,:,:) = en(:,:,:) 
    109          rmxln_25h(:,:,:) = hmxn(:,:,:) 
     109         en_25h   (:,:,:) = en    (:,:,:) 
     110         rmxln_25h(:,:,:) = hmxl_n(:,:,:) 
    110111      ENDIF 
    111112#if defined key_lim3 || defined key_lim2 
     
    156157         ENDIF 
    157158 
    158          tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem) 
    159          sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal) 
    160          sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:) 
    161          un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:) 
    162          vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:) 
    163          wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:) 
    164          avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:) 
    165          avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:) 
    166          IF( ln_zdftke ) THEN 
    167             en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
    168          ENDIF 
    169          IF( ln_zdfgls ) THEN 
    170             en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
    171             rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + hmxn(:,:,:) 
     159         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem) 
     160         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal) 
     161         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:) 
     162         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:) 
     163         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:) 
     164         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:) 
     165         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:) 
     166         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:) 
     167         IF( ln_zdftke ) THEN 
     168            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:) 
     169         ENDIF 
     170         IF( ln_zdfgls ) THEN 
     171            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:) 
     172            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:) 
    172173         ENDIF 
    173174         cnt_25h = cnt_25h + 1 
     
    187188         ENDIF 
    188189         ! 
    189          tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp 
    190          sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp 
    191          sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp 
    192          un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp 
    193          vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp 
    194          wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp 
    195          avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp 
    196          avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp 
    197          IF( ln_zdftke ) THEN 
    198             en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp 
    199          ENDIF 
    200          IF( ln_zdfgls ) THEN 
    201             en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp 
    202             rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp 
     190         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25 
     191         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25 
     192         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25 
     193         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25 
     194         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25 
     195         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25 
     196         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25 
     197         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25 
     198         IF( ln_zdftke ) THEN 
     199            en_25h(:,:,:) = en_25h(:,:,:) * r1_25 
     200         ENDIF 
     201         IF( ln_zdfgls ) THEN 
     202            en_25h   (:,:,:) = en_25h   (:,:,:) * r1_25 
     203            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) * r1_25 
    203204         ENDIF 
    204205         ! 
     
    236237         ! 
    237238         ! After the write reset the values to cnt=1 and sum values equal current value  
    238          tn_25h(:,:,:) = tsn(:,:,:,jp_tem) 
    239          sn_25h(:,:,:) = tsn(:,:,:,jp_sal) 
    240          sshn_25h(:,:) = sshn (:,:) 
    241          un_25h(:,:,:) = un(:,:,:) 
    242          vn_25h(:,:,:) = vn(:,:,:) 
    243          wn_25h(:,:,:) = wn(:,:,:) 
    244          avt_25h(:,:,:) = avt(:,:,:) 
    245          avm_25h(:,:,:) = avm(:,:,:) 
     239         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem) 
     240         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal) 
     241         sshn_25h(:,:)   = sshn(:,:) 
     242         un_25h  (:,:,:) = un  (:,:,:) 
     243         vn_25h  (:,:,:) = vn  (:,:,:) 
     244         wn_25h  (:,:,:) = wn  (:,:,:) 
     245         avt_25h (:,:,:) = avt (:,:,:) 
     246         avm_25h (:,:,:) = avm (:,:,:) 
    246247         IF( ln_zdftke ) THEN 
    247248            en_25h(:,:,:) = en(:,:,:) 
    248249         ENDIF 
    249250         IF( ln_zdfgls ) THEN 
    250             en_25h(:,:,:) = en(:,:,:) 
    251             rmxln_25h(:,:,:) = hmxn(:,:,:) 
     251            en_25h   (:,:,:) = en    (:,:,:) 
     252            rmxln_25h(:,:,:) = hmxl_n(:,:,:) 
    252253         ENDIF 
    253254         cnt_25h = 1 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7990 r7991  
    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 
    2122   USE zdfbfr         ! bottom friction (only for rn_bfrz0) 
    2223   USE sbc_oce        ! surface boundary condition: ocean 
     
    4243 
    4344   ! 
    44    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxn    !: now mixing length 
     45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length 
    4546   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function 
    4647   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
     
    114115      !!                ***  FUNCTION zdf_gls_alloc  *** 
    115116      !!---------------------------------------------------------------------- 
    116       ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    117          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
     117      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustars2(jpi,jpj) ,       & 
     118         &      zwall (jpi,jpj,jpk) , ustarb2(jpi,jpj) ,  STAT= zdf_gls_alloc ) 
    118119         ! 
    119120      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    141142      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves) 
    142143      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before 
    143       REAL(wp), POINTER, DIMENSION(:,:,:) ::   mxlb        ! mixing length at time before 
     144      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hmxl_b      ! mixing length at time before 
    144145      REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear 
    145146      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate 
     
    149150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal 
    150151      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal 
    151       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3du, z3dv  ! u- and v-shears 
     152      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zstt, zstm  ! stability function on tracer and momentum 
    152153      !!-------------------------------------------------------------------- 
    153154      ! 
     
    155156      ! 
    156157      CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro ) 
    157       CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    158       CALL wrk_alloc( jpi,jpj,jpk,   z3du, z3dv ) 
     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 ) 
    159160 
    160161      ! Preliminary computing 
     
    198199      END SELECT 
    199200 
    200       DO jk = 2, jpkm1              !==  Compute shear and dissipation rate  ==! 
     201      !                             !==  Compute shear and dissipation rate  ==! 
     202      CALL zdf_sh2( shear ) 
     203 
     204 
     205      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==! 
    201206         DO jj = 1, jpjm1 
    202207            DO ji = 1, jpim1 
    203                z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   & 
    204                   &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   & 
    205                   &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    206                z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   & 
    207                   &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   & 
    208                   &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
    209                   ! 
    210                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / hmxn(ji,jj,jk) 
     208               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
    211209            END DO 
    212210         END DO 
     
    214212 
    215213      ! Save tke at before time step 
    216       eb  (:,:,:) = en  (:,:,:) 
    217       mxlb(:,:,:) = hmxn(:,:,:) 
     214      eb    (:,:,:) = en    (:,:,:) 
     215      hmxl_b(:,:,:) = hmxl_n(:,:,:) 
    218216 
    219217      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
     
    221219            DO jj = 2, jpjm1  
    222220               DO ji = fs_2, fs_jpim1   ! vector opt. 
    223                   zup   = hmxn(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     221                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 
    224222                  zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 
    225223                  zcoef = ( zup / MAX( zdown, rsmall ) ) 
     
    247245            DO ji = 2, jpim1 
    248246               ! 
    249                ! shear prod. at w-point weightened by mask 
    250                shear(ji,jj,jk) =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    251                   &             + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
    252                ! 
    253                ! stratif. destruction 
    254                buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk) 
    255                ! 
    256                ! shear prod. - stratif. destruction 
    257                diss = eps(ji,jj,jk) 
     247               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction 
     248               ! 
     249               diss = eps(ji,jj,jk)                         ! dissipation 
    258250               ! 
    259251               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     
    344336      ! 
    345337      CASE ( 0 )             ! Dirichlet  
    346          !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin 
     338         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 
    347339         !                      ! Balance between the production and the dissipation terms 
    348340         DO jj = 2, jpjm1 
     
    425417            DO jj = 2, jpjm1 
    426418               DO ji = fs_2, fs_jpim1   ! vector opt. 
    427                   psi(ji,jj,jk)  = eb(ji,jj,jk) * mxlb(ji,jj,jk) 
     419                  psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 
    428420               END DO 
    429421            END DO 
     
    443435            DO jj = 2, jpjm1 
    444436               DO ji = fs_2, fs_jpim1   ! vector opt. 
    445                   psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) ) 
     437                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 
    446438               END DO 
    447439            END DO 
     
    452444            DO jj = 2, jpjm1 
    453445               DO ji = fs_2, fs_jpim1   ! vector opt. 
    454                   psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn  
     446                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn  
    455447               END DO 
    456448            END DO 
     
    563555      ! 
    564556      CASE ( 0 )             ! Dirichlet  
    565          !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0 
     557         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0 
    566558         !                      ! Balance between the production and the dissipation terms 
    567559         DO jj = 2, jpjm1 
     
    686678      ! Limit dissipation rate under stable stratification 
    687679      ! -------------------------------------------------- 
    688       DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time 
     680      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 
    689681         DO jj = 2, jpjm1 
    690682            DO ji = fs_2, fs_jpim1    ! vector opt. 
    691683               ! limitation 
    692                eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    693                hmxn(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
     684               eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
     685               hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    694686               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
    695687               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    696                IF (ln_length_lim) hmxn(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) ) 
     688               IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 
    697689            END DO 
    698690         END DO 
     
    710702               DO ji = fs_2, fs_jpim1   ! vector opt. 
    711703                  ! zcof =  l²/q² 
    712                   zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
     704                  zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 
    713705                  ! Gh = -N²l²/q² 
    714706                  gh = - rn2(ji,jj,jk) * zcof 
     
    719711                  sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 
    720712                  ! 
    721                   ! Store stability function in z3du and z3dv 
    722                   z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    723                   z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     713                  ! Store stability function in zstt and zstm 
     714                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     715                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    724716               END DO 
    725717            END DO 
     
    731723               DO ji = fs_2, fs_jpim1   ! vector opt. 
    732724                  ! zcof =  l²/q² 
    733                   zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
     725                  zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 
    734726                  ! Gh = -N²l²/q² 
    735727                  gh = - rn2(ji,jj,jk) * zcof 
     
    747739                  sh = (rs4 - rs5*gh + rs6*gm) / rcff 
    748740                  ! 
    749                   ! Store stability function in z3du and z3dv 
    750                   z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
    751                   z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
     741                  ! Store stability function in zstt and zstm 
     742                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 
     743                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 
    752744               END DO 
    753745            END DO 
     
    759751      ! Lines below are useless if GOTM style Dirichlet conditions are used 
    760752 
    761       z3dv(:,:,1) = z3dv(:,:,2) 
    762  
     753      zstm(:,:,1) = zstm(:,:,2) 
     754 
     755!!gm should be done for ISF (top boundary cond.) 
    763756      DO jj = 2, jpjm1 
    764757         DO ji = fs_2, fs_jpim1   ! vector opt. 
    765             z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj)) 
     758            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 
    766759         END DO 
    767760      END DO 
     
    772765         DO jj = 2, jpjm1 
    773766            DO ji = fs_2, fs_jpim1   ! vector opt. 
    774                zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxn(ji,jj,jk) 
    775                zav           = zsqen * z3du(ji,jj,jk) 
     767               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
     768               zav           = zsqen * zstt(ji,jj,jk) 
    776769               avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
    777                zav           = zsqen * z3dv(ji,jj,jk) 
     770               zav           = zsqen * zstm(ji,jj,jk) 
    778771               avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
    779772            END DO 
     
    795788      ! 
    796789      CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro ) 
    797       CALL wrk_dealloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 
    798       CALL wrk_dealloc( jpi,jpj,jpk,   z3du, z3dv ) 
     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 ) 
    799792      ! 
    800793      IF( nn_timing == 1 )   CALL timing_stop('zdf_gls') 
     
    11001093 
    11011094      !                                !* read or initialize all required files   
    1102       CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxn) 
     1095      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n) 
    11031096      ! 
    11041097      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init') 
     
    11291122         !                                   ! --------------- 
    11301123         IF( ln_rstart ) THEN                   !* Read the restart file 
    1131             id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. ) 
    1132             id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
    1133             id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
    1134             id4 = iom_varid( numror, 'hmxn' , ldstop = .FALSE. ) 
     1124            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. ) 
     1125            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. ) 
     1126            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. ) 
     1127            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. ) 
    11351128            ! 
    11361129            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist 
    1137                CALL iom_get( numror, jpdom_autoglo, 'en'   , en    ) 
    1138                CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
    1139                CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
    1140                CALL iom_get( numror, jpdom_autoglo, 'hmxn' , hmxn ) 
     1130               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     ) 
     1131               CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k ) 
     1132               CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k ) 
     1133               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 
    11411134            ELSE                         
    1142                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop' 
    1143                en   (:,:,:) = rn_emin 
    1144                hmxn (:,:,:) = 0.05_wp 
    1145                avt_k(:,:,:) = avt(:,:,:) 
    1146                avm_k(:,:,:) = avm(:,:,:) 
     1135               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 
     1136               en    (:,:,:) = rn_emin 
     1137               hmxl_n(:,:,:) = 0.05_wp 
     1138               avt_k (:,:,:) = avt(:,:,:) 
     1139               avm_k (:,:,:) = avm(:,:,:) 
    11471140               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO 
    11481141            ENDIF 
    11491142         ELSE                                   !* Start from rest 
    1150             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values' 
    1151             en  (:,:,:) = rn_emin 
    1152             hmxn(:,:,:) = 0.05_wp 
     1143            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 
     1144            en    (:,:,:) = rn_emin 
     1145            hmxl_n(:,:,:) = 0.05_wp 
    11531146            DO jk = 1, jpk 
    11541147               avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 
     
    11601153         !                                   ! ------------------- 
    11611154         IF(lwp) WRITE(numout,*) '---- gls-rst ----' 
    1162          CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )  
    1163          CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
    1164          CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
    1165          CALL iom_rstput( kt, nitrst, numrow, 'hmxn' , hmxn ) 
     1155         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     )  
     1156         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k ) 
     1157         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k ) 
     1158         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n ) 
    11661159         ! 
    11671160      ENDIF 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r7990 r7991  
    2121   USE oce            ! ocean dynamics and tracers variables 
    2222   USE dom_oce        ! ocean space and time domain variables 
    23    USE zdf_oce        ! ocean vertical physics 
     23   USE zdf_oce        ! vertical physics: variables 
     24   USE zdfsh2         ! vertical physics: shear production term of TKE 
    2425   USE phycst         ! physical constants 
    2526   USE sbc_oce,  ONLY :   taum 
     
    104105      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    105106      LOGICAL  ::   ll_av_weight = .TRUE.      ! local logical 
    106       REAL(wp) ::   zsh2, zcfRi, zav, zustar   ! local scalars 
    107       REAL(wp), DIMENSION(jpi,jpj) ::   zdu2, zdv2, zh_ekm   ! 2D workspace 
     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     - 
    108110      !!---------------------------------------------------------------------- 
    109111      ! 
    110112      IF( nn_timing == 1 )   CALL timing_start('zdf_ric') 
    111113      ! 
    112       DO jk = 2, jpkm1        !==  avm and avt = F(Richardson number)  ==! 
    113          ! 
    114       IF( ll_av_weight ) THEN    !== viscosity weighted shear & stratification terms 
    115          ! 
    116          DO jj = 1, jpjm1           !* viscosity weighted shear term 
     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 
    117121            DO ji = 1, jpim1 
    118                zdu2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) ) * wumask(ji,jj,jk)      & 
    119                   &              * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         & 
    120                   &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    121                zdv2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) ) * wumask(ji,jj,jk)      & 
    122                   &              * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         & 
    123                   &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
    124             END DO 
    125          END DO 
    126          DO jj = 2, jpjm1 
    127             DO ji = 2, jpim1        !* Richardson number dependent coefficient 
    128                !                          ! shear of horizontal velocity 
    129                zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    130                   &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    131                !                          ! coefficient = F(richardson number) 
    132                zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avt(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  ) 
     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 
    133125               ! 
    134                !                    !* Vertical eddy viscosity and diffusivity coefficients 
    135                zav = rn_avmri * zcfRi**nn_ric 
     126               !                          ! avm and avt coefficients 
    136127               avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk) 
    137128               avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk) 
    138129            END DO 
    139130         END DO 
    140       ELSE                    !==  classical Richardson number definition  ==! 
    141          DO jj = 1, jpjm1           !* shear term 
    142             DO ji = 1, jpim1 
    143                zdu2(ji,jj) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         & 
    144                   &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    145                zdv2(ji,jj) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         & 
    146                   &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
    147             END DO 
    148          END DO 
    149          DO jj = 2, jpjm1 
    150             DO ji = 2, jpim1        !* Richardson number dependent coefficient 
    151                !                          ! shear of horizontal velocity 
    152                zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    153                   &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    154                !                          ! coefficient = F(richardson number) 
    155                zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  ) 
    156                ! 
    157                !                    !* Vertical eddy viscosity and diffusivity coefficients 
    158                zav = rn_avmri * zcfRi**nn_ric 
    159                avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk) 
    160                avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk) 
    161             END DO 
    162          END DO 
    163       ENDIF 
    164          ! 
    165131      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 !!! ) 
    166135      ! 
    167136      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==! 
     
    184153            END DO 
    185154         END DO 
    186       END IF 
     155      ENDIF 
    187156      ! 
    188157      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric') 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r7990 r7991  
    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 
    4546   USE zdfmxl         ! vertical physics: mixed layer 
    4647   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    8990   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9091   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation 
    91 #if defined key_c1d 
    92    !                                                                        !!** 1D cfg only  **   ('key_c1d') 
    93    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales 
    94    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers 
    95 #endif 
    9692 
    9793   !! * Substitutions 
     
    108104      !!                ***  FUNCTION zdf_tke_alloc  *** 
    109105      !!---------------------------------------------------------------------- 
    110       ALLOCATE(                                                                    & 
    111 #if defined key_c1d 
    112          &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          & 
    113          &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          & 
    114 #endif 
    115          &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
    116          &      apdlr(jpi,jpj,jpk) ,                                           STAT= zdf_tke_alloc      ) 
     106      ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) ,   STAT= zdf_tke_alloc ) 
    117107         ! 
    118108      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
     
    211201      !! --------------------------------------------------------------------- 
    212202      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
    213 !!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar 
    214 !!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar 
    215       REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
    216       REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
    217       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
    218       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
    219       REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
    220       REAL(wp) ::   ztau  , zdif                    !    -         - 
    221       REAL(wp) ::   zus   , zwlc  , zind            !    -         - 
    222       REAL(wp) ::   zzd_up, zzd_lw                  !    -         - 
    223 !!bfr      REAL(wp) ::   zebot                           !    -         - 
     203!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! local integers 
     204!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          !   -      - 
     205!!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           !   -         - 
    224214      INTEGER , DIMENSION(jpi,jpj) ::   imlc 
    225215      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc 
    226       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 
    227       REAL(wp)                            ::   zri  !   local Richardson number 
     216      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw 
     217      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2 
    228218      !!-------------------------------------------------------------------- 
    229219      ! 
     
    231221      ! 
    232222      CALL wrk_alloc( jpi,jpj,       zhlc )  
    233       CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
     223      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw )  
    234224      ! 
    235225      zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     
    328318      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 
    329319      ! 
    330       DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form) 
    331          DO jj = 1, jpjm1 
    332             DO ji = 1, fs_jpim1   ! vector opt. 
    333                z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   & 
    334                   &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   & 
    335                   &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) & 
    336                   &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 
    337                z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   & 
    338                   &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   & 
    339                   &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) & 
    340                   &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 
    341             END DO 
    342          END DO 
    343       END DO 
    344       ! 
    345       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr 
    346          ! Note that zesh2 is also computed in the next loop. 
    347          ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 
     320      !                           !* Shear production at uw- and vw-points (energy conserving form) 
     321      CALL zdf_sh2( zsh2 ) 
     322      ! 
     323      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri ) 
    348324         DO jk = 2, jpkm1 
    349325            DO jj = 2, jpjm1 
    350                DO ji = fs_2, fs_jpim1   ! vector opt. 
    351                   !                                          ! shear prod. at w-point weightened by mask 
    352                   zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    353                      &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    354                   !                                          ! local Richardson number 
    355                   zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 
     326               DO ji = 2, jpim1 
     327                  !                             ! local Richardson number 
     328                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear ) 
     329                  !                             ! inverse of Prandtl number 
    356330                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
    357                    
    358                END DO 
    359             END DO 
    360          END DO 
    361          ! 
     331               END DO 
     332            END DO 
     333         END DO 
    362334      ENDIF 
    363335      !          
     
    373345                  &          /    ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  ) 
    374346               ! 
    375                !                                   ! shear prod. at w-point weightened by mask 
    376                zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    377                   &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    378                ! 
    379347               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    380348               zd_lw(ji,jj,jk) = zzd_lw 
    381                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
     349               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 
    382350               ! 
    383351               !                                   ! right hand side in en 
    384                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    385                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    386                   &                                 * wmask(ji,jj,jk) 
     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 
     354                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
     355                  &                                ) * wmask(ji,jj,jk) 
    387356            END DO 
    388357         END DO 
     
    470439      ! 
    471440      CALL wrk_dealloc( jpi,jpj,       zhlc )  
    472       CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
     441      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw )  
    473442      ! 
    474443      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
     
    631600      END SELECT 
    632601      ! 
    633 # if defined key_c1d 
    634       e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales 
    635       e_mix(:,:,:) = zmxlm(:,:,:) 
    636 # endif 
    637602 
    638603      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    657622               DO ji = fs_2, fs_jpim1   ! vector opt. 
    658623                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    659 # if defined key_c1d 
    660                   e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number 
    661 # endif 
    662624              END DO 
    663625            END DO 
Note: See TracChangeset for help on using the changeset viewer.