subroutine coa (is, ie, js, je) implicit none integer ie, is, je, js, i, j, nc integer c real time_new &, tarea_earth, area_grid, tarea_coa &, sss,sst,dic_in,co2_in,pHlo,pHhi,pH,sit_in,pt_in,atmpres,zero &, co2star,dco2star,pco2,co3,hplus,ta_1,Omega_c,Omega_a &, conv, Ea_olivine, Rgas, kk, kk25, Rdiss0 &, molvol_olivine, pi, sectim, Rdiss, Rdiss1 &, dm1, Vm1, Mm1, Am1 &, dp1, Vp1, Mp1, Mp1tmp &, coa_tmp,coa_diss_c, count_coa_cells &, boundary_grid, tboundary_coa, time_start logical coa_debug include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "calendar.h" include "csbc.h" include "coord.h" include "grdvar.h" include "tmngr.h" include "switch.h" # if defined O_embm include "cembm.h" include "atm.h" # endif # if defined O_mom include "mw.h" # endif real ttmskxy(imtm2,jmtm2) !----------------------------------------------------------------------- ! coastal ocean alkalinization (COA) !----------------------------------------------------------------------- ! WKs COA code version; // V 20.02.2017 // wkoeve@geomar.de ! The overall approach of COA used here is that of adding a defined amount (batch) ! of olivine to coastal-ocean model grid cells at each segtim and follow this batch ! over time. The batch size may vary with time and location depending on the ! current state of the system, in particular the local omega with respect to ! the mineral aragonite. The code also allows to prescribe an annual global olivine addition. ! ! If you dont use WK's interactive version you can define/change certain value here: coastartyr=2020.0 ! deployment will start at the beginning of yr. 2021 coaendyr=2100.0 ! deployment will end at the end of yr 2100 coa_dia0=100.0 ! default diameter of olivine grains at time of deployment is 100 um omthresh=9.0 ! default omega threshold is 50 (i.e. a high value) coaaddratio=1.0 ! a tuning parameter introduced by Ell coa_do_debug=.true. ! if true you get a lot of additional output into a file coa.log olivine_total = -1. ! if set to a positive values (in Gt olivine!/yr) it forces a constant annual addition ! If you work outside WK's code version please initiate and communicate the following variables in ! common/csbc.h ! logical coa_do_debug ! common /csbc_l coa_do_debug ! integer cumm,coa_unit ! common /csbc_i/ cumm,coa_unit ! real omthresh,coastartyr,coaendyr,coaaddratio,coa_dia0 ! real olivine_total ! common /csbc_r/ omthresh,coastartyr,coaendyr,coaaddratio ! &, coa_dia0, olivine_total ! real coabatchdim,coa_omta,coa_cycles ! real coa_batch,coa_dia ! real coa_batch_cumm ! real coa_olivine,coa_diss ! real coa_grid ! real coa_fT_olivine,coa_fpH_olivine ! real coa_co2target ! common /csbc_r/ coabatchdim,coa_omta,coa_cycles ! &, coa_co2target ! &, coa_batch(imt,jmt,7300), coa_dia(imt,jmt,7300) ! &, coa_olivine(imt,jmt), coa_diss(imt,jmt) ! &, coa_batch_cumm(imt,jmt) ! &, coa_grid(imt,jmt) ! &, coa_fT_olivine(imt,jmt),coa_fpH_olivine(imt,jmt) ! ! Add the following lines to mom/setmom.F ! (before: -- do all consistency checks last -- ! coa_grid(:,:)=0. ! ! Finally call coa() from within gosbc.F, e.g., below the section ! //apply CaCO3 weathering flux proportional to discharge// ! and above the section ! //add in flux from sea level change// ! I.e. at approx line 288. ! ! Avoid changing parameters and code below this line. :-) --------------------- ! Some parameters needed: ! molvol_olivine is the molar volume of olivine (units: m3/mol) (value from Hangx and Spiers, 2009) molvol_olivine=43.02*1e-6 ! pi pi=3.14159 ! sectim is segtim converted to seconds sectim=segtim*24.*3600. ! Ea_olivine is the activation energy of olivine dissolution (J/mol) ! from Hangx and Spiers, 2009, Int. J. Greenhouse Gas Control, 3, 757-767) Ea_olivine = 79.5e3 ! Rgas, ideal age gas constant in cm2 bar K-1 Rgas = 83.1451 coa_unit=20 coa_debug=.false. coa_olivine(:,:) = 0. coa_diss(:,:) = 0. coa_fT_olivine(:,:) = 1. coa_fpH_olivine(:,:) = 1. time_new =year0 + accel_yr0 + (relyr - accel_yr0)*accel if(time_new.eq.coastartyr)then cumm=1 coa_cycles=1 elseif(time_new .gt. coastartyr .and. time_new .lt. coaendyr)then cumm=cumm+1 elseif(time_new .lt. coastartyr) then cumm=0 endif if (time_new .ge. coastartyr) then coa_cycles=coa_cycles + 1. tarea_earth = 0. do j=2,jmtm1 do i=2,imtm1 area_grid = dxt(i)*dyt(j)*cst(j) tarea_earth = tarea_earth + area_grid enddo enddo ! Currently, COA is only carried out north of j=21 and south of j=84. tarea_coa = 0. tboundary_coa = 0. do j=21,84 do i=2,100 ttmskxy(i,j)=tmsk(i+1,j)*tmsk(i,j+1)*tmsk(i-1,j)*tmsk(i,j-1) if ((tmsk(i,j) .ge. 0.5).and. & (ttmskxy(i,j) .eq. 0.))then area_grid = dxt(i)*dyt(j)*cst(j) tarea_coa = tarea_coa + area_grid coa_grid(i,j)=1. endif enddo enddo if (time_new .eq. coastartyr .and. coa_do_debug) then open(unit=coa_unit,file='coa.log',action='write') write (coa_unit,*) '--- gosbc.f, O_coa --------------' write(coa_unit,*) ' time_new: ', time_new write(coa_unit,*) ' coastartyr: ', coastartyr write(coa_unit,*) ' coaendyr: ', coaendyr write(coa_unit,*) ' coa_dia0: ', coa_dia0 write(coa_unit,*) ' omthresh: ', omthresh write(coa_unit,*) ' coaaddratio: ', coaaddratio write(coa_unit,*) & ' tarea_earth: ',tarea_earth, 'in cm2' write(coa_unit,*) & ' tarea_coa: ',tarea_coa, 'in cm2' write(coa_unit,*) & ' tboundary_coa: ', tboundary_coa,' in cm' write(coa_unit,*) & ' co2emit: ',co2emit, ' in g cm-2 s-1' write(coa_unit,*) & '-- starting coa deployment now ----' endif ! The following block deals with the definition of the batch size for the current cumm count_coa_cells=0. coa_debug=.false. do j=21,84 do i=2,100 if (j.eq.21 .and. i.eq.2 .and. coa_do_debug) then ! once each segtime ... write(coa_unit,*) & ' O_coa time_new is', time_new write(coa_unit,*) & " cumm is ", cumm write(coa_unit,*) & ' coaaddratio: ', coaaddratio endif ttmskxy(i,j)=tmsk(i+1,j)*tmsk(i,j+1)*tmsk(i-1,j)*tmsk(i,j-1) if ((tmsk(i,j) .ge. 0.5) .and. (ttmskxy(i,j) .eq. 0.)) then ! if we are in a coastal grid cell if (count_coa_cells .eq. 0. .and. coa_do_debug) & coa_debug=.true. if (count_coa_cells .ne. 0.) coa_debug=.false. count_coa_cells=count_coa_cells+1. sss = 1000.0*sbc(i,j,isss) + 35.0 sst = sbc(i,j,isst) sst = min(35.,max(sst,-2.)) sss = min(45.,max(sss,0.)) dic_in = sbc(i,j,issdic) co2_in = co2ccn pHlo = 6. pHhi = 10. sit_in = 7.6875e-03 !mol/m^3 pt_in = 0.5125e-3 !mol/m^3 atmpres = 1.0 !atm zero = 0. hplus = 8.2 ta_1 = sbc(i,j,issalk) if (ta_1/dic_in .gt. 1.5) then write(coa_unit,*) & 'ta_1/dic_in warning', coa_cycles, j, i, ta_1/dic_in endif call co2calc_SWS (sst, sss, dic_in, ta_1, co2_in, pt_in &, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, hplus, CO3 &, Omega_c, Omega_a) if (time_new .lt. coaendyr) then if (coa_debug) then write (coa_unit,*) & "in deployment period cumm is ", cumm endif !------ Olivine addition -------------------------------------------------- ! Currently, the default mode of olivine addition is related to co2emit, the relative ! area of coastal grid cells and a tuning factor (coaaddratio) which was ! originally hand tuned by EYF in order to restrict omega overshoots. ! This procedure is arbitrary and will be replaced in the future. ! In this version, however, the value of coaaddratio may be modified via the namelist. ! Other variables and parameters: ! - co2emit: anthropogenic CO2 emission, unit "g cm^-2 s^-1". ! - conv: unit conversion from "g" to "umol" ! - dzt(1): unit converion from /cm2 to /cm3; 5000 for a 50m thick surfaca layer ! - coa_batch(i,j,cumm) is the size of the currently added batch. ! (unit conversion: umol cm^-2 s^-1 * () * ()/() * () * s /cm) ! final unit: umol cm¯3 ! I keep memory of the decreasing batch-sizes over the model run. conv=1./12e-6 if(Omega_a .lt. omthresh) then if(olivine_total .lt. 0.) then coa_olivine(i,j)= co2emit*conv*tarea_earth/tarea_coa & * coaaddratio else ! olivine_total is given in Gt olivine /yr ! convert to /s (/365./86400. ! convert to g olivine (*1.e15) ! convert to /cm2 (/tarea_coa) ! convert to mol olivine (/153.31 (g/mol) ! convert to umol *1.e6 ! convert from olivine to TA units (*4) ! final unit umol TA cm-2 s-2 coa_olivine(i,j)=olivine_total/365./86400. & *1.e15/tarea_coa/153.31*1.e6*4. endif coa_batch(i,j,cumm)=coa_olivine(i,j)*sectim/dzt(1) if (coa_debug .and. cumm.eq.1) then write(coa_unit,*) 'coa_olivine(i,j): ', & coa_olivine(i,j) write(coa_unit,*) & 'coa_batch(i,j,cumm): ', coa_batch(i,j,cumm) endif ! coa_dia(i,j,cumm) stores the diameter of olivine particles (unit: m) ! For the shrinking core model I use 'm' units since constants used are in these units. ! Note: The diameter will change/shrink with time (s.b.) coa_dia(i,j,cumm)=coa_dia0/1.e6 else ! coa_olivine(i,j)=0. actually is predefined to be zero coa_batch(i,j,cumm)=0. coa_dia(i,j,cumm)=0. if (coa_debug .and. coa_do_debug) then write(coa_unit,*) & "Omega_a upper thresh has been passed" endif endif else if (coa_debug .and. coa_do_debug) then write(coa_unit,*) 'olivine deployment period is over' endif ! next is end of (time_new .lt. coaendyr) endif ! -------- Olivine dissolution --------------------------------------- ! We compute olivine dissolution from a shrinking core model, including a ! pH-sensitivity and a temperature sensitivity ! ! For the temperature dependence we use an Arrhenius approach ! fT_olivine is the temperature dependence of olivine dissolution, normalised to 25degC kk = exp(-Ea_olivine/(Rgas*(sst+273.15))) kk25 = exp(-Ea_olivine/(Rgas*(25.+273.15))) coa_fT_olivine(i,j) = kk/kk25 # if defined O_coa_nofT coa_fT_olivine(i,j) = 1.0 # endif ! ! For the pH sensitivity we use the truncated pH dependency from Wogelius and Walther, 1991 ! (Geochim Cosmoschim, 55, 943-954) normalized to pH=8.2 coa_fpH_olivine(i,j) = 0.9394+1.623*10.**(3.-0.54*pH) # if defined O_coa_nofpH coa_fpH_olivine(i,j) = 1.0 # endif ! ! Rdiss0 is the dissolution rate constant (for t=25, pH=8.2), taken from Hangx and Spiers, 2009 ! in mol/m2 olivine surface (!!) Rdiss0=1.58e-10 ! Rdiss is the T, pH affected dissolution rate (still in mol/m2 olivine surface) Rdiss = Rdiss0 * coa_fT_olivine(i,j) & * coa_fpH_olivine(i,j) do c=1,cumm ! for all batched already deployed ... ! ------- The shrinking core model --------------------------------------- ! dm1 is the recent diameter (unit: m) ! Vm1 is the recent volume of a characteristic olivine particle (units: m3) ! Mm1 is the mass of a characteristic olivine particle (units: mol) ! Am1 is the recent area of a characteristic olivine particle (units: m2) ! Rdiss1, dissolution rate in mol/s dm1=coa_dia(i,j,c) if (dm1 .gt. 0.) then !--------------- Below is Wolfgang's shrinking core model with explict treatment of olivine volume ! Vm1=1./6.*pi*dm1**3. ! Mm1=Vm1/molvol_olivine ! Am1=pi*dm1**2. ! Rdiss1=Rdiss*Am1 ! Mp1 is the mass of a characteristic olivine particle AFTER segtim ! Vp1 is the volume of a characteristic olivine particle AFTER segtim ! dp1 is the diameter of a characteristic olivine particle AFTER segtim (Units: m) ! if (Rdiss1*sectim .gt. Mm1) then ! Mp1 = 0. ! Vp1 = 0. ! dp1 = 0. ! coa_tmp = coa_batch(i,j,c) ! coa_batch(i,j,c) = 0. ! coa_dia(i,j,c) = 0. ! else ! Mp1tmp = Mm1 - Rdiss1*sectim ! Mp1 = max(0.,Mp1tmp) ! Vp1=Mp1*molvol_olivine ! dp1= (Vp1*6./pi)**(1./3.) ! store dp1 for next passage of gosbc.F ! coa_dia(i,j,c) = max(0.,dp1) ! update the batch size, we simple scale the batchsize by the volume ratio (new/old) ! we make sure that (Vp1 le 0) ! coa_batch(i,j,c)=coa_tmp * max(0.,Vp1)/Vm1 ! endif ! ------------------ End of Wolfgang's model part ! ------------------- Below is Ell's shrinking core model, with only explict treatment to diameter coa_dia(i,j,c)=dm1-molvol_olivine*Rdiss*2*sectim coa_dia(i,j,c) = max(0.,coa_dia(i,j,c)) coa_tmp=coa_batch(i,j,c) coa_batch(i,j,c)=coa_tmp * coa_dia(i,j,c)**3/dm1**3 ! -------------------- End of Ell's model part ! coa_diss is the olivine dissolution, simply the difference of the old and the ! new batch-size; batch size units umol cm^-3 are converted back to umol cm^-2 s^-1 coa_diss_c=(coa_tmp - coa_batch(i,j,c))*dzt(1)/sectim else coa_diss_c = 0. coa_batch(i,j,c) = 0. coa_dia(i,j,c)=0. endif if (coa_debug .and. c.eq.1 .and. & coa_diss_c .gt. 0.) then write (coa_unit,*) & 'diss section c= ',c &, 'coa_cycles=', coa_cycles write (coa_unit,'(A,E)') & 'dm1: ',dm1 write (coa_ unit,'(A,E)') & 'dp1: ',dp1 write (coa_unit,'(A,E)') & 'coa_dia(i,j,c) :',coa_dia(i,j,c) write (coa_unit,'(A,E)') & 'Rdiss1: ',Rdiss1 write (coa_unit,'(A,E)') & 'Vp1/Vm1 :',max(0.,Vp1)/Vm1 write (coa_unit,'(A,E)') & 'coa_diss :',coa_diss_c write (coa_unit,'(A,E)') & 'coa_batch :',coa_batch(i,j,c) write (coa_unit,'(A,E)') & 'fpH :', coa_fpH_olivine(i,j) write (coa_unit,'(A,E)') & 'fT :', coa_fT_olivine(i,j) ! write (coa_unit,'(A,E)') ! write (coa_unit,'(A,E)') endif ! Here we add the olivine dissolution of this batch to the vertial local alkalinity flux. ! Units of both sbc(..,flx) and coa_diss are identical (umol cm^-2 s^-1) coa_diss(i,j) = coa_diss(i,j) + coa_diss_c sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + coa_diss_c enddo ! next line is end of: ((tmsk(i,j) .ge. 0.5) .and. (ttmskxy(i,j) .eq. 0.)) endif enddo enddo ! next: if(time_new .ge. coastartyr.) endif ! print *, 'end of coa gosbc code' return end