Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 15473)
@@ -22,4 +22,6 @@
       
       <field_group id="grid_T" grid_ref="grid_T_2D" >
+         <field id="sst_wl"       long_name="diurnal warm layer SST" standard_name="warm_layer" unit="degC"     grid_ref="grid_T_2D" />
+         <field id="sst_cs"       long_name="diurnal cool skin SST" standard_name="cool_skin" unit="degC"     grid_ref="grid_T_2D" />
          <field id="e3t"          long_name="Ocean Model cell Thickness"   standard_name="cell_thickness"   unit="m"   grid_ref="grid_T_3D"/>
          <field id="e3t_surf"     field_ref="e3t"  long_name="T-cell thickness: surface layer"   standard_name="cell_thickness"   unit="m"  grid_ref="grid_T_surface_extract"/>
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 15473)
@@ -1393,2 +1393,8 @@
    ln_incpc       = .FALSE.
 /
+!-----------------------------------------------------------------------
+&namdiu
+!-----------------------------------------------------------------------
+   ln_diurnal      = .false.
+   ln_diurnal_only = .false.
+/
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 15473)
@@ -57,4 +57,6 @@
    USE timing          ! preformance summary
    USE wrk_nemo        ! working array
+   USE diurnal_bulk    ! diurnal warm layer
+   USE cool_skin       ! Cool skin
 
    IMPLICIT NONE
@@ -144,4 +146,10 @@
          CALL dia_wri_state( 'output.init', kt )
          ninist = 0
+      ENDIF
+      
+      IF(ln_diurnal) THEN
+         CALL iom_put( "sst_wl"   , x_dsst               )    ! warm layer
+         CALL iom_put( "sst_cs"   , x_csdsst             )    ! cool skin
+         !CALL iom_put( "temp_3m"  , tsn(:,:,4,jp_tem)    )    ! temperature at ~3m
       ENDIF
 
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 15473)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 15473)
@@ -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  : timestep 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 :   initalise the coolskin 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 :   Timestep 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                             !Timestep
+     
+      !Local variables
+      REAL(wp), DIMENSION(jpi,jpj) :: z_fv                    !Friction velocity     
+      REAL(wp), DIMENSION(jpi,jpj) :: z_gamma                 !Dimensionless function of windspeed
+      REAL(wp), DIMENSION(jpi,jpj) :: z_lamda                 !Sauders (dimensionless) proportionality constant
+      REAL(wp), DIMENSION(jpi,jpj) :: z_wspd                  !Windspeed (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
+
+
+            !Calulate 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/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 15473)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90	(revision 15473)
@@ -0,0 +1,307 @@
+MODULE diurnal_bulk
+   !!======================================================================
+   !!                    ***  MODULE  diurnal_bulk  ***
+   !!     Takaya model of diurnal warming (Takaya, 2010)
+   !!=====================================================================
+   !! History :        !  11-10  (J. While)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   diurnal_sst_bulk_init  : initalise diurnal model
+   !!   diurnal_sst_bulk_step  : timestep 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 = .FALSE.
+   LOGICAL, PUBLIC :: ln_diurnal_only = .FALSE.
+
+   !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 : Initalise the Takaya diurnal model
+      
+      !!----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+      INTEGER :: ios                 ! Local integer output status for namelist read
+      
+      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
+      
+      !Read the namelist
+      !REWIND( numnam )
+      !READ( numnam, namdiu )
+
+      REWIND( numnam_ref )              ! Namelist namdiu in reference namelist
+      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namdiu in configuration namelist
+      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) )
+         
+         !Initalise the solar fraction
+         x_solfrac=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 :   Timestep the Takaya diurnal model
+      !!
+      !! ** Method :    1) Calculate the Obukhov length
+      !!                2) Calculate the Simularity 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              !timestep
+      
+      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 fuction value 
+      
+      LOGICAL  :: ll_calcfrac
+      
+      INTEGER :: ji,jj
+
+      !Set optional arguments to their defualts
+      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. ) .AND. ( tmask(ji,jj,1) == 1. ) ) &
+               &   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.) 
+         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.
+      ENDWHERE
+      IF( PRESENT(p_hflux_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 ) )  
+      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           !Timestep
+      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. ) 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.) ) 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/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 15473)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/solfrac_mod.F90	(revision 15473)
@@ -0,0 +1,56 @@
+MODULE solfrac_mod
+   !!======================================================================
+   !!                    ***  MODULE  solfrac  ***
+   !!     POSH representation of solar absorbtion (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/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 15473)
@@ -27,4 +27,5 @@
    USE sbc_oce         ! for icesheet freshwater input variables
    USE timing
+   USE diurnal_bulk    ! diurnal SST
 
    IMPLICIT NONE
@@ -136,4 +137,5 @@
       !!----------------------------------------------------------------------
                      IF(nn_timing == 2)  CALL timing_start('iom_rstput')
+                     IF( ln_diurnal ) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    ) 
                      CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       )   ! dynamics time step
                      CALL iom_rstput( kt, nitrst, numrow, 'rdttra1', rdttra(1) )   ! surface tracer time step
@@ -260,4 +262,11 @@
          neuler = 0
       ENDIF
+      IF( ln_diurnal ) THEN
+         IF( iom_varid( numror, 'Dsst', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst  )
+         ELSE
+            x_dsst(:,:) = 0.0
+         ENDIF
+      ENDIF
       !
       CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 15473)
@@ -89,4 +89,5 @@
    USE stopar
    USE stopts
+   USE diurnal_bulk    ! diurnal bulk SST
 
    IMPLICIT NONE
@@ -422,4 +423,7 @@
 
       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
 
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 15473)
@@ -253,4 +253,13 @@
       ENDIF
 
+      ! Cool skin
+      IF ( ln_diurnal ) THEN 
+         IF ( ln_blk_core ) THEN
+            CALL diurnal_sst_coolskin_step(qns(:,:)+(rn_abs*qsr(:,:)), taum, rhop(:,:,1), rdt)
+         ELSE
+            CALL diurnal_sst_coolskin_step(qns, taum, rhop(:,:,1), rdt)
+         ENDIF
+      ENDIF
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs             (ua, va, tsa used as workspace)
@@ -266,4 +275,29 @@
       !
       IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
+
+      !Diurnal warm layer model       
+      IF ( ln_diurnal ) THEN
+         IF ( ln_blk_core ) THEN
+            IF( kstp == nit000 )THEN 
+               CALL diurnal_sst_takaya_step( & 
+                  &    qsr(:,:)-(rn_abs*qsr(:,:)), qns(:,:)+(rn_abs*qsr(:,:)), &
+                  &    taum, rhop(:,:,1), &
+                  &    rdt, ld_calcfrac = .TRUE.) 
+            ELSE 
+               CALL diurnal_sst_takaya_step( & 
+                  &    qsr(:,:)-(rn_abs*qsr(:,:)), qns(:,:)+(rn_abs*qsr(:,:)), &
+                  &    taum, rhop(:,:,1), rdt ) 
+            ENDIF
+         ELSE
+            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
+      ENDIF
 
 #if defined key_top
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 15473)
@@ -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)
 
Index: /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 15472)
+++ /branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 15473)
@@ -96,4 +96,5 @@
       USE lib_mpp,           ONLY: ctl_stop
       USE trcstat,           ONLY: trc_rst_dia_stat 
+      USE cool_skin,         ONLY: x_csdsst
 
    !!* Substitution
@@ -228,5 +229,5 @@
                !!                 chemistry package; note that depth is set to
                !!                 zero in this call
-               CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj),     &
+               CALL mocsy_interface(ztmp(ji,jj)+x_csdsst(ji,jj),zsal(ji,jj),zalk(ji,jj),     &
                                     zdic(ji,jj),zsil(ji,jj),zpho(ji,jj),     &
                                     f_pp0(ji,jj),0.0,                        &
