source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/yoe_spectral_planck.F90

Last change on this file was 5159, checked in by abarral, 3 months ago

Put dimensions.h and paramet.h into modules
File size: 11.2 KB
Line 
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
10MODULE 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
40USE PARKIND1,         ONLY :   JPRB,JPRD,JPIM
41IMPLICIT NONE
42SAVE
43
44!-----------------------------------------------------------------------
45! Type for storing Planck function look-up table
46TYPE 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
66CONTAINS
67  PROCEDURE :: INIT
68  PROCEDURE :: CALC
69  PROCEDURE :: PRINT=> PRINT_SPECTRAL_PLANCK
70  PROCEDURE :: FREE => FREE_SPECTRAL_PLANCK
71
72END TYPE TSPECTRALPLANCK
73
74CONTAINS
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.
81SUBROUTINE 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
236END SUBROUTINE INIT
237
238
239!-----------------------------------------------------------------------
240! Calculate Planck function in spectral intervals from temperature
241SUBROUTINE 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
301END SUBROUTINE CALC
302
303
304!-----------------------------------------------------------------------
305! Print look-up table to a unit
306SUBROUTINE 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
331END SUBROUTINE PRINT_SPECTRAL_PLANCK
332
333
334!-----------------------------------------------------------------------
335! Free allocated memory
336SUBROUTINE 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
350END SUBROUTINE FREE_SPECTRAL_PLANCK
351
352
353END MODULE YOE_SPECTRAL_PLANCK
Note: See TracBrowser for help on using the repository browser.