1 !-----------------------------------------------------------------------
3 ! Contains all subroutines to get time-varying mixing ratios of CO2,
4 ! N2O, CH4, CFC11 and CFC12 into radiation schemes.
5 ! These subroutines enable the user to specify the mixing ratios
6 ! in the file CAMtr_volume_mixing_ratio, giving the user an easy way
7 ! to specify trace gases concentrations.
9 ! The subroutines were first developed by and put together in this module
11 !-----------------------------------------------------------------------
12 MODULE module_ra_clWRF_support
18 INTEGER, PARAMETER :: r8 = 8
19 integer, parameter :: cyr = 1800 ! maximum num. of lines in 'CAMtr_volume_mixing_ratio' file
20 integer, DIMENSION(1:cyr), SAVE :: yrdata
21 real, DIMENSION(1:cyr), SAVE :: juldata
22 real(r8), DIMENSION(1:cyr), SAVE :: co2r, n2or, ch4r, cfc11r, cfc12r
25 ! Same values as in module_ra_cam_support.F
26 integer, parameter :: VARCAM_in_years = 233
27 integer :: yr_CAM(1:VARCAM_in_years) = &
28 (/ 1869, 1870, 1871, 1872, 1873, 1874, 1875, &
29 1876, 1877, 1878, 1879, 1880, 1881, 1882, &
30 1883, 1884, 1885, 1886, 1887, 1888, 1889, &
31 1890, 1891, 1892, 1893, 1894, 1895, 1896, &
32 1897, 1898, 1899, 1900, 1901, 1902, 1903, &
33 1904, 1905, 1906, 1907, 1908, 1909, 1910, &
34 1911, 1912, 1913, 1914, 1915, 1916, 1917, &
35 1918, 1919, 1920, 1921, 1922, 1923, 1924, &
36 1925, 1926, 1927, 1928, 1929, 1930, 1931, &
37 1932, 1933, 1934, 1935, 1936, 1937, 1938, &
38 1939, 1940, 1941, 1942, 1943, 1944, 1945, &
39 1946, 1947, 1948, 1949, 1950, 1951, 1952, &
40 1953, 1954, 1955, 1956, 1957, 1958, 1959, &
41 1960, 1961, 1962, 1963, 1964, 1965, 1966, &
42 1967, 1968, 1969, 1970, 1971, 1972, 1973, &
43 1974, 1975, 1976, 1977, 1978, 1979, 1980, &
44 1981, 1982, 1983, 1984, 1985, 1986, 1987, &
45 1988, 1989, 1990, 1991, 1992, 1993, 1994, &
46 1995, 1996, 1997, 1998, 1999, 2000, 2001, &
47 2002, 2003, 2004, 2005, 2006, 2007, 2008, &
48 2009, 2010, 2011, 2012, 2013, 2014, 2015, &
49 2016, 2017, 2018, 2019, 2020, 2021, 2022, &
50 2023, 2024, 2025, 2026, 2027, 2028, 2029, &
51 2030, 2031, 2032, 2033, 2034, 2035, 2036, &
52 2037, 2038, 2039, 2040, 2041, 2042, 2043, &
53 2044, 2045, 2046, 2047, 2048, 2049, 2050, &
54 2051, 2052, 2053, 2054, 2055, 2056, 2057, &
55 2058, 2059, 2060, 2061, 2062, 2063, 2064, &
56 2065, 2066, 2067, 2068, 2069, 2070, 2071, &
57 2072, 2073, 2074, 2075, 2076, 2077, 2078, &
58 2079, 2080, 2081, 2082, 2083, 2084, 2085, &
59 2086, 2087, 2088, 2089, 2090, 2091, 2092, &
60 2093, 2094, 2095, 2096, 2097, 2098, 2099, &
62 real :: co2r_CAM(1:VARCAM_in_years) = &
63 (/ 289.263, 289.263, 289.416, 289.577, 289.745, 289.919, 290.102, &
64 290.293, 290.491, 290.696, 290.909, 291.129, 291.355, 291.587, 291.824, &
65 292.066, 292.313, 292.563, 292.815, 293.071, 293.328, 293.586, 293.843, &
66 294.098, 294.35, 294.598, 294.842, 295.082, 295.32, 295.558, 295.797, &
67 296.038, 296.284, 296.535, 296.794, 297.062, 297.338, 297.62, 297.91, &
68 298.204, 298.504, 298.806, 299.111, 299.419, 299.729, 300.04, 300.352, &
69 300.666, 300.98, 301.294, 301.608, 301.923, 302.237, 302.551, 302.863, &
70 303.172, 303.478, 303.779, 304.075, 304.366, 304.651, 304.93, 305.206, &
71 305.478, 305.746, 306.013, 306.28, 306.546, 306.815, 307.087, 307.365, &
72 307.65, 307.943, 308.246, 308.56, 308.887, 309.228, 309.584, 309.956, &
73 310.344, 310.749, 311.172, 311.614, 312.077, 312.561, 313.068, 313.599, &
74 314.154, 314.737, 315.347, 315.984, 316.646, 317.328, 318.026, 318.742, &
75 319.489, 320.282, 321.133, 322.045, 323.021, 324.06, 325.155, 326.299, &
76 327.484, 328.698, 329.933, 331.194, 332.499, 333.854, 335.254, 336.69, &
77 338.15, 339.628, 341.125, 342.65, 344.206, 345.797, 347.397, 348.98, &
78 350.551, 352.1, 354.3637, 355.7772, 357.1601, 358.5306, 359.9046, &
79 361.4157, 363.0445, 364.7761, 366.6064, 368.5322, 370.534, 372.5798, &
80 374.6564, 376.7656, 378.9087, 381.0864, 383.2994, 385.548, 387.8326, &
81 390.1536, 392.523, 394.9625, 397.4806, 400.075, 402.7444, 405.4875, &
82 408.3035, 411.1918, 414.1518, 417.1831, 420.2806, 423.4355, 426.6442, &
83 429.9076, 433.2261, 436.6002, 440.0303, 443.5168, 447.06, 450.6603, &
84 454.3059, 457.9756, 461.6612, 465.3649, 469.0886, 472.8335, 476.6008, &
85 480.3916, 484.2069, 488.0473, 491.9184, 495.8295, 499.7849, 503.7843, &
86 507.8278, 511.9155, 516.0476, 520.2243, 524.4459, 528.7127, 533.0213, &
87 537.3655, 541.7429, 546.1544, 550.6005, 555.0819, 559.5991, 564.1525, &
88 568.7429, 573.3701, 578.0399, 582.7611, 587.5379, 592.3701, 597.2572, &
89 602.1997, 607.1975, 612.2507, 617.3596, 622.524, 627.7528, 633.0616, &
90 638.457, 643.9384, 649.505, 655.1568, 660.8936, 666.7153, 672.6219, &
91 678.6133, 684.6945, 690.8745, 697.1569, 703.5416, 710.0284, 716.6172, &
92 723.308, 730.1008, 736.9958, 743.993, 751.0975, 758.3183, 765.6594, &
93 773.1207, 780.702, 788.4033, 796.2249, 804.1667, 812.2289, 820.4118, &
96 PUBLIC :: read_CAMgases
100 SUBROUTINE read_CAMgases(yr, julian, READtrFILE, model, co2vmr, n2ovmr, ch4vmr, cfc11vmr, cfc12vmr)
102 INTEGER, INTENT(IN) :: yr
103 REAL, INTENT(IN) :: julian
104 CHARACTER(LEN=*), INTENT(IN) :: model ! Radiation scheme name
105 REAL(r8), INTENT(OUT) :: co2vmr, n2ovmr, ch4vmr, cfc11vmr, cfc12vmr
106 LOGICAL :: READtrFILE
110 INTEGER :: yearIN, found_yearIN, iyear,yr1,yr2
111 INTEGER :: mondata(1:cyr)
112 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
113 INTEGER, EXTERNAL :: get_unused_unit
115 INTEGER :: istatus, iunit, idata
116 INTEGER, SAVE :: max_years
117 integer :: nyrm, nyrp, njulm, njulp
119 CHARACTER(LEN=256) :: message
120 INTEGER :: monday(13)=(/0,31,28,31,30,31,30,31,31,30,31,30,31/)
121 INTEGER :: mondayi(13)
122 INTEGER :: my1,my2,my3, tot_valid
124 IF ( READtrFILE ) THEN
126 INQUIRE(FILE='CAMtr_volume_mixing_ratio', EXIST=exists)
129 iunit = get_unused_unit()
130 IF ( iunit <= 0 ) THEN
131 IF ( wrf_dm_on_monitor() ) THEN
132 CALL wrf_error_fatal('Error in module_ra_rrtm: could not find a free Fortran unit.')
136 ! Read volume mixing ratio
137 OPEN(UNIT=iunit, FILE='CAMtr_volume_mixing_ratio', FORM='formatted', &
138 STATUS='old', IOSTAT=istatus)
140 IF (istatus == 0) THEN
141 ! Ignore first two lines which constitute a header
142 READ(UNIT=iunit, FMT='(1X)')
143 READ(UNIT=iunit, FMT='(1X)')
147 DO WHILE (istatus == 0)
148 READ(UNIT=iunit, FMT='(I4, 1x, F8.3,1x, 4(F10.3,1x))', IOSTAT=istatus) &
149 yrdata(idata), co2r(idata), n2or(idata), ch4r(idata), cfc11r(idata), &
156 IF ( wrf_dm_on_monitor() ) THEN
157 WRITE(message,*)'Climate GHG input from file from year ',yrdata(1),' to ',yrdata(idata-2)
158 CALL wrf_message( message)
159 WRITE(message,*)'CO2 range = ',co2r (1),co2r (idata-2),' ppm'
160 CALL wrf_message( message)
161 WRITE(message,*)'N2O range = ',n2or (1),n2or (idata-2),' ppb'
162 CALL wrf_message( message)
163 WRITE(message,*)'CH4 range = ',ch4r (1),ch4r (idata-2),' ppb'
164 CALL wrf_message( message)
165 WRITE(message,*)'CFC11 range = ',cfc11r(1),cfc11r(idata-2),' ppt'
166 CALL wrf_message( message)
167 WRITE(message,*)'CFC12 range = ',cfc12r(1),cfc12r(idata-2),' ppt'
168 CALL wrf_message( message)
171 IF (istatus /= -1) THEN
172 WRITE(message,*) " Not normal ending of CAMtr_volume_mixing_ratio file"
173 call wrf_error_fatal( message)
175 max_years = idata - 2
177 ! Calculate the julian day for each month.
180 MY1=MOD(yrdata(idata),4)
181 MY2=MOD(yrdata(idata),100)
182 MY3=MOD(yrdata(idata),400)
183 IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0) mondayi(3)=29
184 juldata(idata) = sum(mondayi(1:mondata(idata)))+real(mondayi(mondata(idata)+1))/2.-0.5 ! 1st Jan 00:00 = 0 julian day
188 max_years = VARCAM_in_years ! Set max number of years to the table size
190 ! For CAM model, recovers original data sets.
191 IF (model .eq. "CAM") THEN
192 yrdata(1:max_years) = yr_CAM
193 co2r(1:max_years) = co2r_CAM
195 CALL wrf_error_fatal("CLWRF: 'CAMtr_volume_mixing_ratio' does not exist")
198 ENDIF ! CAMtr_volume_mixing_ratio exists
200 ENDIF ! File already opened and read
202 ! We have already read the data once. Now we process it to find the valid value
203 ! for this year day combination.
207 DO WHILE (found_yearIN == 0 .and. iyear <= cyr)
208 IF (yrdata(iyear) .GT. yr ) THEN
211 ELSE IF ((yrdata(iyear) .EQ. yr) .AND. (juldata(iyear) .GT. julian)) THEN
218 ! Prevent yr > last year in data
219 IF (iyear .ge. max_years) then
224 IF (found_yearIN .NE. 0 ) THEN
225 if (yearIN .eq. 1) yearIN = yearIN + 1 ! To take 2 first lines of the file
226 nyrm = yrdata(yearIN-1)
227 njulm = juldata(yearIN-1)
228 nyrp = yrdata(yearIN)
229 njulp = juldata(yearIN)
237 if (found_yearIN /= 0) then
238 ! Interpolate data only if we have at least 2 valid concentrations.
239 tot_valid = count(co2r(1:max_years) > 0)
240 IF (tot_valid >= 2 ) THEN
241 CALL valid_years(yearIN, co2r, max_years,yr1, yr2)
243 ! Set nyrm, njulm, nyrp, njulp
249 CALL interpolate_CAMgases(yr, julian, nyrm, njulm, yr1, yr2, nyrp, njulp, max_years, co2r , co2vmr )
250 CALL interpolate_CAMgases(yr, julian, nyrm, njulm, yr1, yr2, nyrp, njulp, max_years, n2or , n2ovmr )
251 CALL interpolate_CAMgases(yr, julian, nyrm, njulm, yr1, yr2, nyrp, njulp, max_years, ch4r , ch4vmr )
252 CALL interpolate_CAMgases(yr, julian, nyrm, njulm, yr1, yr2, nyrp, njulp, max_years, cfc11r, cfc11vmr)
253 CALL interpolate_CAMgases(yr, julian, nyrm, njulm, yr1, yr2, nyrp, njulp, max_years, cfc12r, cfc12vmr)
257 ! Verification of interpolated values. In case of no value
258 ! original values extracted from ghg_surfvals.F90 module
260 IF ( (co2vmr < 0.) .OR. &
263 (cfc11vmr < 0.) .OR. &
264 (cfc12vmr < 0.) .OR. &
265 (found_yearIN == 0) ) THEN
266 CALL orig_val("CO2" ,model,co2vmr )
267 CALL orig_val("N2O" ,model,n2ovmr )
268 CALL orig_val("CH4" ,model,ch4vmr )
269 CALL orig_val("CFC11",model,cfc11vmr)
270 CALL orig_val("CFC12",model,cfc12vmr)
272 ! If extrapolation, need to bound the data to pre-industrial concentrations
273 if (co2vmr < 270.) co2vmr = 270.
274 if (n2ovmr < 270.) n2ovmr = 270.
275 if (ch4vmr < 700.) ch4vmr = 700.
276 co2vmr =co2vmr *1.e-06
277 n2ovmr =n2ovmr *1.e-09
278 ch4vmr =ch4vmr *1.e-09
279 cfc11vmr=cfc11vmr*1.e-12
280 cfc12vmr=cfc12vmr*1.e-12
283 END SUBROUTINE read_CAMgases
285 SUBROUTINE valid_years(yearIN, gas, tot_years, yr1, yr2)
288 INTEGER, INTENT(IN) :: yearIN, tot_years
289 INTEGER, INTENT(OUT) :: yr2, yr1
290 REAL(r8), INTENT(IN) :: gas(:)
293 INTEGER :: yr_loc, idata
300 ! If all valid dates are > yearIN then find the 2 lowest dates with
302 IF (count(gas(1:yr_loc-1) > 0.) == 0) THEN
303 DO idata = yr_loc, tot_years-1
304 IF (gas(idata) > 0.) THEN
309 DO idata = yr1+1, tot_years
310 IF (gas(idata) > 0.) THEN
315 ELSE ! There is at least 1 valid year below yearIN
316 IF (gas(yr_loc) < 0.) THEN
318 ! Find the closest valid year, look for higher year first
319 IF (any(gas(yr_loc:tot_years) > 0)) THEN
320 DO idata=yr_loc+1, tot_years
321 IF (gas(idata) > 0) THEN
326 ELSE ! Need to take lower years and extrapolate data.
327 DO idata=yr_loc-1,2,-1
328 IF (gas(idata) > 0) THEN
336 yr_loc = min(yr_loc-1, yr2-1)
337 IF (gas(yr_loc) < 0.) THEN
339 ! Find the closest valid year, lower than yr1
340 IF (any(gas(1:yr_loc-1) > 0)) THEN
341 DO idata=yr_loc-1,1,-1
342 IF (gas(idata) > 0) THEN
347 ELSE ! Need to take higher years and extrapolate data.
348 print*, 'Problem: this case should never happen'
350 ELSE ! Then yr1 is yr_loc (first valid date < yr2)
354 END SUBROUTINE valid_years
356 SUBROUTINE interpolate_CAMgases(yr, julian, yeari, juli, yr1, yr2, yearf, julf, max_years, gas, interp_gas)
358 ! These subroutine interpolates a trace gas concentration from a non-homogeneously
359 ! distributed gas concentration evolution
360 INTEGER, INTENT (IN) :: yr, yeari, yr1, yr2, yearf, juli, julf, max_years
361 REAL, INTENT (IN) :: julian
362 REAL(r8), DIMENSION(max_years), INTENT(IN):: gas
363 REAL(r8), INTENT (OUT) :: interp_gas
365 REAL(r8) :: yearini, yearend, gas1, gas2
366 REAL :: doymodel, doydatam, doydatap, &
367 deltat, fact1, fact2, x
368 INTEGER :: ny, my1,my2,my3,nday, maxyear, minyear
371 ! Add support for leap-years
373 ! Find smallest and largest year: yearf >= yeari since the file is ordered with increasing dates.
374 minyear = MIN(yr,yeari)
375 maxyear = MAX(yr,yearf)
377 ! Initialise with the day in the year for each date.
382 ! Calculate the julian day (day 0 = 1 Jan minyear at 00:00)
383 DO ny=minyear, maxyear-1
389 IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0) nday=366
403 fact1 = (fact1 - x)/deltat
404 fact2 = (x - fact2)/deltat
406 interp_gas = gas(yr1)*fact1+gas(yr2)*fact2
408 IF (interp_gas .LT. 0. ) THEN
412 END SUBROUTINE interpolate_CAMgases
414 SUBROUTINE interpolate_lin(x1,y1,x2,y2,x0,y)
416 ! Program to interpolate values y=a+b*x with:
421 REAL, INTENT (IN) :: x1,x2,x0
422 REAL(r8), INTENT (IN) :: y1,y2
423 REAL(r8), INTENT (OUT) :: y
438 END SUBROUTINE interpolate_lin
441 SUBROUTINE orig_val(tracer,model,out)
443 CHARACTER(LEN=*), INTENT(IN) :: tracer ! The trace gas name
444 CHARACTER(LEN=*), INTENT(IN) :: model ! The radiation scheme name
445 REAL(r8), INTENT(INOUT) :: out ! The output value
447 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
448 CHARACTER(LEN=256) :: message
451 IF (model .eq. "CAM") THEN
452 IF (tracer .eq. "CO2") THEN
455 ELSE IF (tracer .eq. "N2O") THEN
458 ELSE IF (tracer .eq. "CH4") THEN
461 ELSE IF (tracer .eq. "CFC11") THEN
464 ELSE IF (tracer .eq. "CFC12") THEN
468 IF ( wrf_dm_on_monitor() ) THEN
469 WRITE(message,*) 'CLWRF : Trace gas ',tracer,' not valid for scheme ',model
470 CALL wrf_error_fatal(message)
474 ELSE IF ((model .eq. "RRTM") .or. &
475 (model .eq. "RRTMG")) THEN
476 IF (tracer .eq. "CO2") THEN
479 ELSE IF (tracer .eq. "N2O") THEN
482 ELSE IF (tracer .eq. "CH4") THEN
485 ELSE IF (tracer .eq. "CFC11") THEN
488 ELSE IF (tracer .eq. "CFC12") THEN
492 IF ( wrf_dm_on_monitor() ) THEN
493 WRITE(message,*) 'CLWRF : Trace gas ',tracer,' not valid for scheme ',model
494 CALL wrf_error_fatal(message)
499 IF ( wrf_dm_on_monitor() ) THEN
500 WRITE(message,*) 'CLWRF not implemented for the ',model,' radiative scheme.'
501 CALL wrf_error_fatal(message)
505 END SUBROUTINE orig_val
507 END MODULE module_ra_clWRF_support