Changeset 8093


Ignore:
Timestamp:
2017-05-30T10:13:14+02:00 (3 years ago)
Author:
gm
Message:

#1880 (HPC-09) - step-6: prepare some forthcoming evolutions (ZDF modules mainly)

Location:
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM
Files:
1 added
22 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7990 r8093  
    601601 
    602602!!====================================================================== 
    603 !!                 ***  Bottom boundary condition  *** 
    604 !!====================================================================== 
    605 !!   nambfr        bottom friction 
    606 !!   nambbc        bottom temperature boundary condition 
     603!!                ***  top/Bottom boundary condition  ***             !! 
     604!!====================================================================== 
     605!!   nambfr        bottom friction                                      (default: NONE) 
     606!!   namtfr        top    friction                                      (default: NONE) 
     607!!   nambbc        bottom temperature boundary condition                (default: NO) 
    607608!!   nambbl        bottom boundary layer scheme                         (default: NO) 
    608609!!====================================================================== 
     
    611612&nambfr        !   bottom friction                                      (default: linear) 
    612613!----------------------------------------------------------------------- 
    613    nn_bfr      =    1      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction 
    614                            !                              = 2 : nonlinear friction 
     614   nn_bfr      =    1      !  type of top/bottom drag: free slip (=0), linear drag (=1),  
     615   !                       !  nonlinear drag (=2), nonlinear with logarithmic formulation (=3) 
     616   ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true) 
     617   ln_loglayer = .false.   !  logarithmic formulation (non linear case only) 
     618                           ! 
    615619   rn_bfri1    =    4.e-4  !  bottom drag coefficient (linear case) 
    616620   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T 
     
    619623   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T 
    620624   ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file ) 
    621    rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T) 
     625      rn_bfrien   =    50.    !  local boost factor 
     626   ! 
    622627   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case) 
    623628   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T 
     
    626631   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T 
    627632   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file ) 
    628    rn_tfrien   =   50.     !  local multiplying factor of tfr (ln_tfr2d=T) 
    629  
    630    ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true) 
    631    ln_loglayer = .false.   !  logarithmic formulation (non linear case) 
    632 / 
     633      rn_tfrien   =   50.     !  local boost factor 
     634/ 
     635 
     636!----------------------------------------------------------------------- 
     637&namdrg            !   top/bottom drag coeeficient                      (default: NO selection) 
     638!----------------------------------------------------------------------- 
     639   ln_NONE    = .false.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
     640   ln_lin     = .false.    !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top) 
     641   ln_non_lin = .false.    !  non-linear  drag: Cd = Cd0 |U| 
     642   ln_loglayer= .false.    !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
     643   ! 
     644   ln_drgimp  = .true.     !  implicit top/bottom friction flag 
     645/ 
     646 
     647!----------------------------------------------------------------------- 
     648&namdrg_top        !   TOP friction                                     (ln_isfcav=T) 
     649!----------------------------------------------------------------------- 
     650   rn_Cd0     =  1.e-3     !  drag coefficient [-] 
     651   rn_Uc0     =  0.4       !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     652   rn_Cdmax   =  0.1       !  drag value maximum [-] (logarithmic drag) 
     653   rn_ke0     =  2.5e-3    !  background kinetic energy  [m2/s2] (non-linear cases) 
     654   rn_z0      =  3.0e-3    !  roughness [m] (ln_loglayer=T) 
     655   ln_boost   = .false.    !  =T regional boost of Cd0 ; =F constant 
     656      rn_boost=  50.          !  local boost factor  [-] 
     657/ 
     658 
     659!----------------------------------------------------------------------- 
     660&namdrg_bot        !   BOTTOM friction                                   
     661!----------------------------------------------------------------------- 
     662   rn_Cd0     =  1.e-3    !  drag coefficient [-] 
     663   rn_Uc0     =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0)  
     664   rn_Cdmax   =  0.1      !  drag value maximum [-] (logarithmic drag) 
     665   rn_ke0     =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     666   rn_z0      =  3.e-3    !  roughness [m] (ln_loglayer=T) 
     667   ln_boost   = .false.   !  =T regional boost of Cd0 ; =F constant 
     668      rn_boost=  50.         !  local boost factor  [-] 
     669/ 
     670 
    633671!----------------------------------------------------------------------- 
    634672&nambbc        !   bottom temperature boundary condition                (default: NO) 
     
    935973   rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2] 
    936974   rn_bshear   =   1.e-20  ! background shear (>0) currently a numerical threshold (do not change it) 
     975   nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm) 
    937976   nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom 
    938977                           !                 = 1 bounded by the local vertical scale factor 
    939978                           !                 = 2 first vertical derivative of mixing length bounded by 1 
    940979                           !                 = 3 as =2 with distinct disspipative an mixing length scale 
    941    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm) 
    942980   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F) 
    943981   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value 
     982!!gm   ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE 
    944983   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002) 
    945984   rn_lc       =   0.15    !  coef. associated to Langmuir cells 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90

    r7953 r8093  
    140140      ! Physical and dynamical ocean fields for output or passing to TOP, time-mean fields 
    141141      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE      :: tsn_crs 
    142       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: un_crs, vn_crs, wn_crs, rke_crs 
     142      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: un_crs, vn_crs, wn_crs 
    143143      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: hdivn_crs     
    144144      REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshn_crs     
     
    227227 
    228228 
    229       ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) , vn_crs(jpi_crs,jpj_crs,jpk) , & 
    230          &      wn_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk),& 
    231          &      rke_crs(jpi_crs,jpj_crs,jpk),                                STAT=ierr(11)) 
     229      ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) , vn_crs(jpi_crs,jpj_crs,jpk)   ,     & 
     230         &      wn_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(11)) 
    232231 
    233232     ALLOCATE( sshn_crs(jpi_crs,jpj_crs), emp_crs (jpi_crs,jpj_crs), emp_b_crs(jpi_crs,jpj_crs), & 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90

    r7953 r8093  
    5858      INTEGER  ::   ji, jj, jk        ! dummy loop indices 
    5959      REAL(wp) ::   z2dcrsu, z2dcrsv  ! local scalars 
     60      REAL(wp) ::   zztmp             !   -      - 
    6061      ! 
    6162      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze3t, ze3u, ze3v, ze3w   ! 3D workspace for e3 
    62       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zt, zt_crs 
     63      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zt, zt_crs, z3d 
    6364      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zs, zs_crs   
    6465      !!---------------------------------------------------------------------- 
     
    6970      CALL wrk_alloc( jpi,jpj,jpk,   ze3t, ze3w ) 
    7071      CALL wrk_alloc( jpi,jpj,jpk,   ze3u, ze3v ) 
    71       CALL wrk_alloc( jpi,jpj,jpk,   zt  , zs  ) 
     72      CALL wrk_alloc( jpi,jpj,jpk,   zt  , zs  ,z3d ) 
    7273      ! 
    7374      CALL wrk_alloc( jpi_crs,jpj_crs,jpk,   zt_crs, zs_crs ) 
     
    8687         avs_crs  (:,:,:  ) = 0._wp    ! avt 
    8788         hdivn_crs(:,:,:  ) = 0._wp    ! hdiv 
    88          rke_crs  (:,:,:  ) = 0._wp    ! rke 
    8989         sshn_crs (:,:    ) = 0._wp    ! ssh 
    9090         utau_crs (:,:    ) = 0._wp    ! taux 
     
    158158      CALL iom_put( "voces" , zs_crs )   ! vS 
    159159 
    160       
    161       !  Kinetic energy 
    162       CALL crs_dom_ope( rke, 'VOL', 'T', tmask, rke_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 
    163       CALL iom_put( "eken", rke_crs ) 
    164  
     160      IF( iom_use( "eken") ) THEN     !      kinetic energy 
     161         z3d(:,:,jk) = 0._wp  
     162         DO jk = 1, jpkm1 
     163            DO jj = 2, jpjm1 
     164               DO ji = fs_2, fs_jpim1   ! vector opt. 
     165                  zztmp  = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     166                  z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    & 
     167                     &            un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
     168                     &          + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
     169                     &          + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
     170                     &          + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
     171               END DO 
     172            END DO 
     173         END DO 
     174         CALL lbc_lnk( z3d, 'T', 1. ) 
     175         ! 
     176         CALL crs_dom_ope( z3d, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 
     177         CALL iom_put( "eken", zt_crs ) 
     178      ENDIF 
    165179      !  Horizontal divergence ( following OPA_SRC/DYN/divhor.F90 )  
    166180      DO jk = 1, jpkm1 
     
    175189                   hdivn_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) / crs_volt_wgt(ji,jj,jk)  
    176190               ENDIF 
    177             ENDDO 
    178          ENDDO 
    179       ENDDO 
     191            END DO 
     192         END DO 
     193      END DO 
    180194      CALL crs_lbc_lnk( hdivn_crs, 'T', 1.0 ) 
    181195      ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r7953 r8093  
    2929   USE dynadv, ONLY: ln_dynadv_vec 
    3030   USE zdf_oce         ! ocean vertical physics 
     31   USE zdfdrg          ! ocean vertical physics: top/bottom friction 
    3132   USE ldftra          ! lateral physics: eddy diffusivity coef. 
    3233   USE ldfdyn          ! lateral physics: eddy viscosity   coef. 
     
    119120      !! ** Method  :  use iom_put 
    120121      !!---------------------------------------------------------------------- 
    121       !! 
    122122      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    123123      !! 
    124       INTEGER                      ::   ji, jj, jk              ! dummy loop indices 
    125       INTEGER                      ::   jkbot                   ! 
    126       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !  
    127       !! 
    128       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    129       REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
     124      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     125      INTEGER ::   ikbot            ! local integer 
     126      REAL(wp)::   zztmp , zztmpx   ! local scalar 
     127      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     128      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
    130130      !!---------------------------------------------------------------------- 
    131131      !  
    132132      IF( nn_timing == 1 )   CALL timing_start('dia_wri') 
    133133      !  
    134       CALL wrk_alloc( jpi , jpj      , z2d ) 
    135       CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
    136       ! 
    137134      ! Output the initial state and forcings 
    138135      IF( ninist == 1 ) THEN                        
     
    162159         DO jj = 1, jpj 
    163160            DO ji = 1, jpi 
    164                jkbot = mbkt(ji,jj) 
    165                z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem) 
     161               ikbot = mbkt(ji,jj) 
     162               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 
    166163            END DO 
    167164         END DO 
     
    174171         DO jj = 1, jpj 
    175172            DO ji = 1, jpi 
    176                jkbot = mbkt(ji,jj) 
    177                z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal) 
     173               ikbot = mbkt(ji,jj) 
     174               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 
    178175            END DO 
    179176         END DO 
     
    182179 
    183180      IF ( iom_use("taubot") ) THEN                ! bottom stress 
     181         zztmp = rau0 * 0.25 
    184182         z2d(:,:) = 0._wp 
    185183         DO jj = 2, jpjm1 
    186184            DO ji = fs_2, fs_jpim1   ! vector opt. 
     185!!gm old 
     186!!gm BUG  missing x 0.5 
    187187               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  & 
    188188                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )       
     
    190190                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )  
    191191               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)  
     192!!gm 
     193               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   & 
     194                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   & 
     195                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   & 
     196                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2 
     197               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
     198!!gm new end 
    192199               ! 
    193200            ENDDO 
     
    202209         DO jj = 1, jpj 
    203210            DO ji = 1, jpi 
    204                jkbot = mbku(ji,jj) 
    205                z2d(ji,jj) = un(ji,jj,jkbot) 
     211               ikbot = mbku(ji,jj) 
     212               z2d(ji,jj) = un(ji,jj,ikbot) 
    206213            END DO 
    207214         END DO 
     
    214221         DO jj = 1, jpj 
    215222            DO ji = 1, jpi 
    216                jkbot = mbkv(ji,jj) 
    217                z2d(ji,jj) = vn(ji,jj,jkbot) 
     223               ikbot = mbkv(ji,jj) 
     224               z2d(ji,jj) = vn(ji,jj,ikbot) 
    218225            END DO 
    219226         END DO 
     
    281288      ! 
    282289      IF ( iom_use("eken") ) THEN 
    283          rke(:,:,jk) = 0._wp                               !      kinetic energy  
     290         z3d(:,:,jk) = 0._wp                               !      kinetic energy  
    284291         DO jk = 1, jpkm1 
    285292            DO jj = 2, jpjm1 
    286293               DO ji = fs_2, fs_jpim1   ! vector opt. 
    287                   zztmp   = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
    288                   zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    & 
    289                      &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) )  & 
    290                      &          *  zztmp  
    291                   ! 
    292                   zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)    & 
    293                      &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  & 
    294                      &          *  zztmp  
    295                   ! 
    296                   rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy ) 
    297                   ! 
    298                ENDDO 
    299             ENDDO 
    300          ENDDO 
    301          CALL lbc_lnk( rke, 'T', 1. ) 
    302          CALL iom_put( "eken", rke )            
     294                  zztmp  = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     295                  z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    & 
     296                     &            un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
     297                     &          + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
     298                     &          + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
     299                     &          + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
     300               END DO 
     301            END DO 
     302         END DO 
     303         CALL lbc_lnk( z3d, 'T', 1. ) 
     304         CALL iom_put( "eken", z3d )            
    303305      ENDIF 
    304306      ! 
     
    407409      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2) 
    408410      ! 
    409       CALL wrk_dealloc( jpi , jpj      , z2d ) 
    410       CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
    411       ! 
    412       ! If we want tmb values  
    413  
    414       IF (ln_diatmb) THEN 
     411 
     412      IF (ln_diatmb) THEN      ! If we want tmb values  
    415413         CALL dia_tmb  
    416414      ENDIF  
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r7753 r8093  
    1515   USE dom_oce        ! ocean space and time domain variables  
    1616   USE zdf_oce        ! ocean vertical physics variables 
     17!!gm new 
     18   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     19!!gm old 
    1720   USE zdfbfr         ! ocean bottom friction variables 
     21!!gm 
    1822   USE trd_oce        ! trends: ocean variables 
    1923   USE trddyn         ! trend manager: dynamics 
     24   ! 
    2025   USE in_out_manager ! I/O manager 
    2126   USE prtctl         ! Print control 
     
    5055      INTEGER  ::   ikbu, ikbv   ! local integers 
    5156      REAL(wp) ::   zm1_2dt      ! local scalar 
    52       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     57      REAL(wp) ::   zCdu, zCdv   !   -      - 
     58      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv 
    5359      !!--------------------------------------------------------------------- 
    5460      ! 
     
    6167 
    6268!!gm bug : time step is only rdt (not 2 rdt if euler start !) 
    63         zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     69         zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
    6470 
    65         IF( l_trddyn ) THEN      ! trends: store the input trends 
    66            CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
    67            ztrdu(:,:,:) = ua(:,:,:) 
    68            ztrdv(:,:,:) = va(:,:,:) 
    69         ENDIF 
     71         IF( l_trddyn ) THEN      ! trends: store the input trends 
     72            ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
     73            ztrdu(:,:,:) = ua(:,:,:) 
     74            ztrdv(:,:,:) = va(:,:,:) 
     75         ENDIF 
    7076 
    7177 
    72         DO jj = 2, jpjm1 
    73            DO ji = 2, jpim1 
    74               ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    75               ikbv = mbkv(ji,jj) 
    76               ! 
    77               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    78               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    79               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
    80            END DO 
    81         END DO 
    82         ! 
    83         IF( ln_isfcav ) THEN        ! ocean cavities 
    84            DO jj = 2, jpjm1 
    85               DO ji = 2, jpim1 
    86                  ! (ISF) stability criteria for top friction 
    87                  ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
    88                  ikbv = mikv(ji,jj) 
    89                  ! 
    90                  ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     78         DO jj = 2, jpjm1 
     79            DO ji = 2, jpim1 
     80               ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels 
     81               ikbv = mbkv(ji,jj) 
     82               ! 
     83               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     84!!gm old 
     85               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
     86               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
     87!!gm new 
     88!               zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu) 
     89!               zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv) 
     90!               ! 
     91!               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu) 
     92!               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv) 
     93!!gm 
     94            END DO 
     95         END DO 
     96         ! 
     97         IF( ln_isfcav ) THEN        ! ocean cavities 
     98            DO jj = 2, jpjm1 
     99               DO ji = 2, jpim1 
     100                  ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
     101                  ikbv = mikv(ji,jj) 
     102                  ! 
     103                  ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     104!!gm old 
    91105                 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
    92106                    &             * (1.-umask(ji,jj,1)) 
    93107                 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
    94108                    &             * (1.-vmask(ji,jj,1)) 
    95                  ! (ISF) 
     109!!gm new 
     110!                  zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked 
     111!                  zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv) 
     112!                  ! 
     113!                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu) 
     114!                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv) 
     115!!gm 
    96116              END DO 
    97117           END DO 
    98         END IF 
     118         END IF 
    99119        ! 
    100120        IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics 
     
    102122           ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    103123           CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 
    104            CALL wrk_dealloc( jpi,jpj,jpk,  ztrdu, ztrdv ) 
     124           DEALLOCATE( ztrdu, ztrdv ) 
    105125        ENDIF 
    106126        !                                          ! print mean trends (used for debugging) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r7761 r8093  
    14571457      !!                 ***  ROUTINE interp1  *** 
    14581458      !! 
    1459       !! ** Purpose :   Calculate the first order of deriavtive of 
     1459      !! ** Purpose :   Calculate the first order of derivative of 
    14601460      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3 
    14611461      !! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7831 r8093  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    2929   USE zdf_oce         ! Bottom friction coefts 
     30!!gm new 
     31   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     32!!gm 
    3033   USE sbcisf          ! ice shelf variable (fwfisf) 
    3134   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     
    146149      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    147150      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
    148       REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      - 
     151      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2     !   -      - 
    149152      REAL(wp) ::   zu_spg, zv_spg              !   -      - 
    150       REAL(wp) ::   zhura, zhvra          !   -      - 
    151       REAL(wp) ::   za0, za1, za2, za3    !   -      - 
    152       ! 
    153       REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 
    154       REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
    155       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 
    156       REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
    157       REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    158       REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    159       REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
     153      REAL(wp) ::   zhura, zhvra                !   -      - 
     154      REAL(wp) ::   za0, za1, za2, za3          !   -      - 
     155      REAL(wp) ::   zztmp                       !   -      - 
     156      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e 
     157      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
     158      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zhdiv 
     159      REAL(wp), DIMENSION(jpi,jpj) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
     160      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zsshv_a 
     161      REAL(wp), DIMENSION(jpi,jpj) :: zhf 
     162      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
     163      ! 
     164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy   ! Wetting/Dying gravity filter coef. 
    160165      !!---------------------------------------------------------------------- 
    161166      ! 
    162       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
    163       ! 
    164       !                                         !* Allocate temporary arrays 
    165       CALL wrk_alloc( jpi,jpj,   zsshp2_e, zhdiv ) 
    166       CALL wrk_alloc( jpi,jpj,   zu_trd, zv_trd) 
    167       CALL wrk_alloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc) 
    168       CALL wrk_alloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
    169       CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    170       CALL wrk_alloc( jpi,jpj,   zhf ) 
    171       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
     167      IF( nn_timing == 1 )   CALL timing_start('dyn_spg_ts') 
     168      ! 
     169      IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 
    172170      ! 
    173171      zmdi=1.e+20                               !  missing data indicator for masking 
     
    211209         ! 
    212210      ENDIF 
     211      ! 
     212!!gm old/new 
     213      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities) 
     214         zCdU_u(:,:) = bfrua(:,:) + tfrua(:,:) 
     215         zCdU_v(:,:) = bfrva(:,:) + tfrva(:,:) 
     216      ELSE                    ! bottom friction only 
     217         zCdU_u(:,:) = bfrua(:,:) 
     218         zCdU_v(:,:) = bfrva(:,:) 
     219      ENDIF 
     220!!gm new 
     221!      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities) 
     222!         DO jj = 2, jpjm1 
     223!            DO ji = fs_2, fs_jpim1   ! vector opt. 
     224!               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 
     225!               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 
     226!            END DO 
     227!         END DO 
     228!      ELSE                    ! bottom friction only 
     229!         DO jj = 2, jpjm1 
     230!            DO ji = fs_2, fs_jpim1   ! vector opt. 
     231!               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 
     232!               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 
     233!            END DO 
     234!         END DO 
     235!      ENDIF 
     236!!gm       
    213237      ! 
    214238      ! Set arrays to remove/compute coriolis trend. 
     
    415439                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    416440                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
    417                               &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     441                     &             / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    418442                ELSE 
    419443                  zcpx(ji,jj) = 0._wp 
     
    491515      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    492516      IF( ln_wd ) THEN 
    493         zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
    494         zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     517         zztmp = - 1._wp / rdtbt 
     518         DO jj = 2, jpjm1                           
     519            DO ji = fs_2, fs_jpim1   ! vector opt. 
     520!!gm old 
     521               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * bfrua(ji,jj) , zztmp ) * zwx(ji,jj) 
     522               zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * bfrva(ji,jj) , zztmp ) * zwy(ji,jj) 
     523!!gm new 
     524!               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) 
     525!               zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj) 
     526!!gm 
     527            END DO 
     528         END DO 
    495529      ELSE 
    496         zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
    497         zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 
     530         DO jj = 2, jpjm1                           
     531            DO ji = fs_2, fs_jpim1   ! vector opt. 
     532!!gm old 
     533               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * bfrua(ji,jj) * zwx(ji,jj) 
     534               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * bfrva(ji,jj) * zwy(ji,jj) 
     535!!gm new 
     536!               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 
     537!               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 
     538!!gm 
     539            END DO 
     540         END DO 
    498541      END IF 
    499542      ! 
     
    520563      ! 
    521564      ! Note that the "unclipped" top friction parameter is used even with explicit drag 
    522       zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 
    523       zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 
     565      DO jj = 2, jpjm1               
     566         DO ji = fs_2, fs_jpim1   ! vector opt. 
     567            zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * tfrua(ji,jj) * zwx(ji,jj) 
     568            zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * tfrva(ji,jj) * zwy(ji,jj) 
     569         END DO 
     570      END DO 
    524571      !        
    525572      IF (ln_bt_fw) THEN                        ! Add wind forcing 
    526          zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 
    527          zv_frc(:,:) =  zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 
     573         DO jj = 2, jpjm1 
     574            DO ji = fs_2, fs_jpim1   ! vector opt. 
     575               zu_frc(ji,jj) =  zu_frc(ji,jj) + zraur * utau(ji,jj) * r1_hu_n(ji,jj) 
     576               zv_frc(ji,jj) =  zv_frc(ji,jj) + zraur * vtau(ji,jj) * r1_hv_n(ji,jj) 
     577            END DO 
     578         END DO 
    528579      ELSE 
    529          zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 
    530          zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 
     580         DO jj = 2, jpjm1 
     581            DO ji = fs_2, fs_jpim1   ! vector opt. 
     582               zu_frc(ji,jj) =  zu_frc(ji,jj) + zraur * z1_2 * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 
     583               zv_frc(ji,jj) =  zv_frc(ji,jj) + zraur * z1_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 
     584            END DO 
     585         END DO 
    531586      ENDIF   
    532587      ! 
     
    885940         ENDIF 
    886941         ! 
    887          ! Add bottom stresses: 
    888          zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
    889          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
    890          ! 
    891          ! Add top stresses: 
    892          zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 
    893          zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
     942         DO jj = 2, jpjm1 
     943            DO ji = fs_2, fs_jpim1   ! vector opt. 
     944               ! Add top/bottom stresses: 
     945!!gm old/new 
     946               zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 
     947               zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 
     948!!gm 
     949            END DO 
     950         END DO 
    894951         ! 
    895952         ! Surface pressure trend: 
     
    10911148      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' ) 
    10921149      ! 
    1093       CALL wrk_dealloc( jpi,jpj,   zsshp2_e, zhdiv ) 
    1094       CALL wrk_dealloc( jpi,jpj,   zu_trd, zv_trd ) 
    1095       CALL wrk_dealloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc ) 
    1096       CALL wrk_dealloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
    1097       CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
    1098       CALL wrk_dealloc( jpi,jpj,   zhf ) 
    1099       IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
     1150      IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 
    11001151      ! 
    11011152      IF ( ln_diatmb ) THEN 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r7990 r8093  
    2222   USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form 
    2323   USE zdf_oce        ! ocean vertical physics 
     24!!gm new 
     25   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     26!!gm old 
    2427   USE zdfbfr         ! Bottom friction setup 
     28!!gm 
    2529   ! 
    2630   USE in_out_manager ! I/O manager 
     
    116120               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
    117121               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     122!!gm old 
    118123               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1) 
    119124               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1) 
     125!!gm new 
     126!               avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 
     127!               avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 
     128!!gm 
    120129            END DO 
    121130         END DO 
     
    125134                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
    126135                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     136!!gm old 
    127137                  IF( ikbu >= 2 )   avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu) 
    128138                  IF( ikbv >= 2 )   avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv) 
     139!!gm new 
     140                  ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 
     141!                  avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 
     142!                  avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 
     143!!gm 
    129144               END DO 
    130145            END DO 
     
    148163               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
    149164               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
     165!!gm old 
    150166               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    151167               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
     168!!gm new 
     169!               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     170!               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     171!!gm 
    152172            END DO 
    153173         END DO 
     
    159179                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
    160180                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
     181!!gm old 
    161182                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    162183                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
     184!!gm new 
     185                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     186                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     187!!gm 
    163188               END DO 
    164189            END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r7753 r8093  
    928928               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    929929                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    930                   &            / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 
     930                  &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    931931            END DO 
    932932         END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r7990 r8093  
    4646 
    4747 
    48    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt  , avs     !: vertical T & S  diffusivity coef at w-pt [m2/s] 
    49    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm            !: vertical momentum viscosity coef at w-pt [m2/s] 
     48   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm, avt, avs  !: vertical mixing coefficient (w-point) [m2/s] 
     49!!gm 
    5050   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts  [m2/s] 
    51    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  !: avt, avs computed by turbulent closure alone 
     51!!gm 
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k, avt_k   !: avm, avt computed by turbulent closure alone 
    5253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
     54   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile 
    5355   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt 
    54    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile 
     56!!gm 
    5557   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: bottom friction coefficients 
    5658   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top    friction coefficients 
     59!!gm 
    5760 
    5861   !!---------------------------------------------------------------------- 
    59    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     62   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    6063   !! $Id$  
    6164   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r7990 r8093  
    3939CONTAINS 
    4040 
    41    SUBROUTINE zdf_ddm( kt ) 
     41   SUBROUTINE zdf_ddm( 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   ) ::   kt       ! ocean time-step indexocean time step 
     72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points) 
     73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points) 
     74      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points) 
    7275      ! 
    7376      INTEGER  ::   ji, jj , jk     ! dummy loop indices 
     
    9194!!gm                            and many acces in memory 
    9295          
    93          DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
     96         DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==! 
    9497            DO ji = 1, jpi 
    9598               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
     
    109112         END DO 
    110113 
    111          DO jj = 1, jpj                                     ! indicators: 
     114         DO jj = 1, jpj                !==  indicators  ==! 
    112115            DO ji = 1, jpi 
    113116               ! stability indicator: msks=1 if rn2>0; 0 elsewhere 
     
    154157                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    155158               ! 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 ) 
     159               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds 
     160               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt 
     161               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 
    159162            END DO 
    160163         END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r7990 r8093  
    3838CONTAINS 
    3939 
    40    SUBROUTINE zdf_evd( kt ) 
     40   SUBROUTINE zdf_evd( kt, p_avm, p_avt ) 
    4141      !!---------------------------------------------------------------------- 
    4242      !!                  ***  ROUTINE zdf_evd  *** 
     
    5555      !! ** Action  :   avt, avm   enhanced where static instability occurs 
    5656      !!---------------------------------------------------------------------- 
    57       INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step 
     57      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time-step indexocean time step 
     58      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    5859      ! 
    5960      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    7172      ! 
    7273      ! 
    73       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
     74      zavt_evd(:,:,:) = p_avt(:,:,:)         ! set avt prior to evd application 
    7475      ! 
    7576      SELECT CASE ( nn_evdm ) 
     
    7778      CASE ( 1 )           !==  enhance tracer & momentum Kz  ==!   (if rn2<-1.e-12) 
    7879         ! 
    79          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application 
     80         zavm_evd(:,:,:) = p_avm(:,:,:)      ! set avm prior to evd application 
    8081         ! 
    8182!! change last digits results 
    8283!         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) 
     84!            p_avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
     85!            p_avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
    8586!         END WHERE 
    8687         ! 
     
    8990               DO ji = 2, jpim1 
    9091                  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) 
     92                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
     93                     p_avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
    9394                  ENDIF 
    9495               END DO 
     
    9697         END DO  
    9798         ! 
    98          zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd 
    99          CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
     99         zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd 
     100         CALL iom_put( "avm_evd", zavm_evd )                ! output this change 
    100101         ! 
    101102      CASE DEFAULT         !==  enhance tracer Kz  ==!   (if rn2<-1.e-12)  
    102103!! change last digits results 
    103104!         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) 
     105!            p_avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 
    105106!         END WHERE 
    106107 
     
    109110               DO ji = 2, jpim1 
    110111                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   & 
    111                      avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
     112                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 
    112113               END DO 
    113114            END DO 
     
    116117      END SELECT  
    117118      ! 
    118       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
     119      zavt_evd(:,:,:) = p_avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
    119120      CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
    120121      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7991 r8093  
    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 
    22    USE zdfbfr         ! bottom friction (only for rn_bfrz0) 
     21!!gm old 
     22   USE zdfbfr  , ONLY : rn_tfrz0, rn_bfrz0   ! top/bottom roughness 
     23!!gm new 
     24   USE zdfdrg  , ONLY : r_z0_top , r_z0_bot   ! top/bottom roughness 
     25   USE zdfdrg  , ONLY : rCdU_top , rCdU_bot   ! top/bottom friction 
     26!!gm 
    2327   USE sbc_oce        ! surface boundary condition: ocean 
    2428   USE phycst         ! physical constants 
    2529   USE zdfmxl         ! mixed layer 
    26    USE sbcwave ,  ONLY: hsw   ! significant wave height 
     30   USE sbcwave , ONLY : hsw   ! significant wave height 
    2731   ! 
    2832   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    4549   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length 
    4650   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function 
    47    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_surf !: Squared surface velocity scale at T-points 
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_top  !: Squared top     velocity scale at T-points 
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_bot  !: Squared bottom  velocity scale at T-points 
    4954 
    5055   !                              !! ** Namelist  namzdf_gls  ** 
     
    101106   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        - 
    102107   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        - 
     108   ! 
     109   REAL(wp) ::   r2_3 = 2._wp/3._wp   ! constant=2/3 
    103110 
    104111   !! * Substitutions 
     
    115122      !!                ***  FUNCTION zdf_gls_alloc  *** 
    116123      !!---------------------------------------------------------------------- 
    117       ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustars2(jpi,jpj) ,       & 
    118          &      zwall (jpi,jpj,jpk) , ustarb2(jpi,jpj) ,  STAT= zdf_gls_alloc ) 
     124      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) ,                     & 
     125         &      zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc ) 
    119126         ! 
    120127      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    123130 
    124131 
    125    SUBROUTINE zdf_gls( kt ) 
     132   SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 
    126133      !!---------------------------------------------------------------------- 
    127134      !!                   ***  ROUTINE zdf_gls  *** 
     
    130137      !!              coefficients using the GLS turbulent closure scheme. 
    131138      !!---------------------------------------------------------------------- 
    132       INTEGER, INTENT(in) ::   kt ! ocean time step 
    133       INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments 
    134       REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars 
    135       REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -  
    136       REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      - 
     139      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step 
     140      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
     141      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
     142      ! 
     143      INTEGER  ::   ji, jj, jk    ! dummy loop arguments 
     144      INTEGER  ::   ibot, ibotm1  ! local integers 
     145      INTEGER  ::   itop, itopp1  !   -       - 
     146      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1 , zex2  ! local scalars 
     147      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof, zdir  !   -      -  
     148      REAL(wp) ::   zratio, zrn2, zflxb, sh     , z_en  !   -      - 
    137149      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      - 
    138       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 
     150      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      - 
     151      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep 
     152      REAL(wp), DIMENSION(jpi,jpj)     ::   zkar 
     153      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves  
     154      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
     155      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
     156      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before 
     157      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate 
     158      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi) 
     159      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now 
     160      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_a    ! element of the first  matrix diagonal 
     161      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_b    ! element of the second matrix diagonal 
     162      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_c    ! element of the third  matrix diagonal 
     163      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum 
    153164      !!-------------------------------------------------------------------- 
    154165      ! 
    155166      IF( nn_timing == 1 )   CALL timing_start('zdf_gls') 
    156167      ! 
    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  
    161168      ! Preliminary computing 
    162169 
    163       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    164  
    165       ! restore before values computed GLS alone 
    166       avt(:,:,:) = avt_k(:,:,:) 
    167       avm(:,:,:) = avm_k(:,:,:) 
    168  
    169       ! Compute surface and bottom friction at T-points 
     170      ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp    
     171      ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp 
     172      ustar2_bot (:,:) = 0._wp 
     173       
     174 
     175 
     176      ! Compute surface, top and bottom friction at T-points 
    170177      DO jj = 2, jpjm1           
    171178         DO ji = fs_2, fs_jpim1   ! vector opt.          
    172179            ! 
    173180            ! surface friction 
    174             ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
     181            ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
    175182            !    
     183!!gm old 
    176184            ! bottom friction (explicit before friction)         
    177185            ! Note that we chose here not to bound the friction as in dynbfr)    
     
    180188            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &          
    181189               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )       
    182             ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)          
     190            ustar2_bot(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)          
    183191         END DO          
    184192      END DO     
    185  
     193!!gm new 
     194!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that... 
     195!          ! bottom friction (explicit before friction) 
     196!          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     197!          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0) 
     198!          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 
     199!             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
     200!       END DO 
     201!    END DO 
     202!    IF( ln_isfcav ) THEN       !top friction 
     203!       DO jj = 2, jpjm1 
     204!          DO ji = fs_2, fs_jpim1   ! vector opt. 
     205!             zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
     206!             zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0) 
     207!             ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 
     208!                &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
     209!          END DO 
     210!       END DO 
     211!    ENDIF 
     212!!gm 
    186213      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==! 
    187214      CASE ( 0 )                          ! Constant roughness           
    188215         zhsro(:,:) = rn_hsro 
    189216      CASE ( 1 )             ! Standard Charnock formula 
    190          zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 
     217         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro ) 
    191218      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 
    192219!!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) 
     220!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )      ! Wave age (eq. 10) 
     221         zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )            ! Wave age (eq. 10) 
     222         zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    196223      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
    197224!!gm  BUG missing a multiplicative coefficient.... 
    198225         zhsro(:,:) = hsw(:,:) 
    199226      END SELECT 
    200  
    201       !                             !==  Compute shear and dissipation rate  ==! 
    202       CALL zdf_sh2( shear ) 
    203  
    204  
     227      ! 
    205228      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==! 
    206229         DO jj = 1, jpjm1 
     
    245268            DO ji = 2, jpim1 
    246269               ! 
    247                buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction 
     270               buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction 
    248271               ! 
    249272               diss = eps(ji,jj,jk)                         ! dissipation 
    250273               ! 
    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 
    254                zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
     274               zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     275               ! 
     276               zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term 
     277               zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 
     278!!gm better coding, identical results 
     279!               zesh2 =   p_sh2(ji,jj,jk) + zdir*buoy               ! production term 
     280!               zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 
     281!!gm 
    255282               ! 
    256283               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 
     
    269296               zcof = rfact_tke * tmask(ji,jj,jk) 
    270297               !                                               ! 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) ) 
     298               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) ) 
    272299               !                                               ! 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) ) 
     300               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) ) 
    274301               !                                               ! diagonal 
    275302               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)  
     
    293320      CASE ( 0 )             ! Dirichlet case 
    294321      ! First level 
    295       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     322      en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 
    296323      en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    297324      z_elem_a(:,:,1) = en(:,:,1) 
     
    300327      !  
    301328      ! One level below 
    302       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 
     329      en(:,:,2) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 
    303330         &               / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    304331      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     
    311338      ! 
    312339      ! Dirichlet conditions at k=1 
    313       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     340      en(:,:,1)       = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 
    314341      en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    315342      z_elem_a(:,:,1) = en(:,:,1) 
     
    322349      z_elem_a(:,:,2) = 0._wp 
    323350      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
    324       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
     351      zflxs(:,:)      = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    325352          &                       * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    326353 
     
    340367         DO jj = 2, jpjm1 
    341368            DO ji = fs_2, fs_jpim1   ! vector opt. 
     369!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ???? 
     370!!   With thick deep ocean level thickness, this may be quite large, no ??? 
     371!!   in particular in ocean cavities where top stratification can be large... 
    342372               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    343373               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    344374               ! 
    345                ! Bottom level Dirichlet condition: 
    346                z_elem_a(ji,jj,ibot  ) = 0._wp 
    347                z_elem_c(ji,jj,ibot  ) = 0._wp 
    348                z_elem_b(ji,jj,ibot  ) = 1._wp 
    349                en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    350                ! 
    351                ! Just above last level, Dirichlet condition again 
    352                z_elem_a(ji,jj,ibotm1) = 0._wp 
    353                z_elem_c(ji,jj,ibotm1) = 0._wp 
    354                z_elem_b(ji,jj,ibotm1) = 1._wp 
    355                en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )  
    356             END DO 
    357          END DO 
     375               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
     376               ! 
     377               ! Dirichlet condition applied at:  
     378               !     Bottom level (ibot)      &      Just above it (ibotm1)    
     379               z_elem_a(ji,jj,ibot) = 0._wp   ;   z_elem_a(ji,jj,ibotm1) = 0._wp 
     380               z_elem_c(ji,jj,ibot) = 0._wp   ;   z_elem_c(ji,jj,ibotm1) = 0._wp 
     381               z_elem_b(ji,jj,ibot) = 1._wp   ;   z_elem_b(ji,jj,ibotm1) = 1._wp 
     382               en      (ji,jj,ibot) = z_en    ;   en      (ji,jj,ibotm1) = z_en 
     383            END DO 
     384         END DO 
     385!!gm new 
     386!         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     387!            DO jj = 2, jpjm1 
     388!               DO ji = fs_2, fs_jpim1   ! vector opt. 
     389!                  itop   = mikt(ji,jj)       ! k   top w-point 
     390!                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     391!                  !                                                ! mask at the ocean surface points 
     392!                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
     393!                  ! 
     394!                  ! Dirichlet condition applied at:  
     395!                  !     top level (itop)         &      Just below it (itopp1)    
     396!                  z_elem_a(ji,jj,itop) = 0._wp   ;   z_elem_a(ji,jj,ipotm1) = 0._wp 
     397!                  z_elem_c(ji,jj,itop) = 0._wp   ;   z_elem_c(ji,jj,itopp1) = 0._wp 
     398!                  z_elem_b(ji,jj,itop) = 1._wp   ;   z_elem_b(ji,jj,itopp1) = 1._wp 
     399!                  en      (ji,jj,itop) = zen     ;   en      (ji,jj,itopp1) = z_en 
     400!               END DO 
     401!            END DO 
     402!         ENDIF 
     403!!gm 
    358404         ! 
    359405      CASE ( 1 )             ! Neumman boundary condition 
     
    364410               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    365411               ! 
     412               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 
     413               ! 
    366414               ! Bottom level Dirichlet condition: 
    367                z_elem_a(ji,jj,ibot) = 0._wp 
    368                z_elem_c(ji,jj,ibot) = 0._wp 
    369                z_elem_b(ji,jj,ibot) = 1._wp 
    370                en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    371                ! 
    372                ! Just above last level: Neumann condition 
    373                z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)   ! Remove z_elem_c from z_elem_b 
    374                z_elem_c(ji,jj,ibotm1) = 0._wp 
    375             END DO 
    376          END DO 
     415               !     Bottom level (ibot)      &      Just above it (ibotm1)    
     416               !         Dirichlet            !         Neumann 
     417               z_elem_a(ji,jj,ibot) = 0._wp   !   ! Remove z_elem_c from z_elem_b 
     418               z_elem_b(ji,jj,ibot) = 1._wp   ;   z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) 
     419               z_elem_c(ji,jj,ibot) = 0._wp   ;   z_elem_c(ji,jj,ibotm1) = 0._wp 
     420            END DO 
     421         END DO 
     422!!gm new 
     423!         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity) 
     424!            DO jj = 2, jpjm1 
     425!               DO ji = fs_2, fs_jpim1   ! vector opt. 
     426!                  itop   = mikt(ji,jj)       ! k   top w-point 
     427!                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one 
     428!                  !                                                ! mask at the ocean surface points 
     429!                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 
     430!                  ! 
     431!                  ! Bottom level Dirichlet condition: 
     432!                  !     Bottom level (ibot)      &      Just above it (ibotm1)    
     433!                  !         Dirichlet            !         Neumann 
     434!                  z_elem_a(ji,jj,itop) = 0._wp   !   ! Remove z_elem_c from z_elem_b 
     435!                  z_elem_b(ji,jj,itop) = 1._wp   ;   z_elem_b(ji,jj,itopp1) = z_elem_b(ji,jj,itopp1) + z_elem_c(ji,jj,itopp1) 
     436!                  z_elem_c(ji,jj,itop) = 0._wp   ;   z_elem_c(ji,jj,itopp1) = 0._wp 
     437!               END DO 
     438!            END DO 
     439!         ENDIF 
     440!!gm 
    377441         ! 
    378442      END SELECT 
     
    465529               zratio = psi(ji,jj,jk) / eb(ji,jj,jk)  
    466530               ! 
    467                ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) 
    468                dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 
    469                ! 
    470                rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p 
     531               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 
     532               zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 
     533               ! 
     534               rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 
    471535               ! 
    472536               ! shear prod. - stratif. destruction 
    473                prod = rpsi1 * zratio * shear(ji,jj,jk) 
     537               prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 
    474538               ! 
    475539               ! stratif. destruction 
    476                buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
     540               buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 
    477541               ! 
    478542               ! shear prod. - stratif. destruction 
    479543               diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 
    480544               ! 
    481                dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
    482                ! 
    483                zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term 
    484                zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
     545               zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 
     546               ! 
     547               zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term 
     548               zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
    485549               !                                                         
    486550               ! building the matrix 
    487551               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
    488552               !                                               ! 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) ) 
     553               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) ) 
    490554               !                                               ! 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) ) 
     555               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) ) 
    492556               !                                               ! diagonal 
    493557               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) 
     
    506570      ! 
    507571      CASE ( 0 )             ! Dirichlet boundary conditions 
    508       ! 
    509       ! 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 
    515       ! 
    516       ! 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 
    523       !  
    524       ! 
     572         ! 
     573         ! Surface value 
     574         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic 
     575         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     576         z_elem_a(:,:,1) = psi(:,:,1) 
     577         z_elem_c(:,:,1) = 0._wp 
     578         z_elem_b(:,:,1) = 1._wp 
     579         ! 
     580         ! One level below 
     581         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 
     582         zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 
     583         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     584         z_elem_a(:,:,2) = 0._wp 
     585         z_elem_c(:,:,2) = 0._wp 
     586         z_elem_b(:,:,2) = 1._wp 
     587         !  
    525588      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    526       ! 
    527       ! 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 
    533       ! 
    534       ! 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 
    537       ! 
    538       ! 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) 
    546       !    
    547       ! 
     589         ! 
     590         ! Surface value: Dirichlet 
     591         zdep    (:,:)   = zhsro(:,:) * rl_sf 
     592         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     593         z_elem_a(:,:,1) = psi(:,:,1) 
     594         z_elem_c(:,:,1) = 0._wp 
     595         z_elem_b(:,:,1) = 1._wp 
     596         ! 
     597         ! Neumann condition at k=2 
     598         z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     599         z_elem_a(:,:,2) = 0._wp 
     600         ! 
     601         ! Set psi vertical flux at the surface: 
     602         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     603         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     604         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     605         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
     606            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 
     607         zflxs(:,:)   = zdep(:,:) * zflxs(:,:) 
     608         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 
     609         ! 
    548610      END SELECT 
    549611 
     
    552614      ! 
    553615      SELECT CASE ( nn_bc_bot ) 
    554       ! 
    555616      ! 
    556617      CASE ( 0 )             ! Dirichlet  
     
    597658               ! Set psi vertical flux at the bottom: 
    598659               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) )   & 
     660               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   & 
    600661                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
    601662               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 
     
    730791                  gh = gh * rf6 
    731792                  ! Gm =  M²l²/q² Shear number 
    732                   shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall ) 
     793                  shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 
    733794                  gm = MAX( shr * zcof , 1.e-10 ) 
    734795                  gm = gm * rf6 
     
    765826         DO jj = 2, jpjm1 
    766827            DO ji = fs_2, fs_jpim1   ! vector opt. 
    767                zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
    768                zav           = zsqen * zstt(ji,jj,jk) 
    769                avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
    770                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 
     828               zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 
     829               zavt  = zsqen * zstt(ji,jj,jk) 
     830               zavm  = zsqen * zstm(ji,jj,jk) 
     831               p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 
     832               p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom 
    772833            END DO 
    773834         END DO 
    774835      END DO 
    775836      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. ) 
    778837      ! 
    779838      IF(ln_ctl) THEN 
     
    781840         CALL prt_ctl( tab3d_1=avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk ) 
    782841      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 ) 
    792842      ! 
    793843      IF( nn_timing == 1 )   CALL timing_stop('zdf_gls') 
     
    837887      IF(lwp) THEN                     !* Control print 
    838888         WRITE(numout,*) 
    839          WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme' 
     889         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme' 
    840890         WRITE(numout,*) '~~~~~~~~~~~~' 
    841891         WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters' 
     
    854904         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    855905         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
     906!!gm old 
     907         WRITE(numout,*) '      top    roughness (m) (nambfr namelist)        rn_tfrz0       = ', rn_tfrz0 
    856908         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0 
     909!!gm new 
     910!         WRITE(numout,*) '   Namelist namdrg_top/_bot used values:' 
     911!         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top 
     912!         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot 
     913!!gm 
     914         WRITE(numout,*) 
    857915      ENDIF 
    858916 
     
    861919 
    862920      !                                !* 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' ) 
     921      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' ) 
     922      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' ) 
     923      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' ) 
     924      IF( nn_z0_met == 3  .AND. .NOT.ln_wave     )  CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
     925      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' ) 
     926      IF( nn_clos      < 0 .OR. nn_clos      > 3 )  CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    869927 
    870928      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
     
    872930      CASE( 0 )                              ! k-kl  (Mellor-Yamada) 
    873931         ! 
    874          IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 
     932         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)' 
     933         IF(lwp) WRITE(numout,*) 
    875934         rpp     = 0._wp 
    876935         rmm     = 1._wp 
     
    890949      CASE( 1 )                              ! k-eps 
    891950         ! 
    892          IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 
     951         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen' 
     952         IF(lwp) WRITE(numout,*) 
    893953         rpp     =  3._wp 
    894954         rmm     =  1.5_wp 
     
    908968      CASE( 2 )                              ! k-omega 
    909969         ! 
    910          IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 
     970         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen' 
     971         IF(lwp) WRITE(numout,*) 
    911972         rpp     = -1._wp 
    912973         rmm     =  0.5_wp 
     
    926987      CASE( 3 )                              ! generic 
    927988         ! 
    928          IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 
     989         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen' 
     990         IF(lwp) WRITE(numout,*) 
    929991         rpp     = 2._wp 
    930992         rmm     = 1._wp 
     
    9491011      CASE ( 0 )                             ! Galperin stability functions 
    9501012         ! 
    951          IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 
     1013         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin' 
    9521014         rc2     =  0._wp 
    9531015         rc3     =  0._wp 
     
    9611023      CASE ( 1 )                             ! Kantha-Clayson stability functions 
    9621024         ! 
    963          IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 
     1025         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson' 
    9641026         rc2     =  0.7_wp 
    9651027         rc3     =  0.2_wp 
     
    9731035      CASE ( 2 )                             ! Canuto A stability functions 
    9741036         ! 
    975          IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 
     1037         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A' 
    9761038         rs0 = 1.5_wp * rl1 * rl5*rl5 
    9771039         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 
     
    9971059      CASE ( 3 )                             ! Canuto B stability functions 
    9981060         ! 
    999          IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 
     1061         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B' 
    10001062         rs0 = 1.5_wp * rm1 * rm5*rm5 
    10011063         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 
     
    10521114      IF(lwp) THEN                     !* Control print 
    10531115         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 
     1116         WRITE(numout,*) '   Limit values :' 
     1117         WRITE(numout,*) '      Parameter  m = ', rmm 
     1118         WRITE(numout,*) '      Parameter  n = ', rnn 
     1119         WRITE(numout,*) '      Parameter  p = ', rpp 
     1120         WRITE(numout,*) '      rpsi1    = ', rpsi1 
     1121         WRITE(numout,*) '      rpsi2    = ', rpsi2 
     1122         WRITE(numout,*) '      rpsi3m   = ', rpsi3m 
     1123         WRITE(numout,*) '      rpsi3p   = ', rpsi3p 
     1124         WRITE(numout,*) '      rsc_tke  = ', rsc_tke 
     1125         WRITE(numout,*) '      rsc_psi  = ', rsc_psi 
     1126         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0 
     1127         WRITE(numout,*) '      rc0      = ', rc0 
    10671128         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,*) 
     1129         WRITE(numout,*) '   Shear free turbulence parameters:' 
     1130         WRITE(numout,*) '      rcm_sf   = ', rcm_sf 
     1131         WRITE(numout,*) '      ra_sf    = ', ra_sf 
     1132         WRITE(numout,*) '      rl_sf    = ', rl_sf 
    10731133      ENDIF 
    10741134 
     
    10901150      ! 
    10911151      !                                !* Wall proximity function 
    1092       zwall (:,:,:) = 1._wp * tmask(:,:,:) 
     1152!!gm tmask or wmask ???? 
     1153      zwall(:,:,:) = 1._wp * tmask(:,:,:) 
    10931154 
    10941155      !                                !* read or initialize all required files   
     
    11331194               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 
    11341195            ELSE                         
    1135                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 
     1196               IF(lwp) WRITE(numout,*) 
     1197               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values' 
    11361198               en    (:,:,:) = rn_emin 
    11371199               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 
     1200               ! avt_k, avm_k already set to the background value in zdf_phy_init 
    11411201            ENDIF 
    11421202         ELSE                                   !* Start from rest 
    1143             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 
     1203            IF(lwp) WRITE(numout,*) 
     1204            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values' 
    11441205            en    (:,:,:) = rn_emin 
    11451206            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 
     1207            ! avt_k, avm_k already set to the background value in zdf_phy_init 
    11501208         ENDIF 
    11511209         ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

    r7990 r8093  
    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 * 
     
    7877 
    7978 
    80    SUBROUTINE zdf_iwm( kt ) 
     79   SUBROUTINE zdf_iwm( 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   ) ::   kt             ! ocean time step 
     129      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
     130      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    130131      ! 
    131132      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 
     133      REAL(wp) ::   zztmp        ! scalar workspace 
     134      REAL(wp), DIMENSION(jpi,jpj)    ::  zfact     ! Used for vertical structure 
     135      REAL(wp), DIMENSION(jpi,jpj)    ::  zhdep     ! Ocean depth 
     136      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwkb      ! WKB-stretched height above bottom 
     137      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zweight   ! Weight for high mode vertical distribution 
     138      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     140      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zReb      ! Turbulence intensity parameter 
    140141      !!---------------------------------------------------------------------- 
    141142      ! 
    142143      IF( nn_timing == 1 )   CALL timing_start('zdf_iwm') 
    143144      ! 
    144       CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
    145       CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    146  
    147       !                          ! ----------------------------- ! 
    148       !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
    149       !                          ! ----------------------------- ! 
     145      !                      ! ----------------------------- ! 
     146      !                      !  Internal wave-driven mixing  !  (compute zav_wave) 
     147      !                      ! ----------------------------- ! 
    150148      !                              
    151149      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
     
    162160         emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      & 
    163161            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   & 
     162            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     163 
    164164!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point 
    165             &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     165!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all 
     166 
    166167      END DO 
    167168 
    168169      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
    169170      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
    170       ! 
     171      !                                          ! (NB: N2 is masked, so no use of wmask here) 
    171172      SELECT CASE ( nn_zpyc ) 
    172173      ! 
     
    210211      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    211212      ! 
    212       zwkb(:,:,:) = 0._wp 
    213       zfact(:,:) = 0._wp 
     213      zwkb (:,:,:) = 0._wp 
     214      zfact(:,:)   = 0._wp 
    214215      DO jk = 2, jpkm1 
    215216         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    216217         zwkb(:,:,jk) = zfact(:,:) 
    217218      END DO 
     219!!gm even better: 
     220!      DO jk = 2, jpkm1 
     221!         zwkb(:,:) = zwkb(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) 
     222!      END DO 
     223!      zfact(:,:) = zwkb(:,:,jpkm1) 
     224!!gm or just use zwkb(k=jpk-1) instead of zfact... 
     225!!gm 
    218226      ! 
    219227      DO jk = 2, jpkm1 
     
    221229            DO ji = 1, jpi 
    222230               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) 
    224             END DO 
    225          END DO 
    226       END DO 
    227       zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
     231                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
     232            END DO 
     233         END DO 
     234      END DO 
     235      zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 
    228236      ! 
    229237      zweight(:,:,:) = 0._wp 
     
    247255         emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    248256            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    249       END DO 
    250       ! 
     257!!gm  use of e3t_n just above? 
     258      END DO 
     259      ! 
     260!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    251261      ! Calculate molecular kinematic viscosity 
    252262      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
     
    255265         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    256266      END DO 
     267!!gm end 
    257268      ! 
    258269      ! Calculate turbulence intensity parameter Reb 
     
    285296      ! 
    286297      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    287          ztpc = 0._wp 
     298         zztmp = 0._wp 
    288299!!gm used of glosum 3D.... 
    289300         DO jk = 2, jpkm1 
    290301            DO jj = 1, jpj 
    291302               DO ji = 1, jpi 
    292                   ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
    293                      &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     303                  zztmp = zztmp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
     304                     &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    294305               END DO 
    295306            END DO 
    296307         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  
     308         IF( lk_mpp )   CALL mpp_sum( zztmp ) 
     309         zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    299310         ! 
    300311         IF(lwp) THEN 
     
    303314            WRITE(numout,*) '~~~~~~~ ' 
    304315            WRITE(numout,*) 
    305             WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     316            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztmp * 1.e-12_wp, 'TW' 
    306317         ENDIF 
    307318      ENDIF 
     
    323334         CALL iom_put( "av_ratio", zav_ratio ) 
    324335         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) 
     336            p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
     337            p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
     338            p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    328339         END DO 
    329340         ! 
    330341      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    331342         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) 
     343            p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 
     344            p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
     345            p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    335346         END DO 
    336347      ENDIF 
     
    353364      CALL iom_put( "emix_iwm", emix_iwm ) 
    354365       
    355       CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
    356       CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    357  
    358366      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
    359367      ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

    r7990 r8093  
    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          
     15   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     16!!gm old 
    1517   USE zdfbfr         ! vertical physics: bottom friction 
    16    USE zdfric         ! vertical physics: Richardson vertical mixing    
     18!!gm 
     19   USE zdfsh2         ! vertical physics: shear production term of TKE 
     20   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing    
    1721   USE zdftke         ! vertical physics: TKE vertical mixing 
    1822   USE zdfgls         ! vertical physics: GLS vertical mixing 
     
    4448   INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz 
    4549   INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz 
    46        
     50 
     51   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 
     52 
    4753   !!---------------------------------------------------------------------- 
    4854   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
     
    6167      !!                set horizontal shape and vertical profile of background mixing coef. 
    6268      !!---------------------------------------------------------------------- 
     69      INTEGER ::   jk            ! dummy loop indices 
    6370      INTEGER ::   ioptio, ios   ! local integers 
    6471      !! 
     
    9097         WRITE(numout,*) '      vertical closure scheme' 
    9198         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 
     99         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric 
     100         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke 
     101         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls 
    95102         WRITE(numout,*) '      convection: ' 
    96103         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd 
     
    137144      ENDIF 
    138145      ! 
    139       !                                   ! set 1st/last level Av to zero one for all 
    140       avt(:,:,1) = 0._wp   ;   avs(:,:,1) = 0._wp   ;   avm(:,:,1) = 0._wp 
     146      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k) 
     147         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk) 
     148         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk) 
     149      END DO 
     150!!gm  to be tested only the 1st & last levels 
     151!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp 
     152!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp 
     153!!gm 
     154      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp 
    141155 
    142156      !                          !==  Convection  ==! 
     
    170184      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' ) 
    171185      ENDIF 
     186      !                                ! shear production term flag 
     187      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE. 
     188      ELSE                   ;   l_zdfsh2 = .TRUE. 
     189      ENDIF 
    172190 
    173191      !                          !== gravity wave-driven mixing  ==! 
     
    175193      IF( ln_zdfswm )   CALL zdf_swm_init       ! surface  wave-driven mixing 
    176194 
    177       !                          !== bottom friction  ==! 
     195      !                          !== top/bottom friction  ==! 
     196      CALL zdf_drg_init 
     197!!gm old 
    178198      CALL zdf_bfr_init 
     199!!gm 
    179200      ! 
    180201      !                          !== time-stepping  ==! 
     
    200221      ! 
    201222      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     223      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2   ! shear production 
    202224      !! --------------------------------------------------------------------- 
    203225      ! 
     226!!gm old 
    204227      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic) 
    205       !                        
     228!!gm 
     229      ! 
     230!      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases) 
     231!         ! 
     232!         !                       !* bottom drag 
     233!         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in  
     234!            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   & 
     235!            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s] 
     236!         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities) 
     237!            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in  
     238!               &              r_z0_top,   r_ke0_top,    rCd0_top,   & 
     239!               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s] 
     240!         ENDIF 
     241!      ENDIF 
     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( ub, vb, un, vn, avm_k,   &     ! <<== in 
     247            &                           zsh2    )     ! ==>> out : shear production 
     248      ! 
    206249      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 
    210       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) 
     250      CASE( np_RIC )   ;   CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz 
     251      CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz 
     252      CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz 
     253!     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value) 
     254!         ! avt_k and avm_k set one for all at initialisation phase 
     255!!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
     256!!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 
    213257      END SELECT 
     258      !   
     259      !                          !==  ocean Kz  ==!   (avt, avs, avm) 
     260      ! 
     261      !                                         !* start from turbulent closure values 
     262      avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1) 
     263      avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1) 
    214264      ! 
    215265      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths 
     
    219269      ENDIF 
    220270      ! 
    221       IF( ln_zdfevd )   CALL zdf_evd( kt )      !* convection: enhanced vertical eddy diffusivity 
     271      IF( ln_zdfevd )   CALL zdf_evd( kt, avm, avt )  !* convection: enhanced vertical eddy diffusivity 
    222272      ! 
    223273      !                                         !* double diffusive mixing 
    224274      IF( ln_zdfddm ) THEN                            ! update avt and compute avs 
    225                         CALL zdf_ddm( kt ) 
     275                        CALL zdf_ddm( kt, avm, avt, avs ) 
    226276      ELSE                                            ! same mixing on all tracers 
    227277         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) 
     
    229279      ! 
    230280      !                                         !* 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) 
     281      IF( ln_zdfswm )   CALL zdf_swm( kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)  
     282      IF( ln_zdfiwm )   CALL zdf_iwm( kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017) 
    233283 
    234284 
    235285      !                                         !* 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 
     286      CALL lbc_lnk( avm_k, 'W', 1. )                  ! needed to compute the shear production term 
     287      CALL lbc_lnk( avt_k, 'W', 1. )                  !!gm a priori useless ==>> to be tested 
     288      CALL lbc_lnk( avm  , 'W', 1. )                  ! needed to compute avm at u- and v-points 
     289      CALL lbc_lnk( avt  , 'W', 1. )                  !!gm  a priori only avm_k and avm are required 
     290      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  To be tested 
     291      ! 
    251292 
    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      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file 
     297         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' ) 
     298         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' ) 
     299         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' )  
     300      ENDIF 
     301      ! 
    260302   END SUBROUTINE zdf_phy 
    261303 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r7991 r8093  
    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 
    27    USE eosbn2 ,  ONLY :   neos 
    2827   ! 
    2928   USE in_out_manager ! I/O manager 
    30    USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    31    USE lib_mpp        ! MPP library 
     29   USE iom            ! I/O manager library 
    3230   USE timing         ! Timing 
    3331   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     
    3735   PRIVATE 
    3836 
    39    PUBLIC   zdf_ric         ! called by step.F90 
    40    PUBLIC   zdf_ric_init    ! called by opa.F90 
     37   PUBLIC   zdf_ric         ! called by zdfphy.F90 
     38   PUBLIC   ric_rst         ! called by zdfphy.F90 
     39   PUBLIC   zdf_ric_init    ! called by nemogcm.F90 
    4140 
    4241   !                        !!* Namelist namzdf_ric : Richardson number dependent Kz * 
     
    5453#  include "vectopt_loop_substitute.h90" 
    5554   !!---------------------------------------------------------------------- 
    56    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     55   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    5756   !! $Id$ 
    5857   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5958   !!---------------------------------------------------------------------- 
    6059CONTAINS 
    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  
    16160 
    16261   SUBROUTINE zdf_ric_init 
     
    205104      ENDIF 
    206105      ! 
    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) 
     106      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files 
     107      ! 
     108   END SUBROUTINE zdf_ric_init 
     109 
     110 
     111   SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt ) 
     112      !!---------------------------------------------------------------------- 
     113      !!                 ***  ROUTINE zdfric  *** 
     114      !!                     
     115      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as 
     116      !!                a function of the local richardson number. 
     117      !! 
     118      !! ** Method  :   Local richardson number dependent formulation of the  
     119      !!                vertical eddy viscosity and diffusivity coefficients.  
     120      !!                The eddy coefficients are given by: 
     121      !!                    avm = avm0 + avmb 
     122      !!                    avt = avm0 / (1 + rn_alp*ri) 
     123      !!                with ri  = N^2 / dz(u)**2 
     124      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 
     125      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 
     126      !!                where ri is the before local Richardson number, 
     127      !!                rn_avmri is the maximum value reaches by avm and avt  
     128      !!                and rn_alp, nn_ric are adjustable parameters. 
     129      !!                Typical values : rn_alp=5. and nn_ric=2. 
     130      !! 
     131      !!      As second step compute Ekman depth from wind stress forcing 
     132      !!      and apply namelist provided vertical coeff within this depth. 
     133      !!      The Ekman depth is: 
     134      !!              Ustar = SQRT(Taum/rho0) 
     135      !!              ekd= rn_ekmfc * Ustar / f0 
     136      !!      Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation 
     137      !!      of the above equation indicates the value is somewhat arbitrary; therefore 
     138      !!      we allow the freedom to increase or decrease this value, if the 
     139      !!      Ekman depth estimate appears too shallow or too deep, respectively. 
     140      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the 
     141      !!      namelist 
     142      !!        N.B. the mask are required for implicit scheme, and surface 
     143      !!      and bottom value already set in zdfphy.F90 
     144      !! 
     145      !! ** Action  :   avm, avt  mixing coeff (inner domain values only) 
     146      !! 
     147      !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 
     148      !!              PFJ Lermusiaux 2001. 
     149      !!---------------------------------------------------------------------- 
     150      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step 
     151      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdept          ! depth of t-point  [m] 
     152      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
     153      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points) 
     154      !! 
     155      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
     156      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars 
     157      REAL(wp), DIMENSION(jpi,jpj) ::   zh_ekm  ! 2D workspace 
     158      !!---------------------------------------------------------------------- 
     159      ! 
     160      IF( nn_timing == 1 )   CALL timing_start('zdf_ric') 
     161      ! 
     162      !                       !==  avm and avt = F(Richardson number)  ==! 
     163      DO jk = 2, jpkm1 
     164         DO jj = 1, jpjm1 
     165            DO ji = 1, jpim1              ! coefficient = F(richardson number) (avm-weighted Ri) 
     166               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 ) )  ) 
     167               zav   = rn_avmri * zcfRi**nn_ric 
     168               !                          ! avm and avt coefficients 
     169               p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk) 
     170               p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk) 
     171            END DO 
     172         END DO 
    210173      END DO 
    211174      ! 
    212    END SUBROUTINE zdf_ric_init 
     175!!gm BUG <<<<====  This param can't work at low latitude  
     176!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 
     177      ! 
     178      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==! 
     179         ! 
     180         DO jj = 2, jpjm1        !* Ekman depth 
     181            DO ji = 2, jpim1 
     182               zustar = SQRT( taum(ji,jj) * r1_rau0 ) 
     183               zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth 
     184               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range 
     185            END DO 
     186         END DO 
     187         DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer 
     188            DO jj = 2, jpjm1 
     189               DO ji = 2, jpim1 
     190                  IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 
     191                     p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk) 
     192                     p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk) 
     193                  ENDIF 
     194               END DO 
     195            END DO 
     196         END DO 
     197      ENDIF 
     198      ! 
     199      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric') 
     200      ! 
     201   END SUBROUTINE zdf_ric 
     202 
     203 
     204   SUBROUTINE ric_rst( kt, cdrw ) 
     205      !!--------------------------------------------------------------------- 
     206      !!                   ***  ROUTINE ric_rst  *** 
     207      !!                      
     208      !! ** Purpose :   Read or write TKE file (en) in restart file 
     209      !! 
     210      !! ** Method  :   use of IOM library 
     211      !!                if the restart does not contain TKE, en is either  
     212      !!                set to rn_emin or recomputed  
     213      !!---------------------------------------------------------------------- 
     214      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     215      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     216      ! 
     217      INTEGER ::   jit, jk    ! dummy loop indices 
     218      INTEGER ::   id1, id2   ! local integers 
     219      !!---------------------------------------------------------------------- 
     220      ! 
     221      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     222         !                                   ! --------------- 
     223         !           !* Read the restart file 
     224         IF( ln_rstart ) THEN 
     225            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
     226            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
     227            ! 
     228            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it 
     229               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
     230               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
     231            ENDIF 
     232         ENDIF 
     233         !           !* otherwise Kz already set to the background value in zdf_phy_init 
     234         ! 
     235      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     236         !                                   ! ------------------- 
     237         IF(lwp) WRITE(numout,*) '---- ric-rst ----' 
     238         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
     239         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
     240         ! 
     241      ENDIF 
     242      ! 
     243   END SUBROUTINE ric_rst 
    213244 
    214245   !!====================================================================== 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90

    r7997 r8093  
    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 
     
    3129CONTAINS 
    3230 
    33    SUBROUTINE zdf_sh2(  psh2  )  
     31   SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2  )  
    3432      !!---------------------------------------------------------------------- 
    3533      !!                   ***  ROUTINE zdf_sh2  *** 
     
    4644      !!                NB: wet-point only horizontal averaging of shear 
    4745      !! 
    48       !! ** Action  : - psh2 shear prod. term at w-point (interior ocean domain only) 
    49       !! 
     46      !! ** Action  : - p_sh2 shear prod. term at w-point (inner domain only) 
     47      !!                                                   ***** 
    5048      !! References :   Bruchard, OM 2002 
    5149      !! --------------------------------------------------------------------- 
    52       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   psh2   ! shear production of TKE (w-points) 
     50      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities 
     51      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points) 
     52      REAL(wp), DIMENSION(:,:,:) , INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points) 
    5353      ! 
    5454      INTEGER  ::   ji, jj, jk   ! dummy loop arguments 
    5555      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace 
    5656      !!-------------------------------------------------------------------- 
     57      ! 
    5758      IF( nn_timing == 1 )  CALL timing_start('zdf_sh2') 
    5859      ! 
     
    6061         DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form) 
    6162            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) 
     63               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 
     64                  &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) & 
     65                  &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk) 
     66               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 
     67                  &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) & 
     68                  &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk) 
    6869            END DO 
    6970         END DO 
    7071         DO jj = 2, jpjm1        !* shear production at w-point 
    7172            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) )   ) 
     73               p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   & 
     74                  &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   ) 
    7475            END DO 
    7576         END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90

    r7990 r8093  
    3434CONTAINS 
    3535 
    36    SUBROUTINE zdf_swm( kt ) 
     36   SUBROUTINE zdf_swm( kt, p_avm, p_avt, p_avs ) 
    3737      !!--------------------------------------------------------------------- 
    3838      !!                     ***  ROUTINE zdf_swm *** 
     
    5151      !! reference : Qiao et al. GRL, 2004 
    5252      !!--------------------------------------------------------------------- 
    53       INTEGER, INTENT(in) ::   kt   ! ocean time step 
     53      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     54      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
     55      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    5456      ! 
    5557      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    6365               zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk) 
    6466               ! 
    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 
     67               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zqb 
     68               p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zqb 
     69               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zqb 
    6870            END DO 
    6971         END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r7991 r8093  
    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 
     45!!gm new 
     46   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
     47!!gm 
    4648   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)   
    5549#if defined key_agrif 
    5650   USE agrif_opa_interp 
    5751   USE agrif_opa_update 
    5852#endif 
     53   ! 
     54   USE in_out_manager ! I/O manager 
     55   USE iom            ! I/O manager library 
     56   USE lib_mpp        ! MPP library 
     57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     58   USE prtctl         ! Print control 
     59   USE wrk_nemo       ! work arrays 
     60   USE timing         ! Timing 
     61   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    5962 
    6063   IMPLICIT NONE 
     
    8790   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    8891 
    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 
     92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau) 
     93   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation 
     94   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation 
    9295 
    9396   !! * Substitutions 
     
    112115 
    113116 
    114    SUBROUTINE zdf_tke( kt ) 
     117   SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 
    115118      !!---------------------------------------------------------------------- 
    116119      !!                   ***  ROUTINE zdf_tke  *** 
     
    157160      !!              Bruchard OM 2002 
    158161      !!---------------------------------------------------------------------- 
    159       INTEGER, INTENT(in) ::   kt   ! ocean time step 
     162      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     163      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term 
     164      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points) 
    160165      !!---------------------------------------------------------------------- 
    161166      ! 
     
    165170#endif 
    166171      ! 
    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(:,:,:)  
     172      CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )   ! now tke (en) 
     173      ! 
     174      CALL tke_avn( gdepw_n, e3t_n, e3w_n,        p_avm, p_avt )   ! now avt, avm, dissl 
    176175      ! 
    177176#if defined key_agrif 
     
    179178      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    180179#endif       
    181       !  
    182       IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' ) 
    183180      ! 
    184181  END SUBROUTINE zdf_tke 
    185182 
    186183 
    187    SUBROUTINE tke_tke 
     184   SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2    & 
     185      &                            , p_avm, p_avt ) 
    188186      !!---------------------------------------------------------------------- 
    189187      !!                   ***  ROUTINE tke_tke  *** 
     
    200198      !! ** Action  : - en : now turbulent kinetic energy) 
    201199      !! --------------------------------------------------------------------- 
    202       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
    203 !!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! local integers 
    204 !!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          !   -      - 
    205 !!bfr      REAL(wp) ::   zebot                           ! local scalars 
     200      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! depth of w-points 
     201      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     202      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_sh2          ! shear production term 
     203      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
     204      ! 
     205      INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     206!!bfr      REAL(wp) ::   zebot, zmshu, zmskv      ! local scalars 
    206207      REAL(wp) ::   zrhoa  = 1.22            ! Air density kg/m3 
    207208      REAL(wp) ::   zcdrag = 1.5e-3          ! drag coefficient 
     
    212213      REAL(wp) ::   zus   , zwlc  , zind     !   -         - 
    213214      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 
     215      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
     216      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc 
     217      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    218218      !!-------------------------------------------------------------------- 
    219219      ! 
    220220      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 )  
    224221      ! 
    225222      zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     
    256253      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    257254      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
     255!!gm old 
    258256!!    DO jj = 2, jpjm1 
    259257!!       DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    263261!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 
    264262!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 
    265 !!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 
     263!!          en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
    266264!!       END DO 
    267265!!    END DO 
     266!!gm new 
     267!! 
     268!!    ====>>>>  add below an wet-only calculation of u and v at t-point like in zdfsh2: 
     269!!                   zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     270!!                   zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     271!! 
     272!! 
     273!!    DO jj = 2, jpjm1 
     274!!       DO ji = fs_2, fs_jpim1   ! vector opt. 
     275!!          zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
     276!!          zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
     277!!          !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
     278!!          zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 
     279!!             &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) 
     280!!          en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 
     281!!       END DO 
     282!!    END DO 
     283!!    IF( ln_isfcav ) THEN       !top friction 
     284!!       DO jj = 2, jpjm1 
     285!!          DO ji = fs_2, fs_jpim1   ! vector opt. 
     286!!             zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
     287!!             zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
     288!!             !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
     289!!             zebot = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 
     290!!                &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
     291!!             en(ji,jj,mikt(ji,jj)+1) = MAX( zebot, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface 
     292!!          END DO 
     293!!       END DO 
     294!!    ENDIF 
     295!! 
    268296!!bfr   - end commented area 
    269297      ! 
    270298      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    271       IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002) 
     299      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002) 
    272300         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    273301         ! 
    274302         !                        !* total energy produce by LC : cumulative sum over jk 
    275          zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1) 
     303         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 
    276304         DO jk = 2, jpk 
    277             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk) 
     305            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 
    278306         END DO 
    279307         !                        !* finite Langmuir Circulation depth 
     
    291319         DO jj = 1, jpj  
    292320            DO ji = 1, jpi 
    293                zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj)) 
     321               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 
    294322            END DO 
    295323         END DO 
     
    300328                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    301329                  !                                           ! 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) ) 
     330                  zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 
     331                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    304332                  !                                           ! TKE Langmuir circulation source term 
    305333                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
     
    318346      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 
    319347      ! 
    320       !                           !* Shear production at uw- and vw-points (energy conserving form) 
    321       CALL zdf_sh2( zsh2 ) 
    322       ! 
    323348      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri ) 
    324349         DO jk = 2, jpkm1 
     
    326351               DO ji = 2, jpim1 
    327352                  !                             ! local Richardson number 
    328                   zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear ) 
     353                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 
    329354                  !                             ! inverse of Prandtl number 
    330355                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
     
    340365               !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 
    341366               !                                   ! 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  )  ) 
     367               zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
     368                  &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  ) 
     369               zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
     370                  &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  ) 
    346371               ! 
    347372               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    350375               ! 
    351376               !                                   ! 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 
     377               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
     378                  &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
    354379                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
    355380                  &                                ) * wmask(ji,jj,jk) 
     
    407432            DO jj = 2, jpjm1 
    408433               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) )   & 
     434                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    410435                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    411436               END DO 
     
    416441            DO ji = fs_2, fs_jpim1   ! vector opt. 
    417442               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) )   & 
     443               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    419444                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    420445            END DO 
     
    429454                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    430455                  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) )   & 
     456                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    432457                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    433458               END DO 
     
    435460         END DO 
    436461      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 )  
    442462      ! 
    443463      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
     
    446466 
    447467 
    448    SUBROUTINE tke_avn 
     468   SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 
    449469      !!---------------------------------------------------------------------- 
    450470      !!                   ***  ROUTINE tke_avn  *** 
     
    480500      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 
    481501      !!---------------------------------------------------------------------- 
     502      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! depth (w-points) 
     503      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w   ! level thickness (t- & w-points) 
     504      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
     505      ! 
    482506      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    483507      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
    484508      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
    485509      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
    486       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld 
     510      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmxlm, zmxld 
    487511      !!-------------------------------------------------------------------- 
    488512      ! 
    489513      IF( nn_timing == 1 )  CALL timing_start('tke_avn') 
    490  
    491       CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
    492514 
    493515      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    502524      ! 
    503525      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
     526         zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    504527         DO jj = 2, jpjm1 
    505528            DO ji = fs_2, fs_jpim1 
    506                zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    507529               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    508530            END DO 
     
    529551      ! 
    530552 !!gm Not sure of that coding for ISF.... 
    531       ! where wmask = 0 set zmxlm == e3w_n 
     553      ! where wmask = 0 set zmxlm == p_e3w 
    532554      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    533555         DO jk = 2, jpkm1 
    534556            DO jj = 2, jpjm1 
    535557               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) ) 
     558                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
     559                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 
    538560                  ! 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)) 
     561                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
     562                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 
    541563               END DO 
    542564            END DO 
     
    547569            DO jj = 2, jpjm1 
    548570               DO ji = fs_2, fs_jpim1   ! vector opt. 
    549                   zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) ) 
     571                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    550572                  zmxlm(ji,jj,jk) = zemxl 
    551573                  zmxld(ji,jj,jk) = zemxl 
     
    558580            DO jj = 2, jpjm1 
    559581               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) ) 
     582                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    561583               END DO 
    562584            END DO 
     
    565587            DO jj = 2, jpjm1 
    566588               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) ) 
     589                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    568590                  zmxlm(ji,jj,jk) = zemxl 
    569591                  zmxld(ji,jj,jk) = zemxl 
     
    576598            DO jj = 2, jpjm1 
    577599               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) ) 
     600                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    579601               END DO 
    580602            END DO 
     
    583605            DO jj = 2, jpjm1 
    584606               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) ) 
     607                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    586608               END DO 
    587609            END DO 
     
    609631               zsqen = SQRT( en(ji,jj,jk) ) 
    610632               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) 
     633               p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     634               p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    613635               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    614636            END DO 
     
    621643            DO jj = 2, jpjm1 
    622644               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) 
     645                  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) 
    624646              END DO 
    625647            END DO 
    626648         END DO 
    627649      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) 
    631650 
    632651      IF(ln_ctl) THEN 
     
    634653         CALL prt_ctl( tab3d_1=avm, clinfo1=' tke  - m: ', ovlap=1, kdim=jpk ) 
    635654      ENDIF 
    636       ! 
    637       CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
    638655      ! 
    639656      IF( nn_timing == 1 )  CALL timing_stop('tke_avn') 
     
    737754      !                               !* set vertical eddy coef. to the background value 
    738755      DO jk = 1, jpk 
    739          avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    740          avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     756         avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 
     757         avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 
    741758      END DO 
    742759      dissl(:,:,:) = 1.e-12_wp 
     
    748765 
    749766   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      ! 
     767      !!--------------------------------------------------------------------- 
     768      !!                   ***  ROUTINE tke_rst  *** 
     769      !!                      
     770      !! ** Purpose :   Read or write TKE file (en) in restart file 
     771      !! 
     772      !! ** Method  :   use of IOM library 
     773      !!                if the restart does not contain TKE, en is either  
     774      !!                set to rn_emin or recomputed  
     775      !!---------------------------------------------------------------------- 
     776      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     777      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     778      ! 
     779      INTEGER ::   jit, jk              ! dummy loop indices 
     780      INTEGER ::   id1, id2, id3, id4   ! local integers 
     781      !!---------------------------------------------------------------------- 
     782      ! 
     783      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     784         !                                   ! --------------- 
     785         IF( ln_rstart ) THEN                   !* Read the restart file 
     786            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. ) 
     787            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 
     788            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 
     789            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 
     790            ! 
     791            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist 
     792               CALL iom_get( numror, jpdom_autoglo, 'en', en ) 
     793               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 
     794               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 
     795               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 
     796            ELSE                                          ! start TKE from rest 
     797               IF(lwp) WRITE(numout,*) '   ==>>   previous run without TKE scheme, set en to background values' 
     798               en(:,:,:) = rn_emin * wmask(:,:,:) 
     799               ! avt_k, avm_k already set to the background value in zdf_phy_init 
     800            ENDIF 
     801         ELSE                                   !* Start from rest 
     802            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set en to the background value' 
     803            en(:,:,:) = rn_emin * wmask(:,:,:) 
     804            ! avt_k, avm_k already set to the background value in zdf_phy_init 
     805         ENDIF 
     806         ! 
     807      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     808         !                                   ! ------------------- 
     809         IF(lwp) WRITE(numout,*) '---- tke-rst ----' 
     810         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
     811         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 
     812         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 
     813         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 
     814         ! 
     815      ENDIF 
     816      ! 
    812817   END SUBROUTINE tke_rst 
    813818 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r7646 r8093  
    6363   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   riceload 
    6464 
    65    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rke          !: kinetic energy 
    66  
    6765   !! arrays relating to embedding ice in the ocean. These arrays need to be declared  
    6866   !! even if no ice model is required. In the no ice model or traditional levitating  
     
    9997         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
    10098         ! 
    101       ALLOCATE(rke(jpi,jpj,jpk)  ,                                         & 
    102          &     sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     & 
    103          &     ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     & 
    104          &     vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     & 
    105          &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       & 
    106          &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     & 
    107          &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     & 
    108          &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     & 
    109          &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     & 
    110          &     riceload(jpi,jpj),                             STAT=ierr(2) ) 
     99      ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     & 
     100         &      ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     & 
     101         &      vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     & 
     102         &      spgu  (jpi,jpj)   , spgv(jpi,jpj)                     ,     & 
     103         &      gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts)                ,     & 
     104         &      gru(jpi,jpj)      , grv(jpi,jpj)                      ,     & 
     105         &      gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts)                ,     & 
     106         &      grui(jpi,jpj)     , grvi(jpi,jpj)                     ,     & 
     107         &      riceload(jpi,jpj)                                     , STAT=ierr(2) ) 
    111108         ! 
    112109      ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) ) 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7954 r8093  
    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 
     
    128131                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD) 
    129132 
     133 
     134!!gm  ===>>>   TO BE REMOVED   (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri) 
     135      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
     136         DO jj = 2, jpjm1 
     137            DO ji = 2, jpim1 
     138               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     139               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
     140            END DO 
     141         END DO 
     142      END DO 
     143      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     144!!gm end 
     145 
     146 
    130147      !  LATERAL  PHYSICS 
    131148      ! 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/SETTE/sette.sh

    r7931 r8093  
    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.