chem/module_ghg_fluxes.F

   1 MODULE module_ghg_fluxes
   2    USE module_configure
   3    USE module_state_description
   4
   5    IMPLICIT NONE
   6
   7 ! 04/05/2011- This module contains parameterizations to calculate biospheric greenhouse gas fluxes:
   8 ! CO2 uptake and release by VPRM model (done by Ravan Ahmadov - ravan.ahmadov@noaa.gov)
   9 ! CH4 flux modules: KAPLAN, SOILUPTAKE, TERMITE (done by Veronika Beck - vbeck@bgc-jena.mpg.de)
  10 ! Some references, where VPRM in WRF was used:
  11 ! 1) Ahmadov, R., C. Gerbig, et al. (2007). "Mesoscale covariance of transport and CO2 fluxes:
  12 ! Evidence from observations and simulations using the WRF-VPRM coupled atmosphere-biosphere model." JGR-Atmospheres 112(D22): 14.
  13 ! 2) Ahmadov, R., C. Gerbig, et al. (2009). "Comparing high resolution WRF-VPRM simulations and two global CO2 transport
  14 ! models with coastal tower measurements of CO2." Biogeosciences 6(5): 807-817.
  15 ! 3) Pillai, D., C. Gerbig, et al. (2011). "High-resolution simulations of atmospheric CO(2) over complex terrain -
  16 ! representing the Ochsenkopf mountain tall tower." Atmospheric Chemistry and Physics 11(15): 7445-7464.
  17 ! 4) Pillai, D., C. Gerbig, et al. (2009). "High resolution modeling of CO2 over Europe: implications for representation
  18 ! errors of satellite retrievals." Atmospheric Chemistry and Physics 10(1): 83-94.
  19 !
  20 ! Information on the CH4 flux models are found in:
  21 ! 1) Kaplan, J. O., (2002): Wetlands at the Last Glacial Maximum: Distribution and methane
  22 ! emissions Geophys. Res. Lett., 29, No.6, 1079, doi:10.1029/2001GL013366.
  23 ! Ridgwell, A. J., S. J. Marshall, and K. Gregson, (1999): Consumption of atmospheric
  24 ! methane by soils: A process-based model Global Biochem. Cycles, 13(1), 59-70.
  25 ! Sanderson, M. G., (1996): Biomass of termites and their emissions of methane and carbon
  26 ! dioxode: A global database Global Biochem. Cycles, 10(4), 543-557.
  27 !
  28 ! 2) References where the CH4 flux models in WRF are used:
  29 ! Beck, V., T. Koch, R. Kretschmer, J. Marshall, R. Ahmadov, C. Gerbig, D. Pillai,
  30 ! and M. Heimann, (2011): The WRF Greenhouse Gas Model (WRF-GHG). Technical
  31 ! Report No. 25, Max Planck Institute for Biogeochemistry, Jena, Germany.
  32 ! available online via: http://www.bgc-jena.mpg.de/bgc-systems/index.shtml
  33 ! Beck, V. et al., (2012): WRF-Chem simulations in the Amazon region during wet and
  34 ! dry season transitions: evaluation of methane models and wetland inundation maps, in preparation
  35 !
  36 ! Reference for the WRF_chem module_ghg_fluxes
  37 ! Beck and Ahmadov et al., (2012): module_ghg_fluxes: A new module in WRF-CHEM for the passive tracer
  38 ! transport of greenhouse gases, (in prep.)
  39 !
  40 CONTAINS
  41
  42 SUBROUTINE add_ghg_fluxes (  ids,ide, jds,jde, kds,kde,             &     ! Domain dimensions
  43                              ims,ime, jms,jme, kms,kme,             &     ! Memory dimensions
  44                              its,ite, jts,jte, kts,kte,             &     ! Tile dimensions
  45
  46                              dtstep, dz8w, config_flags, rho_phy,   &
  47                              chem, emis_ant,eghg_bio,ebio_co2oce    )
  48
  49 ! March-28, this subroutine adds all type of greenhouse gases to the chem species
  50
  51 !IMPLICIT NONE
  52
  53 TYPE(grid_config_rec_type), INTENT(IN   )    :: config_flags
  54
  55 INTEGER, INTENT(IN   ) ::   ids,ide, jds,jde, kds,kde,        &
  56                             ims,ime, jms,jme, kms,kme,        &
  57                             its,ite, jts,jte, kts,kte
  58
  59 REAL, INTENT(IN   )    ::   dtstep
  60
  61 REAL, DIMENSION( ims:ime,jms:jme ), INTENT(IN ) :: ebio_co2oce
  62 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(INOUT ) ::      chem
  63 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::   rho_phy, dz8w
  64 REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme, num_emis_ant ), INTENT(IN ) :: emis_ant
  65 REAL, DIMENSION( ims:ime, 1,jms:jme, num_eghg_bio ), INTENT(IN ) :: eghg_bio
  66
  67 INTEGER :: i,j,k
  68 REAL    :: conv_rho
  69
  70 call wrf_debug(15,'add_ghg_fluxes')
  71
  72 ! For both GHG options
  73 DO j=jts,jte
  74    DO i=its,ite
  75       !  3D anthropogenic fluxes
  76       DO k=kts,min(config_flags%kemit,kte)
  77          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]
  78          chem(i,k,j,p_co2_ant)= chem(i,k,j,p_co2_ant) + conv_rho* emis_ant(i,k,j,p_e_co2)
  79          chem(i,k,j,p_co2_tst)= chem(i,k,j,p_co2_tst) + conv_rho* emis_ant(i,k,j,p_e_co2tst)
  80          chem(i,k,j,p_co_ant) = chem(i,k,j,p_co_ant)  + conv_rho* emis_ant(i,k,j,p_e_co)
  81
  82       ! 2D biospheric fluxes:
  83       if (k==1) then
  84          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
  85          chem(i,1,j,p_co2_oce)= chem(i,1,j,p_co2_oce) + conv_rho* ebio_co2oce(i,j)
  86       end if
  87
  88       ENDDO
  89    ENDDO
  90 ENDDO
  91
  92 ! For the GHG_TRACER option only
  93 IF(config_flags%chem_opt==GHG_TRACER) THEN
  94     DO j=jts,jte
  95        DO i=its,ite
  96           !  3D anthropogenic fluxes
  97           DO k=kts,min(config_flags%kemit,kte)
  98              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]
  99              chem(i,k,j,p_ch4_ant)= chem(i,k,j,p_ch4_ant) + conv_rho* emis_ant(i,k,j,p_e_ch4)
 100              chem(i,k,j,p_ch4_tst)= chem(i,k,j,p_ch4_tst) + conv_rho* emis_ant(i,k,j,p_e_ch4tst)
 101              chem(i,k,j,p_co_tst) = chem(i,k,j,p_co_tst)  + conv_rho* emis_ant(i,k,j,p_e_cotst)
 102
 103           ! 2D biospheric fluxes:
 104           if (k==1) then
 105              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) &
 106                                                                        + eghg_bio(i,1,j,p_ebio_ch4term))
 107           end if
 108
 109           ENDDO
 110        ENDDO
 111     ENDDO
 112
 113 END IF
 114
 115 END SUBROUTINE add_ghg_fluxes
 116 !**************************************************************************************************
 117 SUBROUTINE VPRM                 (  ids,ide, jds,jde,                 &
 118                                    ims,ime, jms,jme,                 &
 119                                    its,ite, jts,jte,                 &
 120
 121                                    vprm_in,rad0, lambda,            &
 122                                    alpha, RESP0,                    &
 123                                    T2,RAD, eghg_bio                  )
 124 !IMPLICIT NONE
 125
 126 INTEGER,  INTENT(IN   )   ::  ids,ide, jds,jde, &
 127                               ims,ime, jms,jme, &
 128                               its,ite, jts,jte
 129 INTEGER :: i,j,m
 130
 131 REAL, PARAMETER :: const= 3.6e+3 ! For unit conversion from mmol/m^2/s to mol/km2/hr
 132
 133 REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN)    ::  T2,  RAD
 134
 135 REAL, DIMENSION( ims:ime, 8, jms:jme, num_vprm_in), INTENT(IN)   ::  vprm_in
 136 REAL, DIMENSION( ims:ime, 1,jms:jme, num_eghg_bio ), INTENT(INOUT ) ::  eghg_bio
 137
 138 REAL, DIMENSION(num_vprm_in) :: rad0, lambda, alpha, RESP0
 139 REAL  ::  a1,a2,a3,Tair,Tscale,Wscale,Pscale,GEE_frac,RESP_frac,evithresh,RADscale
 140 ! VPRM vegetation classes:
 141 !1-Evergreen c
 142 !2-Deciduous
 143 !3-Mixed forest
 144 !4-Shrubland
 145 !5-Savanna
 146 !6-Cropland
 147 !7-Grassland
 148 !8-Others
 149
 150 REAL, DIMENSION(8) ::  Tmax,Tmin,Topt
 151 ! These are universal VPRM parameters
 152 DATA Tmin   /0.,0.,0.,2.,2.,5.,2.,0./
 153 DATA Topt   /20.,20.,20.,20.,20.,22.,18.,0./
 154 DATA Tmax   /8*40./
 155
 156 DO j=jts,min(jte,jde-1)
 157 DO i=its,min(ite,ide-1)
 158
 159    GEE_frac= 0.
 160    RESP_frac= 0.
 161
 162    Tair= T2(i,j)-273.15
 163    veg_frac_loop: DO m=1,7
 164
 165         if (vprm_in(i,m,j,p_vegfra_vprm)<1.e-8) CYCLE  ! Then fluxes are zero
 166
 167         a1= Tair-Tmin(m)
 168         a2= Tair-Tmax(m)
 169         a3= Tair-Topt(m)
 170
 171         ! Here a1 or a2 can't be negative
 172         if (a1<0. .OR. a2>0.) then
 173             Tscale= 0.
 174         else
 175             Tscale=a1*a2/(a1*a2 - a3**2)
 176         end if
 177
 178         if (Tscale<0.) then
 179             Tscale=0.
 180         end if
 181
 182        ! modification due to different dependency on ground water
 183         if (m==4 .OR. m==7) then  ! grassland and shrubland are xeric systems
 184             if (vprm_in(i,m,j,p_lswi_max)<1e-7) then  ! in order to avoid NaN for Wscale
 185                 Wscale= 0.
 186             else
 187                 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))
 188             end if
 189         else
 190             Wscale= (1.+vprm_in(i,m,j,p_lswi))/(1.+vprm_in(i,m,j,p_lswi_max))
 191         end if
 192
 193         ! effect of leaf phenology
 194         if (m==1) then  ! evegreen
 195             Pscale= 1.
 196         else if (m==5 .OR. m==7) then  ! savanna or grassland
 197             Pscale= (1.+vprm_in(i,m,j,p_lswi))/2.
 198         else                           ! Other vegetation types
 199             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))
 200             if (vprm_in(i,m,j,p_evi)>=evithresh) then  ! Full canopy period
 201                Pscale= 1.
 202             else
 203                Pscale=(1.+vprm_in(i,m,j,p_lswi))/2.  ! bad-burst to full canopy period
 204             end if
 205         end if
 206
 207         RADscale= 1./(1. + RAD(i,j)/rad0(m))
 208         GEE_frac= lambda(m)*Tscale*Pscale*Wscale*RADscale* vprm_in(i,m,j,p_evi)* RAD(i,j)*vprm_in(i,m,j,p_vegfra_vprm) + GEE_frac
 209
 210         RESP_frac= (alpha(m)*Tair + RESP0(m))*vprm_in(i,m,j,p_vegfra_vprm) + RESP_frac
 211
 212
 213    ENDDO veg_frac_loop
 214
 215    eghg_bio(i,1,j,p_ebio_gee)= min(0.0,const*GEE_frac)
 216    eghg_bio(i,1,j,p_ebio_res)= max(0.0,const*RESP_frac)
 217
 218   ENDDO
 219   ENDDO
 220
 221 END SUBROUTINE VPRM
 222 !****************************************************************************************
 223
 224 !VB: here comes the subroutine for the kaplan wetland inventory (Kaplan, 2002+2006)
 225 !which calculates CH4 wetland emissions from wrf soil temperature and soil moisture fields
 226 !and a global carbon density and a potential wetland map (see also Sitch et al. 2003) for details
 227
 228 SUBROUTINE KAPLAN             ( ids,ide, jds,jde,                        &
 229                                 ims,ime, jms,jme,                        &
 230                                 its,ite, jts,jte,                        &
 231
 232                                 dt, T_soil, Soil_M, wet_in,              &
 233                                 soil_type, T_skin, eghg_bio,             &
 234                                 num_soil_layers,E_f,M_s                  )
 235
 236 !IMPLICIT NONE
 237
 238 INTEGER,  INTENT(IN   )   ::  ids,ide, jds,jde, &
 239                               ims,ime, jms,jme, &
 240                               its,ite, jts,jte, &
 241                               num_soil_layers
 242 INTEGER :: i,j
 243 INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN)  :: soil_type
 244
 245 REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr
 246 REAL, DIMENSION (ims:ime,1,jms:jme, num_wet_in), INTENT(IN)    :: wet_in
 247 REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN)    :: T_skin
 248 REAL, DIMENSION (ims:ime, num_soil_layers, jms:jme), INTENT(IN) :: Soil_M, T_soil
 249 REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) ::  eghg_bio
 250 REAL, INTENT(IN) :: dt,M_s,E_f
 251
 252 REAL, DIMENSION( ims:ime, jms:jme) :: sm_01, sm_02, sm_03, sm_04, g_T, f_SM, &
 253                                       k_r, Carbon, HR, T_flood, sm_tot, T_2s, &
 254                                       T_peat, P_l, soil_sat
 255
 256 REAL  :: tau_10
 257
 258 ! Description of the variables used in the subroutine KAPLAN
 259 ! T_soil: soil temperature from WRF
 260 ! C_pool: carbon content of the fast carbon pool at time 0 - from the LPJ model 2000
 261 ! Wet_map: potential wetland map of Jed Kaplan gives the percentage of wetland area
 262 !          covered by each grid cell
 263 ! ebio_ch4wet: methane wetland emissions - output from kaplan model
 264 ! Soil_M: 4(3) dimensional wrf soil moisture output
 265 ! sm_01: first wrf soil moisture layer
 266 ! sm_02: second wrf soil moisture layer
 267 ! sm_03: third wrf soil moisture layer
 268 ! sm_04: fourth wrf soil moisture layer
 269 ! g_T: gives a modified Arrhenius temperature dependence
 270 ! f_SM: factor accounting for the soil moisture
 271 ! k_r: decomposition rate
 272 ! C_x: carbon pool after a certain time
 273 ! HR: heteorotrophic respiration
 274 ! T_flood: methane emissions for the flooding part - not peat land
 275 ! sm_tot: the actual soil moisture used for the wetland emissions is built up out of the wrf soil moisture
 276 ! T_2s: temperature of first soil layer from surface used for calculation of wetland emissions
 277 ! E_f: correction factor accounting for the peatland wetland emissions
 278 ! T_peat: Ch4 wetland emissions from peat ground
 279 ! P_l: weighting factor for the global peatland and flooded wetland emissions
 280 ! T_a: mean annual surface temperature used for the weighting factor
 281 ! but this mean annual temperature has to come from another source (e.g. GEOS)
 282 ! dt: timestep of the wrf grid in seconds
 283 ! soil_type: soil type for getting saturated soil moisture
 284 ! soil_sat: value of saturated soil moisture
 285 ! T_skin: skin temperature to use instead of soil temperature
 286
 287 ! Initialization of the constant variables
 288
 289   tau_10 = 2.86  !litter soil pool
 290
 291  DO j=jts,min(jte,jde-1)
 292   DO i=its,min(ite,ide-1)
 293
 294 ! zero initialization of the fields to avoid crashing of the program
 295
 296     eghg_bio(i,1,j,p_ebio_ch4wet)=0.
 297     HR(i,j)=0.
 298     k_r(i,j)=0.
 299     g_T(i,j)=0.
 300     T_flood(i,j)=0.
 301     T_peat(i,j)=0.
 302     P_l(i,j)=0.
 303     f_SM(i,j)=0.
 304     sm_01(i,j)=0.
 305     sm_02(i,j)=0.
 306     sm_03(i,j)=0.
 307     sm_04(i,j)=0.
 308     Carbon(i,j)=0.
 309     sm_tot(i,j)=0.
 310     T_2s(i,j)=0.
 311
 312 ! If loop to derive the saturated soil moisture for the specific soil type
 313
 314   IF(soil_type(i,j) == 1)then
 315    soil_sat(i,j) = 0.339
 316   ELSE IF(soil_type(i,j) == 2) then
 317    soil_sat(i,j) = 0.421
 318   ELSE IF(soil_type(i,j) == 3) then
 319    soil_sat(i,j) = 0.434
 320   ELSE IF(soil_type(i,j) == 4) then
 321    soil_sat(i,j) = 0.476
 322   ELSE IF(soil_type(i,j) == 5) then
 323    soil_sat(i,j) = 0.476
 324   ELSE IF(soil_type(i,j) == 6) then
 325    soil_sat(i,j) = 0.439
 326   ELSE IF(soil_type(i,j) == 7) then
 327    soil_sat(i,j) = 0.404
 328   ELSE IF(soil_type(i,j) == 8) then
 329    soil_sat(i,j) = 0.464
 330   ELSE IF(soil_type(i,j) == 9) then
 331    soil_sat(i,j) = 0.465
 332   ELSE IF(soil_type(i,j) == 10) then
 333    soil_sat(i,j) = 0.406
 334   ELSE IF(soil_type(i,j) == 11) then
 335    soil_sat(i,j) = 0.468
 336   ELSE IF(soil_type(i,j) == 12) then
 337    soil_sat(i,j) = 0.468
 338   ELSE IF(soil_type(i,j) == 13) then
 339    soil_sat(i,j) = 0.439
 340   ELSE IF(soil_type(i,j) == 14) then
 341    soil_sat(i,j) = 1.0
 342   ELSE IF(soil_type(i,j) == 15) then
 343    soil_sat(i,j) = 0.20
 344   ELSE IF(soil_type(i,j) == 16) then
 345    soil_sat(i,j) = 0.421
 346   ELSE IF(soil_type(i,j) == 17) then
 347    soil_sat(i,j) = 0.468
 348   ELSE IF(soil_type(i,j) == 18) then
 349    soil_sat(i,j) = 0.200
 350   ELSE IF(soil_type(i,j) == 19) then
 351    soil_sat(i,j) = 0.339
 352   ELSE
 353    soil_sat(i, j) = 0.4
 354   END IF
 355
 356 ! Use mean value of first and 2nd soil layer (up to 30cm depth)
 357 ! and the surface soil temperature
 358
 359   sm_01(i,j)  = Soil_M(i,1,j)
 360   sm_02(i,j)  = Soil_M(i,2,j)
 361   sm_03(i,j)  = Soil_M(i,3,j)
 362   sm_04(i,j)  = Soil_M(i,4,j)
 363
 364   sm_tot(i,j) = 0.5*(sm_01(i,j)+sm_02(i,j))
 365
 366   IF(sm_tot(i,j) < 0.0) then
 367      sm_tot(i,j) = 0.0
 368   END IF
 369
 370   IF ((T_Soil(i,1,j) == 273.16).OR.(T_Soil(i,1,j) == 0.0)) then
 371      T_2s(i, j) = T_skin(i, j) - 273.15
 372   ELSE
 373      T_2s(i,j)= T_Soil(i,1,j)-273.15
 374   END IF
 375
 376 ! calculate temperature dependence all based on sitch (2003)
 377
 378     g_T(i,j) = EXP(308.56*(1.0/56.02 - 1.0/(T_2s(i,j) + 46.02)))
 379
 380 ! calculate soil moisture factor
 381 ! soil moisture should be sm/saturated sm (Folley)
 382
 383     f_SM(i,j) = 0.25 + 0.75 * (sm_tot(i,j)/soil_sat(i,j))
 384
 385 ! calculate decomposition rate per year/per hr 20.04.10
 386
 387     k_r(i,j) = ((1/tau_10)*g_T(i,j)*f_SM(i,j))/(12.0*24.0*30.0)
 388
 389 ! calculate the variable carbon fast pool depending on a month
 390
 391     Carbon(i,j) = wet_in(i,1,j,p_Cpool)*(1.0 - EXP(-k_r(i,j)*1.0))
 392
 393 ! linearize and get the decomposed carbon for one second
 394
 395     Carbon(i,j) = Carbon(i,j)/(3600.0)
 396
 397 ! calculate the heteorotrophic respiration
 398
 399     HR(i,j) =  0.7*Carbon(i,j)
 400
 401 ! calculate the wetland emissions per grid cell for the flooding part
 402
 403   T_flood(i,j) = HR(i,j)*M_s*wet_in(i,1,j,p_wetmap)
 404
 405 ! calculate wetland emissions for the peatland part (higher latitudes)
 406
 407   T_peat(i,j) = HR(i,j)*E_f*wet_in(i,1,j,p_wetmap)
 408
 409 ! the factor for the relation between flood and peat land
 410
 411    P_l(i,j) = EXP((wet_in(i,1,j,p_T_ann) - 303.0)/8.0)
 412    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)
 413
 414   ! units
 415   eghg_bio(i,1,j,p_ebio_ch4wet) = eghg_bio(i,1,j,p_ebio_ch4wet)*const
 416
 417   END DO
 418
 419   END DO
 420
 421 END SUBROUTINE KAPLAN
 422 !***************************************************************************************
 423   SUBROUTINE TERMITE (                    ids,ide, jds,jde,                            &
 424                                           ims,ime, jms,jme,                            &
 425                                           its,ite, jts,jte,                            &
 426
 427                                           dt,eghg_bio, vtype_ter,                      &
 428                                           biomt_par, emit_par                          )
 429
 430 !IMPLICIT NONE
 431
 432 INTEGER,  INTENT(IN   )   ::  ids,ide, jds,jde, &
 433                               ims,ime, jms,jme, &
 434                               its,ite, jts,jte
 435 INTEGER :: i,j,vtype
 436
 437 INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN)    ::  vtype_ter
 438 REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) ::  eghg_bio
 439 REAL, INTENT(IN) :: dt
 440
 441 REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr
 442 REAL, DIMENSION(14), INTENT(IN) :: biomt_par
 443 REAL, DIMENSION(14), INTENT(IN) :: emit_par
 444
 445 CHARACTER*256 :: mminlu_loc
 446
 447  DO j=jts,min(jte,jde-1)
 448   DO i=its,min(ite,ide-1)
 449
 450    CALL nl_get_mminlu(1 , mminlu_loc)
 451
 452     IF(mminlu_loc .EQ. 'USGS')THEN
 453
 454      vtype = INT(vtype_ter(i,j))
 455
 456      SELECT CASE (vtype)
 457
 458        CASE (1)
 459           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 460        CASE (2)
 461           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10)
 462        CASE (3)
 463           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10)
 464        CASE (4)
 465           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11)
 466        CASE (5)
 467           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11)
 468        CASE (6)
 469           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11)
 470        CASE (7)
 471           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(7)*emit_par(7)
 472        CASE (8)
 473           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12)
 474        CASE (9)
 475           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12)
 476        CASE (10)
 477           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(5)*emit_par(5)
 478        CASE (11)
 479           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 480        CASE (12)
 481           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 482        CASE (13)
 483           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(1)*emit_par(1)
 484        CASE (14)
 485           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(2)*emit_par(2)
 486        CASE (15)
 487           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 488        CASE (16)
 489           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 490        CASE (17)
 491           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 492        CASE (18)
 493           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 494        CASE (19)
 495           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(13)*emit_par(13)
 496        CASE (20)
 497           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 498        CASE (21)
 499           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 500        CASE (22)
 501           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 502        CASE (23)
 503           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 504        CASE (24)
 505           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 506        CASE default
 507           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 508        END SELECT
 509
 510
 511     ELSE IF(mminlu_loc .EQ. 'MODIFIED_IGBP_MODIS_NOAH')THEN
 512
 513      vtype = INT(vtype_ter(i,j))
 514
 515      SELECT CASE (vtype)
 516
 517        CASE (1)
 518           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(2)*emit_par(2)
 519        CASE (2)
 520           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(1)*emit_par(1)
 521        CASE (3)
 522           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 523        CASE (4)
 524           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 525        CASE (5)
 526           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(4)*emit_par(4)
 527        CASE (6)
 528           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12)
 529        CASE (7)
 530           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(12)*emit_par(12)
 531        CASE (8)
 532           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(6)*emit_par(6)
 533        CASE (9)
 534           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(5)*emit_par(5)
 535        CASE (10)
 536           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(7)*emit_par(7)
 537        CASE (11)
 538           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 539        CASE (12)
 540           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(10)*emit_par(10)
 541        CASE (13)
 542           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 543        CASE (14)
 544           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(11)*emit_par(11)
 545        CASE (15)
 546           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 547        CASE (16)
 548           eghg_bio(i,1,j,p_ebio_ch4term) = (1.0/3.6)*1.0e-12*biomt_par(13)*emit_par(13)
 549        CASE (17)
 550           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 551        CASE (18)
 552           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 553        CASE (19)
 554           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 555        CASE (20)
 556           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 557        CASE default
 558           eghg_bio(i,1,j,p_ebio_ch4term) = 0.0
 559        END SELECT
 560
 561     END IF
 562
 563     eghg_bio(i,1,j,p_ebio_ch4term) = eghg_bio(i,1,j,p_ebio_ch4term)*const
 564
 565   END DO
 566  END DO
 567
 568
 569 END SUBROUTINE TERMITE
 570
 571 !************************************************************
 572 !subroutine for soil uptake following ridgwell et al. 1999
 573 !************************************************************
 574
 575   SUBROUTINE SOILUPTAKE            ( ids,ide, jds,jde,                      &
 576                                      ims,ime, jms,jme,                      &
 577                                      its,ite, jts,jte,                      &
 578
 579                                      Soil_M, soil_typ, eghg_bio,            &
 580                                      rain_1, rain_2,                        &
 581                                      potevap, sfevap, landuse, T2, dt,      &
 582                                      num_soil_layers, wet_in                )
 583
 584   !IMPLICIT NONE
 585
 586   INTEGER,  INTENT(IN   )   ::  ids,ide, jds,jde, &
 587                                 ims,ime, jms,jme, &
 588                                 its,ite, jts,jte, &
 589                                 num_soil_layers
 590   INTEGER :: i,j
 591
 592   INTEGER, DIMENSION (ims:ime, jms:jme), INTENT(IN)    :: soil_typ
 593
 594
 595   REAL, PARAMETER :: const= 3.6e9/(0.012+4.*0.001) ! For unit conversion from kgCH4/m^2/s to molCH4/km2/hr
 596   REAL, DIMENSION (ims:ime, jms:jme), INTENT(IN)  ::  rain_1, T2,       &
 597                                                       rain_2, landuse,  &
 598                                                       potevap, sfevap
 599
 600   REAL, DIMENSION (ims:ime, num_soil_layers, jms:jme), INTENT(IN) :: Soil_M
 601   REAL, INTENT(IN) :: dt
 602
 603   REAL, DIMENSION (ims:ime, 1,jms:jme, num_eghg_bio), INTENT(INOUT) ::  eghg_bio
 604   REAL, DIMENSION (ims:ime, 1,jms:jme, num_wet_in), INTENT(IN) ::  wet_in
 605
 606   REAL, DIMENSION(ims:ime, jms:jme) :: sm_01, sm_02, sm_03, sm_04,    &
 607                                        i_cult, r_n, r_sm, vero, k_d, &
 608                                        G_soil, G_t, sm_res, dch4,     &
 609                                        phi_soil, eps, b_f, sand_c,    &
 610                                        clay_c, vero1
 611
 612   REAL :: dch4_0, k_0, z_dsoil, conv_Fk, conv_ch4, methane_0
 613
 614   CHARACTER*256 :: mminlu_loc
 615
 616 !Initialization of the constants
 617
 618    dch4_0    =       0.196         !cm2/s  diffusivity of CH4 in free air
 619    k_0       =       0.00087       !1/s   base oxidation rate constant for uncultivated
 620                                    !moist soil at 0 C obtained from experiments
 621    z_dsoil   =       6.0           ! soil depth in cm at which the oxidation activity is assumed
 622    conv_Fk   =       616.9         !mg/ppmv cm CH4  conversion factor fpr giving soil_flux right units of kgCH4/m^2s
 623    conv_ch4  =       (28.97/26.0)*1e6  !ppm to kgCH4/kgair
 624 !   methane_0 =       9.66*1e-7     ! 1750ppb CH4 global average concentration
 625    methane_0 =       1.77     ! 1770ppb CH4 global average concentration
 626
 627  DO j=jts,min(jte,jde-1)
 628   DO i=its,min(ite,ide-1)
 629
 630 ! initialization of all variables
 631
 632   sm_01(i, j)  = 0.0
 633   sm_02(i, j)  = 0.0
 634   sm_03(i, j)  = 0.0
 635   sm_04(i, j)  = 0.0
 636   i_cult(i, j) = 0.0
 637   r_n(i, j)    = 0.0
 638   r_sm(i, j)   = 0.0
 639   vero(i, j)   = 0.0
 640   k_d(i, j)    = 0.0
 641   G_soil(i, j) = 0.0
 642   G_t(i, j)    = 0.0
 643   sm_res(i, j) = 0.0
 644   dch4(i, j)   = 0.0
 645   phi_soil(i, j) = 0.0
 646   eps(i, j)    = 0.0
 647   b_f(i, j)    = 0.0
 648   sand_c(i, j) = 0.0
 649   clay_c(i, j) = 0.0
 650   vero1(i, j)  = 0.0
 651
 652 ! end initialization
 653
 654   sm_01(i,j)  = Soil_M(i,1,j)
 655   sm_02(i,j)  = Soil_M(i,2,j)
 656   sm_03(i,j)  = Soil_M(i,3,j)
 657   sm_04(i,j)  = Soil_M(i,4,j)
 658
 659 !----------------------------------------------------------------------
 660 ! calculation of CH4 diffusion activity
 661 !----------------------------------------------------------------------
 662
 663      G_t(i, j) = 1.0 + 0.0055 * (T2(i, j)-273.15)
 664
 665 !************************
 666 ! here the values for air filled porosity and total por volume have to be filled in
 667 !************************
 668
 669 ! calculation of the soil parameters air filled porosity, total pore volume and b factor
 670 ! depending on the soil type
 671 ! and also clay/sand content - from WRF SOILPARAM.TBL and triangle plot
 672 ! sand and clay content and porosity from Cosby et al. 1984
 673
 674   IF (soil_typ(i, j) == 1) then   !sand
 675          phi_soil(i, j) = 0.339
 676          eps(i, j)      = 0.103
 677          sand_c(i,j)    = 0.92
 678          clay_c(i,j)    = 0.03
 679          b_f(i, j)      = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 680                           *sand_c(i, j)**2 *clay_c(i, j)
 681   ELSE IF (soil_typ(i, j) == 2) then  !loamy sand
 682          phi_soil(i, j) = 0.421
 683          eps(i, j)      = 0.038
 684          sand_c(i,j)    = 0.82
 685          clay_c(i,j)    = 0.06
 686          b_f(i, j)      = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 687                            *sand_c(i, j)**2 *clay_c(i, j)
 688   ELSE IF (soil_typ(i, j) == 3) then !sandy loam
 689          phi_soil(i, j) = 0.434
 690          eps(i, j) = 0.051
 691          sand_c(i,j) = 0.58
 692          clay_c(i,j) = 0.10
 693          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 694                      *sand_c(i, j)**2 *clay_c(i, j)
 695   ELSE IF (soil_typ(i, j) == 4) then ! silt loam
 696          phi_soil(i, j) = 0.476
 697          eps(i, j) = 0.116
 698          sand_c(i,j) = 0.17
 699          clay_c(i,j) = 0.13
 700          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 701                      *sand_c(i, j)**2 *clay_c(i, j)
 702   ELSE IF (soil_typ(i, j) == 5) then !silt
 703          phi_soil(i, j) = 0.476
 704          eps(i, j) = 0.093
 705          sand_c(i,j) = 0.10
 706          clay_c(i,j) = 0.10  !estimated no values
 707          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 708                      *sand_c(i, j)**2 *clay_c(i, j)
 709   ELSE IF (soil_typ(i, j) == 6) then  !loam
 710          phi_soil(i, j) = 0.439
 711          eps(i, j) = 0.11
 712          sand_c(i,j) = 0.43
 713          clay_c(i,j) = 0.18
 714          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 715                      *sand_c(i, j)**2 *clay_c(i, j)
 716   ELSE IF (soil_typ(i, j) == 7) then  !sandy clay loam
 717          phi_soil(i, j) = 0.464
 718          eps(i, j) = 0.09
 719          sand_c(i,j) = 0.58
 720          clay_c(i,j) = 0.27
 721          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 722                      *sand_c(i, j)**2 *clay_c(i, j)
 723   ELSE IF (soil_typ(i, j) == 8) then  !silty clay loam
 724          phi_soil(i, j) = 0.404
 725          eps(i, j) = 0.077
 726          sand_c(i,j) = 0.10
 727          clay_c(i,j) = 0.34
 728          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 729                      *sand_c(i, j)**2 *clay_c(i, j)
 730   ELSE IF (soil_typ(i, j) == 9) then  !clay loam
 731          phi_soil(i, j) = 0.465
 732          eps(i, j) = 0.083
 733          sand_c(i,j) = 0.32
 734          clay_c(i,j) = 0.34
 735          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 736                      *sand_c(i, j)**2 *clay_c(i, j)
 737   ELSE IF (soil_typ(i, j) == 10) then  !sandy clay
 738          phi_soil(i, j) = 0.406
 739          eps(i, j) = 0.068
 740          sand_c(i,j) = 0.52
 741          clay_c(i,j) = 0.42
 742          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 743                      *sand_c(i, j)**2 *clay_c(i, j)
 744   ELSE IF (soil_typ(i, j) == 11) then !silty clay
 745          phi_soil(i, j) = 0.468
 746          eps(i, j) = 0.064
 747          sand_c(i,j) = 0.06
 748          clay_c(i,j) = 0.47
 749          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 750                      *sand_c(i, j)**2 *clay_c(i, j)
 751   ELSE IF (soil_typ(i, j) == 12) then  !clay
 752          phi_soil(i, j) = 0.468
 753          eps(i, j) = 0.056
 754          sand_c(i,j) = 0.22
 755          clay_c(i,j) = 0.58
 756          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 757                      *sand_c(i, j)**2 *clay_c(i, j)
 758   ELSE IF (soil_typ(i, j) == 13) then !organic material
 759          phi_soil(i, j) = 0.439
 760          eps(i, j) = 0.11
 761          sand_c(i,j) = 0.05
 762          clay_c(i,j) = 0.05 !estimated
 763          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 764                      *sand_c(i, j)**2 *clay_c(i, j)
 765   ELSE IF (soil_typ(i, j) == 14) then  !water
 766          phi_soil(i, j) = 1.0
 767          eps(i, j) = 0.0
 768          sand_c(i,j) = 0.60
 769          clay_c(i,j) = 0.40  !estimated
 770          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 771                      *sand_c(i, j)**2 *clay_c(i, j)
 772   ELSE IF (soil_typ(i, j) == 15) then  !bedrock
 773          phi_soil(i, j) = 0.2
 774          eps(i, j) = 0.03
 775          sand_c(i,j) = 0.07
 776          clay_c(i,j) = 0.06
 777          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 778                      *sand_c(i, j)**2 *clay_c(i, j)
 779   ELSE IF (soil_typ(i, j) == 16) then  !other
 780          phi_soil(i, j) = 0.421
 781          eps(i, j) = 0.138
 782          sand_c(i,j) = 0.25
 783          clay_c(i,j) = 0.25
 784          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 785                      *sand_c(i, j)**2 *clay_c(i, j)
 786   ELSE IF (soil_typ(i, j) == 17) then  ! playa
 787          phi_soil(i, j) = 0.468
 788          eps(i, j) = 0.014
 789          sand_c(i,j) = 0.60
 790          clay_c(i,j) = 0.40
 791          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 792                      *sand_c(i, j)**2 *clay_c(i, j)
 793  ELSE IF (soil_typ(i, j) == 18) then  !lava
 794          phi_soil(i, j) = 0.468
 795          eps(i, j) = 0.03
 796          sand_c(i,j) = 0.52
 797          clay_c(i,j) = 0.48
 798          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 799                      *sand_c(i, j)**2 *clay_c(i, j)
 800   ELSE IF (soil_typ(i, j) == 19) then  !white sand
 801          phi_soil(i, j) = 0.339
 802          eps(i, j) = 0.103
 803          sand_c(i,j) = 0.92
 804          clay_c(i,j) = 0.03
 805          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 806                      *sand_c(i, j)**2 *clay_c(i, j)
 807   ELSE
 808          phi_soil(i, j) = 0.339
 809          eps(i, j) = 0.0
 810          sand_c(i,j) = 0.92
 811          clay_c(i,j) = 0.03
 812          b_f(i, j) = -3.140 - 0.00222* clay_c(i, j)**2 - 3.484*10e-5  &
 813                      *sand_c(i, j)**2 *clay_c(i, j)
 814    END IF
 815
 816    G_soil(i, j) = (phi_soil(i, j)**(4.0/3.0))* &
 817                   ((eps(i, j)/phi_soil(i, j))**(1.5 + 3.0/b_f(i, j)))
 818
 819 ! Calculate the Soil CH4 Diffusivity
 820
 821    dch4(i, j) = G_soil(i, j)*G_t(i, j)*dch4_0
 822
 823 !-------------------------------------------------------------!
 824 !  Calculation of CH4 Oxidation Activity                      !
 825 !-------------------------------------------------------------!
 826
 827 ! Calculation of the temperature response to ch4 oxidation rate
 828
 829
 830    IF(T2(i,j)<= 0.0) then
 831
 832    vero(i,j) = 0.0
 833
 834    ELSE IF(T2(i,j) > 0.0) then
 835
 836    vero1(i, j) = (T2(i, j)-273.15)**4.0
 837
 838    vero(i,j) = EXP((0.0693*(T2(i,j)-273.15))-(8.56*10e-7*vero1(i, j)))
 839
 840
 841    END IF
 842
 843
 844 ! Calculation of the cultivation response
 845  CALL nl_get_mminlu(1, mminlu_loc)
 846
 847     IF(mminlu_loc .EQ. 'USGS') THEN
 848
 849        IF(landuse(i, j) <= 5.0) THEN
 850
 851             i_cult(i, j) = 1.0
 852
 853        ELSE IF(landuse(i, j) == 6.0) THEN
 854
 855             i_cult(i, j) = 0.5
 856
 857        ELSE
 858
 859             i_cult(i, j) = 0.0
 860
 861        END IF
 862
 863    ELSE IF(mminlu_loc .EQ. 'MODIFIED_IGBP_MODIS_NOAH')THEN
 864
 865        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
 866
 867              i_cult(i, j) = 1.0
 868
 869        ELSE IF(landuse(i, j) == 14.0 .OR. landuse(i,j) == 31.0) THEN
 870
 871              i_cult(i, j) = 0.5
 872
 873        ELSE
 874
 875              i_cult(i, j) = 0.0
 876
 877        END IF
 878
 879   END IF
 880
 881   r_n(i, j) = 1.0 - (0.75*i_cult(i, j))
 882
 883 ! Calculation of soil moisture response
 884 !precipitation has to be monthly (*30)
 885 ! soil moisture 2nd layer (30cm)
 886
 887
 888    IF(sfevap(i, j) == 0.0) THEN
 889
 890     sm_res(i, j) = 2.0
 891
 892   ELSE
 893
 894      sm_res(i, j) = ((rain_1(i, j) + rain_2(i, j))*30.0*24.*3600./dt*0.001 &
 895                      + sm_02(i, j))/sfevap(i, j)
 896
 897   END IF
 898
 899   IF(sm_res(i, j) > 1.0)THEN
 900
 901     r_sm(i, j) = 1.0
 902
 903   ELSE IF(sm_res(i, j) <= 1.0)THEN
 904
 905     r_sm(i, j) = sm_res(i, j)
 906
 907   END IF
 908
 909   k_d(i, j) = r_n(i, j)*r_sm(i, j)*vero(i, j)*k_0
 910
 911 !******************************************
 912 ! Total soil uptake CH4
 913 !******************************************
 914
 915 ! this gives you the flux in mg CH4/m^2 day
 916
 917   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
 918
 919   eghg_bio(i,1,j,p_ebio_ch4soil) = ((methane_0*dch4(i, j))/z_dsoil)* &
 920                                  (1.0 - dch4(i, j)/(dch4(i, j)+k_d(i, j)*z_dsoil))*conv_Fk
 921
 922 ! flux in mol/km2/hour
 923   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
 924
 925
 926   ELSE
 927
 928   eghg_bio(i,1,j,p_ebio_ch4soil) = 0.0
 929
 930   END IF
 931
 932   IF(eghg_bio(i,1,j,p_ebio_ch4soil) > 0.0) then
 933
 934    eghg_bio(i,1,j,p_ebio_ch4soil) = 0.0
 935
 936   END IF
 937
 938   END DO
 939   END DO
 940
 941   END SUBROUTINE SOILUPTAKE
 942
 943 END MODULE module_ghg_fluxes