Changeset 4489 for LMDZ6/trunk/libf/phylmd/ecrad/radiation_two_stream.F90
- Timestamp:
- Mar 31, 2023, 8:42:57 PM (17 months ago)
- Location:
- LMDZ6/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk
- Property svn:mergeinfo changed
/LMDZ6/branches/LMDZ_ECRad (added) merged: 4175,4177-4183,4188,4192,4200-4203,4355,4366,4387-4388,4390,4444,4482,4486,4488
- Property svn:mergeinfo changed
-
LMDZ6/trunk/libf/phylmd/ecrad/radiation_two_stream.F90
r3908 r4489 18 18 ! 2017-07-26 R Hogan Added calc_frac_scattered_diffuse_sw routine 19 19 ! 2017-10-23 R Hogan Renamed single-character variables 20 ! 2021-02-19 R Hogan Security for shortwave singularity 20 21 21 22 module radiation_two_stream … … 31 32 ! think of acos(1/lw_diffusivity) to be the effective zenith angle 32 33 ! of longwave radiation. 33 real(jprd), parameter :: LwDiffusivity = 1.66_jprd 34 real(jprd), parameter :: LwDiffusivity = 1.66_jprd 35 real(jprb), parameter :: LwDiffusivityWP = 1.66_jprb ! Working precision version 34 36 35 37 ! Shortwave diffusivity factor assumes hemispheric isotropy, assumed … … 87 89 if (lhook) call dr_hook('radiation_two_stream:calc_two_stream_gammas_lw',0,hook_handle) 88 90 #endif 89 91 ! Added for DWD (2020) 92 !NEC$ shortloop 90 93 do jg = 1, ng 91 94 ! Fu et al. (1997), Eq 2.9 and 2.10: … … 136 139 ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to 137 140 ! Atmospheric Physics 53, 147-66) 141 ! Added for DWD (2020) 142 !NEC$ shortloop 138 143 do jg = 1, ng 139 144 ! gamma1(jg) = 2.0_jprb - ssa(jg) * (1.25_jprb + 0.75_jprb*g(jg)) … … 205 210 #endif 206 211 212 ! Added for DWD (2020) 213 !NEC$ shortloop 207 214 do jg = 1, ng 208 215 if (od(jg) > 1.0e-3_jprd) then … … 293 300 #endif 294 301 302 ! Added for DWD (2020) 303 !NEC$ shortloop 295 304 do jg = 1, ng 296 305 k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & … … 359 368 #endif 360 369 370 ! Added for DWD (2020) 371 !NEC$ shortloop 361 372 do jg = 1, ng 362 373 ! Compute upward and downward emission assuming the Planck … … 450 461 integer :: jg 451 462 463 ! Local value of cosine of solar zenith angle, in case it needs to be 464 ! tweaked to avoid near division by zero. This is intentionally in working 465 ! precision (jprb) rather than fixing at double precision (jprd). 466 real(jprb) :: mu0_local 467 452 468 #ifdef DO_DR_HOOK_TWO_STREAM 453 469 real(jprb) :: hook_handle … … 456 472 #endif 457 473 474 ! Added for DWD (2020) 475 !NEC$ shortloop 458 476 do jg = 1, ng 459 od_over_mu0 = max(od(jg) / mu0, 0.0_jprd) 460 ! In the IFS this appears to be faster without testing the value 461 ! of od_over_mu0: 462 if (.true.) then 463 ! if (od_over_mu0 > 1.0e-6_jprd) then 477 464 478 gamma4 = 1.0_jprd - gamma3(jg) 465 479 alpha1 = gamma1(jg)*gamma4 + gamma2(jg)*gamma3(jg) ! Eq. 16 466 480 alpha2 = gamma1(jg)*gamma3(jg) + gamma2(jg)*gamma4 ! Eq. 17 467 481 482 k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & 483 & 1.0e-12_jprd)) ! Eq 18 484 485 ! We had a rare crash where k*mu0 was within around 1e-13 of 1, 486 ! leading to ref_dir and trans_dir_diff being well outside the range 487 ! 0-1. The following approach is appropriate when k_exponent is double 488 ! precision and mu0_local is single precision, although work is needed 489 ! to make this entire routine secure in single precision. 490 mu0_local = mu0 491 if (abs(1.0_jprd - k_exponent*mu0) < 1000.0_jprd * epsilon(1.0_jprd)) then 492 mu0_local = mu0 * (1.0_jprb - 10.0_jprb*epsilon(1.0_jprb)) 493 end if 494 495 od_over_mu0 = max(od(jg) / mu0_local, 0.0_jprd) 496 468 497 ! Note that if the minimum value is reduced (e.g. to 1.0e-24) 469 498 ! then noise starts to appear as a function of solar zenith 470 499 ! angle 471 k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & 472 & 1.0e-12_jprd)) ! Eq 18 473 k_mu0 = k_exponent*mu0 500 k_mu0 = k_exponent*mu0_local 474 501 k_gamma3 = k_exponent*gamma3(jg) 475 502 k_gamma4 = k_exponent*gamma4 … … 482 509 k_2_exponential = 2.0_jprd * k_exponent * exponential 483 510 484 if (k_mu0 == 1.0_jprd) then485 k_mu0 = 1.0_jprd - 10.0_jprd*epsilon(1.0_jprd)486 end if487 488 511 reftrans_factor = 1.0_jprd / (k_exponent + gamma1(jg) + (k_exponent - gamma1(jg))*exponential2) 489 512 … … 498 521 ! to be the flux into a plane perpendicular to the direction of 499 522 ! the sun, not into a horizontal plane 500 reftrans_factor = mu0 * ssa(jg) * reftrans_factor / (1.0_jprd - k_mu0*k_mu0)523 reftrans_factor = mu0_local * ssa(jg) * reftrans_factor / (1.0_jprd - k_mu0*k_mu0) 501 524 502 525 ! Meador & Weaver (1980) Eq. 14, multiplying top & bottom by … … 505 528 & * ( (1.0_jprd - k_mu0) * (alpha2 + k_gamma3) & 506 529 & -(1.0_jprd + k_mu0) * (alpha2 - k_gamma3)*exponential2 & 507 & -k_2_exponential*(gamma3(jg) - alpha2*mu0 )*exponential0)530 & -k_2_exponential*(gamma3(jg) - alpha2*mu0_local)*exponential0) 508 531 509 532 ! Meador & Weaver (1980) Eq. 15, multiplying top & bottom by 510 533 ! exp(-k_exponent*od), minus the 1*exp(-od/mu0) term representing direct 511 534 ! unscattered transmittance. 512 trans_dir_diff(jg) = reftrans_factor * ( k_2_exponential*(gamma4 + alpha1*mu0 ) &535 trans_dir_diff(jg) = reftrans_factor * ( k_2_exponential*(gamma4 + alpha1*mu0_local) & 513 536 & - exponential0 & 514 537 & * ( (1.0_jprd + k_mu0) * (alpha1 + k_gamma4) & 515 538 & -(1.0_jprd - k_mu0) * (alpha1 - k_gamma4) * exponential2) ) 516 539 517 else 518 ! Low optical-depth limit; see equations 19, 20 and 27 from 519 ! Meador & Weaver (1980) 520 trans_diff(jg) = 1.0_jprb - gamma1(jg) * od(jg) 521 ref_diff(jg) = gamma2(jg) * od(jg) 522 trans_dir_diff(jg) = (1.0_jprb - gamma3(jg)) * ssa(jg) * od(jg) 523 ref_dir(jg) = gamma3(jg) * ssa(jg) * od(jg) 524 trans_dir_dir(jg) = 1.0_jprd - od_over_mu0 525 end if 540 ! Final check that ref_dir + trans_dir_diff <= 1 541 ref_dir(jg) = max(0.0_jprb, min(ref_dir(jg), 1.0_jprb)) 542 trans_dir_diff(jg) = max(0.0_jprb, min(trans_dir_diff(jg), 1.0_jprb-ref_dir(jg))) 543 526 544 end do 527 545 … … 587 605 #endif 588 606 607 ! Added for DWD (2020) 608 !NEC$ shortloop 589 609 do jg = 1, ng 590 610 od_over_mu0 = max(gamma0(jg) * depth, 0.0_jprd) … … 699 719 #endif 700 720 721 ! Added for DWD (2020) 722 !NEC$ shortloop 701 723 do jg = 1, ng 702 724 ! Note that if the minimum value is reduced (e.g. to 1.0e-24)
Note: See TracChangeset
for help on using the changeset viewer.