Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/sbcabl.F90	(revision 12154)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/sbcabl.F90	(revision 12179)
@@ -2,6 +2,6 @@
    !!======================================================================
    !!                       ***  MODULE  sbcabl  ***
-   !! Ocean forcing:  momentum, heat and freshwater flux formulation 
-   !!                 derived from an ABL model 
+   !! Ocean forcing:  momentum, heat and freshwater flux formulation
+   !!                 derived from an ABL model
    !!=====================================================================
    !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code
@@ -9,6 +9,6 @@
 
    !!----------------------------------------------------------------------
-   !!   sbc_abl_init  : Initialization of ABL model based on namelist options 
-   !!   sbc_abl       : driver for the computation of momentum, heat and freshwater 
+   !!   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
    !!----------------------------------------------------------------------
@@ -35,5 +35,5 @@
    USE ice    , ONLY : u_ice, v_ice, tm_su, ato_i      ! ato_i = total open water fractional area
    USE sbc_ice, ONLY : wndm_ice, utau_ice, vtau_ice
-#endif   
+#endif
 #if ! defined key_iomput
    USE diawri    , ONLY : dia_wri_alloc_abl
@@ -59,6 +59,6 @@
       !!
       !! ** Purposes :   - read namelist section namsbc_abl
-      !!                 - initialize and check parameter values 
-      !!                 - initialize variables of ABL model 
+      !!                 - initialize and check parameter values
+      !!                 - initialize variables of ABL model
       !!
       !!----------------------------------------------------------------------
@@ -74,5 +74,5 @@
          &                 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, &
-         &                 ln_smth_pblh		 
+         &                 ln_smth_pblh
       !!---------------------------------------------------------------------
 
@@ -87,15 +87,15 @@
       IF(lwm) WRITE( numond, namsbc_abl )
       !
-      ! Check ABL mixing length option      
+      ! 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          
+         &                 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 ' )            
+         &                 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,*)
@@ -106,40 +106,40 @@
             IF(ln_geos_winds) THEN
                ln_geos_winds = .FALSE.
-               WRITE(numout,*) '      ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)'                  
+               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'               
+            WRITE(numout,*) '      ABL -- winds forced by geostrophic winds'
          ELSE
-            WRITE(numout,*) '      ABL -- Geostrophic winds and large-scale pressure gradient are ignored'               
+            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 ( 0 )
+            WRITE(numout,*) '      ABL -- No restoring for ABL winds'
          CASE ( 1 )
-            WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region '             
+            WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region '
          CASE ( 2 )
-            WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere '              
+            WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere '
          END SELECT
          !
-		 IF( ln_smth_pblh )  WRITE(numout,*) '      ABL -- Smoothing of PBL height is activated'   
-		 !
+         IF( ln_smth_pblh )  WRITE(numout,*) '      ABL -- Smoothing of PBL height is activated'
+         !
       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_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 ) 
+      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/) ) ) 
+      jpka   = idimsz(indims - COUNT( (/lluldl/) ) )
       jpkam1 = jpka - 1
 
@@ -152,5 +152,5 @@
 
 #if ! defined key_iomput
-	  IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' )
+      IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' )
 #endif
 
@@ -158,5 +158,5 @@
          WRITE(numout,*)
          WRITE(numout,*) '    sbc_abl_init   : ABL Reference vertical grid'
-         WRITE(numout,*) '    ~~~~~~~'      
+         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 )
@@ -164,6 +164,6 @@
 
       !!---------------------------------------------------------------------
-      !! Check TKE closure parameters 
-      !!--------------------------------------------------------------------- 
+      !! Check TKE closure parameters
+      !!---------------------------------------------------------------------
       rn_Sch  = rn_ce / rn_cm
       mxl_min = (avm_bak / rn_cm) / sqrt( tke_min )
@@ -172,18 +172,18 @@
          WRITE(numout,*)
          WRITE(numout,*) '    abl_zdf_tke   : ABL TKE turbulent closure'
-         WRITE(numout,*) '    ~~~~~~~~~~~'      
+         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                        
+         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
@@ -192,5 +192,5 @@
       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) )  
+         &                   - 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 )
@@ -198,5 +198,5 @@
       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) )       
+         &                   - 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
@@ -229,5 +229,5 @@
          &    + jp_alp1_tra * jp_bmin    + jp_alp0_tra ) * rdt_abl
       zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   &
-         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl 
+         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl
       IF(lwp) THEN
          WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff
@@ -244,10 +244,10 @@
       !!-------------------------------------------------------------------------------------------
       !! 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 
+         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax
          rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8
          !!GS: alternative shape
@@ -260,5 +260,5 @@
       msk_abl(:,:) = tmask(:,:,1)
 
-      !!------------------------------------------------------------------------------------------- 
+      !!-------------------------------------------------------------------------------------------
 
       ! initialize 2D bulk fields AND 3D abl data
@@ -275,14 +275,14 @@
          CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step
 
-          u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:)
-          v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:)
+         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 
+         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
@@ -292,5 +292,5 @@
 
       rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI)
-	  
+
    END SUBROUTINE sbc_abl_init
 
@@ -301,9 +301,9 @@
       !!
       !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
-      !!      the ocean surface from an ABL calculation at each oceanic time step  
-      !!
-      !! ** Method  : 
+      !!      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   
+      !!              - Perform 1 time-step of the ABL model
       !!              - Finalize flux computation in blk_oce_2
       !!
@@ -322,5 +322,5 @@
 #if defined key_si3
       REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi
-#endif  	  
+#endif
       INTEGER                      ::   jbak, jbak_dta, ji, jj
       !!---------------------------------------------------------------------
@@ -329,10 +329,10 @@
       !! 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      ),   &   !   <<= in
@@ -340,6 +340,6 @@
          &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in
          &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in
-         &                zssq, zcd_du, zsen, zevp                          )       !   =>> out
-  
+         &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out
+
 #if defined key_si3
       CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
@@ -347,11 +347,11 @@
          &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
          &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
-#endif  
+#endif
 
       !!-------------------------------------------------------------------------------------------
       !! 3 - Advance ABL variables from now (n) to after (n+1)
-      !!-------------------------------------------------------------------------------------------    
-  
-      CALL abl_stp( kt, sst_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
+      !!-------------------------------------------------------------------------------------------
+
+      CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
          &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
          &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
@@ -360,27 +360,27 @@
          &              zcd_du, zsen, zevp,                                 &   !   <=> in/out
          &              wndm, utau, vtau, taum                              &   !   =>> out
-#if defined key_si3			  
+#if defined key_si3
          &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
          &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
          &            , utau_ice, vtau_ice                                  &   !   =>> out
-#endif			  
+#endif
          &                                                                  )
       !!-------------------------------------------------------------------------------------------
-      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 
+      !! 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( 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                                )
-
-      CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary) 
-      IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file 
+         &            tsk_m, zsen, zevp                                )
+
+      CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)
+      IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file
 
 #if defined key_si3
       ! Avoid a USE abl in icesbc module
       sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
-#endif 
+#endif
 
    END SUBROUTINE sbc_abl
