Thursday, June 14, 2012

5177.txt

date: Wed, 22 Dec 1999 10:35:24 -0700
from: Lisa Butler <lbutleratXYZxyzr.edu>
subject: Message from Wigley
to: m.hulmeatXYZxyz.ac.uk

Dear Mike,

Here is the new MAGICC code (a revision of MAG.FOR). By reading the
header notes you can see what has been done. If the two accompanying
CFG files are "plugged in", and the FORTRAN code is compiled, it should
run with the existing (2.32) shell. No outputs have been changed, so it
should still drive SCENGEN as before. Do you think you could remind me
what Olga is doing to SCENGEN, please? One simple item is that the
default emissions files that come up should be IS92a (reference) and
IS92c (policy). At the moment 92d comes up, which is a little silly
because this is rarely, if ever, used.

Also enclosed are two GAS.EM files. One is IS95A-S (same as IS92a),
which seems to be the only IS92 file I have that has the regional SO2
emissions breakdown. Do you have breakdowns for the other IS92
scenarios??? The second GAS.EM file is WRE550.EM. The CO2 emissions in
this file are set up to reproduce the WRE550 concentration profile.
There is a new switch on the header to this file, explained in the
FORTRAN code header notes. This file goes out to 2500. Gases other
than CO2 follow IS92a to 2100, and then have constant emissions. It
would be better to have these post-2100 emissions set up to achieve
stabilization too. This requires some work, which I hope to do soon.

The CFG files are set up for a SCENGEN run (i.e., 20 loops through the
code), and to use DT=0.5yr. On my (slow) laptop, this takes about 4
seconds to run. Using DT=0.1yr (which is what we usually use for runs
of MAGICC alone, as in the SAR) the code takes about 9 seconds to run.
Apropos this, there are a bunch of comparison files (.TXT) that show how
this version of the code compares with the results used in the SAR (the
differences are line thickness only). Two of these files show what
difference it makes to use DT=0.5yr instead of 0.1yr. You need to
decide on accuracy versus running time in choosing what time step to use
for UNDP. My choice would be 0.5yr, since even 0.1yr doesn't reproduce
the SAR precisely.

For UNDP users, some simplifications to the shell are warranted (or
necessary). You already have some of those underway. I have marked by
an "OK" in MAGUSER.CFG what I think these users should be able to
change. The residual user options from your instructions to Mike Salmon
on what options should be turned *off* in 2.32, should match what I've
indicated should be left *on*.

There are a lot of things that the code can do that UNDP users will not
be able to access at all, since they can only change MAGUSER through the
shell. They could still get at MAGEXTRA from outside of MAGICC, but I
doubt that they would want to do this. Still, for safety's sake, the
MAGEXTRA file that goes along with the software package should have the
explanatory text expunged from it.

I have written a "handbook" for this version, but it has not yet been
word processed. I am afraid I must insist that the Handbook or
"workbook" that goes to UNDP is what *I* have produced. You can modify
it, of course, since there may be aspects about UNDP's specific needs
that I am not aware of. Furthermore, I also insist that *I* am first
author for this publication. The main reason for this should be clear
from my "credit" email of yesterday. Another reason is that I should
take (as well as primary credit), the primary responsibility --- i.e.,
the buck should stop with me. Of course, I appreciate all you have done
and are doing.

My Handbook is much more that just a set of examples (which is what I
recall you were planning to do), although it does have a number of
examples designed to illustrate how MAGICC-SCENGEN can be used. I
figured that it was necessary for *me* to do this, because there are
idiosyncrasies, limitations and strengths in the code that only I know
about.

I can't get this to you until about January 5, because NCAR is closed
until January 3. The handwritten version is too scrappy to send. This
document will still need input from you and/or Mike Salmon, since it
does not yet include diagrams --- just text explaining what these are.
These are all simple dumps from MAGICC-SCENGEN, but you will need to do
these at your end.

Along with this, some updates of the help files in SCENGEN re
necessary. Not least, these could be cross-referenced to the Handbook.
I believe you are doing some of this already; so the most efficient
thing to do is for me to have a go at the next iteration. I've noticed
items that need improvement, but you'll probably cover most of those.

My apologies for being slow with this, and for any confusion regarding
the Handbook. I didn't realize how soon the deadline was. If you need
to call me about any of the enclosed code, my home number is
+303-494-4372. We plan to be home for all except Dec. 30 through
January 2.

Best wishes for Xmas and the New Year. Try to relax!

Cheers,
Tom

--
*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
Lisa Butler
Administrative Assistant
ACACIA Program Office
National Center for Atmospheric Research
P.O. Box 3000
Boulder, CO 80307-3000
(303)497-2697 PH E-mail: lbutleratXYZxyzr.edu
(303)497-2699 FAX http://www.cgd.ucar.edu/cas/ACACIA
*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*51 1
WRE550
WRE550 STABILIZATION FOR CO2, ELSE IS92a
Year Foss CO2 Defo CO2 CH4 N2O SO2,1 SO2,2 SO2,3
(Pg C) (Pg C) (Tg) (Tg N) (Tg S) (Tg S) (Tg S)
1990 6.109 1.300 506.000 12.900 .000 .000 .000
1995 6.604 1.300 525.000 13.300 2.000 .000 .000
2000 7.101 1.300 545.000 13.800 3.000 .000 .000
2005 7.938 1.260 568.000 14.100 10.000 .000 .000
2010 8.679 1.220 588.000 14.500 17.000 .000 .000
2015 9.197 1.180 611.000 14.900 24.000 .000 .000
2020 9.629 1.140 632.000 15.400 34.000 .000 .000
2025 9.938 1.100 659.000 15.800 43.000 .000 .000
2030 10.140 1.040 684.200 15.960 49.800 .000 .000
2035 10.213 0.980 709.400 16.120 56.600 .000 .000
2040 10.177 0.920 734.600 16.280 63.400 .000 .000
2045 10.038 0.860 759.800 16.440 70.200 .000 .000
2050 9.815 0.800 785.000 16.600 77.000 .000 .000
2055 9.599 0.670 797.000 16.620 76.800 .000 .000
2060 9.348 0.540 809.000 16.640 76.600 .000 .000
2065 9.076 0.410 821.000 16.660 76.400 .000 .000
2070 8.790 0.280 833.000 16.680 76.200 .000 .000
2075 8.506 0.150 845.000 16.700 76.000 .000 .000
2080 8.130 0.100 859.400 16.760 75.000 .000 .000
2085 7.772 0.050 873.800 16.820 74.000 .000 .000
2090 7.426 0.000 888.200 16.880 73.000 .000 .000
2095 7.100 -0.050 902.600 16.940 72.000 .000 .000
2100 6.795 0.000 917.000 17.000 71.000 .000 .000
2110 6.187 0.000 917.000 17.000 71.000 .000 .000
2120 5.655 0.000 917.000 17.000 71.000 .000 .000
2130 5.195 0.000 917.000 17.000 71.000 .000 .000
2140 4.795 0.000 917.000 17.000 71.000 .000 .000
2150 4.450 0.000 917.000 17.000 71.000 .000 .000
2160 4.235 0.000 917.000 17.000 71.000 .000 .000
2170 4.047 0.000 917.000 17.000 71.000 .000 .000
2180 3.880 0.000 917.000 17.000 71.000 .000 .000
2190 3.727 0.000 917.000 17.000 71.000 .000 .000
2200 3.588 0.000 917.000 17.000 71.000 .000 .000
2210 3.460 0.000 917.000 17.000 71.000 .000 .000
2220 3.342 0.000 917.000 17.000 71.000 .000 .000
2230 3.232 0.000 917.000 17.000 71.000 .000 .000
2240 3.129 0.000 917.000 17.000 71.000 .000 .000
2250 3.033 0.000 917.000 17.000 71.000 .000 .000
2260 2.943 0.000 917.000 17.000 71.000 .000 .000
2270 2.858 0.000 917.000 17.000 71.000 .000 .000
2280 2.778 0.000 917.000 17.000 71.000 .000 .000
2290 2.703 0.000 917.000 17.000 71.000 .000 .000
2300 2.631 0.000 917.000 17.000 71.000 .000 .000
2325 2.468 0.000 917.000 17.000 71.000 .000 .000
2350 2.324 0.000 917.000 17.000 71.000 .000 .000
2375 2.195 0.000 917.000 17.000 71.000 .000 .000
2400 2.080 0.000 917.000 17.000 71.000 .000 .000
2425 1.975 0.000 917.000 17.000 71.000 .000 .000
2450 1.880 0.000 917.000 17.000 71.000 .000 .000
2475 1.793 0.000 917.000 17.000 71.000 .000 .000
2500 1.713 0.000 917.000 17.000 71.000 .000 .000
12/21/99
IPCC SAR RESULTS (SAR) VERSUS NEW MAG-F
INCLUDES SO2 EMISSIONS
TEMPERATURES FOR DT2X=2.5degC WITH AEROSOLS

YEAR SAR95A NEW95A SAR95B NEW95B SAR95C NEW95C SAR95D NEW95D SAR95E NEW95E SAR95F NEW95F
1990 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000
1995 .0649 .0641 .0650 .0642 .0785 .0779 .0756 .0748 .0657 .0649 .0683 .0676
2000 .1389 .1377 .1418 .1405 .1516 .1509 .1486 .1475 .1348 .1336 .1390 .1378
2005 .2020 .2004 .2065 .2050 .2197 .2193 .2168 .2159 .2038 .2023 .2055 .2040
2010 .2669 .2652 .2719 .2701 .2810 .2812 .2881 .2878 .2709 .2692 .2777 .2760
2015 .3360 .3341 .3380 .3362 .3409 .3418 .3615 .3619 .3393 .3375 .3493 .3475
2020 .4015 .3996 .3998 .3979 .3974 .3992 .4330 .4344 .4112 .4094 .4217 .4198
2025 .4730 .4710 .4623 .4603 .4540 .4568 .5058 .5081 .4802 .4783 .4984 .4966
2030 .5546 .5527 .5411 .5392 .5249 .5289 .5873 .5907 .5637 .5619 .5906 .5889
2035 .6403 .6385 .6253 .6234 .5962 .6013 .6693 .6738 .6548 .6531 .6880 .6863
2040 .7281 .7264 .7116 .7099 .6645 .6706 .7496 .7552 .7512 .7496 .7876 .7861
2045 .8171 .8154 .7990 .7974 .7289 .7360 .8274 .8341 .8517 .8503 .8885 .8872
2050 .9082 .9068 .8884 .8869 .7904 .7986 .9036 .9113 .9568 .9556 .9913 .9902
2055 1.0164 1.0151 .9953 .9940 .8595 .8685 .9807 .9894 1.0791 1.0782 1.1062 1.1053
2060 1.1273 1.1262 1.1052 1.1041 .9264 .9359 1.0556 1.0653 1.2084 1.2077 1.2260 1.2253
2065 1.2378 1.2368 1.2147 1.2137 .9892 .9989 1.1281 1.1385 1.3418 1.3413 1.3485 1.3481
2070 1.3473 1.3466 1.3232 1.3224 1.0478 1.0573 1.1985 1.2095 1.4784 1.4782 1.4731 1.4731
2075 1.4560 1.4554 1.4307 1.4300 1.1024 1.1114 1.2670 1.2782 1.6184 1.6184 1.6000 1.6003
2080 1.5666 1.5662 1.5367 1.5361 1.1502 1.1583 1.3313 1.3427 1.7794 1.7796 1.7449 1.7455
2085 1.6797 1.6794 1.6440 1.6435 1.1940 1.2013 1.3942 1.4056 1.9492 1.9496 1.8964 1.8973
2090 1.7948 1.7948 1.7529 1.7526 1.2349 1.2413 1.4563 1.4675 2.1240 2.1247 2.0514 2.0526
2095 1.9120 1.9121 1.8636 1.8633 1.2731 1.2788 1.5178 1.5288 2.3027 2.3037 2.2089 2.2106
2100 2.0308 2.0311 1.9757 1.9756 1.3089 1.3140 1.5787 1.5895 2.4841 2.4854 2.3686 2.3707
12/21/99
IPCC SAR RESULTS (SAR) VERSUS NEW MAG-F
ABSOLUTE EXTREMES FOR TEMPERATURES
LOW IS IS95C WITH CONSTANT ESO2, DT2X=1.5
A IS IS95A WITH DT2X=2.5, MID MELT
HIGH IS IS95E WITH CONSTANT ESO2, DT2X=4.5
FOR IS95A, WITH INCLUDES ESO2, NONE MEANS CONSTANT ESO2

YEAR C-LOW C-LOW A-WITH A-WITH A-NONE A-NONE E-HIGH E-HIGH
YEAR SAR NEW SAR NEW SAR NEW SAR NEW
1990 .0000 .0000 .0000 .0000 .0000 .0000 .0000 .0000
1995 .0483 .0479 .0649 .0641 .0717 .0710 .0994 .0983
2000 .0975 .0972 .1389 .1377 .1519 .1506 .2155 .2137
2005 .1492 .1491 .2020 .2004 .2397 .2382 .3482 .3460
2010 .2019 .2022 .2669 .2652 .3344 .3327 .4961 .4935
2015 .2551 .2560 .3360 .3341 .4348 .4329 .6584 .6556
2020 .3083 .3099 .4015 .3996 .5404 .5385 .8334 .8305
2025 .3618 .3642 .4730 .4710 .6510 .6491 1.0199 1.0170
2030 .4143 .4175 .5546 .5527 .7655 .7636 1.2171 1.2142
2035 .4637 .4676 .6403 .6385 .8822 .8804 1.4226 1.4198
2040 .5097 .5144 .7281 .7264 1.0004 .9986 1.6350 1.6323
2045 .5523 .5577 .8171 .8154 1.1192 1.1176 1.8528 1.8504
2050 .5922 .5983 .9082 .9068 1.2391 1.2377 2.0761 2.0739
2055 .6284 .6350 1.0164 1.0151 1.3584 1.3571 2.3037 2.3017
2060 .6599 .6668 1.1273 1.1262 1.4753 1.4742 2.5355 2.5339
2065 .6872 .6942 1.2378 1.2368 1.5901 1.5891 2.7715 2.7701
2070 .7107 .7175 1.3473 1.3466 1.7028 1.7021 3.0110 3.0100
2075 .7308 .7371 1.4560 1.4554 1.8139 1.8133 3.2538 3.2531
2080 .7479 .7535 1.5666 1.5662 1.9246 1.9242 3.5003 3.4999
2085 .7625 .7675 1.6797 1.6794 2.0366 2.0364 3.7513 3.7512
2090 .7750 .7794 1.7948 1.7948 2.1503 2.1502 4.0065 4.0068
2095 .7855 .7893 1.9120 1.9121 2.2656 2.2657 4.2655 4.2662
2100 .7940 .7974 2.0308 2.0311 2.3822 2.3825 4.5274 4.5285
12/21/99
MAG-F RESULTS FOR DT=0.1 VS DT=0.5
INCLUDES SO2 EMISSIONS
TEMPERATURES FOR DT2X=2.5degC WITH AEROSOLS

YEAR 95A-0.5 95A-0.1 95C-0.5 95C-0.1 95E-0.5 95E-0.1
1990 .0000 .0000 .0000 .0000 .0000 .0000
1995 .0642 .0641 .0777 .0779 .0651 .0649
2000 .1378 .1377 .1509 .1509 .1338 .1336
2005 .2007 .2004 .2195 .2193 .2025 .2023
2010 .2654 .2652 .2814 .2812 .2695 .2692
2015 .3344 .3341 .3421 .3418 .3377 .3375
2020 .3998 .3996 .3995 .3992 .4096 .4094
2025 .4710 .4710 .4569 .4568 .4784 .4783
2030 .5525 .5527 .5287 .5289 .5617 .5619
2035 .6382 .6385 .6012 .6013 .6527 .6531
2040 .7260 .7264 .6706 .6706 .7491 .7496
2045 .8151 .8154 .7361 .7360 .8497 .8503
2050 .9060 .9068 .7985 .7986 .9546 .9556
2055 1.0140 1.0151 .8682 .8685 1.0768 1.0782
2060 1.1251 1.1262 .9358 .9359 1.2063 1.2077
2065 1.2357 1.2368 .9988 .9989 1.3398 1.3413
2070 1.3455 1.3466 1.0574 1.0573 1.4765 1.4782
2075 1.4543 1.4554 1.1115 1.1114 1.6164 1.6184
2080 1.5650 1.5662 1.1586 1.1583 1.7772 1.7796
2085 1.6782 1.6794 1.2017 1.2013 1.9470 1.9496
2090 1.7935 1.7948 1.2418 1.2413 2.1220 2.1247
2095 1.9107 1.9121 1.2793 1.2788 2.3008 2.3037
2100 2.0297 2.0311 1.3145 1.3140 2.4825 2.4854
12/21/99
IPCC SAR RESULTS (SAR) VERSUS NEW MAG-F
INCLUDES SO2 EMISSIONS
SEA LEVEL FOR DT2X=2.5degC WITH AEROSOLS, MID MELT

YEAR SAR95A NEW95A SAR95B NEW95B SAR95C NEW95C SAR95D NEW95D SAR95E NEW95E SAR95F NEW95F
1990 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00
1995 .95 .92 .95 .92 .97 .94 .96 .94 .95 .92 .95 .93
2000 2.04 1.99 2.05 2.00 2.08 2.03 2.07 2.02 2.04 1.99 2.05 2.00
2005 3.27 3.19 3.28 3.20 3.32 3.24 3.30 3.22 3.28 3.20 3.29 3.21
2010 4.61 4.49 4.62 4.51 4.65 4.55 4.65 4.54 4.66 4.54 4.67 4.55
2015 6.07 5.92 6.08 5.93 6.07 5.93 6.11 5.97 6.16 6.01 6.18 6.03
2020 7.64 7.45 7.64 7.45 7.56 7.39 7.68 7.51 7.80 7.61 7.82 7.63
2025 9.33 9.10 9.29 9.06 9.11 8.92 9.36 9.16 9.55 9.33 9.60 9.37
2030 11.15 10.88 11.07 10.80 10.77 10.55 11.15 10.92 11.46 11.19 11.55 11.28
2035 13.12 12.81 12.99 12.68 12.52 12.28 13.06 12.80 13.53 13.22 13.66 13.35
2040 15.22 14.87 15.05 14.69 14.35 14.09 15.07 14.78 15.77 15.41 15.95 15.58
2045 17.45 17.05 17.23 16.84 16.25 15.96 17.17 16.85 18.17 17.77 18.38 17.97
2050 19.82 19.37 19.55 19.10 18.19 17.89 19.35 19.00 20.75 20.29 20.97 20.50
2055 22.34 21.83 22.02 21.52 20.20 19.87 21.60 21.23 23.52 23.01 23.72 23.20
2060 25.02 24.46 24.65 24.10 22.27 21.91 23.93 23.52 26.51 25.94 26.66 26.08
2065 27.85 27.24 27.43 26.82 24.38 24.00 26.31 25.88 29.70 29.07 29.78 29.14
2070 30.76 30.15 30.32 29.68 26.52 26.11 28.73 28.28 32.93 32.39 32.84 32.33
2075 33.57 33.12 33.08 32.60 28.67 28.23 30.97 30.69 36.21 35.83 35.98 35.62
2080 36.50 36.15 35.94 35.58 30.63 30.33 33.25 33.09 39.71 39.41 39.31 39.01
2085 39.54 39.24 38.91 38.61 32.58 32.37 35.56 35.47 43.45 43.15 42.84 42.54
2090 42.71 42.39 42.00 41.68 34.51 34.38 37.90 37.85 47.24 47.05 46.36 46.18
2095 45.79 45.60 45.03 44.81 36.44 36.34 40.27 40.21 51.16 51.08 50.00 49.93
2100 48.96 48.87 48.10 47.99 38.35 38.26 42.54 42.57 55.29 55.25 53.82 53.79
12/21/99
IPCC SAR RESULTS (SAR) VERSUS NEW MAG-F
ABSOLUTE EXTREMES FOR SEA LEVEL RISE
LOW IS IS95C WITH CONSTANT ESO2, DT2X=1.5, LOW MELT
A IS IS95A WITH DT2X=2.5, MID MELT
HIGH IS IS95E WITH CONSTANT ESO2, DT2X=4.5, HIGH MELT
FOR IS95A, WITH INCLUDES ESO2, NONE MEANS CONSTANT ESO2

YEAR C-LOW C-LOW A-WITH A-WITH A-NONE A-NONE E-HIGH E-HIGH
YEAR SAR NEW SAR NEW SAR NEW SAR NEW
1990 .00 .00 .00 .00 .00 .00 .00 .00
1995 .30 .28 .95 .92 .96 .93 2.01 1.96
2000 .65 .62 2.04 1.99 2.08 2.03 4.39 4.27
2005 1.06 1.03 3.27 3.19 3.38 3.30 7.20 7.00
2010 1.53 1.48 4.61 4.49 4.87 4.75 10.45 10.16
2015 2.05 1.99 6.07 5.92 6.54 6.38 14.18 13.79
2020 2.61 2.54 7.64 7.45 8.40 8.19 18.40 17.90
2025 3.22 3.14 9.33 9.10 10.44 10.19 22.87 22.47
2030 3.87 3.78 11.15 10.88 12.67 12.37 27.53 27.30
2035 4.55 4.44 13.12 12.81 15.09 14.73 32.51 32.35
2040 5.25 5.12 15.22 14.87 17.68 17.26 37.59 37.59
2045 5.95 5.82 17.45 17.05 20.44 19.96 42.92 43.00
2050 6.67 6.52 19.82 19.37 23.36 22.82 48.38 48.55
2055 7.38 7.22 22.34 21.83 26.44 25.83 53.99 54.25
2060 8.09 7.92 25.02 24.46 29.57 28.96 59.79 60.07
2065 8.78 8.59 27.85 27.24 32.58 32.15 65.66 66.02
2070 9.45 9.25 30.76 30.15 35.70 35.37 71.75 72.09
2075 10.10 9.88 33.57 33.12 38.91 38.61 77.90 78.28
2080 10.73 10.49 36.50 36.15 42.15 41.88 84.22 84.58
2085 11.33 11.07 39.54 39.24 45.31 45.17 90.71 91.01
2090 11.90 11.63 42.71 42.39 48.57 48.50 97.29 97.57
2095 12.45 12.16 45.79 45.60 51.93 51.87 104.06 104.26
2100 12.98 12.66 48.96 48.87 55.29 55.28 110.99 111.10
12/21/99
MAG-F RESULTS FOR DT=0.1 VS DT=0.5
INCLUDES SO2 EMISSIONS
SEA LEVEL FOR DT2X=2.5degC WITH AEROSOLS, MID MELT

YEAR 95A-0.5 95A-0.1 95C-0.5 95C-0.1 95E-0.5 95E-0.1
1990 .00 .00 .00 .00 .00 .00
1995 .92 .92 .94 .94 .92 .92
2000 1.98 1.99 2.02 2.03 1.98 1.99
2005 3.18 3.19 3.23 3.24 3.19 3.20
2010 4.48 4.49 4.53 4.55 4.53 4.54
2015 5.90 5.92 5.91 5.93 5.99 6.01
2020 7.43 7.45 7.37 7.39 7.59 7.61
2025 9.07 9.10 8.90 8.92 9.30 9.33
2030 10.85 10.88 10.53 10.55 11.16 11.19
2035 12.77 12.81 12.25 12.28 13.18 13.22
2040 14.83 14.87 14.06 14.09 15.37 15.41
2045 17.01 17.05 15.93 15.96 17.72 17.77
2050 19.32 19.37 17.85 17.89 20.24 20.29
2055 21.78 21.83 19.83 19.87 22.95 23.01
2060 24.41 24.46 21.87 21.91 25.87 25.94
2065 27.18 27.24 23.95 24.00 29.00 29.07
2070 30.08 30.15 26.06 26.11 32.31 32.39
2075 33.05 33.12 28.18 28.23 35.75 35.83
2080 36.08 36.15 30.28 30.33 39.32 39.41
2085 39.16 39.24 32.32 32.37 43.06 43.15
2090 42.31 42.39 34.33 34.38 46.95 47.05
2095 45.51 45.60 36.29 36.34 50.98 51.08
2100 48.78 48.87 38.21 38.26 55.14 55.25
12
IS95A-S
IPCC Scenario 92a modified for 1995 exercise with ESO2 breakdown
Year Foss CO2 Defo CO2 CH4 N2O SO2,1 SO2,2 SO2,3
(Pg C) (Pg C) (Tg) (Tg N) (Tg S) (Tg S) (Tg S)
1990 6.100 1.300 506.000 12.900 0.000 0.000 0.000
1995 6.600 1.300 525.000 13.300 0.900 0.780 0.320
2000 7.100 1.300 545.000 13.800 1.230 1.230 0.540
2005 7.940 1.260 568.000 14.100 3.860 4.180 1.960
2010 8.680 1.220 588.000 14.500 6.150 7.240 3.600
2015 9.420 1.180 611.000 14.900 8.110 10.420 5.470
2020 10.260 1.140 632.000 15.400 10.680 15.030 8.300
2025 11.100 1.100 659.000 15.800 12.470 19.350 11.180
2050 13.700 0.800 785.000 16.600 15.400 37.730 23.870
2075 16.150 0.150 845.000 16.700 13.680 38.000 24.320
2100 20.400 -0.100 917.000 17.000 11.360 36.210 23.430
2200 20.400 -0.100 917.000 17.000 11.360 36.210 23.430
X123451234567890123456789012345678901234567890123456789012345678901234567890
12
IS95A-S
IPCC Scenario 92a modified for 1995 exercise with ESO2 breakdown
Year Foss CO2 Defo CO2 CH4 N2O SO2,1 SO2,2 SO2,3
(Pg C) (Pg C) (Tg) (Tg N) (Tg S) (Tg S) (Tg S)
1990 6.100 1.300 506.000 12.900 0.000 0.000 0.000
1995 6.600 1.300 525.000 13.300 0.900 0.780 0.320
2000 7.100 1.300 545.000 13.800 1.230 1.230 0.540
2005 7.940 1.260 568.000 14.100 3.860 4.180 1.960
2010 8.680 1.220 588.000 14.500 6.150 7.240 3.600
2015 9.420 1.180 611.000 14.900 8.110 10.420 5.470
2020 10.260 1.140 632.000 15.400 10.680 15.030 8.300
2025 11.100 1.100 659.000 15.800 12.470 19.350 11.180
2050 13.700 0.800 785.000 16.600 15.400 37.730 23.870
2075 16.150 0.150 845.000 16.700 13.680 38.000 24.320
2100 20.400 -0.100 917.000 17.000 11.360 36.210 23.430
2200 20.400 -0.100 917.000 17.000 11.360 36.210 23.430
X123451234567890123456789012345678901234567890123456789012345678901234567890
1 * IOLDTZ : IF=1 RETAINS OLD INITIAL OCEAN TEMP PROFILE
0.5 * DT : INTEGRATION TIME STEP FOR CLIMATE MODEL
6.3 * CO2DELQ : QCO2 = CO2DELQ*ALOG(C/C0) : Q2X = CO2DELQ*ALOG(2.0)
1.0 * XKLO : LAND-OCEAN EXCHANGE COEFFICIENT
1.0 * XKNS : N.H.-S.H. EXCHANGE COEFFICIENT
0 * IEXPAN : 0=INCREMENTAL : 1=LUMP EXPANSION ALGORITHM
5 * NOUT : PRINTOUT OPTIONS : SEE BELOW
5 * IEMPRT : PRINTOUT INTERVAL FOR EMISSIONS INPUT
10 * ICO2PRT : PRINTOUT INTERVAL FOR CARBON CYCLE
5 * ICONCPRT : PRINTOUT INTERVAL FOR CONCENTRATIONS
10 * IQDECPRT : PRINTOUT INTERVAL DECADAL FORCING BREAKDOWN
10 * INSLOPRT : PRINTOUT INTERVAL NH/SH, LAND/OCN FORCING BREAKDOWN
5 * IQGASPRT : PRINTOUT INTERVAL FOR GAS BY GAS FORCING BREAKDOWN
5 * IDIS : PRINTOUT INTERVAL FOR *.DIS OUTPUT FILES
0. * TOFFSET : OFFSET FOR TEMPS; MAKES T(KYRREF)=TOFFSET
0.05 * STRATH2O : RATIO DELQH2O/DELQCH4, FOR STRAT H2O FROM CH4 OXIDN
0.02 * DQDCOZ : OZONE DQ/DC; IF =0.0, REMOVES O3 RELATED TO CH4
0.32 * S90OZ : 1990 VALUE OF 'DIRECT' TROPOSPHERIC OZONE FORCING
42.0 * ENAT : NATURAL EMISSIONS IN FORMULA FOR INDIRECT SO4 FORCING
2 * IGHG : IF=0, 1, GE.2 ; GHG-DQ = ZERO, CO2 ONLY, FULL GHG-DQ
0 * ICO2READ : IF.GE.1, READS CO2INPUT.DAT AND OVERWRITES CO2 CONCS
1.00 * CO2SCALE : SCALING FOR QCO2 USED IF ICO2READ=1
12.30 * SO2SCALE : SCALING FOR ESO2 RELT TO ECO2 USED IF ICO2READ=4
0 * ICOMP : IF=1, OVERWRITES QTOT AFTER 1990 WITH COMPOUND DCO2
1.00 * COMPRATE : COMPOUNDING RATE USED POST 1990. USED IF ICOMP=1
0 * IQCONST : IF=1, QTOT(t>1764) = QCONST; IF=2, RAMP TO QCONST
8.734 * QCONST : CONSTANT QTOT USED IF IQCONST.GE.1
0.063 * RAMPDQDT : RATE OF INCREASE IN QTOT USED IF IQCONST=2
160. * TAUSOIL : CH4 SOIL SINK LIFETIME (IPCC VALUE IS 160 YR)
1 * METHHIST : IF=0, OLD CH4 HISTORY; ELSE, NEW (IPCC SAR) HISTORY
0.164 * POWC0 : PARAMETER IN CH4 MODEL (CENTRAL VALUE) * 0.164
0.082 * DELPOWC : PARAMETER IN CH4 MODEL (UNCERTAINTY) * 0.082
0.50 * POWCO0 : PARAMETER IN CH4 MODEL (CENTRAL VALUE) * 0.50
0.25 * DELPOWCO : PARAMETER IN CH4 MODEL (UNCERTAINTY) * 0.25
-0.66 * POWNO0 : PARAMETER IN CH4 MODEL (CENTRAL VALUE) * -0.66
-0.33 * DELPOWNO : PARAMETER IN CH4 MODEL (UNCERTAINTY) * -0.33
-0.21 * POWVO0 : PARAMETER IN CH4 MODEL (CENTRAL VALUE) * -0.21
-0.105 * DELPOWVO : PARAMETER IN CH4 MODEL (UNCERTAINTY) * -0.105
0 * IQREAD : IF=1, QEXTRA ADDED TO ANTHRO DQ; =2 USES QEXTRA ALONE
1371.55 * QOFFSET : VALUE SUBTRACTED FROM ALL QEXTRA.DAT VALUES
0.175 * QFACTOR : FACTOR TO MULTIPLY INPUT FROM QEXTRA.DAT LESS OFFSET
100 * NVALUES : NUMBER OF DIVISIONS IN GSIC MODEL
1234567890 *
NOTES : ICO2READ SWITCH
=1, USES JUST CO2 FORCING AS DETERMINED BY CONCS FROM
CO2INPUT.DAT, SCALED BY CO2SCALE AS SPECIFIED NOW
IN STAG.CFG
=2, USES CO2 CONCS FROM CO2INPUT.DAT AND OTHER GASES
FROM GASEM.$$$
=3, USES CO2 CONCS FROM CO2INPUT.DAT AND AEROSOL
FROM GASEM.$$$. OTHER GASES IGNORED.
=4, USES CO2 CONCS FROM CO2INPUT.DAT AND AEROSOL
BY SCALING EFOSS FROM GASEM.$$$. OTHER GASES IGNORED.
: S90OZ=0.32 BASED ON TOTAL OF 0.40 LESS 0.08W/m**2 DUE TO CH4 BUILDUP
: NOUT OPTIONS, 1 GIVES NH & SH TEMPS
2 GIVES LAND & OCEAN TEMPS AND THEIR RATIO
3 GIVES TEMP DISEQUILM (TE-T) & DEEP OCEAN TEMP
4 AS 3, PLUS GIVES EQUIV CO2 INSTEAD OF DELTA-Q
5 GIVES ABBREVIATED OUTPUT
: CONVERSIONS : 1.0GtC = 3664.13TgCO2
1.0TgN = 1.57113TgN2O
1.0TgS = 1.99809TgSO2
: IF USING QEXTRA.DAT FOR HOYT & SCHATTEN SOLAR INPUT, USE
QOFFSET=1371.55 AND QFACTOR=0.175
C MAG.FOR
C
C Revision history:
C
C 991222 * THIS IS THE VERSION FOR UNDP
C 991221 * THIS VERSION OF THE CODE CAN GIVE ODD ECH4 VALUES, SINCE
C IT USES A SLIGHTLY UNREALISTIC VALUE OF THE 1990 CH4
C LIFETIME IN ORDER TO GET BACK TO THE CH4 CONCENTRATION
C PROJECTIONS USED IN CH. 6 OF THE IPCC SAR. THE BEST
C GUESS ECH4 IN 1990 IS 571.6, SO THE VALUES PRINTED OUT
C HAVE A CORRECTION FACTOR ADDED TO GET TO THIS VALUE.
C METHOFF ADDED TO FIRST LINE OF GAS.EM. IF THIS IS SET
C TO 1, THEN CO2 FROM CH4 OXIDATION IS TURNED OFF SO
C THAT STABILIZATION EMISSIONS SCENARIOS STABILIZE AT THE
C CORRECT LEVELS. IT IS ONLY NECESSARY TO HAVE A METHOFF
C VALUE (METHOFF=1) IN GAS.EM FILES THAT LEAD TO
C STABILIZATION. IF ONE WANTS TO INCLUDE CO2 FROM CH4,
C THEN THE METHOFF VALUE JUST HAS TO BE SET TO A VALUE
C NOT EQUAL TO 1 IN THE GAS.EM FILE.
C 991220 * CODE CORRECTED TO GET BACK TO CH. 6 METHANE RESULTS IN
C ORDER TO REPRODUCE SAR CH. 6/7 TEMP AND MSL RESULTS AS
C NEARLY AS POSSIBLE. DIFFERENCES NON-ZERO, BUT NOT
C DETECTIBLE ON PLOTS.
C 991130 * ERROR IN W SLOWDOWN CORRECTED. ONLY AFFECTED RESULTS
C IF THRESHOLD W VALUE REACHED. WITH 7degC AS THE
C THRESHOLD (AS IN SAR) THIS NEVER HAPPENNED.
C 990727 * QEXTRA FILES NOW ACCESSED THRU MAGEXTRA, PREV MAGUSER.
C NVALUES ALSO ACCESSED THRU MAGEXTRA.
C 990722 * SPECIFICATION OF CH4 MODEL PARAMETERS PUT INTO
C MAGEXTRA.CFG FOR CONVENIENCE.
C METHANE MODEL PARAMETERS MODIFIED TO IMPROVE FIT TO
C IPCC SAR CH. 2. THE CONCENTRATIONS BELOW ARE MIDYEAR
C VALUES BASED ON 1700ppbv AS THE 1990 VALUE.
C SCENARIO 2050 CONCS 2100 CONCS
C IS92A,B 2735 2802 2793 3444 3599 3616
C IS92C 2130 2225 2224 1982 2068 2069
C IS92D 2137 2232 2230 2041 2139 2146
C IS92E 2942 3025 3014 4080 4298 4291
C IS92F 2965 3049 3038 4460 4711 4669
C CH.6 NEW CH.2 CH.6 NEW CH.2
C IN TERMS OF RADIATIVE FORCING, THE 2100 ERRORS ARE
C (W/m**2) : A,B 0.0057
C C 0.0004
C D 0.0030
C E 0.0021
C F 0.0123
C NOTE THAT (E.G.) THE FORCING 'ERROR' BETWEEN CH. 6 AND
C CH. 2 IN THE SAR WAS, FOR IS92a, ABOUT 0.058W/m**2.
C (SEE BELOW FOR MORE ON HOW THIS AROSE.)
C 980217 * FORMAT STATEMENT 226 RETURNED TO ORIGINAL. ALL FILE NAMES
C SET TO LOWER CASE TO ACCORD WITH SCENGEN.
C 980210 * NVALUES REMOVED FROM MAGUSER.CFG AND SET IN MAIN PROGRAM.
C 980208 * REGIONAL SULPHATE AEROSOL PATTERN DRIVERS CORRECTED (OUTPUT
C IN LO*, MID*, HI* AND USRDRIVE.OUT). THE PROBLEM WITH
C CALCULATING THESE IS THAT THE RAW TEMPERATURES (TSO21,2,3)
C CALCULATED BY USING THE SO2 EMISSIONS IN REGIONS 1,2,3
C ALONE ARE INTERNALLY INCONSISTENT BECAUSE OF THE NONLINEAR
C RELATIONSHIP BETWEEN ESO2 AND INDIRECT AEROSOL FORCING.
C IN OTHER WORDS, IN GENERAL, TALL (TEMPS WITH FULL SO2
C EMISSIONS) MINUS TGHG (TEMPS WITH NO SO2 EMISSIONS CHANGES
C AFTER 1990) WILL NOT EQUAL TSO21+TSO22+TSO23. TO CORRECT
C FOR THIS, THE FIRST METHOD USED WAS TO CALCULATE REGION 3
C TEMPS BY SUBTRACTING REGION 1 PLUS REGION 2 RESULTS FROM
C THE GLOBAL EMISSIONS RESULTS. THIS AT LEAST GIVES THE
C CORRECT SUM, BUT CAN LEAD TO RIDICULOUS RESULTS FOR REGION
C 3; FOR EXAMPLE, IF REGION 3 EMISSIONS ARE ZERO, IT WILL
C GIVE NONZERO TEMPS. THIS METHOD HAS THE 'ADVANTAGE' OF
C REQUIRING ONLY 16 NSIM LOOPS, SINCE REGION 3 RESULTS ARE
C NOT DIRECTLY CALCULATED. THE SECOND METHOD TRIED WAS TO
C USE THE REGIONAL BREAKDOWN OF FORCING TO SCALE REGIONAL
C TEMPERATURE RESULTS, AS CALCULATED BY THE MODEL. (THIS
C REQUIRES 4 MORE NSIM LOOPS TO COVER REGION 3.) WHILE THIS
C METHOD AVOIDS GIVING RIDICULOUS RESULTS, IT STILL REQUIRES
C ADDITIONAL CALCULATIONS TO GET THE FORCING BREAKDOWN. THE
C THIRD (AND BEST) METHOD IS TO CALCULATE THE FULL SET OF
C REGIONAL TEMPS (REQUIRING 20 NSIM LOOPS) AND THEN SCALING
C THE TSO2i BY (TALL-TGHG)/(SUM TSO2i) (GIVING XSO2i AS THE
C PATTERN DRIVER WEIGHTS). THIS IS THE METHOD NOW
C IMPLEMENTED. FOR DIAGNOSTIC PURPOSES, THE UNSCALED WEIGHTS
C ARE OUTPUT TO L0*, MID*, HI* AND USRDRIVE.RAW. THE WEIGHTS
C TO USE IN SCENGEN ARE OUTPUT TO LO*, MID*, HI* AND
C USRDRIVE.OUT. ALL THIS INVOLVED SUBSTANTIAL CHANGES TO THE
C CODE, TOO NUMEROUS TO ITEMIZE.
C 980120 * OUTPUT OF FORCING DETAILS (DECADAL BREAKDOWN AND LAND/OCEAN
C SPLIT) MOVED OUT OF MAG.OUT TO QDETAILS.OUT.
C 980118 * 'common /SecOrd/iSecYear,SecCH4' removed : no longer used.
C dcdt90 and emiss removed from 'common /TauMeth/'
C 980110 * METHANE SOIL LIFETIME SPECIFICATION MOVED TO MAGEXTRA.CFG
C INSTANTANEOUS METHANE LIFETIME ADDED TO MAG.OUT OUTPUT.
C NEW 1980-1990 CH4 HISTORY ADDED (ACCESSED BY SETTING
C METHHIST>0 IN MAGEXTRA.CFG) IN WHICH MID-1990 VALUES
C ARE C = 1700 ppbv, dC/dt = 10 ppbv/yr (AS IN IPCC SAR).
C 971227 * CONCENTRATION OUTPUT SECTION IN MAG.OUT CHANGED TO GIVE
C MIDYEAR INSTEAD OF END OF YEAR CONCS.
C 971214 * METHANE PARAMS CHANGED TO GIVE BETTER SIMULATIONS OF IPCC
C SAR RESULTS AS PER TABLE 2.5a IN SAR. 1990 TAU SET TO 9.1:
C PREVIOUSLY 9.08. IF DELTAU IS SET AT +/- 0.8 (CF. IPCC
C RANGE IS +/- 1.8), WHICH MAKES HIGH TAU = 9.9, AND
C ITAUMETH IS SET AT 3 (WHICH GIVES UPPER BOUND CONCS FOR
C THE ORIGINAL OSBORN & WIGLEY MODEL) THEN OUTPUT AGREES
C EXTREMELY WELL PRATHER'S IPCC RESULTS. THIS IMPLIES THAT
C PRATHER'S MODEL MUST HAVE A HIGH TAU AND/OR A TAU
C SENSITIVITY THAT DIFFERS SUBSTANTIALLY FROM THE CENTRAL
C OSBORN & WIGLEY VALUE.
C
C COMPARISON OF PRATHER IPCC (CH.2) AND CURRENT MODEL (IPCC
C VALUES HAVE 17ppbv ADDED BECAUSE THE IPCC VALUE IN 1990
C IS 17ppbv BELOW THE VALUE USED HERE - THIS MAKES THE
C COMPARISON HERE MORE TRUE TO A FORCING COMPARISON).
C
C RESULTS BELOW ARE IN ORDER : OLD VERSION (I.E., NUMBER
C USED IN CH.6 OF SAR); NEW VERSION; IPCC (CH.2). NOTE THAT
C THE MAIN REASON WHY THE CH.6 VALUES ARE SO DIFFERENT FROM
C THE IPCC CH.2 VALUES IS THAT THE PRESENT MODEL WAS TUNED
C TO AN EARLY VERSION OF PRATHER'S RESULTS. (THESE RESULTS
C WERE REVISED ON 12/27/97 BY USING ACTUAL MIDYEAR MODEL
C OUTPUT. PREVIOUSLY HAD INTERPOLATED FROM DECADAL DATA.)
C
C SCENARIO 2050 CONCS 2100 CONCS
C IS92A,B 2752 2821 2810 3461 3603 3633
C IS92C 2147 2249 2241 1999 2092 2086
C IS92D 2154 2256 2247 2058 2163 2163
C IS92E 2959 3041 3031 4097 4289 4308
C IS92F 2982 3066 3055 4477 4694 4686
C CH.6 NEW CH.2 CH.6 NEW CH.2
C
C THE TABLE BELOW COMPARES NEW WITH OLD (IPCC CH.6) T & MSL
C CHANGES OVER 1990-2100.
C
C SCENARIO SIMULATION OLD DT NEW DT OLD DMSL NEW DMSL
C IS92C LOW 0.8622 0.8747 12.63 12.86
C IS92A MID 2.0308 2.0504 48.94 49.28
C IS92E HIGH 3.5260 3.5601 94.75 95.21
C
C 971213 * ERROR IN CORRECTING INPUT METHANE EMISSIONS FOR 1990
C IMBALANCE CORRECTED (SUBSTANTIAL REVISION). THIS WAS FIRST
C CORRECTED IN OPT.FOR ON 961024. IT MAKES NO DIFFERENCE
C TO THE OUTPUT: OUTPUT IS ONLY AFFECTED FOR NON-ZERO DELTAU.
C IN PREVIOUS VERSION OF THE CODE, DELTAU WAS SET TO ZERO.
C 971123 * TWO ERRORS FOUND
C (1) S90IND=0.0 INPUT CAUSES CRASH. CORRECTED TEMPORARILY
C BY RESETTING S90IND=-0.0001
C (2) EQUIVCO2 MISCALCULATED (USED IF NOUT=4 ONLY).
C CODE CORRECTED.
C 970728 * NOTES TO USER :
C EMISSIONS INPUT IS VIA GAS.EM. THIS USED TO BE CALLED
C GASEM.$$$.
C OUTPUT FILES QFORCE, LO* MID* HI* & USRDRIVE NOT FULLY
C IMPLEMENTED IN THIS VERSION. SHOULD USE ISCENGEN=0
C ONLY IN MAGUSER.CFG.
C DATE SUBROUTINE getdat MAY NOT WORK ON ALL MACHINES.
C TO GET AN IDEA HOW THIS PROGRAM WORKS READ THE COMMENTS
C IN THIS FILE AND IN THE .CFG FILES.
C 970721 * ERROR IN TEMP PRINTOUT CORRECTED (AFFECTED NOUT=1 ONLY)
C 960420 * ICOLD OPTION REMOVED FROM CODE
C 960420 * NSIM LOOP CHANGED SO THAT FORCING ONLY CALCULATED
C ON FIRST PASS THROUGH LOOP. FOR FULL SCENGEN RUN
C (NUMSIMS=16), THIS CUTS RUN TIME BY 12%
C 960415 * CHECKED FOR CONSISTENCY WITH STAG.FOR
C NUMSIMS (LOOP FOR FULL MODEL RUNS) MODIFIED TO COVER
C ADDITIONAL AEROSOL RUNS FOR SCENGEN INPUT
C OUTPUT FILES ADDED FOR DATA TRANSFER TO SCENGEN
C 960308 * RUN TIME FOR VARIABLE W REDUCED BY 55% BY PUTTING
C INITIAL TEMP PROFILE CALCULATION IN SUBROUTINE INIT.
C OPTION ADDED TO USE OBSERVED INITIAL TEMP PROFILE
C INSTEAD OF THEORETICAL (EXPONENTIAL) PROFILE.
C 960226 * CO2 HISTORY REMOVED FROM BLOCK DATA TO CO2HIST.IN
C 951205 * OPTION ADDED TO CHANGE W(t) AT DIFF RATES IN NH VS SH
C 951109 * CO2SCALE CHANGED TO APPLY TO ICO2READ = 1, 3 & 4
C 951108 * EXTRA Q INPUT FILE RENAMED QEXTRA.IN
C 951005 * METHANE TEMPERATURE FEEDBACK REVISED BUT NOT ACTIVATED
C 950911 * SEPARATED FROM STAG.FOR FOR NEW MAGICC
C 950814 * MINOR CHANGE MADE TO QCH4OZ TO ENSURE THAT MID 1990
C VALUE IS EXACTLY 0.08. THIS IS NECESSARY IF TOTAL
C O3 FORCING IS TO BE 0.40 IN MID 1990.
C BUG IN 1990 S90OZ VALUE FINALLY TRACKED DOWN!
C 950811 * MINOR CHANGES MADE RE OZONE FORCING AROUND 1990. OUTPUT
C DOESN'T QUITE GIVE S90OZ VALUE IN 1990, BUT CORRECT
C VALUES ARE USED INTERNALLY. BUG STILL NOT QUITE FIXED.
C 950810 * OPTION ADDED TO USE OLD (INPUT=0) OR NEW GASEM.$$$ FILES
C CH4-INDUCED OZONE TERM ADDED AS OUTPUT TO DELQ BREAKDOWN
C 950808 * REVISED IPCC CONCENTRATION HISTORY OF JUNE 1995 ADDED.
C NEW BCO2 FORMULA INSERTED TO ACCOUNT FOR REVISED CONC
C HISTORY.
C CONVOLUTION CONSTANTS UPDATED.
C PSI FOR OCEAN FLUX UPDATED.
C DN(1990) UPDATED.
C CORRECTION MADE TO ALGORITHM FOR ADDING IN METHANE
C OXIDATION CONTIBUTION TO FOSSIL CO2 EMISSIOMS.
C 950804 * GAS-BY-GAS DQ BREAKDOWN CHANGED SO THAT TROP O3 ASSOC
C WITH METHANE IS INCLUDED WITH TROP O3 RATHER THAN CH4.
C 950719 * STAG.CFG INPUT ORDER CHANGED AND RATIONALIZED
C 950717 * QEXTRA INPUT EXPANDED TO ALLOW SEPARATE INPUT FROM ALL
C BOXES (NHO, NHL, SHO AND SHL)
C 950715 * CUMULATIVE ROUNDOFF TIME ERROR CORRECTION ALGORITHM
C IMPROVED
C UNNECESSARY ITEMS RELATED TO STAGOPT.FOR REMOVED
C QEXTRA INPUT ADDED (SET BY IQREAD, DATA FROM QEXTRA.DAT)
C METHOD TO GET SPATIAL FORCING BREAKDOWN ALTERED
C QGLOBE (=MIDYEAR DELTA-Q) ARRAY ADDED
C CALC OF TEQU MOVED TO RUNMOD SUBROUTINE
C 950612 * CARBON CYCLE BCO2, CPART AND AA UPDATED TO ACCOUNT FOR
C NEW INDUSTRIAL EMISSIONS HISTORY. BCO2 FIT SIMPLIFIED.
C 950610 * OPTIONS ADDED TO ADD SINGLE YEAR PULSE OR STEP EMISSIONS
C FOR CO2, CH4, N2O OR SO2. OUTPUT IN STAGPULS.OUT
C OPTION TO CHANGE DELQ(2xCO2) ADDED TO STAG.CFG
C 950609 * ICOMP AND IQCONST ALGORITHMS CORRECTED
C 950529 * INPUT FROM GASEM.$$$ CHANGED TO ALLOW ANY NUMBER OF
C KEY YEARS
C ICO2READ SWITCH MODIFIED :
C =1, USES JUST CO2 FORCING AS DETERMINED BY CONCS FROM
C CO2INPUT.DAT, SCALED BY CO2SCALE AS SPECIFIED NOW
C IN STAG.CFG
C =2, USES CO2 CONCS FROM CO2INPUT.DAT AND OTHER GASES
C FROM GASEM.$$$
C =3, USES CO2 CONCS FROM CO2INPUT.DAT AND AEROSOL
C FROM GASEM.$$$. OTHER GASES IGNORED.
C =4, USES CO2 CONCS FROM CO2INPUT.DAT AND AEROSOL
C BY SCALING EFOSS FROM GASEM.$$$. OTHER GASES IGNORED.
C 950410 * NOUT OPTIONS INCREASED
C 950402 * CARBON CYCLE MODEL PARAMS UPDATED TO IPCC95/6
C (1980S MASS INCREASE = 3.28GtC/yr)
C 950329 * CH4 AND N2O INPUT METHOD CHANGED. CORRECTION APPLIED TO
C ALL VALUES EQUAL TO ERROR IN 1990 BASED ON 1990 VALUE
C BEING CONSISTENT WITH TAU, C AND DC/DT IN 1990
C 950316 * CARBON CYCLE MODEL PARAMS UPDATED TO IPCC94/5
C (1980S MASS INCREASE = 3.24GtC/yr)
c 950315 * CH4 SOIL SINK TAU CHANGED TO 160YRS, FROM PRATHER.
c CORRECTION MADE TO USE TAUINIT CORRECTLY.
C 950315 * OZONE HISTORY CHANGED
C 950312 * NH/SH AND LAND OCEAN SPLIT REVISED
C 950311 * TROP O3 ASSOCIATED WITH CH4 CHANGES ADDED
C CO2 CONC HISTORY CHANGED TO FEB 95 VERSION OF IPCC
c 950309 * TROP O3 SEPARATED FROM AEROSOLS
C QBIO PUT INTO ARRAY
C RAD FORCING OUTPUTS REVISED
C GREENLAND PARAMS RESTORED TO IPCC90 VALUES
C GSIC MODEL PARAMS ALTERED FOR MANYTAU CASE
C ANTARCTIC DB3 CHANGED TO -0.04(0.01)0.06
C 950306 * GLACIER AND SMALL ICE CAP (GSIC) MODEL COMPLETELY REVISED
c 950224 * NEW AEROSOL VALUES AND INPUT METHOD ADDED.
c 950215 * NEW SMALL GLACIER & ICE SHEET MODELS INSERTED.
C 950214 * DIFFERENTIAL CLIMATE SENSITIVITY ADDED.
C VARIABLE UPWELLING RATE ADDED.
C OPTION FOR STEP OR RAMP-TO-CONSTANT DELQ ADDED.
C 941122 * SMALL GLACIER OPTION ADDED (ICEOPT)
C ICEOPT=1, ORIGINAL PARAMETERS
C ICEOPT=2, REVISED Z0 & TAU, ORIG OBS DEL-MSL
C ICEOPT=3, AS 2 WITH OBS DEL-MSL HALVED
C 941122 * GREENLAND AND ANTARCTIC MELT SENSITIVITIES CHANGED
C PER WARRICK EMAIL (G, 0.25+/-0.15, A, -0.30+/-0.15)
C 941118 * OPTIONS ADDED TO REMOVE GHG FORCING (IGHG=0
C IN STAG.CFG) AND TO REMOVE ALL GHG FORCING EXCEPT
C CO2 (IGHG=1 IN STAG.CFG). NOTE THAT GHG FORCING
C HERE INCLUDES BIOMASS AEROSOL FORCING AND TROP O3.
C 941016 * VARIOUS CHANGES MADE TO ACCORD WITH IPCC94. AEROSOL
C FORCING CHANGED TO A BEST GUESS DIRECT VALUE IN 1990
C OF 0.38W/m**2. INDIRECT FORCING MADE EQUAL TO DIRECT.
C PIECEWISE LINEAR BIOMASS BURNING AEROSOL TERM ADDED
C RISING TO -0.16W/m**2 IN 1990, ROUGHLY TRACKING
C TROPICAL GROSS DEFOR. TROP OZONE ADDED IN SULPHATE
C SUBROUTINE SO AS TO BE SPLIT NH/SH AS FOR AEROSOLS.
C FORCING TRACKS SO2 EMISSIONS TO +0.3W/m**2 IN 1990.
C NH/SH SPLIT FOR AEROSOLS CHANGED FROM 9/1 TO 7/1.
c 940918 * DN80S CHANGED TO 0.4(1.1)1.8
C 940715 * N2O HISTORY CHANGED TO ACCORD MORE FULLY WITH IPCC94
c 940715 * LIFETIMES FOR F11, F12, F22, F134A, N2O CHANGED TO
C MATCH IPCC94
C 940713 * MINOR CHANGE TO FORMAT STATEMENT 201
C 940701 * INIT N2O CHANGED TO 270PPBV
C 940617 * EXTRA CALL TO METHANE SUBROUTINE ADDED
C 940616 * METHANE SUBROUTINE CORRECTED
C 940517 * SWITCH ADDED TO ADD IN RANDOM FORCING
C 940508 * NEW (TELLUS) CARBON CYCLE MODEL INSERTED
C 940401 * CFC SUBROUTINE SIMPLIFIED. PRODUCTION INPUT ELIMINATED
C 940401 * EXPANSION ALGORITHM CHANGED TO INCREMENT EXPN AT EACH
C TIME STEP. OLD ALGORITHM RETAINED AND ACCESSIBLE VIA
C STAG.CFG. (THIS HAS NEGLIGIBLE AFFECT ON OUTPUT)
C 940331 * OPTION ADDED TO PRINT OUT LAND/OCEAN OR TOP/BOTTOM
C TEMPS INSTEAD OF NH/SH (NOUT). NOUT=3 ALSO GIVES TEMP
C DISEQUILIBRIUM. NOUT=4 IS AS NOUT=3, BUT GIVES EQUIV
C CO2 INSTEAD OF DELTA-Q
C 940328 * COLD START ALGORITHM CORRECTED TO ACCOUNT FOR AEROSOLS
C 940324 * ICO2READ SWITCH (940117) GENERALIZED TO READ CO2 CONCS
C IF ICO2READ=1 OR 2, AND CALCULATE QTOT BY SCALING UP
C QCO2 OF ICO2READ=1 OR BY ADDING OTHER GASES FROM
C GASEM.$$$ IF ICO2READ=2.
C 940322 * ICOMP SWITCH GENERALIZED. ICOLD=1 SETS DELQ=0.0 TO MID
C 1990. ICOMP=1 OVERWRITES POST-1990 FORCING WITH
C COMPOUNDED CO2 INCREASE AT COMPRATE.
C 940321 * NEW SWITCH ADDED TO DO COMPOUNDED CO2 AND COLD START.
C ICOMP=1 GIVES OBS DELQ TO MID 1990, THEN COMPOUND CO2
C AT COMPRATE. ICOMP=2 DITTO BUT DELQ=0.0 TO MID 1990.
C 940318 * SWITCH ADDED TO DO 1% COMPOUND CO2 INCREASE POST-1990
C ACCESSED BY PUTTING ICO2READ=9. ZERO FORCING TO 1990.
C 940317 * CORRECTION MADE TO 940117 CHANGE
C 940316 * STAG.CFG CHANGED TO ALLOW SPECIFICN OF PI, W, K & H.
C 940117 * N2O HISTORY CHANGED AGAIN
C 940117 * SWITCH ADDED TO STAG.CFG TO ALLOW CO2 CONCS TO BE READ
C FROM CO2INPUT.DAT TO OVER-RIDE CALCULATED CO2 CONCS.
C OTHER GAS FORCINGS ARE ACCOUNTED FOR BY SCALING CO2
C FORCING BY FACTOR SPECIFIED IN CO2INPUT.DAT
c 940104 * ISULPH INPUT IN STAG.CFG CHANGED TO DIRECT INPUT OF
C FACTOR THAT MULTIPLIES BEST GUESS FORCING VALUE (XSULF)
C 940103 * MORE GENERAL LOOP FOR DELTA-T2X PUT INTO STAG.CFG
C 931230 * IPCC CO2 CONC HISTORY INSERTED AND N2O HISTORY
C CHANGED TO HAVE INIT VALUE OF 265 INSTEAD OF 285.
C 931229 * SPLIT FROM STAG92B.FOR. RENAMED STAG.FOR, CFG FILE
C RENAMED STAG.CFG AND OUTPUT RENAMED STAG.OUT
C 931022 * LAST YEAR (LASTYEAR) FOR MODEL RUN ADDED TO CFG FILE
C 930120 * MORE MINOR CHANGES FOR CONSISTENCY WITH MAGMOD.FOR
C 930119 * CONST CH4 TAU OPTION ADDED FOR CONSISTENCY WITH MAGMOD.FOR
C (VIA STAG92B.CFG THRU ITAUMETH=4, AND CHOSEN TAU90CH4)
C 930117 * HALOCARBONS CHANGED TO ACCOUNT FOR ISAKSEN ET AL REPORT
C ON CFC14 AND CFC116, AND REVISED 1990 CONCS FOR CFC13
C AND METH/ENE CHLORIDE. CHANGE IN HISTORY TO GET REVISED
C 1990 EQUIV HFC134a CONC. CHANGE IN SCALING FACTORS TOO.
C 930104 * CH4 & N20 DELTA-Q REVISED TO ALLOW VARIABLE INIT CONCS.
C 921221 * TIME STEP PUT INTO STAG92B.CFG
C 921220 * PRINTOUT (ESCOUT.DAT) MODIFIED AND IMPROVED
C 921219 * HALOCARBON SCALING FACTORS REVISED AGAIN
C * HFC134a CONC REPLACED BY EQUIV HFC134a CONC
C 921216 * HALO COMMON BLOCK ADDED TO MAIN AND DELTAQ
C * IO3FEED AND ISULPH SPECIFICATION PUT INTO STAG92B.CFG
C 921215 * HALOCARBON SCALING FACTORS REVISED
C 921126 * ERROR IN SULPHATE FORCING PRIOR TO 1990 CORRECTED BY
C CHANGING HISTORY (ONLY OCCURRED FOR NSIM=1)
C * ERRORS IN FORMULAE FOR INITIAL OCEAN TEMP PROFILE CORRECTED
C * PRINTOUT OF KYRREF ADDED
C * IY0, IPRT AND KYRREF ADDED TO STAG92.CFG
c 921020 * date/name line
c 920918 * changes to output
C 920916 * .CFG file added
C * changes to carbon cycle INITs
c * OUTPUT changed
c * block data added for RS/6000 compatability
c 920915 * split from STAG92A; new CARBON +init+commons
c ------
c 920915 * minor correction to SULPHATE
c 920911 * revised HISTORY & CFC scaling factors; temps back to 1990;
c changed dcdt90
c 920910 * DT changed from 1.0 to 0.1
c * slight changes to output text
c * XKNS from 1 to 3
c * temps from 1765
c * array errors fixedin forcing
c 920909 * split from STAGGER; heavily cut down
C -----
C
c-------------------------------------------------------------------------------
c Written by Tom Wigley, Sarah Raper & Mike Salmon, Climatic Research Unit, UEA.
c-------------------------------------------------------------------------------
c
PROGRAM CLIMAT
C
C THIS IS THE CLIMATE MODEL MODULE.
C
parameter (iTp =740)
C
common /Limits/KEND
C
INTEGER IY1(100)
DIMENSION DCH4(100),DN2O(100)
DIMENSION DSO2(100),DSO21(100),DSO22(100),DSO23(100)
DIMENSION FOS(100),DEF(100)
DIMENSION QSO2SAVE(0:iTp+1),QDIRSAVE(0:iTp+1)
DIMENSION QRATIO(3,227:iTp+1)
C
DIMENSION TEMUSER(iTp),TEMLO(iTp),TEMMID(iTp),TEMHI(iTp),
&TEMNOSO2(iTp)
DIMENSION SLUSER(iTp),SLLO(iTp),SLMID(iTp),SLHI(iTp),
&SLNOSO2(iTp)
DIMENSION TALL(4,iTp-225),TGHG(4,iTp-225),TSO21(4,iTp-225),
&TSO22(4,iTp-225),TSO23(4,iTp-225),TREF(4)
DIMENSION XSO21(4,iTp-225),XSO22(4,iTp-225),XSO23(4,iTp-225)
C
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
C NOTE THAT EMISSIONS AND DELAY BOX ARRAYS START WITH J=226, =1990.
C
COMMON/CONCS/CH4(iTp),CN2O(iTp),C11(iTp),C12(iTp),C22(iTp),
&C134EQ(iTp),ECH4(226:iTp),EN2O(226:iTp),ECO(226:iTp),
&EVOC(226:iTp),ENOx(226:iTp),ESO2(0:iTp+1),ESO2SUM(226:iTp+1),
&ESO21(226:iTp+1),ESO22(226:iTp+1),ESO23(226:iTp+1)
C
COMMON/COBS/COBS(0:226)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/TANDSL/TEQU(iTp),TGAV(iTp),TNHO(iTp),
&TSHO(iTp),TNHL(iTp),TSHL(iTp),TDEEP(iTp),TNHAV(iTp),TSHAV(iTp),
&TLAND(iTp),TOCEAN(iTp),TOCN(40),TOCNPREV(40),IEXPAN,
&SIP,SGP,SAP,SLI(iTp),SLG(iTp),SLA(iTp),EX(0:iTp),SLT(iTp),
&QTOT(0:iTp),QGH(0:iTp),QOZ(0:iTp),QBIO(0:iTp),
&QSO2(0:iTp+1),QDIR(0:iTp+1)
C
COMMON/CAR/EL1,EL2,EL3,TINV0(5),TINV(4,5),A(3,5),AA(4,5),
&BCO2(4),BTGPP,BTRESP,BTHUM,GAMP,GPP0,RESP0,QA0,U0,C0,B340(4),
&PHI,RG,TAUP,TAUH,TAUS,THP,THS,THH0,THS0,THPL,G1,G2,G3,FACTOR,
&EL21,EL32,XX1,XX2,XX3,XX4,XX5,XX6,DEE1,DEE2,DEE3,DEE4,DEE5,DEE6,
&FL1,FL2,FL3,XL,GAMH,GAMS,QS0,BTSOIL,FERTTYPE,TOTEM,CONVTERP,
&R(4),CPART(4,5),DELMASS(4,226:iTp),ABFRAC(4,226:iTp)
C
common /meth/emeth(226:iTp),imeth,ch4l(225:iTp),ch4b(225:iTp),
+ch4h(225:iTp),ef4(226:iTp)
c
common /radforc/qco2(0:iTp),qm(0:iTp),qn(0:iTp),qcfc(0:iTp)
COMMON /QOZONE/QCH4O3(0:iTp)
c
common /StratH2O/StratH2O
c
common /TauMeth/ITAUMETH,TAU90CH4,TauSoil,DELTAU,
+ TAUINIT,ch4bar90,eno90,eco90,evo90
COMMON /TCH4/TCH4(iTp)
COMMON /CH4HIST/METHHIST
common /powc/powc0,delpowc,powco0,delpowco
common /pown/powno0,delpowno,powvo0,delpowvo
c
common /TauNitr/TAUN2O
c
common /O3feedback/iO3feed
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,enat,es1990,
&dqdcoz,QCH4OZ
COMMON /IGHG/IGHG
COMMON /CO2READ/ICO2READ,XC(226:iTp),CO2SCALE
COMMON /COLDETC/ICOMP,COMPRATE,qtot86
COMMON /QCON/IQCONST,QCONST,RAMPDQDT
COMMON /DSENS/IXLAM,XLAML,XLAMO
common /netdef/ednet(226:iTp),DUSER,FUSER
COMMON /VARW/Z(40),W(2),DW(2),TO0(2),TP0(2),WNH(iTp),WSH(iTp),
&TW0NH,TW0SH,IVARW
COMMON /QSPLIT/QNHO,QNHL,QSHO,QSHL,QGLOBE(0:iTp),
&QQNHO(0:iTp),QQNHL(0:iTp),QQSHO(0:iTp),QQSHL(0:iTp),
&QQQNHO(0:iTp),QQQNHL(0:iTp),QQQSHO(0:iTp),QQQSHL(0:iTp)
C
COMMON /ICE/TAUI,DTSTAR,DTEND,BETAG,BETA1,BETA2,DB3
COMMON /ICEPARA/ZI0,NVALUES,MANYTAU,TAUFACT,ZI0BIT(100)
COMMON /ICEINIT/SIPBIT(100),SIBIT(100)
COMMON /AREAS/FNO,FNL,FSO,FSL
COMMON /TEMEXP/TEMEXP(2,40)
C
COMMON /QADD/IQREAD,JQFIRST,JQLAST,QEX(0:iTp),QEXNH(0:iTp),
&QEXSH(0:iTp),QEXNHO(0:iTp),QEXNHL(0:iTp),QEXSHO(0:iTp),
&QEXSHL(0:iTp)
C
COMMON /OLDTZ/IOLDTZ
COMMON /NSIM/NSIM,NCLIM,ISCENGEN
COMMON /CH4CORR/CORRUSER,CORRMHI,CORRMMID,CORRMLO
C
C ********************************************************************
C
character*20 mnem
character*3 month(12)
data (month(i),i=1,12) /'Jan','Feb','Mar','Apr','May','Jun',
+ 'Jul','Aug','Sep','Oct','Nov','Dec'/
C
C ********************************************************************
C
C READ CO2 CONCENTRATION HISTORY
C
lun = 42 ! spare logical unit no.
open(unit=lun,file='co2hist.IN',status='OLD')
DO ICO2=0,226
READ(LUN,4444)COBS(ICO2)
END DO
CLOSE(lun)
C
C ********************************************************************
C
C READ PARAMETERS FROM MAGUSER.CFG
C
lun = 42 ! spare logical unit no.
open(unit=lun,file='maguser.cfg',status='OLD')
READ(LUN,4240) ISCENGEN
READ(LUN,4240) IUSERGAS
READ(LUN,4241) DT2XUSER
READ(LUN,4241) RLO
READ(LUN,4241) HM
READ(LUN,4241) YK
READ(LUN,4241) W0
READ(LUN,4241) PI
READ(LUN,4240) IVARW
READ(LUN,4241) TW0NH
READ(LUN,4241) TW0SH
READ(LUN,4241) ZI0IN
READ(LUN,4241) TAUIN
READ(LUN,4241) DTSTARIN
READ(LUN,4241) DTENDIN
READ(LUN,4241) BETAGIN
READ(LUN,4241) BETA1IN
READ(LUN,4241) BETA2IN
READ(LUN,4241) DB3IN
READ(LUN,4240) KYRREF
READ(LUN,4240) IY0
READ(LUN,4240) LASTYEAR
READ(LUN,4240) ITEMPRT
READ(LUN,4240) ITAUMETH
READ(LUN,4241) TAU90CH4
READ(LUN,4241) TAUINIT
READ(LUN,4241) DELTAU
READ(LUN,4240) IMETH
READ(LUN,4240) IO3FEED
READ(LUN,4241) S90DIR
READ(LUN,4241) S90IND
READ(LUN,4241) S90BIO
READ(LUN,4241) DUSER
READ(LUN,4241) FUSER
READ(LUN,4241) TAUN2O
close(lun)
C
C INTERIM CORRECTION TO AVOID CRASH IF S90IND SET TO ZERO IN
C MAGUSER.CFG
C
IF(S90IND.EQ.0.0)S90IND=-0.0001
C
XK=YK*3155.76
C
C SELECT USER OR DEFAULT VALUES FOR GAS CYCLE PARAMETERS
C
IF(IUSERGAS.EQ.0)THEN ! DEFAULT
ITAUMETH = 2
TAUINIT = 8.30
TAU90CH4 = 9.08
DELTAU = 0.5
IMETH = 1
IO3FEED = 1
S90DIR = -0.3
S90IND = -0.8
S90BIO = -0.2
DUSER = 1.1
FUSER = 2.0
TAUN2O = 120.
ENDIF
C
C ************************************************************
C
C READ PARAMETERS FROM MAGEXTRA.CFG
C
lun = 42 ! spare logical unit no.
open(unit=lun,file='magextra.cfg',status='OLD')
READ(LUN,4240) IOLDTZ
READ(LUN,4241) DT
READ(LUN,4241) CO2DELQ
READ(LUN,4241) XKLO
READ(LUN,4241) XKNS
READ(LUN,4240) IEXPAN
READ(LUN,4240) NOUT
READ(LUN,4240) IEMPRT
READ(LUN,4240) ICO2PRT
READ(LUN,4240) ICONCPRT
READ(LUN,4240) IQDECPRT
READ(LUN,4240) INSLOPRT
READ(LUN,4240) IQGASPRT
READ(LUN,4240) IDIS
READ(LUN,4241) TOFFSET
READ(LUN,4241) STRATH2O
READ(LUN,4241) DQDCOZ
READ(LUN,4241) S90OZ
READ(LUN,4241) ENAT
READ(LUN,4240) IGHG
READ(LUN,4240) ICO2READ
READ(LUN,4241) CO2SCALE
READ(LUN,4241) SO2SCALE
READ(LUN,4240) ICOMP
READ(LUN,4241) COMPRATE
READ(LUN,4240) IQCONST
READ(LUN,4241) QCONST
READ(LUN,4241) RAMPDQDT
READ(LUN,4241) TAUSOIL
READ(LUN,4240) METHHIST
READ(LUN,4241) POWC0
READ(LUN,4241) DELPOWC
READ(LUN,4241) POWCO0
READ(LUN,4241) DELPOWCO
READ(LUN,4241) POWNO0
READ(LUN,4241) DELPOWNO
READ(LUN,4241) POWVO0
READ(LUN,4241) DELPOWVO
READ(LUN,4240) IQREAD
READ(LUN,4241) QOFFSET
READ(LUN,4241) QFACTOR
READ(LUN,4240) NVALUES
close(lun)
C
C TAU FOR CH4 SOIL SINK CHANGED TO ACCORD WITH IPCC94 (160 yr).
C SPECIFICATION OF TauSoil MOVED TO MAGEXTRA.CFG ON 1/10/97.
C
C ************************************************************
C
IF(ISCENGEN.EQ.1)THEN
NUMSIMS=20
ELSE
NUMSIMS=4
ENDIF
C
C SWITCH TO PRODUCE ONLY 1 SIMULATIONS.
C
IF(ISCENGEN.EQ.9)NUMSIMS=1
C
IXLAM=1
IF(RLO.EQ.1.0) IXLAM=0
C
C ************************************************************
C
IF(ICO2READ.GE.1)IMETH=0
C
IF(DT.GT.1.0)DT=1.0
C
QXX=CO2DELQ
Q2X=QXX*ALOG(2.)
C
IF(IQCONST.NE.0)THEN
S90DIR=0.0
S90IND=0.0
S90BIO=0.0
S90OZ =0.0
ENDIF
C
C *****************************************************************
C
C FO(I) AND FL(I) ARE N.H. AND S.H. OCEAN AND LAND FRACTIONS.
C
FO(1)=1.0-FL(1)
FO(2)=1.0-FL(2)
C
FK=RHO*SPECHT*HTCONS/31.5576
C
C ***************************************************************
C
C TRAP TO CATCH AND OVERWRITE UNREALISTIC D80SIN
C
IF(DUSER.LT.-0.5)THEN
DOLD=DUSER
DUSER=-0.5
WRITE(8,808)DOLD,DUSER
ENDIF
IF(DUSER.GT.3.0)THEN
DOLD=DUSER
DUSER=3.0
WRITE(8,809)DOLD,DUSER
ENDIF
C
C ************************************************************
C
C READ CO2 CONCS DIRECTLY FROM CO2INPUT.DAT IF ICO2READ.GE.1
C
IF(ICO2READ.GE.1)THEN
lun = 42 ! spare logical unit no.
open(unit=lun,file='co2input.dat',status='OLD')
C
C CO2INPUT.DAT MUST HAVE FIRST YEAR = 1990 AND MUST HAVE ANNUAL END
C OF YEAR VALUES. FIRST LINE OF FILE GIVES LAST YEAR OF ARRAY.
C
READ(lun,900)LCO2
ILCO2=LCO2-1764
DO JCO2=226,ILCO2
READ(lun,902)JYEAR,XC(JCO2)
END DO
C
C IF LAST YEAR OF INPUT CO2 DATA LESS THAN LASTYEAR FILL OUT
C ARRAY WITH CONSTANT CO2
C
IF(LASTYEAR.GT.LCO2)THEN
DO JCO2=ILCO2+1,LASTYEAR-1764
XC(JCO2)=XC(ILCO2)
END DO
ENDIF
close(lun)
ENDIF
C
C ************************************************************
C
C READ EXTRA FORCING IF IQREAD=1 OR 2. IF IQREAD=1, FORCING
C IN QEXTRA.IN IS ADDED TO ANTHROPOGENIC FORCING. IF IQREAD=2,
C QEXTRA.IN FORCING IS USED ALONE. QEXTRA.IN HAS A FLAG (NCOLS)
C TO TELL WHETHER THE DATA ARE GLOBAL (ONE Q COLUMN), HEMISPHERIC
C (TWO Q COLUMNS, NH THEN SH) OR FOR ALL BOXES (FOUR Q COLUMNS,
C IN ORDER NHO, NHL, SHO, SHL)
C
IF(IQREAD.GE.1)THEN
lun = 42 ! spare logical unit no.
open(unit=lun,file='qextra.in',status='OLD')
C
READ(LUN,900)NCOLS
READ(lun,901)IQFIRST,IQLAST
JQFIRST=IQFIRST-1764
C
C TRAP IN CASE FIRST YEAR IS BEFORE 1765
C
IF(JQFIRST.LT.1)THEN
DO JQ=JQFIRST,0
IF(NCOLS.EQ.1)READ(lun,902)JYEAR,QQQGL
IF(NCOLS.EQ.2)READ(lun,903)JYEAR,QQQNH,QQQSH
IF(NCOLS.EQ.4)READ(lun,904)JYEAR,QQQNHO,QQQNHL,QQQSHO,QQQSHL
END DO
JQFIRST=1
IQFIRST=1765
ENDIF
C
JQLAST=IQLAST-1764
DO JQ=JQFIRST,JQLAST
IF(NCOLS.EQ.1)THEN
READ(lun,902)JYEAR,QEX(JQ)
QEXNHO(JQ)=QEX(JQ)-QOFFSET
QEXNHL(JQ)=QEX(JQ)-QOFFSET
QEXSHO(JQ)=QEX(JQ)-QOFFSET
QEXSHL(JQ)=QEX(JQ)-QOFFSET
QEXNHO(JQ)=QFACTOR*QEXNHO(JQ)
QEXNHL(JQ)=QFACTOR*QEXNHL(JQ)
QEXSHO(JQ)=QFACTOR*QEXSHO(JQ)
QEXSHL(JQ)=QFACTOR*QEXSHL(JQ)
ENDIF
IF(NCOLS.EQ.2)THEN
READ(lun,903)JYEAR,QEXNH(JQ),QEXSH(JQ)
QEXNHO(JQ)=QEXNH(JQ)-QOFFSET
QEXNHL(JQ)=QEXNH(JQ)-QOFFSET
QEXSHO(JQ)=QEXSH(JQ)-QOFFSET
QEXSHL(JQ)=QEXSH(JQ)-QOFFSET
QEXNHO(JQ)=QFACTOR*QEXNHO(JQ)
QEXNHL(JQ)=QFACTOR*QEXNHL(JQ)
QEXSHO(JQ)=QFACTOR*QEXSHO(JQ)
QEXSHL(JQ)=QFACTOR*QEXSHL(JQ)
ENDIF
IF(NCOLS.EQ.4)THEN
READ(lun,904)JYEAR,QEXNHO(JQ),QEXNHL(JQ),QEXSHO(JQ),
& QEXSHL(JQ)
QEXNHO(JQ)=QEXNHO(JQ)-QOFFSET
QEXNHL(JQ)=QEXNHL(JQ)-QOFFSET
QEXSHO(JQ)=QEXSHO(JQ)-QOFFSET
QEXSHL(JQ)=QEXSHL(JQ)-QOFFSET
QEXNHO(JQ)=QFACTOR*QEXNHO(JQ)
QEXNHL(JQ)=QFACTOR*QEXNHL(JQ)
QEXSHO(JQ)=QFACTOR*QEXSHO(JQ)
QEXSHL(JQ)=QFACTOR*QEXSHL(JQ)
ENDIF
END DO
IF(NCOLS.EQ.1.OR.NCOLS.EQ.2)THEN
QEXNH(JQFIRST-1)=QEXNH(JQFIRST)
QEXSH(JQFIRST-1)=QEXSH(JQFIRST)
ENDIF
IF(NCOLS.EQ.4)THEN
QEXNHO(JQFIRST-1)=QEXNHO(JQFIRST)
QEXNHL(JQFIRST-1)=QEXNHL(JQFIRST)
QEXSHO(JQFIRST-1)=QEXSHO(JQFIRST)
QEXSHL(JQFIRST-1)=QEXSHL(JQFIRST)
ENDIF
close(lun)
ENDIF
C
C ******************************************************************
C
C Read in gas emissions from GAS.EM
C
lun = 42 ! spare logical unit no.
C
open(unit=lun,file='gas.em',status='OLD')
C
C READ HEADER AND NUMBER OR ROWS OF EMISIONS DATA FROM GASEM.$$$
C
read(lun,4244) NVAL,METHOFF
read(lun,'(a)') mnem
read(lun,*) ! skip description
read(lun,*) ! skip column headings
read(lun,*) ! skip units
C
C TURN OFF CO2 FROM CH4 OXIDATION IF METHOFF=1
C
IF(METHOFF.EQ.1)IMETH=0
C
C READ INPUT EMISSIONS DATA FROM GAS.EM
C NOTE THAT, IN CONTRAST TO STAG.FOR AND EARLIER VERSIONS OF MAGICC,
C THE SO2 EMISSIONS (BY REGION) MUST BE INPUT AS CHANGES FROM 1990.
C
C SET 1990 VALUE OF GLOBAL SO2 EMISSIONS
C
es1990=75.0
C
do i=1,NVAL
read(lun,4242) IY1(i),fos(i), def(i), DCH4(i), DN2O(i),
& DSO21(i),DSO22(i),DSO23(i)
C
C ADJUST SO2 EMISSIONS INPUT
C
DSO21(I)= DSO21(I)+es1990
DSO22(I)= DSO22(I)+es1990
DSO23(I)= DSO23(I)+es1990
DSO2(I) = DSO21(I)+DSO22(I)+DSO23(I)-2.0*es1990
C
end do
close(lun)
C
C ************************************************************
C
C TRAP TO CATCH INCONSISTENCY BETWEEN LASTYEAR FROM CFG FILE
C AND LAST YEAR OF EMISSIONS INPUT OR MAX ARRAY SIZE
C
IF(ICO2READ.EQ.0)THEN
IF((LASTYEAR-1764).GT.iTp)LASTYEAR=iTp+1764
IF(LASTYEAR.GT.IY1(NVAL)) LASTYEAR=IY1(NVAL)
ENDIF
IYEND = LASTYEAR
KEND = IYEND-1764
KREF = KYRREF-1764
TEND = FLOAT(KEND-1)
C
C Store 1990 values of CO, VOC and NOx emissions
C
eco90 = 499.0
evo90 = 321.0
eno90 = 55.0
C
C Offset IY1 entries from 1990 : i.e., make IY1(1)=0,
C IY1(2)=IY1(2)-1990, etc.
C
do i=1,NVAL
IY1(i) = IY1(i) - 1990
end do
C
C ***************************************************************
C
C INITIAL (1990) METHANE VALUES: EMISS (LO, MID, HI OR CON) IS THE
C 'CORRECT' 1990 VALUE CALCULATED TO BE CONSISTENT WITH THE
C CORRESPONDING VALUE OF THE 1990 CH4 TAU. IN GENERAL, EMISS
C WILL BE INCONSISTENT WITH THE 1990 INPUT VALUE. THIS IS CORRECTED
C BY OFFSETTING ALL INPUT VALUES BY THE 1990 'ERROR'. SINCE
C THIS ERROR DEPENDS ON TAU, DIFFERENT OFFSETS MUST BE CALCULATED FOR
C EACH 1990 TAU VALUE.
C Note that C and dC/dt must be for mid 1990. VALUES CORRECTED TO
C AGREE WITH IPCC SAR ON 1/10/98. HISTORY CORRECTED TOO. NOTE THAT
C THE OLD C AND dC/dt VALUES DON'T QUITE AGREE WITH THE OLD CONC
C HISTORY (WHICH GAVE C(MID1990)=1717 & dC/dt(MID1990)=10.83).
C
IF(METHHIST.EQ.0)THEN
CH4BAR90 = 1711.5
DCDT90 = 11.22
ELSE
CH4BAR90 = 1700.
DCDT90 = 10.
ENDIF
C
T90LO = TAUINIT-DELTAU
T90MID = TAUINIT
T90HI = TAUINIT+DELTAU
T90CON = TAU90CH4
C
EMISSLO = 2.75*(DCDT90+CH4BAR90*(1./T90LO +1./TAUSOIL))
EMISSMID = 2.75*(DCDT90+CH4BAR90*(1./T90MID+1./TAUSOIL))
EMISSHI = 2.75*(DCDT90+CH4BAR90*(1./T90HI +1./TAUSOIL))
EMISSCON = 2.75*(DCDT90+CH4BAR90*(1./T90CON+1./TAUSOIL))
C
C INITIAL (1990) N2O VALUE: FOLLOWS CH4 CASE, BUT THERE IS ONLY ONE
C CORRECTION FACTOR (emissN). THIS IS CALCULATED to be consistent
C with TAUN2O. Note that C and dC/dt must be for mid 1990.
C
N2Obar90 = 307.6
dcdt90N = 0.8
emissN = 4.81*(dcdt90N+N2Obar90/TAUN2O)
C
C ADD (OR SUBTRACT) CONSTANT TO ALL CH4 AND N2O EMISSIONS TO GIVE
C 1990 VALUE CONSISTENT WITH LIFETIME, CONC AND DC/DT.
C FOR CH4, ONLY THE CORRECTION FOR THE USER-SPECIFIED 1990 LIFETIME
C IS APPLIED (GIVEN BY THE CHOICE OF ITAUMETH).
C
C SPECIFY USER LIFETIME. NOTE, THIS IS RESPECIFIED IN DELTAQ.
C
IF(ITAUMETH.EQ.1) T90USER = T90LO
IF(ITAUMETH.EQ.2) T90USER = T90MID
IF(ITAUMETH.EQ.3) T90USER = T90HI
IF(ITAUMETH.EQ.4) T90USER = T90CON
C
CORRMLO = EMISSLO - DCH4(1)
CORRMMID = EMISSMID - DCH4(1)
CORRMHI = EMISSHI - DCH4(1)
CORRMCON = EMISSCON - DCH4(1)
CORRN2O = EMISSN - DN2O(1)
C
IF(ITAUMETH.EQ.1) CORRUSER = CORRMLO
IF(ITAUMETH.EQ.2) CORRUSER = CORRMMID
IF(ITAUMETH.EQ.3) CORRUSER = CORRMHI
IF(ITAUMETH.EQ.4) CORRUSER = CORRMCON
C
do i=1,NVAL
DCH4(I)=DCH4(I)+CORRUSER
DN2O(I)=DN2O(I)+CORRN2O
end do
C
C ***************************************************************
C
C INTERP(no-of-input-values,start-of-output,key-years,input-values,output)
C INTERP(N,ISTART,IY,X,Y)
C
C SUBROUTINE INTERP TAKES EMISSIONS (E.G. ARRAY X=DCO2) SPECIFIED AT
C KEY YEARS, WITH THE KEY YEARS IDENTIFIED RELATIVE TO 1990=0 IN
C THE IY ARRAY, AND LINEARLY INTERPOLATES TO FIND ANNUAL VALUES
C WHICH ARE OUTPUT TO AN ARRAY (E.G. Y=ECH4) BEGINNING WITH ISTART
C AND ENDING WITH IEND. THUS X(1) AND Y(ISTART) ARE THE SAME, AS
C ARE X(IY(N)) AND Y(IEND). NOTE THAT ONLY ISTART MUST BE SPECIFIED
C (USUALLY AS 226, CORRESP TO 1990), SINCE IEND=ISTART+IY(N).
C HERE, N IS THE SIZE OF THE IY ARRAY.
C
call interp(NVAL,226,IY1,fos,ef)
call interp(NVAL,226,IY1,def,ednet)
C
call interp(NVAL,226,IY1,DCH4,ECH4)
call interp(NVAL,226,IY1,DN2O,EN2O)
C
C NOTE, IF ESO2 WERE BEING INTERPOLATED, WOULD HAVE TO HAVE ESO2(226)
C AS LAST ARGUMENT BECAUSE OF MISMATCH OF ESO2 AND Y ARRAYS IN MAIN
C AND SUBROUTINE INTERP. THIS IS AVOIDED BY USING ESO2SUM.
C
call interp(NVAL,226,IY1,dSO2,ESO2SUM)
C
call interp(NVAL,226,IY1,dSO21,ESO21)
call interp(NVAL,226,IY1,dSO22,ESO22)
call interp(NVAL,226,IY1,dSO23,ESO23)
C
C SET ESO2 ARRAY
C
DO KE=226,KEND
ESO2(KE)=ESO2SUM(KE)
END DO
C
C ****************************************************************
C
CALL INIT
C
C ****************************************************************
C
DO KE=226,KEND
ECO(KE) = eco90
EVOC(KE) = evo90
ENOX(KE) = eno90
END DO
C
C IF ICO2READ=4, TOTAL SO2 EMISSIONS ARE REPLACED BY SCALED EFOSS.
C NUMSIMS MUST BE SET TO 4 SINCE SO2 EM ARE EFFECTIVELY IGNORED.
C
IF(ICO2READ.EQ.4)THEN
NUMSIMS=4
DO IS=226,KEND
ESO2SUM(IS)=ESO2SUM(226)+SO2SCALE*(EF(IS)-EF(226))
END DO
ENDIF
C
C LINEARLY EXTRAPOLATE LAST ESO2 VALUES FOR ONE YEAR
C
ESO2SUM(KEND+1) = 2.*ESO2SUM(KEND)-ESO2SUM(KEND-1)
ESO21(KEND+1) = 2.*ESO21(KEND)-ESO21(KEND-1)
ESO22(KEND+1) = 2.*ESO22(KEND)-ESO22(KEND-1)
ESO23(KEND+1) = 2.*ESO23(KEND)-ESO23(KEND-1)
ESO2(KEND+1) = ESO2SUM(KEND+1)
C
C CALCULATE REGIONAL SO2 EMISSIONS RATIOS FOR LATER APPORTIONING
C OF AEROSOL FORCING
C
DO KE=227,KEND+1
C
IF(ESO2SUM(KE).NE.es1990)THEN
QRATIO(1,KE)=(ESO21(KE)-es1990)/(ESO2SUM(KE)-es1990)
QRATIO(2,KE)=(ESO22(KE)-es1990)/(ESO2SUM(KE)-es1990)
QRATIO(3,KE)=(ESO23(KE)-es1990)/(ESO2SUM(KE)-es1990)
ELSE
QRATIO(1,KE)=0.0
QRATIO(2,KE)=0.0
QRATIO(3,KE)=0.0
ENDIF
C
END DO
C
C ************************************************************
C
C OPEN MAIN OUTPUT FILE (MAG.OUT) AND QDETAILS.OUT FILE
C
OPEN(UNIT=8,FILE='mag.out',STATUS='UNKNOWN')
OPEN(UNIT=88,FILE='qdetails.out',STATUS='UNKNOWN')
C
C ****************************************************************
C
call getdat(myr,imon,iday)
write(8,87) mnem,iday,month(imon),myr
write(88,87) mnem,iday,month(imon),myr
87 format(' Emissions profile: ',a20,20x,' Date: ',i2,1x,a3,1x,i4)
C
C ********************************************************************
C
C WRITE OUT HEADER INFORMATION FOR MAG.OUT
C
C SCALING FACTOR FOR CO2 FORCING :
C (QTOT-Q1990)=(QCO2-QCO2.1990)*CO2SCALE
C
SCAL=100.*(CO2SCALE-1.)
IF(ICO2READ.EQ.1)WRITE(8,871)SCAL
IF(ICO2READ.EQ.2)WRITE(8,872)
IF(ICO2READ.EQ.3)WRITE(8,873)
IF(ICO2READ.EQ.4)WRITE(8,874)SO2SCALE
C
IF(ICO2READ.GE.1.AND.ICOMP.EQ.1)WRITE(8,876)
IF(ICOMP.EQ.1)WRITE(8,877)COMPRATE
C
C ************************************************************
C
IF(IQREAD.EQ.0)WRITE(8,756)
IF(IQREAD.GE.1)THEN
IF(NCOLS.EQ.1)WRITE(8,757)IQFIRST,IQLAST
IF(NCOLS.EQ.2)WRITE(8,758)IQFIRST,IQLAST
IF(NCOLS.EQ.4)WRITE(8,759)IQFIRST,IQLAST
IF(QOFFSET.NE.0.0)WRITE(8,760)QOFFSET
IF(QFACTOR.NE.1.0)WRITE(8,761)QFACTOR
ENDIF
IF(IQREAD.EQ.2)WRITE(8,762)
C
WRITE (8,10) Q2X
C WRITE (8,11) FO(1),FO(2),FL(1),FL(2)
C
C ************************************************************
C
IF(IQREAD.EQ.0)WRITE(88,756)
IF(IQREAD.GE.1)THEN
IF(NCOLS.EQ.1)WRITE(88,757)IQFIRST,IQLAST
IF(NCOLS.EQ.2)WRITE(88,758)IQFIRST,IQLAST
IF(NCOLS.EQ.4)WRITE(88,759)IQFIRST,IQLAST
IF(QOFFSET.NE.0.0)WRITE(88,760)QOFFSET
IF(QFACTOR.NE.1.0)WRITE(88,761)QFACTOR
ENDIF
IF(IQREAD.EQ.2)WRITE(88,762)
C
WRITE (88,10) Q2X
C
C ************************************************************
C
IF(S90DIR.EQ.0.0) write(8,*) 'DIRECT AEROSOL FORCING IGNORED'
IF(ABS(S90DIR).GT.0.0) write(8,60)S90DIR
IF(S90IND.EQ.0.0) write(8,*) 'INDIRECT AEROSOL FORCING IGNORED'
IF(ABS(S90IND).GT.0.0) write(8,61)S90IND
IF(S90BIO.EQ.0.0) write(8,*) 'BIOMASS AEROSOL FORCING IGNORED'
IF(ABS(S90BIO).GT.0.0) write(8,62)S90BIO
IF(S90OZ.EQ.0.0) write(8,*) 'TROPOSPHERIC OZONE FORCING IGNORED'
IF(ABS(S90OZ).GT.0.0) write(8,63)S90OZ
C
if(iO3feed.eq.0) write(8,*)
& 'STRAT OZONE DEPLETION FEEDBACK OMITTED'
if(iO3feed.ne.0) write(8,*)
& 'STRAT OZONE DEPLETION FEEDBACK INCLUDED'
C
write(8,53)STRATH2O
C
C ****************************************************************
C ****************************************************************
C
C Run model NUMSIMS times for different values of DT2X.
C The input parameter ICEOPT determines what ice melt parameter
C values are used in each case.
C
C THE KEY FOR NSIM IS AS FOLLOWS (CASES 17-20 ADDED FEB 7, 1998) ...
C NOTE : IF ISCENGEN=9, ONLY NSIM=1 IS RUN, BUT NCLIM IS SET TO 4.
C
C NSIM CLIM MODEL EMISSIONS NESO2 NCLIM
C 1 LOW ALL 1 1
C 2 MID ALL 1 2
C 3 HIGH ALL 1 3
C 4 USER ALL 1 4
C 5 LOW ESO2 = CONST AFTER 1990 2 1
C 6 MID ESO2 = CONST AFTER 1990 2 2
C 7 HIGH ESO2 = CONST AFTER 1990 2 3
C 8 USER ESO2 = CONST AFTER 1990 2 4
C 9 LOW ESO2 = ESO2(REGION 1) 3 1
C 10 MID ESO2 = ESO2(REGION 1) 3 2
C 11 HIGH ESO2 = ESO2(REGION 1) 3 3
C 12 USER ESO2 = ESO2(REGION 1) 3 4
C 13 LOW ESO2 = ESO2(REGION 2) 4 1
C 14 MID ESO2 = ESO2(REGION 2) 4 2
C 15 HIGH ESO2 = ESO2(REGION 2) 4 3
C 16 USER ESO2 = ESO2(REGION 2) 4 4
C 17 LOW ESO2 = ESO2(REGION 3) 5 1
C 18 MID ESO2 = ESO2(REGION 3) 5 2
C 19 HIGH ESO2 = ESO2(REGION 3) 5 3
C 20 USER ESO2 = ESO2(REGION 3) 5 4
C
C NOTE : NSIM=5-20 ONLY USED IF ISCENGEN=1 (I.E., USER PLANS TO
C GO INTO SCENGEN AFTER MAGICC).
C
DO 1 NSIM=1,NUMSIMS
NESO2=1+INT((NSIM-0.1)/4.0)
NCLIM=NSIM
C
C RE-SET NCLIM FOR SULPHATE PATTERN WEIGHT CASES (NSIM.GE.5).
C
IF(NSIM.GE.5)NCLIM=NSIM-4
IF(NSIM.GE.9)NCLIM=NSIM-8
IF(NSIM.GE.13)NCLIM=NSIM-12
IF(NSIM.GE.17)NCLIM=NSIM-16
C
C RE-SET NCLIM=4 (USER CASE) IF ONLY ONE SIMULATION (ISCENGEN=9).
C
IF(ISCENGEN.EQ.9)NCLIM=4
C
C ****************************************************************
C
IF(NESO2.EQ.1)THEN
DO KE=226,KEND+1
ESO2(KE)=ESO2SUM(KE)
END DO
ENDIF
C
IF(NESO2.EQ.2)THEN
DO KE=226,KEND+1
ESO2(KE)=es1990
END DO
ENDIF
C
IF(NESO2.EQ.3)THEN
DO KE=226,KEND+1
ESO2(KE)=ESO21(KE)
END DO
ENDIF
C
IF(NESO2.EQ.4)THEN
DO KE=226,KEND+1
ESO2(KE)=ESO22(KE)
END DO
ENDIF
C
IF(NESO2.EQ.5)THEN
DO KE=226,KEND+1
ESO2(KE)=ESO23(KE)
END DO
ENDIF
C
C ****************************************************************
C
C SET CLIMATE SENSITIVITY AND ICE MELT PARAMETERS
C
IF(NCLIM.EQ.1)THEN ! LOW
TE = 1.5
ZI0 = 30.0
TAUI = 150.0
DTSTAR= 0.9
DTEND = 4.5
BETAG = 0.01
BETA1 = -0.045
BETA2 = 0.0
DB3 = -0.04
ENDIF
C
IF(NCLIM.EQ.2)THEN ! MID
TE = 2.5
ZI0 = 30.0
TAUI = 100.0
DTSTAR= 0.7
DTEND = 3.0
BETAG = 0.03
BETA1 = -0.03
BETA2 = 0.01
DB3 = 0.01
ENDIF
C
IF(NCLIM.EQ.3)THEN ! HIGH
TE = 4.5
ZI0 = 30.0
TAUI = 50.0
DTSTAR= 0.6
DTEND = 2.5
BETAG = 0.05
BETA1 = -0.015
BETA2 = 0.02
DB3 = 0.06
ENDIF
C
IF(NCLIM.EQ.4)THEN ! USER
TE = DT2XUSER
ZI0 = ZI0IN
TAUI = TAUIN
DTSTAR= DTSTARIN
DTEND = DTENDIN
BETAG = BETAGIN
BETA1 = BETA1IN
BETA2 = BETA2IN
DB3 = DB3IN
ENDIF
C
IF(IXLAM.EQ.1)THEN
CALL LAMCALC(Q2X,FL(1),FL(2),XKLO,XKNS,TE,RLO,XLAMO,XLAML)
ENDIF
C
XLAM=Q2X/TE
C
CALL INIT
C
WRITE (8,179)
WRITE (8,176) NSIM,TE
IF(NESO2.EQ.1)WRITE(8,186)
IF(NESO2.EQ.2)WRITE(8,187)
IF(NESO2.EQ.3)WRITE(8,188)
IF(NESO2.EQ.4)WRITE(8,189)
IF(NESO2.EQ.5)WRITE(8,190)
IF(IVARW.EQ.0)THEN
WRITE (8,122)
ENDIF
IF(IVARW.EQ.1)THEN
WRITE (8,123) TW0NH
WRITE (8,124) TW0SH
ENDIF
IF(IVARW.EQ.2)THEN
WRITE (8,125)
ENDIF
C
WRITE (8,12) XKNS,XKLO
WRITE (8,120) HM,YK
WRITE (8,121) PI,W0
C
WRITE (8,910) DTSTAR,DTEND,taui,zi0
WRITE (8,912) betag
WRITE (8,913) beta1,beta2,db3
IF(IXLAM.EQ.1)THEN
WRITE(8,914) RLO,XLAML,XLAMO
IF(XLAML.LT.0.0)WRITE(8,916)
ELSE
WRITE(8,915) XLAM
ENDIF
C
C ***********************************************************
C
CALL RUNMOD
C
C ***********************************************************
C
C EXTRA CALL TO RUNMOD TO GET FINAL FORCING VALUES FOR K=KEND
C WHEN DT=1.0
C
C IF((K.EQ.KEND).AND.(DT.EQ.1.0))CALL RUNMOD
C
C SAVE SULPHATE AEROSOL FORCINGS IN FIRST PASS THROUGH OF NSIM
C LOOP, WHEN TOTAL SO2 EMISSIONS ARE BEING USED.
C
IF(NSIM.EQ.1)THEN
DO K=1,KEND
QSO2SAVE(K)=QSO2(K)
QDIRSAVE(K)=QDIR(K)
END DO
ENDIF
C
C PRINT OUT RESULTS
C
DT1 = TGAV(226)-TGAV(116)
DMSL1 = SLT(226)-SLT(116)
DTNH1 = TNHAV(226)-TNHAV(116)
DTSH1 = TSHAV(226)-TSHAV(116)
DTLAND = TLAND(226)-TLAND(116)
DTOCEAN= TOCEAN(226)-TOCEAN(116)
DTNHO = TNHO(226)-TNHO(116)
DTSHO = TSHO(226)-TSHO(116)
DTNHL = TNHL(226)-TNHL(116)
DTSHL = TSHL(226)-TSHL(116)
C
WRITE (8,140) DT1,DMSL1
WRITE (8,141) DTNHL,DTNHO,DTSHL,DTSHO
WRITE (8,142) DTNH1,DTSH1,DTLAND,DTOCEAN
WRITE (8,15) KYRREF
WRITE (8,16)
C
IF(IVARW.EQ.1)THEN
WRITE(8,178)TE
ELSE
WRITE(8,177)TE
ENDIF
C
IF(NCLIM.EQ.1)WRITE(8,161)
IF(NCLIM.EQ.2)WRITE(8,162)
IF(NCLIM.EQ.3)WRITE(8,163)
IF(NCLIM.EQ.4)WRITE(8,164)
C
IF(IUSERGAS.EQ.1)THEN
WRITE(8,166)
ELSE
WRITE(8,167)
ENDIF
C
IF(NOUT.EQ.1)WRITE(8,171)
IF(NOUT.EQ.2)WRITE(8,172)
IF(NOUT.EQ.3)WRITE(8,173)
IF(NOUT.EQ.4)WRITE(8,174)
IF(NOUT.EQ.5)THEN
WRITE(8,175)
ENDIF
C
C PRINTOUT OPTIONS
C
IF(NOUT.EQ.1)THEN
COL9 = TNHAV(226)
COL10 = TSHAV(226)
ENDIF
IF(NOUT.EQ.2.OR.NOUT.EQ.5)THEN
COL9 = TLAND(226)
COL10 = TOCEAN(226)
IF(COL10.NE.0.0)THEN
COL11= COL9/COL10
ELSE
COL11= 9.999
ENDIF
ENDIF
IF(NOUT.EQ.3.OR.NOUT.EQ.4)THEN
COL9 = TEQU(226)-TGAV(226)
COL10 = TDEEP(226)
ENDIF
C
TEOUT=TEQU(226)
IF(IQCONST.GE.1)TEOUT=(TE/Q2X)*QCONST
C
IF(NOUT.EQ.1)THEN
WRITE(8,181)QGLOBE(226),TEOUT,TGAV(226),EX(226),SLI(226),
& SLG(226),SLA(226),SLT(226),COL9,COL10,WNH(226),WSH(226)
ENDIF
IF(NOUT.EQ.2) THEN
WRITE(8,182)QGLOBE(226),TEOUT,TGAV(226),EX(226),SLI(226),
& SLG(226),SLA(226),SLT(226),COL9,COL10,COL11,WNH(226),WSH(226)
ENDIF
IF(NOUT.EQ.3)THEN
WRITE(8,183)QGLOBE(226),TEOUT,TGAV(226),EX(226),SLI(226),
& SLG(226),SLA(226),SLT(226),COL9,COL10,WNH(226),WSH(226)
ENDIF
IF(NOUT.EQ.4)THEN
C
C CONVERT W/M**2 TO EQUIV CO2 RELATIVE TO END-1765 CO2 CONC
C
EQUIVCO2=COBS(1)*EXP(QGLOBE(226)/QXX)
WRITE(8,184)EQUIVCO2,TEOUT,TGAV(226),EX(226),SLI(226),
& SLG(226),SLA(226),SLT(226),COL9,COL10,WNH(226),WSH(226)
ENDIF
IF(NOUT.EQ.5)THEN
WRITE(8,185)QGLOBE(226),TGAV(226),COL11,SLT(226),EX(226),
& SLI(226),SLG(226),SLA(226),WNH(226)
ENDIF
C
C ****************************************************************
C
C MAIN PRINT OUT LOOP
C
TE1 = 0.
TT1 = 0.
TN1 = 0.
TS1 = 0.
C
QR = QGLOBE(KREF)
XPRT = FLOAT(ITEMPRT)
NPRT = INT(225./XPRT +0.01)
MPRT = NPRT*ITEMPRT
KYEAR0= 1990-MPRT
C
TREFSUM=0.0
C
DO 987 K=1,KEND
KYEAR=1764+K
QK = QGLOBE(K)
C
IF(IQCONST.EQ.0)THEN
Q1 = QK-QR
TEQUIL = TEQU(K)
TEQUIL0= TEQU(KREF)
ENDIF
IF(IQCONST.EQ.1)THEN
Q1 = QCONST
TEQUIL = (TE/Q2X)*QCONST
TEQUIL0= 0.0
ENDIF
IF(IQCONST.EQ.2)THEN
Q1 = RAMPDQDT*(K-1)
IF(Q1.GT.QCONST)Q1=QCONST
TEQUIL = (TE/Q2X)*Q1
TEQUIL0= 0.0
ENDIF
C
TE1 = TEQUIL-TEQUIL0
TT1 = TGAV(K)-TGAV(KREF)
TN1 = TNHAV(K)-TNHAV(KREF)
TS1 = TSHAV(K)-TSHAV(KREF)
C
C CALCULATE 1961-1990 MEAN TEMPERATURE AS REFERENCE LEVEL FOR
C CALCULATION OF INPUT INTO SCENGEN DRIVER FILES.
C NOTE : TREF DEPENDS ON CLIMATE MODEL PARAMS (I.E. ON NCLIM)
C
IF(K.GE.197.AND.K.LE.226)TREFSUM=TREFSUM+TGAV(K)
IF(K.EQ.226)TREF(NCLIM)=TREFSUM/30.
C
C PRINTOUT OPTIONS
C
IF(NOUT.EQ.1)THEN
COL9 = TN1
COL10 = TS1
ENDIF
C
IF(NOUT.EQ.2.OR.NOUT.EQ.5)THEN
COL9 = TLAND(K)-TLAND(KREF)
COL10 = TOCEAN(K)-TOCEAN(KREF)
IF(TOCEAN(K).NE.0.0)THEN
COL11= TLAND(K)/TOCEAN(K)
ELSE
COL11= 9.999
ENDIF
ENDIF
C
IF(NOUT.EQ.3.OR.NOUT.EQ.4)THEN
COL9 = TE1-TT1
COL10 = TDEEP(K)-TDEEP(KREF)
ENDIF
C
EX1=EX(K) -EX(KREF)
SI1=SLI(K)-SLI(KREF)
SG1=SLG(K)-SLG(KREF)
SA1=SLA(K)-SLA(KREF)
ST1=SLT(K)-SLT(KREF)
C
C PUT TEMPERATURE AND SEA LEVEL RESULTS FOR FULL GLOBAL FORCING
C INTO DISPLAY OUTPUT FILES
C
IF(ISCENGEN.NE.9)THEN
IF(NSIM.EQ.1)THEN
TEMLO(K) = TT1
SLLO(K) = ST1
ENDIF
ENDIF
C
IF(NSIM.EQ.2)THEN
TEMMID(K) = TT1
SLMID(K) = ST1
ENDIF
C
IF(NSIM.EQ.3)THEN
TEMHI(K) = TT1
SLHI(K) = ST1
ENDIF
C
IF((ISCENGEN.EQ.9).OR.(NSIM.EQ.4))THEN
TEMUSER(K)= TT1
SLUSER(K) = ST1
ENDIF
C
C RESULTS FOR ESO2 CONST AFTER 1990 STORED ONLY FOR MID CLIMATE CASE.
C ZERO VALUES STORED IF ISCENGEN=0 OR =9
C
IF(ISCENGEN.EQ.0.OR.ISCENGEN.EQ.9)THEN
TEMNOSO2(K)= 0.0
SLNOSO2(K) = 0.0
ENDIF
C
IF(NSIM.EQ.6)THEN
TEMNOSO2(K)= TT1
SLNOSO2(K) = ST1
ENDIF
C
C IF(NOUT.EQ.1.)THEN
C WRITE (8,191) KYEAR,Q1,TE1,TT1,EX1,SI1,SG1,SA1,ST1,COL9,
C & COL10,WNH(K),WSH(K),KYEAR
C ENDIF
C
C PRINT OUT FLAG IS KKKK=1
C
KKKK=0
C
C ALWAYS PRINT OUT 1765, IY0 AND 1990 VALUES
C
IF(KYEAR.EQ.1764.OR.KYEAR.EQ.IY0.OR.KYEAR.EQ.1990)KKKK=1
C
IF(KYEAR.GE.IY0)THEN
PRIN=(KYEAR-KYEAR0+0.01)/XPRT
BIT=PRIN-INT(PRIN)
IF(PRIN.GT.0.0.AND.BIT.LT.0.02)KKKK=1
IF(KKKK.EQ.1)THEN
C
C ADD CONSTANT TO ALL TEMPS FOR IPCC DETEX TIME FIGURE
C
TT1=TT1+TOFFSET
C
IF(NOUT.EQ.1.)THEN
WRITE (8,191) KYEAR,Q1,TE1,TT1,EX1,SI1,SG1,SA1,ST1,COL9,
& COL10,WNH(K),WSH(K),KYEAR
ENDIF
C
IF(NOUT.EQ.2)THEN
WRITE (8,192) KYEAR,Q1,TE1,TT1,EX1,SI1,SG1,SA1,ST1,COL9,
& COL10,COL11,WNH(K),WSH(K),KYEAR
ENDIF
C
IF(NOUT.EQ.3)THEN
WRITE (8,193) KYEAR,Q1,TE1,TT1,EX1,SI1,SG1,SA1,ST1,COL9,
& COL10,WNH(K),WSH(K),KYEAR
ENDIF
C
IF(NOUT.EQ.4)THEN
EQUIVCO2=COBS(1)*EXP(QK/QXX)
WRITE (8,194) KYEAR,EQUIVCO2,TE1,TT1,EX1,SI1,SG1,SA1,ST1,
& COL9,COL10,WNH(K),WSH(K),KYEAR
ENDIF
C
IF(NOUT.EQ.5)THEN
WRITE(8,195) KYEAR,Q1,TT1,COL11,ST1,EX1,SI1,SG1,SA1,
& WNH(K),KYEAR
ENDIF
C
ENDIF
ENDIF
987 CONTINUE
C
IF(NOUT.EQ.1)WRITE(8,171)
IF(NOUT.EQ.2)WRITE(8,172)
IF(NOUT.EQ.3)WRITE(8,173)
IF(NOUT.EQ.4)WRITE(8,174)
IF(NOUT.EQ.5)WRITE(8,175)
WRITE(8,30)
C
C **************************************************************
C
C DEFINE TEMPERATURE ARRAYS FOR WRITING TO SCENGEN DRIVER FILES.
C ARRAY SUBSCRIPT NCLIM=1,2,3,4 CORRESPONDS TO LO, MID, HIGH
C AND USER CLIMATE MODEL PARAMETER SETS.
C NOTE THAT THESE ARRAYS START WITH KSG=1 IN 1990.
C
C TSO21,2,3 ARE THE RAW TEMPERATURES. THEY ARE INTERNALLY
C INCONSISTENT BECAUSE OF THE NONLINEAR RELATIONSHIP BETWEEN
C ESO2 AND INDIRECT AEROSOL FORCING. IN OTHER WORDS, IN GENERAL,
C TALL MINUS TGHG WILL NOT EQUAL TSO21+TSO22+TSO23. TO CORRECT
C FOR THIS, SCALE TSO2i (GIVING XSO2i) BY (TALL-TGHG)/(SUM TSO2i).
C
DO K=197,KEND
KSG=K-196
IF(NESO2.EQ.1)THEN
TALL(NCLIM,KSG)=TGAV(K)-TREF(NCLIM)
ENDIF
C
IF(NESO2.EQ.2)THEN
TGHG(NCLIM,KSG)=TGAV(K)-TREF(NCLIM)
ENDIF
C
IF(NESO2.EQ.3)THEN
TSO21(NCLIM,KSG)=TGAV(K)-TGHG(NCLIM,KSG)-TREF(NCLIM)
IF(K.LT.226)TSO21(NCLIM,KSG)=0.0
ENDIF
C
IF(NESO2.EQ.4)THEN
TSO22(NCLIM,KSG)=TGAV(K)-TGHG(NCLIM,KSG)-TREF(NCLIM)
IF(K.LT.226)TSO22(NCLIM,KSG)=0.0
ENDIF
C
IF(NESO2.EQ.5)THEN
TSO23(NCLIM,KSG)=TGAV(K)-TGHG(NCLIM,KSG)-TREF(NCLIM)
IF(K.LT.226)TSO23(NCLIM,KSG)=0.0
ENDIF
END DO
C
C end of NSIM loop
C
1 CONTINUE
C
C **************************************************************
C **************************************************************
C
C SCALE TEMPERATURES TO GET DRIVER TEMPERATURES, XSO2i
C
IF(ISCENGEN.EQ.1)THEN
C
DO NCLIM=1,4
DO K=197,KEND
KSG=K-196
XTSUM=TALL(NCLIM,KSG)-TGHG(NCLIM,KSG)
YTSUM=TSO21(NCLIM,KSG)+TSO22(NCLIM,KSG)+TSO23(NCLIM,KSG)
SCALER=1.0
IF(YTSUM.NE.0.0)THEN
SCALER=XTSUM/YTSUM
ENDIF
XSO21(NCLIM,KSG)=TSO21(NCLIM,KSG)*SCALER
XSO22(NCLIM,KSG)=TSO22(NCLIM,KSG)*SCALER
XSO23(NCLIM,KSG)=TSO23(NCLIM,KSG)*SCALER
END DO
END DO
C
C **************************************************************
C
C WRITE TEMPERATURES TO OLD AND NEW SCENGEN DRIVER FILES.
C UNSCALED TEMPS GO TO OLD FILES, SCALED TEMPS TO NEW FILES.
C
OPEN(UNIT=10,FILE='lodrive.raw' ,STATUS='UNKNOWN')
OPEN(UNIT=11,FILE='middrive.raw',STATUS='UNKNOWN')
OPEN(UNIT=12,FILE='hidrive.raw' ,STATUS='UNKNOWN')
OPEN(UNIT=13,FILE='usrdrive.raw',STATUS='UNKNOWN')
C
OPEN(UNIT=14,FILE='lodrive.out' ,STATUS='UNKNOWN')
OPEN(UNIT=15,FILE='middrive.out',STATUS='UNKNOWN')
OPEN(UNIT=16,FILE='hidrive.out' ,STATUS='UNKNOWN')
OPEN(UNIT=17,FILE='usrdrive.out',STATUS='UNKNOWN')
C
DO NCLIM=1,4
DO KSG=1,KEND-196
KYY=KSG+1960
C
IF(NCLIM.EQ.1)THEN
IF(KSG.EQ.1)THEN
WRITE(10,937)mnem
WRITE(10,930)
WRITE(10,934)
WRITE(10,936)TREF(NCLIM)
WRITE(14,937)mnem
WRITE(14,930)
WRITE(14,934)
WRITE(14,936)TREF(NCLIM)
ENDIF
WRITE(10,935)KYY,TGHG(NCLIM,KSG),TSO21(NCLIM,KSG),
& TSO22(NCLIM,KSG),TSO23(NCLIM,KSG),TALL(NCLIM,KSG)
WRITE(14,935)KYY,TGHG(NCLIM,KSG),XSO21(NCLIM,KSG),
& XSO22(NCLIM,KSG),XSO23(NCLIM,KSG),TALL(NCLIM,KSG)
ENDIF
C
IF(NCLIM.EQ.2)THEN
IF(KSG.EQ.1)THEN
WRITE(11,937)mnem
WRITE(11,931)
WRITE(11,934)
WRITE(11,936)TREF(NCLIM)
WRITE(15,937)mnem
WRITE(15,931)
WRITE(15,934)
WRITE(15,936)TREF(NCLIM)
ENDIF
WRITE(11,935)KYY,TGHG(NCLIM,KSG),TSO21(NCLIM,KSG),
& TSO22(NCLIM,KSG),TSO23(NCLIM,KSG),TALL(NCLIM,KSG)
WRITE(15,935)KYY,TGHG(NCLIM,KSG),XSO21(NCLIM,KSG),
& XSO22(NCLIM,KSG),XSO23(NCLIM,KSG),TALL(NCLIM,KSG)
ENDIF
C
IF(NCLIM.EQ.3)THEN
IF(KSG.EQ.1)THEN
WRITE(12,937)mnem
WRITE(12,932)
WRITE(12,934)
WRITE(12,936)TREF(NCLIM)
WRITE(16,937)mnem
WRITE(16,932)
WRITE(16,934)
WRITE(16,936)TREF(NCLIM)
ENDIF
WRITE(12,935)KYY,TGHG(NCLIM,KSG),TSO21(NCLIM,KSG),
& TSO22(NCLIM,KSG),TSO23(NCLIM,KSG),TALL(NCLIM,KSG)
WRITE(16,935)KYY,TGHG(NCLIM,KSG),XSO21(NCLIM,KSG),
& XSO22(NCLIM,KSG),XSO23(NCLIM,KSG),TALL(NCLIM,KSG)
ENDIF
C
IF(NCLIM.EQ.4)THEN
IF(KSG.EQ.1)THEN
WRITE(13,937)mnem
WRITE(13,933)
WRITE(13,934)
WRITE(13,936)TREF(NCLIM)
WRITE(17,937)mnem
WRITE(17,933)
WRITE(17,934)
WRITE(17,936)TREF(NCLIM)
ENDIF
WRITE(13,935)KYY,TGHG(NCLIM,KSG),TSO21(NCLIM,KSG),
& TSO22(NCLIM,KSG),TSO23(NCLIM,KSG),TALL(NCLIM,KSG)
WRITE(17,935)KYY,TGHG(NCLIM,KSG),XSO21(NCLIM,KSG),
& XSO22(NCLIM,KSG),XSO23(NCLIM,KSG),TALL(NCLIM,KSG)
ENDIF
C
END DO
END DO
C
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
CLOSE(14)
CLOSE(15)
CLOSE(16)
CLOSE(17)
C
ENDIF
C
C **************************************************************
C **************************************************************
C
C WRITE TEMPERATURES TO MAG DISPLAY FILE
C
open(unit=9,file='temps.dis',status='UNKNOWN')
C
WRITE (9,213)
C
C PRINTOUT INTERVAL FOR DISPLAY PURPOSES SET BY IDIS IN MAGEXTRA.CFG
C
DO K=1,KEND,IDIS
IYEAR=1764+K
C PRIN=FLOAT(K+4)/IDIS+0.01
C IF((PRIN-INT(PRIN)).LT.0.05)THEN
WRITE (9,226) IYEAR,TEMUSER(K),TEMLO(K),TEMMID(K),
& TEMHI(K),TEMNOSO2(K)
C ENDIF
END DO
C
WRITE (9,213)
C
CLOSE(9)
C
C **************************************************************
C
C
C WRITE SEALEVEL CHANGES TO MAG DISPLAY FILE
C
open(unit=9,file='sealev.dis',status='UNKNOWN')
C
WRITE (9,214)
C
C PRINTOUT INTERVAL FOR DISPLAY PURPOSES SET BY IDIS IN MAGEXTRA.CFG
C
DO K=1,KEND,IDIS
IYEAR=1764+K
C PRIN=FLOAT(K+4)/IDIS+0.01
C IF((PRIN-INT(PRIN)).LT.0.05)THEN
WRITE (9,227) IYEAR,SLUSER(K),SLLO(K),SLMID(K),
& SLHI(K),SLNOSO2(K)
C ENDIF
END DO
C
WRITE (9,214)
C
CLOSE(9)
C
C **************************************************************
C
C PRINT OUT EMISSIONS, CONCS AND FORCING DETAILS
C
C PRINT OUT INPUT EMISSIONS
C
WRITE (8,30)
WRITE (8,31)
WRITE (8,30)
WRITE (8,23)
WRITE (8,231)
WRITE (8,21)
C
C PRINTOUT INTERVAL IS DET BY VALUE OF IEMPRT
C
C NOTE THAT THIS VERSION OF THE CODE CAN GIVE ODD ECH4 VALUES,
C SINCE IT USES A SLIGHTLY UNREALISTIC VALUE OF THE 1990 CH4
C LIFETIME IN ORDER TO GET BACK TO THE CH4 CONCENTRATION
C PROJECTIONS USED IN CH. 6 OF THE IPCC SAR. THE BEST GUESS
C ECH4 IN 1990 IS 571.6, SO THE VALUES PRINTED OUT BELOW HAVE
C A CORRECTION FACTOR ADDED TO GET TO THIS VALUE.
C
CH4CORN=ECH4(226)-571.6
C
DO K=226,KEND
IYEAR=1764+K
PRIN=FLOAT(K+4)/IEMPRT+0.01
IF((PRIN-INT(PRIN)).LT.0.05)THEN
ES1=ESO21(K)-es1990
ES2=ESO22(K)-es1990
ES3=ESO23(K)-es1990
EST=ES1+ES2+ES3
EMCORR=ECH4(K)-CH4CORN
WRITE (8,222) IYEAR,EF(K),EDNET(K),EMCORR,EN2O(K),
& ES1,ES2,ES3,EST,IYEAR
ENDIF
END DO
C
WRITE (8,21)
WRITE (8,30)
WRITE (8,31)
WRITE (8,30)
C
C WRITE EMISSIONS TO MAG DISPLAY FILE
C
open(unit=9,file='emiss.dis',status='UNKNOWN')
C
WRITE (9,212)
C
C PRINTOUT INTERVAL FOR DISPLAY PURPOSES SET BY IDIS IN MAGEXTRA.CFG
C NOTE THAT ESO2(K) IS OVERWRITTEN IN LAST LOOP OF CLIMATE MODEL
C SIMULATIONS, BUT ESO2i(K) REMAIN AS INPUTTED.
C
DO K=226,KEND,IDIS
IYEAR=1764+K
C PRIN=FLOAT(K+4)/IDIS+0.01
C IF((PRIN-INT(PRIN)).LT.0.05)THEN
WRITE (9,225) IYEAR,EF(K),EDNET(K),ECH4(K),EN2O(K),
+ ESO21(K),ESO22(K),ESO23(K)
C ENDIF
END DO
C
WRITE (9,212)
C
CLOSE(9)
C
C **************************************************************
C
C PRINT OUT USER CARBON CYCLE DETAILS
C
WRITE(8,24)
WRITE(8,800)R(1)
WRITE(8,801)R(2)
WRITE(8,802)R(3)
WRITE(8,803)DUSER,R(4)
WRITE(8,804)
WRITE(8,805)
C
if(iMeth.eq.1) then
write(8,806)
else
write(8,807)
endif
C
MID=0
IF(MID.NE.1)WRITE(8,810)
IF(MID.EQ.1)WRITE(8,811)
WRITE(8,812)
C
C PRINTOUT INTERVAL IS DET BY VALUE OF ICO2PRT
C
DO K=226,KEND
IYEAR=1764+K
PRIN=FLOAT(K+4)/ICO2PRT+0.01
IF((PRIN-INT(PRIN)).LT.0.05)THEN
CONCOUT=CCO2(4,K)
IF(MID.EQ.1)CONCOUT=(CCO2(4,K-1)+CCO2(4,K))/2.
C
IF(IMETH.EQ.0)THEN
TOTE=EF(K)+EDNET(K)
ELSE
TOTE=EF(K)+EDNET(K)+EMETH(K)
ENDIF
IF(TOTE.EQ.0.0)THEN
IF(DELMASS(4,K).EQ.0.0)THEN
ABX=1.0
ELSE
ABX=DELMASS(4,K)/ABS(DELMASS(4,K))
ENDIF
ABFRAC(4,K)=ABX*9.999
ELSE
ABFRAC(4,K)=DELMASS(4,K)/TOTE
ENDIF
IF(ABFRAC(4,K).GT.9.999)ABFRAC(4,K)=9.999
IF(ABFRAC(4,K).LT.-9.999)ABFRAC(4,K)=-9.999
C
ECH4OX=EMETH(K)
IF(IMETH.EQ.0)ECH4OX=0.0
WRITE(8,813)IYEAR,TOTE,EF(K),ECH4OX,EDNET(K),EDGROSS(4,K),
& FOC(4,K),ABFRAC(4,K),PL(4,K),HL(4,K),SOIL(4,K),CONCOUT,
& DELMASS(4,K),IYEAR
C
ENDIF
END DO
WRITE(8,812)
C
C **************************************************************
C
C PRINT OUT CONCENTRATIONS
C
WRITE (8,30)
WRITE (8,31)
WRITE (8,30)
WRITE (8,20)
WRITE (8,201)
WRITE (8,210)
C
C PRINTOUT INTERVAL IS DET BY VALUE OF ICONCPRT
C
DO K=1,KEND
IYEAR=1764+K
PRIN=FLOAT(K+4)/ICONCPRT+0.01
IF((PRIN-INT(PRIN)).LT.0.05)THEN
C
C CONVERT END OF YEAR TO MIDYEAR CONCS
C
CO2MID=(CO2(K)+CO2(K-1))/2.
IF(K.GT.1)THEN
CH4MID=(CH4(K)+CH4(K-1))/2.
CN2OMID=(CN2O(K)+CN2O(K-1))/2.
ELSE
CH4MID=CH4(K)
CN2OMID=CN2O(K)
ENDIF
C
IF(K.GE.226)THEN
CH4LMID=(CH4L(K)+CH4L(K-1))/2.
CH4BMID=(CH4B(K)+CH4B(K-1))/2.
CH4HMID=(CH4H(K)+CH4H(K-1))/2.
CO2LMID=(CCO2(1,K)+CCO2(1,K-1))/2.
CO2BMID=(CCO2(2,K)+CCO2(2,K-1))/2.
CO2HMID=(CCO2(3,K)+CCO2(3,K-1))/2.
C
IF(K.EQ.226)THEN
TOR=T90USER
ELSE
TOR=TCH4(K)
ENDIF
C
WRITE (8,220) IYEAR,CO2MID,CH4MID,CN2OMID,
+ ch4lMID,ch4bMID,ch4hMID,
+ CO2LMID,CO2BMID,CO2HMID,IYEAR,TOR
ELSE
WRITE (8,221) IYEAR,CO2MID,CH4MID,CN2OMID,IYEAR
ENDIF
ENDIF
END DO
C
WRITE (8,210)
WRITE (8,201)
C
C WRITE CONCENTRATIONS TO MAG DISPLAY FILE
C
open(unit=9,file='concs.dis',status='UNKNOWN')
C
WRITE (9,211)
C
C PRINTOUT INTERVAL FOR DISPLAY PURPOSES SET BY IDIS IN MAGEXTRA.CFG
C
DO K=1,KEND
IYEAR=1764+K
PRIN=FLOAT(K+4)/IDIS+0.01
IF((PRIN-INT(PRIN)).LT.0.05)THEN
if(k.ge.226) then
WRITE (9,223) IYEAR,CO2(K),cco2(1,k),cco2(2,k),
& cco2(3,k),CH4(K),ch4l(k),ch4b(k),ch4h(k),CN2O(K)
else
WRITE (9,224) IYEAR,CO2(K),CH4(K),CN2O(K)
endif
ENDIF
END DO
C
WRITE (9,211)
C
CLOSE(9)
C
C **************************************************************
C
C PRINT OUT PERCENT RADIATIVE FORCING BREAKDOWN TO QDETAILS.OUT
C
write(88,30)
write(88,31)
WRITE(88,30)
write(88,47)
write(88,48)
write(88,50)
write(88,51)
c
coeff = 100./qtot86
write(88,52) 1765,1850,
& qco2(86),qco2(86)*coeff,qm(86),qm(86)*coeff,
& qn(86),qn(86)*coeff,qcfc(86),qcfc(86)*coeff,
& QSO2SAVE(86),QSO2SAVE(86)*coeff,qtot86

C
C PRINTOUT INTERVAL IS DET BY VALUE OF IQDECPRT
C
do k=1860,IYEND,IQDECPRT
iyr = k-1990+226
coeff = qtot(iyr)-qtot(iyr-10)
if(abs(coeff).lt.0.0001) coeff = 0.0001
coeff = 100./coeff
write(88,52) k-10,k,
& qco2(iyr)-qco2(iyr-10),(qco2(iyr)-qco2(iyr-10))*coeff,
& qm (iyr)-qm (iyr-10),(qm (iyr)-qm (iyr-10))*coeff,
& qn (iyr)-qn (iyr-10),(qn (iyr)-qn (iyr-10))*coeff,
& qcfc(iyr)-qcfc(iyr-10),(qcfc(iyr)-qcfc(iyr-10))*coeff,
& QSO2SAVE(iyr)-QSO2SAVE(iyr-10),
& (QSO2SAVE(iyr)-QSO2SAVE(iyr-10))*coeff,
& qtot(iyr)-qtot(iyr-10)
end do
C
WRITE(88,51)
C
C **************************************************************
C
C PRINTOUT TABLES OF SPATIAL DELTA-Q BREAKDOWN TO QDETAILS.OUT
C
WRITE(88,30)
WRITE(88,31)
WRITE(88,920)
WRITE(88,921)
C
DO K=1770,IYEND,INSLOPRT
IYR = K-1990+226
C
C QQ VALUES ARE RAD FORCING IN W/M**2 FOR EACH BOX. TO GET W/M**2
C COMPONENTS OF TOTAL GLOBAL FORCING, NEED TO MULTIPLY BY RELATIVE
C AREAS.
C
FFNHO=QQQNHO(IYR)*FNO
FFNHL=QQQNHL(IYR)*FNL
FFSHO=QQQSHO(IYR)*FSO
FFSHL=QQQSHL(IYR)*FSL
C
FFTOT=FFNHO+FFNHL+FFSHO+FFSHL
WRITE(88,922)K,QQQNHO(IYR),FFNHO,QQQNHL(IYR),FFNHL,
& QQQSHO(IYR),FFSHO,QQQSHL(IYR),FFSHL,FFTOT
C
END DO
C
WRITE(88,921)
C
C **************************************************************
C
C PRINTOUT TABLES OF DELTA-Q FROM 1990 AND 1765.
C FIRST, DELTA-Q FROM 1990.
C
write(8,30)
write(8,31)
WRITE(8,30)
write(8,55)
write(8,56)
write(8,57)
C
QQQCO2R = (qco2(226)+qco2(225))/2.
QQQMR = (QM(226)+QM(225))/2.
QQQNR = (QN(226)+QN(225))/2.
QQQCFCR = (QCFC(226)+QCFC(225))/2.
QQQSO2R = (QSO2SAVE(226)+QSO2SAVE(225))/2.
QQQDIRR = (QDIRSAVE(226)+QDIRSAVE(225))/2.
C
C NOTE SPECIAL CASE FOR QOZ BECAUSE OF NONLINEAR CHANGE OVER 1990
C
QQQOZR = QOZ(226)
C
QQQBIOR = (QBIO(226)+QBIO(225))/2.
QQQMO3R = (QCH4O3(226)+QCH4O3(225))/2.
C
C QTOTREF=(QTOT(226)+QTOT(225))/2.0
C QTOTR=QTOT(225)+0.5*(QTOT(225)-QTOT(224))
C IF(ICOMP.EQ.1.OR.icold.EQ.1)QTOTREF=QTOTR
C
C PRINTOUT INTERVAL IS DET BY VALUE OF IQGASPRT
C
C FORCING CHANGES FROM MID 1990
C
DO K=1990,IYEND,IQGASPRT
IYR = K-1990+226
IYRP=IYR-1
C
DELQCO2 = (QCO2(IYR)+QCO2(IYRP))/2.-QQQCO2R
DELQM = (QM(IYR)+QM(IYRP))/2. -QQQMR
DELQN = (QN(IYR)+QN(IYRP))/2. -QQQNR
DELQCFC = (QCFC(IYR)+QCFC(IYRP))/2.-QQQCFCR
C
DELQSO2 = (QSO2SAVE(IYR)+QSO2SAVE(IYRP))/2.-QQQSO2R
DELQDIR = (QDIRSAVE(IYR)+QDIRSAVE(IYRP))/2.-QQQDIRR
DELQIND = DELQSO2-DELQDIR
C
C NOTE SPECIAL CASE FOR QOZ BECAUSE OF NONLINEAR CHANGE OVER 1990
C
IF(IYR.EQ.226)THEN
QOZMID= QOZ(IYR)
ELSE
QOZMID= (QOZ(IYR)+QOZ(IYRP))/2.
ENDIF
DELQOZ = QOZMID-QQQOZR
C
DELQBIO = (QBIO(IYR)+QBIO(IYRP))/2.-QQQBIOR
DELQTOT = DELQCO2+DELQM+DELQN+DELQCFC+DELQSO2+DELQBIO+DELQOZ
C
DQCH4O3 = (QCH4O3(IYR)+QCH4O3(IYRP))/2.-QQQMO3R
DELQM = DELQM-DQCH4O3
DELQOZ = DELQOZ+DQCH4O3
C
C DELQM=DELQM+DELQN
C DELQN=DELQTOT-DELQCO2-DELQCFC
C DELQOZ=DELQDIR+DELQIND+DELQBIO
C
WRITE(8,571)K,DELQCO2,DELQM,DELQN,DELQCFC,DELQOZ,
& DELQDIR,DELQIND,DELQBIO,DELQTOT,DQCH4O3
END DO
C
C ************************************************************
C
C WRITE FORCING CHANGES FROM MID-1990 TO MAG DISPLAY FILE
C
open(unit=9,file='forcings.dis',status='UNKNOWN')
C
WRITE (9,215)
C
C PRINTOUT INTERVAL FOR DISPLAY PURPOSES SET BY IDIS IN MAGEXTRA.CFG
C
DO K=1990,IYEND,IDIS
IYR = K-1990+226
IYRP=IYR-1
C
DELQCO2 = (QCO2(IYR)+QCO2(IYRP))/2.-QQQCO2R
DELQM = (QM(IYR)+QM(IYRP))/2. -QQQMR
DELQN = (QN(IYR)+QN(IYRP))/2. -QQQNR
DELQCFC = (QCFC(IYR)+QCFC(IYRP))/2.-QQQCFCR
C
DELQSO2 = (QSO2SAVE(IYR)+QSO2SAVE(IYRP))/2.-QQQSO2R
DELQDIR = (QDIRSAVE(IYR)+QDIRSAVE(IYRP))/2.-QQQDIRR
DELQIND = DELQSO2-DELQDIR
C
C NOTE SPECIAL CASE FOR QOZ BECAUSE OF NONLINEAR CHANGE OVER 1990
C
IF(IYR.EQ.226)THEN
QOZMID= QOZ(IYR)
ELSE
QOZMID= (QOZ(IYR)+QOZ(IYRP))/2.
ENDIF
DELQOZ = QOZMID-QQQOZR
C
DELQBIO = (QBIO(IYR)+QBIO(IYRP))/2.-QQQBIOR
DELQTOT = DELQCO2+DELQM+DELQN+DELQCFC+DELQSO2+DELQBIO+DELQOZ
C
DQCH4O3 = (QCH4O3(IYR)+QCH4O3(IYRP))/2.-QQQMO3R
DELQM = DELQM-DQCH4O3
DELQOZ = DELQOZ+DQCH4O3
C
WRITE(9,228)K,DELQCO2,DELQM,DELQN,DELQCFC,DELQOZ,
& DELQDIR,DELQIND,DELQBIO,DELQTOT
END DO
C
WRITE (9,215)
C
CLOSE(9)
C
C ************************************************************
C
C NOW PRINT OUT FORCING CHANGES FROM MID 1765.
C
WRITE(8,57)
write(8,30)
write(8,31)
WRITE(8,30)
write(8,58)
write(8,56)
write(8,57)
C
DO K=1770,IYEND,IQGASPRT
IYR = K-1990+226
IYRP=IYR-1
C
QQQSO2 = 0.0
QQQDIR = 0.0
IF(K.GT.1860)THEN
QQQSO2 = (QSO2SAVE(IYR)+QSO2SAVE(IYRP))/2.
QQQDIR = (QDIRSAVE(IYR)+QDIRSAVE(IYRP))/2.
ENDIF
QQQIND = QQQSO2-QQQDIR
C
C IF((ICOMP.EQ.1.OR.icold.EQ.1).AND.IYR.EQ.226)THEN
C QQTOT=QTOTR
C ENDIF
C
C if(k.eq.1860) then
C qqtot = qqtot - qqqso2
C endif
C
QQQCO2 = (QCO2(IYR)+QCO2(IYRP))/2.
QQQM = (QM(IYR)+QM(IYRP))/2.
QQQN = (QN(IYR)+QN(IYRP))/2.
QQQCFC = (QCFC(IYR)+QCFC(IYRP))/2.
QQQOZ = (QOZ(IYR)+QOZ(IYRP))/2.
C
C NOTE SPECIAL CASE FOR QOZ BECAUSE OF NONLINEAR CHANGE OVER 1990
C
IF(IYR.EQ.226)QQQOZ=QOZ(IYR)
C
QQQBIO = (QBIO(IYR)+QBIO(IYRP))/2.
QQQTOT = QQQCO2+QQQM+QQQN+QQQCFC+QQQSO2+QQQBIO+QQQOZ
C
QQCH4O3 = (QCH4O3(IYR)+QCH4O3(IYRP))/2.
QQQM = QQQM-QQCH4O3
QQQOZ = QQQOZ+QQCH4O3
C
C QQQM=QQQM+QQQN
C QQQN=QQQTOT-QQQCO2-QQQCFC
C QQQOZ=QQQDIR+QQQIND+QQQBIO
C
WRITE(8,571)K,QQQCO2,QQQM,QQQN,QQQCFC,QQQOZ,QQQDIR,QQQIND,
& QQQBIO,QQQTOT,QQCH4O3
END DO
C
WRITE(8,57)
C
C **************************************************************
C **************************************************************
C
10 FORMAT (/1X,'CO2-DOUBLING FORCING IN W/M**2 =',F6.3)
11 FORMAT (1X,'FNHOC= ',F4.2,' * FSHOC= ',F4.2,' * FNHLAND= ',F4.2,
+' * FSHLAND= ',F4.2)
12 FORMAT (1X,'XKNS=',F4.1,' : XKLO=',F4.1)
120 FORMAT (1X,'HM=',F5.1,'M : XK=',F6.4,'CM**2/SEC')
121 FORMAT (1X,'PI=',F6.4,' : INITIAL W=',F5.2,'M/YR')
122 FORMAT (1X,'CONSTANT W CASE')
123 FORMAT (1X,'VARIABLE W : NH W = ZERO WHEN OCEAN TEMP =',F5.1)
124 FORMAT (1X,'VARIABLE W : SH W = ZERO WHEN OCEAN TEMP =',F5.1)
125 FORMAT (1X,'VARIABLE W : NH AND SH W(t) SPECIFIED IN WINPUT.IN')
140 FORMAT (/1X,'1880-1990 CHANGES : GLOBAL DTEMP =',F7.3,
+' : DMSL =',F7.3)
141 FORMAT (1X,' DTNHL =',F7.3,' : DTNHO =',F7.3,
+' : DTSHL =',f7.3,' : DTSHO =',f7.3)
142 FORMAT (1X,' DTNH =',F7.3,' : DTSH =',F7.3,
+' : DTLAND =',f7.3,' : DTOCEAN =',f7.3)
15 FORMAT (/1X,'** TEMPERATURE AND SEA LEVEL CHANGES FROM',I5,' **')
16 FORMAT (1X,' (FIRST LINE GIVES 1765-1990 CHANGES : ',
+ 'ALL VALUES ARE MID-YEAR TO MID-YEAR)')
161 FORMAT (/1X,'LOW CLIMATE AND SEA LEVEL MODEL PARAMETERS')
162 FORMAT (/1X,'MID CLIMATE AND SEA LEVEL MODEL PARAMETERS')
163 FORMAT (/1X,'HIGH CLIMATE AND SEA LEVEL MODEL PARAMETERS')
164 FORMAT (/1X,'USER CLIMATE AND SEA LEVEL MODEL PARAMETERS')
166 FORMAT (1X,'USER GAS CYCLE MODEL PARAMETERS',/)
167 FORMAT (1X,'DEFAULT GAS CYCLE MODEL PARAMETERS',/)
171 FORMAT (1X,' YEAR DELTAQ TEQU TEMP EXPN GLAC',
+' GREENL ANTAR MSLTOT TNH TSH WNH WSH YEAR')
172 FORMAT (1X,' YEAR DELTAQ TEQU TEMP EXPN GLAC',
+' GREENL ANTAR MSLTOT TLAND TOCN TL/TO WNH WSH YEAR')
173 FORMAT (1X,' YEAR DELTAQ TEQU TEMP EXPN GLAC',
+' GREENL ANTAR MSLTOT TEQ-T TDEEP WNH WSH YEAR')
174 FORMAT (1X,' YEAR EQVCO2 TEQU TEMP EXPN GLAC',
+' GREENL ANTAR MSLTOT TEQ-T TDEEP WNH WSH YEAR')
175 FORMAT (1X,' YEAR DELTAQ TEMP TL/TO MSLTOT EXPN GLAC',
+' GREENL ANTAR WNH YEAR')
176 FORMAT (1X,'NSIM =',I3,' : DELT(2XCO2) =',F6.3,'DEGC')
177 FORMAT (/1X,'DT2X =',F5.2,' : CONSTANT W')
178 FORMAT (/1X,'DT2X =',F5.2,' : VARIABLE W')
179 FORMAT (/1X,'*************************************************')
181 FORMAT ('TO1990 ',F6.3,F7.3,F7.4,5F7.2,2F6.3,2F6.2)
182 FORMAT ('TO1990 ',F6.3,F7.3,F7.4,5F7.2,3F6.3,2F6.2)
183 FORMAT ('TO1990 ',F6.3,F7.3,F7.4,5F7.2,2F6.3,2F6.2)
184 FORMAT ('TO1990 ',F6.1,F7.3,F7.4,5F7.2,2F6.3,2F6.2)
185 FORMAT ('TO1990 ',F6.3,F8.4,F7.3,5F7.2,F6.2)
186 FORMAT (1X,'FULL GLOBAL SO2 EMISSIONS',/)
187 FORMAT (1X,'SO2 EMISSIONS CONSTANT AFTER 1990',/)
188 FORMAT (1X,'REGION 1 SO2 EMISSIONS',/)
189 FORMAT (1X,'REGION 2 SO2 EMISSIONS',/)
190 FORMAT (1X,'REGION 3 SO2 EMISSIONS',/)
191 FORMAT (1X,I5,2F7.3,F7.4,5F7.2,2F6.3,2F6.2,I6)
192 FORMAT (1X,I5,2F7.3,F7.4,5F7.2,3F6.3,2F6.2,I6)
193 FORMAT (1X,I5,2F7.3,F7.4,5F7.2,2F6.3,2F6.2,I6)
194 FORMAT (1X,I5,F7.2,F7.3,F7.4,5F7.2,2F6.3,2F6.2,I6)
195 FORMAT (1X,I5,F7.3,F8.4,F7.3,5F7.2,F6.2,I6)
C
20 FORMAT (1X,'** CONCENTRATIONS (CO2,PPMV : CH4,N2O,PPBV : REST,'
+,'PPTV. MIDYEAR VALUES) **')
201 FORMAT (5X,'<- USER MODEL CONCS ->',
+'<------ CH4 & CO2 MID CONCS & RANGES ------->')
21 FORMAT (1X,'YEAR EFOSS NETDEF CH4 N2O',
+' SO2REG1 SO2REG2 SO2REG3 SO2TOT YEAR')
210 FORMAT (1X,'YEAR CO2 CH4 N2O',
+' CH4LO CH4MID CH4HI CO2LO CO2MID CO2HI YEAR',
+' TAUCH4')
211 FORMAT (1X,'YEAR CO2USER CO2LO CO2MID CO2HI',
& ' CH4USER CH4LO CH4MID CH4HI N2O')
212 FORMAT (1X,'YEAR FOSSCO2 NETDEFOR CH4 N2O',
&' SO2-REG1 SO2-REG2 SO2-REG3')
213 FORMAT (1X,'YEAR TEMUSER TEMLO TEMMID TEMHI TEMNOSO2')
214 FORMAT (1X,'YEAR MSLUSER MSLLO MSLMID MSLHI MSLNOSO2')
215 FORMAT (1X,'YEAR CO2 CH4 N2O HALOS TROPO3',
&' SO4DIR SO4IND BIOAER TOTAL')
220 FORMAT (1X,I4,F7.1,F8.1,F7.1,3F7.0,3F8.1,I6,F7.2)
221 FORMAT (1X,I4,F7.1,F8.1,F7.1,45X,I6)
222 FORMAT (1X,I4,2F7.3,F8.1,F7.2,4F8.2,I6)
223 FORMAT (1X,I4,9F8.1)
224 FORMAT (1X,I4,F8.1,24X,F8.1,24X,F8.1)
225 FORMAT (1X,I4,7F9.2)
226 FORMAT (1X,I4,5F9.3)
227 FORMAT (1X,I4,5F9.1)
228 FORMAT (1X,I4,8F8.3,F9.4)
23 FORMAT (1X,'** INPUT EMISSIONS **')
231 FORMAT (4X,'CO2=EFOSS+NETDEF : SO2 EMISSIONS RELATIVE TO 1990')
24 FORMAT (1X,'** CARBON CYCLE DETAILS **')
C 25 FORMAT (1X,'FEEDBACKS ** TEMPERATURE * GPP :',F7.4,' * RESP :'
C +,F7.4,' * LITT OXDN :',F7.4,' ** FERTIL :',F7.4)
28 FORMAT (1X,F8.1,7F8.2,4f8.1)
C
30 FORMAT (1X,' ')
31 FORMAT (1X,'****************************************************')
47 FORMAT (1X,'** DECADAL CONTRIBUTIONS TO',
& ' GLOBAL RADIATIVE FORCING **')
48 FORMAT (1X,' (DELTA-Q in W/m**2 : PERCENTAGES IN BRACKETS : ',
& 'BASED ON END-OF-YEAR FORCING VALUES)')
C
50 FORMAT (1X,' INTERVAL')
51 FORMAT (1X,'START END CO2 CH4* ',
& 'N2O CFCs SO2 Total')
52 FORMAT (1x,2i5,2x,3(f6.3,' (',f5.0,')'),
& 2(f7.3,' (',f5.0,')'),f8.3)
53 FORMAT (1X,'STRAT H2O FROM CH4 : DELQH2O/DELQCH4 =',F6.3)
55 FORMAT (1X,'** GAS BY GAS DELTA-Q FROM 1990 **')
56 FORMAT (1X,'(BASED ON MID-YEAR FORCING VALUES : QEXTRA',
&' NOT INCLUDED)')
57 FORMAT (1X,'YEAR CO2 CH4tot N2O HALOtot',
& ' TROPO3 SO4DIR SO4IND BIOAER TOTAL CH4-O3')
571 FORMAT (1X,I4,8F8.4,F9.5,F8.4)
58 FORMAT (1X,'** GAS BY GAS DELTA-Q FROM 1765 **')
C
60 FORMAT (1X,'1990 DIRECT AEROSOL FORCING =',F6.3,'W/m**2',
&' (INCLUDES SOOT)')
61 FORMAT (1X,'1990 INDIRECT AEROSOL FORCING =',F6.3,'W/m**2')
62 FORMAT (1X,'1990 BIOMASS AEROSOL FORCING =',F6.3,'W/m**2')
63 FORMAT (1X,'1990 NON-CH4-RELATED TROP O3 FORCING =',F6.3,'W/m**2')
C
756 FORMAT (/1X,'NO EXTRA FORCING ADDED')
757 FORMAT (/1X,'EXTRA GLOBAL MEAN FORCING ADDED FROM',
&' QEXTRA.IN OVER',I5,' TO',I5,' INCLUSIVE')
758 FORMAT (/1X,'EXTRA HEMISPHERIC FORCINGS ADDED FROM',
&' QEXTRA.IN OVER',I5,' TO',I5,' INCLUSIVE')
759 FORMAT (/1X,'EXTRA NHO, NHL, SHO, SHL FORCINGS ADDED FROM',
&' QEXTRA.IN OVER',I5,' TO',I5,' INCLUSIVE')
760 FORMAT (2X,F10.3,' W/m**2 SUBTRACTED FROM ALL VALUES')
761 FORMAT (2X,'FORCING SCALED BY',F7.3,' AFTER OFFSET')
762 FORMAT (/1X,'QEXTRA FORCING USED ALONE')
C
800 FORMAT (1X,'LOW CONC CASE : NETDEF(80s) = 1.80GtC/yr',
&' : GIFFORD FERTILIZATION FACTOR =',F6.3)
801 FORMAT (1X,'MID CONC CASE : NETDEF(80s) = 1.10GtC/yr',
&' : GIFFORD FERTILIZATION FACTOR =',F6.3)
802 FORMAT (1X,'HIGH CONC CASE : NETDEF(80s) = 0.40GtC/yr',
&' : GIFFORD FERTILIZATION FACTOR =',F6.3)
803 FORMAT (1X,'USER CONC CASE : NETDEF(80s) =',F5.2,
&'GtC/yr : GIFFORD FERTILIZATION FACTOR =',F6.3)
804 FORMAT (/1X,'ALL CASES USE 1980s MEAN OCEAN FLUX OF 2.0GtC/yr')
805 FORMAT (1X,'DETAILED CARBON CYCLE OUTPUT IS FOR USER CASE ONLY')
806 FORMAT (1X,'METHANE OXIDATION TERM INCLUDED IN EMISSIONS')
807 FORMAT (1X,'METHANE OXIDATION TERM NOT INCLUDED IN EMISSIONS')
808 FORMAT (/1X,'*** ERROR : D80SIN SET TOO LOW AT',F7.3,
&' : RESET AT'F7.3,' ***')
809 FORMAT (/1X,'*** ERROR : D80SIN SET TOO HIGH AT',F7.3,
&' : RESET AT'F7.3,' ***')
810 FORMAT (/77X,'ENDYEAR')
811 FORMAT (/77X,'MIDYEAR')
812 FORMAT (1X,'YEAR ETOTAL EFOSS CH4OXN NETD GROSSD OFLUX',
&' ABFRAC PLANT C HLITT SOIL CONC DEL-M YEAR')
813 FORMAT (1X,I4,6F7.2,F7.3,F8.1,F6.1,F8.1,F8.2,F7.3,I6)
C
871 FORMAT (/1X,'CO2 CONC INPUT OVERWRITES CO2 EMISSIONS :',
&' Post-1990 CO2 forcing scaled up by',f5.1,'%')
872 FORMAT (/1X,'CO2 CONC INPUT OVERWRITES CO2 EMISSIONS :',
&' Other gas emissions as specified in GASEM.$$$')
873 FORMAT (/1X,'CO2 CONC INPUT OVERWRITES CO2 EMISSIONS :',
&' SO2 emissions as specified in GASEM.$$$')
874 FORMAT (/1X,'CO2 CONC INPUT OVERWRITES CO2 EMISSIONS :',
&' ESO2=ESO2(1990)+',F6.3,'*(EFOSS-EFOSS(1990))')
876 FORMAT (/1X,'Both CO2 conc scenario and compounded CO2 increase',
&' specified. Conc scenario overwritten.')
877 FORMAT (/1X,'Profile replaced by compounded CO2',
&' conc increase at',f5.2,'% per yr from mid1990')
C
900 FORMAT (1X,I5)
901 FORMAT (1X,2I5)
902 FORMAT (1X,I5,F10.0)
903 FORMAT (1X,I5,2F10.0)
904 FORMAT (1X,I5,4F10.0)
C
910 FORMAT (1X,'GSIC MODEL PARAMS : DTSTAR =',F6.3,
&' : DTEND =',F6.3,' : TAU =',F6.1,' : INIT VOL =',F6.1)
912 FORMAT (1X,'GREENLAND PARAMS : BETA =',F6.3)
913 FORMAT (1X,'ANTARCTIC PARAMS : BETA1 =',F6.3,
&' : BETA2=',F6.3,' : DB3 =',F6.3)
914 FORMAT (/1X,'DIFF/L SENSITIVITY CASE : RLO =',F6.3,' : XLAML =',
&F10.4,' : XLAMO =',F10.4)
915 FORMAT (/1X,'GLOBAL SENSITIVITY CASE : INITIAL XLAM =',F10.4)
916 FORMAT (1X,' ** WARNING, XLAML<0.0 : USE SMALLER XKLO **')
920 FORMAT (/1X,' LAND/OCEAN, NH/SH DQ SPLIT : (..) GIVES COMPONENT',
&' OF GLOBAL TOTAL')
921 FORMAT (2X,'YEAR QNHO QNHL',
&' QSHO QSHL QTOTAL')
922 FORMAT (1X,I5,4(F7.3,'(',F6.3,')'),F9.4)
930 FORMAT (/1X,'LOW CLIMATE MODEL PARAMETERS')
931 FORMAT (/1X,'MID CLIMATE MODEL PARAMETERS')
932 FORMAT (/1X,'HIGH CLIMATE MODEL PARAMETERS')
933 FORMAT (/1X,'USER CLIMATE MODEL PARAMETERS')
934 FORMAT (/2X,'YEAR GHG ESO21 ESO22 ESO23',
&' SUM')
935 FORMAT (1X,I5,5F10.4)
936 FORMAT (1X,'REF T',40X,F10.4)
937 format (1x,'PROFILE: ',A20)
C
4240 FORMAT (I10)
4241 FORMAT (F10.0)
4242 FORMAT (1X,I5,7F10.0)
C FOSSCO2,DEFCO2,CH4,N2O,SO2-1,SO2-2,SO2-3
4243 FORMAT (I2)
4244 FORMAT (I2,I5)
4444 FORMAT(6X,F10.0)
C
END
C
C********************************************************************
C
BLOCK DATA
c
parameter (iTp =740)
c
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
COMMON/CAR/EL1,EL2,EL3,TINV0(5),TINV(4,5),A(3,5),AA(4,5),
&BCO2(4),BTGPP,BTRESP,BTHUM,GAMP,GPP0,RESP0,QA0,U0,C0,B340(4),
&PHI,RG,TAUP,TAUH,TAUS,THP,THS,THH0,THS0,THPL,G1,G2,G3,FACTOR,
&EL21,EL32,XX1,XX2,XX3,XX4,XX5,XX6,DEE1,DEE2,DEE3,DEE4,DEE5,DEE6,
&FL1,FL2,FL3,XL,GAMH,GAMS,QS0,BTSOIL,FERTTYPE,TOTEM,CONVTERP,
&R(4),CPART(4,5),DELMASS(4,226:iTp),ABFRAC(4,226:iTp)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/COBS/COBS(0:226)
C
DATA FL(1)/0.420/,FL(2)/0.210/
C
DATA RHO/1.026/,SPECHT/0.9333/,HTCONS/4.1856/
C
C CO2 CONCS FROM END OF 1764 TO END OF 1990.
C UPDATED TO FEB 95 VERSION OF IPCC ON 3.11.95
C UPDATED TO MAR 95 VERSION OF IPCC ON 4.02.95
C UPDATED TO JUN 95 VERSION OF IPCC ON 8.08.95
C REMOVED FROM BLOCK DATA TO CO2HIST.IN ON 2.26.96
C
C INITIALISE CARBON CYCLE MODEL PARAMETERS.
C FIRST SPECIFY CUMULATIVE EMISSIONS TRANSITION POINTS.
C
DATA EL1/141./,EL2/565./,EL3/2500./
C
DATA DEE1/0.25/,DEE2/0.5/,DEE3/1.0/,DEE4/2.0/
&,DEE5/4.0/,DEE6/8.0/
C
C THESE ARE THE INVERSE DECAY TIMES AND MULTIPLYING CONSTANTS
C
DATA (TINV0(J),J=1,5)/0.0,0.0030303,0.0125,0.05,0.625/
DATA (A(1,J),J=1,5)/0.131,0.216,0.261,0.294,0.098/
DATA (A(2,J),J=1,5)/0.142,0.230,0.335,0.198,0.095/
DATA (A(3,J),J=1,5)/0.166,0.363,0.304,0.088,0.079/
C
end
C
C********************************************************************
C
SUBROUTINE INIT
C
parameter (iTp =740)
c
common /Limits/KEND
C
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
COMMON/CONCS/CH4(iTp),CN2O(iTp),C11(iTp),C12(iTp),C22(iTp),
&C134EQ(iTp),ECH4(226:iTp),EN2O(226:iTp),ECO(226:iTp),
&EVOC(226:iTp),ENOx(226:iTp),ESO2(0:iTp+1),ESO2SUM(226:iTp+1),
&ESO21(226:iTp+1),ESO22(226:iTp+1),ESO23(226:iTp+1)
C
COMMON/TANDSL/TEQU(iTp),TGAV(iTp),TNHO(iTp),
&TSHO(iTp),TNHL(iTp),TSHL(iTp),TDEEP(iTp),TNHAV(iTp),TSHAV(iTp),
&TLAND(iTp),TOCEAN(iTp),TOCN(40),TOCNPREV(40),IEXPAN,
&SIP,SGP,SAP,SLI(iTp),SLG(iTp),SLA(iTp),EX(0:iTp),SLT(iTp),
&QTOT(0:iTp),QGH(0:iTp),QOZ(0:iTp),QBIO(0:iTp),
&QSO2(0:iTp+1),QDIR(0:iTp+1)
C
common /radforc/qco2(0:iTp),qm(0:iTp),qn(0:iTp),qcfc(0:iTp)
COMMON /QOZONE/QCH4O3(0:iTp)
C
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,ENAT,es1990,
&dqdcoz,QCH4OZ
COMMON /VARW/Z(40),W(2),DW(2),TO0(2),TP0(2),WNH(iTp),WSH(iTp),
&TW0NH,TW0SH,IVARW
COMMON /ICEPARA/ZI0,NVALUES,MANYTAU,TAUFACT,ZI0BIT(100)
COMMON /ICEINIT/SIPBIT(100),SIBIT(100)
COMMON /AREAS/FNO,FNL,FSO,FSL
COMMON /TEMEXP/TEMEXP(2,40)
COMMON /NSIM/NSIM,NCLIM,ISCENGEN
C
C THE PROGRAM USES VARIOUS COUNTERS TO KEEP TRACK OF TIME.
C IC BEGINS WITH IC=1 TO IDENTIFY THE YEAR 1765. THUS IC
C =226 IS THE YEAR 1990, IC=336 IS THE YEAR 2100, ETC.
C CONC AND FORCING ARRAYS GIVE VALUES AT THE END OF THE
C YEAR. THUS CONC(1) IS THE VALUE AT THE END OF 1765, ETC.
C TIME (T) IS COUNTED FROM T=0 AT THE MIDDLE OF 1765, SO
C THAT THE MIDDLE OF 1990 IS T=225.0. TEMP AND SEALEVEL
C OUTPUT VALUES ARE AVERAGES OVER CALENDAR YEARS WITH THE
C VALUES BEING THOSE CALCULATED AT THE MIDPOINTS OF YEARS.
C TEMP(1) IS THEREFORE THE VALUE FOR THE MIDDLE OF 1765
C CORRESPONDING TO T=0.0 AND IC=1. EMISSIONS ARE TOTALS
C OVER CALENDAR YEARS, SO THE E(1) WOULD BE THE TOTAL FOR
C 1765, E(226) THE TOTAL FOR 1990, ETC.
C
C ****************************************************************
C
FNL=FL(1)/2.0
FNO=(1.0-FL(1))/2.0
FSL=FL(2)/2.0
FSO=(1.0-FL(2))/2.0
FLSUM=2.0*(FNL+FSL)
FOSUM=2.0*(FNO+FSO)
C
IP=0
IC=1
KC=1
T=0.0
DZ=100.
DZ1=DZ/2.
DTH=DT/HM
DTZ=DT/DZ
XXX=XK/DZ1
YYY=XK/DZ
AL=-DTZ*YYY
C
C INITIALIZE TO(I,L) = OCEAN TEMP CHANGE IN HEMISPHERE "I"
C AT LEVEL "L", AND TOCN(L) = AREA WEIGHTED MEAN OF HEMIS TEMPS.
C
DO L=1,40
TOCN(L)=0.0
DO I=1,2
TO(I,L)=0.0
END DO
END DO
C
Y(1)=0.0
Y(2)=0.0
Y(3)=0.0
Y(4)=0.0
HEM(1)=0.0
HEM(2)=0.0
C
C Z(I) DEPTH FROM BOTTOM OF MIXED LAYER FOR VARIABLE W
C
DO I=2,40
Z(I)=(I-2)*DZ+0.5*DZ
ENDDO
C
C DEFINE INITIAL TEMP PROFILE (TEM(I)) AND PRESSURE (P(I)).
C
ZED=HM/2.
P(1)=0.0098*(0.1005*ZED+10.5*(EXP(-ZED/3500.)-1.))
TEM(1)=20.98-13.12/HM-0.04025*HM
DO I=2,40
ZED=HM+(I-1)*DZ-DZ/2.
P(I)=0.0098*(0.1005*ZED+10.5*(EXP(-ZED/3500.)-1.))
ZZ=ZED/100.
IF(ZED.LE.130.)THEN
TEM(I)=18.98-4.025*ZZ
ELSE IF(ZED.LE.500.0)THEN
TEM(I)=17.01-2.829*ZZ+0.228*ZZ*ZZ
ELSE IF(ZED.LT.2500.)THEN
TEM(I)=EXP(EXP(1.007-0.0665*ZZ))
ELSE
TEM(I)=2.98-0.052*ZZ
ENDIF
END DO
C
C DEFINE THEORETICAL INITIAL TEMP PROFILE
C
TO0(1)=17.2 ! INITIAL MIXED LAYER TEMPE
TO0(2)=17.2 ! DITTO
TP0(1)=1.0 ! INITIAL TEMP OF POLAR SINKING WATER
TP0(2)=1.0 ! DITTO
C
DO I=1,2
DO L=2,40
TEMEXP(I,L)=TP0(I)+(TO0(I)-TP0(I))*EXP(-W0*Z(L)/XK)
END DO
END DO
C
C SET INITIAL VALUES FOR USE WITH VARIABLE W
C
W(1)=W0
W(2)=W0
DW(1)=0.0
DW(2)=0.0
C
C DEFINE INITIAL TEMP AND SEA LEVEL COMPONENTS.
C
TGAV(1)=0.0
TDEEP(1)=0.0
SLI(1)=0.0
SLG(1)=0.0
SLA(1)=0.0
SLT(1)=0.0
EX(0)=0.0
c
C SET ICE PARAMS AND INITIAL VALUES
C
MANYTAU=1
TAUFACT=0.3
DO N=1,NVALUES
ZI0BIT(N)=ZI0/NVALUES
SIBIT(N)=0.0
SIPBIT(N)=0.0
END DO
SIP=0.0
SGP=0.0
SAP=0.0
C
C CALCULATE NEW RADIATIVE FORCINGS EVERY FOURTH LOOP (I.E., WHEN
C NSIM=1,5,9,13,17) WHEN NEW SO2 EMISSIONS ARE USED.
C
IF(NCLIM.EQ.1.OR.ISCENGEN.EQ.9)THEN
CALL DELTAQ
C
C INITIALISE QTOT ETC AT START OF 1765.
C THIS ENSURES THAT ALL FORCINGS ARE ZERO AT THE MIDPOINT OF 1765.
C
QTOT(0) = -qtot(1)
qso2(0) = -qso2(1)
qdir(0) = -qdir(1)
qco2(0) = -qco2(1)
qm(0) = -qm(1)
qn(0) = -qn(1)
qcfc(0) = -qcfc(1)
QCH4O3(0)= -QCH4O3(1)
qgh(0) = -qgh(1)
QBIO(0) = -QBIO(1)
ENDIF
C
RETURN
END
C
C *******************************************************************
C
SUBROUTINE TSLCALC(N)
C
parameter (iTp =740)
C
common /Limits/KEND
C
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
COMMON/CONCS/CH4(iTp),CN2O(iTp),C11(iTp),C12(iTp),C22(iTp),
&C134EQ(iTp),ECH4(226:iTp),EN2O(226:iTp),ECO(226:iTp),
&EVOC(226:iTp),ENOx(226:iTp),ESO2(0:iTp+1),ESO2SUM(226:iTp+1),
&ESO21(226:iTp+1),ESO22(226:iTp+1),ESO23(226:iTp+1)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/TANDSL/TEQU(iTp),TGAV(iTp),TNHO(iTp),
&TSHO(iTp),TNHL(iTp),TSHL(iTp),TDEEP(iTp),TNHAV(iTp),TSHAV(iTp),
&TLAND(iTp),TOCEAN(iTp),TOCN(40),TOCNPREV(40),IEXPAN,
&SIP,SGP,SAP,SLI(iTp),SLG(iTp),SLA(iTp),EX(0:iTp),SLT(iTp),
&QTOT(0:iTp),QGH(0:iTp),QOZ(0:iTp),QBIO(0:iTp),
&QSO2(0:iTp+1),QDIR(0:iTp+1)
C
COMMON /VARW/Z(40),W(2),DW(2),TO0(2),TP0(2),WNH(iTp),WSH(iTp),
&TW0NH,TW0SH,IVARW
COMMON /QSPLIT/QNHO,QNHL,QSHO,QSHL,QGLOBE(0:iTp),
&QQNHO(0:iTp),QQNHL(0:iTp),QQSHO(0:iTp),QQSHL(0:iTp),
&QQQNHO(0:iTp),QQQNHL(0:iTp),QQQSHO(0:iTp),QQQSHL(0:iTp)
C
COMMON /ICE/TAUI,DTSTAR,DTEND,BETAG,BETA1,BETA2,DB3
COMMON /ICEPARA/ZI0,NVALUES,MANYTAU,TAUFACT,ZI0BIT(100)
COMMON /ICEINIT/SIPBIT(100),SIBIT(100)
C
QQNHO(N) = QNHO
QQNHL(N) = QNHL
QQSHO(N) = QSHO
QQSHL(N) = QSHL
TNHO(N) = Y(1)
TNHL(N) = Y(3)
TSHO(N) = Y(2)
TSHL(N) = Y(4)
TNHAV(N) = FO(1)*Y(1)+FL(1)*Y(3)
TSHAV(N) = FO(2)*Y(2)+FL(2)*Y(4)
TLAND(N) = (FL(1)*Y(3)+FL(2)*Y(4))/FLSUM
TOCEAN(N)= (FO(1)*Y(1)+FO(2)*Y(2))/FOSUM
TGAV(N) = (TNHAV(N)+TSHAV(N))/2.
C
C VARIABLE W TERMS
C
C IF(IVARW.EQ.0)THEN
C WNH(N)=W0
C WSH(N)=W0
C ENDIF
C
C IF(IVARW.EQ.1)THEN
C WTHRESH=0.1
C WNH(N)=W0-W0*Y(1)/TW0NH
C IF(WNH(N).LT.WTHRESH)THEN
C WNH(N)=WTHRESH
C ENDIF
C WSH(N)=W0-W0*Y(2)/TW0SH
C IF(WSH(N).LT.WTHRESH)THEN
C WSH(N)=WTHRESH
C ENDIF
C ENDIF
C
C CALCULATE MEAN OCEAN TEMP CHANGE PROFILE AND ITS INCREMENTAL
C CHANGE OVER ONE YEAR
C
DO L=1,40
TOCNPREV(L)=TOCN(L)
TOCN(L)=(FO(1)*TO(1,L)+FO(2)*TO(2,L))/FOSUM
END DO
TDEEP(N)=TOCN(40)
C
C THERMAL EXPANSION CONTRIBUTION TO SEA LEVEL CHANGE. FIRST PICK
C INCREMENTAL OR LUMP EXPANSION ALGORITHM. (FORMER PREFERED)
C
IF(IEXPAN.EQ.0)THEN
EXPAN=EX(N-1)
ELSE
EXPAN=0.0
ENDIF
C
ZLAYER=HM
DO I=1,40
C
IF(IEXPAN.EQ.0)THEN
DELTOC=TOCN(I)-TOCNPREV(I)
DELTBAR=(TOCN(I)+TOCNPREV(I))/2.0
ELSE
DELTOC=TOCN(I)
DELTBAR=TOCN(I)/2.0
ENDIF
C
TL1=TEM(I)+DELTBAR
TL2=TL1*TL1
TL3=TL2*TL1/6000.
PPL=P(I)
COEFE=52.55+28.051*PPL-0.7108*PPL*PPL+TL1*(12.9635-1.0833*PPL)-
+TL2*(0.1713-0.019263*PPL)+TL3*(10.41-1.1338*PPL)
DLR=COEFE*DELTOC*ZLAYER/10000.
ZLAYER=DZ
EXPAN=EXPAN+DLR
END DO
EX(N)=EXPAN
C
C ICE MELT CONTRIBUTIONS TO SEA LEVEL RISE.
C NSTEADY IS THE N VALUE FOR BEGINNING ICE MELT AT STEADY STATE
C
NSTEADY=116
IF(N.LE.NSTEADY)THEN
DO NNN=1,NVALUES
SIBIT(NNN)=0.0
END DO
SI=0.0
SG=0.0
SA=0.0
ELSE
DD=TGAV(N)-TGAV(NSTEADY)
CALL ICEMELT(DD,SIPBIT,SGP,SAP,SIBIT,SI,SG,SA)
ENDIF
C
C ABOVE GIVES END OF YEAR VALUES FOR ICE MELT. CALC MID YEAR
C VALUES FOR CONSISTENCY WITH TEMP AND EXPANSION AND THEN RESET
C FOR NEXT CALL TO SUBROUTINNE.
C
SLI(N)=(SIP+SI)/2.0
SLG(N)=(SGP+SG)/2.0
SLA(N)=(SAP+SA)/2.0
SLT(N)=SLI(N)+SLG(N)+SLA(N)+EX(N)
SIP=SI
SGP=SG
SAP=SA
DO NNN=1,NVALUES
SIPBIT(NNN)=SIBIT(NNN)
END DO
C
RETURN
END
C
C *******************************************************************
C
SUBROUTINE SPLIT(QGLOBE,A,BN,BS,QNO,QNL,QSO,QSL)
C
COMMON /AREAS/FNO,FNL,FSO,FSL
C
C Q VALUES ARE FORCINGS OVER AREAS IN W/M**2, F VALUES ARE THESE
C MULTIPLIED BY AREA FRACTIONS. THE SUM OF THE F MUST BE THE
C GLOBAL Q.
C
C FIRST SPLIT QGLOBE INTO NH AND SH
C
Q=QGLOBE
QS=2.*Q/(A+1.)
QN=A*QS
C
C NOW SPLIT NH AND SH INTO LAND AND OCEAN
C
FFN=2.0*(BN*FNL+FNO)
QNO=QN/FFN
QNL=BN*QNO
C
FFS=2.0*(BS*FSL+FSO)
QSO=QS/FFS
QSL=BS*QSO
C
RETURN
END
C
C *******************************************************************
C
SUBROUTINE RUNMOD
C
parameter (iTp =740)
C
common /Limits/KEND
C
DIMENSION AA(40),BB(40),A(40),B(40),C(40),D(40)
DIMENSION XLLGLOBE(2),XLLDIFF(2)
C
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
COMMON/CONCS/CH4(iTp),CN2O(iTp),C11(iTp),C12(iTp),C22(iTp),
&C134EQ(iTp),ECH4(226:iTp),EN2O(226:iTp),ECO(226:iTp),
&EVOC(226:iTp),ENOx(226:iTp),ESO2(0:iTp+1),ESO2SUM(226:iTp+1),
&ESO21(226:iTp+1),ESO22(226:iTp+1),ESO23(226:iTp+1)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/TANDSL/TEQU(iTp),TGAV(iTp),TNHO(iTp),
&TSHO(iTp),TNHL(iTp),TSHL(iTp),TDEEP(iTp),TNHAV(iTp),TSHAV(iTp),
&TLAND(iTp),TOCEAN(iTp),TOCN(40),TOCNPREV(40),IEXPAN,
&SIP,SGP,SAP,SLI(iTp),SLG(iTp),SLA(iTp),EX(0:iTp),SLT(iTp),
&QTOT(0:iTp),QGH(0:iTp),QOZ(0:iTp),QBIO(0:iTp),
&QSO2(0:iTp+1),QDIR(0:iTp+1)
C
COMMON /CO2READ/ICO2READ,XC(226:iTp),CO2SCALE
COMMON /COLDETC/ICOMP,COMPRATE,qtot86
COMMON /DSENS/IXLAM,XLAML,XLAMO
C
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,ENAT,es1990,
&dqdcoz,QCH4OZ
COMMON /VARW/Z(40),W(2),DW(2),TO0(2),TP0(2),WNH(iTp),WSH(iTp),
&TW0NH,TW0SH,IVARW
COMMON /QSPLIT/QNHO,QNHL,QSHO,QSHL,QGLOBE(0:iTp),
&QQNHO(0:iTp),QQNHL(0:iTp),QQSHO(0:iTp),QQSHL(0:iTp),
&QQQNHO(0:iTp),QQQNHL(0:iTp),QQQSHO(0:iTp),QQQSHL(0:iTp)
COMMON /AREAS/FNO,FNL,FSO,FSL
C
COMMON /QADD/IQREAD,JQFIRST,JQLAST,QEX(0:iTp),QEXNH(0:iTp),
&QEXSH(0:iTp),QEXNHO(0:iTp),QEXNHL(0:iTp),QEXSHO(0:iTp),
&QEXSHL(0:iTp)
C
COMMON /OLDTZ/IOLDTZ
COMMON /TEMEXP/TEMEXP(2,40)
COMMON /NSIM/NSIM,NCLIM,ISCENGEN
C
11 CONTINUE
C
C INCREMENT COUNTER AND ADD TIME STEP.
C T CORRESPONDS TO THE TIME AT WHICH NEW VALUES ARE TO
C BE CALCULATED
C
IP=IC
T=T+DT
C
IC=INT(T+1.49)
C
C PROTECT AGAINST CUMULATIVE ROUND-OFF ERROR IN T.
C (IF T IS VERY NEAR AN INTEGER, ROUND UP OR DOWN TO INTEGER.)
C
DIFF=ABS(T-INT(T))
IF(DIFF.LE.0.5)THEN
XT=FLOAT(INT(T))
ELSE
DIFF=1.0-DIFF
XT=FLOAT(INT(T)+1)
ENDIF
IF(DIFF.LT.0.01)T=XT
C
C AS SOON AS A NEW YEAR IS ENCOUNTERED (I.E. IC IS
C INCREMENTED), CALCULATE NEW END OF YEAR VALUES FOR
C CONCENTRATIONS AND RADIATIVE FORCING.
C THIS IS ONLY NECESSARY FOR FIRST PASS THROUGH NUMSIMS LOOP
C WHICH WILL SET ALL VALUES OF FORCING COMPONENT ARRAYS.
C
IF(NCLIM.EQ.1.OR.ISCENGEN.EQ.9)THEN
IF(IC.GT.IP) CALL DELTAQ
ENDIF
C
C INTERPOLATE FORCING FROM VALUES AT ENDS OF YEARS TO
C VALUE AT MIDPOINT OF TIME STEP.
C
C ***********************************************************
C ***********************************************************
C
C REVISED ALGORITHM FOR INTERPOLATION (15 APR, 1994)
C
C FIRST CALCULATE FRACTION OF YEAR THAT TIME CORRESPONDS TO AND
C IDENTIFY THE INTEGER (JC) THAT IS THE END OF THE YEAR IN
C WHICH THE TIME LIES.
C
T1=T-DT/2.0
JC=INT(T1+1.49)
FRAC=T1+0.5-INT(T1+0.5)
IF(FRAC.LT.0.001)FRAC=1.0
C
C CALCULATE GREENHOUSE GAS AND BIOMASS AEROSOL COMPONENTS OF
C FORCING AT START (0) AND END (1) OF YEAR, AND AT START OF
C PREVIOUS YEAR (P).
C
JPREV=0
IF(JC.GE.2)JPREV=JC-2
C
QGHP = QGH(JPREV)
QGH0 = QGH(JC-1)
QGH1 = QGH(JC)
QBIOP = QBIO(JPREV)
QBIO0 = QBIO(JC-1)
QBIO1 = QBIO(JC)
C
C RELABEL DIRECT AEROSOL AND OZONE FORCINGS
C
QDIRP = QDIR(JPREV)
QDIR0 = QDIR(JC-1)
QDIR1 = QDIR(JC)
QOZP = QOZ(JPREV)
QOZ0 = QOZ(JC-1)
QOZ1 = QOZ(JC)
C
C CALCULATE INDIRECT AEROSOL COMPONENT.
C
QINDP =QSO2(JPREV)-QDIRP
QIND0 =QSO2(JC-1)-QDIR0
QIND1 =QSO2(JC)-QDIR1
C
C ***********************************************************
C
C SPLIT AEROSOL & TROP O3 FORCING INTO LAND AND OCEAN IN NH AND SH
C
C A IS THE NH/SH FORCING RATIO
C BN IS THE LAND/OCEAN FORCING RATIO IN THE NH
C BS IS THE LAND/OCEAN FORCING RATIO IN THE SH
C
ADIR = 4.0
AIND = 2.0
AOZ =99.0
BNDIR= 9.0
BNIND= 9.0
BNOZ = 9.0
BSDIR= 9.0
BSIND= 9.0
BSOZ = 9.0
C
CALL SPLIT(QINDP,AIND,BNIND,BSIND,QINDNOP,QINDNLP,
&QINDSOP,QINDSLP)
CALL SPLIT(QIND0,AIND,BNIND,BSIND,QINDNO0,QINDNL0,
&QINDSO0,QINDSL0)
CALL SPLIT(QIND1,AIND,BNIND,BSIND,QINDNO1,QINDNL1,
&QINDSO1,QINDSL1)
C
CALL SPLIT(QDIRP,ADIR,BNDIR,BSDIR,QDIRNOP,QDIRNLP,
&QDIRSOP,QDIRSLP)
CALL SPLIT(QDIR0,ADIR,BNDIR,BSDIR,QDIRNO0,QDIRNL0,
&QDIRSO0,QDIRSL0)
CALL SPLIT(QDIR1,ADIR,BNDIR,BSDIR,QDIRNO1,QDIRNL1,
&QDIRSO1,QDIRSL1)
C
CALL SPLIT(QOZP,AOZ,BNOZ,BSOZ,QOZNOP,QOZNLP,
&QOZSOP,QOZSLP)
CALL SPLIT(QOZ0,AOZ,BNOZ,BSOZ,QOZNO0,QOZNL0,
&QOZSO0,QOZSL0)
CALL SPLIT(QOZ1,AOZ,BNOZ,BSOZ,QOZNO1,QOZNL1,
&QOZSO1,QOZSL1)
C
C CATCH QOZ AT 1990
C
IF((QOZ1.EQ.S90OZ).AND.(QOZ0.LT.S90OZ))THEN
OZBITNO=(QOZNO1-QOZNO0)/2.0
OZBITNL=(QOZNL1-QOZNL0)/2.0
OZBITSO=(QOZSO1-QOZSO0)/2.0
OZBITSL=(QOZSL1-QOZSL0)/2.0
ENDIF
C
C ***********************************************************
C
C NOW ADD
C
QSNHOP =QDIRNOP+QINDNOP+QOZNOP
QSNHLP =QDIRNLP+QINDNLP+QOZNLP
QSSHOP =QDIRSOP+QINDSOP+QOZSOP
QSSHLP =QDIRSLP+QINDSLP+QOZSLP
C
QSNHO0 =QDIRNO0+QINDNO0+QOZNO0
QSNHL0 =QDIRNL0+QINDNL0+QOZNL0
QSSHO0 =QDIRSO0+QINDSO0+QOZSO0
QSSHL0 =QDIRSL0+QINDSL0+QOZSL0
C
QSNHO1 =QDIRNO1+QINDNO1+QOZNO1
QSNHL1 =QDIRNL1+QINDNL1+QOZNL1
QSSHO1 =QDIRSO1+QINDSO1+QOZSO1
QSSHL1 =QDIRSL1+QINDSL1+QOZSL1
C
C ***********************************************************
C
C IF EXTRA FORCINGS ADDED THROUGH QEXTRA.IN, CALC CORRESP
C 'P', '0' AND '1' COMPONENTS FOR THESE. NOTE THAT THESE
C DATA ARE INPUT AS ANNUAL-MEAN VALUES, SO THEY MUST
C BE APPLIED EQUALLY AT THE START AND END OF THE YEAR,
C I.E., Q0=Q1=Q(JC).
C
QEXNHOP=0.0
QEXSHOP=0.0
QEXNHLP=0.0
QEXSHLP=0.0
QEXNHO0=0.0
QEXSHO0=0.0
QEXNHL0=0.0
QEXSHL0=0.0
QEXNHO1=0.0
QEXSHO1=0.0
QEXNHL1=0.0
QEXSHL1=0.0
C
IF(IQREAD.GE.1)THEN
IF((JC.GE.JQFIRST).AND.(JC.LE.JQLAST))THEN
QEXNHOP=QEXNHO(JC-1)
QEXSHOP=QEXSHO(JC-1)
QEXNHLP=QEXNHL(JC-1)
QEXSHLP=QEXSHL(JC-1)
QEXNHO0=QEXNHO(JC)
QEXSHO0=QEXSHO(JC)
QEXNHL0=QEXNHL(JC)
QEXSHL0=QEXSHL(JC)
QEXNHO1=QEXNHO(JC)
QEXSHO1=QEXSHO(JC)
QEXNHL1=QEXNHL(JC)
QEXSHL1=QEXSHL(JC)
ENDIF
ENDIF
C
C IF EXTRA NH AND SH FORCING INPUT THROUGH QEXTRA.IN, ADD
C TO AEROSOL FORCING
C
C IF IQREAD=2, USE EXTRA FORCING ONLY
C
IQR=1
IF(IQREAD.EQ.2)IQR=0
C
QSNHOP=IQR*QSNHOP+QEXNHOP
QSSHOP=IQR*QSSHOP+QEXSHOP
QSNHLP=IQR*QSNHLP+QEXNHLP
QSSHLP=IQR*QSSHLP+QEXSHLP
QGHP =IQR*QGHP
QBIOP =IQR*QBIOP
C
QSNHO0=IQR*QSNHO0+QEXNHO0
QSSHO0=IQR*QSSHO0+QEXSHO0
QSNHL0=IQR*QSNHL0+QEXNHL0
QSSHL0=IQR*QSSHL0+QEXSHL0
QGH0 =IQR*QGH0
QBIO0 =IQR*QBIO0
C
QSNHO1=IQR*QSNHO1+QEXNHO1
QSSHO1=IQR*QSSHO1+QEXSHO1
QSNHL1=IQR*QSNHL1+QEXNHL1
QSSHL1=IQR*QSSHL1+QEXSHL1
QGH1 =IQR*QGH1
QBIO1 =IQR*QBIO1
C
C CALCULATE FORCING COMPONENTS FOR MIDPOINT OF INTERVAL.
C
QSNHLM =(QSNHL0+QSNHL1)/2.0
QSSHLM =(QSSHL0+QSSHL1)/2.0
QSNHOM =(QSNHO0+QSNHO1)/2.0
QSSHOM =(QSSHO0+QSSHO1)/2.0
QGHM =(QGH0 +QGH1 )/2.0
QBIOM =(QBIO0 +QBIO1 )/2.0
C
C CORRECT FOR NONLINEAR CHANGE IN QOZ AT 1990
C
IF((QOZ1.EQ.S90OZ).AND.(QOZ0.LT.S90OZ))THEN
QSNHOM = QSNHOM+OZBITNO
QSNHLM = QSNHLM+OZBITNL
QSSHOM = QSSHOM+OZBITSO
QSSHLM = QSSHLM+OZBITSL
ENDIF
C
C ************************************************************
C
C OVERWRITE MIDPOINT VALUES FOR SPECIAL CASE OF ICOMP=1 FOR 1990
C
IF(JC.EQ.226)THEN
C
IF(ICOMP.EQ.1)THEN
QSNHLM =QSNHL0+(QSNHL0-QSNHLP)/2.0
QSSHLM =QSSHL0+(QSSHL0-QSSHLP)/2.0
QSNHOM =QSNHO0+(QSNHO0-QSNHOP)/2.0
QSSHOM =QSSHO0+(QSSHO0-QSSHOP)/2.0
QGHM =QGH0 +(QGH0 -QGHP )/2.0
QBIOM =QBIO0 +(QBIO0 -QBIOP )/2.0
ENDIF
ENDIF
C
C *************************************************************
C
C PUT TOTAL FORCINGS AT MIDPOINT OF YEAR JC INTO ARRAYS. NOTE
C THAT THIS GIVES CORRECT 1990 RESULTS EVEN IF ICOMP = 1.
C
QQQNHL(JC)=QSNHLM+QGHM+QBIOM
QQQNHO(JC)=QSNHOM+QGHM+QBIOM
QQQSHL(JC)=QSSHLM+QGHM+QBIOM
QQQSHO(JC)=QSSHOM+QGHM+QBIOM
FNHL=QQQNHL(JC)*FNL
FNHO=QQQNHO(JC)*FNO
FSHL=QQQSHL(JC)*FSL
FSHO=QQQSHO(JC)*FSO
QGLOBE(JC)=FNHL+FNHO+FSHL+FSHO
TEQU(JC)=TE*QGLOBE(JC)/Q2X
C
C CALCULATE FORCING INCREMENTS OVER YEAR IN WHICH TIME STEP LIES.
C DONE FOR 1ST AND 2ND HALVES OF YEAR SEPARATELY.
C
DSNHL0M =(QSNHLM-QSNHL0)*2.0
DSSHL0M =(QSSHLM-QSSHL0)*2.0
DSNHO0M =(QSNHOM-QSNHO0)*2.0
DSSHO0M =(QSSHOM-QSSHO0)*2.0
DGH0M =(QGHM -QGH0 )*2.0
DBIO0M =(QBIOM -QBIO0 )*2.0
C
DSNHLM1 =(QSNHL1-QSNHLM)*2.0
DSSHLM1 =(QSSHL1-QSSHLM)*2.0
DSNHOM1 =(QSNHO1-QSNHOM)*2.0
DSSHOM1 =(QSSHO1-QSSHOM)*2.0
DGHM1 =(QGH1 -QGHM )*2.0
DBIOM1 =(QBIO1 -QBIOM )*2.0
C
C NOW CALCULATE THE FORCING VALUES AT THE MIDPOINT OF THE TIME STEP.
C
IF(FRAC.LE.0.5)THEN
QSNHL=QSNHL0+FRAC*DSNHL0M
QSSHL=QSSHL0+FRAC*DSSHL0M
QSNHO=QSNHO0+FRAC*DSNHO0M
QSSHO=QSSHO0+FRAC*DSSHO0M
QGHG =QGH0 +FRAC*DGH0M
QBIOG=QBIO0 +FRAC*DBIO0M
ELSE
QSNHL=QSNHLM+(FRAC-0.5)*DSNHLM1
QSSHL=QSSHLM+(FRAC-0.5)*DSSHLM1
QSNHO=QSNHOM+(FRAC-0.5)*DSNHOM1
QSSHO=QSSHOM+(FRAC-0.5)*DSSHOM1
QGHG =QGHM +(FRAC-0.5)*DGHM1
QBIOG=QBIOM +(FRAC-0.5)*DBIOM1
ENDIF
C
C COMBINE GREENHOUSE, (SO4 AEROSOL + TROP O3 + QEXTRA) AND BIO
C AEROSOL FORCINGS
C
QNHL=QGHG+QSNHL+QBIOG
QNHO=QGHG+QSNHO+QBIOG
QSHL=QGHG+QSSHL+QBIOG
QSHO=QGHG+QSSHO+QBIOG
C
C **********************************************************
C
DO 10 II=1,2
C
C ***** START OF DIFFERENTIAL SENSITIVITY TERMS *****
C
IF(IXLAM.EQ.1)THEN
WWW=FO(II)*(XKLO+FL(II)*XLAML)
XLLDIFF(II)=XLAMO+XLAML*XKLO*FL(II)/WWW
XLL=XLLDIFF(II)
ELSE
WWW=FO(II)*(XKLO+FL(II)*XLAM)
XLLGLOBE(II)=XLAM*(1.+XKLO*FL(II)/WWW)
XLL=XLLGLOBE(II)
ENDIF
C
C ***** END OF DIFFERENTIAL SENSITIVITY TERMS *****
C
CL=-DTZ*(YYY+W(II))
BL=1.-AL-CL
C
A(1)=1.+DTH*(XXX+W(II)*PI+XLL/FK)
B(1)=-DTH*(XXX+W(II))
A(2)=-DTZ*XXX
C(2)=CL
B(2)=1.-A(2)-C(2)
D(2)=TO(II,2)
C
C TERMS FOR VARIABLE W
C
DTZW=DTZ*DW(II)
C
IF(IVARW.GE.1)THEN
IF(IOLDTZ.EQ.1)THEN
D(2)=D(2)+DTZW*(TEMEXP(II,3)-TEMEXP(II,2))
ELSE
D(2)=D(2)+DTZW*(TEM(3)-TEM(2))
ENDIF
ENDIF
C
DO L=3,39
A(L)=AL
B(L)=BL
C(L)=CL
D(L)=TO(II,L)
C
C TERMS FOR VARIABLE W
C
IF(IVARW.GE.1)THEN
IF(IOLDTZ.EQ.1)THEN
D(L)=D(L)+DTZW*(TEMEXP(II,L+1)-TEMEXP(II,L))
ELSE
D(L)=D(L)+DTZW*(TEM(L+1)-TEM(L))
ENDIF
ENDIF
C
END DO
C
A(40)=AL
B(40)=1.-CL
D(40)=TO(II,40)+TO(II,1)*PI*DTZ*W(II)
C
C TERMS FOR VARIABLE W
C
IF(IVARW.GE.1)THEN
IF(IOLDTZ.EQ.1)THEN
D(40)=D(40)+DTZW*(TP0(II)-TEMEXP(II,40))
ELSE
D(40)=D(40)+DTZW*(TP0(II)-TEM(40))
ENDIF
ENDIF
C
C FORL is the land forcing term
C
if(ii.eq.1) then
forl = qnhl*xklo*fl(ii)/www
heat = (qnho+hem(ii)+forl)*dth/fk
else
forl = qshl*xklo*fl(ii)/www
heat = (qsho+hem(ii)+forl)*dth/fk
endif
C
D(1)=TO(II,1)+HEAT
C
C TERMS FOR VARIABLE W
C
IF(IVARW.GE.1)THEN
IF(IOLDTZ.EQ.1)THEN
D(1)=D(1)+DTH*DW(II)*(TEMEXP(II,2)-TP0(II))
ELSE
D(1)=D(1)+DTH*DW(II)*(TEM(2)-TP0(II))
ENDIF
ENDIF
C
C THIS IS THE OLD BAND SUBROUTINE
C
AA(1)=-B(1)/A(1)
BB(1)=D(1)/A(1)
DO L=2,39
VV=A(L)*AA(L-1)+B(L)
AA(L)=-C(L)/VV
BB(L)=(D(L)-A(L)*BB(L-1))/VV
END DO
TO(II,40)=(D(40)-A(40)*BB(39))/(A(40)*AA(39)+B(40))
DO I=1,39
L=40-I
TO(II,L)=AA(L)*TO(II,L+1)+BB(L)
END DO
C
10 CONTINUE
C
C Y(1,2,3,4) ARE NH OCEAN, SH OCEAN, NH LAND & SH LAND TEMPS.
C
Y(1)=TO(1,1)
Y(2)=TO(2,1)
C
C ***** START OF DIFFERENTIAL SENSITIVITY TERMS *****
C
IF(IXLAM.EQ.1)THEN
Y(3)=(FL(1)*qnhl+XKLO*Y(1))/(FL(1)*XLAML+XKLO)
Y(4)=(FL(2)*qshl+XKLO*Y(2))/(FL(2)*XLAML+XKLO)
ELSE
Y(3)=(FL(1)*qnhl+XKLO*Y(1))/(FL(1)*XLAM+XKLO)
Y(4)=(FL(2)*qshl+XKLO*Y(2))/(FL(2)*XLAM+XKLO)
ENDIF
C
C ***** END OF DIFFERENTIAL SENSITIVITY TERMS *****
C
HEM(1)=(XKNS/FO(1))*(Y(2)-Y(1))
HEM(2)=(XKNS/FO(2))*(Y(1)-Y(2))
C
C VARIABLE W TERMS
C
IF(IVARW.EQ.0)THEN
W(1)=W0
DW(1)=0.0
W(2)=W0
DW(2)=0.0
ENDIF
C
IF(IVARW.GE.1)THEN
WTHRESH=0.1
OCEANT=(FO(1)*Y(1)+FO(2)*Y(2))/(FO(1)+FO(2))
C
DW(1)=-W0*Y(1)/TW0NH
IF(IVARW.EQ.2)DW(1)=-W0+WNH(JC)
W(1)=W0+DW(1)
IF(IVARW.EQ.1.AND.W(1).LT.WTHRESH)THEN
W(1)=WTHRESH
DW(1)=WTHRESH-W0
ENDIF
C
DW(2)=-W0*Y(2)/TW0SH
IF(IVARW.EQ.2)DW(2)=-W0+WSH(JC)
W(2)=W0+DW(2)
IF(IVARW.EQ.1.AND.W(2).LT.WTHRESH)THEN
W(2)=WTHRESH
DW(2)=WTHRESH-W0
ENDIF
ENDIF
C
C IF(IVARW.LE.1)THEN
WNH(JC)=W(1)
WSH(JC)=W(2)
C ENDIF
C
C HAVING CALCULATED VALUES AT NEW TIME 'T', DECIDE WHETHER OR NOT
C TEMPS ARE ANNUAL VALUES (I.E., CORRESPOND TO T=MIDPOINT OF
C YEAR). IF SO, GO TO TSLCALC, CALCULATE GLOBAL MEAN TEMP (TGAV)
C AND INSERT INTO ARRAY, AND CALCULATE SEA LEVEL RISE COMPONENTS.
C NOTE THAT, BECAUSE FORCING VALUES ARE GIVEN AT ENDS OF YEARS
C AND TIMES ARE INTEGERS AT MIDPOINTS OF YEARS, A ONE YEAR TIME
C STEP WILL STILL GIVE MID-YEAR VALUES FOR CALCULATED TEMPERATURES.
C
KP=KC
KC=INT(T+1.01)
IF(KC.GT.KP)CALL TSLCALC(KC)
C
IF(T.GE.TEND)RETURN
GO TO 11
END
C
C *******************************************************************
C
SUBROUTINE DELTAQ
C
parameter (iTp =740)
C
common /Limits/KEND
C
COMMON/CLIM/IC,IP,KC,DT,DZ,FK,HM,Q2X,QXX,PI,T,TE,TEND,W0,XK,XKLO,
+XKNS,XLAM,FL(2),FO(2),FLSUM,FOSUM,HEM(2),P(40),TEM(40),TO(2,40),
+AL,BL,CL,DTH,DTZ,DZ1,XLL,WWW,XXX,YYY,RHO,SPECHT,HTCONS,Y(4)
C
COMMON/CONCS/CH4(iTp),CN2O(iTp),C11(iTp),C12(iTp),C22(iTp),
&C134EQ(iTp),ECH4(226:iTp),EN2O(226:iTp),ECO(226:iTp),
&EVOC(226:iTp),ENOx(226:iTp),ESO2(0:iTp+1),ESO2SUM(226:iTp+1),
&ESO21(226:iTp+1),ESO22(226:iTp+1),ESO23(226:iTp+1)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/TANDSL/TEQU(iTp),TGAV(iTp),TNHO(iTp),
&TSHO(iTp),TNHL(iTp),TSHL(iTp),TDEEP(iTp),TNHAV(iTp),TSHAV(iTp),
&TLAND(iTp),TOCEAN(iTp),TOCN(40),TOCNPREV(40),IEXPAN,
&SIP,SGP,SAP,SLI(iTp),SLG(iTp),SLA(iTp),EX(0:iTp),SLT(iTp),
&QTOT(0:iTp),QGH(0:iTp),QOZ(0:iTp),QBIO(0:iTp),
&QSO2(0:iTp+1),QDIR(0:iTp+1)
C
COMMON/CAR/EL1,EL2,EL3,TINV0(5),TINV(4,5),A(3,5),AA(4,5),
&BCO2(4),BTGPP,BTRESP,BTHUM,GAMP,GPP0,RESP0,QA0,U0,C0,B340(4),
&PHI,RG,TAUP,TAUH,TAUS,THP,THS,THH0,THS0,THPL,G1,G2,G3,FACTOR,
&EL21,EL32,XX1,XX2,XX3,XX4,XX5,XX6,DEE1,DEE2,DEE3,DEE4,DEE5,DEE6,
&FL1,FL2,FL3,XL,GAMH,GAMS,QS0,BTSOIL,FERTTYPE,TOTEM,CONVTERP,
&R(4),CPART(4,5),DELMASS(4,226:iTp),ABFRAC(4,226:iTp)
c
common /meth/emeth(226:iTp),imeth,ch4l(225:iTp),ch4b(225:iTp),
+ch4h(225:iTp),ef4(226:iTp)
c
common /radforc/qco2(0:iTp),qm(0:iTp),qn(0:iTp),qcfc(0:iTp)
COMMON /QOZONE/QCH4O3(0:iTp)
c
common /StratH2O/StratH2O
c
common /TauMeth/ITAUMETH,TAU90CH4,TauSoil,DELTAU,
+ TAUINIT,ch4bar90,eno90,eco90,evo90
COMMON /TCH4/TCH4(iTp)
c
common /O3feedback/iO3feed
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,ENAT,es1990,
&dqdcoz,QCH4OZ
COMMON /IGHG/IGHG
COMMON /CO2READ/ICO2READ,XC(226:iTp),CO2SCALE
COMMON /COLDETC/ICOMP,COMPRATE,qtot86
COMMON /QCON/IQCONST,QCONST,RAMPDQDT
common /netdef/ednet(226:iTp),DUSER,FUSER
COMMON /CH4CORR/CORRUSER,CORRMHI,CORRMMID,CORRMLO
C
C THIS SUBROUTINE IS ONLY ENTERED WHEN THE IC YEAR COUNT
C INCREMENTS. CONCENTRATIONS AND RADIATIVE FORCINGS ARE
C THEN CALCULATED FOR THE END OF THE IC YEAR. SINCE THE
C IC INCREMENT MAY BE GREATER THAN 1, A LOOP IS NEEDED TO
C FILL IN THE MISSING IC VALUES.
C
DO 10 J=IP+1,IC
C
C *******************************************************
C *******************************************************
C
C START OF LONG IF STATEMENT
C
IF(J.LE.226)THEN
C
CALL HISTORY(J,CO2(J),CH4(J),CN2O(J),C11(J),C12(J),
+ C22(J),C134EQ(J),eso2(J),eso2(j+1))
c
if(j.eq.226) then
ch4l(225) = ch4(225)
ch4b(225) = ch4(225)
ch4h(225) = ch4(225)
ch4l(226) = ch4(226)
ch4b(226) = ch4(226)
ch4h(226) = ch4(226)
c
c Calculate additional CO2 emissions due to CH4 oxidation for 1990
c First specify fossil fuel fraction of total CH4 emission
c
fffrac = 0.2
emeth(226) = fffrac*0.0020625*(ch4bar90-790.)/TAUINIT
if(imeth.eq.1) then
ef4(226) = ef(226)+emeth(226)
else
ef4(226) = ef(226)
endif
endif
C
ELSE
C
C FOR CONCS BEYOND THE END OF 1990, CALL THE VARIOUS EMISSIONS
C TO CONCENTRATIONS MODELS. NOTE THAT THE INPUT EMISSIONS
C ARRAYS ARE NUMBERED AS FOR OTHER VARIABLES ; I.E. E(226)
C IS THE 1990 VALUE. EMISSIONS VALUES UP TO AND INCLUDING E(225)
C ARE NOT USED.
C
C FOR CALCULATING FEEDBACKS IN CARBON AND METHANE, NEED TO USE
C EXTRAPOLATED VALUES FOR TEMPERATURE AND CO2 CONCENTRATION.
C RELEVANT VALUES ARE FROM 1765/1770. FOR TEMP, THE MODELLED
C CHANGE IS USED TO 1990 FOR CONSISTENCY WITH POST-1990 TEMPS.
C FOR METHANE THIS TEMP IS CALCULATED BELOW. FOR THE CARBON CYCLE
C MODEL, THE EXTRAPOLATION IS DONE IN SUBROUTINE CARBON.
C
TX=2.0*TGAV(J-1)-TGAV(J-2)
IF(J.EQ.226)DELT90=TX
if(j.eq.227) then
t90lo=TAUINIT-DELTAU
t90mid=TAUINIT
t90hi=TAUINIT+DELTAU
if(ITAUMETH.eq.1)t90user=t90lo
if(ITAUMETH.eq.2)t90user=t90mid
if(ITAUMETH.eq.3)t90user=t90hi
if(ITAUMETH.eq.4)t90user=TAU90CH4
tlm2 = t90lo
tlm1 = t90lo
tbm2 = t90mid
tbm1 = t90mid
thm2 = t90hi
thm1 = t90hi
tm2 = t90user
tm1 = t90user
endif
C
C *******************************************************
C
C FOR METHANE CONC PROJECTIONS, NEED TO USE EMISSIONS CONSISTENT WITH
C THE LIFETIME. THIS IS DONE USING CORRECTION FACTORS CALCULATED
C IN THE MAIN PROGRAM. NOTE THAT THE INPUT EMISSIONS HAVE ALREADY
C BEEN OFFSET BY THE AMOUNT APPROPRIATE TO THE USER-SPECIFIED
C LIFETIME (I.E., BY CORRUSER).
C THIS WAS CORRECTED ON 97/12/13.
C
C LOW LIFETIME
C
EECH4 = ECH4(J)-CORRUSER+CORRMLO

CALL METHANE(1,CH4L(J-2),CH4L(J-1),eech4,ENOx(J),ECO(J),
& EVOc(J),tlm2,tlm1,tx,CH4L(J),TauLo)
tlm2 = tlm1
tlm1 = TauLo
C
C MID (BEST) LIFETIME
C
EECH4 = ECH4(J)-CORRUSER+CORRMMID

CALL METHANE(2,CH4B(J-2),CH4B(J-1),eech4,ENOx(J),ECO(J),
& EVOc(J),tbm2,tbm1,tx,CH4B(J),TauBest)
tbm2 = tbm1
tbm1 = TauBest
C
C HIGH LIFETIME
C
EECH4 = ECH4(J)-CORRUSER+CORRMHI

CALL METHANE(3,CH4H(J-2),CH4H(J-1),eech4,ENOx(J),ECO(J),
& EVOc(J),thm2,thm1,tx,CH4H(J),TauHi)
thm2 = thm1
thm1 = TauHi
C
C USER LIFETIME (ONE OF ABOVE, OR CONSTANT AT SPECIFIED 1990 VALUE)
C
EECH4 = ECH4(J)
C
CALL METHANE(ITAUMETH,CH4(J-2),CH4(J-1),eech4,ENOx(J),ECO(J),
& EVOc(J),tm2,tm1,tx,CH4(J),TauCH4)
tm2 = tm1
tm1 = TauCH4
C
C SAVE USER-MODEL METHANE LIFETIME. TCH4(J) = CHEMICAL (OH)
C LIFETIME. THIS IS THE SAME AS ......
C TCH4EFF(J)=CH4BAR/(ECH4(J)/2.75-DELCH4-SOILSINK)
C
TCH4(J)=TauCH4
C
C Methane oxidation source: based on user methane projection only.
C User projection determined by choice of ITAUMETH (1,2,3,or 4).
C
CH4BAR=(CH4(J-1)+CH4(J))/2
EMETH(J) = FFFRAC*0.0020625*(CH4BAR-790.)/TAUCH4
IF(IMETH.EQ.1)THEN
EF4(J) = EF(J)+EMETH(J)
ELSE
EF4(J) = EF(J)
ENDIF
C
C *******************************************************
C
C CARBON CYCLE MODEL CALL
C NC=1, BCO2 BASED ON D80S=1.8 TO GIVE LOWER BOUND CONCS
C NC=2, BCO2 BASED ON D80S=1.1 TO GIVE BEST GUESS CONCS
C NC=3, BCO2 BASED ON D80S=0.4 TO GIVE UPPER BOUND CONCS
C NC=4, BCO2 BASED ON USER SELECTED D80S
C ALL CASES USE F80S =2.0
C
DO NC = 1,4
C
C CALL INITCAR TO INITIALIZE CARBON CYCLE MODEL
C
IF(J.EQ.227)THEN
FIN=2.0
DIN=2.5-0.7*NC
IF(NC.EQ.4)THEN
FIN=FUSER
DIN=DUSER
ENDIF
CALL INITCAR(NC,DIN,FIN)
ENDIF
C
C IF IDRELAX.NE.O, OVERWRITE EDNET(J) FOR 1990 TO 1990+DRELAX WITH
C A LINEAR INTERP BETWEEN THE 1990 VALUE BASED ON BALANCING THE
C 1980S-MEAN CARBON BUDGET TO THE APPROPRIATE VALUE OBTAINED FROM
C THE EMISSIONS INPUT FILE.
C
IDRELAX=10
IF(IDRELAX.NE.0)THEN
EDNET(226)=EDNET90(NC)
JDEND=226+IDRELAX
DELED=(EDNET(JDEND)-EDNET(226))/FLOAT(IDRELAX)
DO JD=227,JDEND-1
EDNET(JD)=EDNET(226)+DELED*(JD-226)
END DO
ENDIF
C
c Note: for temp feedback on CO2, default or user temp is used in all
c four nc cases. Strictly in (eg) the upper bound case should
c use the corresponding temp. However the upper bound CO2 is not
c generally passed to the climate model so this temp is not
c calculated. The error in making this approx must be small
c as it is only a second order effect.
c Note: this also applies to the methane model.
C
TEMP=TGAV(J)+0.6
CMINUS3=CCO2(NC,J-3)
CMINUS2=CCO2(NC,J-2)
CMINUS1=CCO2(NC,J-1)
CALL CARBON(NC,TEMP,EF4(J),EDNET(J),CMINUS3,CMINUS2,CMINUS1,
& PL(NC,J-1),HL(NC,J-1),SOIL(NC,J-1),REGROW(NC,J-1),ETOT(NC,J-1),
& PL(NC,J) ,HL(NC,J) ,SOIL(NC,J) ,REGROW(NC,J) ,ETOT(NC,J) ,
& ESUM(J),FOC(NC,J),DELMASS(NC,J),EDGROSS(NC,J),CCO2(NC,J))
C
end do
co2(J) = cco2(4,j)
C
C OVERWRITE CARBON CYCLE CO2 CONCS FOR YEARS.GE.1990 WITH INPUT
C DATA FROM CO2INPUT.DAT IF ICO2READ.GE.1. THIS INPUT IS
C OVERWRITTEN BY COMPOUND CO2 CONC INCREASE IF ICOMP=1
C
IF(ICO2READ.GE.1) co2(J)=xc(J)
C IF(ICOMP.EQ.1)THEN
C IF(J.EQ.225)CO289=CO2(J)
C IF(J.EQ.226)THEN
C CO290=CO2(J)
C CO2M =(CO289+CO290)/2.0
C ENDIF
C IF(J.GE.226)THEN
C CO2(J)=CO2M*EXP((J-225.5)*COMPRATE/100.)
C ENDIF
C ENDIF
C
C N2O CONCs
C
CALL NITROUS(CN2O(J-1),EN2O(J),CN2O(J))
ENDIF
C
C **************************************************
C
C END OF LONG IF STATEMNT
C
C *******************************************************
C *******************************************************
C
C HALOCARBON FORCING
C
IF(J.LE.186)QCFC(J)=0.0
IF(J.LE.215.AND.J.GT.186)QCFC(J)=0.0884*(J-186)/29.
IF(J.LE.226.AND.J.GT.215)QCFC(J)=0.0884+0.0233*(J-215)/11.
IF(J.GT.226)THEN
XJ=(J-226)/10.
ZJ=XJ-5.5
IF(J.LT.266)DQCFC=XJ*.032
IF(J.GE.266.AND.J.LT.281)DQCFC=0.128+(XJ-4.)*0.04133
IF(J.GT.281)DQCFC=0.190+0.043*ZJ-0.003719*ZJ*ZJ
QCFC(J)=0.1117+DQCFC
IF(QCFC(J).LT.0.0)QCFC(J)=0.0
ENDIF
C
C METHANE FORCING
C
QCH4=0.036*(SQRT(CH4(J))-SQRT(CH4(1)))
C
C QH2O IS THE ADDITIONAL INDIRECT CH4 FORCING DUE TO PRODUCTION
C OF STRAT H2O FORM CH4 OXIDATION. QCH4OZ IS THE ENHANCEMENT
C OF QCH4 DUE TO CH4-INDUCED OZONE PRODUCTION.
C THE DENOMINATOR IN QCH4OZ HAS BEEN CHOSEN TO GIVE 0.08W/m**2
C FORCING IN MID 1990. IT DIFFERS SLIGHTLY FROM THE VALUE
C ORIGINALLY ADVISED BY PRATHER (WHICH WAS 353). THE 0.0025
C CORRECTION TERM IS TO MAKE THE CHANGE RELATIVE TO MID 1765.
C
QH2O=STRATH2O*QCH4
QCH4OZ=DQDCOZ*1.5*(CH4(J)-(CH4(1)-.0025))/347.62
C
XM=CH4(J)/1000.
WW=CN2O(1)/1000.
ZZ=XM*WW
AB=0.636*ZZ**0.75+0.007*XM*ZZ**1.52
F=0.47*ALOG(1.+AB)
XM0=CH4(1)/1000.
ZZ0=XM0*WW
AB0=0.636*ZZ0**0.75+0.007*XM0*ZZ0**1.52
F0=0.47*ALOG(1.+AB0)
QLAP=-F+F0
QMeth=QCH4+QLAP
QM(J) = qMeth+qH2O+QCH4OZ
QCH4O3(J)=QCH4OZ
c
C NITROUS OXIDE FORCING
c
QN2O=0.14*(SQRT(CN2O(J))-SQRT(CN2O(1)))
c
XM=XM0
WW=CN2O(J)/1000.
ZZ=XM*WW
AB=0.636*ZZ**0.75+0.007*XM*ZZ**1.52
F=0.47*ALOG(1.+AB)
QLAP=-F+F0
QN(J)=QN2O+QLAP
c
C TOTAL GREENHOUSE FORCING
C
QCO2(J)=QXX*ALOG(CO2(J)/CO2(1))
QGH(J)=QCO2(J)+QM(J)+QN(J)+QCFC(J)
C
C ADD IN BIOMASS BURNING TERM
C
IF(J.LE.136)THEN
QBIO(J)=0.0
ELSE IF(J.LE.211)THEN
QBIO(J)=-0.04*(J-136)/75.
ELSE IF(J.LE.221)THEN
QBIO(J)=-0.04-(J-211)*0.012
ELSE
QBIO(J)=-0.16
ENDIF
QBIO(J)=QBIO(J)*(-S90BIO/0.16)
C
IF(IGHG.EQ.0)QGH(J)=0.0
IF(IGHG.EQ.1)QGH(J)=QCO2(J)
C
C GLOBAL SULPHATE AEROSOL FORCING
C
CALL SULPHATE(J,ESO2(J),ESO2(J+1),QSO2(J),QDIR(J),QOZ(J))
C
C TOTAL GLOBAL FORCING INCLUDING AEROSOLS.
C NOTE THAT AEROSOL FORCING IS RESET IN SULPHATE SUBROUTINE
C IF ICOMP=1 OR ICO2READ=1 SINCE THIS IS USED IN
C RUNMOD TO CALCULATE SPATIALLY DISAGGREGATED FORCING.
C 950530 : JUST NOTICED THAT THIS IS NO LONGER VALID
C
qtot(J) = QGH(J) + qso2(J) + qoz(J) + QBIO(J)
if(j.eq.86)qtot86=qtot(J)
C
C ***************************************************************
C
C SWITCH TO GIVE CONSTANT FORCING AT QCONST IF IQCONST=1 OR
C RAMP UP TO QCONST IF IQCONST=2 (RATE = RAMPDQDT (W/m**2/YEAR)).
C
IF(IQCONST.GE.1)THEN
QSO2(J)=0.0
QOZ(J)=0.0
QBIO(J)=0.0
QDIR(J)=0.0
ENDIF
IF(IQCONST.EQ.1)THEN
QTOT(J)=QCONST
QGH(J)=QCONST
ENDIF
IF(IQCONST.EQ.2)THEN
QTOT(J)=RAMPDQDT*(J-1)
QGH(J)=RAMPDQDT*(J-1)
IF(QTOT(J).GT.QCONST)THEN
QTOT(J)=QCONST
QGH(J)=QCONST
ENDIF
ENDIF
C
C SWITCH TO OVERWRITE CO2 CONCS CALCULATED BY CARBON CYCLE MODEL.
C ICO2READ=1, THEN DELTA-QTOT FROM 1990 IS A DIRECT MULTIPLE
C OF THE CO2 FORCING WITH SCALING FACTOR FROM CFG FILE.
C ICO2READ=2, QTOT IS THE SUM OF THE NEW CO2 FORCING AND OTHER
C FORCINGS AS DETERMINED BY THE GASEM.$$$ EMISSIONS INPUTS.
C ICO2READ=3, QTOT IS THE SUM OF THE NEW CO2 FORCING AND SULPHATE
C FORCING AS DETERMINED BY THE EMISSIONS INPUTS.
C ICO2READ=4, QTOT IS THE SUM OF THE NEW CO2 FORCING AND SULPHATE
C FORCING, WITH THE LATTER DETERMINED BY ESO2 FROM SCALED
C FOSSIL EMISSIONS.
C
IF(ICO2READ.EQ.1.AND.J.GT.226)THEN
QTOT(J)=QTOT(226)+CO2SCALE*(QXX*ALOG(CO2(J)/CO2(226)))
QSO2(J)=QSO2(226)
QDIR(J)=QDIR(226)
QOZ(J) =QOZ(226)
QBIO(J)=QBIO(226)
QGH(J) =QGH(226)+CO2SCALE*(QXX*ALOG(CO2(J)/CO2(226)))
ENDIF
C
C IF ICO2READ=2, THE CHANGE REQUIRED IS ONLY FOR CO2 AND THIS HAS
C ALREADY BEEN MADE. IF ICO2READ=4, THE SO2 SCALING HAS ALREADY
C BEEN DONE JUST AFTER READING GASEM.$$$
C
IF(ICO2READ.GE.3.AND.J.GT.226)THEN
QTOT(J)=QTOT(226)+CO2SCALE*(QXX*ALOG(CO2(J)/CO2(226)))
QTOT(J)=QTOT(J)+(QSO2(J)-QSO2(226))
QOZ(J) =QOZ(226)
QBIO(J)=QBIO(226)
QGH(J) =QGH(226)+CO2SCALE*(QXX*ALOG(CO2(J)/CO2(226)))
ENDIF
C
IF(J.EQ.225)THEN
QTOT89=QTOT(J)
QSO289=QSO2(J)
QDIR89=QDIR(J)
QOZ89 =QOZ(J)
QBIO89=QBIO(J)
QGH89 =QGH(J)
ENDIF
C
IF(J.EQ.226)THEN
QTOT90=QTOT(J)
QSO290=QSO2(J)
QDIR90=QDIR(J)
QOZ90 =QOZ(J)
QBIO90=QBIO(J)
QGH90 =QGH(J)
QTOTM =(QTOT89+QTOT90)/2.0
QSO2M =(QSO289+QSO290)/2.0
QDIRM =(QDIR89+QDIR90)/2.0
QOZM =(QOZ89+QOZ90)/2.0
QBIOM =(QBIO89+QBIO90)/2.0
QGHM =(QGH89+QGH90)/2.0
ENDIF
C
C SWITCH FOR SIMULATING COMPOUND CO2 GROWTH RATE BEGINNING IN MID
C 1990. GROWTH RATE SPECIFIED BY COMPRATE.
C
IF(ICOMP.EQ.1.AND.J.GE.226)THEN
QTOT(J)=QTOTM+QXX*(J-225.5)*(COMPRATE/100.)
QSO2(J)=QSO2M
QDIR(J)=QDIRM
QOZ(J) =QOZM
QBIO(J)=QBIOM
QGH(J) =QGHM+QXX*(J-225.5)*(COMPRATE/100.)
ENDIF
C
10 CONTINUE
C
RETURN
END
C
C *******************************************
C
SUBROUTINE INITCAR(NN,D80,F80)
C
parameter (iTp =740)
C
INTEGER FERTTYPE,TOTEM,CONVTERP
C
COMMON/COBS/COBS(0:226)
C
COMMON/CARB/CCO2(4,224:iTp),EDGROSS(4,226:iTp),EF(226:iTp),
&REGROW(4,226:iTp),PL(4,226:iTp),HL(4,226:iTp),SOIL(4,226:iTp),
&TTT(226:iTp),ESUM(226:iTp),ETOT(4,226:iTp),EDNET90(4),
&FOC(4,226:iTp),co2(0:iTp)
C
COMMON/CAR/EL1,EL2,EL3,TINV0(5),TINV(4,5),A(3,5),AA(4,5),
&BCO2(4),BTGPP,BTRESP,BTHUM,GAMP,GPP0,RESP0,QA0,U0,C0,B340(4),
&PHI,RG,TAUP,TAUH,TAUS,THP,THS,THH0,THS0,THPL,G1,G2,G3,FACTOR,
&EL21,EL32,XX1,XX2,XX3,XX4,XX5,XX6,DEE1,DEE2,DEE3,DEE4,DEE5,DEE6,
&FL1,FL2,FL3,XL,GAMH,GAMS,QS0,BTSOIL,FERTTYPE,TOTEM,CONVTERP,
&R(4),CPART(4,5),DELMASS(4,226:iTp),ABFRAC(4,226:iTp)
C
C FIRST INITIALISE PARAMETERS THAT DEPEND ON MM
C
C BCO2 ***************************
C
C COMPLETELY REVISED ON 6/12/95
C REVISED AGAIN ON 8/8/95. CUBIC RETAINED AS BASIC APPROXIMATION
C FORMULA FOR BCO2, BUT CORRECTION ADDED FOR D80<0.4 OR >1.8.
C FORMULA RAPIDLY LOSES ACCURACY OUTSIDE RANGE SHOWN BELOW.
C FORMULA ADDED TO ACCOUNT FOR CASES WHERE F80S.NE.2.0 USES
C FACT THAT BCO2 FOR (F80),(D80) IS APPROX THE SAME AS FOR
C (F80-DEL),(D80-DEL)
C
DD =(D80-1.1)-(F80-2.0)
DD2=DD*DD
DD3=DD*DD2
BCO2(NN)=0.27628+DD*0.301304+DD2*0.060673+DD3*0.012383
XD=ABS(DD)-0.7
IF(XD.GT.0.0) BCO2(NN)=BCO2(NN)-(0.0029*XD-.022*XD*XD)
C
C CORRECTION FOR F80.NE.2.0. ERROR IN BCO2 LESS THAN 0.0014
C FOR 1.0<F80<3.0
C
BCO2(NN)=BCO2(NN)-0.0379*(F80-2.0)
C
C COMPARISON OF FIT TO TRUE BCO2 FOR F80=2.0
C
C D80 BCO2(TRUE) BCO2(FIT) TRUE-FIT
C 0.0 0.00365 0.00414 -0.00049
C 0.1 0.02438 0.02438 0.00000
C 0.4 0.09085 0.09085 0.00000
C 0.8 0.19092 0.19102 -0.00010
C 1.1 0.27628 0.27628 0.00000
C 1.4 0.37237 0.37247 -0.00010
C 1.8 0.52117 0.52117 0.00000
C 2.2 0.69994 0.69997 -0.00003
C 2.6 0.91868 0.91830 0.00038
C 3.0 1.19350 1.18092 0.01258
C
C EDNET ETC. **********************
C
C LATEST (8/8/95) VALUES. NOT YET INSERTED.
C
C YEAR EDNET REGROW PLANT HLITT SOIL ETOT
C 1989 .724 1.251 707.895 84.195 1416.528 284.230
C 1990 .762 1.254 708.285 84.300 1416.655 289.707
C
C EDNET90 CORRECTED ON 8/8/95. LEADING TERMS ONLY IN OTHER ITEMS
C CORRECTED.
C
EDNET90(NN) =0.762+1.1185*DD-0.0020*DD2
C
REGROW(NN,226) =1.254+0.584*DD-0.0173*DD2
EDGROSS(NN,226)=EDNET90(NN)+REGROW(NN,226)
C
PL(NN,226) = 708.285-6.131*DD+0.137*DD2
HL(NN,226) = 84.300+4.569*DD-0.039*DD2
SOIL(NN,226) =1416.655+1.564*DD-0.099*DD2
C
C NEXT THREE ITEMS NOT UPDATED
C
FOC(NN,226) = 2.24
DELMASS(NN,226)= 3.609
ABFRAC(NN,226) = 0.514-0.084*(D80-1.0)+0.013*(D80-1.0)**2
C
C ******************************
C
C SPECIFY INIT, END-1988, END-1989 AND END-1990 CO2 CONCS.
C
C0 =COBS(0)
CCO2(NN,224)=COBS(224)
CCO2(NN,225)=COBS(225)
CCO2(NN,226)=COBS(226)
C
ETOT(NN,226) = 290.797
C
PL0 = 750.0
HL0 = 80.0
SOIL0=1450.0
C
FERTTYPE=2
TOTEM =1
CONVTERP=1
FACTOR =2.123
C
C ADJUST INVERSE DECAY TIMES ACCORDING TO 1980S-MEAN OCEAN FLUX.
C PSI CONTROLS THE MAGNITUDE OF THE FLUX INTO THE OCEAN.
C (IF PSI IS LESS THAN 1, FLUX IS LESS THAN IN THE ORIGINAL
C MAIER-REIMER AND HASSELMANN MODEL.)
C A SIMPLE BUT ACCURATE APPROXIMATE EMPIRICAL EXPRESSION IS USED
C TO ESTIMATE PSI AS A FUNCTION OF F80SIN. THE PSI-F80SIN RELATION
C DEPENDS ON THE ASSUMED HISTORY OF OBSERVED CONCENTRATION
C CHANGES. THUS, DIFFERENT PSI-F80SIN RELATIONSHIPS APPLY TO
C DIFFERENT CO2 HISTORIES.
C
C REVISED ON 3/16/95, BUT APSI ONLY (HENCE OK FOR F80=2.0 ONLY)
C REVISED AGAIN ON 4/2/95
C REVISED AGAIN ON 6/12/95
C REVISED AGAIN ON 8/8/95
C
FF =F80-2.0
FF2 =FF*FF
APSI=1.029606
BPSI=0.873692
CPSI=0.165084
PSI =APSI+BPSI*FF+CPSI*FF2
C
C PSI COMPARISON
C
C F80 PSITRUE PSIEST
C 1.0 .320730 .320998
C 1.5 .634031 .634031
C 2.0 1.029606 1.029606
C 2.5 1.507723 1.502733
C 3.0 2.074348 2.068382
C
DO J=1,5
TINV(NN,J)=TINV0(J)*PSI
END DO
C
C SPECIFY 1990 (J=226) PARTIAL CONCS AND CONVOLUTION CONSTANTS
C REVISED ON 3/16/95
C REVISED AGAIN ON 4/2/95
C REVISED AGAIN ON 6/12/95 (AA NOT CHANGED, CPART VERY MINOR)
C REVISED AGAIN ON 8/8/95 (AA AND CPART CHANGED)
C
CPART(NN,1)=18.40290
CPART(NN,2)=25.51023
CPART(NN,3)=22.51496
CPART(NN,4)= 9.64078
CPART(NN,5)= 0.38709
C
AA(NN,1)= 0.13486
AA(NN,2)= 0.22091
AA(NN,3)= 0.28695
AA(NN,4)= 0.26033
AA(NN,5)= 0.09695
C
C SPECIFY OR CALCULATE OTHER MAIN MODEL PARAMETERS
C
C GPP IS THE PART OF GPP THAT IS NOT IMMEDIATELY RESPIRED BY
C LEAVES AND GROUND VEG
C RESP IS THE PART OF RESPIRATION THAT COMES FROM THE TREES
C RG=RESP0/GPP0
C
GPP0 =76.0
RESP0=14.0
PHI =0.98
XL =0.05
G1 =0.35
G2 =0.60
GAMP =0.70
GAMH =0.05
C
RG =RESP0/GPP0
G3 =1.0-G1-G2
GAMS=1.0-GAMP-GAMH
C
C SPECIFY TEMPERATURE FEEDBACK TERMS.
C
BTGPP =0.00
BTRESP=0.00
BTHUM =0.00
BTSOIL=0.00
C
THPL=G1*GPP0-RESP0
TAUP=PL0/THPL
THP =1./(2.*TAUP)
TAUH=HL0/(G2*GPP0+PHI*THPL)
THH0=0.5/TAUH
TAUS=SOIL0/(GPP0-RESP0-(1.0-XL)*HL0/TAUH)
THS0=0.5/TAUS
C
C QA0 IS THE INITIAL HLITT DECOMPOSITION FLUX TO ATMOSPHERE.
C
QA0=(1.0-XL)*HL0/TAUH
C
C QS0 IS THE INITIAL HLITT FLUX TO SOIL : U0 DITTO SOIL TO ATMOS
C
QS0=XL*HL0/TAUH
U0 =SOIL0/TAUS
C
C FOC(NN,0) =0.0
C EF(0) =0.0
C ESUM(0) =0.0
C ETOT(NN,0)=0.0
C
FL1 =EL1/100.
FL2 =EL2/100.
FL3 =EL3/100.
EL21=FL2-FL1
EL32=FL3-FL2
XX1 =FL1-DEE1
XX2 =FL1+DEE2
XX3 =FL2-DEE3
XX4 =FL2+DEE4
XX5 =FL3-DEE5
XX6 =FL3+DEE6
C
RETURN
END
C
C *******************************************
c
SUBROUTINE CARBON(MM,TEM,EFOSS,ENETDEF,CPP,CPREV,C,
&PL ,HU ,SO ,REGRO ,ETOT ,
&PL1,HU1,SO1,REGRO1,ETOT1,
&SUMEM1,FLUX,DELM,EGROSSD,C1)
C
parameter (iTp =740)
C
DIMENSION A1(4,5)
C
INTEGER FERTTYPE,TOTEM,CONVTERP
C
COMMON/CAR/EL1,EL2,EL3,TINV0(5),TINV(4,5),A(3,5),AA(4,5),
&BCO2(4),BTGPP,BTRESP,BTHUM,GAMP,GPP0,RESP0,QA0,U0,C0,B340(4),
&PHI,RG,TAUP,TAUH,TAUS,THP,THS,THH0,THS0,THPL,G1,G2,G3,FACTOR,
&EL21,EL32,XX1,XX2,XX3,XX4,XX5,XX6,DEE1,DEE2,DEE3,DEE4,DEE5,DEE6,
&FL1,FL2,FL3,XL,GAMH,GAMS,QS0,BTSOIL,FERTTYPE,TOTEM,CONVTERP,
&R(4),CPART(4,5),DELMASS(4,226:iTp),ABFRAC(4,226:iTp)
C
C TERRESTRIAL CARBON CYCLE SECTION
C
C FOR CALCULATING FEEDBACKS, NEED TO USE EXTRAPOLATED CO2
C CONCENTRATION FOR THE MIDPOINT OF YEAR 'J' AFTER 1990.
C FORMULA BELOW GIVES QUADRATIC EXTRAPOLATED RESULT.
C
C CBAR=C+0.5*(C-CPREV)
CBAR=(3.0*CPP-10.0*CPREV+15.0*C)/8.0
C
C NOW CALCULATE GROSS DEFOR CORRESP TO INPUT NET DEFOR
C
REGRO1=(REGRO*(TAUP-0.5)+GAMP*ENETDEF)/(TAUP+0.5-GAMP)
EGROSSD=ENETDEF+REGRO1
C
C TEMPERATURE TERMS
C
FG=EXP(BTGPP*TEM)
FR=EXP(BTRESP*TEM)
FH=EXP(BTHUM*TEM)
FS=EXP(BTSOIL*TEM)
C
C DEFORESTATION TERMS
C
GDP=GAMP*EGROSSD
GDH=GAMH*EGROSSD
GDS=GAMS*EGROSSD
C
C CO2 FERTILIZATION TERMS. LOG FORM FIRST THEN FORM USED BY ENTING,
C THEN FORM USED BY GIFFORD. ALSO CALCULATE BCO2 AT C=CCC=340PPMV
C
IF(FERTTYPE.EQ.1)THEN
Y1=BCO2(MM)*ALOG(CBAR/C0)+1.0
B340(MM)=BCO2(MM)
ELSE
C
C CONCBASE=80.
C GINF=2.4
C BEE=C0*(GINF-1.0)-GINF*CONCBASE
C GEE=GINF*(CBAR-CONCBASE)/(CBAR+BEE)
C Y1=BCO2(MM)*(GEE-1.0)+1.0
C
CB=31.
R(MM)=(1.+BCO2(MM)*ALOG(680./C0))/(1.+BCO2(MM)*ALOG(340./C0))
AR=680.-CB
BR=340.-CB
BEE=999.9
IF(R(MM).NE.1.)BEE=(AR/BR-R(MM))/(R(MM)-1.)/AR
DR=CBAR-CB
CR=C0-CB
Y1=1.0
CCC=340.
B340(MM)=0.0
IF(R(MM).NE.1.)THEN
Y1=(1./CR+BEE)/(1./DR+BEE)
B340(MM)=(1./CR+BEE)*CCC/(1.+BEE*(CCC-CB))**2
ENDIF
ENDIF
C
GPP=GPP0*Y1*FG
C
PGPP=GPP*G1
DPGPP=PGPP-GPP0*G1
HGPP=GPP*G2
DHGPP=HGPP-GPP0*G2
SGPP=GPP*G3
DSGPP=SGPP-GPP0*G3
RESP=RESP0*Y1*FR
C
DRESP=RESP-RESP0
C
C NEW PLANT MASS
C
PTERM=PGPP-RESP-GDP
PL1=(PL*(1.0-THP)+PTERM)/(1.0+THP)
C
C NEW HLITT MASS
C
THH=THH0*FH
Y2=THP*(PL+PL1)
HTERM=HGPP+PHI*Y2-GDH
HU1=(HU*(1.0-THH)+HTERM)/(1.0+THH)
C
C NEW SOIL MASS
C
THS=THS0*FS
Y3=THH*(HU+HU1)
STERM=SGPP+(1.0-PHI)*Y2+XL*Y3-GDS
SO1=(SO*(1.0-THS)+STERM)/(1.0+THS)
C
C FERTILIZATION FEEDBACK FLUX
C
BFEED1=GPP0-RESP0-(GPP-RESP)
C
C FEEDBACK FLUX DUE TO ACTIVE TERRESTRIAL BIOMASS
C
HUBAR=(HU+HU1)/2.0
SOBAR=(SO+SO1)/2.0
QA1=FH*(1.0-XL)*HUBAR/TAUH
U1=FS*SOBAR/TAUS
AFEED1=QA1+U1-(QA0+U0)
C
C TOTAL BIOMASS CHANGE
C
DELSOIL=HU1+SO1-HU-SO
DELB=PL1-PL+DELSOIL
C
C ANN SUM OF EMISSIONS, OCEAN FLUX AND CUM SUM OF EMISSIONS.
C
SUMEM1=EFOSS-DELB
FLUX=SUMEM1-FACTOR*DELC
IF(TOTEM.EQ.1)ETOT1=ETOT+SUMEM1
IF(TOTEM.NE.1)ETOT1=ETOT+EFOSS+EGROSSD
C IF(TOTEM.NE.1)ETOT1=ETOT+EFOSS-DELB
C
C CALCULATE CONVOLUTION CONSTANTS AT THE END OF YEAR I. THIS
C MAY BE DONE IN TWO WAYS, DETERMINED BY CONVTERP. NORMALLY
C CONVTERP=1 SHOULD BE USED.
C
EE=ETOT1/100.
C
IF(CONVTERP.EQ.1)THEN
C
DO K=1,5
D21=(A(2,K)-A(1,K))/EL21
D32=(A(3,K)-A(2,K))/EL32
IF(EE.LE.XX1)THEN
A1(MM,K)=A(1,K)
C
ELSE IF(EE.LE.XX2)THEN
DY=D21*DEE2
X1=DEE1+DEE2
X12=X1*X1
AX=(D21-2.0*DY/X1)/X12
BX=(3.*DY/X1-D21)/X1
U=EE-XX1
U2=U*U
A1(MM,K)=A(1,K)+BX*U2+AX*U2*U
C
ELSE IF(EE.LE.XX3)THEN
U=EE-FL1
A1(MM,K)=A(1,K)+D21*U
C
ELSE IF(EE.LE.XX4)THEN
Y0=A(2,K)-D21*DEE3
DY=D21*DEE3+D32*DEE4
DD=D21+D32
X1=DEE3+DEE4
X12=X1*X1
AX=(DD-2.0*DY/X1)/X12
BX=(3.*DY/X1-DD-D21)/X1
U=EE-XX3
U2=U*U
A1(MM,K)=Y0+D21*U+BX*U2+AX*U2*U
C
ELSE IF(EE.LE.XX5)THEN
U=EE-FL2
A1(MM,K)=A(2,K)+D32*U
C
ELSE IF(EE.LE.XX6)THEN
Y0=A(3,K)-D32*DEE5
DY=D32*DEE5
X1=DEE5+DEE6
X12=X1*X1
AX=(D32-2.0*DY/X1)/X12
BX=(3.*DY/X1-2.0*D32)/X1
U=EE-XX5
U2=U*U
A1(MM,K)=Y0+D32*U+BX*U2+AX*U2*U
C
ELSE
A1(MM,K)=A(3,K)
ENDIF
C
END DO
C
ELSE
C
XXX1=0.75*FL1
XXX2=0.75*FL2
XXX3=1.25*FL2
XXX4=1.25*FL3
DX12=XXX2-XXX1
DX23=XXX4-XXX3
C
DO K=1,5
IF(EE.LE.XXX1)THEN
A1(MM,K)=A(1,K)
ELSE IF(EE.LE.XXX2)THEN
Z1=EE-XXX1
Z2=EE-XXX2
DY12=A(2,K)-A(1,K)
A1(MM,K)=A(1,K)+(DY12/DX12**3)*Z1*Z1*(Z1-3.0*Z2)
ELSE IF(EE.LE.XXX3)THEN
A1(MM,K)=A(2,K)
ELSE IF(EE.LE.XXX4)THEN
Z1=EE-XXX3
Z2=EE-XXX4
DY23=A(3,K)-A(2,K)
A1(MM,K)=A(2,K)+(DY23/DX23**3)*Z1*Z1*(Z1-3.0*Z2)
ELSE
A1(MM,K)=A(3,K)
ENDIF
END DO
C
ENDIF
C
C ******************************************************************
C
C CALCULATE NEW PARTIAL CONCS AND NEW CONCENTRATION.
C
DELC=0.0
DO J=1,5
DELA=A1(MM,J)-AA(MM,J)
ABAR=(A1(MM,J)+AA(MM,J))/2.0
Z=(TINV(MM,J)-DELA/ABAR)/2.0
DEL=(ABAR*SUMEM1/FACTOR-2.0*Z*CPART(MM,J))/(1.0+Z)
CPART(MM,J)=CPART(MM,J)+DEL
DELC=DELC+DEL
FLUX=SUMEM1-FACTOR*DELC
AA(MM,J)=A1(MM,J)
END DO
C1=C+DELC
CBAR=(C1+C)/2.0
C
C CALCULATE ANNUAL CHANGES IN ATMOSPHERIC MASS
C
DELM=FACTOR*(C1-C)
C
RETURN
END
C
C *******************************************
C
SUBROUTINE HISTORY(JJJ,CO2,CH4,CN2O,C11,C12,C22,C134EQ,eso2,eso21)
C
C THIS SUBROUTINE CALCULATES THE CONCS UP TO AND INCLUDING
C 1990 USING FITS TO OBSERVED DATA.
c Concs are end of year values. Conc(1) is end of 1765 etc.
c SO2 emissions values are for whole year, assumed to apply to midpoint
c of year.
C
COMMON/COBS/COBS(0:226)
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,ENAT,es1990,
&dqdcoz,QCH4OZ
COMMON /CH4HIST/METHHIST
C
C Y CORRESPS TO MID YEAR. THUS, IF JJJ=186 SAY (= END OF 1950)
C THEN Y=1950.5 (I.E., Y=1950.0 IS THE MIDDLE OF 1950)
C
Y=JJJ+1764.5
Y1=Y-1765.
C
C CO2 : END OF YEAR CONCS ARE SPECIFIED IN A DATA STATEMENT IN
C BLOCK DATA, IN ARRAY COBS. THIS ARRAY IS THE IPCC DATA SET
C GENERATED BY ENTING AND WIGLEY FOR THE IPCC CONC STABILIZATION
C EXERCISE. IMPLEMENTED ON DEC 30, 1993, UPDATED MAR 11, 1995
C
CO2=COBS(JJJ)
c
c new CH4 (9209)
C
if(y.lt.1940.) then
CH4=790.0+0.0101*Y1*Y1
else if(y.lt.1955.) then
CH4=1067.32+0.09765*(Y-1921.9)**2
else if(y.lt.1980.) then
Y6 = Y-2000.45
CH4=1804.68-0.7731*Y6**2-0.010296*Y6**3
else
C
C TWO VERSIONS OF 1980-1990 CONCS ....
C METHHIST
C (0) OLD VERSION, C(MID1990)=1717.0, DC/DT(MID1990)=10.828
C (1) NEW (1/10/98), C(MID1990)=1700.0, DC/DT(MID1990)=10.000
C
IF(METHHIST.EQ.0)THEN
CH4=1791.44-0.39375*(Y-2003.75)**2
ELSE
Y9=Y-1990.
YY9=Y9*Y9
YYY9=YY9*Y9
CH4=1700.+10.*Y9-0.0474*YY9+0.02584*YYY9
ENDIF
C
endif
C
C N2O REVISED ON JULY 15, 1994 TO ACCORD WITH IPCC94. CONCS
C HERE ARE MID YEAR VALUES. FIT IS LINEAR TO 1950, QUADRATIC TO
C 1965, LINEAR TO 1985 AND LINEAR WITH A DIFFERENT SLOPE TO
C 310.4 IN MID-1993. THE IPCC94 MID-1993 VALUE IS 310.0, A
C MINOR COMPROMISE IN ORDER TO GET A SIMPLE ANALYTIC FIT TO
C EARLIER PRESCRIBED VALUES.
C
YR0=1765.0
YR1=1950.0
YR2=1965.0
YR3=1985.0
CN0=270.0
C CN0=285.0
CN1=290.0
CN2=294.0
CN3=304.0
C
IF(Y.LT.YR1)THEN
CN2O=CN0+(CN1-CN0)*(Y-YR0)/(YR1-YR0)
ELSE IF(Y.LT.YR2)THEN
CN2O=288.667+(4.0/675.0)*(Y-1935.0)**2
ELSE IF(Y.LT.YR3)THEN
CN2O=CN2+0.5*(Y-YR2)
ELSE
CN2O=CN3+0.8*(Y-YR3)
ENDIF
C
C HALOCARBONS (revised 9209)
C
IF(Y.LT.1950.)THEN
C11=0.0
C134EQ=0.0
ELSE
CORR = 1.+(1990.-Y)/100.
C11=0.175*(Y-1950.)**2*CORR
ENDIF
C12=1.73*C11
factor = 0.33+0.08*(Y-1979.)/11.
C22=factor*C11
C
C EQUIV HFC134a CONC ADDED 921219, BASED ON LINEAR F14+F116 DELQ
C DELQ IS OFFSET BY 0.5 YEARS TO GET END OF YEAR CONC
C
IF(Y.GE.1950.)THEN
Q134EQ=0.00948*(Y-1950.)/40.
C134EQ = Q134EQ/0.000169
ENDIF
C
c SO2 EMISSIONS : NEED TO CALC JJJ & JJJ+1 VALUES.
c eso2 IS THE VALUE FOR YEAR J, NOMINALLY A MID-YEAR VALUE.
C eso21 IS THE VALUE FOR YEAR J+1, NEEDED TO CALCULATE AN
C EFFECTIVE END OF YEAR VALUE IN DETERMINING FORCING.
C es1990 IS THE FIRST INPUT VALUE FROM GASEM.$$$.
C
DO I=0,1
J=jjj+I
ymid = J+1764.
if(ymid.lt.1860.) then
ee = 0.0
else if(ymid.lt.1953) then
ee = 35.0*(ymid-1860.)/93.0
else if(ymid.lt.1973) then
ee = 35.0+33.0*(ymid-1953.)/20.
else
ee = 68.0+(es1990-68.0)*(ymid-1973.)/17.
endif
IF(I.EQ.0)eso2=ee
IF(I.EQ.1)eso21=ee
END DO
C
RETURN
END
c
c--------------------------
c
c New methane model, Aug92
c
subroutine methane(ilmh,cm1,c,em,eno,eco,evo,
+ taum2,taum1,DELT,c1,tau)
c
common /TauMeth/ITAUMETH,TAU90CH4,TauSoil,DELTAU,
+ TAUINIT,ch4bar90,eno90,eco90,evo90
common /powc/powc0,delpowc,powco0,delpowco
common /pown/powno0,delpowno,powvo0,delpowvo
c
c Temperature feedback term, parameters specified below
c
Q10=1.0
BIOTA0=100.0
DELT90=0.6
C
BIOTA = BIOTA0*Q10**(DELT/10.)
BIOTA90 = BIOTA0*Q10**(DELT90/10.)
DELBIOTA = BIOTA-BIOTA90
emm = em+DELBIOTA
c
if(ilmh.eq.1) then
POWC = powc0-delpowc
POWCO = powco0-delpowco
POWNO = powno0-delpowno
POWVO = powvo0-delpowvo
tau90 = TAUINIT-DELTAU
else if(ilmh.eq.2) then
POWC = powc0
POWCO = powco0
POWNO = powno0
POWVO = powvo0
tau90 = TAUINIT
else if(ilmh.eq.3) then
POWC = powc0+delpowc
POWCO = powco0+delpowco
POWNO = powno0+delpowno
POWVO = powvo0+delpowvo
tau90 = TAUINIT+DELTAU
endif
c
b = 2.75
cosource = 270.
vosource = 184.

if(ilmh.eq.4) then
tau = TAU90CH4
yy = 0.
else
xcos90 = vosource+cosource
c
c guessc = first order estimate of mid-year value in year J (=0.5*(C1+C))
c guesst = first order estimate of tau = mean lifetime over year during
c which conc changes from C to C1
c
guessc = 1.5*c - 0.5*cm1
guesst = 1.5*taum1 - 0.5*taum2
xcos = cosource*(guessc/ch4bar90)*(tau90/guesst)
+ + vosource*(evo/evo90)
c
tau = tau90*((c/ch4bar90)**powc)
+ *((eno/eno90)**powno)
+ *(((eco+xcos)/(eco90+xcos90))**powco)
+ *((evo/evo90)**powvo)
yy = powc/(2.*tau)
endif

y = 1./tau + 1./TauSoil
aa = emm/b - c*y
bb = 1. + y/2. - yy
c1 = c + aa/bb
c
return
end
C
C ********************************
C
SUBROUTINE NITROUS(C,EM,C1)
common /TauNitr/TAUN2O
BB=4.81
TT=1.0/(2.0*TAUN2O)
C1=(C*(1.0-TT)+EM/BB)/(1.0+TT)
RETURN
END
C
C ***************************************
C
SUBROUTINE ICEMELT(DTEM,S0BIT,G0,A0,S1BIT,S1,G1,A1)
c
c S0 = Small glaciers input S1 = 1 year increment
c G0 = Greenland etc
c A0 = Antartica
c
C NOTE : S0 NOT USED DIRECTLY - INPUT IS S0BIT
C
COMMON /ICE/TAUI,DTSTAR,DTEND,BETAG,BETA1,BETA2,DB3
COMMON /ICEPARA/ZI0,NVALUES,MANYTAU,TAUFACT,ZI0BIT(100)
DIMENSION DTST(100),S0BIT(100),S1BIT(100)
C
TAU=TAUI
DTSTBIT=(DTEND-DTSTAR)/(NVALUES-1)
DTSTMIN=DTSTAR
DTSTMAX=DTEND
TAUMIN=(1.0-TAUFACT)*TAU
TAUMAX=(1.0+TAUFACT)*TAU
C
S1=0.0
C
DO 1 N=1,NVALUES
C
DTST(N)=DTSTAR+(N-1)*DTSTBIT
IF(MANYTAU.EQ.1)THEN
DELTAU=(TAUMAX-TAUMIN)/(NVALUES-1)
TAU=TAUMIN+(N-1)*DELTAU
ENDIF
C
S1BIT(N)=(S0BIT(N)*(TAU-0.5)+ZI0BIT(N)*DTEM/DTST(N))
&/(TAU+0.5)
IF(S1BIT(N).GE.ZI0BIT(N))THEN
S1BIT(N)=ZI0BIT(N)
ENDIF
C
S1=S1+S1BIT(N)
C
1 CONTINUE
C
G1=G0+BETAG*DTEM*1.5
A1=A0+DB3+(BETA1+BETA2)*DTEM
C
RETURN
END
C
C ****************************************
C
SUBROUTINE INTERP(N,ISTART,IY,X,Y)
c
parameter (iTp =740)
c
common /Limits/KEND
c
DIMENSION IY(100),X(100),Y(226:iTp)
IEND=ISTART+IY(N)
DO I=0,IY(N)-1
DO K=1,N
IF(I.GE.IY(K).AND.I.LT.IY(K+1))THEN
J=I+ISTART
Y(J)=X(K)+(I-IY(K))*(X(K+1)-X(K))/(IY(K+1)-IY(K))
ENDIF
END DO
END DO
Y(IEND)=X(N)
c
c If last year in profile (relative to 1990) not KEND, then assume
c constant emissions from last year specified in emissions profile
c to KEND.
c
if(iy(n).lt.KEND) then
do i=iend+1,KEND
y(i) = x(n)
end do
end if
c
RETURN
END
C
C*************************************************
C
SUBROUTINE SULPHATE(JY,ESO2,ESO21,QSO2,QDIR,QOZ)
c
parameter (iTp =740)
C
COMMON /IGHG/IGHG
common /Sulph/S90DIR,S90IND,S90BIO,S90OZ,ENAT,es1990,
&dqdcoz,QCH4OZ
c
c Tall stack effect factor
c
if(jy.lt.186) then
f = 0.7
else if(jy.lt.206) then
f = 0.7 + 0.3*(jy-186)/20.
else
f = 1.0
endif
c
ky = jy + 1
if(ky.lt.186) then
f1 = 0.7
else if(ky.lt.206) then
f1 = 0.7 + 0.3*(ky-186)/20.
else
f1 = 1.0
endif
c
c Calculate end of year emissions and ditto corrected for tall
c stack effect
c
eraw=(eso2+eso21)/2.0
e = (f*eso2+f1*eso21)/2.
c
c Calculate global forcing at end of year.
c
c Initialise SO2 parameters. ORIGINALLY USED VALUES BASED ON
C CHARLSON ET AL. ALTERED TO FIT IPCC94 ON OCT 16, 1994.
c FURTHER MODIFIED ON FEB 24, 1995 FOR IPCC95
c MOVED FROM SUBROUTINE INIT TO SUBROUTINE SULPHATE ON 3/10/95
C WITH ASO2 AND BSO2 ELIMINATED
c
qdir = e*s90dir/es1990
qindir = s90ind*(alog(1.0+e/ENAT))/(alog(1.0+es1990/ENAT))
C
c Calculate global tropospheric ozone forcing. (Separate from that
c due to ch4 changes.)
C
C REVISED OZONE HISTORY (parallels fossil fuel emissions)
C NOTE THAT THE KEY COZ VALUES ARE SPECIFIED AS MID YEAR VALUES, BUT
C THIS SUBROUTINE REQUIRES END OF YEAR OUTPUT FOR QOZ. HENCE, 0.5
C YEAR OFFSET NEEDED.
C
IY=JY+1764
IF(IY.LT.1850)THEN
COZ=25.0
ELSE IF(IY.LT.1890)THEN
COZ=25.0+21.0*((IY-1850+0.5)/40.0)*(0.3/6.1)
ELSE IF(IY.LT.1945)THEN
COZ=25.0+21.0*(0.3/6.1)+((IY-1890+0.5)/55.0)*(1.0/6.1)*21.0
ELSE IF(IY.LT.1990)THEN
COZ=25.0+21.0*(1.3/6.1)+((IY-1945+0.5)/45.0)*(4.8/6.1)*21.0
ELSE
COZ=46.0
ENDIF
QOZ=(S90OZ/21.0)*(COZ-25.0)
C
IF(IGHG.LE.1)QOZ=0.0
qso2 = qdir+qindir
C
return
end
C
C ******************************************************************
C
SUBROUTINE LAMCALC(Q,FNHL,FSHL,XK,XKH,DT2X,A,LAMOBEST,LAMLBEST)
C
C Revision history:
C 950215 : CONVERTED TO SUBROUTINE FOR STAG.FOR
C 950208 : ITERATION ALGORITHM IMPROVED YET AGAIN
C 950206 : ITERATION ALGORITHM IMPROVED AGAIN
C 950205 : MATRIX INVERSION CORRECTED
C 950204 : ITERATION ALGORITHM IMPROVED
C 950203 : FIRST VERSION OF PROGRAM WRITTEN
C
C THIS SUBROUTINE CALCULATES LAND AND OCEAN FEEDBACK
C PARAMETER VALUES GIVEN THE GLOBAL DT2X AND THE LAND TO
C OCEAN EQUILIBRIUM WARMING RATIO (A).
C BOTH THE INPUT VALUES OF FNHL AND FSHL, AND THE INPUT XK
C (XKLO IN MAIN) AND XKH (XKNS IN MAIN) ARE DOUBLE WHAT
C ARE USED HERE.
C
DIMENSION AEST(100),DIFF(100)
REAL LAMO(100),LAML(100),LAMOBEST,LAMLBEST,LAM,KLO,KNS
C
KLO=XK/2.0
KNS=XKH/2.0
FNL=FNHL/2.0
FSL=FSHL/2.0
C
IMAX=40
DLAMO=1.0
DIFFLIM=0.001
C
FNO=0.5-FNL
FSO=0.5-FSL
FL=FNL+FSL
FO=FNO+FSO
FRATIO=FO/FL
C
DT2XO=DT2X/(FO+A*FL)
DT2XL=A*DT2XO
LAM=Q/DT2X
LAMO(1)=LAM
LAMO(2)=LAM+DLAMO
C
IFLAG=0
DO 1 I=1,IMAX
C
LAML(I)=LAM+FRATIO*(LAM-LAMO(I))/A
C
C SOLVE FOR NH/SH OCEAN/LAND TEMPS
C FIRST SPECIFY COEFFICIENT MATRIX : A(I,J)
C
A11= FNO*LAMO(I)+KLO+KNS
A12=-KLO
A13=-KNS
A14= 0.0
A22= FNL*LAML(I)+KLO
A23= 0.0
A24= 0.0
A33= FSO*LAMO(I)+KLO+KNS
A34= A12
A44= FSL*LAML(I)+KLO
C
C CALCULATE INVERSE OF COEFFICIENT MATRIX : B(I,J)
C FIRST DETERMINE DETERMINANT OF A(I,J) MATRIX
C
C1 = A11*A22-A12*A12
C2 = A33*A44-A12*A12
C3 = A22*A13
C4 = A44*A13
DET= C1*C2-C3*C4
C
B11= A22*C2/DET
B12=-A12*C2/DET
B13=-A44*C3/DET
B14= A12*C3/DET
B22= (A11*C2-A13*C4)/DET
B23= A12*C4/DET
B24=-A12*A12*A13/DET
B33= A44*C1/DET
B34=-A12*C1/DET
B44= (A33*C1-A13*C3)/DET
C
C CALCULATE ESTIMATED NH/SH OCEAN/LAND EQUILIBRIUM TEMPS
C
TNO=(B11*FNO+B12*FNL+B13*FSO+B14*FSL)*Q
TNL=(B12*FNO+B22*FNL+B23*FSO+B24*FSL)*Q
TSO=(B13*FNO+B23*FNL+B33*FSO+B34*FSL)*Q
TSL=(B14*FNO+B24*FNL+B34*FSO+B44*FSL)*Q
C
C CALCULATE ESTIMATED OCEAN-MEAN AND LAND-MEAN TEMPS
C
DT2XLE=(TNL*FNL+TSL*FSL)/FL
DT2XOE=(TNO*FNO+TSO*FSO)/FO
C
C CALCULATE ESTIMATED N.H. AND S.H. TEMPERATURES
C
TNH=(TNL*FNL+TNO*FNO)/0.5
TSH=(TSL*FSL+TSO*FSO)/0.5
C
C CALCULATE ESTIMATED VALUE OF A
C
AEST(I)=DT2XLE/DT2XOE
AAAA=AEST(I)
DIFF(I)=A-AEST(I)
C
C TEST DIFF TO DECIDE WHETHER TO END ITERATION LOOP
C
IF(ABS(DIFF(I)).LT.DIFFLIM)GO TO 2
C
IF(I.GE.2)THEN
DD=DIFF(I)*DIFF(I-1)
C
IF(DD.LT.0.0)THEN
IFLAG=1
ELSE
IF(ABS(DIFF(I)).GT.ABS(DIFF(I-1)))DLAMO=-DLAMO
LAMO(I+1)=LAMO(I)+DLAMO
ENDIF
C
IF(IFLAG.EQ.1)THEN
IF(DD.LT.0.0)THEN
RATIO=(LAMO(I)-LAMO(I-1))/(DIFF(I)-DIFF(I-1))
LAMO(I+1)=LAMO(I)-RATIO*DIFF(I)
ELSE
RATIO=(LAMO(I)-LAMO(I-2))/(DIFF(I)-DIFF(I-2))
LAMO(I+1)=LAMO(I)-RATIO*DIFF(I)
ENDIF
ENDIF
C
ENDIF
C
1 CONTINUE
2 CONTINUE
C
LAMOBEST=LAMO(I)
LAMLBEST=LAML(I)
C
RETURN
END
1 * ISCENGEN : IF=1, RUNS 20 SIMULATIONS TO DRIVE SCENGEN : SEE BELOW
0 * IUSERGAS : IF=0, USES DEFAULT VALUES FOR GAS CYCLE PARAMETERS.
2.5 *OK DT2XUSER : USER VALUE OF DT2X IN CLIMATE MODEL LOOP
1.3 * RLO : EQUIL LAND/OCEAN TEMP RATIO (USED IF IXLAM=1)
90.0 * HM : MIXED LAYER DEPTH (M)
1.0 * YK : VERTICAL DIFFUSIVITY (CM**2/SEC)
4.0 * W0 : INITIAL UPWELLING RATE (M/YR) : CAN USE W=4.0*YK
0.2 * PI : RATIO, SINKING WATER DELT TO GLOBAL DELT
1 * IVARW : =1, VARIABLE UPWELLING: =2, CONST W=WTHRESH AFTER TW0
7.0 * TW0NH : TEMP AT WHICH LINEAR DECLINE IN WNH REACHES ZERO
7.0 * TW0SH : TEMP AT WHICH LINEAR DECLINE IN WSH REACHES ZERO
30.0 * ZI0IN : USER VALUE OF ZI0
100.0 * TAUIN : USER VALUE OF GSIC TAU
0.7 * DTSTARIN : USER VALUE OF DTSTAR (LOWER BOUND VALUE)
3.0 * DTENDIN : USER VALUE OF DTEND (UPPER BOUND FOR DTSTAR)
0.03 * BETAGIN : USER VALUE OF GREENLAND SENSITIVITY
-0.03 * BETA1IN : USER VALUE OF ANTARCTIC SENSITIVITY (1)
0.01 * BETA2IN : USER VALUE OF ANTARCTIC SENSITIVITY (2)
0.01 * DB3IN : USER VALUE OF ANTARCTIC INITIAL MASS BALANCE
1990 *OK KYRREF : REFERENCE YEAR FOR CLIMATE MODEL OUTPUT
1990 *OK IY0 : FIRST YEAR FOR CLIMATE MODEL OUTPUT
2100 *OK LASTYEAR : LAST YEAR FOR CLIMATE MODEL RUN
5 *OK? ITEMPRT : PRINTOUT INTERVAL FOR CLIMATE MODEL
2 * ITAUMETH : 4=CONST TAU.
9.08 * TAU90CH4 : USED ONLY IF ITAUMETH=4.
8.30 * TAUINIT : value of ch4 tau90 (in variable tau formula)
0.5 * DELTAU : TAUINIT RANGE FOR LO, MID, HI CONCS (FOR SAR USE 0.8).
1 * IMETH : IF=1 ADDS IN CH4 OXIDATION TO CO2 EMISSIONS
1 * IO3FEED : IF=1 INCLUDES STRAT OZONE DEPLETION FEEDBACK
-0.30 * S90DIR : 1990 VALUE OF DIRECT AEROSOL FORCING (INCL SOOT)
-0.800 * S90IND : 1990 VALUE OF INDIRECT AEROSOL FORCING
-0.2 * S90BIO : 1990 VALUE OF BIOMASS AEROSOL FORCING
1.1 * DUSER : USER INPUT VALUE FOR D80S IN CARBON CYCLE MODEL
2.0 * FUSER : USER INPUT VALUE FOR F80S IN CARBON CYCLE MODEL
120.0 * TAUN2O : N20 LIFETIME
1234567890 *
: ISCENGEN : =0 RUNS 4 CASES FOR MAGICC ALONE, =1 RUNS 20 CASES TO GET
SULPHATE PATTERN WEIGHTS, =9 RUNS USER CASE ONLY
: METHANE : IF IUSERGAS=0, DEFAULT TAUINIT (=9.08) IS USED AND RESULTS
AGREE WELL WITH IPCC SAR CH.2 (PRATHER). USE THIS WITH MAG-F2.FOR.
IF IUSERGAS=1 AND TAUINIT IS SET AT 8.30, THEN RESULTS
WILL AGREE REASONABLY WELL WITH IPCC SAR CH.6 (WIGLEY). USE THIS
WITH MAG-F AND 232 VERSION OF MAGICC.
1764 277.951
1765 277.980
1766 278.011
1767 278.045
1768 278.082
1769 278.121
1770 278.164
1771 278.210
1772 278.258
1773 278.310
1774 278.364
1775 278.422
1776 278.483
1777 278.546
1778 278.613
1779 278.683
1780 278.756
1781 278.833
1782 278.912
1783 278.995
1784 279.081
1785 279.171
1786 279.263
1787 279.359
1788 279.459
1789 279.561
1790 279.668
1791 279.777
1792 279.890
1793 280.007
1794 280.127
1795 280.250
1796 280.377
1797 280.507
1798 280.640
1799 280.776
1800 280.914
1801 281.054
1802 281.195
1803 281.337
1804 281.479
1805 281.622
1806 281.764
1807 281.906
1808 282.047
1809 282.185
1810 282.322
1811 282.457
1812 282.588
1813 282.717
1814 282.841
1815 282.962
1816 283.078
1817 283.190
1818 283.297
1819 283.401
1820 283.502
1821 283.600
1822 283.696
1823 283.790
1824 283.883
1825 283.974
1826 284.066
1827 284.157
1828 284.249
1829 284.343
1830 284.438
1831 284.535
1832 284.634
1833 284.738
1834 284.845
1835 284.956
1836 285.072
1837 285.193
1838 285.321
1839 285.454
1840 285.594
1841 285.737
1842 285.881
1843 286.026
1844 286.168
1845 286.309
1846 286.448
1847 286.583
1848 286.716
1849 286.845
1850 286.971
1851 287.094
1852 287.214
1853 287.330
1854 287.442
1855 287.552
1856 287.658
1857 287.763
1858 287.868
1859 287.973
1860 288.078
1861 288.185
1862 288.294
1863 288.407
1864 288.524
1865 288.646
1866 288.773
1867 288.905
1868 289.044
1869 289.188
1870 289.338
1871 289.495
1872 289.659
1873 289.830
1874 290.009
1875 290.196
1876 290.390
1877 290.592
1878 290.801
1879 291.017
1880 291.241
1881 291.470
1882 291.704
1883 291.944
1884 292.189
1885 292.437
1886 292.688
1887 292.943
1888 293.200
1889 293.457
1890 293.715
1891 293.971
1892 294.225
1893 294.475
1894 294.721
1895 294.963
1896 295.201
1897 295.439
1898 295.677
1899 295.917
1900 296.160
1901 296.408
1902 296.663
1903 296.926
1904 297.198
1905 297.477
1906 297.764
1907 298.056
1908 298.353
1909 298.655
1910 298.958
1911 299.265
1912 299.573
1913 299.884
1914 300.196
1915 300.509
1916 300.823
1917 301.137
1918 301.451
1919 301.766
1920 302.080
1921 302.394
1922 302.707
1923 303.018
1924 303.326
1925 303.630
1926 303.929
1927 304.222
1928 304.510
1929 304.792
1930 305.069
1931 305.343
1932 305.613
1933 305.880
1934 306.147
1935 306.413
1936 306.680
1937 306.950
1938 307.225
1939 307.506
1940 307.794
1941 308.092
1942 308.400
1943 308.720
1944 309.054
1945 309.403
1946 309.766
1947 310.146
1948 310.542
1949 310.956
1950 311.388
1951 311.841
1952 312.314
1953 312.809
1954 313.328
1955 313.870
1956 314.439
1957 315.035
1958 315.659
1959 316.310
1960 316.983
1961 317.673
1962 318.378
1963 319.106
1964 319.872
1965 320.692
1966 321.574
1967 322.517
1968 323.525
1969 324.594
1970 325.716
1971 326.883
1972 328.086
1973 329.311
1974 330.555
1975 331.834
1976 333.163
1977 334.544
1978 335.965
1979 337.415
1980 338.885
1981 340.371
1982 341.880
1983 343.420
1984 344.992
1985 346.601
1986 348.193
1987 349.768
1988 351.335
1989 352.865
1990 354.407
YEAR CONC
END OF YEAR CONCS AS USED IN IPCC95

No comments:

Post a Comment