c#define NoWeatheringFeedback (crashes the world) #define CarbonateLandAreaFactor program geocarb implicit none integer np, n_record parameter(np=2) integer i_alk,i_tc,it,ip,nsteps,is,i_print real*8 co2_weat(np),f_land(np),f_plants(np),co2_degas(np) real*8 co2_slug(np), MOLpPPM, MOLpMM3 real*8 ocn_conc(2), ocn_inv(2) real*8 mean_latitude, geol_year, dt2x,GEOG,weat_carb_k,weat_sil_k, $ time_step,ocn_vol,GEpPPM,pco2,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,GE 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 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 do ip=1,np co2_slug(ip) = co2_slug(ip) * 1000. / 12. ! gton c to e12 mol enddo c parameters weat_carb_k = .00261 * 3000. weat_sil_k = 7. c DT2x = 3. time_step = 50. ! years MOLpPPM = 2. * 1.e3 / 12. ocn_vol = 3.e2 * 4000 ! e12 m^3 MOLpMM3 = ocn_vol GEpPPM = 1./300. * MOLpPPM ! mol/yr i_alk = 1 i_tc = 2 c initial condition pco2 = 200. 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,ocn_temp_0,pco2_test,co3) if(pco2_test .GT. pco2) then ocn_tco2_top = ocn_conc(i_tc) else ocn_tco2_bot = ocn_conc(i_tc) endif enddo time = 0. do ip = 1,np if(ip .EQ. 1) then nsteps = 5.e6 / time_step else nsteps = 2.0e6 / time_step endif do it = 1,nsteps time = (ip-2)*nsteps*time_step $ + (it-1)*time_step co2_doublings = LOG(pco2/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 / 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/579)) 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 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,atm_temp,pco2_sat,co3) GE = GEpPPM * (pco2_sat - pco2) ! + up c caco3 precipitation BC = 3.E-3 ! from Archer neut_long $ * (co3-150.) ! 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,2 ocn_inv(is) = ocn_conc(is) * ocn_vol ! mol enddo atm_inv = pco2 * MOLpPPM 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 atm_inv = atm_inv $ + ( GE $ + co2_degas(ip) $ - WEAT_SIL $ ) * time_step do is=1,2 ocn_conc(is) = ocn_inv(is) / ocn_vol enddo pco2 = atm_inv / MOLpPPM i_print = 0 if(ip .EQ. 1) then ! spinup period if(time .GT. -10000) then i_print = 1 elseif(time .GT. -10000) 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 if(i_print .EQ. 1) then write(6,'(f15.5,4f12.3,7f8.3,3f12.3)')time, ! 1 $ ocn_conc(1),ocn_conc(2), ! 2,3 $ pco2,co3, ! 4,5 $ WEAT_CARB,WEAT_SIL, ! 6,7 $ WEAT_CARB+WEAT_SIL,BC, ! 8,9 $ co2_degas(ip), ! 10 $ atm_temp,ocn_temp ! 11,12 endif enddo ! time steps atm_inv = atm_inv + co2_slug(ip) pco2 = atm_inv / MOLpPPM enddo end subroutine calc_co2(conc,temp,pco2,co3) implicit none real*8 conc(2), 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(conc(1)/1.e3, conc(2)/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