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