[6] | 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 |
---|