1 | SUBROUTINE forcing |
---|
2 | |
---|
3 | !!------------------------------------------------------------------ |
---|
4 | !! *** ROUTINE forcing *** |
---|
5 | !! ** Purpose : |
---|
6 | !! This routine computes the model forcing |
---|
7 | !! forc_swi = 0 -> read |
---|
8 | !! 1 -> computed |
---|
9 | !! 99 -> prescribed |
---|
10 | !! |
---|
11 | !! ** Method : |
---|
12 | !! |
---|
13 | !! ** Arguments : |
---|
14 | !! |
---|
15 | !! ** Inputs / Ouputs : (global commons) |
---|
16 | !! |
---|
17 | !! ** External : |
---|
18 | !! |
---|
19 | !! ** References : |
---|
20 | !! |
---|
21 | !! ** History : |
---|
22 | !! |
---|
23 | !!------------------------------------------------------------------ |
---|
24 | !! * Arguments |
---|
25 | |
---|
26 | USE lib_fortran |
---|
27 | |
---|
28 | INCLUDE 'type.com' |
---|
29 | INCLUDE 'para.com' |
---|
30 | INCLUDE 'const.com' |
---|
31 | INCLUDE 'ice.com' |
---|
32 | INCLUDE 'thermo.com' |
---|
33 | INCLUDE 'forcing.com' |
---|
34 | |
---|
35 | INTEGER :: |
---|
36 | & ji , ! : index for space |
---|
37 | & jk , ! : index for ice layers |
---|
38 | & jf , ! : index for forcing field |
---|
39 | & numforc= 600 ! : reference number for bio.param |
---|
40 | |
---|
41 | CHARACTER(len=10) :: |
---|
42 | & filenc='forcing.nc' |
---|
43 | |
---|
44 | REAL(4) zforc(1), zforc2(2) ! forcing field dummy array |
---|
45 | DIMENSION ws(96),zmue(96),zalcnp(96) ! for solar flux formula |
---|
46 | |
---|
47 | REAL(8) :: |
---|
48 | & zhour |
---|
49 | |
---|
50 | DIMENSION budyko(19) |
---|
51 | DATA budyko /1.00,0.98,0.95,0.92,0.89,0.86,0.83,0.80,0.78,0.75, |
---|
52 | & 0.72,0.69,0.67,0.64,0.61,0.58,0.56,0.53,0.50/ |
---|
53 | |
---|
54 | LOGICAL :: |
---|
55 | & ln_write_forc |
---|
56 | |
---|
57 | ln_write_forc = .TRUE. |
---|
58 | |
---|
59 | WRITE(numout,*) ' * forcing_nc : ' |
---|
60 | WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' |
---|
61 | |
---|
62 | ! Control parameters number of steps in the integration of shortwave rad |
---|
63 | nintsr = 24 ! 24 ideally |
---|
64 | zsolar = 1368. ! solar constant |
---|
65 | ji = 1 |
---|
66 | |
---|
67 | ! |
---|
68 | !-----------------------------------------------------------------------------! |
---|
69 | ! 1) Names of the forcing variables |
---|
70 | !-----------------------------------------------------------------------------! |
---|
71 | ! |
---|
72 | forc_nam(1) = 'fsw' |
---|
73 | forc_nam(2) = 'flw' |
---|
74 | forc_nam(3) = 'par' |
---|
75 | forc_nam(4) = 'tair' |
---|
76 | forc_nam(5) = 'pres' |
---|
77 | forc_nam(6) = 'qair' |
---|
78 | forc_nam(7) = 'wspd' |
---|
79 | forc_nam(8) = 'cld' |
---|
80 | forc_nam(9) = 'foce' |
---|
81 | forc_nam(10)= 'sfal' |
---|
82 | forc_nam(11)= 'albe' |
---|
83 | ! |
---|
84 | !-----------------------------------------------------------------------------! |
---|
85 | ! 2) Reads namelist |
---|
86 | !-----------------------------------------------------------------------------! |
---|
87 | ! |
---|
88 | IF ( numit .EQ. nstart ) THEN |
---|
89 | |
---|
90 | IF ( ln_write_forc ) THEN |
---|
91 | WRITE(numout,*) ' Forcing parameters ... ' |
---|
92 | WRITE(numout,*) |
---|
93 | WRITE(numout,*) ' forc_nam : ',( forc_nam(i), i = 1, n_forc) |
---|
94 | ENDIF |
---|
95 | |
---|
96 | OPEN( unit = numforc, file='forcing.param', status='old' ) |
---|
97 | READ(numforc,*) |
---|
98 | READ(numforc,*) |
---|
99 | READ(numforc,*) |
---|
100 | READ(numforc,*) |
---|
101 | READ(numforc,*) |
---|
102 | READ(numforc,*) ts_forc ! Forcing time step |
---|
103 | READ(numforc,*) |
---|
104 | READ(numforc,*) n0_forc ! Number of the first time step in the forcing file |
---|
105 | READ(numforc,*) |
---|
106 | READ(numforc,*) n1_forc ! Number of the last time step in the forcing file |
---|
107 | IF ( ln_write_forc ) WRITE(numout,*) ' ts_forc : ', ts_forc |
---|
108 | |
---|
109 | READ(numforc,*) |
---|
110 | READ(numforc,*) |
---|
111 | READ(numforc,*) |
---|
112 | READ(numforc,*) |
---|
113 | READ(numforc,*) |
---|
114 | READ(numforc,*) |
---|
115 | READ(numforc,*) |
---|
116 | |
---|
117 | IF ( ln_write_forc ) WRITE(numout,*) ' forc_swi, forc_uni, |
---|
118 | & forc_val : ' |
---|
119 | DO i = 1, n_forc |
---|
120 | READ(numforc,*) idum, dummy1, dummy2 |
---|
121 | READ(numforc,*) |
---|
122 | forc_swi(i) = idum |
---|
123 | forc_uni(i) = dummy1 |
---|
124 | forc_val(i) = dummy2 |
---|
125 | IF ( ln_write_forc ) WRITE(numout,*) forc_nam(i), forc_swi(i), |
---|
126 | & forc_uni(i), forc_val(i) |
---|
127 | END DO |
---|
128 | |
---|
129 | READ(numforc,*) |
---|
130 | READ(numforc,*) |
---|
131 | READ(numforc,*) |
---|
132 | READ(numforc,*) opt_dept |
---|
133 | WRITE(numout,*) 'opt_dept : ', opt_dept |
---|
134 | |
---|
135 | CALL CF_OPEN (filenc,id) ! open forcing file |
---|
136 | |
---|
137 | IF ( ln_write_forc ) WRITE(numout,*) |
---|
138 | |
---|
139 | ! number of days in a year |
---|
140 | yeaday = 365.0 |
---|
141 | |
---|
142 | ENDIF ! numit=nstart |
---|
143 | ! |
---|
144 | !-----------------------------------------------------------------------------! |
---|
145 | ! 3) Calendar parameters |
---|
146 | !-----------------------------------------------------------------------------! |
---|
147 | ! |
---|
148 | IF ( ln_write_forc ) THEN |
---|
149 | WRITE(numout,*) |
---|
150 | WRITE(numout,*) ' Calendar ' |
---|
151 | WRITE(numout,*) |
---|
152 | ENDIF |
---|
153 | |
---|
154 | ! Number of days in the year |
---|
155 | ! and number of the year |
---|
156 | IF ( num_d .EQ. INT( yeaday+1 ) ) THEN |
---|
157 | nyear = nyear + 1 |
---|
158 | yeaday = 365.0 |
---|
159 | IF ( MOD(nyear , 4) .EQ. 0 ) THEN |
---|
160 | yeaday = 366.0 |
---|
161 | ENDIF |
---|
162 | ENDIF |
---|
163 | |
---|
164 | ! day of year |
---|
165 | doy = REAL(numit) * ddtb / 86400. ! this does not reset xjour for the new year |
---|
166 | zz_years = INT ( doy / 365.) + 1 ! number of years in the forcing |
---|
167 | |
---|
168 | doy = doy - ( zz_years - 1 ) * 365. ! this is not yet perfect for leap years but error is very small |
---|
169 | num_d = INT(doy) |
---|
170 | |
---|
171 | WRITE(numout,*) ' zz_years: ', zz_years |
---|
172 | WRITE(numout,*) ' ddtb : ', ddtb |
---|
173 | WRITE(numout,*) ' nyear : ', nyear |
---|
174 | WRITE(numout,*) ' doy : ', doy |
---|
175 | WRITE(numout,*) ' num_d : ', num_d |
---|
176 | ! |
---|
177 | !-----------------------------------------------------------------------------! |
---|
178 | ! 3) Read the variables |
---|
179 | !-----------------------------------------------------------------------------! |
---|
180 | ! |
---|
181 | |
---|
182 | IF ( numit .EQ. nstart ) THEN |
---|
183 | n_fofr = INT( ts_forc / ddtb ) ! forcing frequency |
---|
184 | n_forc_min = n_fofr / 2 + nstart |
---|
185 | n_forc_max = ( nitrun - INT( FLOAT(n_fofr) / 2. ) ) + nstart |
---|
186 | WRITE(numout,*) ' n_fofr : ', n_fofr |
---|
187 | WRITE(numout,*) ' nstart : ', nstart |
---|
188 | WRITE(numout,*) ' n_forc_min : ', n_forc_min |
---|
189 | WRITE(numout,*) ' n_forc_max : ', n_forc_max |
---|
190 | WRITE(numout,*) ' nend : ', nend |
---|
191 | WRITE(numout,*) ' n0_forc : ', n0_forc |
---|
192 | i_forc_ts = nstart / n_fofr - n0_forc + 1 |
---|
193 | WRITE(numout,*) ' i_forc_ts : ', i_forc_ts |
---|
194 | ENDIF |
---|
195 | |
---|
196 | i_forc = MOD( ( numit - nstart ) - n_fofr / 2 , n_fofr) ! indicates if forcing is to be read |
---|
197 | WRITE(numout,*) ' numit - nstart : ', numit - nstart |
---|
198 | WRITE(numout,*) ' n_fofr : ', n_fofr |
---|
199 | WRITE(numout,*) ' i_forc : ', i_forc |
---|
200 | |
---|
201 | !------------------------------------- |
---|
202 | ! Time steps at which values are READ |
---|
203 | !------------------------------------- |
---|
204 | IF ( ( numit .EQ. nstart ) .OR. ( i_forc .EQ. 0.0 ) ) THEN |
---|
205 | i_forc_count = 0 |
---|
206 | i_forc_ts = i_forc_ts + 1 |
---|
207 | |
---|
208 | IF ( i_forc_ts .GT. n1_forc ) i_forc_ts = 1 |
---|
209 | |
---|
210 | WRITE(numout,*) ' Forcing read at this time step, numit : ', |
---|
211 | & numit |
---|
212 | WRITE(numout,*) ' Forcing step, i_forc_ts : ', i_forc_ts |
---|
213 | WRITE(numout,*) ' n_forc : ', n_forc |
---|
214 | |
---|
215 | DO i = 1, n_forc |
---|
216 | |
---|
217 | IF ( forc_swi(i) .EQ. 0 ) THEN |
---|
218 | IF ( numit .GT. nstart ) forc_arr_old(i) =forc_arr_new(i) |
---|
219 | IF ( numit .LT. n_forc_max ) THEN |
---|
220 | CALL CF_READ1D ( filenc, forc_nam(i), |
---|
221 | & i_forc_ts , 1, zforc ) |
---|
222 | |
---|
223 | forc_arr_new(i) = REAL(zforc(1)) * forc_uni(i) |
---|
224 | ENDIF |
---|
225 | IF ( numit .EQ. nstart ) forc_arr_old(i) =forc_arr_new(i) |
---|
226 | forc_coeff(i) = ( forc_arr_new(i) - forc_arr_old(i) ) / |
---|
227 | & FLOAT(n_fofr) |
---|
228 | ENDIF |
---|
229 | |
---|
230 | END DO |
---|
231 | |
---|
232 | ENDIF |
---|
233 | |
---|
234 | !----------------------------------------- |
---|
235 | ! Time steps at which values are COMPUTED |
---|
236 | !----------------------------------------- |
---|
237 | i_forc_count = i_forc_count + 1 |
---|
238 | DO i = 1, n_forc |
---|
239 | IF ( forc_swi(i) .EQ. 0 ) forc_arr(i) = forc_arr_old(i) |
---|
240 | & + forc_coeff(i) * FLOAT(i_forc_count) |
---|
241 | END DO |
---|
242 | |
---|
243 | ! Recover arrays (dirty patch to remove in the end) |
---|
244 | IF ( forc_swi(1) .EQ. 0 ) zfsw = forc_arr(1) |
---|
245 | IF ( forc_swi(2) .EQ. 0 ) zflw = forc_arr(2) |
---|
246 | IF ( forc_swi(3) .EQ. 0 ) zpar = forc_arr(3) |
---|
247 | IF ( forc_swi(4) .EQ. 0 ) ztair= forc_arr(4) |
---|
248 | IF ( forc_swi(5) .EQ. 0 ) zpres= forc_arr(5) |
---|
249 | IF ( forc_swi(6) .EQ. 0 ) zqair= forc_arr(6) |
---|
250 | IF ( forc_swi(7) .EQ. 0 ) zwspd= forc_arr(7) |
---|
251 | IF ( forc_swi(8) .EQ. 0 ) zcld = forc_arr(8) |
---|
252 | IF ( forc_swi(9) .EQ. 0 ) zfoce= forc_arr(9) |
---|
253 | IF ( forc_swi(10).EQ. 0 ) zsfal= forc_arr(10) |
---|
254 | |
---|
255 | IF ( forc_swi(1) .EQ. 99 ) zfsw = forc_val(1) * forc_uni(1) |
---|
256 | IF ( forc_swi(2) .EQ. 99 ) zflw = forc_val(2) * forc_uni(2) |
---|
257 | IF ( forc_swi(3) .EQ. 99 ) zpar = forc_val(3) * forc_uni(3) |
---|
258 | IF ( forc_swi(4) .EQ. 99 ) ztair = forc_val(4) * forc_uni(4) |
---|
259 | IF ( forc_swi(5) .EQ. 99 ) zpres = forc_val(5) * forc_uni(5) |
---|
260 | IF ( forc_swi(6) .EQ. 99 ) zqair = forc_val(6) * forc_uni(6) |
---|
261 | IF ( forc_swi(7) .EQ. 99 ) zwspd = forc_val(7) * forc_uni(7) |
---|
262 | IF ( forc_swi(8) .EQ. 99 ) zcld = forc_val(8) * forc_uni(8) |
---|
263 | IF ( forc_swi(9) .EQ. 99 ) zfoce = forc_val(9) * forc_uni(9) |
---|
264 | IF ( forc_swi(10).EQ. 99 ) zsfal = forc_val(10) * forc_uni(10) |
---|
265 | IF ( forc_swi(11).EQ. 99 ) zalbe = forc_val(11) * forc_uni(11) |
---|
266 | |
---|
267 | IF ( forc_swi(11).EQ. 0 ) THEN |
---|
268 | zalbe = forc_arr(11) |
---|
269 | zfswn = ( 1. - zalbe ) * zfsw |
---|
270 | ENDIF |
---|
271 | ! snowfall has to be adjusted |
---|
272 | zsfal = zsfal * ddtb / ts_forc |
---|
273 | ! ! sensitivity test |
---|
274 | ! ztair = ztair + 2. |
---|
275 | ! zsfal = zsfal + zsfal |
---|
276 | ! zcld = zcld - zcld / 10. |
---|
277 | ! |
---|
278 | !-----------------------------------------------------------------------------! |
---|
279 | ! 4) Compute some of the forcing fields if required |
---|
280 | !-----------------------------------------------------------------------------! |
---|
281 | ! |
---|
282 | !---------------------- |
---|
283 | ! 4.1) FLW computation |
---|
284 | !---------------------- |
---|
285 | i = 1 ; j = 1 |
---|
286 | IF ( ( forc_swi(2) .GT. 0 ) .AND. ( forc_swi(2) .LT. 99 ) ) THEN |
---|
287 | ze = zqair / ( 1. - zqair ) * zpres / 100. ! wvp in hPa |
---|
288 | & / ( 0.622 - zqair / ( 1. - zqair ) ) |
---|
289 | zta4 = ztair * ztair * ztair * ztair |
---|
290 | WRITE(numout,*) ' Efimova a' |
---|
291 | WRITE(numout,*) ' zta4 : ', zta4 |
---|
292 | WRITE(numout,*) ' ze : ', ze |
---|
293 | WRITE(numout,*) ' zcld : ', zcld |
---|
294 | ENDIF |
---|
295 | |
---|
296 | ! Efimova (61) and Key et al (96) |
---|
297 | !--------------------------------- |
---|
298 | IF ( forc_swi(2) == 1 ) THEN |
---|
299 | ! zflw = emig * stefan * zta4 * |
---|
300 | ! & ( 0.746 + 0.0066 * ze ) + ( 1 + 0.26 * zcld ) |
---|
301 | WRITE(numout,*) ' Efimova b' |
---|
302 | WRITE(numout,*) ' zta4 : ', zta4 |
---|
303 | WRITE(numout,*) ' ze : ', ze |
---|
304 | WRITE(numout,*) ' zcld : ', zcld |
---|
305 | ze = zqair / ( 1. - zqair ) * zpres / 100. ! wvp in hPa |
---|
306 | & / ( 0.622 - zqair / ( 1. - zqair ) ) |
---|
307 | zflw = 0.97 * stefan * zta4 * |
---|
308 | & ( 0.746 + 0.0066 * ze ) * ( 1. + 0.26 * zcld ) |
---|
309 | ENDIF |
---|
310 | |
---|
311 | ! Berliand |
---|
312 | !---------- |
---|
313 | IF ( forc_swi(2) == 2 ) THEN |
---|
314 | ! covrai(i,j)=sin(rlat*radian) |
---|
315 | ! alat = asin(covrai(i,j))/radian |
---|
316 | ! alat = ASIN ( SIN ( rlat * radian ) ) / radian ! Is it absolute or not |
---|
317 | clat = (90.0-rlat)/10.0 |
---|
318 | indx = 1+int(clat) |
---|
319 | zflw = stefan*zta4* |
---|
320 | & (1.0-(0.39-0.05*sqrt(ze))*(1.0-budyko(indx)* |
---|
321 | & zcld*zcld)) |
---|
322 | zcorr_fac_lw = 0.1 |
---|
323 | zflw = zflw + zflw * zcorr_fac_lw |
---|
324 | ENDIF ! forc_swi(2) |
---|
325 | |
---|
326 | ! Lab |
---|
327 | !---------- |
---|
328 | IF ( forc_swi(2) == 3 ) THEN |
---|
329 | zflw = emig * stefan * zta4 |
---|
330 | ENDIF |
---|
331 | |
---|
332 | !------------------------- |
---|
333 | ! 4.2) Albedo computation |
---|
334 | !------------------------- |
---|
335 | ! Shine |
---|
336 | !------- |
---|
337 | IF ( forc_swi(11) == 1 ) THEN |
---|
338 | CALL shine(tfsn, tfsg, t_su_b(ji), t_s_b(ji,1), |
---|
339 | & ht_i_b(ji), ht_s_b(ji), |
---|
340 | & alb_c, alb_o) |
---|
341 | zalbe = ( ( 1. - zcld ) * alb_c + zcld * alb_o ) |
---|
342 | ENDIF ! |
---|
343 | |
---|
344 | ! Albedo from observed values , ISPOL |
---|
345 | !------------------------------------- |
---|
346 | IF ( forc_swi(11) == 2 ) THEN |
---|
347 | IF ( ht_s_b(ji) .LE. 0.0 ) THEN |
---|
348 | zalbe = 0.50 !bare ice albedo |
---|
349 | ELSE |
---|
350 | IF ((t_su_b(ji) .GE. 273.15 ).OR.(t_s_b(ji,1) .GE. 273.15 )) |
---|
351 | & THEN |
---|
352 | zalbe = 0.65 ! melting snow albedo |
---|
353 | zalbe = 0.72 ! melting snow albedo (corrected from reprocessed data) |
---|
354 | ELSE |
---|
355 | zalbe = 0.80 ! dry snow albedo |
---|
356 | ENDIF |
---|
357 | ENDIF |
---|
358 | alb_c = zalbe |
---|
359 | alb_o = zalbe + 0.06 |
---|
360 | zalbe = ( ( 1. - zcld ) * alb_c + zcld * alb_o ) |
---|
361 | ENDIF ! |
---|
362 | |
---|
363 | ! net SW flux in case SW is read in a file |
---|
364 | !------------------------------------------ |
---|
365 | IF ( ( forc_swi(11) .GT. 0 ) .AND. ( forc_swi(11) .LT. 99 ) .AND. |
---|
366 | & ( forc_swi(1) .EQ. 0 ) ) zfswn = zfsw * ( 1.0 - zalbe ) |
---|
367 | |
---|
368 | !---------------------- |
---|
369 | ! 4.3) FSW computation |
---|
370 | !---------------------- |
---|
371 | IF ( ( forc_swi(1) .GT. 0 ) .AND. ( forc_swi(1) .LT. 99 ) ) THEN |
---|
372 | zeps0 = 1.d-13 |
---|
373 | dpi = 2*pi |
---|
374 | indaet = 1 |
---|
375 | |
---|
376 | dec = pdecli(indaet,num_d) * radian |
---|
377 | sdec = sin(dec) |
---|
378 | cdec = cos(dec) |
---|
379 | DO j = 1, 1 |
---|
380 | DO i = 1, 1 |
---|
381 | ! geometric factors |
---|
382 | !------------------- |
---|
383 | slat = SIN ( rlat * radian ) |
---|
384 | zps = slat*sdec |
---|
385 | zpc = COS(ASIN(slat))*cdec |
---|
386 | |
---|
387 | zljour = ACOS(-SIGN(one,zps) |
---|
388 | & * MIN(one,SIGN(one,zps)*(zps/zpc))) |
---|
389 | |
---|
390 | dws = (2.0*zljour)/REAL(nintsr) |
---|
391 | zlmidi = ASIN(( zps +zpc ))/radian |
---|
392 | zalcnq = 0.0 |
---|
393 | DO k = 1, nintsr |
---|
394 | ws(k) = zljour-(REAL(k)-0.5)*dws |
---|
395 | zmue(k) = MAX(c_zero,zps+zpc*COS(ws(k))) |
---|
396 | zalcnp(k) = 0.05/(1.1*zmue(k)**1.4+0.15) |
---|
397 | zalcnq = zalcnp(k) |
---|
398 | END DO |
---|
399 | zalcnq = zalcnq/MAX(2.0*zljour,zeps0) |
---|
400 | zmudum = 0.4 |
---|
401 | |
---|
402 | ! Irradiance |
---|
403 | !------------ |
---|
404 | ze = zqair / ( 1. - zqair ) * zpres ! wvp in Pa |
---|
405 | & / ( 0.622 - zqair / ( 1. - zqair ) ) |
---|
406 | ztc = ztair - 273.15 |
---|
407 | zesw = 611.*EXP(17.269*ztc/(ztc+237.3)) ! Pa |
---|
408 | ztdew = ze / zesw ! no units |
---|
409 | |
---|
410 | !------- |
---|
411 | ! Shine |
---|
412 | !------- |
---|
413 | IF ( forc_swi(1) == 1 ) THEN |
---|
414 | !----------------- |
---|
415 | ! DAILY time step |
---|
416 | !----------------- |
---|
417 | IF ( ddtb .EQ. 86400.0 ) THEN |
---|
418 | frsdrg = 0.0 |
---|
419 | frsdfg = 0.0 |
---|
420 | frsdro = 0.0 |
---|
421 | frsdfo = 0.0 |
---|
422 | ! opt_dept = 16.297 ! to put in the namelist |
---|
423 | ! opt_dept = 5.6000 ! to put in the namelist |
---|
424 | DO k = 1 , nintsr ! integrate over the whole day |
---|
425 | frsdrg = frsdrg+dws* ! clear |
---|
426 | & (zsolar*zmue(k)*zmue(k)*(1.0-alb_c))/ |
---|
427 | & (1.2*zmue(k)+(1.0+zmue(k))*ze*1.0e-05+0.0455) |
---|
428 | frsdfg = frsdfg+dws* ! overcast |
---|
429 | & ( ( 53.5 + 1274.5*zmue(k) ) * SQRT(zmue(k) ) |
---|
430 | & * ( 1.0 - 0.996*alb_o ) ) |
---|
431 | & / (1.0+0.139*(1.0-0.9435*alb_o)*opt_dept) |
---|
432 | END DO |
---|
433 | ! net solar heat flux (1-a)FSW |
---|
434 | zfswn = ( ( 1.0 - zcld ) * frsdrg + zcld * frsdfg ) / |
---|
435 | & dpi |
---|
436 | ENDIF ! ddtb |
---|
437 | |
---|
438 | !-------------------- |
---|
439 | ! SUBDAILY time step |
---|
440 | !-------------------- |
---|
441 | IF ( ddtb .LT. 86400.00 ) THEN |
---|
442 | ! WRITE(numout,*) ' Shine formula, subdaily time step ' |
---|
443 | ! WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' |
---|
444 | |
---|
445 | ! Solar angle |
---|
446 | zdecl = pdecli(indaet,num_d) * radian |
---|
447 | zsdec = SIN(zdecl) |
---|
448 | zcdec = COS(zdecl) |
---|
449 | zslat = SIN ( rlat * radian ) |
---|
450 | zps = zslat*zsdec |
---|
451 | zpc = COS(ASIN(slat))*cdec |
---|
452 | zhour = ( doy - FLOAT(num_d) ) * 24. |
---|
453 | |
---|
454 | ! Hour angle |
---|
455 | zljour = ACOS(-SIGN(one,zps) |
---|
456 | & * MIN(one,SIGN(one,zps)*(zps/zpc))) |
---|
457 | zhourang = - pi + 2.*pi / 24. * zhour |
---|
458 | ! Cosine of solar angle |
---|
459 | zcosz = 0.0 |
---|
460 | IF ( ( zhourang .GT. -zljour ) .AND. |
---|
461 | & ( zhourang .LT. zljour ) ) |
---|
462 | & zcosz = MAX(0.0,zps+zpc*COS(zhourang)) |
---|
463 | zcosz2 = zcosz * zcosz |
---|
464 | |
---|
465 | ! Irradiance |
---|
466 | zqsr_clear = zsolar * zcosz2 |
---|
467 | & / ( 1.2 * zcosz |
---|
468 | & + ( 1.0 + zcosz ) * ze * 1.0e-5 |
---|
469 | & + 0.0455 ) |
---|
470 | zqsr_cloud = ( 53.5 + 1274.5 * zcosz ) * |
---|
471 | & SQRT( zcosz ) / ( 1.0 + 0.139 * |
---|
472 | & ( 1.0 - 0.9435 * zalbe ) * opt_dept ) |
---|
473 | zfsw = ( 1. - zcld ) * zqsr_clear |
---|
474 | & + zcld * zqsr_cloud |
---|
475 | zfswn = ( 1. - zcld ) * ( 1. - alb_c ) * zqsr_clear |
---|
476 | & + zcld * ( 1. - alb_o ) * zqsr_cloud |
---|
477 | |
---|
478 | WRITE(numout,*) ' zqsr_clear : ', zqsr_clear |
---|
479 | WRITE(numout,*) ' zqsr_cloud : ', zqsr_cloud |
---|
480 | |
---|
481 | ENDIF ! ddtb |
---|
482 | |
---|
483 | ENDIF ! Shine |
---|
484 | |
---|
485 | !---------- |
---|
486 | ! Zillmann |
---|
487 | !---------- |
---|
488 | IF ( forc_swi(1) == 2 ) THEN |
---|
489 | IF ( ddtb .EQ. 86400.0 ) THEN |
---|
490 | frsdtg = 0.0 |
---|
491 | frsdto = 0.0 |
---|
492 | DO k = 1 , nintsr |
---|
493 | albo = (1.0-zcld)*zalcnp(k)+zcld*zalbe |
---|
494 | frsdtg = frsdtg+dws*(1.0-zalbe)* |
---|
495 | & (zsolar*zmue(k)*zmue(k))/ |
---|
496 | & ((zmue(k)+2.7)*ze*1.0e-05+1.085*zmue(k)+0.10) |
---|
497 | END DO |
---|
498 | ! rewrite this next line correctly |
---|
499 | zfswn = 0.9*min(one,(1-.62*zcld+.0019*zcld))* |
---|
500 | & frsdtg/dpi |
---|
501 | ENDIF |
---|
502 | ENDIF ! |
---|
503 | END DO |
---|
504 | END DO |
---|
505 | ENDIF ! forc_swi(1) |
---|
506 | |
---|
507 | IF ( forc_swi(3) == 1 ) THEN |
---|
508 | !----------------- |
---|
509 | ! 4.2 PAR formula |
---|
510 | !----------------- |
---|
511 | zpar = 0.43 |
---|
512 | ENDIF |
---|
513 | ! |
---|
514 | !-----------------------------------------------------------------------------! |
---|
515 | ! 5) Case of prescribed values |
---|
516 | !-----------------------------------------------------------------------------! |
---|
517 | ! |
---|
518 | IF ( forc_swi(1) .EQ. 99 ) zfsw = forc_val(1) * forc_uni(1) |
---|
519 | IF ( forc_swi(2) .EQ. 99 ) zflw = forc_val(2) * forc_uni(2) |
---|
520 | IF ( forc_swi(3) .EQ. 99 ) zpar = forc_val(3) * forc_uni(3) |
---|
521 | IF ( forc_swi(4) .EQ. 99 ) ztair = forc_val(4) * forc_uni(4) |
---|
522 | IF ( forc_swi(5) .EQ. 99 ) zpres = forc_val(5) * forc_uni(5) |
---|
523 | IF ( forc_swi(6) .EQ. 99 ) zqair = forc_val(6) * forc_uni(6) |
---|
524 | IF ( forc_swi(7) .EQ. 99 ) zwspd = forc_val(7) * forc_uni(7) |
---|
525 | IF ( forc_swi(8) .EQ. 99 ) zcld = forc_val(8) * forc_uni(8) |
---|
526 | IF ( forc_swi(9) .EQ. 99 ) zfoce = forc_val(9) * forc_uni(9) |
---|
527 | IF ( forc_swi(10).EQ. 99 ) zsfal = forc_val(10) * forc_uni(10) |
---|
528 | IF ( forc_swi(11).EQ. 99 ) THEN |
---|
529 | zalbe = forc_val(11) * forc_uni(11) |
---|
530 | zfswn = ( 1. - zalbe ) * zfsw |
---|
531 | ENDIF |
---|
532 | ! |
---|
533 | !-----------------------------------------------------------------------------! |
---|
534 | ! 6) Temporary dirty plug |
---|
535 | !-----------------------------------------------------------------------------! |
---|
536 | ! |
---|
537 | |
---|
538 | fsolgb(ji) = zfswn |
---|
539 | ratbqb(ji) = zflw |
---|
540 | ! par is not yet included, should be |
---|
541 | tabqb(ji) = ztair |
---|
542 | psbqb(ji) = zpres |
---|
543 | qabqb(ji) = zqair |
---|
544 | vabqb(ji) = zwspd |
---|
545 | cldqb(ji) = zcld |
---|
546 | oce_flx = zfoce |
---|
547 | hnpbqb(ji) = zsfal |
---|
548 | albgb(ji) = zalbe |
---|
549 | ze = zqair / ( 1. - zqair ) * zpres ! wvp in Pa |
---|
550 | & / ( 0.622 - zqair / ( 1. - zqair ) ) |
---|
551 | ztc = ztair - 273.15 |
---|
552 | zesw = 611.*EXP(17.269*ztc/(ztc+237.3)) ! Pa |
---|
553 | ztdew = ze / zesw ! no units |
---|
554 | tdewb(ji) = ztdew |
---|
555 | |
---|
556 | IF ( ln_write_forc ) THEN |
---|
557 | WRITE(numout,*) ' Forcing fields ... ' |
---|
558 | WRITE(numout,*) |
---|
559 | WRITE(numout,*) ' fsolgb : ', fsolgb(ji) |
---|
560 | WRITE(numout,*) ' ratbqg : ', ratbqb(ji) |
---|
561 | WRITE(numout,*) ' tabqb : ', tabqb(ji) |
---|
562 | WRITE(numout,*) ' psbqb : ', psbqb(ji) |
---|
563 | WRITE(numout,*) ' qabqb : ', qabqb(ji) |
---|
564 | WRITE(numout,*) ' vabqb : ', vabqb(ji) |
---|
565 | WRITE(numout,*) ' cldqb : ', cldqb(ji) |
---|
566 | WRITE(numout,*) ' oce_flx : ', oce_flx |
---|
567 | WRITE(numout,*) ' hnpbqb : ', hnpbqb(ji) |
---|
568 | WRITE(numout,*) ' albgb : ', albgb(ji) |
---|
569 | WRITE(numout,*) |
---|
570 | ENDIF |
---|
571 | |
---|
572 | |
---|
573 | !------------------------------------------------------------------------------ |
---|
574 | ! end of forcing |
---|
575 | RETURN |
---|
576 | END |
---|