MODULE module_ghg_fluxes USE module_configure USE module_state_description IMPLICIT NONE ! 04/05/2011- This module contains parameterizations to calculate biospheric greenhouse gas fluxes: ! CO2 uptake and release by VPRM model (done by Ravan Ahmadov - ravan.ahmadov@noaa.gov) ! CH4 flux modules: KAPLAN, SOILUPTAKE, TERMITE (done by Veronika Beck - vbeck@bgc-jena.mpg.de) ! Some references, where VPRM in WRF was used: ! 1) Ahmadov, R., C. Gerbig, et al. (2007). "Mesoscale covariance of transport and CO2 fluxes: ! Evidence from observations and simulations using the WRF-VPRM coupled atmosphere-biosphere model." JGR-Atmospheres 112(D22): 14. ! 2) Ahmadov, R., C. Gerbig, et al. (2009). "Comparing high resolution WRF-VPRM simulations and two global CO2 transport ! models with coastal tower measurements of CO2." Biogeosciences 6(5): 807-817. ! 3) Pillai, D., C. Gerbig, et al. (2011). "High-resolution simulations of atmospheric CO(2) over complex terrain - ! representing the Ochsenkopf mountain tall tower." Atmospheric Chemistry and Physics 11(15): 7445-7464. ! 4) Pillai, D., C. Gerbig, et al. (2009). "High resolution modeling of CO2 over Europe: implications for representation ! errors of satellite retrievals." Atmospheric Chemistry and Physics 10(1): 83-94. ! ! Information on the CH4 flux models are found in: ! 1) Kaplan, J. O., (2002): Wetlands at the Last Glacial Maximum: Distribution and methane ! emissions Geophys. Res. Lett., 29, No.6, 1079, doi:10.1029/2001GL013366. ! Ridgwell, A. J., S. J. Marshall, and K. Gregson, (1999): Consumption of atmospheric ! methane by soils: A process-based model Global Biochem. Cycles, 13(1), 59-70. ! Sanderson, M. G., (1996): Biomass of termites and their emissions of methane and carbon ! dioxode: A global database Global Biochem. Cycles, 10(4), 543-557. ! ! 2) References where the CH4 flux models in WRF are used: ! Beck, V., T. Koch, R. Kretschmer, J. Marshall, R. Ahmadov, C. Gerbig, D. Pillai, ! and M. Heimann, (2011): The WRF Greenhouse Gas Model (WRF-GHG). Technical ! Report No. 25, Max Planck Institute for Biogeochemistry, Jena, Germany. ! available online via: http://www.bgc-jena.mpg.de/bgc-systems/index.shtml ! Beck, V. et al., (2012): WRF-Chem simulations in the Amazon region during wet and ! dry season transitions: evaluation of methane models and wetland inundation maps, in preparation ! ! Reference for the WRF_chem module_ghg_fluxes ! Beck and Ahmadov et al., (2012): module_ghg_fluxes: A new module in WRF-CHEM for the passive tracer ! transport of greenhouse gases, (in prep.) ! CONTAINS SUBROUTINE add_ghg_fluxes ( ids,ide, jds,jde, kds,kde, & ! Domain dimensions ims,ime, jms,jme, kms,kme, & ! Memory dimensions its,ite, jts,jte, kts,kte, & ! Tile dimensions dtstep, dz8w, config_flags, rho_phy, & chem, emis_ant,eghg_bio,ebio_co2oce ) ! March-28, this subroutine adds all type of greenhouse gases to the chem species !IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, INTENT(IN ) :: dtstep REAL, DIMENSION( ims:ime,jms:jme ), INTENT(IN ) :: ebio_co2oce REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: rho_phy, dz8w REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme, num_emis_ant ), INTENT(IN ) :: emis_ant REAL, DIMENSION( ims:ime, 1,jms:jme, num_eghg_bio ), INTENT(IN ) :: eghg_bio INTEGER :: i,j,k REAL :: conv_rho call wrf_debug(15,'add_ghg_fluxes') ! For both GHG options DO j=jts,jte DO i=its,ite ! 3D anthropogenic fluxes DO k=kts,min(config_flags%kemit,kte) conv_rho=8.0461e-6/rho_phy(i,k,j)*dtstep/dz8w(i,k,j) ! 8.0461e-6=molar_mass(air)/3600, [g/mol/s] chem(i,k,j,p_co2_ant)= chem(i,k,j,p_co2_ant) + conv_rho* emis_ant(i,k,j,p_e_co2) chem(i,k,j,p_co2_tst)= chem(i,k,j,p_co2_tst) + conv_rho* emis_ant(i,k,j,p_e_co2tst) chem(i,k,j,p_co_ant) = chem(i,k,j,p_co_ant) + conv_rho* emis_ant(i,k,j,p_e_co) ! 2D biospheric fluxes: if (k==1) then chem(i,1,j,p_co2_bio)= chem(i,1,j,p_co2_bio) + conv_rho* (eghg_bio(i,1,j,p_ebio_gee) + eghg_bio(i,1,j,p_ebio_res)) ! both uptake and release chem(i,1,j,p_co2_oce)= chem(i,1,j,p_co2_oce) + conv_rho* ebio_co2oce(i,j) end if ENDDO ENDDO ENDDO ! For the GHG_TRACER option only IF(config_flags%chem_opt==GHG_TRACER) THEN DO j=jts,jte DO i=its,ite ! 3D anthropogenic fluxes DO k=kts,min(config_flags%kemit,kte) conv_rho=8.0461e-6/rho_phy(i,k,j)*dtstep/dz8w(i,k,j) ! 8.0461e-6=molar_mass(air)/3600, [g/mol/s] chem(i,k,j,p_ch4_ant)= chem(i,k,j,p_ch4_ant) + conv_rho* emis_ant(i,k,j,p_e_ch4) chem(i,k,j,p_ch4_tst)= chem(i,k,j,p_ch4_tst) + conv_rho* emis_ant(i,k,j,p_e_ch4tst) chem(i,k,j,p_co_tst) = chem(i,k,j,p_co_tst) + conv_rho* emis_ant(i,k,j,p_e_cotst) ! 2D biospheric fluxes: if (k==1) then chem(i,1,j,p_ch4_bio)= chem(i,1,j,p_ch4_bio) + conv_rho* (eghg_bio(i,1,j,p_ebio_ch4wet) + eghg_bio(i,1,j,p_ebio_ch4soil) & + eghg_bio(i,1,j,p_ebio_ch4term)) end if ENDDO ENDDO ENDDO END IF END SUBROUTINE add_ghg_fluxes !************************************************************************************************** SUBROUTINE VPRM ( ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & vprm_in,rad0, lambda, & alpha, RESP0, & T2,RAD, eghg_bio ) !IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte INTEGER :: i,j,m REAL, PARAMETER :: const= 3.6e+3 ! For unit conversion from mmol/m^2/s to mol/km2/hr REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: T2, RAD REAL :: PARadiation ! Hilton use PAR, consistent with rad0 REAL, DIMENSION( ims:ime, 8, jms:jme, num_vprm_in), INTENT(IN) :: vprm_in REAL, DIMENSION( ims:ime, 1,jms:jme, num_eghg_bio ), INTENT(INOUT ) :: eghg_bio REAL, DIMENSION(num_vprm_in) :: rad0, lambda, alpha, RESP0 REAL :: a1,a2,a3,Tair,Tscale,Wscale,Pscale,GEE_frac,RESP_frac,evithresh,RADscale ! VPRM vegetation classes: !1-Evergreen c !2-Deciduous !3-Mixed forest !4-Shrubland !5-Savanna !6-Cropland !7-Grassland !8-Others REAL, DIMENSION(8) :: Tmax,Tmin,Topt ! These are universal VPRM parameters DATA Tmin /0.,0.,0.,2.,2.,5.,2.,0./ DATA Topt /20.,20.,20.,20.,20.,22.,18.,0./ DATA Tmax /8*40./ DO j=jts,min(jte,jde-1) DO i=its,min(ite,ide-1) GEE_frac= 0. RESP_frac= 0. Tair= T2(i,j)-273.15 veg_frac_loop: DO m=1,7 if (vprm_in(i,m,j,p_vegfra_vprm)<1.e-8) CYCLE ! Then fluxes are zero a1= Tair-Tmin(m) a2= Tair-Tmax(m) a3= Tair-Topt(m) ! Here a1 or a2 can't be negative if (a1<0. .OR. a2>0.) then Tscale= 0. else Tscale=a1*a2/(a1*a2 - a3**2) end if if (Tscale<0.) then Tscale=0. end if ! modification due to different dependency on ground water if (m==4 .OR. m==7) then ! grassland and shrubland are xeric systems ! if (vprm_in(i,m,j,p_lswi_max)<1e-7 .or. abs(vprm_in(i,m,j,p_lswi_max)-vprm_in(i,m,j,p_lswi_min)).lt.1e-7 ) then ! in order to avoid NaN for Wscale if (vprm_in(i,m,j,p_lswi_max)<1e-7 .or. abs(vprm_in(i,m,j,p_lswi_max)-vprm_in(i,m,j,p_lswi_min)).lt.0.0000327 ) then ! in order to avoid NaN for Wscale Wscale= 0. else Wscale= (vprm_in(i,m,j,p_lswi)-vprm_in(i,m,j,p_lswi_min))/(vprm_in(i,m,j,p_lswi_max)-vprm_in(i,m,j,p_lswi_min)) end if else Wscale= (1.+vprm_in(i,m,j,p_lswi))/(1.+vprm_in(i,m,j,p_lswi_max)) end if ! effect of leaf phenology if (m==1) then ! evegreen Pscale= 1. else if (m==5 .OR. m==7) then ! savanna or grassland Pscale= (1.+vprm_in(i,m,j,p_lswi))/2. else ! Other vegetation types evithresh= vprm_in(i,m,j,p_evi_min) + 0.55*(vprm_in(i,m,j,p_evi_max)-vprm_in(i,m,j,p_evi_min)) if (vprm_in(i,m,j,p_evi)>=evithresh) then ! Full canopy period Pscale= 1. else Pscale=(1.+vprm_in(i,m,j,p_lswi))/2. ! bad-burst to full canopy period end if end if PARadiation = RAD(i,j)/0.505 RADscale= 1./(1. + PARadiation/rad0(m)) GEE_frac= lambda(m)*Tscale*Pscale*Wscale*RADscale* vprm_in(i,m,j,p_evi)* PARadiation*vprm_in(i,m,j,p_vegfra_vprm) + GEE_frac RESP_frac= (alpha(m)*Tair + RESP0(m))*vprm_in(i,m,j,p_vegfra_vprm) + RESP_frac ENDDO veg_frac_loop eghg_bio(i,1,j,p_ebio_gee)= min(0.0,const*GEE_frac) eghg_bio(i,1,j,p_ebio_res)= max(0.0,const*RESP_frac) ENDDO ENDDO END SUBROUTINE VPRM !**************************************************************************************** !VB: here comes the subroutine for the kaplan wetland inventory (Kaplan, 2002+2006) !which calculates CH4 wetland emissions from wrf soil temperature and soil moisture fields !and a global carbon density and a potential wetland map (see also Sitch et al. 2003) for details SUBROUTINE KAPLAN ( ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & dt, T_soil, Soil_M, wet_in, & soil_type, T_skin, eghg_bio, & num_soil_layers,E_f,M_s ) !IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & num_soil_layers INTEGER :: i,j INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: soil_type REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr REAL, DIMENSION (ims:ime,1,jms:jme, num_wet_in), INTENT(IN) :: wet_in REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: T_skin REAL, DIMENSION (ims:ime, num_soil_layers, jms:jme), INTENT(IN) :: Soil_M, T_soil REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) :: eghg_bio REAL, INTENT(IN) :: dt,M_s,E_f REAL, DIMENSION( ims:ime, jms:jme) :: sm_01, sm_02, sm_03, sm_04, g_T, f_SM, & k_r, Carbon, HR, T_flood, sm_tot, T_2s, & T_peat, P_l, soil_sat REAL :: tau_10 ! Description of the variables used in the subroutine KAPLAN ! T_soil: soil temperature from WRF ! C_pool: carbon content of the fast carbon pool at time 0 - from the LPJ model 2000 ! Wet_map: potential wetland map of Jed Kaplan gives the percentage of wetland area ! covered by each grid cell ! ebio_ch4wet: methane wetland emissions - output from kaplan model ! Soil_M: 4(3) dimensional wrf soil moisture output ! sm_01: first wrf soil moisture layer ! sm_02: second wrf soil moisture layer ! sm_03: third wrf soil moisture layer ! sm_04: fourth wrf soil moisture layer ! g_T: gives a modified Arrhenius temperature dependence ! f_SM: factor accounting for the soil moisture ! k_r: decomposition rate ! C_x: carbon pool after a certain time ! HR: heteorotrophic respiration ! T_flood: methane emissions for the flooding part - not peat land ! sm_tot: the actual soil moisture used for the wetland emissions is built up out of the wrf soil moisture ! T_2s: temperature of first soil layer from surface used for calculation of wetland emissions ! E_f: correction factor accounting for the peatland wetland emissions ! T_peat: Ch4 wetland emissions from peat ground ! P_l: weighting factor for the global peatland and flooded wetland emissions ! T_a: mean annual surface temperature used for the weighting factor ! but this mean annual temperature has to come from another source (e.g. GEOS) ! dt: timestep of the wrf grid in seconds ! soil_type: soil type for getting saturated soil moisture ! soil_sat: value of saturated soil moisture ! T_skin: skin temperature to use instead of soil temperature ! Initialization of the constant variables tau_10 = 2.86 !litter soil pool DO j=jts,min(jte,jde-1) DO i=its,min(ite,ide-1) ! zero initialization of the fields to avoid crashing of the program eghg_bio(i,1,j,p_ebio_ch4wet)=0. HR(i,j)=0. k_r(i,j)=0. g_T(i,j)=0. T_flood(i,j)=0. T_peat(i,j)=0. P_l(i,j)=0. f_SM(i,j)=0. sm_01(i,j)=0. sm_02(i,j)=0. sm_03(i,j)=0. sm_04(i,j)=0. Carbon(i,j)=0. sm_tot(i,j)=0. T_2s(i,j)=0. ! If loop to derive the saturated soil moisture for the specific soil type IF(soil_type(i,j) == 1)then soil_sat(i,j) = 0.339 ELSE IF(soil_type(i,j) == 2) then soil_sat(i,j) = 0.421 ELSE IF(soil_type(i,j) == 3) then soil_sat(i,j) = 0.434 ELSE IF(soil_type(i,j) == 4) then soil_sat(i,j) = 0.476 ELSE IF(soil_type(i,j) == 5) then soil_sat(i,j) = 0.476 ELSE IF(soil_type(i,j) == 6) then soil_sat(i,j) = 0.439 ELSE IF(soil_type(i,j) == 7) then soil_sat(i,j) = 0.404 ELSE IF(soil_type(i,j) == 8) then soil_sat(i,j) = 0.464 ELSE IF(soil_type(i,j) == 9) then soil_sat(i,j) = 0.465 ELSE IF(soil_type(i,j) == 10) then soil_sat(i,j) = 0.406 ELSE IF(soil_type(i,j) == 11) then soil_sat(i,j) = 0.468 ELSE IF(soil_type(i,j) == 12) then soil_sat(i,j) = 0.468 ELSE IF(soil_type(i,j) == 13) then soil_sat(i,j) = 0.439 ELSE IF(soil_type(i,j) == 14) then soil_sat(i,j) = 1.0 ELSE IF(soil_type(i,j) == 15) then soil_sat(i,j) = 0.20 ELSE IF(soil_type(i,j) == 16) then soil_sat(i,j) = 0.421 ELSE IF(soil_type(i,j) == 17) then soil_sat(i,j) = 0.468 ELSE IF(soil_type(i,j) == 18) then soil_sat(i,j) = 0.200 ELSE IF(soil_type(i,j) == 19) then soil_sat(i,j) = 0.339 ELSE soil_sat(i, j) = 0.4 END IF ! Use mean value of first and 2nd soil layer (up to 30cm depth) ! and the surface soil temperature sm_01(i,j) = Soil_M(i,1,j) sm_02(i,j) = Soil_M(i,2,j) sm_03(i,j) = Soil_M(i,3,j) sm_04(i,j) = Soil_M(i,4,j) sm_tot(i,j) = 0.5*(sm_01(i,j)+sm_02(i,j)) IF(sm_tot(i,j) < 0.0) then sm_tot(i,j) = 0.0 END IF IF ((T_Soil(i,1,j) == 273.16).OR.(T_Soil(i,1,j) == 0.0)) then T_2s(i, j) = T_skin(i, j) - 273.15 ELSE T_2s(i,j)= T_Soil(i,1,j)-273.15 END IF ! calculate temperature dependence all based on sitch (2003) g_T(i,j) = EXP(308.56*(1.0/56.02 - 1.0/(T_2s(i,j) + 46.02))) ! calculate soil moisture factor ! soil moisture should be sm/saturated sm (Folley) f_SM(i,j) = 0.25 + 0.75 * (sm_tot(i,j)/soil_sat(i,j)) ! calculate decomposition rate per year/per hr 20.04.10 k_r(i,j) = ((1/tau_10)*g_T(i,j)*f_SM(i,j))/(12.0*24.0*30.0) ! calculate the variable carbon fast pool depending on a month Carbon(i,j) = wet_in(i,1,j,p_Cpool)*(1.0 - EXP(-k_r(i,j)*1.0)) ! linearize and get the decomposed carbon for one second Carbon(i,j) = Carbon(i,j)/(3600.0) ! calculate the heteorotrophic respiration HR(i,j) = 0.7*Carbon(i,j) ! calculate the wetland emissions per grid cell for the flooding part T_flood(i,j) = HR(i,j)*M_s*wet_in(i,1,j,p_wetmap) ! calculate wetland emissions for the peatland part (higher latitudes) T_peat(i,j) = HR(i,j)*E_f*wet_in(i,1,j,p_wetmap) ! the factor for the relation between flood and peat land P_l(i,j) = EXP((wet_in(i,1,j,p_T_ann) - 303.0)/8.0) eghg_bio(i,1,j,p_ebio_ch4wet) = P_l(i,j)*T_flood(i,j)+(1.0 - P_l(i,j))*T_peat(i,j) ! units eghg_bio(i,1,j,p_ebio_ch4wet) = eghg_bio(i,1,j,p_ebio_ch4wet)*const END DO END DO END SUBROUTINE KAPLAN !*************************************************************************************** SUBROUTINE TERMITE ( ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & dt,eghg_bio, vtype_ter, & biomt_par, emit_par ) !IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte INTEGER :: i,j,vtype INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: vtype_ter REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) :: eghg_bio REAL, INTENT(IN) :: dt REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr REAL, DIMENSION(14), INTENT(IN) :: biomt_par REAL, DIMENSION(14), INTENT(IN) :: emit_par CHARACTER*256 :: mminlu_loc DO j=jts,min(jte,jde-1) DO i=its,min(ite,ide-1) CALL nl_get_mminlu(1 , mminlu_loc) IF(mminlu_loc .EQ. 'USGS')THEN vtype = INT(vtype_ter(i,j)) SELECT CASE (vtype) CASE (1) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (2) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10) CASE (3) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10) CASE (4) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11) CASE (5) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11) CASE (6) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11) CASE (7) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(7)*emit_par(7) CASE (8) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12) CASE (9) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12) CASE (10) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(5)*emit_par(5) CASE (11) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (12) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (13) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(1)*emit_par(1) CASE (14) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(2)*emit_par(2) CASE (15) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (16) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (17) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (18) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (19) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(13)*emit_par(13) CASE (20) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (21) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (22) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (23) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (24) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE default eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 END SELECT ELSE IF(mminlu_loc .EQ. 'MODIFIED_IGBP_MODIS_NOAH')THEN vtype = INT(vtype_ter(i,j)) SELECT CASE (vtype) CASE (1) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(2)*emit_par(2) CASE (2) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(1)*emit_par(1) CASE (3) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (4) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (5) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4) CASE (6) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12) CASE (7) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12) CASE (8) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(6)*emit_par(6) CASE (9) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(5)*emit_par(5) CASE (10) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(7)*emit_par(7) CASE (11) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (12) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10) CASE (13) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (14) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11) CASE (15) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (16) eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(13)*emit_par(13) CASE (17) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (18) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (19) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE (20) eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 CASE default eghg_bio(i,1,j,p_ebio_ch4term) = 0.0 END SELECT END IF eghg_bio(i,1,j,p_ebio_ch4term) = eghg_bio(i,1,j,p_ebio_ch4term)*const END DO END DO END SUBROUTINE TERMITE !************************************************************ !subroutine for soil uptake following ridgwell et al. 1999 !************************************************************ SUBROUTINE SOILUPTAKE ( ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & Soil_M, soil_typ, eghg_bio, & rain_1, rain_2, & potevap, sfevap, landuse, T2, dt, & num_soil_layers, wet_in ) !IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, & ims,ime, jms,jme, & its,ite, jts,jte, & num_soil_layers INTEGER :: i,j INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: soil_typ REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN) :: rain_1, T2, & rain_2, landuse, & potevap, sfevap REAL, DIMENSION (ims:ime, num_soil_layers, jms:jme), INTENT(IN) :: Soil_M REAL, INTENT(IN) :: dt REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) :: eghg_bio REAL, DIMENSION (ims:ime, 1,jms:jme, num_wet_in), INTENT(IN) :: wet_in REAL, DIMENSION(ims:ime, jms:jme) :: sm_01, sm_02, sm_03, sm_04, & i_cult, r_n, r_sm, vero, k_d, & G_soil, G_t, sm_res, dch4, & phi_soil, eps, b_f, sand_c, & clay_c, vero1 REAL :: dch4_0, k_0, z_dsoil, conv_Fk, conv_ch4, methane_0 CHARACTER*256 :: mminlu_loc !Initialization of the constants dch4_0 = 0.196 !cm2/s diffusivity of CH4 in free air k_0 = 0.00087 !1/s base oxidation rate constant for uncultivated !moist soil at 0 C obtained from experiments z_dsoil = 6.0 ! soil depth in cm at which the oxidation activity is assumed conv_Fk = 616.9 !mg/ppmv cm CH4 conversion factor fpr giving soil_flux right units of kgCH4/m^2s conv_ch4 = (28.97/26.0)*1e6 !ppm to kgCH4/kgair ! methane_0 = 9.66*1e-7 ! 1750ppb CH4 global average concentration methane_0 = 1.77 ! 1770ppb CH4 global average concentration DO j=jts,min(jte,jde-1) DO i=its,min(ite,ide-1) ! initialization of all variables sm_01(i, j) = 0.0 sm_02(i, j) = 0.0 sm_03(i, j) = 0.0 sm_04(i, j) = 0.0 i_cult(i, j) = 0.0 r_n(i, j) = 0.0 r_sm(i, j) = 0.0 vero(i, j) = 0.0 k_d(i, j) = 0.0 G_soil(i, j) = 0.0 G_t(i, j) = 0.0 sm_res(i, j) = 0.0 dch4(i, j) = 0.0 phi_soil(i, j) = 0.0 eps(i, j) = 0.0 b_f(i, j) = 0.0 sand_c(i, j) = 0.0 clay_c(i, j) = 0.0 vero1(i, j) = 0.0 ! end initialization sm_01(i,j) = Soil_M(i,1,j) sm_02(i,j) = Soil_M(i,2,j) sm_03(i,j) = Soil_M(i,3,j) sm_04(i,j) = Soil_M(i,4,j) !---------------------------------------------------------------------- ! calculation of CH4 diffusion activity !---------------------------------------------------------------------- G_t(i, j) = 1.0 + 0.0055 * (T2(i, j)-273.15) !************************ ! here the values for air filled porosity and total por volume have to be filled in !************************ ! calculation of the soil parameters air filled porosity, total pore volume and b factor ! depending on the soil type ! and also clay/sand content - from WRF SOILPARAM.TBL and triangle plot ! sand and clay content and porosity from Cosby et al. 1984 IF (soil_typ(i, j) == 1) then !sand phi_soil(i, j) = 0.339 eps(i, j) = 0.103 sand_c(i,j) = 0.92 clay_c(i,j) = 0.03 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 2) then !loamy sand phi_soil(i, j) = 0.421 eps(i, j) = 0.038 sand_c(i,j) = 0.82 clay_c(i,j) = 0.06 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 3) then !sandy loam phi_soil(i, j) = 0.434 eps(i, j) = 0.051 sand_c(i,j) = 0.58 clay_c(i,j) = 0.10 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 4) then ! silt loam phi_soil(i, j) = 0.476 eps(i, j) = 0.116 sand_c(i,j) = 0.17 clay_c(i,j) = 0.13 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 5) then !silt phi_soil(i, j) = 0.476 eps(i, j) = 0.093 sand_c(i,j) = 0.10 clay_c(i,j) = 0.10 !estimated no values b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 6) then !loam phi_soil(i, j) = 0.439 eps(i, j) = 0.11 sand_c(i,j) = 0.43 clay_c(i,j) = 0.18 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 7) then !sandy clay loam phi_soil(i, j) = 0.464 eps(i, j) = 0.09 sand_c(i,j) = 0.58 clay_c(i,j) = 0.27 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 8) then !silty clay loam phi_soil(i, j) = 0.404 eps(i, j) = 0.077 sand_c(i,j) = 0.10 clay_c(i,j) = 0.34 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 9) then !clay loam phi_soil(i, j) = 0.465 eps(i, j) = 0.083 sand_c(i,j) = 0.32 clay_c(i,j) = 0.34 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 10) then !sandy clay phi_soil(i, j) = 0.406 eps(i, j) = 0.068 sand_c(i,j) = 0.52 clay_c(i,j) = 0.42 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 11) then !silty clay phi_soil(i, j) = 0.468 eps(i, j) = 0.064 sand_c(i,j) = 0.06 clay_c(i,j) = 0.47 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 12) then !clay phi_soil(i, j) = 0.468 eps(i, j) = 0.056 sand_c(i,j) = 0.22 clay_c(i,j) = 0.58 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 13) then !organic material phi_soil(i, j) = 0.439 eps(i, j) = 0.11 sand_c(i,j) = 0.05 clay_c(i,j) = 0.05 !estimated b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 14) then !water phi_soil(i, j) = 1.0 eps(i, j) = 0.0 sand_c(i,j) = 0.60 clay_c(i,j) = 0.40 !estimated b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 15) then !bedrock phi_soil(i, j) = 0.2 eps(i, j) = 0.03 sand_c(i,j) = 0.07 clay_c(i,j) = 0.06 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 16) then !other phi_soil(i, j) = 0.421 eps(i, j) = 0.138 sand_c(i,j) = 0.25 clay_c(i,j) = 0.25 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 17) then ! playa phi_soil(i, j) = 0.468 eps(i, j) = 0.014 sand_c(i,j) = 0.60 clay_c(i,j) = 0.40 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 18) then !lava phi_soil(i, j) = 0.468 eps(i, j) = 0.03 sand_c(i,j) = 0.52 clay_c(i,j) = 0.48 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE IF (soil_typ(i, j) == 19) then !white sand phi_soil(i, j) = 0.339 eps(i, j) = 0.103 sand_c(i,j) = 0.92 clay_c(i,j) = 0.03 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) ELSE phi_soil(i, j) = 0.339 eps(i, j) = 0.0 sand_c(i,j) = 0.92 clay_c(i,j) = 0.03 b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5 & *sand_c(i, j)**2 *clay_c(i, j) END IF G_soil(i, j) = (phi_soil(i, j)**(4.0/3.0))* & ((eps(i, j)/phi_soil(i, j))**(1.5 + 3.0/b_f(i, j))) ! Calculate the Soil CH4 Diffusivity dch4(i, j) = G_soil(i, j)*G_t(i, j)*dch4_0 !-------------------------------------------------------------! ! Calculation of CH4 Oxidation Activity ! !-------------------------------------------------------------! ! Calculation of the temperature response to ch4 oxidation rate IF(T2(i,j)<= 0.0) then vero(i,j) = 0.0 ELSE IF(T2(i,j) > 0.0) then vero1(i, j) = (T2(i, j)-273.15)**4.0 vero(i,j) = EXP((0.0693*(T2(i,j)-273.15))-(8.56*10e-7*vero1(i, j))) END IF ! Calculation of the cultivation response CALL nl_get_mminlu(1, mminlu_loc) IF(mminlu_loc .EQ. 'USGS') THEN IF(landuse(i, j) <= 5.0) THEN i_cult(i, j) = 1.0 ELSE IF(landuse(i, j) == 6.0) THEN i_cult(i, j) = 0.5 ELSE i_cult(i, j) = 0.0 END IF ELSE IF(mminlu_loc .EQ. 'MODIFIED_IGBP_MODIS_NOAH')THEN IF(landuse(i, j) == 12.0 .OR. landuse(i,j) == 32.0 .OR. landuse(i,j) == 33.0 .OR. landuse(i,j) == 13.0) THEN i_cult(i, j) = 1.0 ELSE IF(landuse(i, j) == 14.0 .OR. landuse(i,j) == 31.0) THEN i_cult(i, j) = 0.5 ELSE i_cult(i, j) = 0.0 END IF END IF r_n(i, j) = 1.0 - (0.75*i_cult(i, j)) ! Calculation of soil moisture response !precipitation has to be monthly (*30) ! soil moisture 2nd layer (30cm) IF(sfevap(i, j) == 0.0) THEN sm_res(i, j) = 2.0 ELSE sm_res(i, j) = ((rain_1(i, j) + rain_2(i, j))*30.0*24.*3600./dt*0.001 & + sm_02(i, j))/sfevap(i, j) END IF IF(sm_res(i, j) > 1.0)THEN r_sm(i, j) = 1.0 ELSE IF(sm_res(i, j) <= 1.0)THEN r_sm(i, j) = sm_res(i, j) END IF k_d(i, j) = r_n(i, j)*r_sm(i, j)*vero(i, j)*k_0 !****************************************** ! Total soil uptake CH4 !****************************************** ! this gives you the flux in mg CH4/m^2 day IF(dch4(i,j) /= 0.0 .AND. wet_in(i,1,j,p_wetmap) < 0.10 .AND. soil_typ(i,j) /= 14 .AND. soil_typ(i,j) /=1) then eghg_bio(i,1,j,p_ebio_ch4soil) = ((methane_0*dch4(i, j))/z_dsoil)* & (1.0 - dch4(i, j)/(dch4(i, j)+k_d(i, j)*z_dsoil))*conv_Fk ! flux in mol/km2/hour eghg_bio(i,1,j,p_ebio_ch4soil) = -eghg_bio(i,1,j,p_ebio_ch4soil) * 1.0e-6/(24.0*3600.0)*const !negative soil flux ELSE eghg_bio(i,1,j,p_ebio_ch4soil) = 0.0 END IF IF(eghg_bio(i,1,j,p_ebio_ch4soil) > 0.0) then eghg_bio(i,1,j,p_ebio_ch4soil) = 0.0 END IF END DO END DO END SUBROUTINE SOILUPTAKE END MODULE module_ghg_fluxes