Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2452)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2497)
@@ -28,5 +28,4 @@
    USE in_out_manager ! I/O manager
    USE iom            ! I/O manager library
-   USE zdfbfr, ONLY : rn_hbro, wbotu, wbotv ! bottom roughness and bottom stresses
 
    IMPLICIT NONE
@@ -61,4 +60,5 @@
 
    REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness
+   REAL(wp) ::   hbro          =  0.003_wp    ! Bottom roughness (m)
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
    REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
@@ -143,21 +143,21 @@
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
-      ! Compute wind stress at T-points
+      ! Compute surface and bottom friction at T-points
 !CDIR NOVERRCHK
       DO jj = 2, jpjm1
 !CDIR NOVERRCHK
-        DO ji = fs_2, fs_jpim1   ! vector opt.
-          ! 
-          ! wind stress
-          ! squared surface velocity scale
-          ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1)
-          !
-          ! bottom friction 
-          ztx2 = ( wbotu(ji-1,jj)*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1) )   &
-            &  / MAX( 1.,         umask(ji-1,jj,1) +              umask(ji,jj,1) )
-          zty2 = ( wbotv(ji,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1) )   &
-            &  / MAX( 1.,         vmask(ji,jj-1,1) +              vmask(ji,jj,1) )
-          ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
-        END DO
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ! 
+            ! surface friction 
+            ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1)
+            !
+            ! bottom friction (explicit before friction)
+            ! Note that we chose here not to bound the friction as in dynbfr)
+            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
+               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
+            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
+               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
+            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+         END DO
       END DO  
 
@@ -617,5 +617,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro
+         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro
          !                      ! Balance between the production and the dissipation terms
 !CDIR NOVERRCHK
@@ -625,5 +625,5 @@
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
-               zdep(ji,jj) = vkarmn * rn_hbro
+               zdep(ji,jj) = vkarmn * hbro
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                z_elem_a(ji,jj,ibot) = 0._wp
@@ -632,5 +632,5 @@
                !
                ! Just above last level, Dirichlet condition again (GOTM like)
-               zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1))
+               zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) )
                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
                z_elem_a(ji,jj,ibotm1) = 0._wp
@@ -650,5 +650,5 @@
                !
                ! Bottom level Dirichlet condition:
-               zdep(ji,jj) = vkarmn * rn_hbro
+               zdep(ji,jj) = vkarmn * hbro
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                !
@@ -662,5 +662,5 @@
                !
                ! Set psi vertical flux at the bottom:
-               zdep(ji,jj) = rn_hbro + 0.5_wp*fse3t(ji,jj,ibotm1)
+               zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1)
                zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
@@ -906,6 +906,4 @@
          WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin
          WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin
-         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro
-         WRITE(numout,*) '      Bottom roughness (m)                          rn_hbro       = ', rn_hbro
          WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim
          WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp
@@ -920,5 +918,7 @@
          WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
-         WRITE(numout,*)
+         WRITE(numout,*) '   Hard coded parameters'
+         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro
+         WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro
       ENDIF
 
@@ -1195,6 +1195,4 @@
             id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
             id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
-            id7 = iom_varid( numror, 'wbotu', ldstop = .FALSE. )
-            id8 = iom_varid( numror, 'wbotv', ldstop = .FALSE. )
             !
             IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN        ! all required arrays exist
@@ -1205,40 +1203,14 @@
                CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
                CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
-               CALL iom_get( numror, jpdom_autoglo, 'wbotu' , wbotu  )
-               CALL iom_get( numror, jpdom_autoglo, 'wbotv' , wbotv  )
             ELSE                        
                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
-               IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' 
                en  (:,:,:) = rn_emin
                mxln(:,:,:) = 0.001        
-               ! Initialize bottom stresses
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ikbu = mbku(ji,jj)         ! bottom ocean level of u-point
-                     ikbv = mbkv(ji,jj)
-                     cbx  = avmu(ji,jj,ikbu+1) / fse3uw(ji,jj,ikbu+1)
-                     cby  = avmv(ji,jj,ikbv+1) / fse3vw(ji,jj,ikbv+1)
-                     wbotu(ji,jj) = -cbx * un(ji,jj,ikbu) * umask(ji,jj,1)
-                     wbotv(ji,jj) = -cby * vn(ji,jj,ikbv) * vmask(ji,jj,1)
-                  END DO
-               END DO
                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
             ENDIF
          ELSE                                   !* Start from rest
             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
-            IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated'
             en  (:,:,:) = rn_emin
             mxln(:,:,:) = 0.001       
-            ! Initialize bottom stresses
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ikbu = mbku(ji,jj)            ! bottom ocean level of u-point
-                  ikbv = mbkv(ji,jj)
-                  cbx  = avmu(ji,jj,ikbu+1) / fse3uw(ji,jj,ikbu+1)
-                  cby  = avmv(ji,jj,ikbv+1) / fse3vw(ji,jj,ikbv+1)
-                  wbotu(ji,jj) = -cbx * un(ji,jj,ikbu) * umask(ji,jj,1)
-                  wbotv(ji,jj) = -cby * vn(ji,jj,ikbv) * vmask(ji,jj,1)
-               END DO
-            END DO
          ENDIF
          !
@@ -1252,6 +1224,4 @@
          CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
          CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln  )
-         CALL iom_rstput( kt, nitrst, numrow, 'wbotu' , wbotu )
-         CALL iom_rstput( kt, nitrst, numrow, 'wbotv' , wbotv )
          !
       ENDIF
