Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/abl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/abl.F90	(revision 11305)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/abl.F90	(revision 11305)
@@ -0,0 +1,76 @@
+MODULE abl
+   !!======================================================================
+   !!                      ***  MODULE  abl  ***
+   !! Abl        :  ABL dynamics and active tracers defined in memory 
+   !!======================================================================
+   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
+   !!----------------------------------------------------------------------
+   USE par_abl        ! abl parameters
+   USE lib_mpp        ! MPP library
+   USE dom_oce, ONLY: glamt, gphit               ! latitude/longitude
+   USE dom_oce, ONLY: e1t, e1u, e1v, e1f         ! scale factors for horizontal grid
+   USE dom_oce, ONLY: e2t, e2u, e2v, e2f         ! 
+   USE dom_oce, ONLY: rdt                        ! oceanic time-step
+   USE sbc_oce, ONLY: ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka   ! scale factors and altitudes of ABL grid points in the vertical 
+  
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC abl_alloc ! routine called by nemo_init in nemogcm.F90
+
+   !! ABL dynamics and tracer fields                            ! 
+   !! --------------------------                            ! 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   u_abl        !: i-horizontal velocity   [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   v_abl        !: j-horizontal velocity   [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:,:) ::   tq_abl       !: 4D T-q fields           [Kelvin,kg/kg] 
+   
+   !! ABL TKE closure scheme                            ! 
+   !! --------------------------        
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avm_abl      !: turbulent viscosity   [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avt_abl      !: turbulent diffusivity [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxl_abl      !: mixing length         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   tke_abl      !: turbulent kinetic energy [m2/s2]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   fft_abl      !: Coriolis parameter    [1/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   pblh         !: PBL height            [m]
+
+   !! ABL Land/sea mask and restoring term                           ! 
+   !! --------------------------     
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   msk_abl
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   rest_eq
+   !
+   INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2)
+   ! 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION abl_alloc()
+      !!----------------------------------------------------------------------
+      !!                   ***  FUNCTION abl_alloc  ***
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( u_abl  (1:jpi,1:jpj,1:jpka,jptime), &
+         &      v_abl  (1:jpi,1:jpj,1:jpka,jptime), &
+         &      tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), &
+         &      avm_abl(1:jpi,1:jpj,1:jpka), &
+         &      avt_abl(1:jpi,1:jpj,1:jpka), &
+         &      mxl_abl(1:jpi,1:jpj,1:jpka), &         
+         &      tke_abl(1:jpi,1:jpj,1:jpka,jptime), &
+         &      fft_abl(1:jpi,1:jpj), &
+         &      pblh   (1:jpi,1:jpj), &         
+         &      msk_abl(1:jpi,1:jpj), &
+         &      rest_eq(1:jpi,1:jpj), &         
+         &      e3t_abl(1:jpka), e3w_abl(1:jpka), ght_abl(1:jpka), ghw_abl(1:jpka),                 STAT=ierr )
+         !
+      abl_alloc = ierr
+      IF( abl_alloc /= 0 )   CALL ctl_warn('abl_alloc: failed to allocate arrays')
+      !
+   END FUNCTION abl_alloc
+
+   !!======================================================================
+END MODULE abl
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11305)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90	(revision 11305)
@@ -0,0 +1,889 @@
+MODULE ablmod
+   !!======================================================================
+   !!                       ***  MODULE  ablmod  ***
+   !! Surface module :  ABL computation to provide atmospheric data 
+   !!                   for surface fluxes computation
+   !!======================================================================
+   !! History :  3.6  ! 2019-03  (F. Lemarié & G. Samson)  Original code
+   !!----------------------------------------------------------------------
+   
+   !!----------------------------------------------------------------------
+   !!   abl_stp       : ABL single column model
+   !!   abl_zdf_tke   : atmospheric vertical closure
+   !!----------------------------------------------------------------------
+   USE abl            ! ABL dynamics and tracers
+   USE par_abl        ! ABL constants
+
+   USE phycst         ! physical constants
+   USE dom_oce, ONLY  : tmask  
+   USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1
+   USE sbcblk         ! use some physical constants for flux computation
+   !
+   USE prtctl         ! Print control                    (prt_ctl routine)
+   USE iom            ! IOM library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+
+   PUBLIC   abl_stp   ! called by sbcabl.F90
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+
+CONTAINS
+
+
+!===================================================================================================
+   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in
+              &            pu_dta, pv_dta, pt_dta, pq_dta, &
+              &            pslp_dta, pgu_dta, pgv_dta, &
+              &            pcd_du, psen, pevp,     &                            ! in/out
+              &            pwndm, ptaui, ptauj, ptaum )
+!---------------------------------------------------------------------------------------------------
+
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE abl_stp ***
+      !!
+      !! ** Purpose :   Time-integration of the ABL model 
+      !!
+      !! ** Method  :   Compute atmospheric variables : vertical turbulence 
+      !!                             + Coriolis term + newtonian relaxation
+      !!                
+      !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
+      !!              - Advance tracers to time n+1 (Euler backward scheme)
+      !!              - Compute Coriolis term with forward-backward scheme (possibly with geostrophic guide)
+      !!              - Advance u,v to time n+1 (Euler backward scheme)
+      !!              - Apply newtonian relaxation on the dynamics and the tracers
+      !!              - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER  , INTENT(in   )                   ::   kt         ! time step index
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point) 
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pv_dta     ! large-scale windj
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgu_dta    ! large-scale hpgi
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgv_dta    ! large-scale hpgj
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pt_dta     ! large-scale pot. temp.
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pq_dta     ! large-scale humidity
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp_dta   ! sea-level pressure
+      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du     ! Cd x Du (T-point)
+      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
+      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
+      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||      
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau|| 
+      !
+      REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zrhoa, zwnd_i, zwnd_j
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   zFC
+      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF    
+      REAL(wp), DIMENSION(1:jpi,  jptq  )        ::   zBC
+      !
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c
+      !
+      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka )   ::   z_cft  !--FL--to be removed after the test phase
+      !
+      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
+      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1
+      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv
+      LOGICAL             ::   ln_old_coriolis  = .FALSE.              ! possibility to switch off Coriolis term       
+      !
+      !!---------------------------------------------------------------------      
+      !
+      IF(lwp .AND. kt == nit000) THEN                  ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'abl_stp : ABL time stepping'
+         WRITE(numout,*) '~~~~~~'
+      ENDIF 
+      !
+      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
+      !! Compute ustar squared as Cd || Uatm-Uoce ||^2 
+      !! needed for surface boundary condition of TKE 
+      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj)
+         END DO
+      END DO  
+      !
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 
+         
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  2 *** Advance tracers to time n+1
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      
+      !-------------
+      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
+      !-------------    
+         ! Compute matrix elements for interior points 
+         DO jk = 3, jpkam1
+            DO ji = 1, jpi   ! vector opt.
+               z_elem_a( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal
+               z_elem_c( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal       
+               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal           
+            END DO
+         END DO  
+         ! Boundary conditions  
+         DO ji = 1, jpi   ! vector opt. 
+            ! Neumann at the bottom          
+            z_elem_a( ji,     2              ) = 0._wp
+            z_elem_c( ji,     2              ) = - rdt * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                         
+            ! Homogeneous Neumann at the top
+            z_elem_a( ji,     jpka           ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
+            z_elem_c( ji,     jpka           ) = 0._wp
+            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )
+         END DO            
+                   
+         DO jtra = 1,jptq  ! loop on active tracers
+               
+            DO jk = 3, jpkam1
+               DO ji = fs_2, fs_jpim1                 
+                  tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
+               END DO
+            END DO
+
+            IF(jtra == jp_ta) THEN
+               DO ji = 1,jpi  ! surface boundary condition for temperature              
+                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj) 
+                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
+                    &                               + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )                
+                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
+               END DO 
+            ELSE   
+               DO ji = 1,jpi  ! surface boundary condition for humidity              
+                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj) 
+                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   &
+                    &                               + rdt * pevp(ji, jj) * pssq(ji, jj)               
+                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
+               END DO                                     
+            END IF
+            !!
+            !! Matrix inversion
+            !! ----------------------------------------------------------
+            DO ji = 1,jpi            
+               zcff                       =  1._wp / z_elem_b( ji, 2 )
+               zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 )
+               tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra)
+            END DO             
+
+            DO jk = 3, jpka            
+               DO ji = 1,jpi            
+                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
+                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
+                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
+                  &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) )
+               END DO
+            END DO
+!!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...         
+            DO jk = jpkam1,2,-1            
+               DO ji = 1,jpi  
+                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
+                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
+               END DO
+            END DO
+         
+         END DO   !<-- loop on tracers                    
+         !! 
+      !-------------
+      END DO             ! end outer loop 
+      !-------------   
+      
+      
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  3 *** Compute Coriolis term with geostrophic guide
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
+      !-------------        
+      DO jk = 2, jpka    ! outer loop
+      !-------------
+         !             
+         ! Advance u_abl & v_abl to time n+1
+         DO jj = 1, jpj
+            DO ji = 1, jpi           
+               zcff = ( ff_t(ji,jj) * rdt )*( ff_t(ji,jj) * rdt )  ! (f dt)**2
+      
+               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  &
+                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    &
+                  &                 +  rdt * ff_t(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  &
+                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff)
+                  
+               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  &
+                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   &
+                  &                 -  rdt * ff_t(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) &
+                  &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)                
+            END DO
+         END DO 
+         !                                   
+      !-------------
+      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
+      !-------------                
+      !
+      IF( ln_geos_winds ) THEN
+         DO jj = 1, jpj    ! outer loop       
+            DO jk = 1, jpka
+               DO ji = 1, jpi 
+                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
+                     &                      - rdt * e3t_abl(jk) * ff_t(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
+                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
+                     &                      + rdt * e3t_abl(jk) * ff_t(ji, jj  ) * pgu_dta(ji  ,jj  ,jk)
+               END DO
+            END DO
+         END DO      
+      END IF 
+      !-------------            
+      !
+      IF( ln_hpgls_frc ) THEN
+         DO jj = 1, jpj    ! outer loop       
+            DO jk = 1, jpka
+               DO ji = 1, jpi 
+                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)
+                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk)  
+               ENDDO
+            ENDDO
+         ENDDO 
+      END IF
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  4 *** Advance u,v to time n+1 
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+      !
+      !  Vertical diffusion for u_abl     
+      !-------------
+      DO jj = 1, jpj    ! outer loop
+      !-------------    
+
+         DO jk = 3, jpkam1
+            DO ji = 1, jpi  
+               z_elem_a( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal
+               z_elem_c( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                
+               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal
+            END DO
+         END DO  
+            
+         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)            
+            !++ Surface boundary condition
+            z_elem_a( ji,     2    ) = 0._wp
+            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
+            zcff  = rdt * pcd_du(ji, jj)
+            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
+            u_abl( ji, jj,    2, nt_a ) = u_abl( ji, jj,    2, nt_a ) + 0.5_wp * zcff * ( pssu(ji-1, jj) + pssu(ji,jj) ) 
+            !++ Top Neumann B.C.
+            !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 
+            !z_elem_c( ji,     jpka ) = 0._wp
+            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
+            !++ Top Dirichlet B.C.
+            z_elem_a( ji,     jpka )       = 0._wp
+            z_elem_c( ji,     jpka )       = 0._wp
+            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
+            u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)              
+         END DO 
+         !!
+         !! Matrix inversion
+         !! ----------------------------------------------------------
+         DO ji = 2, jpi          
+            zcff                 =   1._wp / z_elem_b( ji, 2 )
+            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
+            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a)
+         END DO             
+
+         DO jk = 3, jpka            
+            DO ji = 2, jpi              
+               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
+               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
+               u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
+               &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
+            END DO
+         END DO
+            
+         DO jk = jpkam1,2,-1            
+            DO ji = 2, jpi
+               u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
+            END DO
+         END DO
+                                          
+      !-------------
+      END DO             ! end outer loop 
+      !-------------   
+
+      !
+      !  Vertical diffusion for v_abl     
+      !-------------
+      DO jj = 2, jpj   ! outer loop
+      !-------------    
+         !
+         DO jk = 3, jpkam1
+            DO ji = 1, jpi   
+               z_elem_a( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
+               z_elem_c( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal              
+               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal
+            END DO
+         END DO
+
+         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)          
+            !++ Surface boundary condition
+            z_elem_a( ji,     2    ) = 0._wp
+            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                       
+            zcff  = rdt * pcd_du(ji, jj)
+            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff 
+            v_abl( ji, jj,    2, nt_a ) =  v_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssv(ji, jj) + pssv(ji, jj-1) )  
+            !++ Top Neumann B.C.            
+            !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
+            !z_elem_c( ji,     jpka ) = 0._wp
+            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                               
+            !++ Top Dirichlet B.C.
+            z_elem_a( ji,     jpka )       = 0._wp
+            z_elem_c( ji,     jpka )       = 0._wp
+            z_elem_b( ji,     jpka )       = e3t_abl( jpka )
+            v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)             
+         END DO 
+         !!
+         !! Matrix inversion
+         !! ----------------------------------------------------------
+         DO ji = 1, jpi              
+            zcff                 =  1._wp / z_elem_b( ji, 2 )
+            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
+            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
+         END DO             
+
+         DO jk = 3, jpka            
+            DO ji = 1, jpi                   
+               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
+               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
+               v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
+               &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
+            END DO
+         END DO
+            
+         DO jk = jpkam1,2,-1            
+            DO ji = 1, jpi    
+               v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
+            END DO
+         END DO
+         !                                          
+      !-------------
+      END DO             ! end outer loop 
+      !-------------   
+  
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  5 *** Apply nudging on the dynamics and the tracers 
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+      z_cft(:,:,:) = 0._wp      
+      
+      IF( nn_dyn_restore > 0  ) THEN
+         !------------- 
+         DO jk = 2, jpka    ! outer loop
+         !-------------       
+            DO jj = 2, jpj
+               DO ji = 2, jpi
+                  zcff1  = pblh( ji, jj )
+                  zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                        
+                  zsig  =            MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
+                  zmsk  = msk_abl(ji,jj)
+                  zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   &
+                     &  + jp_alp1_dyn * zsig    + jp_alp0_dyn
+                  zcff  = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points
+                  
+                  u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   &
+                     &                               + zcff   * pu_dta( ji, jj, jk       )                      
+                  v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   &
+                     &                               + zcff   * pv_dta( ji, jj, jk       )
+                END DO
+            END DO   
+         !-------------
+         END DO             ! end outer loop
+         !-------------               
+      END IF
+
+      !------------- 
+      DO jk = 2, jpka    ! outer loop
+      !-------------       
+         DO jj = 1,jpj
+            DO ji = 1,jpi 
+               zcff1  = pblh( ji, jj )
+               zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
+               zsig   =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
+               zmsk   = msk_abl(ji,jj)
+               zcff2  = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   &
+                  &   + jp_alp1_tra * zsig    + jp_alp0_tra
+               zcff   = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points      
+               tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   &
+                  &                                       + zcff   * pt_dta( ji, jj, jk )
+               
+               tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
+                  &                                       + zcff   * pq_dta( ji, jj, jk )
+               
+               z_cft( ji, jj, jk ) = zcff
+            END DO
+         END DO
+      !-------------
+      END DO             ! end outer loop
+      !------------- 
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  6 *** MPI exchanges  
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), kfillmode = jpfillnothing  )
+      CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing  )   ! ++++ this should not be needed...
+      !
+      CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      ))
+      CALL iom_put ( 'v_abl',  v_abl(:,:,2:jpka,nt_a      ))
+      CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta))
+      CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa))
+!      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           ))
+      CALL iom_put (  'pblh',   pblh(:,:                  ))
+      !
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  7 *** Finalize flux computation
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
+            zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
+            zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ...
+               &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  )
+            psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )
+            pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) )
+            zrhoa( ji, jj )   = zcff              
+         END DO
+      END DO
+      
+      DO jj = 2, jpj
+         DO ji = 2, jpi   ! vect. opt.
+            zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) )  
+            zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) 
+         END DO
+      END DO
+      ! 
+      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. )
+      !
+      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
+      DO jj = 1, jpj
+         DO ji = 1, jpi      
+            zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
+               &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  ) * msk_abl(ji,jj)
+            zztmp         = zrhoa(ji,jj) * pcd_du(ji,jj)
+            
+            pwndm (ji,jj) =         zcff
+            ptaum (ji,jj) = zztmp * zcff
+            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+         END DO
+      END DO   
+      ! ... utau, vtau at U- and V_points, resp.
+      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
+            zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
+            ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
+            zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
+            zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
+            ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
+         END DO
+      END DO
+      !
+      CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1., ptauj(:,:), 'V', -1. )
+
+      CALL iom_put( "taum_oce", ptaum )
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab2d_1=pwndm  , clinfo1=' abl_stp: wndm   : ' )
+         CALL prt_ctl( tab2d_1=ptaui  , clinfo1=' abl_stp: utau   : ' )
+         CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' )
+      ENDIF
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  8 *** Swap time indices for the next timestep
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+      nt_n = 1 + MOD( kt  , 2)
+      nt_a = 1 + MOD( kt+1, 2)
+      !
+!---------------------------------------------------------------------------------------------------
+   END SUBROUTINE abl_stp
+!===================================================================================================
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+!===================================================================================================
+   SUBROUTINE abl_zdf_tke( )
+!---------------------------------------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE abl_zdf_tke  ***
+      !!
+      !! ** Purpose :   Time-step Turbulente Kinetic Energy (TKE) equation
+      !!
+      !! ** Method  : - source term due to shear
+      !!              - source term due to stratification
+      !!              - resolution of the TKE equation by inverting
+      !!                a tridiagonal linear system
+      !!
+      !! ** Action  : - en : now turbulent kinetic energy)
+      !!              - avmu, avmv : production of TKE by shear at u and v-points
+      !!                (= Kz dz[Ub] * dz[Un] )
+      !! ---------------------------------------------------------------------
+      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
+      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
+      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
+      REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2
+      REAL(wp)                                ::   zdU2,zdV2      
+      REAL(wp)                                ::   zwndi,zwndj
+      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
+      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF       
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
+      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
+      LOGICAL                                 ::   ln_Patankar    = .FALSE.
+      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
+      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
+      !
+      tind  = nt_n
+      ziRic = 1._wp / rn_Ric
+      ! if tind = nt_a it is required to apply lbc_lnk on u_abl(nt_a) and v_abl(nt_a)
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  Advance TKE equation to time n+1
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !-------------
+      DO jj = 1, jpj    ! outer loop
+      !-------------
+         !
+         ! Compute vertical shear         
+         DO jk = 2, jpkam1
+            DO ji = 1,jpi  
+               zcff        = 1.0_wp / e3w_abl( jk )**2                
+               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                           
+               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
+               zsh2(ji,jk) = zdU+zdV
+            END DO
+         END DO   
+         !
+         ! Compute brunt-vaisala frequency
+         DO jk = 2, jpkam1
+            DO ji = 1,jpi 
+               zcff  = grav * itvref / e3w_abl( jk )  
+               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
+               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
+                  &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
+               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
+            END DO
+         END DO 
+         !
+         ! Terms for the tridiagonal problem
+         DO jk = 2, jpkam1
+            DO ji = 1,jpi 
+               zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point
+               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation          
+               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
+               
+               z_elem_a( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal
+               z_elem_c( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal          
+               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE
+                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
+                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal       
+                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) )             ! right-hand-side
+               ELSE
+                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   &
+                     &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal    
+                     &                     - e3w_abl(jk) * rdt * zbuoy   
+                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt *  zshear )             ! right-hand-side                     
+               END IF
+            END DO
+         END DO  
+                            
+         DO ji = 1,jpi    ! vector opt.            
+            zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )           
+            zetop   = tke_min             
+            z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                   
+            z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp            
+            tke_abl( ji, jj,    1, nt_a ) = zesrf
+            tke_abl( ji, jj, jpka, nt_a ) = zetop 
+            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
+            zsh2(ji,      1) = zsh2( ji,      2)                                      
+            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
+            zsh2(ji,   jpka) = zsh2( ji  , jpkam1)
+         END DO        
+         !!
+         !! Matrix inversion
+         !! ----------------------------------------------------------
+         DO ji = 1,jpi
+            zcff                  =  1._wp / z_elem_b( ji, 1 )
+            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
+            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
+         END DO             
+
+         DO jk = 2, jpka            
+            DO ji = 1,jpi
+               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
+               zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
+               tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
+               &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
+            END DO
+         END DO
+            
+         DO jk = jpkam1,1,-1            
+            DO ji = 1,jpi
+               tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
+            END DO
+         END DO
+         
+!!FL should not be needed because of Patankar procedure  
+         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
+
+         !!
+         !! Diagnose PBL height
+         !! ----------------------------------------------------------
+         
+         
+         !                                                       
+         ! arrays zRH, zFC and zCF are available at this point
+         ! and zFC(:, 1 ) = 0.
+         ! diagnose PBL height based on zsh2 and zbn2
+         zFC (  :  ,1) = 0._wp
+         ikbl( 1:jpi ) = 0 
+         
+         DO jk = 2,jpka
+            DO ji = 1, jpi 
+               zcff  = ghw_abl( jk-1 )
+               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
+               zcff  = ghw_abl( jk   )
+               zFC( ji, jk ) = zFC( ji, jk-1) + 0.5_wp * e3t_abl( jk )*(                 &
+                               zcff2 * ( zsh2( ji, jk  ) - ziRic * zbn2( ji, jj, jk   ) &
+                           - rn_Cek  * ( ff_t( ji, jj  ) * ff_t( ji, jj ) ) ) &
+                             + zcff1 * ( zsh2( ji, jk-1) - ziRic * zbn2( ji, jj, jk-1 ) &
+                           - rn_Cek  * ( ff_t( ji, jj  ) * ff_t( ji, jj ) ) ) &
+                           &                                                 )
+               IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
+            END DO
+         END DO
+         !
+         ! finalize the computation of the PBL height
+         DO ji = 1, jpi
+            jk = ikbl(ji)
+            IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
+               pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
+                  &              -  ghw_abl( jk   ) * zFC( ji, jk-1 )       &
+                  &              ) / ( zFC( ji, jk   ) - zFC( ji, jk-1 ) )
+            ELSE IF( jk==2 ) THEN
+               pblh( ji, jj ) = ghw_abl(2   )
+            ELSE
+               pblh( ji, jj ) = ghw_abl(jpka)
+            END IF 
+         END DO
+         ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
+         !
+      !-------------      
+      END DO      
+      !------------- 
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  Diagnostic mixing length computation
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      SELECT CASE ( nn_amxl )
+      !
+      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom         
+#   define zlup zRH      
+#   define zldw zFC 
+         DO jj = 1, jpj     ! outer loop
+            !
+            DO ji = 1, jpi
+               mxl_abl ( ji, jj,    1 )  = 0._wp
+               mxl_abl ( ji, jj, jpka )  = mxl_min
+               zldw( ji,        1 )  = 0._wp              
+               zlup( ji,     jpka )  = 0._wp
+            END DO
+            !
+            DO jk = 2, jpkam1   
+               DO ji = 1, jpi  
+                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )  
+                  mxl_abl( ji, jj, jk ) = MAX( mxl_min,  &
+                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
+               END DO
+            END DO   
+            !
+            ! Limit mxl
+            DO jk = jpkam1,1,-1   
+               DO ji = 1, jpi 
+                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) )
+               END DO
+            END DO   
+            !
+            DO jk = 2, jpka
+               DO ji = 1, jpi  
+                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )   
+               END DO
+            END DO 
+            !
+            DO jk = 1, jpka
+               DO ji = 1, jpi
+                  mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )   
+               END DO
+            END DO                         
+            !
+         END DO
+#   undef zlup      
+#   undef zldw  
+         !
+         !
+      CASE ( 1 )             ! length-scale computed as the distance to the PBL height
+         DO jj = 1,jpj      ! outer loop
+            !
+            DO ji = 1, jpi   ! vector opt.
+               zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl              
+               DO jk = 1, jpka              
+                  zsig  = MIN( zcff * ghw_abl( jk ), 1. )    
+                  zcff1 = pblh( ji, jj )                 
+                  mxl_abl( ji, jj, jk ) =  mxl_min                           &
+                     &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) &
+                     &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) &
+                     &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) &
+                     &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) &
+                     &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min )
+               END DO
+            END DO 
+            !            
+         END DO            
+         !
+      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale
+         !
+#   define zlup zRH      
+#   define zldw zFC           
+! zCF is used for matrix inversion
+!             
+       DO jj = 1, jpj      ! outer loop
+         
+         DO ji = 1, jpi           
+            zlup( ji,    1 ) = mxl_min
+            zldw( ji,    1 ) = mxl_min               
+            zlup( ji, jpka ) = mxl_min
+            zldw( ji, jpka ) = mxl_min                
+         END DO            
+         
+         DO jk = 2,jpka-1
+            DO ji = 1, jpi
+               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
+               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)        
+            END DO
+         END DO         
+         !!
+         !! BL89 search for lup
+         !! ----------------------------------------------------------           
+         DO jk=2,jpka-1 
+            !
+            DO ji = 1, jpi
+               zCF(ji,1:jpka) = 0._wp
+               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
+               ln_foundl(ji ) = .false.
+            END DO   
+            !           
+            DO jkup=jk+1,jpka-1
+               DO ji = 1, jpi
+                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
+                     &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) &
+                     &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )
+                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
+                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
+                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                   
+                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
+                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
+                     zlup(ji,jk) = zcff                
+                     ln_foundl(ji) = .true.
+                  END IF
+               END DO
+            END DO
+            !
+         END DO   
+         !!
+         !! BL89 search for ldwn
+         !! ----------------------------------------------------------          
+         DO jk=2,jpka-1         
+            !
+            DO ji = 1, jpi
+               zCF(ji,1:jpka) = 0._wp
+               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
+               ln_foundl(ji ) = .false.
+            END DO  
+            !   
+            DO jkdwn=jk-1,1,-1
+               DO ji = 1, jpi             
+                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
+                     &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
+                                       + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )   
+                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
+                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
+                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )              
+                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
+                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
+                     zldw(ji,jk) = zcff 
+                     ln_foundl(ji) = .true.               
+                  END IF                                   
+               END DO
+            END DO
+            !      
+         END DO
+
+         DO jk = 1, jpka
+            DO ji = 1, jpi  
+               mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )  
+            END DO
+         END DO  
+
+      END DO
+#   undef zlup      
+#   undef zldw           
+         !
+      END SELECT      
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  Finalize the computation of turbulent visc./diff.
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      
+      !-------------
+      DO jj = 1, jpj     ! outer loop
+      !-------------
+         DO jk = 1, jpka   
+            DO ji = 1, jpi  ! vector opt.
+               zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  &
+                  &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) ) 
+               zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z)
+               zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
+!!FL: MAX function probably useless because of the definition of mxl_min             
+               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
+               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                              
+            END DO
+         END DO
+      !-------------         
+      END DO      
+      !-------------
+
+!---------------------------------------------------------------------------------------------------
+   END SUBROUTINE abl_zdf_tke
+!===================================================================================================
+
+
+!!======================================================================
+END MODULE ablmod
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90	(revision 11305)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90	(revision 11305)
@@ -0,0 +1,90 @@
+MODULE par_abl
+   !!======================================================================
+   !!                        ***  par_abl  ***
+   !! ABL :   set the Atmospheric Boundary Layer parameters
+   !!======================================================================
+   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
+   !!----------------------------------------------------------------------
+   USE par_kind          ! kind parameters
+   USE par_oce           ! to access horizontal domain size & mpp decomposition
+
+   IMPLICIT NONE
+   PUBLIC
+
+   !!---------------------------------------------------------------------
+   !! Active tracer parameters
+   !!---------------------------------------------------------------------
+   INTEGER , PUBLIC, PARAMETER ::   jptq   = 2     !: Number of active tracers (=2, i.e. T & q )
+   INTEGER , PUBLIC, PARAMETER ::   jp_ta  = 1     !: indice for temperature
+   INTEGER , PUBLIC, PARAMETER ::   jp_qa  = 2     !: indice for humidity
+   INTEGER , PUBLIC, PARAMETER ::   jptime = 2     !: number of time indices stored in memory
+   !!---------------------------------------------------------------------
+   !! namABL Namelist options  
+   !!---------------------------------------------------------------------   
+   INTEGER , PUBLIC            ::   nn_amxl        !: mixing length option 
+   INTEGER , PUBLIC            ::   nn_dyn_restore !: restoring option for dynamical ABL variables
+   LOGICAL , PUBLIC            ::   ln_geos_winds  !: large-scale restoring of ABL winds toward geostrophic winds
+   LOGICAL , PUBLIC            ::   ln_hpgls_frc   !: forcing of ABL winds by large-scale pressure gradient 
+
+   !!---------------------------------------------------------------------
+   !! ABL parameters for TKE turbulent closure
+   !!---------------------------------------------------------------------
+   REAL(wp), PUBLIC, PARAMETER ::   tke_min   = 1.e-6_wp          !: minimum TKE
+   REAL(wp), PUBLIC, PARAMETER ::   avm_bak   = 1.e-4_wp          !: background viscosity
+   REAL(wp), PUBLIC, PARAMETER ::   avt_bak   = 1.e-5_wp          !: background diffusion
+   REAL(wp), PUBLIC, PARAMETER ::   itvref    = 1.0_wp / 288.0_wp !: inverse of reference virtual temperature     
+   !++ TKE closure parameters
+   REAL(wp), PUBLIC, PARAMETER ::   rn_phimax = (1._wp - 2.2_wp) / 2.2_wp !: maximum value for Ri * mxl^2 * N^2 / tke in phiz computation
+   REAL(wp), PUBLIC, PARAMETER ::   rn_Cek    = 258._wp                   !: Ekman constant for Richardson number 
+   REAL(wp), PUBLIC, PARAMETER ::   rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp )
+   REAL(wp), PUBLIC            ::   rn_ceps                       !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_cm                         !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_ct                         !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_ce                         !: namelist parameter 
+   REAL(wp), PUBLIC            ::   rn_Rod                        !: namelist parameter   
+   REAL(wp), PUBLIC            ::   rn_Sch    
+   REAL(wp), PUBLIC            ::   mxl_min    
+   REAL(wp), PUBLIC            ::   rn_ldyn_min                   !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_ldyn_max                   !: namelist parameter  
+   REAL(wp), PUBLIC            ::   rn_ltra_min                   !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_ltra_max                   !: namelist parameter
+   REAL(wp), PUBLIC            ::   rn_Ric                        !: critical Richardson number 
+
+   !!---------------------------------------------------------------------
+   !! ABL parameters for the vertical profile of the restoring term 
+   !!---------------------------------------------------------------------
+   REAL(wp), PUBLIC, PARAMETER ::   jp_bmin    =   0.5_wp 
+   REAL(wp), PUBLIC, PARAMETER ::   jp_bmax    =   1.5_wp
+   REAL(wp), PUBLIC            ::   jp_alp0_tra
+   REAL(wp), PUBLIC            ::   jp_alp1_tra   
+   REAL(wp), PUBLIC            ::   jp_alp2_tra   
+   REAL(wp), PUBLIC            ::   jp_alp3_tra   
+   REAL(wp), PUBLIC            ::   jp_alp0_dyn
+   REAL(wp), PUBLIC            ::   jp_alp1_dyn   
+   REAL(wp), PUBLIC            ::   jp_alp2_dyn   
+   REAL(wp), PUBLIC            ::   jp_alp3_dyn  
+   REAL(wp), PUBLIC            ::   jp_pblh_min
+   REAL(wp), PUBLIC            ::   jp_pblh_max      
+   ! parameter for the semi-implicit treatment of Coriolis term  
+   REAL(wp), PUBLIC, PARAMETER ::   gamma_Cor  = 0.55_wp
+
+   !!---------------------------------------------------------------------
+   !! ABL parameters for the diagnostic mixing length option nn_amxl = 1 
+   !!---------------------------------------------------------------------   
+   REAL(wp), PUBLIC, PARAMETER ::   amx1 =    4.3995604393911742_wp
+   REAL(wp), PUBLIC, PARAMETER ::   amx2 =  -18.159100102732943_wp
+   REAL(wp), PUBLIC, PARAMETER ::   amx3 =   40.241226956967239_wp
+   REAL(wp), PUBLIC, PARAMETER ::   amx4 =  -43.603409583363678_wp
+   REAL(wp), PUBLIC, PARAMETER ::   amx5 =   17.121715347554314_wp
+   REAL(wp), PUBLIC, PARAMETER ::   bmx1 =  -16.262675447730114_wp
+   REAL(wp), PUBLIC, PARAMETER ::   bmx2 =   85.088728134110781_wp
+   REAL(wp), PUBLIC, PARAMETER ::   bmx3 = -193.46548261141191_wp
+   REAL(wp), PUBLIC, PARAMETER ::   bmx4 =  196.71548261141191_wp
+   REAL(wp), PUBLIC, PARAMETER ::   bmx5 =  -72.076052686380677_wp 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/ABL 4.0 , NEMO Consortium (2018)
+   !! $Id: sbc_oce.F90 10882 2019-04-17 15:40:17Z clem $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!======================================================================
+END MODULE par_abl
Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11305)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90	(revision 11305)
@@ -0,0 +1,343 @@
+MODULE sbcabl
+   !!======================================================================
+   !!                       ***  MODULE  sbcabl  ***
+   !! Ocean forcing:  momentum, heat and freshwater flux formulation 
+   !!                 derived from an ABL model 
+   !!=====================================================================
+   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_abl_init  : Initialization of ABL model based on namelist options 
+   !!   sbc_abl       : driver for the computation of momentum, heat and freshwater 
+   !!                   fluxes over ocean via the ABL model
+   !!----------------------------------------------------------------------
+   USE abl            ! ABL
+   USE par_abl        ! abl parameters
+   USE ablmod
+
+   USE phycst         ! physical constants
+   USE fldread        ! read input fields
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbcblk         ! Surface boundary condition: bulk formulae
+   USE dom_oce, ONLY  : tmask
+   !
+   USE iom            ! I/O manager library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distribued memory computing library
+   USE lib_fortran    ! to use key_nosignedzero
+   USE timing         ! Timing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_abl_init       ! routine called in sbcmod module
+   PUBLIC   sbc_abl            ! routine called in sbcmod module
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO-consortium (2014)
+   !! $Id: sbcabl.F90 6416 2016-04-01 12:22:17Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_abl_init
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_abl_init  ***
+      !!
+      !! ** Purposes :   - read namelist section namsbc_abl
+      !!                 - initialize and check parameter values 
+      !!                 - initialize variables of ABL model 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER            ::   ji, jj, jk, jbak, jbak_dta       ! dummy loop indices
+      INTEGER            ::   ios, ierror, ioptio              ! Local integer
+      INTEGER            ::   inum, indims, idimsz(4), id
+      CHARACTER(len=100) ::   cn_dir, cn_dom                   ! Atmospheric grid directory
+      REAL(wp)           ::   zcff,zcff1
+      LOGICAL            ::   lluldl
+      NAMELIST/namsbc_abl/ cn_dir , cn_dom, ln_hpgls_frc, ln_geos_winds,          &
+         &                 nn_dyn_restore,                                        &
+         &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   &
+         &                 nn_amxl, rn_cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric 
+      !!---------------------------------------------------------------------
+
+      REWIND( numnam_ref )              ! Namelist namsbc_abl in reference namelist : ABL parameters
+      READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namsbc_abl in configuration namelist : ABL parameters
+      READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist', lwp )
+      !
+      IF(lwm) WRITE( numond, namsbc_abl )
+      !
+      ! Check ABL mixing length option      
+      IF( nn_amxl  < 0   .OR.  nn_amxl  > 2 )   &
+         &                 CALL ctl_stop( 'abl_init : bad flag, nn_amxl must be  0, 1 or 2 ' )         
+      !
+      ! Check ABL dyn restore option          
+      IF( nn_dyn_restore  < 0   .OR.  nn_dyn_restore  > 2 )   &
+         &                 CALL ctl_stop( 'abl_init : bad flag, nn_dyn_restore must be  0, 1 or 2 ' )            
+      !
+      !!---------------------------------------------------------------------
+      !! Control prints
+      !!--------------------------------------------------------------------- 
+      IF(lwp) THEN                              ! Control print (other namelist variable)
+         WRITE(numout,*)
+         WRITE(numout,*) '      ABL -- cn_dir         = ', cn_dir
+         WRITE(numout,*) '      ABL -- cn_dom         = ', cn_dom
+         IF( ln_hpgls_frc ) THEN
+            WRITE(numout,*) '      ABL -- winds forced by large-scale pressure gradient'
+            IF(ln_geos_winds) THEN
+               ln_geos_winds = .FALSE.
+               WRITE(numout,*) '      ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)'                  
+            END IF
+         ELSE IF( ln_geos_winds ) THEN
+            WRITE(numout,*) '      ABL -- winds forced by geostrophic winds'               
+         ELSE
+            WRITE(numout,*) '      ABL -- Geostrophic winds and large-scale pressure gradient are ignored'               
+         END IF
+         !
+         SELECT CASE ( nn_dyn_restore )
+         CASE ( 0 )   
+            WRITE(numout,*) '      ABL -- No restoring for ABL winds'   
+         CASE ( 1 )
+            WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region '             
+         CASE ( 2 )
+            WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere '              
+         END SELECT
+         !
+      ENDIF
+
+      !!---------------------------------------------------------------------
+      !! Convert nudging coefficient from hours to 1/sec
+      !!---------------------------------------------------------------------         
+      zcff        = 1._wp / 3600._wp
+      rn_ldyn_min = zcff / rn_ldyn_min
+      rn_ldyn_max = zcff / rn_ldyn_max      
+      rn_ltra_min = zcff / rn_ltra_min
+      rn_ltra_max = zcff / rn_ltra_max   
+
+      !!---------------------------------------------------------------------
+      !! ABL grid initialization 
+      !!---------------------------------------------------------------------      
+      CALL iom_open( TRIM(cn_dir)//TRIM(cn_dom), inum ) 
+      id     = iom_varid( inum, 'e3t_abl', kdimsz=idimsz, kndims=indims, lduld=lluldl )
+      jpka   = idimsz(indims - COUNT( (/lluldl/) ) ) 
+      jpkam1 = jpka - 1
+
+      IF( abl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' )
+      CALL iom_get( inum, jpdom_unknown, 'e3t_abl', e3t_abl(:) )
+      CALL iom_get( inum, jpdom_unknown, 'e3w_abl', e3w_abl(:) )
+      CALL iom_get( inum, jpdom_unknown, 'ght_abl', ght_abl(:) )
+      CALL iom_get( inum, jpdom_unknown, 'ghw_abl', ghw_abl(:) )
+      CALL iom_close( inum )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '    sbc_abl_init   : ABL Reference vertical grid'
+         WRITE(numout,*) '    ~~~~~~~'      
+         WRITE(numout, "(9x,'  level ght_abl   ghw_abl   e3t_abl   e3w_abl  ')" )
+         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, ght_abl(jk), ghw_abl(jk), e3t_abl(jk), e3w_abl(jk), jk = 1, jpka )
+      END IF
+
+      !!---------------------------------------------------------------------
+      !! Check TKE closure parameters 
+      !!--------------------------------------------------------------------- 
+      rn_Sch  = rn_ce / rn_cm
+      mxl_min = (avm_bak / rn_cm) / sqrt( tke_min )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '    abl_zdf_tke   : ABL TKE turbulent closure'
+         WRITE(numout,*) '    ~~~~~~~~~~~'      
+         IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale '
+         IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height '
+         WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2'
+         WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m'
+         WRITE(numout,*) ' Constant for turbulent viscosity           = ',rn_Cm         
+         WRITE(numout,*) ' Constant for turbulent diffusivity         = ',rn_Ct         
+         WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch            
+         WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps                        
+      END IF
+
+      !!-------------------------------------------------------------------------------------------
+      !! Compute parameters to build the vertical profile for the nudging term (used in abl_stp())
+      !!-------------------------------------------------------------------------------------------       
+      zcff1       = 1._wp / ( jp_bmax - jp_bmin )**3
+      ! for active tracers
+      jp_alp3_tra = -2._wp * zcff1                       * ( rn_ltra_max - rn_ltra_min )
+      jp_alp2_tra =  3._wp * zcff1 * (jp_bmax + jp_bmin) * ( rn_ltra_max - rn_ltra_min )
+      jp_alp1_tra = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ltra_max - rn_ltra_min )
+      jp_alp0_tra = zcff1 * (  rn_ltra_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      &
+         &                   - rn_ltra_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) )  
+      ! for dynamics
+      jp_alp3_dyn = -2._wp * zcff1                       * ( rn_ldyn_max - rn_ldyn_min )
+      jp_alp2_dyn =  3._wp * zcff1 * (jp_bmax + jp_bmin) * ( rn_ldyn_max - rn_ldyn_min )
+      jp_alp1_dyn = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ldyn_max - rn_ldyn_min )
+      jp_alp0_dyn = zcff1 * (  rn_ldyn_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      &
+         &                   - rn_ldyn_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) )       
+
+      jp_pblh_min = ghw_abl(     4) / jp_bmin  !<-- at least 3 grid points at the bottom have value rn_ltra_min
+      jp_pblh_max = ghw_abl(jpka-3) / jp_bmax  !<-- at least 3 grid points at the top    have value rn_ltra_max
+
+      ! Check parameters for dynamics
+      zcff  = jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   &
+         &  + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn 
+      zcff1 = jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   &
+         &  + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn      
+      IF(lwp) THEN
+         IF(nn_dyn_restore > 0) THEN
+            WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff  * rdt
+            WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 * rdt
+            ! Check that restoring coefficients are between 0 and 1
+            IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
+               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
+            IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
+               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' )
+            IF( zcff  * rdt > zcff1 * rdt ) &
+               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' )
+         END IF
+      END IF
+
+      ! Check parameters for active tracers
+      zcff  = jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   &
+         &  + jp_alp1_tra * jp_bmin    + jp_alp0_tra 
+      zcff1 = jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   &
+         &  + jp_alp1_tra * jp_bmax    + jp_alp0_tra  
+      IF(lwp) THEN
+         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff  * rdt
+         WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 * rdt
+         ! Check that restoring coefficients are between 0 and 1
+         IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   &
+            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
+         IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   &
+            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' )
+         IF( zcff  * rdt > zcff1  * rdt ) &
+            &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' )
+      END IF
+
+      !!-------------------------------------------------------------------------------------------
+      !! Initialize Coriolis frequency, equatorial restoring and land/sea mask
+      !!-------------------------------------------------------------------------------------------  
+      fft_abl(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
+
+      ! Equatorial restoring
+      IF( nn_dyn_restore == 1 ) THEN
+         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax 
+         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff)/zcff ) )**8
+      ELSE
+         rest_eq(:,:) = 1._wp
+      END IF
+      ! T-mask
+      msk_abl(:,:) = tmask(:,:,1)
+
+      !!------------------------------------------------------------------------------------------- 
+
+      ! initialize 2D bulk fields AND 3D abl data
+      CALL sbc_blk_init
+
+      ! initialize ABL from data or restart
+      !!GS  disabled for now
+      !IF( ln_rstart ) THEN
+      !  CALL ctl_stop( 'STOP', 'sbc_abl_init: restart mode not supported yet' )
+      !ELSE
+
+      CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step
+
+      ! Initialize the time index for now time (nt_n) and after time (nt_a)
+      nt_n = 1 + MOD( nit000  , 2)
+      nt_a = 1 + MOD( nit000+1, 2)
+
+       u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:)
+       v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:)
+      tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:)
+      tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:)
+
+      tke_abl(:,:,:,nt_n     ) = tke_min
+      avm_abl(:,:,:          ) = avm_bak
+      avt_abl(:,:,:          ) = avt_bak
+      mxl_abl(:,:,:          ) = mxl_min 
+      pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points 
+      u_abl  (:,:,:,nt_a     ) = 0._wp
+      v_abl  (:,:,:,nt_a     ) = 0._wp
+      tq_abl (:,:,:,nt_a,:   ) = 0._wp
+      tke_abl(:,:,:,nt_a     ) = 0._wp
+      !ENDIF
+      !!GS restart case not supported
+
+   END SUBROUTINE sbc_abl_init
+
+  
+
+   SUBROUTINE sbc_abl( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_abl  ***
+      !!
+      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
+      !!      the ocean surface from an ABL calculation at each oceanic time step  
+      !!
+      !! ** Method  : 
+      !!              - Pre-compute part of turbulent fluxes in blk_oce_1
+      !!              - Perform 1 time-step of the ABL model   
+      !!              - Finalize flux computation in blk_oce_2
+      !!
+      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
+      !!              - vtau    : j-component of the stress at V-point  (N/m2)
+      !!              - taum    : Wind stress module at T-point         (N/m2)
+      !!              - wndm    : Wind speed module at T-point          (m/s)
+      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
+      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
+      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER ,         INTENT(in) ::   kt   ! ocean time step
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
+      INTEGER                      ::   jbak, jbak_dta, ji, jj
+      !!---------------------------------------------------------------------
+      !
+      !!-------------------------------------------------------------------------------------------
+      !! 1 - Read Atmospheric 3D data for large-scale forcing
+      !!-------------------------------------------------------------------------------------------
+      
+      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
+     
+      !!-------------------------------------------------------------------------------------------
+      !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
+      !!-------------------------------------------------------------------------------------------  
+         
+      CALL blk_oce_1( kt, u_abl(:,:,2,nt_n), v_abl(:,:,2,nt_n), tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),             &
+              &           sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, &
+              &              zssq, zcd_du, zsen, zevp )
+  
+      !!-------------------------------------------------------------------------------------------
+      !! 3 - Advance ABL variables from now (n) to after (n+1)
+      !!-------------------------------------------------------------------------------------------    
+  
+      CALL abl_stp( kt, sst_m, ssu_m, ssv_m, zssq, &                            ! in
+              &         sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &
+              &         sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), sf(jp_slp)%fnow(:,:,1),  &
+              &         sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),                          &
+              &         zcd_du, zsen, zevp,    &                                ! in/out
+              &         wndm, utau, vtau, taum )                                ! out
+
+      !!-------------------------------------------------------------------------------------------
+      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 
+      !!                                                                time swap is done in abl_stp
+      !!-------------------------------------------------------------------------------------------   
+
+      CALL blk_oce_2( kt, tq_abl(:,:,2,nt_n,jp_ta), sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),     &
+              &           sf(jp_prec)%fnow(:,:,1),  sf(jp_snow)%fnow(:,:,1),   &
+              &           sst_m , zsen, zevp )
+
+   END SUBROUTINE sbc_abl
+
+
+   !!======================================================================
+END MODULE sbcabl
