Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
MODULE diu_bulk
!!======================================================================
!! *** MODULE diu_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
PRIVATE
! Namelist parameters
LOGICAL, PUBLIC :: ln_diurnal = .false. ! force definition if diurnal_sst_bulk_init is not called
LOGICAL, PUBLIC :: ln_diurnal_only = .false. ! force definition if diurnal_sst_bulk_init is not called
! 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
PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE diurnal_sst_bulk_init
!!----------------------------------------------------------------------
!! *** ROUTINE diurnal_sst_init ***
!!
!! ** Purpose : Initialise the Takaya diurnal model
!!----------------------------------------------------------------------
INTEGER :: ios ! local integer
!!
NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
!!----------------------------------------------------------------------
! Read the namelist
READ ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist' )
READ ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist' )
!
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( x_dsst(A2D(0)), x_solfrac(A2D(0)) )
!
x_solfrac = 0._wp ! Initialise the solar fraction
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(kt, psolflux, pqflux, ptauflux, prho, p_rdt, &
& pla, pthick, pcoolthick, pmu, &
& 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
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
REAL(wp), DIMENSION(A2D(0)) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(A2D(0)) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(A2D(0)) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(A2D(0)) :: z_fla ! Langmuir function value
!!----------------------------------------------------------------------
! 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 done already, calculate the solar fraction
IF ( kt==nit000 ) THEN
DO_2D( 0, 0, 0, 0 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( smask0(ji,jj) == 1._wp ) ) &
& x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
END_2D
ENDIF
! convert solar flux and heat flux to absorbed flux
WHERE ( smask0(:,:) == 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. ( smask0(:,:) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(A2D(0)) )
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 ( smask0(:,:) == 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(A2D(0)) :: t_imp
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(A2D(0)), 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
IF ( smask0(ji,jj) /= 1._wp ) THEN
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
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_2D
END FUNCTION t_imp
END MODULE diu_bulk