[4728] | 1 | ! (C) Copyright 2019- ECMWF. |
---|
| 2 | ! |
---|
| 3 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
| 4 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 5 | ! |
---|
| 6 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 7 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
| 8 | ! nor does it submit to any jurisdiction. |
---|
| 9 | |
---|
| 10 | MODULE YOE_SPECTRAL_PLANCK |
---|
| 11 | |
---|
| 12 | ! YOE_SPECTRAL_PLANCK |
---|
| 13 | ! |
---|
| 14 | ! PURPOSE |
---|
| 15 | ! ------- |
---|
| 16 | ! Calculate Planck function integrated across user-specified |
---|
| 17 | ! spectral intervals, used in RADHEATN by approximate longwave |
---|
| 18 | ! update scheme to modify longwave fluxes to account for the |
---|
| 19 | ! spectral emissivity on the high-resolution model grid (rather than |
---|
| 20 | ! the lower resolution grid seen by the radiation scheme). |
---|
| 21 | ! |
---|
| 22 | ! INTERFACE |
---|
| 23 | ! --------- |
---|
| 24 | ! Call the INIT member routine to configure the look-up table of the |
---|
| 25 | ! TSPECRALPLANCK type, followed by any number of CALC calls with the |
---|
| 26 | ! temperatures at which the Planck function is required. FREE then |
---|
| 27 | ! deallocates memory. |
---|
| 28 | ! |
---|
| 29 | ! AUTHOR |
---|
| 30 | ! ------ |
---|
| 31 | ! Robin Hogan, ECMWF |
---|
| 32 | ! Original: 2019-02-04 |
---|
| 33 | ! |
---|
| 34 | ! MODIFICATIONS |
---|
| 35 | ! ------------- |
---|
| 36 | ! A Dawson 2019-08-05 avoid single precision overflow in INIT |
---|
| 37 | |
---|
| 38 | !----------------------------------------------------------------------- |
---|
| 39 | |
---|
| 40 | USE PARKIND1, ONLY : JPRB,JPRD,JPIM |
---|
| 41 | IMPLICIT NONE |
---|
| 42 | SAVE |
---|
| 43 | |
---|
| 44 | !----------------------------------------------------------------------- |
---|
| 45 | ! Type for storing Planck function look-up table |
---|
| 46 | TYPE TSPECTRALPLANCK |
---|
| 47 | ! Number of intervals over which the integrated Planck function is |
---|
| 48 | ! required. Note that an interval need not be contiguous in |
---|
| 49 | ! wavelength. |
---|
| 50 | INTEGER(KIND=JPIM) :: NINTERVALS |
---|
| 51 | |
---|
| 52 | ! Number of temperatures in look-up table |
---|
| 53 | INTEGER(KIND=JPIM) :: NTEMPS |
---|
| 54 | |
---|
| 55 | ! Start temperature and temperature spacing of look-up table |
---|
| 56 | REAL(KIND=JPRB) :: TEMP1, DTEMP |
---|
| 57 | |
---|
| 58 | ! Integrated Planck functions in look-up table, dimensioned |
---|
| 59 | ! (NINTERVALS,NTEMPS) |
---|
| 60 | REAL(KIND=JPRB), ALLOCATABLE :: PLANCK_LUT(:,:) |
---|
| 61 | |
---|
| 62 | ! Store interval data |
---|
| 63 | REAL(KIND=JPRB), ALLOCATABLE :: WAVLEN_BOUND(:) |
---|
| 64 | INTEGER(KIND=JPIM), ALLOCATABLE :: INTERVAL_MAP(:) |
---|
| 65 | |
---|
| 66 | CONTAINS |
---|
| 67 | PROCEDURE :: INIT |
---|
| 68 | PROCEDURE :: CALC |
---|
| 69 | PROCEDURE :: PRINT=> PRINT_SPECTRAL_PLANCK |
---|
| 70 | PROCEDURE :: FREE => FREE_SPECTRAL_PLANCK |
---|
| 71 | |
---|
| 72 | END TYPE TSPECTRALPLANCK |
---|
| 73 | |
---|
| 74 | CONTAINS |
---|
| 75 | |
---|
| 76 | !----------------------------------------------------------------------- |
---|
| 77 | ! Generate a Planck function look-up table consisting of KINTERVALS |
---|
| 78 | ! spectral intervals (which need not be contiguous in wavelength), |
---|
| 79 | ! whose wavelength bounds are defined by PWAVLEN_BOUND and mapping |
---|
| 80 | ! on to KINTERVALS described by KINTERVAL_MAP. |
---|
| 81 | SUBROUTINE INIT(SELF, KINTERVALS, PWAVLEN_BOUND, KINTERVAL_MAP) |
---|
| 82 | |
---|
| 83 | USE YOMCST, ONLY : RPI, RKBOL, RHPLA, RCLUM |
---|
| 84 | USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK |
---|
| 85 | USE YOMLUN, ONLY : NULOUT |
---|
| 86 | |
---|
| 87 | CLASS(TSPECTRALPLANCK), INTENT(INOUT) :: SELF |
---|
| 88 | INTEGER(KIND=JPIM) , INTENT(IN) :: KINTERVALS |
---|
| 89 | REAL(KIND=JPRB) , INTENT(IN) :: PWAVLEN_BOUND(:) |
---|
| 90 | INTEGER(KIND=JPIM) , INTENT(IN) :: KINTERVAL_MAP(:) |
---|
| 91 | |
---|
| 92 | ! Current temperature (K) |
---|
| 93 | REAL(KIND=JPRB) :: ZTEMP |
---|
| 94 | |
---|
| 95 | ! Combinations of constants in the Planck function |
---|
| 96 | REAL(KIND=JPRB) :: ZCOEFF1, ZCOEFF2 |
---|
| 97 | |
---|
| 98 | ! Wavelengths at start and end of range |
---|
| 99 | REAL(KIND=JPRB) :: ZWAVLEN1, ZWAVLEN2, DWAVLEN |
---|
| 100 | |
---|
| 101 | ! Wavelength, wavelength squared |
---|
| 102 | REAL(KIND=JPRB) :: ZWAVLEN, ZWAVLEN_SQR |
---|
| 103 | |
---|
| 104 | ! Sum of Planck values, integration weight |
---|
| 105 | REAL(KIND=JPRB) :: ZSUM, ZWEIGHT |
---|
| 106 | |
---|
| 107 | ! A double-precision temporary to hold the exponential term in Planck's law. |
---|
| 108 | ! A single precision float can overflow during this calculation. |
---|
| 109 | REAL(KIND=JPRD) :: ZPLANCKEXP |
---|
| 110 | |
---|
| 111 | ! Number of wavelength ranges represented by PWAVLEN_BOUND and |
---|
| 112 | ! KINTERVAL_MAP |
---|
| 113 | INTEGER(KIND=JPIM) :: NRANGES, NWAVLEN |
---|
| 114 | |
---|
| 115 | INTEGER(KIND=JPIM) :: JT, JI, JW, JR |
---|
| 116 | |
---|
| 117 | REAL(KIND=JPHOOK) :: ZHOOK_HANDLE |
---|
| 118 | |
---|
| 119 | #include "abor1.intfb.h" |
---|
| 120 | |
---|
| 121 | IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',0,ZHOOK_HANDLE) |
---|
| 122 | |
---|
| 123 | IF (KINTERVALS == 1) THEN |
---|
| 124 | ! We can use Stefan-Boltzmann law without a look-up table |
---|
| 125 | WRITE(NULOUT,'(a)') 'YOE_SPECTRAL_PLANCK: Single-band look-up table requested: use Stefan-Boltzmann law' |
---|
| 126 | SELF%NINTERVALS = KINTERVALS |
---|
| 127 | CALL SELF%FREE |
---|
| 128 | ELSE |
---|
| 129 | ! Full look-up table required |
---|
| 130 | ZCOEFF1 = 2.0_JPRB * RPI * RHPLA * RCLUM * RCLUM |
---|
| 131 | ZCOEFF2 = RHPLA * RCLUM / RKBOL |
---|
| 132 | |
---|
| 133 | NRANGES = SIZE(KINTERVAL_MAP,1) |
---|
| 134 | IF (SIZE(PWAVLEN_BOUND,1) /= NRANGES-1) THEN |
---|
| 135 | CALL ABOR1('YOS_SPECTRAL_PLANCK:INIT: PWAVLEN_BOUND must have one fewer elements than KINTERVAL_MAP') |
---|
| 136 | ENDIF |
---|
| 137 | |
---|
| 138 | CALL SELF%FREE |
---|
| 139 | |
---|
| 140 | ALLOCATE(SELF%WAVLEN_BOUND(NRANGES-1)) |
---|
| 141 | ALLOCATE(SELF%INTERVAL_MAP(NRANGES)) |
---|
| 142 | SELF%WAVLEN_BOUND(1:NRANGES-1) = PWAVLEN_BOUND(1:NRANGES-1) |
---|
| 143 | SELF%INTERVAL_MAP(1:NRANGES) = KINTERVAL_MAP(1:NRANGES) |
---|
| 144 | |
---|
| 145 | SELF%NINTERVALS = KINTERVALS |
---|
| 146 | ! Temperature in 1-K intervals from 150 K to 350 K |
---|
| 147 | SELF%TEMP1 = 150.0_JPRB |
---|
| 148 | SELF%DTEMP = 1.0_JPRB |
---|
| 149 | SELF%NTEMPS = 1 + NINT((350.0_JPRB - SELF%TEMP1) / SELF%DTEMP) |
---|
| 150 | |
---|
| 151 | ALLOCATE(SELF%PLANCK_LUT(SELF%NINTERVALS,SELF%NTEMPS)) |
---|
| 152 | SELF%PLANCK_LUT(:,:) = 0.0_JPRB |
---|
| 153 | |
---|
| 154 | ! Print the properties of the look-up table |
---|
| 155 | WRITE(NULOUT,'(a,i0,a,f5.1,a,f5.1,a)') & |
---|
| 156 | & 'YOE_SPECTRAL_PLANCK: Generating Planck look-up table with ', & |
---|
| 157 | & SELF%NTEMPS, ' temperatures from ', & |
---|
| 158 | & SELF%TEMP1, ' to ', SELF%TEMP1+SELF%DTEMP*(SELF%NTEMPS-1), ' K:' |
---|
| 159 | DO JI = 1,SELF%NINTERVALS |
---|
| 160 | WRITE(NULOUT,'(a,i0,a)',advance='no') ' Band ', JI, ':' |
---|
| 161 | DO JR = 1,NRANGES |
---|
| 162 | IF (KINTERVAL_MAP(JR) == JI) THEN |
---|
| 163 | IF (JR == 1) THEN |
---|
| 164 | WRITE(NULOUT,'(a,f0.2)',advance='no') ' 0.00-', & |
---|
| 165 | & PWAVLEN_BOUND(1)*1.0e6_JPRB |
---|
| 166 | ELSEIF (JR == NRANGES) THEN |
---|
| 167 | WRITE(NULOUT,'(a,f0.2,a)',advance='no') ' ', & |
---|
| 168 | & PWAVLEN_BOUND(JR-1)*1.0e6_JPRB, '-Inf' |
---|
| 169 | ELSE |
---|
| 170 | WRITE(NULOUT,'(a,f0.2,a,f0.2)',advance='no') ' ', & |
---|
| 171 | & PWAVLEN_BOUND(JR-1)*1.0e6_JPRB, '-', & |
---|
| 172 | & PWAVLEN_BOUND(JR)*1.0e6_JPRB |
---|
| 173 | ENDIF |
---|
| 174 | ENDIF |
---|
| 175 | ENDDO |
---|
| 176 | WRITE(NULOUT,'(a)') ' microns' |
---|
| 177 | ENDDO |
---|
| 178 | |
---|
| 179 | ! Create the look-up table |
---|
| 180 | DO JT = 1,SELF%NTEMPS |
---|
| 181 | |
---|
| 182 | ZTEMP = SELF%TEMP1 + (JT-1) * SELF%DTEMP |
---|
| 183 | |
---|
| 184 | DO JI = 1,NRANGES |
---|
| 185 | |
---|
| 186 | IF (JI == 1) THEN |
---|
| 187 | ZWAVLEN1 = MIN(1.0E-6_JPRB, 0.8_JPRB * PWAVLEN_BOUND(1)) |
---|
| 188 | ZWAVLEN2 = PWAVLEN_BOUND(1) |
---|
| 189 | ELSEIF (JI == NRANGES) THEN |
---|
| 190 | ZWAVLEN1 = PWAVLEN_BOUND(NRANGES-1) |
---|
| 191 | ! Simulate up to at least 200 microns wavelength |
---|
| 192 | ZWAVLEN2 = MAX(200.0E-6_JPRB, PWAVLEN_BOUND(NRANGES-1)+20.0E-6_JPRB) |
---|
| 193 | ELSE |
---|
| 194 | ZWAVLEN1 = PWAVLEN_BOUND(JI-1) |
---|
| 195 | ZWAVLEN2 = PWAVLEN_BOUND(JI) |
---|
| 196 | ENDIF |
---|
| 197 | |
---|
| 198 | NWAVLEN = 100 |
---|
| 199 | DWAVLEN = (ZWAVLEN2 - ZWAVLEN1) / NWAVLEN |
---|
| 200 | ZSUM = 0.0_JPRB |
---|
| 201 | DO JW = 0,NWAVLEN |
---|
| 202 | ZWAVLEN = ZWAVLEN1 + DWAVLEN*JW |
---|
| 203 | ! Weights for trapezoidal rule |
---|
| 204 | !IF (JW > 0 .AND. JW < NWAVLEN) THEN |
---|
| 205 | ! ZWEIGHT = 2.0_JPRB |
---|
| 206 | !ELSE |
---|
| 207 | ! ZWEIGHT = 1.0_JPRB |
---|
| 208 | !ENDIF |
---|
| 209 | ! Weights for Simpson's rule |
---|
| 210 | IF (JW > 0 .AND. JW < NWAVLEN) THEN |
---|
| 211 | ZWEIGHT = 2.0_JPRB + 2.0_JPRB * MOD(JW,2) |
---|
| 212 | ELSE |
---|
| 213 | ZWEIGHT = 1.0_JPRB |
---|
| 214 | ENDIF |
---|
| 215 | ! Planck's law |
---|
| 216 | ! |
---|
| 217 | ! The exponential term is computed in double precision to avoid |
---|
| 218 | ! overflow. The final result should still be in the range of a single |
---|
| 219 | ! precision float. |
---|
| 220 | ZWAVLEN_SQR = ZWAVLEN*ZWAVLEN |
---|
| 221 | ZPLANCKEXP = EXP(REAL(ZCOEFF2, JPRD) & |
---|
| 222 | & / (REAL(ZWAVLEN, JPRD) * REAL(ZTEMP, JPRD))) |
---|
| 223 | ZSUM = ZSUM + ZWEIGHT / (ZWAVLEN_SQR*ZWAVLEN_SQR*ZWAVLEN & |
---|
| 224 | & * (ZPLANCKEXP - 1.0_JPRB)) |
---|
| 225 | ENDDO |
---|
| 226 | SELF%PLANCK_LUT(KINTERVAL_MAP(JI),JT) = SELF%PLANCK_LUT(KINTERVAL_MAP(JI),JT) & |
---|
| 227 | & + ZCOEFF1 * ZSUM * DWAVLEN / 3.0_JPRB |
---|
| 228 | ENDDO |
---|
| 229 | |
---|
| 230 | ENDDO |
---|
| 231 | |
---|
| 232 | ENDIF |
---|
| 233 | |
---|
| 234 | IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',1,ZHOOK_HANDLE) |
---|
| 235 | |
---|
| 236 | END SUBROUTINE INIT |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | !----------------------------------------------------------------------- |
---|
| 240 | ! Calculate Planck function in spectral intervals from temperature |
---|
| 241 | SUBROUTINE CALC(SELF, KIDIA, KFDIA, KLON, PTEMPERATURE, PPLANCK) |
---|
| 242 | |
---|
| 243 | USE YOMCST, ONLY : RSIGMA |
---|
| 244 | USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK |
---|
| 245 | |
---|
| 246 | CLASS(TSPECTRALPLANCK), INTENT(IN) :: SELF |
---|
| 247 | ! Process columns KIDIA-KFDIA from total of KLON columns |
---|
| 248 | INTEGER(KIND=JPIM) , INTENT(IN) :: KIDIA, KFDIA, KLON |
---|
| 249 | ! Temperature in Kelvin |
---|
| 250 | REAL(KIND=JPRB) , INTENT(IN) :: PTEMPERATURE(KLON) |
---|
| 251 | ! Integrated Planck function as an irradiance, in W m-2 |
---|
| 252 | REAL(KIND=JPRB) , INTENT(OUT) :: PPLANCK(KLON,SELF%NINTERVALS) |
---|
| 253 | |
---|
| 254 | ! Column loop counter, index to temperature interval |
---|
| 255 | INTEGER(KIND=JPRB) :: JL, ITEMP |
---|
| 256 | |
---|
| 257 | ! Interpolation weight, highest temperature in look-up table |
---|
| 258 | REAL(KIND=JPRB) :: ZWEIGHT, ZTEMP2 |
---|
| 259 | |
---|
| 260 | REAL(KIND=JPHOOK) :: ZHOOK_HANDLE |
---|
| 261 | |
---|
| 262 | IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',0,ZHOOK_HANDLE) |
---|
| 263 | |
---|
| 264 | IF (SELF%NINTERVALS == 1) THEN |
---|
| 265 | ! Stefan-Boltzmann law |
---|
| 266 | PPLANCK(KIDIA:KFDIA,1) = RSIGMA * PTEMPERATURE(KIDIA:KFDIA)**4 |
---|
| 267 | |
---|
| 268 | ELSE |
---|
| 269 | ! Look-up table |
---|
| 270 | ZTEMP2 = SELF%TEMP1 + SELF%DTEMP * (SELF%NTEMPS - 1) |
---|
| 271 | |
---|
| 272 | DO JL = KIDIA,KFDIA |
---|
| 273 | |
---|
| 274 | IF (PTEMPERATURE(JL) <= SELF%TEMP1) THEN |
---|
| 275 | ! Cap the Planck function at the low end |
---|
| 276 | ITEMP = 1 |
---|
| 277 | ZWEIGHT = 0.0_JPRB |
---|
| 278 | ELSEIF (PTEMPERATURE(JL) < ZTEMP2) THEN |
---|
| 279 | ! Linear interpolation |
---|
| 280 | ZWEIGHT = 1.0_JPRB + (PTEMPERATURE(JL) - SELF%TEMP1) / SELF%DTEMP |
---|
| 281 | ITEMP = NINT(ZWEIGHT) |
---|
| 282 | ZWEIGHT = ZWEIGHT - ITEMP |
---|
| 283 | ELSE |
---|
| 284 | ! Linear extrapolation at high temperatures off the scale |
---|
| 285 | ITEMP = SELF%NTEMPS-1 |
---|
| 286 | ZWEIGHT = 1.0_JPRB + (PTEMPERATURE(JL) - SELF%TEMP1) / SELF%DTEMP - ITEMP |
---|
| 287 | ENDIF |
---|
| 288 | |
---|
| 289 | PPLANCK(JL,:) = SELF%PLANCK_LUT(:,ITEMP) & |
---|
| 290 | & + ZWEIGHT * (SELF%PLANCK_LUT(:,ITEMP+1) - SELF%PLANCK_LUT(:,ITEMP)) |
---|
| 291 | |
---|
| 292 | ! Force sum to equal Stefan-Boltzmann law |
---|
| 293 | PPLANCK(JL,:) = PPLANCK(JL,:) * RSIGMA * PTEMPERATURE(JL)**4 / SUM(PPLANCK(JL,:),1) |
---|
| 294 | |
---|
| 295 | ENDDO |
---|
| 296 | |
---|
| 297 | ENDIF |
---|
| 298 | |
---|
| 299 | IF (LHOOK) CALL DR_HOOK('YOE_SPECTRAL_PLANCK:INIT',1,ZHOOK_HANDLE) |
---|
| 300 | |
---|
| 301 | END SUBROUTINE CALC |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | !----------------------------------------------------------------------- |
---|
| 305 | ! Print look-up table to a unit |
---|
| 306 | SUBROUTINE PRINT_SPECTRAL_PLANCK(SELF, IUNIT) |
---|
| 307 | |
---|
| 308 | CLASS(TSPECTRALPLANCK), INTENT(IN) :: SELF |
---|
| 309 | INTEGER(KIND=JPIM), INTENT(IN) :: IUNIT |
---|
| 310 | |
---|
| 311 | INTEGER(KIND=JPIM) :: JT |
---|
| 312 | |
---|
| 313 | CHARACTER(len=24) :: MY_FORMAT |
---|
| 314 | |
---|
| 315 | IF (SELF%NINTERVALS == 1) THEN |
---|
| 316 | |
---|
| 317 | WRITE(IUNIT,'(A)') 'Spectral Planck in only one interval: using Stefan-Boltzmann law' |
---|
| 318 | |
---|
| 319 | ELSE |
---|
| 320 | |
---|
| 321 | WRITE(IUNIT,'(A,I0,A)') 'Spectral Planck look-up table defined in ', & |
---|
| 322 | & SELF%NINTERVALS, ' intervals:' |
---|
| 323 | WRITE(MY_FORMAT,'(A,I0,A)') '(f7.2,', SELF%NINTERVALS, 'e15.5)' |
---|
| 324 | DO JT = 1,SELF%NTEMPS |
---|
| 325 | WRITE(IUNIT,TRIM(MY_FORMAT)) SELF%TEMP1 + (JT-1) * SELF%DTEMP, & |
---|
| 326 | & SELF%PLANCK_LUT(:,JT) |
---|
| 327 | ENDDO |
---|
| 328 | |
---|
| 329 | ENDIF |
---|
| 330 | |
---|
| 331 | END SUBROUTINE PRINT_SPECTRAL_PLANCK |
---|
| 332 | |
---|
| 333 | |
---|
| 334 | !----------------------------------------------------------------------- |
---|
| 335 | ! Free allocated memory |
---|
| 336 | SUBROUTINE FREE_SPECTRAL_PLANCK(SELF) |
---|
| 337 | |
---|
| 338 | CLASS(TSPECTRALPLANCK), INTENT(INOUT) :: SELF |
---|
| 339 | |
---|
| 340 | IF (ALLOCATED(SELF%PLANCK_LUT)) THEN |
---|
| 341 | DEALLOCATE(SELF%PLANCK_LUT) |
---|
| 342 | ENDIF |
---|
| 343 | IF (ALLOCATED(SELF%WAVLEN_BOUND)) THEN |
---|
| 344 | DEALLOCATE(SELF%WAVLEN_BOUND) |
---|
| 345 | ENDIF |
---|
| 346 | IF (ALLOCATED(SELF%INTERVAL_MAP)) THEN |
---|
| 347 | DEALLOCATE(SELF%INTERVAL_MAP) |
---|
| 348 | ENDIF |
---|
| 349 | |
---|
| 350 | END SUBROUTINE FREE_SPECTRAL_PLANCK |
---|
| 351 | |
---|
| 352 | |
---|
| 353 | END MODULE YOE_SPECTRAL_PLANCK |
---|