[10989] | 1 | MODULE diu_bulk |
---|
[5676] | 2 | !!====================================================================== |
---|
[10989] | 3 | !! *** MODULE diu_bulk *** |
---|
[5676] | 4 | !! Takaya model of diurnal warming (Takaya, 2010) |
---|
| 5 | !!===================================================================== |
---|
| 6 | !! History : ! 11-10 (J. While) Original code |
---|
| 7 | !!---------------------------------------------------------------------- |
---|
| 8 | |
---|
| 9 | !!---------------------------------------------------------------------- |
---|
| 10 | !! diurnal_sst_bulk_init : initialise diurnal model |
---|
| 11 | !! diurnal_sst_bulk_step : time-step the diurnal model |
---|
| 12 | !!---------------------------------------------------------------------- |
---|
| 13 | USE par_kind |
---|
| 14 | USE phycst |
---|
| 15 | USE dom_oce |
---|
| 16 | USE lib_mpp |
---|
| 17 | USE solfrac_mod |
---|
| 18 | USE in_out_manager |
---|
| 19 | |
---|
| 20 | IMPLICIT NONE |
---|
[9168] | 21 | PRIVATE |
---|
[5676] | 22 | |
---|
| 23 | ! Namelist parameters |
---|
[13942] | 24 | LOGICAL, PUBLIC :: ln_diurnal = .false. ! force definition if diurnal_sst_bulk_init is not called |
---|
| 25 | LOGICAL, PUBLIC :: ln_diurnal_only = .false. ! force definition if diurnal_sst_bulk_init is not called |
---|
[5676] | 26 | |
---|
| 27 | ! Parameters |
---|
| 28 | REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp |
---|
| 29 | REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp |
---|
| 30 | REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp |
---|
| 31 | |
---|
| 32 | ! Key variables |
---|
| 33 | REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst ! Delta SST |
---|
| 34 | REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac ! Fraction of |
---|
| 35 | ! ! absorbed radiation |
---|
| 36 | |
---|
| 37 | PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step |
---|
[12340] | 38 | !! * Substitutions |
---|
| 39 | # include "do_loop_substitute.h90" |
---|
[5676] | 40 | |
---|
[9168] | 41 | !!---------------------------------------------------------------------- |
---|
| 42 | CONTAINS |
---|
[5676] | 43 | |
---|
| 44 | SUBROUTINE diurnal_sst_bulk_init |
---|
| 45 | !!---------------------------------------------------------------------- |
---|
| 46 | !! *** ROUTINE diurnal_sst_init *** |
---|
| 47 | !! |
---|
[9168] | 48 | !! ** Purpose : Initialise the Takaya diurnal model |
---|
| 49 | !!---------------------------------------------------------------------- |
---|
| 50 | INTEGER :: ios ! local integer |
---|
| 51 | !! |
---|
[5676] | 52 | NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only |
---|
[9168] | 53 | !!---------------------------------------------------------------------- |
---|
| 54 | |
---|
[5676] | 55 | ! Read the namelist |
---|
[9168] | 56 | READ ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 ) |
---|
[11822] | 57 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist' ) |
---|
[9168] | 58 | READ ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 ) |
---|
[11822] | 59 | 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist' ) |
---|
[9168] | 60 | ! |
---|
| 61 | IF( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN |
---|
[5676] | 62 | CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" ) |
---|
| 63 | ENDIF |
---|
| 64 | |
---|
[9168] | 65 | IF( ln_diurnal ) THEN |
---|
| 66 | ! |
---|
[5676] | 67 | ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) ) |
---|
[9168] | 68 | ! |
---|
| 69 | x_solfrac = 0._wp ! Initialise the solar fraction |
---|
| 70 | x_dsst = 0._wp |
---|
| 71 | ! |
---|
| 72 | IF( ln_diurnal_only ) THEN |
---|
[5676] | 73 | CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" ) |
---|
| 74 | ENDIF |
---|
| 75 | ENDIF |
---|
| 76 | |
---|
| 77 | END SUBROUTINE diurnal_sst_bulk_init |
---|
[9168] | 78 | |
---|
| 79 | |
---|
| 80 | SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt, & |
---|
[6010] | 81 | & pla, pthick, pcoolthick, pmu, & |
---|
[5676] | 82 | & p_fvel_bkginc, p_hflux_bkginc) |
---|
| 83 | !!---------------------------------------------------------------------- |
---|
| 84 | !! *** ROUTINE diurnal_sst_takaya_step *** |
---|
| 85 | !! |
---|
| 86 | !! ** Purpose : Time-step the Takaya diurnal model |
---|
| 87 | !! |
---|
| 88 | !! ** Method : 1) Calculate the Obukhov length |
---|
| 89 | !! 2) Calculate the Similarity function |
---|
| 90 | !! 2) Calculate the increment to dsst |
---|
| 91 | !! 3) Apply the increment |
---|
| 92 | !! ** Reference : Refinements to a prognostic scheme of skin sea surface |
---|
| 93 | !! temperature, Takaya et al, JGR, 2010 |
---|
| 94 | !!---------------------------------------------------------------------- |
---|
[9168] | 95 | INTEGER , INTENT(in) :: kt ! time step |
---|
| 96 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: psolflux ! solar flux (Watts) |
---|
| 97 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts) |
---|
| 98 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2) |
---|
| 99 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3) |
---|
| 100 | REAL(wp) , INTENT(in) :: p_rdt ! time-step |
---|
| 101 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pLa ! Langmuir number |
---|
| 102 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m) |
---|
| 103 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m) |
---|
| 104 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pmu ! mu parameter |
---|
| 105 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux |
---|
| 106 | REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity |
---|
| 107 | ! |
---|
| 108 | INTEGER :: ji,jj |
---|
| 109 | LOGICAL :: ll_calcfrac |
---|
[5676] | 110 | REAL(wp), DIMENSION(jpi,jpj) :: z_fvel ! friction velocity |
---|
| 111 | REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla |
---|
| 112 | REAL(wp), DIMENSION(jpi,jpj) :: z_abflux ! absorbed flux |
---|
| 113 | REAL(wp), DIMENSION(jpi,jpj) :: z_fla ! Langmuir function value |
---|
[9168] | 114 | !!---------------------------------------------------------------------- |
---|
[5676] | 115 | |
---|
| 116 | ! Set optional arguments to their defaults |
---|
[9168] | 117 | IF( .NOT. PRESENT( pthick ) ) THEN ; zthick(:,:) = 3._wp |
---|
| 118 | ELSE ; zthick(:,:) = pthick(:,:) |
---|
[5676] | 119 | ENDIF |
---|
[9168] | 120 | IF( .NOT. PRESENT(pcoolthick) ) THEN ; zcoolthick(:,:) = 0._wp |
---|
| 121 | ELSE ; zcoolthick(:,:) = pcoolthick(:,:) |
---|
[5676] | 122 | ENDIF |
---|
[9168] | 123 | IF( .NOT. PRESENT( pmu ) ) THEN ; zmu(:,:) = 0.3_wp |
---|
| 124 | ELSE ; zmu(:,:) = pmu(:,:) |
---|
[5676] | 125 | ENDIF |
---|
[9168] | 126 | IF( .NOT. PRESENT(pla) ) THEN ; zla(:,:) = 0.3_wp |
---|
| 127 | ELSE ; zla(:,:) = pla(:,:) |
---|
[5676] | 128 | ENDIF |
---|
| 129 | |
---|
| 130 | ! If not done already, calculate the solar fraction |
---|
[6010] | 131 | IF ( kt==nit000 ) THEN |
---|
[13295] | 132 | DO_2D( 1, 1, 1, 1 ) |
---|
[12340] | 133 | IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) & |
---|
| 134 | & x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) ) |
---|
| 135 | END_2D |
---|
[5676] | 136 | ENDIF |
---|
| 137 | |
---|
| 138 | ! convert solar flux and heat flux to absorbed flux |
---|
| 139 | WHERE ( tmask(:,:,1) == 1._wp) |
---|
| 140 | z_abflux(:,:) = ( x_solfrac(:,:) * psolflux (:,:)) + pqflux(:,:) |
---|
| 141 | ELSEWHERE |
---|
| 142 | z_abflux(:,:) = 0._wp |
---|
| 143 | ENDWHERE |
---|
| 144 | IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc ! Optional increment |
---|
| 145 | WHERE ( ABS( z_abflux(:,:) ) < rsmall ) |
---|
| 146 | z_abflux(:,:) = rsmall |
---|
| 147 | ENDWHERE |
---|
| 148 | |
---|
| 149 | ! Calculate the friction velocity |
---|
| 150 | WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) ) |
---|
| 151 | z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) ) |
---|
| 152 | ELSEWHERE |
---|
| 153 | z_fvel(:,:) = 0._wp |
---|
| 154 | ENDWHERE |
---|
| 155 | IF( PRESENT(p_fvel_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc ! Optional increment |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | ! Calculate the Langmuir function value |
---|
| 160 | WHERE ( tmask(:,:,1) == 1.) |
---|
| 161 | z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) ) |
---|
| 162 | ELSEWHERE |
---|
| 163 | z_fla(:,:) = 0._wp |
---|
| 164 | ENDWHERE |
---|
| 165 | |
---|
| 166 | ! Increment the temperature using the implicit solution |
---|
[9168] | 167 | x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:), & |
---|
| 168 | & z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) ) |
---|
| 169 | ! |
---|
[5676] | 170 | END SUBROUTINE diurnal_sst_takaya_step |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, & |
---|
| 174 | p_fla, pmu, pthick, prho ) |
---|
| 175 | |
---|
| 176 | IMPLICIT NONE |
---|
| 177 | |
---|
| 178 | ! Function definition |
---|
| 179 | REAL(wp), DIMENSION(jpi,jpj) :: t_imp |
---|
| 180 | ! Dummy variables |
---|
| 181 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst ! Delta SST |
---|
| 182 | REAL(wp), INTENT(IN) :: p_rdt ! Time-step |
---|
| 183 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux ! Heat forcing |
---|
| 184 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel ! Friction velocity |
---|
| 185 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla ! Langmuir number |
---|
| 186 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu ! Structure parameter |
---|
| 187 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick ! Layer thickness |
---|
| 188 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho ! Water density |
---|
| 189 | |
---|
| 190 | ! Local variables |
---|
| 191 | REAL(wp) :: z_olength ! Obukhov length |
---|
| 192 | REAL(wp) :: z_sigma, z_sigma2 |
---|
| 193 | REAL(wp) :: z_term1, z_term2 |
---|
| 194 | REAL(wp) :: z_stabfunc ! stability function value |
---|
| 195 | REAL(wp) :: z_fvel |
---|
| 196 | |
---|
| 197 | CHARACTER(200) :: warn_string |
---|
| 198 | |
---|
| 199 | INTEGER :: ji,jj |
---|
| 200 | |
---|
[13295] | 201 | DO_2D( 1, 1, 1, 1 ) |
---|
[12340] | 202 | |
---|
| 203 | ! Only calculate outside tmask |
---|
| 204 | IF ( tmask(ji,jj,1) /= 1._wp ) THEN |
---|
| 205 | t_imp(ji,jj) = 0._wp |
---|
| 206 | CYCLE |
---|
| 207 | END IF |
---|
| 208 | |
---|
| 209 | IF (p_fvel(ji,jj) < pp_min_fvel) THEN |
---|
| 210 | z_fvel = pp_min_fvel |
---|
| 211 | WRITE(warn_string,*) "diurnal_sst_takaya step: "& |
---|
| 212 | &//"friction velocity < minimum\n" & |
---|
| 213 | &//"Setting friction velocity =",pp_min_fvel |
---|
| 214 | CALL ctl_warn(warn_string) |
---|
[5676] | 215 | |
---|
[12340] | 216 | ELSE |
---|
| 217 | z_fvel = p_fvel(ji,jj) |
---|
| 218 | ENDIF |
---|
| 219 | |
---|
| 220 | ! Calculate the Obukhov length |
---|
| 221 | IF ( (z_fvel < pp_veltol ) .AND. & |
---|
| 222 | & (p_dsst(ji,jj) > 0._wp) ) THEN |
---|
| 223 | z_olength = z_fvel / & |
---|
| 224 | & SQRT( p_dsst(ji,jj) * vkarmn * grav * & |
---|
| 225 | & pp_alpha / ( 5._wp * pthick(ji,jj) ) ) |
---|
| 226 | ELSE |
---|
| 227 | z_olength = & |
---|
| 228 | & ( prho(ji,jj) * rcp * z_fvel**3._wp ) / & |
---|
| 229 | & ( vkarmn * grav * pp_alpha *& |
---|
| 230 | & p_abflux(ji,jj) ) |
---|
| 231 | ENDIF |
---|
| 232 | |
---|
| 233 | ! Calculate the stability function |
---|
| 234 | z_sigma = pthick(ji,jj) / z_olength |
---|
| 235 | z_sigma2 = z_sigma * z_sigma |
---|
| 236 | IF ( z_sigma >= 0. ) THEN |
---|
| 237 | z_stabfunc = 1._wp + & |
---|
| 238 | & ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / & |
---|
| 239 | & ( 1._wp + 3._wp * z_sigma + 0.25_wp * & |
---|
| 240 | & z_sigma2 ) ) |
---|
| 241 | ELSE |
---|
| 242 | z_stabfunc = 1._wp / & |
---|
| 243 | & SQRT( 1._wp - 16._wp * z_sigma ) |
---|
| 244 | ENDIF |
---|
[5676] | 245 | |
---|
[12340] | 246 | ! Calculate the T increment |
---|
| 247 | z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp) / & |
---|
| 248 | & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) ) |
---|
| 249 | |
---|
[5676] | 250 | |
---|
[12340] | 251 | z_term2 = -( ( pmu(ji,jj) + 1._wp) * & |
---|
| 252 | & ( vkarmn * z_fvel * p_fla(ji,jj) ) / & |
---|
| 253 | & ( pthick(ji,jj) * z_stabfunc ) ) |
---|
| 254 | |
---|
| 255 | t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / & |
---|
| 256 | ( 1._wp - p_rdt * z_term2 ) |
---|
[5676] | 257 | |
---|
[12340] | 258 | END_2D |
---|
[5676] | 259 | |
---|
| 260 | END FUNCTION t_imp |
---|
| 261 | |
---|
[10989] | 262 | END MODULE diu_bulk |
---|