Changeset 833 for trunk/NEMO/OPA_SRC/SBC/flxblk.F90
- Timestamp:
- 2008-03-07T14:51:35+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SBC/flxblk.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SBC/flxblk.F90
r789 r833 4 4 !! Ocean forcing: bulk thermohaline forcing of the ocean (or ice) 5 5 !!===================================================================== 6 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily6 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily 7 7 !!---------------------------------------------------------------------- 8 8 !! 'key_flx_bulk_monthly' or MONTHLY bulk … … 25 25 USE albedo 26 26 USE prtctl ! Print control 27 27 #if defined key_lim3 28 USE par_ice 29 USE ice 30 #elif defined key_lim2 31 USE ice_2 32 #endif 28 33 IMPLICIT NONE 29 34 PRIVATE … … 79 84 80 85 CONTAINS 86 #if defined key_lim3 87 88 SUBROUTINE flx_blk( psst ) 89 !!--------------------------------------------------------------------------- 90 !! *** ROUTINE flx_blk *** 91 !! 92 !! ** Purpose : Computation of the heat fluxes at ocean and snow/ice 93 !! surface the solar heat at ocean and snow/ice surfaces and the 94 !! sensitivity of total heat fluxes to the SST variations 95 !! 96 !! ** Method : The flux of heat at the ice and ocean surfaces are derived 97 !! from semi-empirical ( or bulk ) formulae which relate the flux to 98 !! the properties of the surface and of the lower atmosphere. Here, we 99 !! follow the work of Oberhuber, 1988 100 !! 101 !! ** Action : call flx_blk_albedo to compute ocean and ice albedo 102 !! computation of snow precipitation 103 !! computation of solar flux at the ocean and ice surfaces 104 !! computation of the long-wave radiation for the ocean and sea/ice 105 !! computation of turbulent heat fluxes over water and ice 106 !! computation of evaporation over water 107 !! computation of total heat fluxes sensitivity over ice (dQ/dT) 108 !! computation of latent heat flux sensitivity over ice (dQla/dT) 109 !! 110 !! History : 111 !! 8.0 ! 97-06 (Louvain-La-Neuve) Original code 112 !! 8.5 ! 02-09 (C. Ethe , G. Madec ) F90: Free form and module 113 !!---------------------------------------------------------------------- 114 !! * Arguments 115 REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) :: & 116 & psst ! Sea Surface Temperature 117 118 !! * Local variables 119 INTEGER :: & 120 ji, jj, jl, jt , & ! dummy loop indices 121 indaet , & ! = -1, 0, 1 for odd, normal and leap years resp. 122 iday , & ! integer part of day 123 indxb , & ! index for budyko coefficient 124 indxc ! index for cloud depth coefficient 125 126 REAL(wp) :: & 127 zalat , zclat , & ! latitude in degrees 128 zmt1, zmt2, zmt3 , & ! tempory air temperatures variables 129 ztatm3, ztatm4 , & ! power 3 and 4 of air temperature 130 z4tatm3 , & ! 4 * ztatm3 131 zcmue , & ! cosine of local solar altitude 132 zcmue2 , & ! root of zcmue1 133 zscmue , & ! square-root of zcmue1 134 zpcmue , & ! zcmue1**1.4 135 zdecl , & ! solar declination 136 zsdecl , zcdecl , & ! sine and cosine of solar declination 137 zalbo , & ! albedo of sea-water 138 zalbi , & ! albedo of ice 139 ztamr , & ! air temperature minus triple point of water (rtt) 140 ztaevbk , & ! part of net longwave radiation 141 zevi , zevo , & ! vapour pressure of ice and ocean 142 zind1,zind2,zind3 , & ! switch for testing the values of air temperature 143 zinda , & ! switch for testing the values of sea ice cover 144 zpis2 , & ! pi / 2 145 z2pi ! 2 * pi 146 147 REAL(wp) :: & 148 zxday , & ! day of year 149 zdist , & ! distance between the sun and the earth during the year 150 zdaycor , & ! corr. factor to take into account the variation of 151 ! ! zday when calc. the solar rad. 152 zesi, zeso , & ! vapour pressure of ice and ocean at saturation 153 zesi2 , & ! root of zesi 154 zqsato , & ! humidity close to the ocean surface (at saturation) 155 zqsati , & ! humidity close to the ice surface (at saturation) 156 zqsati2 , & ! root of zqsati 157 zdesidt , & ! derivative of zesi, function of ice temperature 158 zdteta , & ! diff. betw. sst and air temperature 159 zdeltaq , & ! diff. betw. spec. hum. and hum. close to the surface 160 ztvmoy, zobouks , & ! tempory scalars 161 zpsims, zpsihs, zpsils, zobouku, zxins, zpsimu , & 162 zpsihu, zpsilu, zstab,zpsim, zpsih, zpsil , & 163 zvatmg, zcmn, zchn, zcln, zcmcmn, zdenum , & 164 zdtetar, ztvmoyr, zlxins, zcmn2, zchcm, zclcm , zcoef 165 166 REAL(wp) :: & 167 zrhova , & ! air density per wind speed 168 zcsho , zcleo , & ! transfer coefficient over ocean 169 zcshi , zclei , & ! transfer coefficient over ice-free 170 zrhovacleo , & ! air density per wind speed per transfer coef. 171 zrhovacsho, zrhovaclei, zrhovacshi, & 172 ztice3 , & ! power 3 of ice temperature 173 zticemb, zticemb2 , & ! tempory air temperatures variables 174 zdqlw_ice , & ! sensitivity of long-wave flux over ice 175 zdqsb_ice , & ! sensitivity of sensible heat flux over ice 176 zdqla_ice , & ! sensitivity of latent heat flux over ice 177 zdl, zdr ! fractionnal part of latitude 178 REAL(wp), DIMENSION(jpi,jpj) :: & 179 zpatm , & ! atmospheric pressure 180 zqatm , & ! specific humidity 181 zes , & ! vapour pressure at saturation 182 zev, zevsqr , & ! vapour pressure and his square-root 183 zrhoa , & ! air density 184 ztatm , & ! air temperature in Kelvins 185 zfrld , & ! fraction of sea ice cover 186 zcatm1 , & ! fraction of cloud 187 zcldeff ! correction factor to account cloud effect 188 REAL(wp), DIMENSION(jpi,jpj) :: & 189 zalbocsd , & ! albedo of ocean 190 zalboos , & ! albedo of ocean under overcast sky 191 zalbomu , & ! albedo of ocean when zcmue is 0.4 192 zqsro , & ! solar radiation over ocean 193 zqsrics , & ! solar radiation over ice under clear sky 194 zqsrios , & ! solar radiation over ice under overcast sky 195 zcldcor , & ! cloud correction 196 zlsrise, zlsset , & ! sunrise and sunset 197 zlmunoon , & ! local noon solar altitude 198 zdlha , & ! length of the ninstr segments of the solar day 199 zps , & ! sine of latitude per sine of solar decli. 200 zpc ! cosine of latitude per cosine of solar decli. 201 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 202 zalbics , & ! albedo of ice under clear sky 203 zalbios ! albedo of ice under overcast sky 204 205 REAL(wp), DIMENSION(jpi,jpj) :: & 206 zqlw_oce , & ! long-wave heat flux over ocean 207 zqla_oce , & ! latent heat flux over ocean 208 zqsb_oce ! sensible heat flux over ocean 209 210 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 211 zqlw_ice , & ! long-wave heat flux over ice 212 zqla_ice , & ! latent heat flux over ice 213 zqsb_ice ! sensible heat flux over ice 214 215 REAL(wp), DIMENSION(jpi,jpj,jpintsr) :: & 216 zlha , & ! local hour angle 217 zalbocs , & ! tempory var. of ocean albedo under clear sky 218 zsqsro , & ! tempory var. of solar rad. over ocean 219 zsqsrics , & ! temp. var. of solar rad. over ice under clear sky 220 zsqsrios ! temp. var. of solar rad. over ice under overcast sky 221 !!--------------------------------------------------------------------- 222 223 ! Determine cloud optical depths as a function of latitude (Chou et al., 1981). 224 ! and the correction factor for taking into account the effect of clouds 225 !------------------------------------------------------ 226 IF( lbulk_init ) THEN 227 DO jj = 1, jpj 228 DO ji = 1 , jpi 229 zalat = ( 90.e0 - ABS( gphit(ji,jj) ) ) / 5.e0 230 zclat = ( 95.e0 - gphit(ji,jj) ) / 10.e0 231 indxb = 1 + INT( zalat ) 232 ! correction factor to account for the effect of clouds 233 sbudyko(ji,jj) = budyko(indxb) 234 indxc = 1 + INT( zclat ) 235 zdl = zclat - INT( zclat ) 236 zdr = 1.0 - zdl 237 stauc(ji,jj) = zdr * tauco( indxc ) + zdl * tauco( indxc + 1 ) 238 END DO 239 END DO 240 IF( nleapy == 1 ) THEN 241 yearday = 366.e0 242 ELSE IF( nleapy == 0 ) THEN 243 yearday = 365.e0 244 ELSEIF( nleapy == 30) THEN 245 yearday = 360.e0 246 ENDIF 247 lbulk_init = .FALSE. 248 ENDIF 249 250 zqlw_oce(:,:) = 0.e0 251 zqla_oce(:,:) = 0.e0 252 zqsb_oce(:,:) = 0.e0 253 zqlw_ice(:,:,:) = 0.e0 254 zqla_ice(:,:,:) = 0.e0 255 zqsb_ice(:,:,:) = 0.e0 256 257 zpis2 = rpi / 2. 258 z2pi = 2. * rpi 259 260 !CDIR NOVERRCHK 261 DO jj = 1, jpj 262 !CDIR NOVERRCHK 263 DO ji = 1, jpi 264 265 ztatm (ji,jj) = 273.15 + tatm (ji,jj) ! air temperature in Kelvins 266 zcatm1(ji,jj) = 1.0 - catm (ji,jj) ! fractional cloud cover 267 zfrld (ji,jj) = 1.0 - freeze(ji,jj) ! fractional sea ice cover 268 zpatm(ji,jj) = 101000. ! pressure 269 270 ! Computation of air density, obtained from the equation of state for dry air. 271 zrhoa(ji,jj) = zpatm(ji,jj) / ( 287.04 * ztatm(ji,jj) ) 272 273 ! zes : Saturation water vapour 274 ztamr = ztatm(ji,jj) - rtt 275 zmt1 = SIGN( 17.269, ztamr ) 276 zmt2 = SIGN( 21.875, ztamr ) 277 zmt3 = SIGN( 28.200, -ztamr ) 278 zes(ji,jj) = 611.0 * EXP ( ABS( ztamr ) * MIN ( zmt1, zmt2 ) & 279 & / ( ztatm(ji,jj) - 35.86 + MAX( zzero, zmt3 ) ) ) 280 281 ! zev : vapour pressure (hatm is relative humidity) 282 zev(ji,jj) = hatm(ji,jj) * zes(ji,jj) 283 ! square-root of vapour pressure 284 !CDIR NOVERRCHK 285 zevsqr(ji,jj) = SQRT( zev(ji,jj) * 0.01 ) 286 ! zqapb : specific humidity 287 zqatm(ji,jj) = 0.622 * zev(ji,jj) / ( zpatm(ji,jj) - 0.378 * zev(ji,jj) ) 288 289 290 !---------------------------------------------------- 291 ! Computation of snow precipitation (Ledley, 1985) | 292 !---------------------------------------------------- 293 294 zmt1 = 253.0 - ztatm(ji,jj) 295 zmt2 = ( 272.0 - ztatm(ji,jj) ) / 38.0 296 zmt3 = ( 281.0 - ztatm(ji,jj) ) / 18.0 297 zind1 = MAX( zzero, SIGN( zone, zmt1 ) ) 298 zind2 = MAX( zzero, SIGN( zone, zmt2 ) ) 299 zind3 = MAX( zzero, SIGN( zone, zmt3 ) ) 300 ! total precipitation 301 tprecip(ji,jj) = watm(ji,jj) 302 ! solid (snow) precipitation 303 sprecip(ji,jj) = tprecip(ji,jj) * & 304 & ( zind1 & 305 & + ( 1.0 - zind1 ) * ( zind2 * ( 0.5 + zmt2 ) + ( 1.0 - zind2 ) * zind3 * zmt3 ) ) 306 END DO 307 END DO 308 309 !---------------------------------------------------------- 310 ! Computation of albedo (need to calculates heat fluxes)| 311 !----------------------------------------------------------- 312 313 CALL flx_blk_albedo( zalbios, zalboos, zalbics, zalbomu ) 314 315 !------------------------------------- 316 ! Computation of solar irradiance. | 317 !---------------------------------------- 318 indaet = 1 319 ! compution of the day of the year at which the fluxes have to be calculate 320 !--The date corresponds to the middle of the time step. 321 zxday=nday_year + rdtbs2/rday 322 323 iday = INT( zxday ) 324 325 IF(ln_ctl) CALL prt_ctl_info('declin : iday ', ivar1=iday, clinfo2=' nfbulk= ', ivar2=nfbulk) 326 327 ! computation of the solar declination, his sine and his cosine 328 CALL flx_blk_declin( indaet, iday, zdecl ) 329 330 zdecl = zdecl * rad 331 zsdecl = SIN( zdecl ) 332 zcdecl = COS( zdecl ) 333 334 ! correction factor added for computation of shortwave flux to take into account the variation of 335 ! the distance between the sun and the earth during the year (Oberhuber 1988) 336 zdist = zxday * z2pi / yearday 337 zdaycor = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 338 339 !CDIR NOVERRCHK 340 DO jj = 1, jpj 341 !CDIR NOVERRCHK 342 DO ji = 1, jpi 343 ! product of sine of latitude and sine of solar declination 344 zps (ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl 345 ! product of cosine of latitude and cosine of solar declination 346 zpc (ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl 347 ! computation of the both local time of sunrise and sunset 348 zlsrise (ji,jj) = ACOS( - SIGN( zone, zps(ji,jj) ) * MIN( zone, SIGN( zone, zps(ji,jj) ) & 349 & * ( zps(ji,jj) / zpc(ji,jj) ) ) ) 350 zlsset (ji,jj) = - zlsrise(ji,jj) 351 ! dividing the solar day into jpintsr segments of length zdlha 352 zdlha (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jpintsr ) 353 ! computation of the local noon solar altitude 354 zlmunoon(ji,jj) = ASIN ( ( zps(ji,jj) + zpc(ji,jj) ) ) / rad 355 356 ! cloud correction taken from Reed (1977) (imposed lower than 1) 357 zcldcor (ji,jj) = MIN( zone, ( 1.e0 - 0.62 * catm(ji,jj) + 0.0019 * zlmunoon(ji,jj) ) ) 358 END DO 359 END DO 360 361 ! Computation of solar heat flux at each time of the day between sunrise and sunset. 362 ! We do this to a better optimisation of the code 363 !------------------------------------------------------ 364 DO jl = 1, jpl 365 366 !CDIR NOVERRCHK 367 DO jt = 1, jpintsr 368 zcoef = FLOAT( jt ) - 0.5 369 !CDIR NOVERRCHK 370 DO jj = 1, jpj 371 !CDIR NOVERRCHK 372 DO ji = 1, jpi 373 ! local hour angle 374 zlha (ji,jj,jt) = COS ( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) ) 375 376 ! cosine of local solar altitude 377 zcmue = MAX ( zzero , zps(ji,jj) + zpc(ji,jj) * zlha (ji,jj,jt) ) 378 zcmue2 = 1368.0 * zcmue * zcmue 379 zscmue = SQRT ( zcmue ) 380 zpcmue = zcmue**1.4 381 ! computation of sea-water albedo (Payne, 1972) 382 zalbocs(ji,jj,jt) = 0.05 / ( 1.1 * zpcmue + 0.15 ) 383 zalbo = zcatm1(ji,jj) * zalbocs(ji,jj,jt) + catm(ji,jj) * zalboos(ji,jj) 384 ! solar heat flux absorbed at ocean surfaces (Zillman, 1972) 385 zevo = zev(ji,jj) * 1.0e-05 386 zsqsro(ji,jj,jt) = ( 1.0 - zalbo ) * zdlha(ji,jj) * zcmue2 & 387 / ( ( zcmue + 2.7 ) * zevo + 1.085 * zcmue + 0.10 ) 388 ! solar heat flux absorbed at sea/ice surfaces 389 ! Formulation of Shine and Crane, 1984 adapted for high albedo surfaces 390 391 ! For clear sky 392 zevi = zevo 393 zalbi = zalbics(ji,jj,jl) 394 zsqsrics(ji,jj,jt) = ( 1.0 - zalbi ) * zdlha(ji,jj) * zcmue2 & 395 & / ( ( 1.0 + zcmue ) * zevi + 1.2 * zcmue + 0.0455 ) 396 397 ! For overcast sky 398 zalbi = zalbios(ji,jj,jl) 399 zsqsrios(ji,jj,jt) = zdlha(ji,jj) * & 400 & ( ( 53.5 + 1274.5 * zcmue ) * zscmue * ( 1.0 - 0.996 * zalbi ) ) & 401 & / ( 1.0 + 0.139 * stauc(ji,jj) * ( 1.0 - 0.9435 * zalbi ) ) 402 END DO 403 END DO 404 END DO 405 406 407 ! Computation of daily (between sunrise and sunset) solar heat flux absorbed 408 ! at the ocean and snow/ice surfaces. 409 !-------------------------------------------------------------------- 410 411 zalbocsd(:,:) = 0.e0 412 zqsro (:,:) = 0.e0 413 zqsrics (:,:) = 0.e0 414 zqsrios (:,:) = 0.e0 415 416 DO jt = 1, jpintsr 417 # if defined key_vectopt_loop 418 DO ji = 1, jpij 419 zalbocsd(ji,1) = zalbocsd(ji,1) + zdlha (ji,1) * zalbocs(ji,1,jt) & 420 & / MAX( 2.0 * zlsrise(ji,1) , zeps0 ) 421 zqsro (ji,1) = zqsro (ji,1) + zsqsro (ji,1,jt) 422 zqsrics (ji,1) = zqsrics (ji,1) + zsqsrics(ji,1,jt) 423 zqsrios (ji,1) = zqsrios (ji,1) + zsqsrios(ji,1,jt) 424 END DO 425 # else 426 DO jj = 1, jpj 427 DO ji = 1, jpi 428 zalbocsd(ji,jj) = zalbocsd(ji,jj) + zdlha(ji,jj) * zalbocs(ji,jj,jt) & 429 & / MAX( 2.0 * zlsrise(ji,jj) , zeps0 ) 430 zqsro (ji,jj) = zqsro (ji,jj) + zsqsro (ji,jj,jt) 431 zqsrics(ji,jj) = zqsrics (ji,jj) + zsqsrics(ji,jj,jt) 432 zqsrios(ji,jj) = zqsrios (ji,jj) + zsqsrios(ji,jj,jt) 433 END DO 434 END DO 435 # endif 436 END DO 437 438 DO jj = 1, jpj 439 DO ji = 1, jpi 440 441 !------------------------------------------- 442 ! Computation of shortwave radiation. 443 !------------------------------------------- 444 445 ! the solar heat flux absorbed at ocean and snow/ice surfaces 446 !------------------------------------------------------------ 447 448 ! For snow/ice 449 qsr_ice(ji,jj,jl) = ( zcatm1(ji,jj) * zqsrics(ji,jj) + catm(ji,jj) * zqsrios(ji,jj) ) / z2pi 450 451 ! Taking into account the ellipsity of the earth orbit 452 !----------------------------------------------------- 453 454 qsr_ice(ji,jj,jl) = qsr_ice(ji,jj,jl) * zdaycor 455 !--------------------------------------------------------------------------- 456 ! Computation of long-wave radiation ( Berliand 1952 ; all latitudes ) 457 !--------------------------------------------------------------------------- 458 459 ! tempory variables 460 ztatm3 = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj) 461 ztatm4 = ztatm3 * ztatm(ji,jj) 462 z4tatm3 = 4. * ztatm3 463 zcldeff(ji,jj) = 1.0 - sbudyko(ji,jj) * catm(ji,jj) * catm(ji,jj) 464 ztaevbk = ztatm4 * zcldeff(ji,jj) * ( 0.39 - 0.05 * zevsqr(ji,jj) ) 465 466 ! Long-Wave for Ice 467 !---------------------- 468 zqlw_ice(ji,jj,jl) = - emic * stefan * ( ztaevbk + z4tatm3 * ( t_su(ji,jj,jl) - ztatm(ji,jj) ) ) 469 470 END DO !ji 471 END DO !jj 472 473 END DO !jl 474 475 DO jj = 1, jpj 476 DO ji = 1, jpi 477 478 ! fraction of net shortwave radiation which is not absorbed in the 479 ! thin surface layer and penetrates inside the ice cover 480 ! ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 481 !------------------------------------------------------------------ 482 483 fr1_i0(ji,jj) = 0.18 * zcatm1(ji,jj) + 0.35 * catm(ji,jj) 484 fr2_i0(ji,jj) = 0.82 * zcatm1(ji,jj) + 0.65 * catm(ji,jj) 485 486 ! the solar heat flux absorbed at ocean and snow/ice surfaces 487 !------------------------------------------------------------ 488 ! For ocean 489 qsr_oce(ji,jj) = srgamma * zcldcor(ji,jj) * zqsro(ji,jj) / z2pi 490 zinda = SIGN( zone , -( -0.5 - zfrld(ji,jj) ) ) 491 zinda = 1.0 - MAX( zzero , zinda ) 492 qsr_oce(ji,jj) = ( 1.- zinda ) * qsr_oce(ji,jj) 493 494 ! Taking into account the ellipsity of the earth orbit 495 !----------------------------------------------------- 496 qsr_oce(ji,jj) = qsr_oce(ji,jj) * zdaycor 497 498 !--------------------------------------------------------------------------- 499 ! Computation of long-wave radiation ( Berliand 1952 ; all latitudes ) 500 !--------------------------------------------------------------------------- 501 502 ! tempory variables 503 ztatm3 = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj) 504 ztatm4 = ztatm3 * ztatm(ji,jj) 505 z4tatm3 = 4. * ztatm3 506 zcldeff(ji,jj) = 1.0 - sbudyko(ji,jj) * catm(ji,jj) * catm(ji,jj) 507 ztaevbk = ztatm4 * zcldeff(ji,jj) * ( 0.39 - 0.05 * zevsqr(ji,jj) ) 508 509 ! Long-Wave for Ocean 510 !----------------------- 511 zqlw_oce(ji,jj) = - emic * stefan * ( ztaevbk + z4tatm3 * ( psst (ji,jj) - ztatm(ji,jj) ) ) 512 513 END DO 514 END DO 515 516 !---------------------------------------- 517 ! Computation of turbulent heat fluxes ( Latent and sensible ) 518 !---------------------------------------- 519 !CDIR NOVERRCHK 520 DO jj = 2 , jpjm1 521 !CDIR NOVERRCHK 522 DO ji = 1, jpi 523 524 ! Turbulent heat fluxes over water 525 !---------------------------------- 526 527 ! zeso : vapour pressure at saturation of ocean 528 ! zqsato : humidity close to the ocean surface (at saturation) 529 zeso = 611.0 * EXP ( 17.2693884 * ( psst(ji,jj) - rtt ) * tmask(ji,jj,1) / ( psst(ji,jj) - 35.86 ) ) 530 zqsato = ( 0.622 * zeso ) / ( zpatm(ji,jj) - 0.378 * zeso ) 531 532 ! Drag coefficients from Large and Pond (1981,1982) 533 !--------------------------------------------------- 534 535 ! Stability parameters 536 zdteta = psst(ji,jj) - ztatm(ji,jj) 537 zdeltaq = zqatm(ji,jj) - zqsato 538 ztvmoy = ztatm(ji,jj) * ( 1. + 2.2e-3 * ztatm(ji,jj) * zqatm(ji,jj) ) 539 zdenum = MAX( vatm(ji,jj) * vatm(ji,jj) * ztvmoy, zeps ) 540 zdtetar = zdteta / zdenum 541 ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum 542 543 ! For stable atmospheric conditions 544 zobouks = -70.0 * 10. * ( zdtetar + 3.2e-3 * ztvmoyr ) 545 zobouks = MAX( zzero , zobouks ) 546 zpsims = -7.0 * zobouks 547 zpsihs = zpsims 548 zpsils = zpsims 549 550 ! For unstable atmospheric conditions 551 zobouku = -100.0 * 10.0 * ( zdtetar + 2.2e-3 * ztvmoyr ) 552 zobouku = MIN( zzero , zobouku ) 553 zxins = ( 1. - 16. * zobouku )**0.25 554 zlxins = LOG( ( 1. + zxins * zxins ) / 2. ) 555 zpsimu = 2. * LOG( ( 1 + zxins ) / 2. ) + zlxins - 2. * ATAN( zxins ) + zpis2 556 zpsihu = 2. * zlxins 557 zpsilu = zpsihu 558 559 ! computation of intermediate values 560 zstab = MAX( zzero , SIGN( zone , zdteta ) ) 561 zpsim = zstab * zpsimu + (1.0 - zstab ) * zpsims 562 zpsih = zstab * zpsihu + (1.0 - zstab ) * zpsihs 563 zpsil = zpsih 564 565 zvatmg = MAX( 0.032 * 1.5e-3 * vatm(ji,jj) * vatm(ji,jj) / grav, zeps ) 566 567 zcmn = vkarmn / LOG ( 10. / zvatmg ) 568 zcmn2 = zcmn * zcmn 569 zchn = 0.0327 * zcmn 570 zcln = 0.0346 * zcmn 571 zcmcmn = 1 / ( 1 - zcmn * zpsim / vkarmn ) 572 zchcm = zcmcmn / ( 1 - zchn * zpsih / ( vkarmn * zcmn ) ) 573 zclcm = zchcm 574 575 576 ! Transfer cofficient zcsho and zcleo over ocean according to Large and Pond (1981,1982) 577 !-------------------------------------------------------------- 578 zcsho = zchn * zchcm 579 zcleo = zcln * zclcm 580 581 582 ! Computation of sensible and latent fluxes over Ocean 583 !---------------------------------------------------------------- 584 585 ! computation of intermediate values 586 zrhova = zrhoa(ji,jj) * vatm(ji,jj) 587 zrhovacsho = zrhova * zcsho 588 zrhovacleo = zrhova * zcleo 589 590 ! sensible heat flux 591 zqsb_oce(ji,jj) = zrhovacsho * 1004.0 * ( psst(ji,jj) - ztatm(ji,jj) ) 592 593 ! latent heat flux 594 zqla_oce(ji,jj) = MAX(0.e0, zrhovacleo * 2.5e+06 * ( zqsato - zqatm(ji,jj) ) ) 595 596 ! Calculate evaporation over water. (kg/m2/s) 597 !------------------------------------------------- 598 evap(ji,jj) = zqla_oce(ji,jj) / cevap 599 600 END DO !ji 601 END DO !jj 602 603 DO jl = 1, jpl 604 !CDIR NOVERRCHK 605 DO jj = 2 , jpjm1 606 !CDIR NOVERRCHK 607 DO ji = 1, jpi 608 609 ! Turbulent heat fluxes over snow/ice. 610 !-------------------------------------------------- 611 612 ! zesi : vapour pressure at saturation of ice 613 ! zqsati : humidity close to the ice surface (at saturation) 614 zesi = 611.0 * EXP ( 21.8745587 * tmask(ji,jj,1) & ! tmask needed to avoid overflow in the exponential 615 & * ( t_su(ji,jj,jl) - rtt )/ ( t_su(ji,jj,jl)- 7.66 ) ) 616 zqsati = ( 0.622 * zesi ) / ( zpatm(ji,jj) - 0.378 * zesi ) 617 618 ! computation of intermediate values 619 zticemb = t_su(ji,jj,jl) - 7.66 620 zticemb2 = zticemb * zticemb 621 ztice3 = t_su(ji,jj,jl) * t_su(ji,jj,jl) * t_su(ji,jj,jl) 622 zqsati2 = zqsati * zqsati 623 zesi2 = zesi * zesi 624 zdesidt = zesi * ( 9.5 * LOG( 10.0 ) * ( rtt - 7.66 ) / zticemb2 ) 625 626 ! Transfer cofficient zcshi and zclei over ice. Assumed to be constant Parkinson 1979 ; Maykut 1982 627 !-------------------------------------------------------------------- 628 zcshi = 1.75e-03 629 zclei = zcshi 630 631 ! Computation of sensible and latent fluxes over ice 632 !---------------------------------------------------------------- 633 634 ! computation of intermediate values 635 zrhova = zrhoa(ji,jj) * vatm(ji,jj) 636 zrhovaclei = zrhova * zcshi * 2.834e+06 637 zrhovacshi = zrhova * zclei * 1004.0 638 639 ! sensible heat flux 640 zqsb_ice(ji,jj,jl) = zrhovacshi * ( t_su(ji,jj,jl) - ztatm(ji,jj) ) 641 642 ! latent heat flux 643 zqla_ice(ji,jj,jl) = zrhovaclei * ( zqsati - zqatm(ji,jj) ) 644 qla_ice (ji,jj,jl) = MAX(0.e0, zqla_ice(ji,jj,jl) ) 645 646 ! Computation of sensitivity of non solar fluxes (dQ/dT) 647 !--------------------------------------------------------------- 648 649 ! computation of long-wave, sensible and latent flux sensitivity 650 zdqlw_ice = 4.0 * emic * stefan * ztice3 651 zdqsb_ice = zrhovacshi 652 zdqla_ice = zrhovaclei * ( zdesidt * ( zqsati2 / zesi2 ) * ( zpatm(ji,jj) / 0.622 ) ) 653 654 ! total non solar sensitivity 655 dqns_ice(ji,jj,jl) = -( zdqlw_ice + zdqsb_ice + zdqla_ice ) 656 657 ! latent flux sensitivity 658 dqla_ice(ji,jj,jl) = zdqla_ice 659 660 END DO 661 END DO 662 END DO !jl 663 664 ! total non solar heat flux over ice 665 qnsr_ice(:,:,:) = zqlw_ice(:,:,:) - zqsb_ice(:,:,:) - zqla_ice(:,:,:) 666 ! total non solar heat flux over water 667 qnsr_oce(:,:) = zqlw_oce(:,:) - zqsb_oce(:,:) - zqla_oce(:,:) 668 669 ! solid precipitations ( kg/m2/day -> kg/m2/s) 670 tprecip(:,:) = tprecip (:,:) / rday 671 ! snow precipitations ( kg/m2/day -> kg/m2/s) 672 sprecip(:,:) = sprecip (:,:) / rday 673 674 CALL lbc_lnk( qsr_oce (:,:) , 'T', 1. ) 675 CALL lbc_lnk( qnsr_oce(:,:) , 'T', 1. ) 676 CALL lbc_lnk( fr1_i0 (:,:) , 'T', 1. ) 677 CALL lbc_lnk( fr2_i0 (:,:) , 'T', 1. ) 678 CALL lbc_lnk( tprecip (:,:) , 'T', 1. ) 679 CALL lbc_lnk( sprecip (:,:) , 'T', 1. ) 680 CALL lbc_lnk( evap (:,:) , 'T', 1. ) 681 DO jl = 1, jpl 682 CALL lbc_lnk( qsr_ice (:,:,jl) , 'T', 1. ) 683 CALL lbc_lnk( qnsr_ice(:,:,jl) , 'T', 1. ) 684 CALL lbc_lnk( dqns_ice(:,:,jl) , 'T', 1. ) 685 CALL lbc_lnk( qla_ice (:,:,jl) , 'T', 1. ) 686 CALL lbc_lnk( dqla_ice(:,:,jl) , 'T', 1. ) 687 END DO 688 689 qsr_oce (:,:) = qsr_oce (:,:)*tmask(:,:,1) 690 qnsr_oce(:,:) = qnsr_oce(:,:)*tmask(:,:,1) 691 DO jl = 1, jpl 692 qsr_ice (:,:,jl) = qsr_ice (:,:,jl)*tmask(:,:,1) 693 qnsr_ice(:,:,jl) = qnsr_ice(:,:,jl)*tmask(:,:,1) 694 qla_ice (:,:,jl) = qla_ice (:,:,jl)*tmask(:,:,1) 695 dqns_ice(:,:,jl) = dqns_ice(:,:,jl)*tmask(:,:,1) 696 dqla_ice(:,:,jl) = dqla_ice(:,:,jl)*tmask(:,:,1) 697 END DO 698 fr1_i0 (:,:) = fr1_i0 (:,:)*tmask(:,:,1) 699 fr2_i0 (:,:) = fr2_i0 (:,:)*tmask(:,:,1) 700 tprecip (:,:) = tprecip (:,:)*tmask(:,:,1) 701 sprecip (:,:) = sprecip (:,:)*tmask(:,:,1) 702 evap (:,:) = evap (:,:)*tmask(:,:,1) 703 704 705 END SUBROUTINE flx_blk 706 707 708 #else 81 709 82 710 SUBROUTINE flx_blk( psst ) … … 215 843 ! Initilization ! 216 844 !--------------------- 217 #if ! defined key_ ice_lim845 #if ! defined key_lim2 218 846 tn_ice(:,:) = psst(:,:) 219 847 #endif … … 680 1308 681 1309 END SUBROUTINE flx_blk 1310 #endif 682 1311 683 1312
Note: See TracChangeset
for help on using the changeset viewer.
