Thursday, March 22, 2012

2726.txt

cc: Mike Hulme <m.hulmeatXYZxyz.ac.uk>
date: 27 Feb 1998 14:40:37 -0700
from: Tom Wigley <Tom_WigleyatXYZxyzte.ucar.edu>
subject: MAG.F and EXE (2 of 2)
to: Mike Salmon <m.salmonatXYZxyz.ac.uk>


REGARDING MAG.F and EXE (2 of 2)


C MAG.FOR

C

C Revision history:

C

C 980223 * NVALUES PUT BACK INTO MAGUSER.CFG.

C 980216 * FORMAT STATEMENT 226 RETURNED TO ORIGINAL VERSION. ALL FILE

C NAMES 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 =550)

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

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,4240) NVALUES

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,4240) IQREAD

READ(LUN,4241) QOFFSET

READ(LUN,4241) QFACTOR

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

TAU90CH4 = 9.08

TAUINIT = 9.08

DELTAU = 0.0

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

close(lun)

C

C SPECIFICATION OF CH4 MODEL PARAMS MOVED TO HERE FOR CONVENIENCE.

C

POWC0 = 0.030

DELPOWC= 0.015

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,4243) NVAL

read(lun,'(a)') mnem

read(lun,*) ! skip description

read(lun,*) ! skip column headings

read(lun,*) ! skip units

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 TAUINIT. 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

IF(ITAUMETH.EQ.1)T90USER=T90LO

IF(ITAUMETH.EQ.1)T90USER=T90MID

IF(ITAUMETH.EQ.1)T90USER=T90HI

IF(ITAUMETH.EQ.1)T90USER=T90CON

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

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

WRITE (8,222) IYEAR,EF(K),EDNET(K),ECH4(K),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)

4444 FORMAT(6X,F10.0)

C

END

C

C********************************************************************

C

BLOCK DATA

c

parameter (iTp =550)

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 =550)

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 /TauMeth/ITAUMETH,TAU90CH4,TauSoil,DELTAU,

+ TAUINIT,ch4bar90,eno90,eco90,evo90

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 =550)

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 =550)

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)=0.0

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)=0.0

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 =550)

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 =550)

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 =550)

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

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

powno = -0.083

powco = 0.229

powvo = -0.026

tau90 = TAUINIT-DELTAU

else if(ilmh.eq.2) then

POWC = powc0

powno = -0.166

powco = 0.457

powvo = -0.051

tau90 = TAUINIT

else if(ilmh.eq.3) then

POWC = powc0+delpowc

powno = -0.249

powco = 0.686

powvo = -0.077

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 =550)

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 =550)

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

No comments:

Post a Comment