Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5676)
@@ -36,5 +36,8 @@
          <field id="sstgrad2"     long_name="square of module of sst gradient"                                                        unit="degC2/m2" />
          <field id="sbt"          long_name="sea bottom temperature"                                                                  unit="degC"     />
-	 
+         <field id="sst_wl"       long_name="Delta SST of warm layer"                                                                 unit="degC"     />
+         <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                  unit="degC"     />
+	 <field id="temp_3m"      long_name="temperature at 3m"                                                                       unit="degC"     />
+         
          <field id="sss"          long_name="sea surface salinity"             standard_name="sea_surface_salinity"   unit="1e-3" />
          <field id="sss2"         long_name="square of sea surface salinity"                                          unit="1e-6"  > sss * sss </field >
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5676)
@@ -1293,2 +1293,8 @@
    rn_htrmax         =  200.0   ! max. depth of transition range
 /
+!-----------------------------------------------------------------------
+&namdiu !   Cool skin and warm layer modesl
+!-----------------------------------------------------------------------
+   ln_diurnal = .false.
+   ln_diurnal_only = .false.
+/
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 5676)
@@ -55,4 +55,6 @@
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
+   USE diurnal_bulk    ! diurnal warm layer
+   USE cool_skin       ! Cool skin
    USE wrk_nemo        ! working array
 
@@ -145,4 +147,11 @@
       ENDIF
 
+      IF(ln_diurnal) THEN
+         CALL iom_put( "sst_wl"   , x_dsst               )    ! warm layer
+         CALL iom_put( "sst_cs"   , x_csdsst             )    ! cool skin
+      ENDIF
+       
+      IF( ln_diurnal_only ) RETURN 
+      
       IF( .NOT.lk_vvl ) THEN
          CALL iom_put( "e3t" , fse3t_n(:,:,:) )
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 5676)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 5676)
@@ -0,0 +1,149 @@
+MODULE cool_skin
+   !!======================================================================
+   !!                    ***  MODULE  cool_skin  ***
+   !!     Cool skin thickness and delta T correction using Artele et al. (2002)
+   !!     [see also Tu and Tsuang (2005)]
+   !!
+   !!=====================================================================
+   !! History :        !  2012-01  (P. Sykes)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_coolskin_step  : time-step the cool skin corrections
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE in_out_manager
+   USE sbc_oce
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   
+   IMPLICIT NONE
+   
+   ! Namelist parameters
+
+   ! Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_k = 0.596_wp          ! Thermal conductivity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_v = 1.05e-6_wp        ! Kinematic viscosity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_C = 86400             ! seconds [see Tu and Tsuang (2005)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_cw = 3993._wp         ! specific heat capacity of seawater
+   REAL(wp), PRIVATE, PARAMETER :: pp_h = 10._wp            ! reference depth [using 10m from Artale et al. (2002)]
+   REAL(wp), PRIVATE, PARAMETER :: pp_rhoa = 1.20421_wp     ! density of air (at 20C)
+   REAL(wp), PRIVATE, PARAMETER :: pp_cda = 1.45e-3_wp      ! assumed air-sea drag coefficient for calculating wind speed
+   
+   ! Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csdsst    ! Cool skin delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_csthick   ! Cool skin thickness
+
+   PRIVATE
+   PUBLIC diurnal_sst_coolskin_step, diurnal_sst_coolskin_init
+
+      !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   
+   CONTAINS 
+   
+   SUBROUTINE diurnal_sst_coolskin_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_coolskin_init ***
+      !!
+      !! ** Purpose :   initialise the cool skin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference :
+      !! 
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
+      x_csdsst = 0.
+      x_csthick = 0.
+      
+   END SUBROUTINE diurnal_sst_coolskin_init
+ 
+   SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Time-step the Artale cool skin model
+      !!
+      !! ** Method : 
+      !!
+      !! ** Reference : 
+      !!----------------------------------------------------------------------
+     
+      IMPLICIT NONE
+     
+      ! Dummy variables
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux     ! Heat (non-solar)(Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux   ! Wind stress (kg/ m s^2)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho       ! Water density (kg/m^3)
+      REAL(wp), INTENT(IN) :: rdt                             ! Time-step
+     
+      ! Local variables
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fv                    ! Friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: z_gamma                 ! Dimensionless function of wind speed
+      REAL(wp), DIMENSION(jpi,jpj) :: z_lamda                 ! Sauders (dimensionless) proportionality constant
+      REAL(wp), DIMENSION(jpi,jpj) :: z_wspd                  ! Wind speed (m/s)
+      REAL(wp) :: z_ztx                                       ! Temporary u wind stress
+      REAL(wp) :: z_zty                                       ! Temporary v wind stress
+      REAL(wp) :: z_zmod                                      ! Temporary total wind stress
+     
+      INTEGER :: ji,jj
+      
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            
+            ! Calcualte wind speed from wind stress and friction velocity
+            IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
+               z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
+               z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
+            ELSE
+               z_fv(ji,jj) = 0.
+               z_wspd(ji,jj) = 0.     
+            ENDIF
+
+ 
+            ! Calculate gamma function which is dependent upon wind speed
+            IF( tmask(ji,jj,1) == 1. ) THEN
+               IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
+               IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
+               IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
+            ENDIF
+
+
+            ! Calculate lamda function
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
+               z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
+            ELSE
+               z_lamda(ji,jj) = 0.
+            ENDIF
+
+
+
+            ! Calculate the cool skin thickness - only when heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
+                x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
+            ELSE
+                x_csthick(ji,jj) = 0.
+            ENDIF
+
+
+
+            ! Calculate the cool skin correction - only when the heat flux is out of the ocean
+            IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
+               x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
+             ELSE
+               x_csdsst(ji,jj) = 0.
+            ENDIF
+
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE diurnal_sst_coolskin_step
+
+
+END MODULE cool_skin
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 5676)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 5676)
@@ -0,0 +1,305 @@
+MODULE diurnal_bulk
+   !!======================================================================
+   !!                    ***  MODULE  diurnal_bulk  ***
+   !!     Takaya model of diurnal warming (Takaya, 2010)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_bulk_init  : initialise diurnal model
+   !!   diurnal_sst_bulk_step  : time-step the diurnal model
+   !!----------------------------------------------------------------------
+   USE par_kind
+   USE phycst
+   USE dom_oce
+   USE lib_mpp
+   USE solfrac_mod
+   USE in_out_manager
+   
+   IMPLICIT NONE
+   
+   ! Namelist parameters
+   LOGICAL, PUBLIC :: ln_diurnal
+   LOGICAL, PUBLIC :: ln_diurnal_only
+
+   ! Parameters
+   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
+   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
+   
+   ! Key variables
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     ! Delta SST
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  ! Fraction of 
+   !                                                           ! absorbed radiation
+
+   PRIVATE
+   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
+   
+   CONTAINS 
+   
+   SUBROUTINE diurnal_sst_bulk_init
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_init ***
+      !!
+      !! ** Purpose : Initialise the Takaya diurnal model
+      
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      INTEGER :: ios
+      
+      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
+      
+      ! Read the namelist
+      REWIND( numnam_ref )
+      READ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
+      REWIND( numnam_cfg )
+      READ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )      
+      
+      IF ( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
+         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
+      ENDIF
+      
+      IF ( ln_diurnal ) THEN      
+         
+         ! Allocate arrays
+         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
+         
+         ! Initialise the solar fraction
+         x_solfrac = 0._wp
+         x_dsst = 0._wp
+
+         IF ( ln_diurnal_only ) THEN
+            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
+         ENDIF
+      ENDIF
+      
+   END SUBROUTINE diurnal_sst_bulk_init
+   
+   SUBROUTINE diurnal_sst_takaya_step(psolflux, pqflux, ptauflux, prho, p_rdt,&
+            &                  pla, pthick, pcoolthick, pmu, ld_calcfrac, &
+            &                  p_fvel_bkginc, p_hflux_bkginc)
+      !!----------------------------------------------------------------------
+      !! *** ROUTINE diurnal_sst_takaya_step ***
+      !!
+      !! ** Purpose :   Time-step the Takaya diurnal model
+      !!
+      !! ** Method :    1) Calculate the Obukhov length
+      !!                2) Calculate the Similarity function
+      !!                2) Calculate the increment to dsst
+      !!                3) Apply the increment
+      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
+      !!                temperature, Takaya et al, JGR, 2010 
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+      ! Dummy variables
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psolflux  ! solar flux (Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pqflux    ! heat (non-solar) 
+      !                                                     ! flux (Watts)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: ptauflux  ! wind stress 
+      !                                                     ! (kg/ m s^2)
+      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: prho      ! water density 
+      !                                                     ! (kg/m^3)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pLa    
+      !                                                     ! Langmuir number
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pthick 
+      !                                                     ! warm layer thickness (m)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pcoolthick 
+      !                                                     ! cool skin thickness (m)
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pmu    
+      !                                                     ! mu parameter
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_hflux_bkginc  
+      !                                                     ! optional increment to the
+      !                                                     ! heat flux
+      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_fvel_bkginc  
+      !                                                     ! optional increment to the
+      !                                                     ! friction velocity
+      REAL(wp), INTENT(IN) :: p_rdt                         ! time-step
+      
+      LOGICAL, OPTIONAL, &
+      &     INTENT(IN)  :: ld_calcfrac   ! Set TRUE to recalculate the 
+                                         ! solar fraction
+      
+      ! Local variables 
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
+      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               ! Langmuir function value 
+      
+      LOGICAL  :: ll_calcfrac
+      
+      INTEGER :: ji,jj
+
+      ! Set optional arguments to their defaults
+      IF ( .NOT. PRESENT(pthick) ) THEN 
+         zthick(:,:) = 3._wp
+      ELSE
+         zthick(:,:) = pthick(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pcoolthick) ) THEN 
+         zcoolthick(:,:) = 0._wp
+      ELSE
+         zcoolthick(:,:) = pcoolthick(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pmu) ) THEN
+         zmu(:,:) = 0.3_wp
+      ELSE
+         zmu(:,:) = pmu(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(pla) ) THEN
+         zla(:,:) = 0.3_wp
+      ELSE
+         zla(:,:) = pla(:,:)
+      ENDIF
+      IF ( .NOT. PRESENT(ld_calcfrac) ) THEN
+         ll_calcfrac = .FALSE.
+      ELSE
+         ll_calcfrac = ld_calcfrac
+      ENDIF     
+      
+      ! If not done already, calculate the solar fraction
+      IF (ll_calcfrac ) THEN
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
+               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
+            END DO
+         END DO   
+      ENDIF   
+
+      ! convert solar flux and heat flux to absorbed flux   
+      WHERE ( tmask(:,:,1) == 1._wp) 
+         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
+      ELSEWHERE
+         z_abflux(:,:) = 0._wp
+      ENDWHERE
+      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc   ! Optional increment
+      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
+         z_abflux(:,:) = rsmall
+      ENDWHERE 
+      
+      ! Calculate the friction velocity
+      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
+         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
+      ELSEWHERE
+         z_fvel(:,:) = 0._wp
+      ENDWHERE
+      IF( PRESENT(p_fvel_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc   ! Optional increment
+      
+       
+       
+      ! Calculate the Langmuir function value
+      WHERE ( tmask(:,:,1) == 1.)
+         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )  
+      ELSEWHERE
+         z_fla(:,:) = 0._wp
+      ENDWHERE     
+      
+      ! Increment the temperature using the implicit solution
+      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:), &
+                                 z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:))
+      
+                                
+      
+   END SUBROUTINE diurnal_sst_takaya_step
+
+   
+   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
+                          p_fla, pmu, pthick, prho )
+                          
+      IMPLICIT NONE
+      
+      ! Function definition
+      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
+      ! Dummy variables
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst     ! Delta SST
+      REAL(wp), INTENT(IN)                     :: p_rdt      ! Time-step
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux   ! Heat forcing
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel     ! Friction velocity
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla      ! Langmuir number
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu        ! Structure parameter
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick     ! Layer thickness
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho       ! Water density
+   
+      ! Local variables
+      REAL(wp) :: z_olength          ! Obukhov length
+      REAL(wp) :: z_sigma, z_sigma2
+      REAL(wp) :: z_term1, z_term2     
+      REAL(wp) :: z_stabfunc         ! stability function value
+      REAL(wp) :: z_fvel      
+      
+      CHARACTER(200) :: warn_string
+       
+      INTEGER :: ji,jj
+                     
+      DO jj = 1, jpj
+         DO ji = 1, jpi   
+            
+            ! Only calculate outside tmask
+            IF ( tmask(ji,jj,1) /= 1._wp ) THEN
+               t_imp(ji,jj) = 0._wp
+               CYCLE   
+            END IF
+            
+            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
+               z_fvel = pp_min_fvel
+               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
+               &//"friction velocity < minimum\n" &
+               &//"Setting friction velocity =",pp_min_fvel
+               CALL ctl_warn(warn_string)
+               
+            ELSE
+               z_fvel = p_fvel(ji,jj)
+            ENDIF
+                 
+            ! Calculate the Obukhov length
+            IF ( (z_fvel < pp_veltol ) .AND. &
+            &    (p_dsst(ji,jj) > 0._wp) ) THEN
+               z_olength =  z_fvel  / &
+               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
+               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
+            ELSE
+               z_olength = &
+               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
+               &   ( vkarmn * grav * pp_alpha *&
+               &     p_abflux(ji,jj) )
+            ENDIF
+             
+            ! Calculate the stability function          
+            z_sigma = pthick(ji,jj) / z_olength
+            z_sigma2 = z_sigma * z_sigma
+  
+            IF ( z_sigma >= 0. ) THEN
+               z_stabfunc = 1._wp + &
+               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
+               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
+               &    z_sigma2 ) )
+            ELSE
+               z_stabfunc = 1._wp / &
+               &                   SQRT( 1._wp - 16._wp * z_sigma )
+            ENDIF
+
+            ! Calculate the T increment
+            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
+            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
+            
+             
+            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
+            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
+            &      ( pthick(ji,jj) * z_stabfunc ) )     
+          
+            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
+                           ( 1._wp - p_rdt * z_term2 )
+
+         END DO
+      END DO
+      
+      END FUNCTION t_imp
+
+END MODULE diurnal_bulk
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 5676)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 5676)
@@ -0,0 +1,56 @@
+MODULE solfrac_mod
+   !!======================================================================
+   !!                    ***  MODULE  solfrac  ***
+   !!     POSH representation of solar absorption (Gntermann, 2009)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   solfrac  : function to calculate the solar fraction
+   !!----------------------------------------------------------------------
+   
+   USE par_kind
+   IMPLICIT NONE
+     
+   ! Parameters 
+   REAL(wp), PRIVATE, PARAMETER, DIMENSION(9) :: &
+   &                                     pp_wgt = (/0.2370, 0.36,  0.1790, &
+   &                                                0.087,  0.08,  0.025,  &
+   &                                                0.025,  0.007, 0.0004/)
+   REAL(wp), PRIVATE, PARAMETER, DIMENSION(9) :: &
+   &                                    pp_len = (/34.84,   2.266,   0.0315,  &
+   &                                               0.0055,  8.32e-4, 1.26e-4, &
+   &                                                3.13e-4, 7.82e-4, 1.44e-5/)
+   
+   PUBLIC solfrac
+                                                   
+CONTAINS
+
+   REAL(dp) FUNCTION solfrac(ptop,pbottom)
+       !!----------------------------------------------------------------------
+      !! *** ROUTINE solfrac ***
+      !!
+      !! ** Purpose :   Calculate the solar fraction absorbed between two 
+      !!                layers
+      !!
+      !! ** Reference : POSH a model of diurnal warming, Gentemann et al, 
+      !!                 JGR, 2009 
+      !!----------------------------------------------------------------------
+      
+      ! Dummy variabes
+      REAL(wp), INTENT(IN) :: ptop, pbottom   ! Top and bottom of layer
+      
+      ! local variables
+      INTEGER :: jt
+      
+      ! Calculate the solar fraction absorbed between the two layers
+      solfrac = 0._wp
+      DO jt = 1, 9
+           solfrac = solfrac + pp_wgt(jt) * ( exp ( -ptop / pp_len(jt) ) &
+            &                                 - exp ( -pbottom / pp_len(jt) ) )
+      END DO
+      
+   END FUNCTION
+   
+END MODULE solfrac_mod
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 5676)
@@ -24,5 +24,6 @@
    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
    USE divcur          ! hor. divergence and curl      (div & cur routines)
-
+   USE diurnal_bulk
+   
    IMPLICIT NONE
    PRIVATE
@@ -125,5 +126,5 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       )   ! dynamics time step
                      CALL iom_rstput( kt, nitrst, numrow, 'rdttra1', rdttra(1) )   ! surface tracer time step
-
+      IF ( .NOT. ln_diurnal_only ) THEN
                      CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub        )     ! before fields
                      CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb        )
@@ -145,4 +146,8 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
 #endif
+      ENDIF
+      
+      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    )  
+      
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -221,5 +226,16 @@
          IF( zrdttra1 /= rdttra(1) )   neuler = 0
       ENDIF
-      ! 
+       
+      
+      ! Diurnal DSST 
+      IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst  ) 
+      IF ( ln_diurnal_only ) THEN 
+         IF(lwp) WRITE( numout, * ) &
+         &   "rst_read:- ln_diurnal_only set, setting rhop=rau0" 
+         rhop = rau0
+         CALL iom_get( numror, jpdom_autoglo, 'tn'     , tsn(:,:,1,jp_tem) ) 
+         RETURN 
+      ENDIF  
+      
       IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub      )   ! before fields
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 5676)
@@ -55,4 +55,7 @@
    USE bdy_par          ! Require lk_bdy
 
+   USE diurnal_bulk, ONLY: &
+      & ln_diurnal_only 
+
    IMPLICIT NONE
    PRIVATE
@@ -404,6 +407,9 @@
       IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
 
-      IF( nn_closea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain 
-      !                                                           ! (update freshwater fluxes)
+      ! treatment of closed sea in the model domain 
+      ! (update freshwater fluxes)
+      ! Should not be ran if ln_diurnal_only
+      IF( .NOT.(ln_diurnal_only) .AND. (nn_closea == 1) )   CALL sbc_clo( kt )   
+
 !RBbug do not understand why see ticket 667
 !clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5676)
@@ -76,4 +76,5 @@
 #endif
    USE lib_mpp         ! distributed memory computing
+   USE diurnal_bulk    ! diurnal bulk SST 
 #if defined key_iomput
    USE xios
@@ -169,5 +170,9 @@
             CALL stp                         ! AGRIF: time stepping
 #else
-            CALL stp( istp )                 ! standard time stepping
+            IF ( .NOT. ln_diurnal_only ) THEN 
+               CALL stp( istp )                 ! standard time stepping 
+            ELSE 
+               CALL stp_diurnal( istp )        ! time step only the diurnal SST 
+            ENDIF 
 #endif
             istp = istp + 1
@@ -400,5 +405,24 @@
 
       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
-
+      
+      CALL diurnal_sst_bulk_init            ! diurnal sst
+      IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin   
+      
+      ! IF ln_diurnal_only, then we only want a subset of the initialisation routines
+      IF ( ln_diurnal_only ) THEN
+         CALL  istate_init   ! ocean initial state (Dynamics and tracers)
+         CALL     sbc_init   ! Forcings : surface module
+         CALL tra_qsr_init   ! penetrative solar radiation qsr
+         IF( lk_diaobs     ) THEN                  ! Observation & model comparison
+            CALL dia_obs_init            ! Initialize observational data
+            CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
+         ENDIF     
+         !                                     ! Assimilation increments
+         IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
+                 
+         IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
+         RETURN
+      ENDIF
+      
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5676)
@@ -37,5 +37,5 @@
    PRIVATE
 
-   PUBLIC   stp   ! called by opa.F90
+   PUBLIC   stp, stp_diurnal   ! called by opa.F90
 
    !! * Substitutions
@@ -227,4 +227,16 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! cool skin
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
+      IF ( ln_diurnal ) THEN 
+         IF ( .NOT. ln_blk_core ) THEN
+            CALL ctl_stop("step.f90: diurnal flux processing not implemented"//&
+            &             " for direct forcing")
+         ENDIF
+         CALL diurnal_sst_coolskin_step( & 
+                    qns, taum, rhop(:,:,1), rdt)
+      ENDIF
+      
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs             (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -238,5 +250,20 @@
       !
       IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
-
+      
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! diurnal warm layer
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF ( ln_diurnal ) THEN 
+         IF( kstp == nit000 )THEN 
+            CALL diurnal_sst_takaya_step( & 
+            &    qsr, qns, taum, rhop(:,:,1), &
+            &    rdt, ld_calcfrac = .TRUE.) 
+         ELSE 
+            CALL diurnal_sst_takaya_step( & 
+            &    qsr, qns, taum, rhop(:,:,1), rdt ) 
+         ENDIF
+      ENDIF     
+      
+      
 #if defined key_top
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -386,5 +413,72 @@
       !
    END SUBROUTINE stp
-
+   
+   SUBROUTINE stp_diurnal( kstp ) 
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
+      !!---------------------------------------------------------------------- 
+      !!                     ***  ROUTINE stp  *** 
+      !!                       
+      !! ** Purpose : - Time stepping of diurnal SST model only 
+      !!   
+      !! ** Method  : -1- Update forcings and data   
+      !!              -2- Update ocean physics   
+      !!              -3- Compute the t and s trends   
+      !!              -4- Update t and s   
+      !!              -5- Compute the momentum trends 
+      !!              -6- Update the horizontal velocity 
+      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w) 
+      !!              -8- Outputs and diagnostics 
+      !!---------------------------------------------------------------------- 
+      INTEGER ::   jk       ! dummy loop indices
+      INTEGER ::   indic    ! error indicator if < 0 
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel_bkginc, z_hflux_bkginc     
+      !! --------------------------------------------------------------------- 
+ 
+      indic = 1                                 ! reset to no error condition 
+ 
+      IF( kstp /= nit000 )   CALL day( kstp )   ! Calendar (day was already called at nit000 in day_init) 
+ 
+      CALL iom_setkt( kstp - nit000 + 1, cxios_context )   ! tell iom we are at time step kstp
+      IF( ln_crs ) THEN
+         CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell iom we are at time step kstp
+      ENDIF
+       
+      CALL sbc    ( kstp )                      ! Sea Boundary Conditions 
+     
+      ! Cool skin
+      IF ( .NOT.ln_diurnal ) CALL ctl_stop(  "stp_diurnal: ln_diurnal not set"  )
+         
+      IF ( .NOT. ln_blk_core ) THEN
+         CALL ctl_stop("step.f90: diurnal flux processing not implemented"//&
+         &             " for direct forcing") 
+      ENDIF
+
+      CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rdt)
+
+      CALL dia_wri( kstp )                            ! ocean model: outputs 
+       
+      ! Diurnal warm layer model       
+           
+      IF( kstp == nit000 )THEN 
+         CALL diurnal_sst_takaya_step( & 
+         &    qsr, qns, taum, rhop(:,:,1), &
+         &    rdt, ld_calcfrac = .TRUE.) 
+      ELSE 
+         CALL diurnal_sst_takaya_step( & 
+         &    qsr, qns, taum, rhop(:,:,1), rdt ) 
+      ENDIF
+
+      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+     
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
+      ! Control and restarts 
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
+      IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file 
+      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
+     
+      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset 
+       
+   END SUBROUTINE stp_diurnal  
+   
    !!======================================================================
 END MODULE step
Index: /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5675)
+++ /branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5676)
@@ -85,4 +85,8 @@
    USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
 
+   USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine) 
+   USE cool_skin       ! diurnal cool skin correction (diurnal_sst_coolskin routine)   
+   USE sbc_oce         ! surface fluxes  
+   
    USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
 
