1 2! KGEN-generated Fortran source file 3! 4! Filename : mo_lrtm_solver.f90 5! Generated at: 2015-02-19 15:30:36 6! KGEN version: 0.4.4 7 8 9 10 MODULE mo_lrtm_solver 11 USE mo_kind, ONLY: wp 12 USE mo_math_constants, ONLY: pi 13 USE mo_rrtm_params, ONLY: nbndlw 14 USE mo_rad_fastmath, ONLY: tautrans 15 USE mo_rad_fastmath, ONLY: transmit 16 IMPLICIT NONE 17 REAL(KIND=wp), parameter :: fluxfac = 2.0e+04_wp * pi 18 CONTAINS 19 20 ! read subroutines 21 ! ------------------------------------------------------------------------------- 22 23 SUBROUTINE lrtm_solver(kproma, kbdim, klev, tau, layplnk, levplnk, weights, secdiff, surfplanck, surfemis, fluxup, fluxdn) 24 ! 25 ! Compute IR (no scattering) radiative transfer for a set of columns 26 ! Based on AER code RRTMG_LW_RTNMC, including approximations used there 27 ! Layers are ordered from botton to top (i.e. tau(1) is tau in lowest layer) 28 ! Computes all-sky RT given a total optical thickness in each layer 29 ! 30 INTEGER, intent(in) :: kbdim 31 INTEGER, intent(in) :: klev 32 INTEGER, intent(in) :: kproma 33 !< Number of columns 34 !< Maximum number of columns as declared in calling (sub)program 35 !< number of layers (one fewer than levels) 36 REAL(KIND=wp), intent(in) :: tau(kbdim,klev) 37 REAL(KIND=wp), intent(in) :: layplnk(kbdim,klev) 38 REAL(KIND=wp), intent(in) :: weights(kbdim,klev) !< dimension (kbdim, klev) 39 !< Longwave optical thickness 40 !< Planck function at layer centers 41 !< Fraction of total Planck function for this g-point 42 REAL(KIND=wp), intent(in) :: levplnk(kbdim, 0:klev) 43 !< Planck function at layer edges, level i is the top of layer i 44 REAL(KIND=wp), intent(in) :: secdiff(kbdim) 45 REAL(KIND=wp), intent(in) :: surfemis(kbdim) 46 REAL(KIND=wp), intent(in) :: surfplanck(kbdim) !< dimension (kbdim) 47 !< Planck function at surface 48 !< Surface emissivity 49 !< secant of integration angle - depends on band, column water vapor 50 REAL(KIND=wp), intent(out) :: fluxup(kbdim, 0:klev) 51 REAL(KIND=wp), intent(out) :: fluxdn(kbdim, 0:klev) !< dimension (kbdim, 0:klev) 52 !< Fluxes at the interfaces 53 ! ----------- 54 INTEGER :: jk 55 !< Loop index for layers 56 REAL(KIND=wp) :: odepth(kbdim,klev) 57 REAL(KIND=wp) :: tfn(kbdim) 58 REAL(KIND=wp) :: dplnkup(kbdim,klev) 59 REAL(KIND=wp) :: dplnkdn(kbdim,klev) 60 REAL(KIND=wp) :: bbup(kbdim,klev) 61 REAL(KIND=wp) :: bbdn(kbdim,klev) 62 REAL(KIND=wp) :: trans(kbdim,klev) 63 !< Layer transmissivity 64 !< TFN_TBL 65 !< Tau transition function; i.e. the transition of the Planck 66 !< function from that for the mean layer temperature to that for 67 !< the layer boundary temperature as a function of optical depth. 68 !< The "linear in tau" method is used to make the table. 69 !< Upward derivative of Planck function 70 !< Downward derivative of Planck function 71 !< Interpolated downward emission 72 !< Interpolated upward emission 73 !< Effective IR optical depth of layer 74 REAL(KIND=wp) :: rad_dn(kbdim,0:klev) 75 REAL(KIND=wp) :: rad_up(kbdim,0:klev) 76 !< Radiance down at propagation angle 77 !< Radiance down at propagation angle 78 ! This secant and weight corresponds to the standard diffusivity 79 ! angle. The angle is redefined for some bands. 80 REAL(KIND=wp), parameter :: wtdiff = 0.5_wp 81 ! ----------- 82 ! 83 ! 1.0 Initial preparations 84 ! Weight optical depth by 1/cos(diffusivity angle), which depends on band 85 ! This will be used to compute layer transmittance 86 ! ----- 87 !IBM* ASSERT(NODEPS) 88 DO jk = 1, klev 89 odepth(1:kproma,jk) = max(0._wp, secdiff(1:kproma) * tau(1:kproma,jk)) 90 END DO 91 ! 92 ! 2.0 Radiative transfer 93 ! 94 ! ----- 95 ! 96 ! Plank function derivatives and total emission for linear-in-tau approximation 97 ! 98 !IBM* ASSERT(NODEPS) 99 DO jk = 1, klev 100 tfn(1:kproma) = tautrans(odepth(:,jk), kproma) 101 dplnkup(1:kproma,jk) = levplnk(1:kproma,jk) - layplnk(1:kproma,jk) 102 dplnkdn(1:kproma,jk) = levplnk(1:kproma,jk-1) - layplnk(1:kproma,jk) 103 bbup(1:kproma,jk) = weights(1:kproma,jk) * (layplnk(1:kproma,jk) + dplnkup(1:kproma,jk) * tfn(1:kproma)) 104 bbdn(1:kproma,jk) = weights(1:kproma,jk) * (layplnk(1:kproma,jk) + dplnkdn(1:kproma,jk) * tfn(1:kproma)) 105 END DO 106 ! ----- 107 ! 2.1 Downward radiative transfer 108 ! 109 ! Level 0 is closest to the ground 110 ! 111 rad_dn(:, klev) = 0. ! Upper boundary condition - no downwelling IR 112 DO jk = klev, 1, -1 113 trans(1:kproma,jk) = transmit(odepth(:,jk), kproma) 114 ! RHS is a rearrangment of rad_dn(:,jk) * (1._wp - trans(:,jk)) + trans(:,jk) * bbdn(:) 115 rad_dn(1:kproma,jk-1) = rad_dn(1:kproma,jk) + (bbdn(1:kproma,jk) - rad_dn(1:kproma,jk)) * trans(1:kproma,jk) 116 END DO 117 ! 118 ! 2.2 Surface contribution, including reflection 119 ! 120 rad_up(1:kproma, 0) = weights(1:kproma, 1) * surfemis(1:kproma) * surfplanck(1:kproma) + (1._wp - & 121 surfemis(1:kproma)) * rad_dn(1:kproma, 0) 122 ! 123 ! 2.3 Upward radiative transfer 124 ! 125 DO jk = 1, klev 126 rad_up(1:kproma,jk) = rad_up(1:kproma,jk-1) * (1._wp - trans(1:kproma,jk)) + trans(1:kproma,jk) * bbup(1:kproma,& 127 jk) 128 END DO 129 ! 130 ! 3.0 Covert intensities at diffusivity angles to fluxes 131 ! 132 ! ----- 133 fluxup(1:kproma, 0:klev) = rad_up(1:kproma,:) * wtdiff * fluxfac 134 fluxdn(1:kproma, 0:klev) = rad_dn(1:kproma,:) * wtdiff * fluxfac 135 END SUBROUTINE lrtm_solver 136 ! ------------------------------------------------------------------------------- 137 138 elemental FUNCTION find_secdiff(iband, pwvcm) 139 INTEGER, intent(in) :: iband 140 !< RRTMG LW band number 141 REAL(KIND=wp), intent(in) :: pwvcm 142 !< Precipitable water vapor (cm) 143 REAL(KIND=wp) :: find_secdiff 144 ! Compute diffusivity angle for Bands 2-3 and 5-9 to vary (between 1.50 145 ! and 1.80) as a function of total column water vapor. The function 146 ! has been defined to minimize flux and cooling rate errors in these bands 147 ! over a wide range of precipitable water values. 148 REAL(KIND=wp), dimension(nbndlw), parameter :: a0 = (/ 1.66_wp, 1.55_wp, 1.58_wp, 1.66_wp, 1.54_wp, 1.454_wp, & 149 1.89_wp, 1.33_wp, 1.668_wp, 1.66_wp, 1.66_wp, 1.66_wp, 1.66_wp, 1.66_wp, 1.66_wp, 1.66_wp /) 150 REAL(KIND=wp), dimension(nbndlw), parameter :: a1 = (/ 0.00_wp, 0.25_wp, 0.22_wp, 0.00_wp, 0.13_wp, 0.446_wp, & 151 -0.10_wp, 0.40_wp, -0.006_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp /) 152 REAL(KIND=wp), dimension(nbndlw), parameter :: a2 = (/ 0.00_wp, -12.0_wp, -11.7_wp, 0.00_wp, -0.72_wp,-0.243_wp, & 153 0.19_wp,-0.062_wp, 0.414_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp /) 154 IF (iband == 1 .or. iband == 4 .or. iband >= 10) THEN 155 find_secdiff = 1.66_wp 156 ELSE 157 find_secdiff = max(min(a0(iband) + a1(iband) * exp(a2(iband)*pwvcm), 1.8_wp), 1.5_wp) 158 END IF 159 END FUNCTION find_secdiff 160 ! ------------------------------------------------------------------------------- 161 END MODULE mo_lrtm_solver 162