cc: Mike Hulme date: 27 Feb 1998 14:40:37 -0700 from: Tom Wigley subject: MAG.F and EXE (2 of 2) to: Mike Salmon 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), ®ROW(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), ®ROW(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), ®ROW(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), ®ROW(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), ®ROW(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), ®ROW(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