#define NoWeatheringFeedback #define CarbonateLandAreaFactor #define Carbon13 c#define SustainedCO2Slug c#define LTMIP c#define LTMIP_Header c#define LTMIP_Hist program geocarb implicit none integer np, n_record parameter(np=2) integer i_alk,i_tc,i_tc_13,it,ip,nsteps,is,i_print,i_blew_already real*8 co2_weat(np),f_land(np),f_plants(np),co2_degas(np) real*8 co2_slug(2), MOLpPPM, MOLpMM3 #ifdef Carbon13 real*8 ocn_conc(3), ocn_inv(3) #else real*8 ocn_conc(2), ocn_inv(2) #endif #ifdef LTMIP character*3 instring #endif #ifdef LTMIP_Hist real year_ipcc(250), pco2_ipcc(250) integer i_ipcc #endif real*8 mean_latitude, geol_year, dt2x,GEOG,weat_carb_k,weat_sil_k, $ time_step,ocn_vol,GEpPPM,pco2(2),atm_temp_0,ocn_temp_0, $ atm_temp,ocn_temp,ocn_tco2_bot,ocn_tco2_top, $ pco2_test,co3,time,co2_doublings,ocn_temp_target,r_co2, $ Solar_Const,RUN,tau,Z,f_bt_co2,WEAT_CARB,weat_sil,pco2_sat, $ BC,f_caco3_T, $ atm_inv(2),GE, ge13, pco2_13_eq,co2_slug_d13c read(5,*) geol_year, ! 1 position $ co2_degas, ! 2,3 $ mean_latitude, ! 4 $ f_plants, ! 5,6 $ dt2x, ! 7 $ f_land, ! 8,9 $ co2_slug(1), ! 10 $ co2_slug_d13c c . , co2_slug(1,2) co2_slug(2) = co2_slug(1) * ( 1.d0 + co2_slug_d13c/1.d3 ) #ifdef LTMIP call getarg(1,instring) #endif c#ifdef LTMIP_Header write(6,'(a15,5a12,7a8,3a12)') "year", $ "tco2", "alk", $ "d13Cocn", $ "pCO2", "d13Catm","CO3", $ "WeatC", "WeatS", $ "TotW","BurC", $ "Degas", "Tatm","Tocn" c#endif c write(6,*)year,co2_degas,mean_latitude,f_plants,dt2x, c $ f_land, co2_slug GEOG = (cos(mean_latitude/360*2*3.14159)-0.5) $ * 10. - 3. ! T offset from mean abs latitude of continents ! 0 at 30 N or S co2_slug(1) = co2_slug(1) * 1000. / 12. ! gton c to e12 mol #ifdef Carbon13 co2_slug(2) = co2_slug(2) * 1000. / 12. ! gton c to e12 mol #endif c parameters weat_carb_k = .00261 * 3000. weat_sil_k = 7. c DT2x = 3. #ifdef LTMIP if(instring .EQ. "___" .OR. instring .EQ. "_s_") then dt2x = 0. endif #endif time_step = 50. ! years #ifdef LTMIP time_step = 10. ! years #endif #ifdef LTMIP_Hist time_step = 1. #endif MOLpPPM = 2. * 1.e3 / 12. ocn_vol = 3.e2 * 4000 ! e12 m^3 MOLpMM3 = ocn_vol GEpPPM = 1./300. * MOLpPPM ! mol/yr i_blew_already = 0 i_alk = 1 i_tc = 2 #ifdef Carbon13 i_tc_13 = 3 #endif c initial condition pco2(1) = 200. #ifdef Carbon13 pco2(2) = 200. #endif atm_temp_0 = 15. ocn_temp_0 = 4. ocn_temp = ocn_temp_0 ocn_conc(i_alk) = 2.3 ! mol/m^3 ocn_tco2_top = 5. ocn_tco2_bot = 1. do it = 1,10 ocn_conc(i_tc) = (ocn_tco2_top+ocn_tco2_bot) / 2 call calc_co2(ocn_conc(i_tc),ocn_conc(i_alk), $ ocn_temp_0,pco2_test,co3) if(pco2_test .GT. pco2(1)) then ocn_tco2_top = ocn_conc(i_tc) else ocn_tco2_bot = ocn_conc(i_tc) endif enddo #ifdef Carbon13 ocn_conc(i_tc_13) = ocn_conc(i_tc) #endif time = 0. do ip = 1,np if(ip .EQ. 1) then nsteps = 5.0e6 / time_step else nsteps = 2.0e6 / time_step endif WEAT_SIL = 0. do it = 1,nsteps time = (ip-2)*nsteps*time_step $ + (it-1)*time_step #ifdef LTMIP if(time .LT. 0) then pco2 = 278. if(instring .EQ. "csw") then co2_degas = WEAT_SIL endif endif #endif #ifdef SustainedCO2Slug if(time .GE. 0 .AND. time .LE. 30000.) then atm_inv = atm_inv + co2_slug / 30000. * time_step pco2 = atm_inv / MOLpPPM endif c#else if(time .GE. 0. .AND. i_blew_already .EQ. 0) then atm_inv(:) = atm_inv(:) + co2_slug(:) pco2 = atm_inv / MOLpPPM i_blew_already = 1 endif #endif #ifdef LTMIP_Hist if(ip .EQ. 2) then if(it .EQ. 1) then OPEN(UNIT=11,FILE='pco2.ipcc.a.dat') do i_ipcc = 1, 250 read(11,*) year_ipcc(i_ipcc), $ pco2_ipcc(i_ipcc) enddo elseif(time .GE. 1750 .AND. time .LE. 2000) then do i_ipcc = 1, 250 if(year_ipcc(i_ipcc) .LE. time) then pco2 = pco2_ipcc(i_ipcc) endif enddo endif endif #endif co2_doublings = LOG(pco2(1)/280.)/LOG(2.) atm_temp = atm_temp_0 $ + DT2x * co2_doublings ocn_temp_target = atm_temp - 10. ocn_temp = ocn_temp $ + (ocn_temp_target-ocn_temp) $ / 1000. * time_step r_co2 = pco2(1) / 280. Solar_Const = -7.4 ! p32 RUN = 0.045 ! B+K2001, Fig 3, cold periods tau = DT2x Z = 0.09 f_bt_co2 = r_co2**(Z*tau) f_bt_co2 = f_bt_co2 $ * exp(z*GEOG + z*Solar_Const*(geol_year/570)) f_bt_co2 = f_bt_co2 $ * ( 1 $ + RUN * ( $ tau * log(r_co2) $ + Solar_Const * ( geol_year / 570 ) $ + GEOG $ ) $ )**0.65 if(f_plants(ip) .LT. 0.5) then f_bt_co2 = f_bt_co2 * 0.8 endif #ifdef LTMIP if(instring .NE. "csw") then f_bt_co2 = 1.0 endif #endif WEAT_CARB = weat_carb_k * f_bt_co2 #ifdef CarbonateLandAreaFactor $ * f_land(ip) ! mol/yr #endif WEAT_SIL = weat_sil_k * f_bt_co2 $ * f_land(ip) c gas exchange call calc_co2(ocn_conc(i_tc),ocn_conc(i_alk), $ atm_temp,pco2_sat,co3) GE = GEpPPM * (pco2_sat - pco2(1)) ! + up #ifdef Carbon13 pco2_13_eq = pco2_sat $ / ocn_conc(i_tc) * ocn_conc(i_tc_13) $ * ( 1.d0 - 7.d-3 ) c call calc_co2(ocn_conc(i_tc_13),ocn_conc(i_alk), c $ atm_temp,pco2_sat,co3) GE13 = GEpPPM * (pco2_13_eq - pco2(2)) ! + up #endif c caco3 precipitation BC = 3.E-3 ! from Archer neut_long $ * (co3-170.) ! offset between ocn mean and deep $ * 1.e3 / 12. ! e12 mol/yr f_caco3_T = 0. if(f_caco3_T .GT. 0.) then BC = BC * exp((atm_temp-atm_temp_0)*f_caco3_T) endif do is=1,3 ocn_inv(is) = ocn_conc(is) * ocn_vol ! mol enddo do is=1,2 atm_inv(is) = pco2(is) * MOLpPPM enddo #ifdef LTMIP if(instring .EQ. "___" .OR. instring .EQ. "c__") then WEAT_CARB = 0. BC = 0. endif if(instring .NE. "csw") then WEAT_SIL = 0. co2_degas = 0. endif #endif ocn_inv(i_alk) = ocn_inv(i_alk) $ + ( WEAT_CARB * 2. $ + WEAT_SIL * 2. $ - BC * 2. $ ) * time_step ocn_inv(i_tc) = ocn_inv(i_tc) $ + ( WEAT_CARB $ + WEAT_SIL $ - BC $ - GE $ ) * time_step #ifdef Carbon13 ocn_inv(i_tc_13) = ocn_inv(i_tc_13) $ + ( WEAT_CARB $ + WEAT_SIL $ - BC * ocn_inv(i_tc_13) / ocn_inv(i_tc) $ - GE13 $ ) * time_step #endif atm_inv(1) = atm_inv(1) $ + ( GE $ + co2_degas(ip) $ - WEAT_SIL $ ) * time_step #ifdef Carbon13 atm_inv(2) = atm_inv(2) $ + ( GE13 $ + co2_degas(ip) $ - WEAT_SIL $ ) * time_step #endif do is=1,3 ocn_conc(is) = ocn_inv(is) / ocn_vol enddo pco2(1) = atm_inv(1) / MOLpPPM pco2(2) = atm_inv(2) / MOLpPPM i_print = 0 if(ip .EQ. 1) then ! spinup period if(time .GT. -1000) then i_print = 1 elseif(time .GT. -1000) then if(MOD(it,100) .EQ. 0) then i_print = 1 ENDIF endif else if(time .LE. 1001.) then i_print = 1 elseif(time .LE. 25001.) then if(MOD(it,10).eq. 1) then i_print = 1 endif elseif(time .LE. 100001.) then if(MOD(it,100).eq. 1) then i_print = 1 endif elseif(MOD(it,1000) .EQ. 1) then i_print = 1 endif endif #ifdef LTMIP i_print = 0 if(ip .EQ. 1) then if(time .GT. -100.1) then i_print = 1 endif else if(time .LE. 1000) then i_print = 1 else if(time .LE. 10000 .AND. MOD(it,10) .EQ. 1) then i_print = 1 endif endif #endif c i_print = 1 if(i_print .EQ. 1) then write(6,'(f15.5,5f15.5,8f15.5,3f15.5)')time, ! 1 $ ocn_conc(i_tc),ocn_conc(i_alk), ! 2,3 $ (ocn_conc(3)/ocn_conc(2)-1)*1000., ! 4 $ pco2(1), $ (pco2(2)/pco2(1)-1)*1000., $ co3, ! 5,6 $ WEAT_CARB,WEAT_SIL, ! 7,8 $ WEAT_CARB+WEAT_SIL,BC, ! 9,10 $ co2_degas(ip), ! 11, 12 $ atm_temp,ocn_temp ! 13, 14 endif enddo ! time steps c#ifdef BoundaryCO2Slug atm_inv = atm_inv + co2_slug(:) pco2 = atm_inv / MOLpPPM c#endif enddo end subroutine calc_co2(tc,alk,temp,pco2,co3) implicit none real*8 tc,alk, k1,k2,khco2,temp,pco2,co3,tk,sal, $ co2,hco3 sal = 35. tk = temp + 273.15 c write(6,*) "in calc_co2", conc,temp,tk k1 = 13.7201 - 0.031334 * tk * - 3235.76 / tk * - 1.3E-5 * sal * tk * + 0.1032 * sal**(0.5) k1 = 10**(k1) k2 = - 5371.9645 * - 1.671221 * tk * + 128375.28 / tk * + 2194.3055 * LOG10( tk ) * - 0.22913 * sal * - 18.3802 * LOG10( sal ) * + 8.0944E-4 * sal * tk * + 5617.11 * LOG10( sal ) / tk * - 2.136 * sal / tk k2 = 10**(k2) c write(6,*) "before newt", conc,sal,temp,k1,k2 call newt(alk/1.e3, tc/1.e3, $ sal,temp,k1,k2, $ co2,hco3,co3) khco2 = exp ( * -60.2409 + 9345.17 / tk * + 23.3585 * log (tk/100.) * + sal * ( * 0.023517 - 2.3656e-4 * tk * + 4.7036e-7 * tk * tk * ) * ) PCO2 = CO2 / KHCO2 * 1.e6 co3 = co3 * 1.e6 c write(6,*) "after newt", co2,hco3,co3 c write(6,*) "end of newt", khco2,pco2,co3 return End #ifdef DamagedNewt SUBROUTINE NEWT(ALK,TCO2,SAL,TEMP, * K1,K2, * CO2,HCO3,CO3) implicit none REAL*8 K1,K2,KW,tbor,tkt,sal,fh,c1,c2,c4,ah1,wv,retard,a, $ alk,tco2,temp,tk,aht,wm,x,co2,hco3,co3 integer icnt TBOR = 0. TKT = TEMP+273 TK = TKT/100. KW = EXP( 148.9802 - 13847.26/tkt * - 23.6521*LOG(tkt) - 0.019813*sal * + sal**0.5 * ( * - 79.2447 + 3298.72/tkt * + 12.0408*LOG(tkt) * ) * ) FH = 1.29 - 0.00204*tkt + 4.6*1.E-4*sal**2 * - 1.48*1.E-6*sal**2*tkt C1 = K1/2.0 C2 = 1.0 - 4.0*K2/K1 C4 = 0. AHT= 1.E-8 DO 100 ICNT=1,100 WM = ( (KW*FH/AHT) - (AHT/FH) ) c A = ALK - BM - SiM - PM - WM retard = 3. ! 1 would be nothing A = ( (retard-1.)* A + ALK - WM) / retard c write(6,*) "in newt", icnt, a X = A/TCO2 AH1 = C1/X*(1.0-X+SQRT(1.0+C2*X*(-2.+X))) IF(0.5E-4.GE.ABS(1.-AHT/AH1)) GOTO 200 AHT=AH1 100 CONTINUE 200 CONTINUE CO3 = (A-TCO2)/(1.0-(AH1*AH1)/(K1*K2)) HCO3 = TCO2/(1.+AH1/K1+K2/AH1) CO2 = TCO2/(1.0+K1/AH1+K1*K2/(AH1*AH1)) RETURN END #endif SUBROUTINE NEWT(ALK,TCO2,SAL,TEMP, * K1,K2, * CO2,HCO3,CO3) implicit none REAL*8 K1,K2,KB,KSi,KP2,KP3,KW,TSi,TP,TBOR,alk,tco2,temp,sal, $ tkt,tk,fh,c1,c2,c4,aht,bm,sim,pm,wm,a,x,ah1,co2,hco3,co3 integer icnt kb=1.e-9 TSi = 0. TP = 0. TBOR = 4.106E-4*SAL/35. TKT = TEMP+273 TK = TKT/100. KSi = 1.E-10 KP2 = EXP( -9.039 - 1450./tkt) KP3 = EXP( 4.466 - 7276./tkt) KW = EXP( 148.9802 - 13847.26/tkt * - 23.6521*LOG(tkt) - 0.019813*sal * + sal**0.5 * ( * - 79.2447 + 3298.72/tkt * + 12.0408*LOG(tkt) * ) * ) FH = 1.29 - 0.00204*tkt + 4.6*1.E-4*sal**2 * - 1.48*1.E-6*sal**2*tkt C1 = K1/2.0 C2 = 1.0 - 4.0*K2/K1 C4 = TBOR*KB AHT= 1.E-8 DO 100 ICNT=1,100 BM = TBOR * KB / (AHT + KB) SiM = TSi * 4 * 1E-10 / (AHT + 4 * 1.E-10) PM = TP * (1 / ( * 1 + KP2/AHT + KP2*KP3 / AHT**2 * ) * + 2 / ( * 1 + AHT/KP2 + KP3/AHT * ) * + 3 / ( * 1 + AHT/KP3 + AHT**2 / (KP3*KP2) * ) * ) WM = ( (KW*FH/AHT) - (AHT/FH) ) A = ALK - BM - SiM - PM - WM c write(6,*) "in newt", icnt,a X = A/TCO2 AH1 = C1/X*(1.0-X+SQRT(1.0+C2*X*(-2.+X))) IF(0.5E-4.GE.ABS(1.-AHT/AH1)) GOTO 200 AHT=AH1 100 CONTINUE 200 CONTINUE CO3 = (A-TCO2)/(1.0-(AH1*AH1)/(K1*K2)) HCO3 = TCO2/(1.+AH1/K1+K2/AH1) CO2 = TCO2/(1.0+K1/AH1+K1*K2/(AH1*AH1)) RETURN END