diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index a6c2d8f771..3f0673227f 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -2913,8 +2913,6 @@ package ctsmscheme sf_surface_physics==6 - - package pxlsmscheme sf_surface_physics==7 - state:t2_ndg_new,q2_ndg_new,t2_ndg_old,q2_ndg_old,t2obs,q2obs,vegf_px,imperv,canfra,lai_px,wwlt_px,wfc_px,wsat_px,clay_px,csand_px,fmsand_px package ssibscheme sf_surface_physics==8 - state:ssib_fm,ssib_fh,ssib_cm,ssibxdd,ssib_br,ssib_lhf,ssib_shf,ssib_ghf,ssib_egs,ssib_eci,ssib_ect,ssib_egi,ssib_egt,ssib_sdn,ssib_sup,ssib_ldn,ssib_lup,ssib_wat,ssib_shc,ssib_shg,ssib_lai,ssib_vcf,ssib_z00,ssib_veg,isnow,swe,snowden,snowdepth,tkair,dzo1,wo1,tssn1,tssno1,bwo1,bto1,cto1,fio1,flo1,bio1,blo1,ho1,dzo2,wo2,tssn2,tssno2,bwo2,bto2,cto2,fio2,flo2,bio2,blo2,ho2,dzo3,wo3,tssn3,tssno3,bwo3,bto3,cto3,fio3,flo3,bio3,blo3,ho3,dzo4,wo4,tssn4,tssno4,bwo4,bto4,cto4,fio4,flo4,bio4,blo4,ho4 -package noahmpoptcrop2 opt_crop==2 - state:gecros_state - package noahmosaicscheme sf_surface_mosaic==1 - state:TSK_mosaic,QSFC_mosaic,TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic,CANWAT_mosaic,SNOW_mosaic,SNOWH_mosaic,SNOWC_mosaic,ALBEDO_mosaic,ALBBCK_mosaic,EMISS_mosaic,EMBCK_mosaic,ZNT_mosaic,Z0_mosaic,HFX_mosaic,QFX_mosaic,LH_mosaic,GRDFLX_mosaic,SNOTIME_mosaic,TR_URB2D_mosaic,TB_URB2D_mosaic,TG_URB2D_mosaic,TC_URB2D_mosaic,TS_URB2D_mosaic,TS_RUL2D_mosaic,QC_URB2D_mosaic,UC_URB2D_mosaic,TRL_URB3D_mosaic,TBL_URB3D_mosaic,TGL_URB3D_mosaic,SH_URB2D_mosaic,LH_URB2D_mosaic,G_URB2D_mosaic,RN_URB2D_mosaic,mosaic_cat_index,landusef2,smcrel,RS_mosaic,LAI_mosaic # IRRIGATION package channel sf_surf_irr_alloc==1 - state:irrigation,irr_rand_field diff --git a/Registry/Registry.NMM b/Registry/Registry.NMM index bd9849bc9c..0316f0854a 100644 --- a/Registry/Registry.NMM +++ b/Registry/Registry.NMM @@ -1211,8 +1211,6 @@ state real season_gdd ij - 1 - i02rhd=(interp_mask_land state integer cropcat ij - 1 - rh "cropcat" "dominant crop category" "category" state integer pgsxy ij - 1 - rh "pgs" "pgs" "" -state real gecros_state i{gecros}j - 1 Z rhd=(interp_mask_land_field:lu_index)u=(UpNear) "gecros" "gecros state array" "-" - # added state for etampnew microphysics (needed for restarts) state real mp_restart_state p misc 1 - r "MP_RESTART_STATE" "STATE VECTOR FOR MICROPHYSICS RESTARTS" state real tbpvs_state p misc 1 - r "TBPVS_STATE" "STATE FOR ETAMPNEW MICROPHYSICS" @@ -1846,7 +1844,7 @@ package slabscheme sf_surface_physics==1 - - package lsmscheme sf_surface_physics==2 - state:flx4,fvb,fbur,fgsn package ruclsmscheme sf_surface_physics==3 - state:ACRUNOFF -package noahmpscheme sf_surface_physics==4 - state:isnowxy,tvxy,tgxy,canliqxy,canicexy,eahxy,tahxy,cmxy,chxy,fwetxy,sneqvoxy,alboldxy,qsnowxy,wslakexy,zwtxy,waxy,wtxy,tsnoxy,zsnsoxy,snicexy,snliqxy,lfmassxy,rtmassxy,stmassxy,woodxy,stblcpxy,fastcpxy,xsaixy,taussxy,t2mvxy,t2mbxy,q2mvxy,q2mbxy,tradxy,neexy,gppxy,nppxy,fvegxy,qinxy,runsfxy,runsbxy,ecanxy,edirxy,etranxy,fsaxy,firaxy,aparxy,psnxy,savxy,sagxy,rssunxy,rsshaxy,bgapxy,wgapxy,tgvxy,tgbxy,chvxy,chbxy,shgxy,shcxy,shbxy,evgxy,evbxy,ghvxy,ghbxy,irgxy,ircxy,irbxy,trxy,evcxy,chleafxy,chucxy,chv2xy,chb2xy,chstarxy,smoiseq,smcwtdxy,rechxy,deeprechxy,fdepthxy,areaxy,rivercondxy,riverbedxy,eqzwt,pexpxy,qrfxy,qrfsxy,qspringxy,qspringsxy,qslatxy,stepwtd,gddxy,grainxy,croptype,planting,harvest,season_gdd,cropcat,pgsxy,soilcomp,soilcl1,soilcl2,soilcl3,soilcl4,gecros_state +package noahmpscheme sf_surface_physics==4 - state:isnowxy,tvxy,tgxy,canliqxy,canicexy,eahxy,tahxy,cmxy,chxy,fwetxy,sneqvoxy,alboldxy,qsnowxy,wslakexy,zwtxy,waxy,wtxy,tsnoxy,zsnsoxy,snicexy,snliqxy,lfmassxy,rtmassxy,stmassxy,woodxy,stblcpxy,fastcpxy,xsaixy,taussxy,t2mvxy,t2mbxy,q2mvxy,q2mbxy,tradxy,neexy,gppxy,nppxy,fvegxy,qinxy,runsfxy,runsbxy,ecanxy,edirxy,etranxy,fsaxy,firaxy,aparxy,psnxy,savxy,sagxy,rssunxy,rsshaxy,bgapxy,wgapxy,tgvxy,tgbxy,chvxy,chbxy,shgxy,shcxy,shbxy,evgxy,evbxy,ghvxy,ghbxy,irgxy,ircxy,irbxy,trxy,evcxy,chleafxy,chucxy,chv2xy,chb2xy,chstarxy,smoiseq,smcwtdxy,rechxy,deeprechxy,fdepthxy,areaxy,rivercondxy,riverbedxy,eqzwt,pexpxy,qrfxy,qrfsxy,qspringxy,qspringsxy,qslatxy,stepwtd,gddxy,grainxy,croptype,planting,harvest,season_gdd,cropcat,pgsxy,soilcomp,soilcl1,soilcl2,soilcl3,soilcl4 package clmscheme sf_surface_physics==5 - - package gfdlslab sf_surface_physics==88 - - diff --git a/Registry/registry.dimspec b/Registry/registry.dimspec index a680a6fce7..71fba341a5 100644 --- a/Registry/registry.dimspec +++ b/Registry/registry.dimspec @@ -93,7 +93,6 @@ dimspec lake_intl 2 constant=16 z interface_level dimspec lake_sll 2 constant=10 z soil_levels_or_lake_levels dimspec crop 2 constant=5 z crop_types dimspec soilc 2 constant=8 z soil_composition_layers -dimspec gecros 2 constant=60 z gecros # Set up the dust erosion dimension irrespective of WRF-Chem, but namelist erosion_dim=3 mandatory when not using WRF-Chem (gthompsn, 2017Jun23) dimspec . 3 namelist=erosion_dim z dust_erosion_dimension diff --git a/Registry/registry.noahmp b/Registry/registry.noahmp index daf0511f41..71c00dae62 100644 --- a/Registry/registry.noahmp +++ b/Registry/registry.noahmp @@ -108,5 +108,3 @@ state real season_gdd ij - 1 - i012rd=(interp_mask_fiel state integer cropcat ij - 1 - rh "cropcat" "dominant crop category" "category" state integer pgsxy ij - 1 - rh "pgs" "pgs" "" -# Gecros crop model only -state real gecros_state i{gecros}j - 1 Z rhd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "gecros" "gecros state array" "-" diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 51d77ec230..4e9c5a235d 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -821,7 +821,6 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,riverbedxy=grid%riverbedxy, eqzwt=grid%eqzwt, pexpxy=grid%pexpxy & & ,qrfxy=grid%qrfxy, qspringxy=grid%qspringxy, qslatxy=grid%qslatxy, qrfsxy=grid%qrfsxy & & ,qspringsxy=grid%qspringsxy, smoiseq=grid%smoiseq, wtddt=config_flags%wtddt, stepwtd=grid%stepwtd & - & ,gecros_state=grid%gecros_state & ! Optional gecros crop ! Noah UA changes & ,ua_phys=config_flags%ua_phys,flx4=grid%flx4,fvb=grid%fvb & & ,fbur=grid%fbur,fgsn=grid%fgsn & diff --git a/dyn_em/start_em.F b/dyn_em/start_em.F index 08f6f39bbf..c254e232ea 100644 --- a/dyn_em/start_em.F +++ b/dyn_em/start_em.F @@ -1107,7 +1107,6 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & grid%GRAINXY, grid%GDDXY, & ! Optional Noah-MP grid%CROPTYPE, grid%CROPCAT, & ! Optional Noah-MP config_flags%opt_crop, & - grid%gecros_state, & ! Optional gecros crop grid%XSAIXY,grid%LAI, & ! Optional Noah-MP grid%T2MVXY, grid%T2MBXY, grid%CHSTARXY, & ! Optional Noah-MP grid%SMOISEQ ,grid%SMCWTDXY ,grid%RECHXY, grid%DEEPRECHXY, grid%AREAXY, & ! Optional Noah-MP diff --git a/dyn_nmm/module_PHYSICS_CALLS.F b/dyn_nmm/module_PHYSICS_CALLS.F index 65e07ace5a..cc413a580f 100644 --- a/dyn_nmm/module_PHYSICS_CALLS.F +++ b/dyn_nmm/module_PHYSICS_CALLS.F @@ -1592,7 +1592,6 @@ SUBROUTINE TURBL(NTSD,DT,NPHS,RESTRT & & , chv2xy=grid%chv2xy , chb2xy=grid%chb2xy , chstarxy=grid%chstarxy & & , smoiseq=grid%smoiseq, smcwtdxy=grid%smcwtdxy, rechxy=grid%rechxy & & , deeprechxy=grid%deeprechxy & - & ,gecros_state=grid%gecros_state & ! Optional gecros crop & ,coszen=grid%czen,xlat_urb2d=grid%hlat & & ,rs=rs,xlaidyn=xlaidyn & ! mosaic tiling for Noah diff --git a/dyn_nmm/start_domain_nmm.F b/dyn_nmm/start_domain_nmm.F index 53025d959c..b1110d4c48 100644 --- a/dyn_nmm/start_domain_nmm.F +++ b/dyn_nmm/start_domain_nmm.F @@ -2199,7 +2199,6 @@ END SUBROUTINE med_set_egrid_locs & smoiseq=grid%smoiseq, smcwtdxy=grid%smcwtdxy, rechxy=grid%rechxy, & & deeprechxy=grid%deeprechxy, & & iopt_crop=config_flags%opt_crop, & - & gecros_state=grid%gecros_state, & ! Optional gecros crop & lakedepth2d=grid%lakedepth2d, savedtke12d=grid%savedtke12d, snowdp2d=grid%snowdp2d, h2osno2d=grid%h2osno2d, & !lake & snl2d= grid%snl2d, t_grnd2d=grid%t_grnd2d, t_lake3d=grid%t_lake3d, lake_icefrac3d=grid%lake_icefrac3d, & !lake & z_lake3d=grid%z_lake3d, dz_lake3d=grid%dz_lake3d, t_soisno3d=grid%t_soisno3d, h2osoi_ice3d=grid%h2osoi_ice3d, & !lake diff --git a/main/depend.common b/main/depend.common index e5a58206c4..3528c990f0 100644 --- a/main/depend.common +++ b/main/depend.common @@ -463,14 +463,12 @@ module_sf_noahmpdrv.o: module_sf_noahmplsm.o \ module_data_gocart_dust.o \ module_sf_noahmp_glacier.o \ module_sf_noahmp_groundwater.o \ - module_sf_gecros.o \ ../share/module_model_constants.o \ module_sf_urban.o module_sf_bep.o module_sf_bep_bem.o module_sf_noahlsm_glacial_only.o: module_sf_noahlsm.o module_sf_noahmplsm.o module_sf_noahmplsm.o: ../share/module_model_constants.o \ - module_sf_gecros.o \ module_sf_myjsfc.o module_sf_noahmp_groundwater.o: module_sf_noahmplsm.o diff --git a/phys/Makefile b/phys/Makefile index ee46f37441..47583dbb2b 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -138,7 +138,6 @@ MODULES = \ module_sf_noah_seaice_drv.o \ module_sf_noahlsm_glacial_only.o \ module_sf_noahmp_groundwater.o \ - module_sf_gecros.o \ module_sf_noahmpdrv.o \ module_sf_noahmplsm.o \ module_sf_noahmp_glacier.o \ diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index cc37d864b8..23c56d8aed 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -120,7 +120,6 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & GRAINXY, GDDXY, & ! Optional Noah-MP croptype, cropcat, & ! Noah-MP Crop model iopt_crop, & - gecros_state, & ! Optional gecros crop XSAIXY, LAI, & ! Optional Noah-MP T2MVXY, T2MBXY, CHSTARXY , & ! Optional Noah-MP SMOISEQ ,SMCWTDXY ,RECHXY, DEEPRECHXY, AREAXY, & ! Optional Noah-MP @@ -540,7 +539,6 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & INTEGER, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: CROPCAT REAL, OPTIONAL, DIMENSION(ims:ime,5,jms:jme), INTENT(IN ) :: CROPTYPE INTEGER, INTENT(IN ) :: iopt_crop - REAL, OPTIONAL, DIMENSION(ims:ime,60,jms:jme) :: gecros_state ! Optional gecros crop ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, INTENT(INOUT ) :: qnn_curr @@ -1350,7 +1348,6 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & GRAINXY, GDDXY, & ! Noah-MP Crop model croptype, cropcat, & ! Noah-MP Crop model iopt_crop, & - gecros_state, & ! Optional gecros crop XSAIXY, LAI, & SMOISEQ, SMCWTDXY, RECHXY, DEEPRECHXY, AREAXY, & WTDDT, STEPWTD, QRFSXY ,QSPRINGSXY ,QSLATXY, & @@ -2357,7 +2354,6 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & GRAINXY, GDDXY, & ! Noah-MP Crop model croptype, cropcat, & ! Noah-MP Crop model iopt_crop, & - gecros_state, & ! Optional gecros crop XSAIXY, LAI, & SMOISEQ, SMCWTDXY, RECHXY, DEEPRECHXY, AREAXY, & WTDDT, STEPWTD,QRFSXY ,QSPRINGSXY ,QSLATXY, & @@ -2680,7 +2676,6 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & INTEGER, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: CROPCAT REAL, OPTIONAL, DIMENSION(ims:ime,5,jms:jme), INTENT(IN ) :: CROPTYPE INTEGER , INTENT(IN) :: iopt_crop - REAL, OPTIONAL, DIMENSION(ims:ime,60,jms:jme) :: gecros_state ! Optional gecros crop INTEGER, OPTIONAL, INTENT(IN) :: ISICE INTEGER , OPTIONAL, INTENT(OUT) :: STEPWTD REAL , OPTIONAL, INTENT(IN) :: WTDDT @@ -3197,8 +3192,7 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & ,smoiseq ,smcwtdxy ,rechxy ,deeprechxy, areaxy ,dx, dy, msftx, msfty,& wtddt ,stepwtd ,dt ,qrfsxy ,qspringsxy ,qslatxy, & fdepthxy ,ht ,riverbedxy ,eqzwt ,rivercondxy ,pexpxy, & - rechclim, & - gecros_state & ! Optional gecros crop + rechclim & #endif ) diff --git a/phys/module_sf_gecros.F b/phys/module_sf_gecros.F deleted file mode 100644 index aa02408401..0000000000 --- a/phys/module_sf_gecros.F +++ /dev/null @@ -1,2411 +0,0 @@ -!*********************************************************************** -!* GECROS for early and late covering crops * -!* Genotype-by-Environment interaction on CROp growth Simulator * -!* * -!* Author: Xinyou YIN * -!* Crop and Weed Ecology Group * -!* Wageningen University & Research Centre * -!* PO Box 430, 6700 AK Wageningen, Netherlands * -!* * -!* Modified and extended for winter crops * -!* by Joachim INGWERSEN * -!* Biogeophysics Group * -!* University of Hohenheim * -!* * -!*********************************************************************** - -MODULE module_sf_gecros - -implicit none - -!Runtime constant variables defined within this subroutine -!Crop-specific parameters (former *.tbl file) -REAL, PARAMETER :: EG = 1. !Efficiency of germination -REAL, PARAMETER :: CFV = 0.48 !Carbon fraction in vegetative organs -REAL, PARAMETER :: YGV = 0.81 !Grwoth efficiency for vegetative growth -REAL, PARAMETER :: FFAT = 0.02 !Fraction of fat in storage organs -REAL, PARAMETER :: FLIG = 0.06 !Fraction of lignin in storage organs -REAL, PARAMETER :: FOAC = 0.02 !Fraction of organic acids in storage organs -REAL, PARAMETER :: FMIN = 0.02 !Fraction of minerals in storage organs -REAL, PARAMETER :: LNCI = 0.03586 !Initial value of LNC (g N g-1) -REAL, PARAMETER :: TBD = 0.0 !Base temperature of phenology (C) -REAL, PARAMETER :: TOD = 22.5 !Optimum temperature of phenology (C) -REAL, PARAMETER :: TCD = 37.0 !Ceiling temperature of phenology (C) -REAL, PARAMETER :: TSEN = 1.0 !Curvature of temperature response (1) -REAL, PARAMETER :: SPSP = 0.2 !DS for start of photoperiod-sensitive phase -REAL, PARAMETER :: EPSP = 0.7 !DS for end of photoperiod-sensitive phase -REAL, PARAMETER :: CDMHT = 492.6 !Stem dry weight per unit plant height (g m-2 m2) -REAL, PARAMETER :: PMEH = 0.6468 !Fraction of sigmoid curve inflexion in plant height growth period (1) -REAL, PARAMETER :: ESDI = 1.35 !ESD for indeterminate crops -REAL, PARAMETER :: PMES = 0.50 !Fraction of sigmoid curve inflexion in seed growth growth period (1) -REAL, PARAMETER :: TBDV = -1.3 !Base temperature of vernalization (C) -REAL, PARAMETER :: TODV = 4.9 !Optimum temperature of vernalization (C) -REAL, PARAMETER :: TCDV = 15.7 !Ceiling temperature of vernalization (C) -REAL, PARAMETER :: NUPTX = 0.5/24./3600. !Maximum crop nitrogen uptake (g N m-2 s-1) -REAL, PARAMETER :: SLA0 = 0.0237 !Specific leaf area constant (m2 leaf g-1) -REAL, PARAMETER :: SLNMIN = 0.35 !Base SLN for photosynthesis (g N m-2 leaf) -REAL, PARAMETER :: RNCMIN = 0.005 !Minimum N concentration in roots (g N g-1) -REAL, PARAMETER :: STEMNC = 0.01 !Minimum N concentration in stem (g N g-1) -REAL, PARAMETER :: RLVDS = 0.0904 !Rate of C turnover from dead leaves to litter (d-1) -REAL, PARAMETER :: DSCRIT = 0.225 !DS at which winter dormancy ends (1) -REAL, PARAMETER :: EAJMAX = 48270.!Energy of activation for JMAX (J mol-1) -REAL, PARAMETER :: XVN = 24.96 !Slope of linear relationship between JMAX and leaf N (mu-mol CO2 s-1 g-1 N) -REAL, PARAMETER :: XJN = 49.92 !Slope of linear relationship between VCMX and leaf N (mu-mol e- s-1 g-1 N) -REAL, PARAMETER :: NPL = 390.0 !Number of plants per m2 -REAL, PARAMETER :: SEEDW = 0.0475 !Seed weight (g seed-1) -REAL, PARAMETER :: SEEDNC = 0.020 !Seed N concentration (g N g-1) -REAL, PARAMETER :: BLD = 25.58 !Leaf angle (from horizontal) (degree) -REAL, PARAMETER :: HTMX = 1.1 !Maximum plant height (m) -REAL, PARAMETER :: MTDV = 46.12 !Minimum thermal days for vegetative growth phase (d) -REAL, PARAMETER :: MTDR = 40.22 !Minimum thermal days for reproductive phase (d) -REAL, PARAMETER :: PSEN = -.104 !Photoperiod sensitivity of phenological development (h-1) -REAL, PARAMETER :: C3C4 = 1. !C3C4 = 1. for C3 crops; C3C4 = -1. for C4 crops. - -!Runtime-constant general parameters -REAL, PARAMETER :: LEGUME = -1. !LEGUME = 1. for leguminous crops; = -1. for non-leguminous crops. -REAL, PARAMETER :: DETER = 1. !DETER = 1. for determinate crops; = -1. for indeterminate crops. -REAL, PARAMETER :: SLP = -1. !SLP = 1. for short-day crops; = -1. for long-day crops. -REAL, PARAMETER :: NSUP1 = 0. !NH4-N supply in g/m2/s -REAL, PARAMETER :: NSUP2 = 1./172800. !NO3-N supply in g/m2/s -REAL, PARAMETER :: CO2A=350. !Atmospherric CO2 concentration (ppm) -REAL, PARAMETER :: FCRSH=0.5 !Initial fraction of C in shoot (-) -REAL, PARAMETER :: FNRSH=0.63 !Initial fraction of N in shoot (-) -REAL, PARAMETER :: PNPRE=0.7 !Proportion of seed N that comes from non-structural stems during seed fill -REAL, PARAMETER :: CB=0.75 !Factor for initial N concentration of seed fill -REAL, PARAMETER :: CX=1.00 !Factor for final N concentration of seed fill -REAL, PARAMETER :: TM=1.5 !DS when transition from CB to CX is fastest -REAL, PARAMETER :: RSS=100. !Soil resistance to evaporation (s/m) -REAL, PARAMETER :: LS=0. !Lodging sverity (value between zero and unity) -REAL, PARAMETER :: WRB=0.25 !Critical root weight density (g/m2/cm depth) -REAL, PARAMETER :: THETA=0.7 !Convexity for light response curve of electron transport (J2) in photosynthesis -REAL, PARAMETER :: PNLS=1. !Fraction of dead leaf N incorporated into soil litter N -REAL, PARAMETER :: INSP=-2. !Inclination of sun angle (degree) -REAL, PARAMETER :: CCFIX=6. !Carbon cost of symbiotic N fixation (g C/g N) -REAL, PARAMETER :: RDMX=130. !Maximum rooting depth (m) -REAL, PARAMETER :: TCP=86400. !Conversion factor from day into sec -REAL, PARAMETER :: Z1244=0.272727272727 !12./44. -REAL, PARAMETER :: LOG005=-2.995732274 !LOG(0.05) - -!Runtime constant variables computed in driver.f -REAL :: YGO !Growth efficiency of storage organs (g C/g C) -REAL :: CFO !Carbon fraction in storage organs (g C/g) -REAL :: LNCMIN !Minimum N concentration in leaves (g N/g) -REAL :: CLVI !Initial value of CLV -REAL :: CRTI !Initial value of CRT -REAL :: NLVI !Initial value of NLV -REAL :: NRTI !Initial value of NRT -REAL :: HTI !Initial value of HT -REAL :: RDI !Initial value of RD - -!*** Debugging on/off -LOGICAL :: debugging=.false. - -CONTAINS - -SUBROUTINE gecros (DOY, DT, CROP, RB, RT, RTS, FB, SNOWH, & !I - WN, SFCTMP, EAIR, RSD, RLD, PRCP, WUL, WLL, WCMIN, LWIDTH, & !I & !I - STATE_GECROS, & !H - ATRJC, ATRJS, FSR, FRSU, ARSWSU, ARSWSH) !O - -IMPLICIT NONE - -! character(len=19), INTENT(IN) :: nowdate ! string of current date e.g. 2012-30-10_12:00:00 - INTEGER, INTENT(IN) :: CROP !CROP=1 -> early-covering crop, CROP=2 -> later covering crop - REAL, INTENT(IN) :: DOY !Julian day of the current time step - REAL, INTENT(IN) :: DT !Integration time step (s) - REAL, INTENT(IN) :: RB !Leaf boundary layer resistance (s/m) - REAL, INTENT(IN) :: RT !Aerodynamic canopy resistance (s/m) - REAL, INTENT(IN) :: RTS !Aerodynamic ground resistance (s/m) - REAL, INTENT(IN) :: FB !Fraction of vegetation cover by snow - REAL, INTENT(IN) :: SNOWH !Snow height (cm) - REAL, INTENT(IN) :: WN !Wind speed (m/s) - REAL, INTENT(IN) :: SFCTMP !Air temperature (K) - REAL, INTENT(IN) :: EAIR !Vapour pressure (Pa) - REAL, INTENT(IN) :: RSD !Downwelling shortwave radiation (W/m2) - REAL, INTENT(IN) :: RLD !Downwelling longwave raidation (W/m2) - REAL, INTENT(IN) :: PRCP !Precipiration (kg/m2/s) - REAL, INTENT(IN) :: WUL !Soil water in the upper layer (L/m2 or mm) - REAL, INTENT(IN) :: WLL !Soil water in the lower layer (L/m2 or mm) - REAL, INTENT(IN) :: WCMIN !Minimum soil watter content (m3/m3) - REAL, INTENT(IN) :: LWIDTH !Leaf width (m) read in from MPTABLE.TBL - REAL, INTENT(OUT) :: ATRJC !Absorbed global radiation by canopy (W/m2) - REAL, INTENT(OUT) :: ATRJS !Absorbed global radiation by soil (W/m2) - REAL, INTENT(OUT) :: FSR !Reflected shortwave radiation (W/m2) - REAL, INTENT(OUT) :: FRSU !Fraction of sunlit leaves in the canopy - REAL, INTENT(OUT) :: ARSWSU !Actual stomatal resistance of sunlit leaves (s/m) - REAL, INTENT(OUT) :: ARSWSH !Actual stomatal resistance of shaded leaves (s/m) - REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS !Array with Gecros state variables - - !*** Gecros state variables - REAL :: DS !Development stage - REAL :: LAI !Green leaf area index - REAL :: TLAI !Total leaf area index - REAL :: CTDU !Cumulative thermal day unit (d) - REAL :: CVDU !Cumulative thermal day unit of vernalization (d) - REAL :: CLV !Carbon in living leaves (g C/m2) - REAL :: CLVD !Carbon in dead leaves (g C/m2) - REAL :: CSST !Carbon in structural stems (g C/m2) - REAL :: CSO !Carbon in storage organs (g C/m2) - REAL :: CSRT !Carbon in structural roots (g C/m2) - REAL :: CRTD !Carbon in dead roots (g C/m2) - REAL :: CLVDS !Carbon in dead leaves that have become litter in soil (g C/m2) - REAL :: NRT !Nitrogen in living roots (g N/m2) - REAL :: NST !Nitrogen in stems (g N/m2) - REAL :: NLV !Nitrogen in living leaves (g N/m2) - REAL :: NSO !Nitrogen in storage organs (g N/m2) - REAL :: TNLV !Total leaf nitrogen (including N in senesced leaves) (g N/m2) - REAL :: NLVD !Nitrogen in dead leaves (g N/m2) - REAL :: NRTD !Nitrogen in dead roots (g N/m2) - REAL :: CRVS !Carbon in stem reserves (g C/m2) - REAL :: CRVR !Carbon in root reserves (g C/m2) - REAL :: NREOE !NRES accumulated till the end of seed number determining period (g N/m2) - REAL :: NREOF !NRES accumulated till the moment at which seed fill starts (g N/m2) - REAL :: DCDSR !Short fall of carbon demand for seed fill in previous time steps (g C/m2) - REAL :: DCDTR !Short fall of carbon demand for structural stems in previous time steps (g C/m2) - REAL :: SLNB !SLN in bottom leaves of canopy (g N/m2) - REAL :: LAIC !Carbon determined LAI (m2/m2) - REAL :: RMUL !Total of RMUN+RMUA+RMUS+RMLD - REAL :: NDEMP !Crop nitrogen demand of the previous time step (g N/m2/s) - REAL :: NSUPP !Nitrogen supply of the previous time step (g N/m2/s) - REAL :: NFIXT !Total symbiotically fixed nitrogen during growth (g N/m2) - REAL :: NFIXR !Reserve pool of symbiotically fixed nitrogen (g N/m2) - REAL :: DCDTP !Carbon demand for structural stem growth of the previous time step (g C/m2/s) - REAL :: HT !Canopy height (m) - REAL :: TPCAN !Cumulative canopy photosynthesis over growth period (g CO2/m2) - REAL :: TRESP !Total crop respiratory cost during growth (g CO2/m2) - REAL :: TNUPT !Total crop nitrogen uptake during growth (g N/m2) - REAL :: LITNT !Total litter nitrogen entering soil during growth (g N/m2) - REAL :: WSO !Yield (g/m2) - REAL :: WSTRAW !Straw (g/m2) - REAL :: GrainNC !Nitrogen in grain (kg N/ha) - REAL :: StrawNC !Nitrogen in straw (kg N/ha) - REAL :: APCAN !Actual gross canopy photosynthesis (g CO2/m2/d) - - character(Len=12) :: inputstring - INTEGER :: nowday, minutes - REAL :: SC !Solar constant (W/m2) - REAL :: SINLD !Seasonal offset of sine of solar height (-) - REAL :: COSLD !Amplitude of sine of solar height (-) - REAL :: DAYL !Day length (h) - REAL :: DDLP !Day length for photoperiodism (h) - REAL :: DSINBE !DSINB to correct for lower atmospheric transmission at lower solar elevation (s/d) - REAL :: DVP !Vapour pressure (kPa) - REAL :: WNM !Wind speed (m/s) - REAL :: TAIR !Air temperature (C) - REAL :: ROOTD !Rooting depth (m) - REAL :: WLV !Dry weight of living leaves (g/m2) - REAL :: WST !Dry weight of stems (g/m2) - REAL :: WRT !Dry weight of roots (g/m2) - REAL :: WSH !Dry weight of shoot (g/m2) - REAL :: WLVD !Dry weight of dead leaves (g/m2) - REAL :: WRTD !Dry weight of dead roots (g/m2) - REAL :: CRT !Carbon in living roots (g C/m2) - REAL :: CSH !Carbon in living shoot (g C/m2) - REAL :: NSH !Nitrogen in living shoot (g N/m2) - REAL :: DLAI !LAI of dead leaves still hanging on stem (m2/m2) - REAL :: ESD !DS for end of seed-number determing phase (-) - REAL :: KCRN !EXtinction coefficient of root nitrogen (m2/g C) - REAL :: NFIXD !Crop demand-determined NFIX (g N/2/s) - REAL :: KR !Extinction coefficient of root weight density over soil septh (1/cm) - REAL :: HNCCR !Critical shoot N concentration ( g N/g) - REAL :: FVPD !Slope of linear effect of VPD of intercelluar to ambient CO2 ratio (1/kPa) - REAL :: NRETS !Total crop-residue N returned to soil (g N/m2) - REAL :: WCUL !Soil water content in uppler layer (m3/m3) - REAL :: DWSUP !Water supply for evapotranspiration (mm/s) - REAL :: CSRTN !N determined CSRT (g C/m2) - REAL :: NRES !Estimated vegetative-organ N remobilizable for seed fill (g N/m2) - REAL :: ONC !N concentration in storage organ (g N/g) - REAL :: RNC !N concentration in roots (g N/g) - REAL :: LNC !N concentration in living leaves (g N/g) - REAL :: KL !Extinction coefficient diffuse component for PAE (m2/m2 leaf) - REAL :: CTOT !Total C in living shoots and roots (g C/m2) - REAL :: NSHH !N in shoots (excluding dead leaves incorporated into soil litter (g N/m2) - REAL :: NTOT !Total N in living shoots and roots (g N/m2) - REAL :: WSHH !Dry weight of shoots organs (excluding shadded leaves) (g/m2) - REAL :: TSN !Total seed number (seeds/m2) - REAL :: HNC !Actural N concentration in living shoot (g N/g) - REAL :: PSO !Protein content in storage organs (g protein/m2) - REAL :: KLN !Intermediate variable to compute KN (g N/m2 leaf) - REAL :: NBK !Intermediate variable to compute KN (g N/m2 leaf) - REAL :: KW !Wind speed extinction coefficient in canopy (m2/m2 leaf) - REAL :: WTOT !Dry weight of total living plant parts (g/m2) - REAL :: CCHKIN !C in crop accumulated since start of simulation - REAL :: NCHKIN !N in crop accumulated since start of simulation - REAL :: LCRT !Rate of C loss in rooots because of senescence (g C/m2/s) - REAL :: TSW !Thousand seed weight (g) - REAL :: PNC !N concentration in living plant material (g N/m2) - REAL :: NDEMD !Defiency-driven N demand (g N/m2/s) - REAL :: FCRVR !Fraction of new root C partitioned to root reserves (g C/g C) - REAL :: KN !Leaf N extinction coefficient in the canopy (m2/m2 leaf) - REAL :: CCHK !Diff. between C added to the crop since start and net total C fluxes rel. to CCHKIN (%) - REAL :: HI !Harvest index (g/g) - REAL :: LWRT !Rate of loss of root biomass because of senescence (g/m2/s) - REAL :: NCHK !As CCHK but for N - REAL :: LAIN !N-determined LAI (m2/m2) - REAL :: LNRT !Rate of loss of N because of senescence (g N/m2/s) - REAL :: LWLVM !Intermediate variable to compute LWLV (g/m2/s) - REAL :: SLNNT !Value of SLNT with small plant-N increment (g N/m2 leaf) - REAL :: SLNBC !SLNB calculated from exponential N profile in canopy (g N/m2 leaf) - REAL :: SLN !Specific leaf N content (g N/m2 leaf) - REAL :: SLA !Specific leaf area (m2 leaf/g) - REAL :: LWLV !Rate of loss of leaf weight because of leaf senescence (g/m2/s) - REAL :: AESOIL !Actual soil evaporation (mm/s) - REAL :: APCANN !APCANS with small plant-N increment (g CO2/m2/s) - REAL :: APCANS !Actual standing canopy CO2 assimilation (g CO2/m2/s) - REAL :: ATCAN !Actual canopy transpiration (mm/s) - REAL :: DAPAR !PAR absorbed by canopy (J/m2/s) - REAL :: DIFS !Daytime average soil-air temperature difference (C) - REAL :: DIFSH !Daytime average shaded leaf-air temperature difference (C) - REAL :: DIFSU !Daytime average sunlit leaf-air temperature difference (C) - REAL :: FRSH !Fraction of shaded leaves in the canopy - REAL :: HOD !Hour of the day - REAL :: PESOIL !Potential soil evaporation (mm/s) - REAL :: PPCAN !Potential canopy assimilation (g CO2/m2/s) - REAL :: PTCAN !Potential canopy transpiration (mm/s) - REAL :: RCAN !Canopy resistance (s/m) - REAL :: RSLNB !Rate of change in SLNB (g N/m2 leaf/s) - REAL :: LNLV !Rate of loss of leaf N because of senescence (g N/m2/s) - REAL :: LCLV !Rate of loss of leaf C because of senescence (g C/m2/s) - REAL :: RMRE !Residual maintenance respiration (g CO2/m2/s) - REAL :: TAVSS !Soil surface temperature (C) - REAL :: RMN !RM calculated with a small increment in plant N (g CO2/m2/s) - REAL :: VDU !Rate of Vernalization day unit increase (d/s) - REAL :: TDU !Rate of thermal day unit increase (d/s) - REAL :: RM !Non-growth components of respiration, excluding the cost of N fixation (g CO2/m2/s) - REAL :: LVDS !Rate of transfer of C from dead leaves to litter (g C/m2/s) - REAL :: NFIXE !Available energy-determined NFIX (g N/m2/s) - REAL :: DVR !Phenological development rate (1/s) - REAL :: RX !Respiratory cost of N fixation (g CO2/m2/s) - REAL :: RNSUPP !Rate of change in NSUPP (g N/m2/s) - REAL :: NSUP !N supply to crop (g N/m2/s) - REAL :: NFIX !Symbioticall fixed N (g N/m2/s) - REAL :: LITN !Litter N entering soil (g N/m2/s) - REAL :: LITC !Litter C entering soil (g C/m2/s) - REAL :: FDH !Expected relative growth rate of plant height (1/s) - REAL :: FDS !Expected relative growth rate of storage organs (1/s) - REAL :: SHSAN !SHSA calculated with a small increment in plant-N (g C/g C/s) - REAL :: SHSA !Relative shoot activity (g C/g C/s) - REAL :: RMUS !Respiratory cost of mineral uptake (g CO2/m2/s) - REAL :: NDEMAD !Intermediate variable related to NDEM (g N/m2/s) - REAL :: NDEMA !Activity-driven NDEM (g N/m2/s) - REAL :: NCR !Intermediate variable (g N/g C) - REAL :: DERI !First order-derivative of SHSA with respect to crop N/C ratio (g C/g N/s) - REAL :: ASSA !Assimilates available from current photosynthesis for growth (g CO2/m2/s) - REAL :: RNFIXR !Rate of change in NFIXR (g N/m2/s) - REAL :: RNDEMP !Rate of change in NDEMP (g N/m2/s) - REAL :: RMLD !Respiration due to phloem loading of C assimilates to root (g CO2/m2/s) - REAL :: RCSRT !Rate of change in CSRT (g C/m2/s) - REAL :: NUPTN !Nitrate-N uptake by the crop (g N/m2/s) - REAL :: NUPTA !Ammonium-N uptake by the crop (g N/m2/s) - REAL :: NDEM !Crop N demand (g N/m2/s)SUBROUTINE ENERGY - REAL :: FNSH !Fraction of newly absorbed N partitioned to shoot (g N/g N) - REAL :: FCSH !Fraction of new V partitioned to shoot (g C/g C) - REAL :: DCSR !C supply from current photosynthesis for root growth (g C/m2/s) - REAL :: SLNT !SLN for top leaves in canopy (g N/m2 leaf) - REAL :: RMUN !Respiratory cost of nitrate-N uptake (g CO2/m2/s) - REAL :: RMUA !Respiratory cost of ammonium-N uptake (g CO2/m2/s) - REAL :: NUPT !Crop N uptake (g N/m2/s) - REAL :: DCDSC !C demand for seed filling at current time step (g C/m2/s) - REAL :: DCDS !C demand for filling of storage organgs at current time step (g C/m2/s) - REAL :: FLWCS !Flow of current assimilated C to storage organs (g C/m2/s) - REAL :: DCSS !C supply from current photosynthesis for shoot growth (g C/m2/s) - REAL :: DCST !C suppyl from current photosynthesis for structural stem growth (g C/m2/s) - REAL :: FCSO !Fraction of new C partitioned to storage organs (g C/m2/s) - REAL :: RRMUL !Rate of change in RMUL (g CO2/m2/s) - REAL :: RHT !Rate of change in HT (m/s) - REAL :: IFSH !Integral factor of stresses on plant height growth (-) - REAL :: GAP !Gap betweeb C supply and C demand for seed growth (g C/m2/s) - REAL :: CREMSI !Intermediate variable to compute CREMS (g C/m2/s) - REAL :: CREMS !C remobilized from stem reserves to storage organs (g C/m2/s) - REAL :: CREMRI !Intermediate variable to compute CREMR (g C/m2/s) - REAL :: CREMR !C remobilized from root reserves to storage organs (g C/m2/s) - REAL :: RWSO !Rate of change in storage organs(g/m2/s) - REAL :: RWRT !Rate of change in roots (g/m2/s) - REAL :: RRD !Rate of change in RD (cm/s) - REAL :: RDCDTP !Rate of change in DCDTP (g C/m2/s) - REAL :: RDCDSR !Rate of change in DCDSR (g C/m2/s) - REAL :: RCSST !Rate of change in CSST (g C/m2/s) - REAL :: RCSO !Rate of change in CSO (g C/m2/s) - REAL :: RCRVR !Rate of change in CRVR (g C/m2/s) - REAL :: FLWCT !Flow of assimilated C to structural stems (g C/m2/s) - REAL :: FCSST !Fraction of new shoot C partitioned to structural stems (g C/m2/s) - REAL :: FCLV !Fraction of new shoot C partitioned to leaves (g C/m2/s) - REAL :: DCDTC !C demand of structural stem growth at current time step (g C/m2/s) - REAL :: DCDT !C demand for the growth of structural stems (g C/m2/s) - REAL :: FCRVS !Fraction of new C paritioned to stem reserves (g C/m2/s) - REAL :: RCLV !Rate of change in CLV (g C/m2/s) - REAL :: RCRVS !Rate of change in CRVS (g C/m2/s) - REAL :: RDCDTR !Rate of change in DCDTR (g C/m2/s) - REAL :: RESTOT !Total respiratory cost (g CO2/m2/s) - REAL :: RG !Growth respiration (g CO2/m2/s) - REAL :: RLAI !Rate of change in LAI (m2 leaf/m2/s) - REAL :: RNLV !Rate of change in NLV (g N/m2/s) - REAL :: RNREOE !Rate of change in NREOE (g N/m2/s) - REAL :: RNREOF !Rate of change in NREOF (g N/m2/s) - REAL :: RNRES !Rate of change in NRES (g N/m2/s) - REAL :: RNRT !Rate of change in NRT (g N/m2/s) - REAL :: RNSO !Rate of change in NSO (g N/m2/s) - REAL :: RNST !Rate of change in NSZ (g N/m2/s) - REAL :: RRP !Respiration/photosynthesis ratio (-) - REAL :: RTNLV !Rate of change in TNLV (g N/m2/s) - REAL :: RWLV !Rate of change in WLV (g/m2/s) - REAL :: RWST !Rate of change in WST (g/m2/s) - REAL :: GLAT !Latitude (degree) - REAL :: SD1 !Thickness of upper evaporative soil layer (cm) (equals RDI) - - !*** write STATE_GECROS array into Gecros variables - DS = STATE_GECROS(1) - CTDU = STATE_GECROS(2) - CVDU = STATE_GECROS(3) - CLV = STATE_GECROS(4) - CLVD = STATE_GECROS(5) - CSST = STATE_GECROS(6) - CSO = STATE_GECROS(7) - CSRT = STATE_GECROS(8) - CRTD = STATE_GECROS(9) - CLVDS = STATE_GECROS(10) - NRT = STATE_GECROS(11) - NST = STATE_GECROS(12) - NLV = STATE_GECROS(13) - NSO = STATE_GECROS(14) - TNLV = STATE_GECROS(15) - NLVD = STATE_GECROS(16) - NRTD = STATE_GECROS(17) - CRVS = STATE_GECROS(18) - CRVR = STATE_GECROS(19) - NREOE = STATE_GECROS(20) - NREOF = STATE_GECROS(21) - DCDSR = STATE_GECROS(22) - DCDTR = STATE_GECROS(23) - SLNB = STATE_GECROS(24) - LAIC = STATE_GECROS(25) - RMUL = STATE_GECROS(26) - NDEMP = STATE_GECROS(27) - NSUPP = STATE_GECROS(28) - NFIXT = STATE_GECROS(29) - NFIXR = STATE_GECROS(30) - DCDTP = STATE_GECROS(31) - HT = STATE_GECROS(32) - ROOTD = STATE_GECROS(33) - TPCAN = STATE_GECROS(34) - TRESP = STATE_GECROS(35) - TNUPT = STATE_GECROS(36) - LITNT = STATE_GECROS(37) - GLAT = STATE_GECROS(44) - WSO = STATE_GECROS(45) - WSTRAW = STATE_GECROS(46) - GrainNC = STATE_GECROS(47) - StrawNC = STATE_GECROS(48) - LAI = STATE_GECROS(49) - TLAI = STATE_GECROS(50) - SD1 = STATE_GECROS(52) - - ! Used for debugging -! if (nowdate(9:12).eq.'2500') then -! write(*,*) nowdate, ' 1 ', DS -! write(*,*) nowdate, ' 2 ', CTDU -! write(*,*) nowdate, ' 3 ', CVDU -! write(*,*) nowdate, ' 4 ', CLV -! write(*,*) nowdate, ' 5 ', CLVD -! write(*,*) nowdate, ' 6 ', CSST -! write(*,*) nowdate, ' 7 ', CSO -! write(*,*) nowdate, ' 8 ', CSRT -! write(*,*) nowdate, ' 9 ', CRTD -! write(*,*) nowdate, ' 10 ', CLVDS -! write(*,*) nowdate, ' 11 ', NRT -! write(*,*) nowdate, ' 12 ', NST -! write(*,*) nowdate, ' 13 ', NLV -! write(*,*) nowdate, ' 14 ', NSO -! write(*,*) nowdate, ' 15 ', TNLV -! write(*,*) nowdate, ' 16 ', NLVD -! write(*,*) nowdate, ' 17 ', NRTD -! write(*,*) nowdate, ' 18 ', CRVS -! write(*,*) nowdate, ' 19 ', CRVR -! write(*,*) nowdate, ' 20 ', NREOE -! write(*,*) nowdate, ' 21 ', NREOF -! write(*,*) nowdate, ' 22 ', DCDSR -! write(*,*) nowdate, ' 23 ', DCDTR -! write(*,*) nowdate, ' 24 ', SLNB -! write(*,*) nowdate, ' 25 ', LAIC -! write(*,*) nowdate, ' 26 ', RMUL -! write(*,*) nowdate, ' 27 ', NDEMP -! write(*,*) nowdate, ' 28 ', NSUPP -! write(*,*) nowdate, ' 29 ', NFIXT -! write(*,*) nowdate, ' 30 ', NFIXR -! write(*,*) nowdate, ' 31 ', DCDTP -! write(*,*) nowdate, ' 32 ', HT -! write(*,*) nowdate, ' 33 ', ROOTD -! write(*,*) nowdate, ' 34 ', TPCAN -! write(*,*) nowdate, ' 35 ', TRESP -! write(*,*) nowdate, ' 36 ', TNUPT -! write(*,*) nowdate, ' 37 ', LITNT -! read(*,*) -! endif - - PPCAN=0. - APCANS=0. - APCANN=0. - APCAN=0. - PTCAN=0. - ATCAN=0. - PESOIL=0. - AESOIL=0. - DIFS=0. - DIFSU=0. - DIFSH=0. - DAPAR=0. - RCAN=0. - DVR=0. - - nowday = INT(DOY) - HOD = float(nint((DOY-int(DOY))*86400.))/3600. - - ! Conversion from K to C - TAIR = SFCTMP - 273.15 - - ! Conversion of rel. humidity into VP (kPa) - DVP = EAIR*0.001 !Converts EAIR from Pa to kPa - WNM = MAX (0.1, WN) - - ! Photoperiod, solar constant and daily extraterrestrial radiation - CALL ASTRO(aint(DOY),GLAT,INSP,SC,SINLD,COSLD,DAYL,DDLP,DSINBE) - - ! Plant weights (g/m2) - WLV = CLV / CFV - WST = CSST / CFV + CRVS/0.444 - WSO = CSO / CFO - WRT = CSRT / CFV + CRVR/0.444 - WLVD = CLVD / CFV - WRTD = CRTD / CFV - - ! Carbon in shoot and root (g C/m2) - CSH = CLV + CSST + CRVS + CSO - CRT = CSRT + CRVR - - ! Nitrogen in shoot (g C/m2) - NSH = NST + NLV + NSO - - ! Extinction coefficient of root nitrogen (m2/g C) - KCRN = -LOG005/6.3424/CFV/WRB/RDMX - - ! DS for end of seed number determining period - ESD = INSW(DETER, ESDI, 1.) - - ! Dead leaves still hanging on the plant (m2/m2 - DLAI = (CLVD-CLVDS)/CFV*SLA0 - - ! Total leaf area index (dead plus living leaves) - TLAI = MAX(0.01,LAI + DLAI) - - ! Crop demand-determined NFIX - NFIXD = MAX(0., NDEMP - NSUPP) - - ! Critical shoot nitrogen concentration g N g-1 - HNCCR = LNCI*EXP(-.4*DS) - - !Extinction coefficient of root weight density over the soil depth cm-1 - KR = -LOG005/RDMX !cm-1 - - ! Slope of linear effect of VPD on intercelluar to ambient CO2 ratio (1/kPa) - FVPD = INSW (C3C4, 0.195127, 0.116214) - - ! Total crop-residue nitrogen returned to soil (g N/m2) - NRETS = LITNT+INSW(DS-2.,0.,NLV+NST+NRT+NFIXR+(CLVD-CLVDS)/ & - CFV*LNCMIN*(1.+PNLS)/2.) - - !if (nowdate(1:12).eq.'201110061930') then - !write(*,*) 'Bis hier her bin ich gekommen', ROOTD - !endif - - ! Soil water content of the upper and lower layer (m3/m3) - WCUL = (WUL+WCMIN*10.*ROOTD)/10./ROOTD - - ! Daily water supply for evapotranspiration (mm/s) - DWSUP = MAX(.1/TCP,WUL/TCP+.1/TCP) - - ! Nitrogen-determined CSRT (g C/m2) - CSRTN = 1./KCRN*LOG(1.+KCRN*MAX(0.,(NRT*CFV-CRVR*RNCMIN))/RNCMIN) - - ! Straw weight (g/m2) - WSTRAW = WLV + WST + (WLVD - CLVDS/CFV) - - ! Shoot weight (g/m2) - WSH = WLV + WST + WSO - - ! Estimated vegetative-organ N remobilizable for seed growth (g N/m2) - NRES = NREOF + (NREOE-NREOF)*(ESD-1.)/NOTNUL(MIN(DS,ESD)-1.) - - ! Nitrogen concentration in living leaves LAI(g N/m2) - LNC = NLV / NOTNUL(WLV) - - ! Nitrogen concentration in roots (g N/m2) - RNC = NRT / NOTNUL(WRT) - - ! Nitrogen concentration seeds (g N/m2) - ONC = INSW(-WSO, NSO/NOTNUL(WSO), 0.) - - ! Nitrogen in grain and straw in kg N per ha - GrainNC = NSO*10. - StrawNC = (NLV+NST)*10. - - ! Extinction coefficient of nitrogen and wind - CALL KDIFF (TLAI,BLD*3.141592654/180.,0.2, KL) - - CTOT = CSH + CRT - NSHH = NSH +(WLVD-CLVDS/CFV)*LNCMIN - NTOT = NSH + NRT - WSHH = WSH + (WLVD-CLVDS/CFV) - TSN = NRES/PNPRE/SEEDNC/SEEDW - HNC = NSH / NOTNUL(WSH) - - ! Amount of seed protein - PSO = 6.25*WSO*ONC - - KLN = KL*(TNLV-SLNMIN*TLAI) - - NBK = SLNMIN*(1.-EXP(-KL*TLAI)) - KW = KL - WTOT = WSH + WRT - - ! Crop carbon balance check - CCHKIN = CTOT + CLVD + CRTD -CLVI-CRTI - - ! Crop nitrogen balance check - NCHKIN = NTOT + NLVD + NRTD -NLVI-NRTI - - LCRT = MAX(MIN(CSRT-1.E-4,CSRT-MIN(CSRTN,CSRT)),0.)/TCP !Eq. 43, p. 38, DELT ok, gC m-2 s-1 - - FCRVR = INSW(CSRTN-CSRT, 1., 0.) - PNC = NTOT / NOTNUL(WTOT) - - KN = 1./TLAI*LOG(MAX(1.001,(KLN+NBK)/(KLN*EXP(-KL*TLAI)+NBK))) - - TSW = WSO/NOTNUL(TSN)*1000. - - NDEMD = INSW(DS-1., WSH*(HNCCR-HNC)*(1.+NRT/MAX(1E-2,NSH))/TCP, 0.) !Eq. 20, p.25, DELT ok, g N m-2 s-1 - - ! Biomass formation - HI = WSO / NOTNUL(WSHH) - CCHK = (CCHKIN-(TPCAN-TRESP)*Z1244)/NOTNUL(CCHKIN)*100. - NCHK = (NCHKIN-TNUPT)/NOTNUL(TNUPT)*100. - LWRT = LCRT/CFV !g m-2 s-1 - - ! Leaf area development - LAIN = LOG(1.+ KN*MAX(0.,NLV)/SLNMIN)/KN - LNRT = LWRT*RNCMIN !g N m-2 s-1 - - ! Green leaves still hanging on the plant (m2/m2) - LAI = MAX(0.01,MIN(LAIN, LAIC)) - - ! Leaf senescence - ! The equation differs somewhat from Yin. The right hand side is not divided by deltaT - ! This is done by computing LWLV - ! LWLVM = (LAIC-MIN(LAIC,LAIN))/SLA0/DELT !Eq. 42, p.36, DELT ok, g m-2 s-1 - - LWLVM = (LAIC-MIN(LAIC,LAIN))/SLA0 !Eq. 42, p.36, DELT ok, g m-2 s-1 - - ! Specific leaf nitrogen and its profile in the canopy - SLN = NLV/LAI - SLNT = NLV*KN/(1.-EXP(-KN*LAI)) - SLNBC = NLV*KN*EXP(-KN*LAI)/(1.-EXP(-KN*LAI)) - SLNNT = (NLV+0.001*NLV)*KN /(1.-EXP(-KN*LAI)) - SLA = LAI/NOTNUL(WLV) - - ! ji, 1402, Stay-green maize -> no senescense - LWLV = MIN((WLV-1.E-5)/TCP, (LWLVM+REANOR(ESD-DS,LWLVM)*0.03*WLV)/TCP) !g m-2 s-1, gecheckt ok - - !write(*,*) nowdate, DWSUP,CO2A,LS,EAJMAX,XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS - !read(*,*) - ! Call TOTPT: Computes daily total photosynthesis and stomatal resistance - CALL TOTPT(HOD,DS,SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR,DVP,WNM,C3C4,LAI, & - TLAI,HT,LWIDTH,ROOTD,SD1,RSS,BLD,KN,KW,SLN,SLNT,SLNNT,SLNMIN,FB,& - DWSUP,CO2A,LS,EAJMAX,XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS, & - PPCAN,APCANS,APCANN,APCAN,PTCAN,ATCAN,PESOIL,AESOIL,DIFS,DIFSU, & - DIFSH,DAPAR,RCAN,ATRJC,ATRJS,FSR,FRSU,FRSH,ARSWSU,ARSWSH) - - RSLNB = (SLNBC-SLNB)/TCP !!Eq. 38, p. 35, DELT ok, g N m-2 leaf s-1 - - LNLV = MIN(LWLV,LWLVM)*LNCMIN + (LWLV-MIN(LWLV,LWLVM))*LNC !g m-2 s-1 - LCLV = LWLV*CFV !g m-2 s-1 - - ! Residual maintenance respiration g m-2 s-1 - RMRE = MAX(MIN(44./12.*0.218*(NTOT-WSH*LNCMIN-WRT*RNCMIN)/TCP & - ,APCAN-(1.E-5+RMUL)/TCP), 0.) !0.218 has the unit g C g-1 N d-1! - - TAVSS = TAIR + DIFS - - ! Call TUNIT: Computes cumulative thermal units (CTDU) - CALL TUNIT (1.*CROP, DS,TAIR,MAX(0.,DIFS),DAYL,TBD,TOD,TCD,TBDV,TODV,TCDV,TSEN,TDU,VDU) - - RMN = MAX(0., MIN(APCAN-1.E-5/TCP,RMUL/TCP) + MAX(MIN(44./12.*0.218* & - (1.001*NTOT-WSH*LNCMIN-WRT*RNCMIN)/TCP,APCAN-(1.E-5+RMUL)/TCP), 0.)) - - RM = MAX(0., MIN(APCAN-1.E-5/TCP,RMUL/TCP) + RMRE) !gC m-2 s-1 - - ! Daily and total C and N returns from crop to soil - ! ji: RLVDS eingefuehrt. Standardmaessig war der Wert von RLVDS auf 0.1 hard-gecoded - ! CLVD: Amount of carbon in dead leaves - ! CLVDS: Amount of carbon in dead leaves that have become litter in soil - - IF(DS.lt.0.25) THEN - LVDS = (CLVD-CLVDS)/TCP !gC m-2 s-1 - ELSE - LVDS = RLVDS*(CLVD-CLVDS)*(TAVSS-TBD)/(TOD-TBD)/TCP !gC m-2 s-1 - ENDIF - - CALL PHENO (1.*CROP,DS,SLP,DDLP,SPSP,EPSP,PSEN,MTDV,MTDR,TDU,CVDU,DVR) - - NFIXE = MAX(0., APCAN-(1.E-5+RM)/TCP)/CCFIX*Z1244 - - ! Daily carbon flow for seed filling - CALL BETAF(DVR,1.,PMES,LIMIT(1.,2.,DS)-1., FDS) - - LITC = LCRT + LVDS !gC m-2 s-1 - LITN = LNRT + LVDS/CFV *LNCMIN*PNLS !gN m-2 s-1 - - CALL BETAF(DVR,(1.+ESD)/2.,PMEH*(1.+ESD)/2.,MIN((1.+ESD)/2.,DS), FDH) - - NSUP = NSUP1 + NSUP2 !gN m-2 s-1 - - NFIX = INSW (LEGUME, 0., MIN(NFIXE, NFIXD)) !N fixation - - RNSUPP = (NSUP - NSUPP)/DT !RNSUPP g N m-2 s-2 - RX = 44./12.*(CCFIX*NFIX) - -!*** Current photo-assimilates (g CO2 m-2 s-1) for growth, and R/P ratio - ASSA = APCAN - RM - RX - -!*** Crop nitrogen demand and uptake (g N m-2 s-1) - - SHSA = Z1244 * YGV*MAX(1E-16,APCAN -RM -RX)/ MAX(0.1,CSH) - SHSAN = Z1244 * YGV*(APCANN-RMN-RX)/ MAX(0.1,CSH) - DERI = MAX(0.,(SHSAN - SHSA)/(0.001*MAX(0.01,NTOT)/MAX(0.1,CTOT))) - - RMUS = 0.06*0.05/0.454*YGV*ASSA - NDEMA = CRT * SHSA**2/NOTNUL(DERI) - NCR = INSW(SLNT-SLNMIN,0.,MIN(NUPTX,NDEMA))/(YGV*MAX(1E-16,APCANS-RM-RX)*Z1244) - NDEMAD = INSW(LNC-1.5*LNCI, MAX(NDEMA, NDEMD), 0.) - -!*** Nitrogen partitioning between shoots and roots - FNSH = 1./(1.+NCR*DERI/SHSA*CSH/MAX(1E-2,CRT)*NRT/MAX(1E-2,NSH)) - -!*** Carbon partitioning among organs and reserve pools - FCSH = 1./(1.+NCR*DERI/SHSA) - -!*** ji 10.02.14, winter dormancy -!*** ji 16.06.15, ecc/lcc switch - IF (CROP==1) THEN - NDEM = INSW(DS-DSCRIT,INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,.01*NDEMAD), 0.), & - INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,NDEMAD), 0.)) - ELSE - NDEM = INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,NDEMAD), 0.) - ENDIF - - DCSR = Z1244*(1.-FCSH)*ASSA - NUPTN = MIN(NSUP2, NSUP2/NOTNUL(NSUP)*MAX(0.,NDEM-NFIXR/TCP)) - RCSRT = Z1244*ASSA*(1.-FCSH)*(1.-FCRVR)*YGV - LCRT - RNFIXR = NFIX - MIN(NDEM,NFIXR/TCP) - RMLD = 0.06*(1.-FCSH)*ASSA - RNDEMP = (NDEM - NDEMP)/DT !as for RNSUPP, DELT ok - NUPTA = MIN(NSUP1, NSUP1/NOTNUL(NSUP)*MAX(0.,NDEM-NFIXR/TCP)) - - ! Carbon supply from current photo-assimilates for shoot & root growth - DCSS = Z1244* FCSH *ASSA - NUPT = MAX(0., NUPTA + NUPTN + MIN(NDEM, NFIXR/TCP)) - RMUA = 44./12.*0.17*NUPTA - - CALL SINKG(DS,1.,TSN*SEEDW*CFO,YGO,FDS,DCDSR,DCSS,DT,& - DCDSC,DCDS,FLWCS) - - ! Maintenance and total respiration (g CO2 m-2 d-1) - RMUN = 44./12.*2.05*NUPTN - - ! Daily carbon flow for structural stem growth - DCST = DCSS - FLWCS - FCSO = FLWCS/DCSS - - RRMUL = (RMUN+RMUA+RMUS+RMLD-RMUL)/DT - - GAP = MAX(0., DCDS-DCSS) - - !*** ji, 16.06.15, ecc/lcc switch - IF (CROP==1) THEN - IFSH = INSW(DCST-1E-11, 1., LIMIT(0.,1.,DCST/NOTNUL(DCDTP))) - ELSE - IFSH = LIMIT(0.,1.,DCST/NOTNUL(DCDTP)) - ENDIF - - CREMSI = MIN(0.94*CRVS, CRVS/NOTNUL(CRVS+CRVR)*GAP)/0.94 - CREMRI = MIN(0.94*CRVR, CRVR/NOTNUL(CRVS+CRVR)*GAP)/0.94 - CREMS = INSW(DCDS-DCSS, 0., CREMSI) - CREMR = INSW(DCDS-DCSS, 0., CREMRI) - - IF (CROP==1) THEN - RHT = MIN(HTMX-HT, FDH*HTMX*IFSH) - ELSE - RHT = MIN(HTMX-HT, FDH*HTMX*INSW(DCST-1E-4, 1., LIMIT(0.,1.,DCST/NOTNUL(DCDTP)))) - ENDIF - - IF (CROP==1) THEN - CALL SINKG(DS,.1,CDMHT*HTMX*CFV,YGV,FDH*IFSH,DCDTR,DCST,DT, & - DCDTC,DCDT,FLWCT) - ELSE - CALL SINKG(DS,.0,CDMHT*HTMX*CFV,YGV,FDH*IFSH,DCDTR,DCST,DT, & - DCDTC,DCDT,FLWCT) - ENDIF - - RCRVR = FCRVR*DCSR - CREMR - - IF (CROP==1) THEN - RDCDTP = (DCDTC-DCDTP)/DT - ELSE - DCDTP = DCDTC - RDCDTP = 0. - ENDIF - - !ji, the MIN function avoids that a situation occurs during which - !no assimilates are partioned to the leaves. FCSST max. 85% - IF (CROP==1) THEN - FCSST = MIN(.85,INSW(DS-(ESD+0.2), FLWCT/MAX(1E-16,DCSS), 0.)) - !FCSST = INSW(DS-(ESD+0.2), FLWCT/DCSS, 0.) original version - ELSE - FCSST = MIN(.85,INSW(DS-(ESD+0.2), FLWCT/MAX(1E-16,DCSS), 0.)) - ENDIF - - RCSO = Z1244*ASSA*FCSH*FCSO*YGO + 0.94*(CREMS+CREMR)*YGO - RWSO = RCSO / CFO - RDCDSR = MAX(0., (DCDSC-RCSO/YGO))-(FLWCS-MIN(DCDSC,DCSS)) - RWRT = RCSRT/CFV + RCRVR/0.444 - RCSST = Z1244*ASSA* FCSH * FCSST *YGV - - FCLV = REAAND(LAIN-LAIC,ESD-DS)*(1.-FCSO-FCSST) - - RRD = INSW(ROOTD-RDMX, MIN((RDMX-ROOTD)/TCP,(RWRT+LWRT)/(WRB+KR* & - (WRT+WRTD))), 0.) - - RCLV = Z1244*ASSA* FCSH * FCLV * YGV - LCLV - - FCRVS = 1. - FCLV - FCSO - FCSST - - RDCDTR = MAX(0., (DCDTC-RCSST/YGV))-(FLWCT-MIN(DCDTC,DCST)) - RCRVS = FCRVS*DCSS - CREMS - RWLV = RCLV / CFV - - RNRES = NUPT-(LNCMIN*(RCLV+LCLV)+RNCMIN*(RCSRT+LCRT)+STEMNC* & - RCSST)/CFV - - RG = 44./12.*((1.-YGV)/YGV*(RCLV+RCSST+RCSRT+LCLV+LCRT)+ & - (1.-YGO)/YGO* RCSO) - - RWST = RCSST/ CFV + RCRVS/0.444 - RNREOF = INSW (DS-1.0, RNRES, 0.) - RNREOE = INSW (DS-ESD, RNRES, 0.) - - RESTOT = RM+RX+RG + 44./12.*0.06*(CREMS+CREMR) - - CALL RNACC (FNSH,NUPT,RWST,STEMNC,LNCMIN,RNCMIN,LNC,RNC,NLV,NRT,WLV,WRT, & - DT,CB,CX,TM,DS,SEEDNC,RWSO,LNLV,LNRT, RNRT,RNST,RNLV,RTNLV,RNSO) - - RRP = RESTOT / APCAN - - CALL RLAIC(1.*CROP,DS,SLA0,RWLV,LAIC,KN,NLV,RNLV,SLNB,RSLNB, RLAI) - - ! Used for debugging -! if (nowdate(9:12).eq.'2500') then -! write(*,*) nowdate, ' 1 ', DVR -! write(*,*) nowdate, ' 2 ', TDU -! write(*,*) nowdate, ' 3 ', VDU -! write(*,*) nowdate, ' 4 ' -! write(*,*) nowdate, ' 5 ', RCLV -! write(*,*) nowdate, ' 6 ', LCLV -! write(*,*) nowdate, ' 7 ', RCSST -! write(*,*) nowdate, ' 8 ', RCSO -! write(*,*) nowdate, ' 9 ', RCSRT -! write(*,*) nowdate, ' 10 ',LCRT -! write(*,*) nowdate, ' 11 ',LVDS -! write(*,*) nowdate, ' 12 ',RNRT -! write(*,*) nowdate, ' 13 ',RNST -! write(*,*) nowdate, ' 14 ',RNLV -! write(*,*) nowdate, ' 15 ',RNSO -! write(*,*) nowdate, ' 16 ',RTNLV -! write(*,*) nowdate, ' 17 ',LNLV -! write(*,*) nowdate, ' 18 ',LNRT -! write(*,*) nowdate, ' 19 ',RCRVS -! write(*,*) nowdate, ' 20 ',RCRVR -! write(*,*) nowdate, ' 21 ',RNREOE -! write(*,*) nowdate, ' 22 ',RNREOF -! write(*,*) nowdate, ' 23 ',RDCDSR -! write(*,*) nowdate, ' 24 ',RDCDTR -! write(*,*) nowdate, ' 25 ',RSLNB -! write(*,*) nowdate, ' 26 ',RLAI -! write(*,*) nowdate, ' 27 ',RRMUL -! write(*,*) nowdate, ' 27 ',RRMUL -! write(*,*) nowdate, ' 28 ',RNDEMP -! write(*,*) nowdate, ' 29 ',RNSUPP -! write(*,*) nowdate, ' 30 ',NFIX -! write(*,*) nowdate, ' 31 ',RNFIXR -! write(*,*) nowdate, ' 32 ',RDCDTP -! write(*,*) nowdate, ' 33 ',RHT -! write(*,*) nowdate, ' 34 ',RRD -! write(*,*) nowdate, ' 35 ',APCAN -! write(*,*) nowdate, ' 36 ',RESTOT -! write(*,*) nowdate, ' 37 ',NUPT -! write(*,*) nowdate, ' 38 ',LITN -! read(*,*) -! endif - - ! Integration of ODEs - DS = MAX(0., INTGRL(DS, DVR, DT)) - CTDU = MAX(0., INTGRL(CTDU, TDU, DT)) - CVDU = MAX(0., INTGRL(CVDU, VDU, DT)) - CLV = MAX(0., INTGRL (CLV, RCLV, DT)) - CLVD = MAX(0., INTGRL (CLVD, LCLV, DT)) - CSST = MAX(0., INTGRL (CSST, RCSST, DT)) - CSO = MAX(0., INTGRL (CSO, RCSO, DT)) - CSRT = MAX(0., INTGRL (CSRT, RCSRT, DT)) - CRTD = MAX(0., INTGRL (CRTD, LCRT, DT)) - CLVDS = MAX(0.,INTGRL (CLVDS, LVDS, DT)) - NRT = MAX(0.,INTGRL (NRT, RNRT, DT)) - NST = MAX(0.,INTGRL (NST, RNST, DT)) - NLV = MAX(0.,INTGRL (NLV, RNLV, DT)) - NSO = MAX(0.,INTGRL (NSO, RNSO, DT)) - TNLV = MAX(0.,INTGRL (TNLV, RTNLV, DT)) - NLVD = MAX(0.,INTGRL (NLVD, LNLV, DT)) - NRTD = MAX(0.,INTGRL (NRTD, LNRT, DT)) - CRVS = MAX(0.,INTGRL (CRVS, RCRVS, DT)) - CRVR = MAX(0.,INTGRL (CRVR, RCRVR, DT)) - NREOE = MAX(0.,INTGRL(NREOE, RNREOE, DT)) - NREOF = MAX(0.,INTGRL(NREOF, RNREOF, DT)) - DCDSR = MAX(0.,INTGRL(DCDSR, RDCDSR, DT)) - DCDTR = MAX(0.,INTGRL(DCDTR, RDCDTR, DT)) - SLNB = MAX(0.,INTGRL(SLNB, RSLNB, DT)) - LAIC = MAX(0.,INTGRL(LAIC, RLAI, DT)) - RMUL = MAX(0.,INTGRL(RMUL, RRMUL, DT)) - NDEMP = MAX(0.,INTGRL(NDEMP, RNDEMP, DT)) - NSUPP = MAX(0.,INTGRL(NSUPP, RNSUPP, DT)) - NFIXT = MAX(0.,INTGRL(NFIXT, NFIX, DT)) - NFIXR = MAX(0.,INTGRL(NFIXR, RNFIXR, DT)) - DCDTP = MAX(0.,INTGRL(DCDTP, RDCDTP, DT)) - HT = MAX(0.,INTGRL(HT, RHT, DT)) - ROOTD = MAX(1.,INTGRL(ROOTD, RRD, DT)) - TPCAN = MAX(0.,INTGRL(TPCAN, APCAN, DT)) - TRESP = MAX(0.,INTGRL(TRESP, RESTOT, DT)) - TNUPT = MAX(0.,INTGRL(TNUPT, NUPT, DT)) - LITNT = MAX(1e-3,INTGRL(LITNT, LITN, DT)) - - ! Write updated Gecros variables into STATE_GECROS array - STATE_GECROS(1) = DS - STATE_GECROS(2) = CTDU - STATE_GECROS(3) = CVDU - STATE_GECROS(4) = CLV - STATE_GECROS(5) = CLVD - STATE_GECROS(6) = CSST - STATE_GECROS(7) = CSO - STATE_GECROS(8) = CSRT - STATE_GECROS(9) = CRTD - STATE_GECROS(10) = CLVDS - STATE_GECROS(11) = NRT - STATE_GECROS(12) = NST - STATE_GECROS(13) = NLV - STATE_GECROS(14) = NSO - STATE_GECROS(15) = TNLV - STATE_GECROS(16) = NLVD - STATE_GECROS(17) = NRTD - STATE_GECROS(18) = CRVS - STATE_GECROS(19) = CRVR - STATE_GECROS(20) = NREOE - STATE_GECROS(21) = NREOF - STATE_GECROS(22) = DCDSR - STATE_GECROS(23) = DCDTR - STATE_GECROS(24) = SLNB - STATE_GECROS(25) = LAIC - STATE_GECROS(26) = RMUL - STATE_GECROS(27) = NDEMP - STATE_GECROS(28) = NSUPP - STATE_GECROS(29) = NFIXT - STATE_GECROS(30) = NFIXR - STATE_GECROS(31) = DCDTP - STATE_GECROS(32) = HT - STATE_GECROS(33) = ROOTD - STATE_GECROS(34) = TPCAN - STATE_GECROS(35) = TRESP - STATE_GECROS(36) = TNUPT - STATE_GECROS(37) = LITNT - STATE_GECROS(45) = WSO - STATE_GECROS(46) = WSTRAW - STATE_GECROS(47) = GrainNC - STATE_GECROS(48) = StrawNC - STATE_GECROS(49) = LAI - STATE_GECROS(50) = TLAI - - ! Used for debugging -! if(debugging) then -! if (nowdate(9:12).eq.'1200') then -! write(*,*) nowdate, RSD -! write(*,*) "DS: ", DS -! write(*,*) "CTDU: ",CTDU -! write(*,*) "CVDU: ",CVDU -! write(*,*) "Carbon in living leaves CLV: ",CLV -! write(*,*) "Carbon in dead leaves CLVD: ",CLVD -! write(*,*) "Carbon in structural stems CSST: ",CSST -! write(*,*) "Carbon in storage organs CSO: ",CSO -! write(*,*) "Carbon in structural roots CSRT: ",CSRT -! write(*,*) "Carbon in living roots CRT: ",CRT -! write(*,*) "Carbon in dead roots CRTD: ",CRTD -! write(*,*) "Rooting depth ROOTD: ",ROOTD -! write(*,*) "Carbon litter on soil CLVDS: ",CLVDS -! write(*,*) "Carbon in stem reserves CRVS: ",CRVS -! write(*,*) "Carbon in root reserves CRVR: ",CRVR -! write(*,*) "Total leaf area index TLAI: ",TLAI -! write(*,*) "Dead leaf area index DLAI: ",DLAI -! write(*,*) "Living leaf area index GLAI: ",LAI, LAIC, LAIN -! write(*,*) "Specific leaf area index SLA: ",SLA -! write(*,*) "Carbon balance check CCHK: ",CCHK -! write(*,*) "Nitrogen balance check NCHK: ",NCHK -! read(*,*) -! endif -! endif - -! ---------------------------------------------------------------------- -END SUBROUTINE gecros -! ---------------------------------------------------------------------- - -!******************** SUBROUTINES FOR CROP SIMULATION ******************* -!*----------------------------------------------------------------------* -!* SUBROUTINE TUNIT * -!* Purpose: This subroutine calculates the daily amount of thermal day * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DS R4 Development stage - I * -!* TMAX R4 Daily maximum temperature oC I * -!* TMIN R4 Daily minimum temperature oC I * -!* DIF R4 Daytime plant-air temperature differential oC I * -!* DAYL R4 Astronomic daylength (base = 0 degrees) h I * -!* TBD R4 Base temperature for phenology oC I * -!* TOD R4 Optimum temperature for phenology oC I * -!* TCD R4 Ceiling temperature for phenology oC I * -!* TSEN R4 Curvature for temperature response - I * -!* TDU R4 Rate of thermal days increase d s-1 O * -!* VDU R4 Rate of vernalization days increase d s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE TUNIT(CROP,DS,TAIR,DIF,DAYL,TBD,TOD,TCD,TBDV,TODV,TCDV,TSEN,TDU,VDU) - IMPLICIT REAL (A-Z) - INTEGER I - -!*---assuming development rate at supra-optimum temperatures during -!* the reproductive phase equals that at the optimum temperature - IF (DS.GT.1.) THEN - TAIR = MIN (TAIR,TOD) - ELSE - TAIR = TAIR - ENDIF - -!*---vernalization response (Lenz 2007, p. 26) -!*---Ingwersen 29.6.2010 - IF (TAIR.LT.TBDV .OR. TAIR.GT.TCDV) THEN - TUV = 0. - ELSE - TUV = (((TCDV-TAIR)/(TCDV-TODV))*((TAIR-TBDV)/(TODV-TBDV))**((TODV-TBDV)/(TCDV-TODV)))**TSEN - ENDIF - -!*---instantaneous thermal unit based on bell-shaped temperature response - IF (TAIR.LT.TBD .OR. TAIR.GT.TCD) THEN - TU = 0. - ELSE - TU = (((TCD-TAIR)/(TCD-TOD))*((TAIR-TBD)/(TOD-TBD))**((TOD-TBD)/(TCD-TOD)))**TSEN - ENDIF - -!*---daily thermal unit -!*** ji, 16.6.15, ecc/lcc switch - IF (CROP==1) THEN - TDU = TU/TCP - ELSE - TDU = INSW(DS-2.,TU/TCP,0.) - ENDIF - VDU = TUV/TCP - - RETURN - END SUBROUTINE TUNIT - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PHENO * -!* Purpose: This subroutine calculates phenological development rate. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DS R4 Development stage - I * -!* SLP R4 Crop type(1. for short-day,-1. for long-day) - I * -!* DDLP R4 Daylength for photoperiodism h I * -!* SPSP R4 DS for start of photoperiod-sensitive phase - I * -!* EPSP R4 DS for end of photoperiod-sensitive phase - I * -!* PSEN R4 Photoperiod sensitivity (+ for SD, - for LD) h-1 I * -!* MTDV R4 Minimum thermal days for vegetative phase d I * -!* MTDR R4 Minimum thermal days for reproductive phase d I * -!* TDU R4 Thermal unit - I * -!* CVDU R4 Cumulative vernalization unit (ji, 29.6) - I * -!* DVR R4 Development rate s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE PHENO (CROP,DS,SLP,DDLP,SPSP,EPSP,PSEN,MTDV,MTDR,TDU,CVDU,DVR) - IMPLICIT REAL (A-Z) - -!*---determining if it is for short-day or long-day crop - IF (SLP.LT.0.) THEN - MOP = 18. !minimum optimum photoperiod for long-day crop - DLP = MIN(MOP,DDLP) - ELSE - MOP = 11. !maximum optimum photoperiod for short-day crop - DLP = MAX(MOP,DDLP) - ENDIF - -!*---effect of photoperiod on development rate - IF (DS.LT.SPSP .OR. DS.GT.EPSP) THEN - EFP = 1. - ELSE - EFP = MAX(0., 1.-PSEN*(DLP-MOP)) - ENDIF - -!*---effect of vernalization (ji, 21.6.10) -!*** ji, 16.6.15, ecc/lcc switch - IF (CROP==1) THEN - EFV = CVDU**5./(22.5**5. + CVDU**5.) - ELSE - EFV = 1.0 - ENDIF - -!*---development rate of vegetative and reproductive phases -!*---extended for vernalization according to Lenz (2007); ji: 21.6.2010 - IF (DS.LE.0.4) THEN - DVR = 1./MTDV*TDU*EFP*EFV - ENDIF - - IF (DS.GT.0.4 .AND. DS.LE.1.0) THEN - DVR = 1./MTDV*TDU*EFP - ENDIF - - IF (DS.GT.1.0) THEN - DVR = 1./MTDR*TDU - ENDIF - - RETURN - END SUBROUTINE PHENO - - -!*----------------------------------------------------------------------* -!* SUBROUTINE RNACC * -!* Purpose: This subroutine calculates rate of N accumulation in organs* -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* FNSH R4 Fraction of new N partitioned to shoot - I * -!* NUPT R4 Nitrogen uptake at a time step gN/m2/s I * -!* RWST R4 Rate of stem weight g/m2/s I * -!* STEMNC R4 Nitrogen concentration in stem gN/g I * -!* LNCMIN R4 Minimum N concentration in leaf gN/g I * -!* RNCMIN R4 Minimum N concentration in root gN/g I * -!* LNC R4 Nitrogen concentration in leaf gN/g I * -!* RNC R4 Nitrogen concentration in root gN/g I * -!* NLV R4 Canopy (green)leaf N content gN/m2 I * -!* NRT R4 (living)root N content gN/m2 I * -!* WLV R4 Canopy (green)leaf weight g/m2 I * -!* WRT R4 (living)Root weight g/m2 I * -!* DELT R4 Time step of simulation s I * -!* CB R4 Factor for initial N concent. of seed-fill - I * -!* CX R4 Factor for final N concent. of seed-fill - I * -!* TM R4 DS when transition from CB to CX is fastest - I * -!* DS R4 Development stage - I * -!* SEEDNC R4 Standard seed N concentration gN/g I * -!* RWSO R4 growth rate of seed g/m2/s I * -!* LNLV R4 Loss rate of NLV due to senescence gN/m2/s I * -!* LNRT R4 Loss rate of NRT due to senescence gN/m2/s I * -!* RNRT R4 rate of N accumulation in root gN/m2/s O * -!* RNST R4 rate of N accumulation in stem gN/m2/s O * -!* RNLV R4 rate of N accumulation in leaf gN/m2/s O * -!* RTNLV R4 Positive value of RNLV gN/m2/s O * -!* RNSO R4 rate of N accumulation in seed(storage organ)gN/m2/s O * -!*----------------------------------------------------------------------* - SUBROUTINE RNACC (FNSH,NUPT,RWST,STEMNC,LNCMIN,RNCMIN,LNC,RNC, & - NLV,NRT,WLV,WRT,DELT,CB,CX,TM,DS,SEEDNC, & - RWSO,LNLV,LNRT, RNRT,RNST,RNLV,RTNLV,RNSO) - IMPLICIT REAL (A-Z) - -!*---amount of N partitioned to shoot - NSHN = FNSH * NUPT - -!*---leaf N (NLVA) or root N (NRTA) available for remobilization within a day - NLVA = INSW(LNCMIN-LNC, NLV-WLV*LNCMIN, 0.)/TCP - NRTA = INSW(RNCMIN-RNC, NRT-WRT*RNCMIN, 0.)/TCP - - NTA = NLVA + NRTA - -!*---rate of N accumulation in stem - RNST = RWST * INSW(-NTA,STEMNC,0.) - -!*---expected N dynamics during seed(storage organ) filling - CDS = CB+(CX-CB)*(4.-TM-DS)/(2.-TM)*(DS-1.)**(1./(2.-TM)) - ENSNC = LIMIT(CB,CX,CDS) * SEEDNC - -!*---rate of N accumulation in seed - NGS = NSHN - RNST - ENSNC*RWSO - NONC = MAX(0.,INSW(NTA+NGS,(NTA+NSHN-RNST)/NOTNUL(RWSO),ENSNC)) - RNSO = RWSO*NONC - -!*---rate of N accumulation in leaf! - - NLVN = INSW(NTA+NGS,-NLVA-LNLV,-NLVA/NOTNUL(NTA)*(-NGS)-LNLV) - GNLV = INSW(NGS, NLVN, NSHN-RNST-RNSO-LNLV) - RNLV = MAX (-NLV+1.E-7, GNLV) - RTNLV = MAX(0., RNLV) - -!*---rate of N accumulation in root - NRTN = INSW(NTA+NGS, NUPT-NSHN-NRTA-LNRT, & - NUPT-NSHN-NRTA/NOTNUL(NTA)*(-NGS)-LNRT) - GNRT = INSW(NGS, NRTN, NUPT-NSHN-LNRT) - RNRT = MAX (-NRT+5.E-8, GNRT) - - RETURN - END SUBROUTINE RNACC - - -!*----------------------------------------------------------------------* -!* SUBROUTINE RLAIC * -!* Purpose: This subroutine calculates the daily increase of leaf * -! area index (m2 leaf/m2 ground/day). * -! * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DS R4 Development stage - I * -!* SLA0 R4 Specific leaf area constant m2 g-1 I * -!* RWLV R4 Rate of increment in leaf weight g m-2 s-1 I * -!* LAI R4 Leaf area index m2 m-2 I * -!* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I * -!* NLV R4 Total leaf nitrogen content in a canopy g m-2 I * -!* RNLV R4 Rate of increment in NLV g m-2 s-1 I * -!* SLNB R4 Nitrogen content of bottom leaves g m-2 I * -!* RSLNB R4 Rate of increment in SLNB g m-2 s-1 I * -!* RLAI R4 Rate of increment in leaf area index m2 m-2s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE RLAIC(CROP,DS,SLA0,RWLV,LAI,KN,NLV,RNLV,SLNB,RSLNB, RLAI) - IMPLICIT REAL (A-Z) - - SLNB = MAX(1E-2, SLNB) - !*---rate of LAI driven by carbon supply - RLAI = INSW(RWLV, MAX(-LAI+1.E-5,SLA0*RWLV), SLA0*RWLV) - - !*---rate of LAI driven by nitrogen during juvenile phase - !*** ji, 16.6.15, ecc/lcc switch - IF ((CROP==2) .AND. (LAI.LT.1.5) .AND. (DS.LT.0.75)) THEN - RLAI = MAX(0.,(SLNB*RNLV-NLV*RSLNB)/SLNB/(SLNB+KN*NLV)) - ENDIF - - - RETURN - END SUBROUTINE RLAIC - - -!*----------------------------------------------------------------------* -!* SUBROUTINE BETAF * -!* Purpose: This subroutine calculates the dynamics of expected growth * -!* of sinks, based on the beta sigmoid growth equation * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DVR R4 Development rate s-1 I * -!* TE R4 Stage at which sink growth stops - I * -!* TX R4 Stage at which sink growth rate is maximal - I * -!* TI R4 Stage of a day considered - I * -!* FD R4 Relative expected growth of a sink at a day s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE BETAF(DVRX,TE,TX,TI, FD) - - REAL, INTENT(IN) :: DVRX, TE, TX, TI - REAL, INTENT(OUT) :: FD - - FD = DVRX*(2.*TE-TX)*(TE-TI)/TE/(TE-TX)**2*(TI/TE)**(TX/(TE-TX)) - !FD = DVRX*(TE-TI)/(TE-TX)*(TI/TX)**(TX/(TE-TX)) !Eq. a 1a Yin et al 2003 - - END SUBROUTINE BETAF - - -!*----------------------------------------------------------------------* -!* SUBROUTINE SINKG * -!* Purpose: This subroutine calculates carbon demand for sink growth. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DS R4 Development stage - I * -!* SSG R4 Stage at which sink growth starts - I * -!* TOTC R4 Total carbon in a sink at end of its growth g C/m2 I * -!* YG R4 Growth efficiency g C/g C I * -!* FD R4 Relative expected growth of a sink at a day s-1 I * -!* DCDR R4 Shortfall of C demand in previous days g C/m2 I * -!* DCS R4 Daily C supply for sink growth g C/m2/s I * -!* DELT R4 Time step of integration s I * -!* DCDC R4 C demand of the current day g C/m2/s O * -!* DCD R4 Daily C demand for sink growth g C/m2/s O * -!* FLWC R4 Flow of current assimilated C to sink g C/m2/s O * -!*----------------------------------------------------------------------* - SUBROUTINE SINKG(DS,SSG,TOTC,YG,FD,DCDR,DCS,DELT,DCDC,DCD,FLWC) - IMPLICIT REAL (A-Z) - -!*---expected demand for C of the current time step - DCDC = INSW (DS-SSG, 0., TOTC/YG*FD) - -!*---total demand for C at the time step considered - DCD = DCDC + MAX(0.,DCDR)/DELT - -!*---flow of current assimilated carbon to sink - FLWC = MIN(DCD, DCS) - - RETURN - END SUBROUTINE SINKG - - -!*----------------------------------------------------------------------* -!* SUBROUTINE ASTRO (from the SUCROS model) * -!* Purpose: This subroutine calculates astronomic daylength, * -!* diurnal radiation characteristics such as the daily * -!* integral of sine of solar elevation and solar constant. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* DOY R4 Daynumber (Jan 1st = 1) - I * -!* LAT R4 Latitude of the site degree I * -!* INSP R4 Inclination of sun angle for computing DDLP degree I * -!* SC R4 Solar constant J m-2 s-1 O * -!* SINLD R4 Seasonal offset of sine of solar height - O * -!* COSLD R4 Amplitude of sine of solar height - O * -!* DAYL R4 Astronomic daylength (base = 0 degrees) h O * -!* DDLP R4 Photoperiodic daylength h O * -!* DSINBE R4 Daily total of effective solar height s d-1 O * -!* * -!* FATAL ERROR CHECKS (execution terminated, message) * -!* condition: LAT > 67, LAT < -67 * -!* * -!* FILE usage : none * -!*----------------------------------------------------------------------* - SUBROUTINE ASTRO (DOY,LAT,INSP,SC,SINLD,COSLD,DAYL,DDLP,DSINBE) - IMPLICIT REAL (A-Z) - -!*---PI and conversion factor from degrees to radians - PI = 3.141592654 - RAD = PI/180. -!*---check on input range of parameters - IF (LAT.GT.67.) STOP 'ERROR IN ASTRO: LAT> 67' - IF (LAT.LT.-67.) STOP 'ERROR IN ASTRO: LAT>-67' - -!*---declination of the sun as function of daynumber (DOY) - DEC = -ASIN (SIN (23.45*RAD)*COS (2.*PI*(DOY+10.)/365.)) - -!*---SINLD, COSLD and AOB are intermediate variables - SINLD = SIN (RAD*LAT)*SIN (DEC) - COSLD = COS (RAD*LAT)*COS (DEC) - AOB = SINLD/COSLD - -!*---daylength (DAYL) - DAYL = 12.0*(1.+2.*ASIN (AOB)/PI) - DDLP = 12.0*(1.+2.*ASIN((-SIN(INSP*RAD)+SINLD)/COSLD)/PI) - - DSINB = 3600.*(DAYL*SINLD+24.*COSLD*SQRT (1.-AOB*AOB)/PI) - DSINBE = 3600.*(DAYL*(SINLD+0.4*(SINLD*SINLD+COSLD*COSLD*0.5))+ & - 12.0*COSLD*(2.0+3.0*0.4*SINLD)*SQRT (1.-AOB*AOB)/PI) - -!*---solar constant (SC) - SC = 1367.*(1.+0.033*COS(2.*PI*(DOY-10.)/365.)) - - RETURN - END SUBROUTINE ASTRO - - -!*----------------------------------------------------------------------* -!* SUBROUTINE TOTPT * -!* Purpose: This subroutine calculates daily total gross photosynthesis* -!* and transpiration by performing a Gaussian integration * -!* over time. At five different times of the day, temperature * -!* and radiation are computed to determine assimilation * -!* and transpiration whereafter integration takes place. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* HOD R4 Hour of the day h I * -!* SC R4 Solar constant J m-2 s-1 I * -!* SINLD R4 Seasonal offset of sine of solar height - I * -!* COSLD R4 Amplitude of sine of solar height - I * -!* DAYL R4 Astronomic daylength (base = 0 degrees) h I * -!* DSINBE R4 Daily total of effective solar height s d-1 I * -!* RSD R4 Global radiation J m-2 s-1 I * -!* TAIR R4 Air temperature oC I * -!* DVP R4 Vapour pressure kPa I * -!* WNM R4 daily average wind speed (>=0.1 m/s) m s-1 I * -!* C3C4 R4 Crop type (=1 for C3, -1 for C4 crops) - I * -!* LAI R4 (green)Leaf area index m2 m-2 I * -!* TLAI R4 Total Leaf area index m2 m-2 I * -!* HT R4 Plant height m I * -!* LWIDTH R4 Leaf width m I * -!* RD R4 Rooting depth cm I * -!* SD1 R4 Depth of evaporative upper soil layer cm I * -!* RSS R4 Soil resistance,equivalent to leaf stomata s m-1 I * -!* BLD R4 Leaf angle from horizontal degree I * -!* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I * -!* KW R4 Windspeed extinction coefficient in canopy m2 m-2 I * -!* SLN R4 Average leaf nitrogen content in canopy g m-2 I * -!* SLNT R4 Top-leaf nitrogen content g m-2 I * -!* SLNN R4 Value of SLNT with small plant-N increment g m-2 I * -!* SLNMIN R4 Minimum or base SLNT for photosynthesis g m-2 I * -!* DWSUP R4 Daily water supply for evapotranspiration mm s-1 I * -!* CO2A R4 Ambient CO2 concentration ml m-3 I * -!* LS R4 Lodging severity - I * -!* EAJMAX R4 Energy of activation for Jmax J mol-1 I * -!* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I * -!* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I * -!* THETA R4 Convexity for light response of e-transport - I * -!* WCUL R4 Water content of the upper soil layer m3 m-3 I * -!* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I * -!* PPCAN R4 Potential canopy CO2 assimilation g m-2 s-1 O * -!* APCANS R4 Actual standing-canopy CO2 assimilation g m-2 s-1 O * -!* APCANN R4 APCANS with small plant-N increment g m-2 s-1 O * -!* APCAN R4 Actual canopy CO2 assimilation g m-2 s-1 O * -!* PTCAN R4 Potential canopy transpiration mm s-1 O * -!* ATCAN R4 Actual canopy transpiration mm s-1 O * -!* PESOIL R4 Potential soil evaporation mm s-1 O * -!* AESOIL R4 Actual soil evaporation mm s-1 O * -!* DIFS R4 Soil-air temp. difference oC O * -!* DIFSU R4 Sunlit leaf-air temp. diff. oC O * -!* DIFSH R4 Shaded leaf-air temp. diff. oC O * -!* DAPAR R4 PAR absorbed by crop canopy J m-2 s-1 O * -!* RCAN R4 Canopy resistance s/m 0 * -!*----------------------------------------------------------------------* -SUBROUTINE TOTPT(HOD,DS, SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR, DVP, & - WNM,C3C4,LAI,TLAI,HT,LWIDTH,RD,SD1,RSS,BLD,KN,KW, & - SLN,SLNT,SLNN,SLNMIN,FB,DWSUP,CO2A,LS,EAJMAX, & - XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS,PPCAN,APCANS, & - APCANN,APCAN,PTCAN,ATCAN,PESOIL,AESOIL,DIFS,DIFSU,DIFSH,DAPAR, & - RCAN, ATRJC, ATRJS, FSR, FRSU, FRSH, ARSWSU, ARSWSH) - - REAL, INTENT(IN) :: HOD,DS,SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR, DVP, & - WNM,C3C4,LAI,TLAI,HT,LWIDTH,RD,SD1,RSS,BLD,KN,KW, & - SLN,SLNT,SLNN,SLNMIN,FB,DWSUP,CO2A,LS,EAJMAX, & - XVN,XJN,THETA,WCUL,FVPD,RT,RTS !RB - - REAL, INTENT(INOUT) :: PPCAN,APCANS,APCANN,APCAN,PTCAN,ATCAN,PESOIL, & - AESOIL,DIFS,DIFSU,DIFSH,DAPAR,RCAN,ATRJC,ATRJS,FSR,FRSU,FRSH,ARSWSU,ARSWSH - REAL :: ACO2I,ADIFS,ADIFSH,ADIFSU,ANIRSH,ANIRSU,ANRAD,APAR,APARSH - REAL :: APARSU,ASVP,ATMTR,ATRJSH,ATRJSU,ATSH,ATSU,AV_Albedo - REAL :: AV_RSWSH,AV_RSWSU,BL,CUMRSD,DATRJC,DATRJS,DNRAD,FRDF - REAL :: GBHC,GBHLF,GBHSH,GBHSU,IAE,IAP,IAPL,IAPN,IAPNN,IAPS,IAT,IPE - REAL :: IPH,IPHSOIL,IPP,IPPL,IPT,IRDL,KB,KBPNIR,KBPPAR,KDPNIR,KDPPAR - REAL :: NIR,NIRDF,NIRDR,NPSH,NPSHN,NPSU,NPSUN,NRADS,NRADSH,PANSH,PANSU - REAL :: PAR,NRADSU,PARDF,PARDR,PASSH,PASSU,PCBNIR,PCBPAR,PCDNIR,PCDPAR - REAL :: PHCAN,PHSH,PHSOIL,PHSU,PI,PLFSH,PLFSU,PSNIR,PSPAR,PT1,PTSH,PTSU - REAL :: RBHS,RBHSH,RBHSU,RBWS,RBWSH,RBWSU,RSWSH,RSWSU,SCPNIR - REAL :: SCPPAR,SINB,SLOPSH,SLOPSU,WND,WSUP,WSUP1,Albedo,RB, DSsw - REAL :: TLEAFSH, TLEAFSU - - PI = 3.141592654 - -!*---output-variables set to 0. and five different times of a day(HOUR) - PPCAN = 0. - APCANS = 0. - APCANN = 0. - APCAN = 0. - PTCAN = 0. - PHCAN = 0. - ATCAN = 0. - PESOIL = 0. - AESOIL = 0. - PHSOIL = 0. - DIFS = 0. - DIFSU = 0. - DIFSH = 0. - DAPAR = 0. - DNRAD = 0. - DATRJC = 0. - DATRJS = 0. - AV_RSWSU = 0. - AV_RSWSH = 0. - AV_Albedo = 0. - CUMRSD = 0. - -!*---sine of solar elevation - SINB = MAX (.01, SINLD+COSLD*COS(2.*PI*(HOD-12.)/24.)) - -!*---daytime course of water supply - WSUP = DWSUP - WSUP1 = WSUP*SD1/RD - -!*---daytime course of wind speed - WND = WNM !m s-1 - -!*---total incoming PAR and NIR - PAR = (1.-FB)*0.5*RSD !J m-2 s-1 - NIR = (1.-FB)*0.5*RSD !J m-2 s-1 - -!*---diffuse light fraction (FRDF) from atmospheric transmission (ATMTR) - ATMTR = PAR/(0.5*SC*SINB) ! unitless, fraction - - IF (ATMTR.LE.0.22) THEN - FRDF = 1. - ELSE IF (ATMTR.GT.0.22 .AND. ATMTR.LE.0.35) THEN - FRDF = 1.-6.4*(ATMTR-0.22)**2 - ELSE - FRDF = 1.47-1.66*ATMTR - ENDIF - - FRDF = MAX (FRDF, 0.15+0.85*(1.-EXP (-0.1/SINB))) - -!*---incoming diffuse PAR (PARDF) and direct PAR (PARDR) - PARDF = PAR * FRDF !J m-2 s-1 - PARDR = PAR - PARDF !J m-2 s-1 - -!*---incoming diffuse NIR (NIRDF) and direct NIR (NIRDR) - NIRDF = NIR * FRDF !J m-2 s-1 - NIRDR = NIR - NIRDF !J m-2 s-1 - -!*---extinction and reflection coefficients - BL = BLD*PI/180. !leaf angle, conversion to radians - CALL KBEAM (SINB,BL,KB) - - SCPPAR = 0.2 !leaf scattering coefficient for PAR - SCPNIR = 0.8 !leaf scattering coefficient for NIR - CALL KDIFF (TLAI,BL,SCPPAR, KDPPAR) - CALL KDIFF (TLAI,BL,SCPNIR, KDPNIR) - - CALL REFL (SCPPAR,KB, KBPPAR,PCBPAR) - CALL REFL (SCPNIR,KB, KBPNIR,PCBNIR) - - PCDPAR = 0.057 !canopy diffuse PAR reflection coefficient - PCDNIR = 0.389 !canopy diffuse NIR reflection coefficient - -!*---fraction of sunlit and shaded components in canopy -!* ji 4.7.11 LAI -> TLAI - FRSU = 1./KB/TLAI*(1.-EXP(-KB*TLAI)) - FRSH = 1.-FRSU - -!*---leaf boundary layer conductance for canopy, sunlit and shaded leaves !m s-1 - GBHLF = 0.01*SQRT(WND/LWIDTH) - RBHSU = RB/(FRSU*LAI) !boundary layer resistance to heat,sunlit part !s m-1 - RBWSU = RB/(FRSU*LAI) !boundary layer resistance to H2O, sunlit part - RBHSH = RB/(FRSH*LAI) !boundary layer resistance to heat,shaded part - RBWSH = RB/(FRSH*LAI) !boundary layer resistance to H2O, shaded part - RCAN = WNM*10. - -!*---boundary layer resistance for soil !s m-1 - RBHS = 172.*SQRT(0.05/MAX(0.1,WND*EXP(-KW*TLAI))) - RBWS = 0.93*RBHS - -!*---photosynthetically active nitrogen for sunlit and shaded leaves - CALL PAN (SLNT,SLNMIN,LAI,KN,KB, NPSU,NPSH) - CALL PAN (SLNN,SLNMIN,LAI,KN,KB, NPSUN,NPSHN) - -!*---absorbed PAR and NIR by sunlit leaves and shaded leaves -!* ji 4.7.11 LAI->TLAI -!* Ansonsten fuehrt es dazu, dass nach dem Abreifen die Albedo auf >0.7 ansteigt und vom Bestand viel -!* zu wenig Energie mehr absorbiert wird - CALL LIGAB (SCPPAR,KB,KBPPAR,KDPPAR,PCBPAR,PCDPAR,PARDR,PARDF,TLAI,APARSU,APARSH) - CALL LIGAB (SCPNIR,KB,KBPNIR,KDPNIR,PCBNIR,PCDNIR,NIRDR,NIRDF,TLAI,ANIRSU,ANIRSH) - APAR = APARSU+APARSH !J m-2 s-1 - -!*---absorbed total radiation (PAR+NIR) by sunlit and shaded leaves - ATRJSU = APARSU+ANIRSU !J m-2 s-1 - ATRJSH = APARSH+ANIRSH !J m-2 s-1 - ATRJC = ATRJSH + ATRJSU !J m-2 s-1 - -!*---absorbed total radiation (PAR+NIR) by soil - PSPAR = 0.1 !soil PAR reflection - PSNIR = INSW(WCUL-0.5, 0.52-0.68*WCUL, 0.18) !soil NIR reflection - ATRJS=(1.-PSPAR)*(PARDR*EXP(-KBPPAR*TLAI)+PARDF*EXP(-KDPPAR*TLAI)) & - +(1.-PSNIR)*(NIRDR*EXP(-KBPNIR*TLAI)+NIRDF*EXP(-KDPNIR*TLAI)) - -!* ji 4.7.11 Berechnung der Albedo eingefuegt. Nachts wird die Albedo auf einen Durchschnittswert gesetzt (.2) - FSR = RSD-ATRJC-ATRJS - - if (RSD.gt.0) then - Albedo = (RSD-ATRJC-ATRJS)/RSD - else - Albedo = .2 - endif - -!*---instantaneous potential photosynthesis and transpiration - CALL PPHTR(FRSU,TAIR,DVP,CO2A,C3C4,FVPD,APARSU,NPSU,RBWSU,RBHSU, & - RT*FRSU,ATRJSU,ATMTR,EAJMAX,XVN,XJN,THETA,PLFSU, & - PTSU,PHSU,RSWSU,NRADSU,SLOPSU) - - CALL PPHTR(FRSH,TAIR,DVP,CO2A,C3C4,FVPD,APARSH,NPSH,RBWSH,RBHSH, & - RT*FRSH,ATRJSH,ATMTR,EAJMAX,XVN,XJN,THETA,PLFSH, & - PTSH,PHSH,RSWSH,NRADSH,SLOPSH) - - IPP = PLFSU + PLFSH !gCO2/m2/s - IPT = PTSU + PTSH !mm s-1 - IPH = PHSU + PHSH !J m-2 s-1 - ANRAD = NRADSU + NRADSH !J m-2 s-1 - -!*--- PT1: Potential transpiration using water from the upper evaporative layer (mm s-1) -!*--- SD1: thickness of upper evaporative layer (cm); default value: 5 cm -!*--- RD: Rooting depth (cm) - PT1 = IPT * SD1/RD - -!*---instantaneous potential soil evaporation - CALL PEVAP (TAIR,DVP,RSS,RTS,RBWS,RBHS,ATRJS,ATMTR, & - PT1,WSUP1,IPE,IPHSOIL,NRADS) - -!*---instantaneous actual soil evaporation, actual canopy -!* transpiration and photosynthesis - IAE = MIN (IPE,IPE/(PT1+IPE)*WSUP1) !actual soil evaporation mm s-1 - IAT = MIN (IPT,PT1/(PT1+IPE)*WSUP1+WSUP-WSUP1) !actual transpiration mm s-1 - ATSU = PTSU/IPT*IAT !actual transpiration of sunlit leaves mm s-1 - ATSH = PTSH/IPT*IAT !actual transpiration of shaded leaves mm s-1 - - CALL DIFLA (NRADS,IAE,RBHS,RTS, ADIFS) - - CALL APHTR (TAIR,APARSU,DVP,CO2A,C3C4,FVPD,NRADSU,ATSU,PTSU, & - RT*FRSU,RBHSU,RBWSU,RSWSU,SLOPSU,NPSU,NPSUN, & - EAJMAX,XVN,XJN,THETA,PASSU,PANSU,ADIFSU,ARSWSU) - - CALL APHTR (TAIR,APARSH,DVP,CO2A,C3C4,FVPD,NRADSH,ATSH,PTSH, & - RT*FRSH,RBHSH,RBWSH,RSWSH,SLOPSH,NPSH,NPSHN, & - EAJMAX,XVN,XJN,THETA,PASSH,PANSH,ADIFSH,ARSWSH) - - IAPS = PASSU + PASSH !actual canopy photosynthesis gCO2 m-2 s-1 - IAPN = PANSU + PANSH !actual canopy photosynthesis with a small N increment gCO2 m-2 s-1 - IAP = IAPS - IAPNN = IAPN - -!*---integration of assimilation and transpiration to a daily total - PPCAN = IPP - APCANS = IAPS - APCANN = IAPNN - APCAN = IAP - PTCAN = IPT - PHCAN = IPH - ATCAN = IAT - PESOIL = IPE - AESOIL = IAE - PHSOIL = IPHSOIL - DIFS = ADIFS - DIFSU = ADIFSU - DIFSH = ADIFSH - DAPAR = APAR - DNRAD = ANRAD !net absorbed radiation by canopy - DATRJC = ATRJC !absorbed radiation by canopy - DATRJS = ATRJS !absorbed radiation by soil - CUMRSD = RSD !incoming solar radiation - AV_RSWSU = RSWSU - AV_RSWSH = RSWSH - AV_Albedo = Albedo - TLEAFSU = TAIR + DIFSU - TLEAFSH = TAIR + DIFSH - - RETURN - END SUBROUTINE TOTPT - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PPHTR * -!* Purpose: This subroutine calculates potential leaf photosynthesis * -!* and transpiration. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* FRAC R4 Fraction of leaf classes (sunlit vs shaded) - I * -!* TAIR R4 Air temperature oC I * -!* DVP R4 Vapour pressure kPa I * -!* CO2A R4 Ambient CO2 concentration ml m-3 I * -!* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I * -!* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I * -!* PAR R4 Absorbed photosynth. active radiation J m-2 s-1 I * -!* NP R4 Photosynthetically active N content g m-2 I * -!* RBW R4 Leaf boundary layer resistance to water s m-1 I * -!* RBH R4 Leaf boundary layer resistance to heat s m-1 I * -!* RT R4 Turbulence resistance s m-1 I * -!* ATRJ R4 Absorbed global radiation J m-2 s-1 I * -!* ATMTR R4 Atmospheric transmissivity - I * -!* EAJMAX R4 Energy of activation for Jmax J mol-1 I * -!* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I * -!* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I * -!* THETA R4 Convexity for light response of e-transport - I * -!* PLF R4 Potential leaf photosynthesis gCO2/m2/s O * -!* PT R4 Potential leaf transpiration mm s-1 O * -!* PH R4 Potential leaf sensible heat flux J m-2 s-1 O * -!* RSW R4 Potential stomatal resistance to water s m-1 O * -!* NRADC R4 Net leaf absorbed radiation J m-2 s-1 O * -!* SLOPEL R4 Slope of saturated vapour pressure curve kPa oC-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE PPHTR(FRAC,TAIR,DVP,CO2A,C3C4,FVPD,PAR,NP,RBW,RBH,RT, & - ATRJ,ATMTR,EAJMAX,XVN,XJN,THETA,PLF,PT,PH,RSW,NRADC,SLOPEL) - IMPLICIT REAL (A-Z) - -!*---first-round calculation to determine leaf temperature - CALL ICO2 (TAIR,DVP,FVPD,CO2A,C3C4, SVP,FCO2I) - - CALL PHOTO(C3C4,PAR,TAIR,FCO2I,NP,EAJMAX,XVN,XJN,THETA,FPLF, & - FLRD) - - VPD = MAX (0., SVP- DVP) - SLOPE = 4158.6 * SVP/(TAIR + 239.)**2 - CALL GCRSW(FPLF,FLRD,TAIR,CO2A,FCO2I,RBW,RT, FRSW) - - CALL PTRAN(FRSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TAIR,DVP, & - SLOPE, VPD, FPT, FPH, FNRADC) - - CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF) - - TLEAF = TAIR + FDIF - -!*---second-round calculation to determine potential photosynthesis -!* and transpiration - CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I) - CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD) - - SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR) - - CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW) - CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, & - SLOPEL,VPD, PT, PH, NRADC) - - CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF) - - TLEAF = TAIR + FDIF - -!*---third-round calculation to determine potential photosynthesis -!* and transpiration - CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I) - CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD) - - SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR) - - CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW) - CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, & - SLOPEL,VPD, PT, PH, NRADC) - - CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF) - - TLEAF = TAIR + FDIF - -!*---fourth-round calculation to determine potential photosynthesis -!* and transpiration - CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I) - CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD) - - SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR) - - CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW) - CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, & - SLOPEL,VPD, PT, PH, NRADC) - - RETURN - END SUBROUTINE PPHTR - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PEVAP * -!* Purpose: This subroutine calculates potential soil evaporation. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* TAIR R4 Air temperature oC I * -!* DVP R4 Vapour pressure kPa I * -!* RSS R4 Soil resistance,equivalent to leaf stomata s m-1 I * -!* RTS R4 Turbulence resistance for soil s m-1 I * -!* RBWS R4 Soil boundary layer resistance to water s m-1 I * -!* RBHS R4 Soil boundary layer resistance to heat s m-1 I * -!* ATRJS R4 Absorbed global radiation by soil J m-2 s-1 I * -!* ATMTR R4 Atmospheric transmissivity - I * -!* PT1 R4 Potential leaf transpiration using water mm s-1 I * -!* from upper evaporative soil layer * -!* WSUP1 R4 Water supply from upper evaporative soil mm s-1 I * -!* layer for evapotranspiration * -!* PESOIL R4 Potential soil evaporation mm s-1 O * -!* PHSOIL R4 Potential soil sensible heat flux J m-2 s-1 O * -!* NRADS R4 Net soil absorbed radiation J m-2 s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE PEVAP (TAIR,DVP,RSS,RTS,RBWS,RBHS,ATRJS,ATMTR, & - PT1,WSUP1,PESOIL,PHSOIL,NRADS) - IMPLICIT REAL (A-Z) - -!*--- first-round calculation to estimate soil surface temperature (TAVS) - SVP = 0.611*EXP(17.4*TAIR/(TAIR+239.)) - VPD = MAX (0., SVP-DVP) - SLOPE = 4158.6 * SVP/(TAIR + 239.)**2 - CALL PTRAN(RSS,RTS,RBWS,RBHS,ATRJS,ATMTR,1.,TAIR,DVP, & - SLOPE,VPD, FPE, FPH, FNRADS) - FPESOL = MAX(0., FPE) - FAESOL = MIN(FPESOL,FPESOL/(PT1+FPESOL)*WSUP1) - CALL DIFLA (FNRADS,FAESOL,RBHS,RTS, FDIFS) - TAVS = TAIR + FDIFS - -!*---second-round calculation to estimate potential soil evaporation - SVPS = 0.611*EXP(17.4*TAVS/(TAVS+239.)) - SLOPES = (SVPS-SVP)/NOTNUL(FDIFS) - - CALL PTRAN(RSS,RTS,RBWS,RBHS,ATRJS,ATMTR,1.,TAVS,DVP, & - SLOPES,VPD, PE, PH, NRADS) - PESOIL = MAX(0., PE) - PHSOIL = PH - - RETURN - END SUBROUTINE PEVAP - -!*----------------------------------------------------------------------* -!* SUBROUTINE APHTR * -!* Purpose: This subroutine calculates actual leaf photosynthesis when * -!* water stress occurs. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* TAIR R4 Air temperature oC I * -!* PAR R4 Absorbed photosynth. active radiation J m-2 s-1 I * -!* DVP R4 Vapour pressure kPa I * -!* CO2A R4 Ambient CO2 concentration ml m-3 I * -!* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I * -!* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I * -!* NRADC R4 Net leaf absorbed radiation J m-2 s-1 I * -!* AT R4 Actual leaf transpiration mm s-1 I * -!* PT R4 Potential leaf transpiration mm s-1 I * -!* RT R4 Turbulence resistance s m-1 I * -!* RBH R4 Leaf boundary layer resistance to heat s m-1 I * -!* RBW R4 Leaf boundary layer resistance to water s m-1 I * -!* RSW R4 Potential stomatal resistance to water s m-1 I * -!* SLOPEL R4 Slope of saturated vapour pressure curve kPa oC-1 I * -!* NP R4 Photosynthet. active leaf N content g m-2 I * -!* NPN R4 NP with small plant-N increment g m-2 I * -!* EAJMAX R4 Energy of activation for Jmax J mol-1 I * -!* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I * -!* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I * -!* THETA R4 Convexity for light response of e-transport - I * -!* PLFAS R4 Actual leaf photosynthesis gCO2/m2/s O * -!* PLFAN R4 PLFAS with small plant-N increment gCO2/m2/s O * -!* ADIF R4 Actual leaf-air temperature difference oC O * -!* ARSW R4 Actual stomatal resistance to water s m-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE APHTR(TAIR,PAR,DVP,CO2A,C3C4,FVPD,NRADC,AT,PT,RT,RBH, & - RBW,RSW,SLOPEL,NP,NPN,EAJMAX,XVN,XJN,THETA, & - PLFAS,PLFAN,ADIF,ARSW) - IMPLICIT REAL (A-Z) - - PSYCH = 0.067 !psychrometric constant (kPa/oC) - -!*---leaf temperature if water stress occurs - CALL DIFLA (NRADC,AT,RBH,RT, ADIF) - ATLEAF = TAIR + ADIF - -!*---stomatal resistance to water if water stress occurs - ARSW = (PT-AT)*(SLOPEL*(RBH+RT)+PSYCH*(RBW+RT))/AT/PSYCH+PT/AT*RSW - -!*---potential photosynthesis at the new leaf temperature - CALL ICO2 (ATLEAF,DVP,FVPD,CO2A,C3C4, SVPA,ACO2I) - CALL PHOTO(C3C4,PAR,ATLEAF,ACO2I,NPN,EAJMAX,XVN,XJN,THETA,APLFN,ARDN) - CALL PHOTO(C3C4,PAR,ATLEAF,ACO2I,NP,EAJMAX,XVN,XJN,THETA,APLF,ARD) - -!*---actual photosynthesis under water stress condition - PLFAS = (1.6*RSW+1.3*RBW+RT)/(1.6*ARSW+1.3*RBW+RT)*(APLF-ARD)+ARD - PLFAN = (1.6*RSW+1.3*RBW+RT)/(1.6*ARSW+1.3*RBW+RT)*(APLFN-ARDN)+ARDN - - RETURN - END SUBROUTINE APHTR - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PTRAN * -!* Purpose: This subroutine calculates leaf transpiration, using the * -!* Penman-Monteith equation * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* RSW R4 Potential stomatal resistance to water s m-1 I * -!* RT R4 Turbulence resistance s m-1 I * -!* RBW R4 Leaf boundary layer resistance to water s m-1 I * -!* RBH R4 Leaf boundary layer resistance to heat s m-1 I * -!* ATRJ R4 Absorbed global radiation J m-2 s-1 I * -!* ATMTR R4 Atmospheric transmissivity - I * -!* FRAC R4 Fraction of leaf classes (sunlit vs shaded)- I * -!* TLEAF R4 Leaf temperature oC I * -!* DVP R4 Vapour pressure kPa I * -!* SLOPE R4 Slope of saturated vapour pressure curve kPa oC-1 I * -!* VPD R4 Saturation vapour pressure deficit of air kPa I * -!* PT R4 Potential leaf transpiration mm s-1 O * -!* PH R4 Potential leaf sensible heat flux J m-2 s-1 O * -!* NRADC R4 Net leaf absorbed radiation J m-2 s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE PTRAN (RSW,RT,RBW,RBH,ATRJ, & - ATMTR,FRAC,TLEAF,DVP,SLOPE,VPD,PT,PH,NRADC) - IMPLICIT REAL (A-Z) - -!*---some physical constants - BOLTZM = 5.668E-8 !Stefan-Boltzmann constant(J/m2/s/K4) - LHVAP = 2.4E6 !latent heat of water vaporization(J/kg) - VHCA = 1200. !volumetric heat capacity (J/m3/oC) - PSYCH = 0.067 !psychrometric constant (kPa/oC) - -!*---net absorbed radiation - CLEAR = MAX(0., MIN(1., (ATMTR-0.25)/0.45)) !sky clearness - BBRAD = BOLTZM*(TLEAF +273.)**4 - RLWN = BBRAD*(0.56-0.079*SQRT(DVP*10.))*(0.1+0.9*CLEAR)*FRAC - NRADC = ATRJ - RLWN - -!*---intermediate variable related to resistances - PSR = PSYCH*(RBW+RT+RSW)/(RBH+RT) - -!*---Compute PT -!*---radiation-determined term - PTR = NRADC*SLOPE/(SLOPE+PSR)/LHVAP - -!*---vapour pressure-determined term - PTD = (VHCA*VPD/(RBH+RT))/(SLOPE+PSR)/LHVAP - -!*---potential evaporation or transpiration - PT = MAX(1.E-10,PTR+PTD) - -!*---ji -!*---Compute PH in W/m2 (see Bolan p. 202) -!*---radiation-determined term - PHR = NRADC*PSR/(SLOPE+PSR) - -!*---vapour pressure-determined term - PHD = (VHCA*VPD/(RBH+RT))/(SLOPE+PSR) - -!*---potential evaporation or transpiration - PH = PHR-PHD - - RETURN - END SUBROUTINE PTRAN - - -!*----------------------------------------------------------------------* -!* SUBROUTINE DIFLA * -!* Purpose: This subroutine calculates leaf(canopy)-air temperature * -!* differential. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* NRADC R4 Net leaf absorbed radiation J m-2 s-1 I * -!* PT R4 Potential leaf transpiration mm s-1 I * -!* RBH R4 Leaf boundary layer resistance to heat s m-1 I * -!* RT R4 Turbulence resistance s m-1 I * -!* DIF R4 Leaf-air temperature difference oC O * -!*----------------------------------------------------------------------* - SUBROUTINE DIFLA (NRADC,PT,RBH,RT, DIF) - IMPLICIT REAL (A-Z) - - LHVAP = 2.4E6 !latent heat of water vaporization(J/kg) - VHCA = 1200. !volumetric heat capacity (J/m3/oC) - - DIF = LIMIT (-25., 25., (NRADC-LHVAP*PT)*(RBH+RT)/VHCA) - - RETURN - END SUBROUTINE DIFLA -!*----------------------------------------------------------------------* -!* SUBROUTINE ICO2 * -!* Purpose: This subroutine calculates the internal CO2 concentration * -!* as affected by vapour pressure deficit. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* TLEAF R4 Leaf temperature oC I * -!* DVP R4 Vapour pressure kPa I * -!* FVPD R4 Slope for linear effect of VPDL on Ci/Ca (kPa)-1 I * -!* CO2A R4 Ambient CO2 concentration ml m-3 I * -!* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I * -!* SVPL R4 Saturated vapour pressure of leaf kPa O * -!* CO2I R4 intercellular CO2 concentration ml m-3 O * -!*----------------------------------------------------------------------* - SUBROUTINE ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I) - IMPLICIT REAL (A-Z) - -!*---air-to-leaf vapour pressure deficit - SVPL = 0.611 * EXP(17.4 * TLEAF / (TLEAF + 239.)) - VPDL = MAX (0., SVPL - DVP) - -!*---Michaelis-Menten const. for CO2 at 25oC (umol/mol) - KMC25 = INSW(C3C4, 650., 404.9) !greater KMC25 for C4 than C3 - -!*---Michaelis-Menten const. for O2 at 25oC (mmol/mol) - KMO25 = INSW(C3C4, 450., 278.4) !greater KMO25 for C4 than C3 - -!*---CO2 compensation point in absence of dark respiration (GAMMAX) - O2 = 210. !oxygen concentration(mmol/mol) - EAVCMX = 65330. !energy of activation for Vcmx(J/mol) - EAKMC = 79430. !energy of activation for KMC (J/mol) - EAKMO = 36380. !energy of activation for KMO (J/mol) - EARD = 46390. !energy of activation for dark respiration(J/mol) - RDVX25 = 0.0089 !ratio of dark respiration to Vcmax at 25oC - TO = 298.15 - - KMC = KMC25*EXP((1./TO-1./(TLEAF+273.))*EAKMC/8.314) - KMO = KMO25*EXP((1./TO-1./(TLEAF+273.))*EAKMO/8.314) - GAMMAX = 0.5*EXP(-3.3801+5220./(TLEAF+273.)/8.314)*O2*KMC/KMO - -!*---CO2 compensation point (GAMMA) - RDVCX = RDVX25*EXP((1./TO-1./(TLEAF+273.))*(EARD-EAVCMX)/8.314) - GAMMA0 = (GAMMAX+RDVCX*KMC*(1.+O2/KMO))/(1.-RDVCX) - GAMMA_ = INSW (C3C4, GAMMA0/10., GAMMA0) - -!*---internal/ambient CO2 ratio, based on data of Morison & Gifford (1983) - RCICA = 1.-(1.-GAMMA_/CO2A)*(0.14+FVPD*VPDL) - -!*---intercellular CO2 concentration - CO2I = RCICA * CO2A - - RETURN - END SUBROUTINE ICO2 - - -!*----------------------------------------------------------------------* -!* SUBROUTINE GCRSW * -!* Purpose: This subroutine calculates overall leaf conductance * -!* for CO2 (GC) and the stomatal resistance to water (RSW). * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* PLEAF R4 Gross leaf photosynthesis gCO2/m2/s I * -!* RDLEAF R4 Leaf dark respiration gCO2/m2/s I * -!* TLEAF R4 Leaf temperature oC I * -!* CO2A R4 Ambient CO2 concentration ml m-3 I * -!* CO2I R4 Internal CO2 concentration ml m-3 I * -!* RT R4 Turbulence resistance s m-1 I * -!* RBW R4 Leaf boundary layer resistance to water s m-1 I * -!* RSW R4 Potential stomatal resistance to water s m-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE GCRSW (PLEAF,RDLEAF,TLEAF,CO2A,CO2I,RBW,RT, RSW) - IMPLICIT REAL (A-Z) - -!*---potential conductance for CO2 -!*ji 13.7.11 MAX routine faengt negative Leitfaehigkeiten ab - GC = MAX(1E-6,(PLEAF-RDLEAF)*(273.+TLEAF)/0.53717/(CO2A-CO2I)) - -!*---potential stomatal resistance to water - RSW = MAX(1E-30, 1./GC - RBW*1.3 - RT)/1.6 - - RETURN - END SUBROUTINE GCRSW - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PAN * -!* Purpose: This subroutine calculates photosynthetically active * -!* nitrogen content for sunlit and shaded parts of canopy. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* SLNT R4 Top-leaf nitrogen content g m-2 I * -!* SLNMIN R4 Minimum or base SLNT for photosynthesis g m-2 I * -!* LAI R4 (green)Leaf area index m2 m-2 I * -!* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I * -!* KB R4 Direct beam radiation extinction coeff. m2 m-2 I * -!* NPSU R4 Photosynthet. active N for sunlit leaves g m-2 O * -!* NPSH R4 Photosynthet. active N for shaded leaves g m-2 O * -!*----------------------------------------------------------------------* - SUBROUTINE PAN(SLNT,SLNMIN,LAI,KN,KB, NPSU,NPSH) - IMPLICIT REAL (A-Z) - -!*---total photosynthetic nitrogen in canopy - NPC = SLNT*(1.-EXP(-KN*LAI))/KN-SLNMIN*LAI - -!*---photosynthetic nitrogen for sunlit and shaded parts of canopy - NPSU = SLNT*(1.-EXP(-(KN+KB)*LAI))/(KN+KB)-SLNMIN*(1.-EXP(-KB*LAI))/KB - NPSH = NPC-NPSU - - RETURN - END SUBROUTINE PAN - - -!*----------------------------------------------------------------------* -!* SUBROUTINE PHOTO * -!* Purpose: This subroutine calculates leaf photosynthesis and dark * -!* respiration, based on a renewed Farquhar biochemistry * -!* (cf Yin et al.2004. Plant, Cell & Environment 27:1211-1222)* -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* C3C4 R4 Crop type (=1. for C3, -1. for C4 crops) - I * -!* PAR R4 Leaf absorbed photosynth. active radiance J m-2 s-1 I * -!* TLEAF R4 Leaf temperature oC I * -!* CO2I R4 Intercellular CO2 concentration ml m-3 I * -!* NP R4 Photosynthetically active leaf N content g m-2 I * -!* EAJMAX R4 Energy of activation for Jmax J mol-1 I * -!* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I * -!* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I * -!* THETA R4 Convexity for light response of e-transport - I * -!* PLEAF R4 Gross leaf photosynthesis gCO2/m2/s O * -!* RDLEAF R4 Leaf dark respiration gCO2/m2/s O * -!*----------------------------------------------------------------------* - SUBROUTINE PHOTO(C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN, & - THETA,PLEAF,RDLEAF) - IMPLICIT REAL (A-Z) - -!*---Michaelis-Menten constants for CO2 and O2 at 25oC - IF (C3C4.LT.0.) THEN - KMC25 = 650. !greater KMC25 for C4 than C3; unit:(umol/mol) - KMO25 = 450. !greater KMO25 for C4 than C3; unit:(mmol/mol) - ELSE - KMC25 = 404.9 !unit:(umol/mol) - KMO25 = 278.4 !unit:(mmol/mol) - ENDIF - -!*---other constants related to the Farquhar-type photosynthesis model - O2 = 210. !oxygen concentration(mmol/mol) - EAVCMX = 65330. !energy of activation for Vcmx(J/mol) - EAKMC = 79430. !energy of activation for KMC (J/mol) - EAKMO = 36380. !energy of activation for KMO (J/mol) - EARD = 46390. !energy of activation for dark respiration(J/mol) - DEJMAX = 200000. !energy of deactivation for JMAX (J/mol) - SJ = 650. !entropy term in JT equation (J/mol/K) - PHI2M = 0.85 !maximum electron transport efficiency of PS II - HH = 3. !number of protons required to synthesise 1 ATP - KTMP = 1.0 !Factor for reducing photosynthesis in case of T<5C - JTMAX = 3.12 - TO = 298.15 - -!*---PAR photon flux in umol/m2/s absorbed by leaf photo-systems - UPAR = 4.56*PAR !4.56 conversion factor in umol/J - -!*---Michaelis-Menten constants for CO2 and O2 respectively - KMC = KMC25*EXP((1./TO-1./(TLEAF+273.))*EAKMC/8.314) - KMO = KMO25*EXP((1./TO-1./(TLEAF+273.))*EAKMO/8.314) - -!*---CO2 compensation point in the absence of dark respiration - GAMMAX = 0.5*EXP(-3.3801+5220./(TLEAF+273.)/8.314)*O2*KMC/KMO - -!*---Arrhenius function for the effect of temperature on carboxylation - VCT = EXP((1./TO-1./(TLEAF+273.))*EAVCMX/8.314) - -!*---function for the effect of temperature on electron transport - JT = EXP((1./TO-1./(TLEAF+273.))*EAJMAX/8.314)* & - (1.+EXP(SJ/8.314-DEJMAX/TO/8.314))/ & - (1.+EXP(SJ/8.314-1./(TLEAF+273.) *DEJMAX/8.314)) - - VCMX = XVN*VCT*NP - JMAX = XJN*JT *NP - -!*---CO2 concentration at carboxylation site & electron pathways and -!* their stoichiometries - FPSEUD = 0. !assuming no pseudocyclic e- transport - IF (C3C4.LT.0.) THEN - ZZ = 0.2 !CO2 leakage from bundle-sheath to mesophyll - CC = 10.*CO2I !to mimic C4 CO2 concentrating mechanism - SF = 2.*(CC-GAMMAX)/(1.-ZZ) - FQ = 1.- FPSEUD- 2.*(4.*CC+8.*GAMMAX)/HH/(SF+3.*CC+7.*GAMMAX) - FCYC = FQ - ELSE - CC = CO2I - SF = 0. - FQ = 0. - FCYC = 1.-(FPSEUD*HH*(SF+3.*CC+7.*GAMMAX)/(4.*CC+8.*GAMMAX)+1.)/ & - (HH*(SF+3.*CC+7.*GAMMAX)/(4.*CC+8.*GAMMAX)-1.) - ENDIF - -!*--- electron transport rate in dependence on PAR photon flux - ALPHA2 = (1.-FCYC)/(1.+(1.-FCYC)/PHI2M) - X = ALPHA2*UPAR/MAX(1.E-10,JMAX) - J2 = JMAX*(1.+X-((1.+X)**2.-4.*X*THETA)**0.5)/2./THETA - -!*---rates of carboxylation limited by Rubisco and electron transport - VC = VCMX * CC/(CC + KMC*(O2/KMO+1.)) - VJ = J2 * CC*(2.+FQ-FCYC)/HH/(SF+3.*CC+7.*GAMMAX)/(1.-FCYC) - -!*---gross rate of leaf photosynthesis - ALF = (1.-GAMMAX/CC)*MIN(VC,VJ) - PLEAF = MAX(1.E-10, (1.E-6)*44.*ALF) - -!*---rate of leaf dark respiration - RDVX25 = 0.0089 !ratio of dark respiration to Vcmax at 25oC - RDT = EXP((1./TO-1./(TLEAF+273.))*EARD/8.314) - RDLEAF = (1.E-6)*44. *RDVX25*(XVN*NP) * RDT - - RETURN - END SUBROUTINE PHOTO - - -!*----------------------------------------------------------------------* -!* SUBROUTINE REFL * -!* Purpose: This subroutine calculates reflection coefficients. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* SCP R4 Leaf scattering coefficient - I * -!* KB R4 Direct beam radiation extinction coeff. m2 m-2 I * -!* KBP R4 Scattered beam radiation extinction coeff. m2 m-2 O * -!* PCB R4 Canopy beam radiation reflection coeff. - O * -!*----------------------------------------------------------------------* - SUBROUTINE REFL (SCP,KB, KBP,PCB) - IMPLICIT REAL (A-Z) - -!*--- scattered beam radiation extinction coefficient - KBP = KB*SQRT(1.-SCP) - -!*---canopy reflection coefficient for horizontal leaves - PH = (1.-SQRT(1.-SCP))/(1.+SQRT(1.-SCP)) - -!*---Canopy beam radiation reflection coefficient - PCB = 1.-EXP(-2.*PH*KB/(1.+KB)) - - RETURN - END SUBROUTINE REFL - - -!*----------------------------------------------------------------------* -!* SUBROUTINE LIGAB * -!* Purpose: This subroutine calculates absorbed light for sunlit and * -!* shaded leaves. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* SCP R4 Leaf scattering coefficient - I * -!* KB R4 Direct beam radiation extinction coeff. m2 m-2 I * -!* KBP R4 Scattered beam radiation extinction coeff. m2 m-2 I * -!* KDP R4 Diffuse radiation extinction coefficient m2 m-2 I * -!* PCB R4 Canopy beam radiation reflection coeff. - I * -!* PCD R4 Canopy diffuse radiation reflection coeff. - I * -!* IB0 R4 Incident direct-beam radiation J m-2 s-1 I * -!* ID0 R4 Incident diffuse radiation J m-2 s-1 I * -!* LAI R4 (green)Leaf area index m2 m-2 I * -!* ISU R4 Absorbed radiation by sunlit leaves J m-2 s-1 O * -!* ISH R4 Absorbed radiation by shaded leaves J m-2 s-1 O * -!*----------------------------------------------------------------------* - SUBROUTINE LIGAB (SCP,KB,KBP,KDP,PCB,PCD,IB0,ID0,LAI, ISU,ISH) - IMPLICIT REAL (A-Z) - -!*---total absorbed light by canopy - IC = (1.-PCB)*MAX(1e-30,IB0)*(1.-EXP(-KBP*LAI))+ & - (1.-PCD)*MAX(1e-30,ID0)*(1.-EXP(-KDP*LAI)) - -!*---absorbed light by sunlit and shaded fractions of canopy - ISU = (1.-SCP)*MAX(1e-30,IB0)*(1.-EXP(-KB *LAI))+(1.-PCD)*MAX(1e-30,ID0)/(KDP+KB)* & - KDP*(1.-EXP(-(KDP+KB)*LAI))+MAX(1e-30,IB0)*((1.-PCB)/(KBP+KB)*KBP* & - (1.-EXP(-(KBP+KB)*LAI))-(1.-SCP)*(1.-EXP(-2.*KB*LAI))/2.) - - ISH = IC-ISU - - RETURN - END SUBROUTINE LIGAB - - -!*----------------------------------------------------------------------* -!* SUBROUTINE KBEAM * -!* Purpose: This subroutine calculates extinction coefficient for * -!* direct beam radiation. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* SINB R4 Sine of solar elevation - I * -!* BL R4 Leaf angle (from horizontal) radians I * -!* KB R4 Direct beam radiation extinction coeff. m2 m-2 O * -!*----------------------------------------------------------------------* - SUBROUTINE KBEAM (SINB,BL, KB) - IMPLICIT REAL (A-Z) - -!*---solar elevation in radians - B = ASIN(SINB) - -!*---average projection of leaves in the direction of a solar beam - IF (SINB.GE.SIN(BL)) THEN - OAV = SINB*COS(BL) - ELSE - OAV = 2./3.141592654*(SINB*COS(BL)*ASIN(TAN(B)/TAN(BL)) & - +((SIN(BL))**2-SINB**2)**0.5) - ENDIF - -!*---beam radiation extinction coefficient - KB = OAV/SINB - - RETURN - END SUBROUTINE KBEAM - - -!*----------------------------------------------------------------------* -!* SUBROUTINE KDIFF * -!* Purpose: This subroutine calculates extinction coefficient for * -!* diffuse radiation. * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* LAI R4 Total leaf area index m2 m-2 I * -!* BL R4 Leaf angle (from horizontal) radians I * -!* SCP R4 Leaf scattering coefficient - I * -!* KDP R4 Diffuse radiation extinction coefficient m2 m-2 O * -!*----------------------------------------------------------------------* - SUBROUTINE KDIFF (LAI,BL,SCP, KDP) - IMPLICIT REAL (A-Z) - - PI = 3.141592654 - -!*---extinction coefficient of beam lights from 15, 45 and 75 elevations - CALL KBEAM (SIN(15.*PI/180.),BL, KB15) - CALL KBEAM (SIN(45.*PI/180.),BL, KB45) - CALL KBEAM (SIN(75.*PI/180.),BL, KB75) - -!*---diffuse light extinction coefficient - KDP = -1./LAI*LOG(0.178*EXP(-KB15*(1.-SCP)**0.5*LAI) & - +0.514*EXP(-KB45*(1.-SCP)**0.5*LAI) & - +0.308*EXP(-KB75*(1.-SCP)**0.5*LAI)) - - RETURN - END SUBROUTINE KDIFF - -!*----------------------------------------------------------------------* -!* FUNCTION INTGRL * -!* Purpose: This function integrates a differential equation * -!* using the Euler method. Substitutes the intrinsic FST function * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* POOL R4 Currenmt stock size g X m-2 I * -!* RATE R4 Rate of change 1/d I * -!* DT R4 Time step d I * -!* INTGRL R4 Returns the stock size in next time step g X m-2 0 * -FUNCTION INTGRL(POOL, RATE, DT) -IMPLICIT NONE -real:: INTGRL, POOL, RATE, DT - -INTGRL = POOL + RATE*DT - -RETURN -END FUNCTION INTGRL - -!*----------------------------------------------------------------------* -!* FUNCTION NOTNUL * -!* Y = NOTNUL(X) * -!* Y is equal to X but 1.0 in case of x=0.0. Note that X is * -!* evaluated without any tolerance interval * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* X R4 Variable to be checked I * -!*----------------------------------------------------------------------* -FUNCTION NOTNUL(X) -IMPLICIT NONE -real:: NOTNUL, X - -if (X.ne.0) then - NOTNUL = X -else - NOTNUL = 1 -endif - -RETURN -END FUNCTION NOTNUL - -!*----------------------------------------------------------------------* -!* FUNCTION LIMIT * -!* Y = LIMIT(XL, XH, X) * -!* Y is equal to X but limited between XL and XH * -!* Y - Returned as X bounded on [Xl,XH] * -!* XL- Lower bound of X * -!* XH- Upper bound of X * -!* ---- ---- ------- * -!* X R4 Variable to be checken I * -!*----------------------------------------------------------------------* -FUNCTION LIMIT(XL, XH, X) -IMPLICIT NONE -real:: LIMIT, XL, XH, X - -if (X.ge.XL.and.X.le.XH) then - LIMIT = X -else - if (X.gt.XH) then - LIMIT = XH - else - LIMIT = XL - endif -endif - -RETURN -END FUNCTION LIMIT - -!*----------------------------------------------------------------------* -!* FUNCTION INSW * -!* Y = INSW(X, Y1, Y2) * -!* Input switch. Y is set equal to Y1 orY2 depending on the value of X * -!* Y - Returned as either Y1 or Y2 * -!* X - Control variable * -!* Y1- Returned value of Y if X<0 * -!* Y2- Returned value of Y if X>=0 * -!* ---- ---- ------- * -!* X R4 Variable to be checken I * -!*----------------------------------------------------------------------* -FUNCTION INSW(X, Y1, Y2) -implicit none -real :: INSW, X, Y1, Y2 - -if (X.lt.0) then - INSW = Y1 -else - INSW = Y2 -endif - -RETURN -END FUNCTION INSW - - -!*----------------------------------------------------------------------* -!* FUNCTION REAAND * -!* Y = REAAND(X1, X2) * -!* Returns 1.0 if both input variables are positive, otherwise Y=0.0 * -!* ---- ---- ------- * -!* X1 R4 1. variable to be checked I * -!* X2 R4 2. variable to be checked I * -!*----------------------------------------------------------------------* -FUNCTION REAAND(X1, X2) -implicit none -real :: REAAND, X1, X2 - -if (X1.gt.0.and.X2.gt.0) then - REAAND = 1. -else - REAAND = 0. -endif - -RETURN -END FUNCTION REAAND - - -!*----------------------------------------------------------------------* -!* FUNCTION REAnOR * -!* Y = REANOR(X1, X2) * -!* Returns 1.0 if both input variables are less than or equal to 0., * -!* otherwise Y=0. -!* ---- ---- ------- * -!* X1 R4 1. variable to be checked I * -!* X2 R4 2. variable to be checked I * -!*----------------------------------------------------------------------* -FUNCTION REANOR(X1, X2) -implicit none -real :: REANOR, X1, X2 - -if (X1.le.0.and.X2.le.0) then - REANOR = 1. -else - REANOR = 0. -endif - -RETURN -END FUNCTION REANOR - -END MODULE module_sf_gecros diff --git a/phys/module_sf_noahmpdrv.F b/phys/module_sf_noahmpdrv.F index f77b959845..14a7fb5959 100644 --- a/phys/module_sf_noahmpdrv.F +++ b/phys/module_sf_noahmpdrv.F @@ -32,7 +32,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN ZSNSOXY, SNICEXY, SNLIQXY, LFMASSXY, RTMASSXY, STMASSXY, & ! IN/OUT Noah MP only WOODXY, STBLCPXY, FASTCPXY, XLAIXY, XSAIXY, TAUSSXY, & ! IN/OUT Noah MP only SMOISEQ, SMCWTDXY,DEEPRECHXY, RECHXY, GRAINXY, GDDXY,PGSXY, & ! IN/OUT Noah MP only - GECROS_STATE, & ! IN/OUT gecros model T2MVXY, T2MBXY, Q2MVXY, Q2MBXY, & ! OUT Noah MP only TRADXY, NEEXY, GPPXY, NPPXY, FVEGXY, RUNSFXY, & ! OUT Noah MP only RUNSBXY, ECANXY, EDIRXY, ETRANXY, FSAXY, FIRAXY, & ! OUT Noah MP only @@ -100,7 +99,7 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN INTEGER, INTENT(IN ) :: IOPT_RSF ! surface resistance (1->Sakaguchi/Zeng; 2->Seller; 3->mod Sellers; 4->1+snow) INTEGER, INTENT(IN ) :: IOPT_SOIL ! soil configuration option INTEGER, INTENT(IN ) :: IOPT_PEDO ! soil pedotransfer function option - INTEGER, INTENT(IN ) :: IOPT_CROP ! crop model option (0->none; 1->Liu et al.; 2->Gecros) + INTEGER, INTENT(IN ) :: IOPT_CROP ! crop model option (0->none; 1->Liu et al. 2016) INTEGER, INTENT(IN ) :: IZ0TLND ! option of Chen adjustment of Czil (not used) INTEGER, INTENT(IN ) :: sf_urban_physics ! urban physics option REAL, DIMENSION( ims:ime, 8, jms:jme ), INTENT(IN ) :: SOILCOMP ! soil sand and clay percentage @@ -135,8 +134,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: GDDXY ! growing degree days XING (based on 10C) INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PGSXY -! gecros model - REAL, DIMENSION( ims:ime, 60,jms:jme ), INTENT(INOUT) :: gecros_state ! gecros crop #ifdef WRF_HYDRO REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfcheadrt,INFXSRT,soldrain ! for WRF-Hydro @@ -420,10 +417,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN ! Local terms - REAL, DIMENSION(1:60) :: gecros1d ! gecros crop - REAL :: gecros_dd ,gecros_tbem,gecros_emb ,gecros_ema, & - gecros_ds1,gecros_ds2 ,gecros_ds1x,gecros_ds2x - REAL :: FSR ! total reflected solar radiation (w/m2) REAL, DIMENSION(-2:0) :: FICEOLD ! snow layer ice fraction [] REAL :: CO2PP ! CO2 partial pressure [Pa] @@ -632,34 +625,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN RECH = 0. DEEPRECH = 0. - if(iopt_crop == 2) then ! gecros crop model - - gecros1d(1:60) = gecros_state(I,1:60,J) ! Gecros variables 2D -> local - - if(croptype == 1) then - gecros_dd = 2.5 - gecros_tbem = 2.0 - gecros_emb = 10.2 - gecros_ema = 40.0 - gecros_ds1 = 2.1 !BBCH 92 - gecros_ds2 = 2.0 !BBCH 90 - gecros_ds1x = 0.0 - gecros_ds2x = 10.0 - end if - - if(croptype == 2) then - gecros_dd = 5.0 - gecros_tbem = 8.0 - gecros_emb = 15.0 - gecros_ema = 6.0 - gecros_ds1 = 1.78 !BBCH 85 - gecros_ds2 = 1.63 !BBCH 80 - gecros_ds1x = 0.0 - gecros_ds2x = 14.0 - end if - - end if - SLOPETYP = 1 ! set underground runoff slope term IST = 1 ! MP surface type: 1 = land; 2 = lake SOILCOLOR = 4 ! soil color: assuming a middle color category ????????? @@ -864,7 +829,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN CM , CH , TAUSS , & ! IN/OUT : GRAIN , GDD , PGS , & ! IN/OUT SMCWTD ,DEEPRECH , RECH , & ! IN/OUT : - GECROS1D, & ! IN/OUT : Z0WRF , & FSA , FSR , FIRA , FSH , SSOIL , FCEV , & ! OUT : FGEV , FCTR , ECAN , ETRAN , ESOIL , TRAD , & ! OUT : @@ -1018,19 +982,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN GDDXY (I,J) = GDD !XING PGSXY (I,J) = PGS - if(iopt_crop == 2) then ! gecros crop model - - !*** Check for harvest - if ((gecros1d(1) >= gecros_ds1).and.(gecros1d(42) < 0)) then - if (checkIfHarvest(gecros_state, DT, gecros_ds1, gecros_ds2, gecros_ds1x, & - gecros_ds2x) == 1) then - - call gecros_reinit(gecros1d) - endif - endif - - gecros_state (i,1:60,j) = gecros1d(1:60) - end if ENDIF ! endif of land-sea test @@ -1400,12 +1351,9 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT smoiseq ,smcwtdxy ,rechxy ,deeprechxy, areaxy, dx, dy, msftx, msfty,& ! Optional groundwater wtddt ,stepwtd ,dt ,qrfsxy ,qspringsxy , qslatxy , & ! Optional groundwater fdepthxy ,ht ,riverbedxy ,eqzwt ,rivercondxy ,pexpxy , & ! Optional groundwater - rechclim, & ! Optional groundwater - gecros_state) ! Optional gecros crop + rechclim ) ! Optional groundwater USE NOAHMP_TABLES - use module_sf_gecros, only: seednc,sla0,slnmin,ffat,flig,foac,fmin,npl,seedw,eg,fcrsh,seednc,lnci,cfv - IMPLICIT NONE @@ -1499,8 +1447,6 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT), OPTIONAL :: PEXPXY !factor for river conductance REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) , OPTIONAL :: rechclim - REAL, DIMENSION(ims:ime,60,jms:jme), INTENT(INOUT), OPTIONAL :: gecros_state ! Optional gecros crop - INTEGER, INTENT(OUT) , OPTIONAL :: STEPWTD REAL, INTENT(IN) , OPTIONAL :: DT, WTDDT @@ -1517,10 +1463,6 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT REAL :: BEXP, SMCMAX, PSISAT REAL :: FK, masslai,masssai -! gecros local variables - REAL :: hti,rdi,fpro,lncmin,fcar,cfo,clvi,crti,ygo,nlvi,laii,nrti,slnbi - - REAL, PARAMETER :: BLIM = 5.5 REAL, PARAMETER :: HLICE = 3.335E5 REAL, PARAMETER :: GRAV = 9.81 @@ -1750,34 +1692,6 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT end if end if -! Initialize cropcat for gecros crop model - - if(iopt_crop == 2) then - cropcat (i,j) = 0 - if(croptype(i,5,j) >= 0.5) then - if(croptype(i,3,j) > 0.0) cropcat(i,j) = 1 ! if any wheat, set to wheat - if(croptype(i,1,j) > croptype(i,3,j)) cropcat(i,j) = 2 ! change to maize - end if - - hti = 0.01 - rdi = 10. - fpro = 6.25*seednc - lncmin = sla0*slnmin - fcar = 1.-fpro-ffat-flig-foac-fmin - cfo = 0.444*fcar+0.531*fpro+0.774*ffat+0.667*flig+0.368*foac - clvi = npl * seedw * cfo * eg * fcrsh - crti = npl * seedw * cfo * eg * (1.-fcrsh) - ygo = cfo/(1.275*fcar+1.887*fpro+3.189*ffat+2.231*flig+0.954* & - foac)*30./12. - nlvi = min(0.75 * npl * seedw * eg * seednc, lnci * clvi/cfv) - laii = clvi/cfv*sla0 - nrti = npl * seedw * eg * seednc - nlvi - slnbi = nlvi/laii - - call gecros_init(xlat(i,j),hti,rdi,clvi,crti,nlvi,laii,nrti,slnbi,gecros_state(i,:,j)) - - end if - END IF enddo @@ -2288,189 +2202,6 @@ SUBROUTINE EQSMOISTURE(NSOIL , ZSOIL , SMCMAX , SMCWLT, DWSAT , DKSAT ,BEXP , END SUBROUTINE EQSMOISTURE -! gecros initialization routines - -SUBROUTINE gecros_init(xlat,hti,rdi,clvi,crti,nlvi,laii,nrti,slnbi,state_gecros) -implicit none -REAL, INTENT(IN) :: HTI -REAL, INTENT(IN) :: RDI -REAL, INTENT(IN) :: CLVI -REAL, INTENT(IN) :: CRTI -REAL, INTENT(IN) :: NLVI -REAL, INTENT(IN) :: LAII -REAL, INTENT(IN) :: NRTI -REAL, INTENT(IN) :: SLNBI -REAL, INTENT(IN) :: XLAT -REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS - - !Inititalization of Gecros variables - STATE_GECROS(1) = 0. !DS - STATE_GECROS(2) = 0. !CTDURDI, HTI, CLVI, CRTI, NLVI, LAII, NRTI, SLNBI, - STATE_GECROS(3) = 0. !CVDU - STATE_GECROS(4) = CLVI !CLV - STATE_GECROS(5) = 0. !CLVD - STATE_GECROS(6) = 0. !CSST - STATE_GECROS(7) = 0. !CSO - STATE_GECROS(8) = CRTI !CSRT - STATE_GECROS(9) = 0. !CRTD - STATE_GECROS(10) = 0. !CLVDS - STATE_GECROS(11) = NRTI !NRT - STATE_GECROS(12) = 0. !NST - STATE_GECROS(13) = NLVI !NLV - STATE_GECROS(14) = 0. !NSO - STATE_GECROS(15) = NLVI !TNLV - STATE_GECROS(16) = 0. !NLVD - STATE_GECROS(17) = 0. !NRTD - STATE_GECROS(18) = 0. !CRVS - STATE_GECROS(19) = 0. !CRVR - STATE_GECROS(20) = 0. !NREOE - STATE_GECROS(21) = 0. !NREOF - STATE_GECROS(22) = 0. !DCDSR - STATE_GECROS(23) = 0. !DCDTR - STATE_GECROS(24) = SLNBI !SLNB - STATE_GECROS(25) = LAII !LAIC - STATE_GECROS(26) = 0. !RMUL - STATE_GECROS(27) = 0. !NDEMP - STATE_GECROS(28) = 0. !NSUPP - STATE_GECROS(29) = 0. !NFIXT - STATE_GECROS(30) = 0. !NFIXR - STATE_GECROS(31) = 0. !DCDTP - STATE_GECROS(32) = 0.01 !HTI - STATE_GECROS(33) = RDI !RDI - STATE_GECROS(34) = 0. !TPCAN - STATE_GECROS(35) = 0. !TRESP - STATE_GECROS(36) = 0. !TNUPT - STATE_GECROS(37) = 0. !LITNT - STATE_GECROS(38) = 0. !daysSinceDS1 - STATE_GECROS(39) = 0. !daysSinceDS2 - STATE_GECROS(40) = -1. !drilled -1:false, 1:true - STATE_GECROS(41) = -1. !emerged -1:false, 1:true - STATE_GECROS(42) = -1. !harvested -1:false, 1:true - STATE_GECROS(43) = 0. !TTEM - STATE_GECROS(44) = XLAT !GLAT - STATE_GECROS(45) = 0. !WSO - STATE_GECROS(46) = 0. !WSTRAW - STATE_GECROS(47) = 0. !GrainNC - STATE_GECROS(48) = 0. !StrawNC - STATE_GECROS(49) = 0.01 !GLAI - STATE_GECROS(50) = 0.01 !TLAI - STATE_GECROS(51) = HTI !Fields 51-58 set for reinitialization - STATE_GECROS(52) = RDI - STATE_GECROS(53) = CLVI - STATE_GECROS(54) = CRTI - STATE_GECROS(55) = NRTI - STATE_GECROS(56) = NLVI - STATE_GECROS(57) = SLNBI - STATE_GECROS(58) = LAII - -END SUBROUTINE gecros_init - -SUBROUTINE gecros_reinit(STATE_GECROS) -implicit none -REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS - - !Re-inititalization of Gecros variables after harvest - STATE_GECROS(1) = 0. !DS - STATE_GECROS(2) = 0. !CTDU - STATE_GECROS(3) = 0. !CVDU - STATE_GECROS(4) = STATE_GECROS(53) !CLV - STATE_GECROS(5) = 0. !CLVD - STATE_GECROS(6) = 0. !CSST - STATE_GECROS(7) = 0. !CSO - STATE_GECROS(8) = STATE_GECROS(54) !CRT - STATE_GECROS(9) = 0. !CRTD - STATE_GECROS(10) = 0. !CLVDS - STATE_GECROS(11) = STATE_GECROS(55)!NRT - STATE_GECROS(12) = 0. !NST - STATE_GECROS(13) = STATE_GECROS(56)!NLV - STATE_GECROS(14) = 0. !NSO - STATE_GECROS(15) = STATE_GECROS(56)!TNLV - STATE_GECROS(16) = 0. !NLVD - STATE_GECROS(17) = 0. !NRTD - STATE_GECROS(18) = 0. !CRVS - STATE_GECROS(19) = 0. !CRVR - STATE_GECROS(20) = 0. !NREOE - STATE_GECROS(21) = 0. !NREOF - STATE_GECROS(22) = 0. !DCDSR - STATE_GECROS(23) = 0. !DCDTR - STATE_GECROS(24) = STATE_GECROS(57)!SLNB - STATE_GECROS(25) = STATE_GECROS(58)!LAIC - STATE_GECROS(26) = 0. !RMUL - STATE_GECROS(27) = 0. !NDEMP - STATE_GECROS(28) = 0. !NSUPP - STATE_GECROS(29) = 0. !NFIXT - STATE_GECROS(30) = 0. !NFIXR - STATE_GECROS(31) = 0. !DCDTP - STATE_GECROS(32) = STATE_GECROS(51)!HT - STATE_GECROS(33) = STATE_GECROS(52)!ROOTD - STATE_GECROS(34) = 0. !TPCAN - STATE_GECROS(35) = 0. !TRESP - STATE_GECROS(36) = 0. !TNUPT - STATE_GECROS(37) = 0. !LITNT - STATE_GECROS(38) = 0. !daysSinceDS1 - STATE_GECROS(39) = 0. !daysSinceDS2 - STATE_GECROS(40) = -1. !drilled -1:false, 1:true - STATE_GECROS(41) = -1. !emerged -1:false, 1:true - STATE_GECROS(42) = 1. !harvested -1:false, 1:true - STATE_GECROS(43) = 0. !TTEM - STATE_GECROS(45) = 0. !WSO - STATE_GECROS(46) = 0. !WSTRAW - STATE_GECROS(47) = 0. !GrainNC - STATE_GECROS(48) = 0. !StrawNC - STATE_GECROS(49) = 0.01 !GLAI - STATE_GECROS(50) = 0.01 !TLAI - -END SUBROUTINE gecros_reinit - -!***Function for HARVEST DATES: - -!Determine if crop is to be harvested today -!function to be called once a day -!return codes: 0 - no, 1- yes -!requires two counters 'daysSinceDS2', 'daysSinceDS1' , zero-initialized to be maintained within caller -!STATE_GECROS(1) = current DS -!STATE_GECROS(38)=daysSinceDS1 -!STATE_GECROS(39)=daysSinceDS2 - -function checkIfHarvest(STATE_GECROS, DT, harvestDS1, harvestDS2, harvestDS1ExtraDays, harvestDS2ExtraDays) -implicit none -real :: DT, harvestDS1, harvestDS2 -real :: daysSinceDS1, daysSinceDS2 -real :: harvestDS1ExtraDays, harvestDS2ExtraDays -integer :: checkIfHarvest -REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS - - - !***check whether maturity (DS1) has been reached - if (STATE_GECROS(1) >= harvestDS1) then - - if (STATE_GECROS(38) >= harvestDS1ExtraDays) then - checkIfHarvest=1 - !if we are > DS1, but not over the limit, increase the counter of days - else - STATE_GECROS(38) = STATE_GECROS(38) + DT/86400. - endif - else - - !if maturity has not been reached, but we are close (> DS2) - !check the number of days for which we have been > DS2 - !and harvest in case we are over the limit given for that stage - !(in case that maturity will not be reached at all) - - checkIfHarvest=0 - if (STATE_GECROS(1) >= harvestDS2 ) then - - if (STATE_GECROS(39) >= harvestDS2ExtraDays) then - checkIfHarvest=1 - else !if we are > DS2, but not over the limit, increase the counter of days - STATE_GECROS(39) = STATE_GECROS(39) + DT/86400. - checkIfHarvest=0 - endif - endif - endif - return -end function checkIfHarvest - !------------------------------------------------------------------------------------------ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, & ! IN : Model configuration diff --git a/phys/module_sf_noahmplsm.F b/phys/module_sf_noahmplsm.F index 9686a403c6..320b3c20ab 100644 --- a/phys/module_sf_noahmplsm.F +++ b/phys/module_sf_noahmplsm.F @@ -20,8 +20,6 @@ MODULE MODULE_SF_NOAHMPLSM - use module_sf_gecros, only : gecros - IMPLICIT NONE public :: noahmp_options @@ -170,7 +168,6 @@ MODULE MODULE_SF_NOAHMPLSM INTEGER :: OPT_CROP ! options for crop model ! **0 -> No crop model, will run default dynamic vegetation ! 1 -> Liu, et al. 2016 - ! 2 -> Gecros (Genotype-by-Environment interaction on CROp growth Simulator) Yin and van Laar, 2005 !------------------------------------------------------------------------------------------! ! Physical Constants: ! @@ -400,7 +397,6 @@ SUBROUTINE NOAHMP_SFLX (parameters, & CM , CH , TAUSS , & ! IN/OUT : GRAIN , GDD , PGS , & ! IN/OUT SMCWTD ,DEEPRECH , RECH , & ! IN/OUT : - GECROS1D, & ! IN/OUT : Z0WRF , & FSA , FSR , FIRA , FSH , SSOIL , FCEV , & ! OUT : FGEV , FCTR , ECAN , ETRAN , EDIR , TRAD , & ! OUT : @@ -508,8 +504,6 @@ SUBROUTINE NOAHMP_SFLX (parameters, & REAL, INTENT(INOUT) :: DEEPRECH !recharge to or from the water table when deep [m] REAL, INTENT(INOUT) :: RECH !recharge to or from the water table when shallow [m] (diagnostic) - REAL, DIMENSION(1:60) , INTENT(INOUT) :: gecros1d ! gecros crop - ! output REAL , INTENT(OUT) :: Z0WRF !combined z0 sent to coupled model REAL , INTENT(OUT) :: FSA !total absorbed solar radiation (w/m2) @@ -674,8 +668,6 @@ SUBROUTINE NOAHMP_SFLX (parameters, & LOGICAL :: dveg_active ! flag to run dynamic vegetation LOGICAL :: crop_active ! flag to run crop model - REAL :: FB - ! INTENT (OUT) variables need to be assigned a value. These normally get assigned values ! only if DVEG == 2. nee = 0.0 @@ -781,8 +773,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & FSRG ,RSSUN ,RSSHA ,ALBSND ,ALBSNI, BGAP ,WGAP,TGV,TGB,& Q1 ,Q2V ,Q2B ,Q2E ,CHV ,CHB , & !out EMISSI ,PAH , & - SHG,SHC,SHB,EVG,EVB,GHV,GHB,IRG,IRC,IRB,TR,EVC,CHLEAF,CHUC,CHV2,CHB2,& - JULIAN, SWDOWN, PRCP, FB, GECROS1D ) + SHG,SHC,SHB,EVG,EVB,GHV,GHB,IRG,IRC,IRB,TR,EVC,CHLEAF,CHUC,CHV2,CHB2 ) !jref:end SICE(:) = MAX(0.0, SMC(:) - SH2O(:)) @@ -1532,8 +1523,7 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in T2MV ,T2MB ,FSRV , & FSRG ,RSSUN ,RSSHA ,ALBSND ,ALBSNI,BGAP ,WGAP,TGV,TGB,& Q1 ,Q2V ,Q2B ,Q2E ,CHV ,CHB, EMISSI,PAH ,& - SHG,SHC,SHB,EVG,EVB,GHV,GHB,IRG,IRC,IRB,TR,EVC,CHLEAF,CHUC,CHV2,CHB2, & - JULIAN, SWDOWN, PRCP, FB, GECROS1D ) + SHG,SHC,SHB,EVG,EVB,GHV,GHB,IRG,IRC,IRB,TR,EVC,CHLEAF,CHUC,CHV2,CHB2 ) !jref:end ! -------------------------------------------------------------------------------------------------- @@ -1783,9 +1773,6 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in REAL,INTENT(OUT) :: CHB2 !sensible heat conductance, canopy air to ZLVL air (m/s) REAL :: noahmpres - REAL, INTENT(IN) :: JULIAN, SWDOWN, PRCP, FB - REAL,DIMENSION(1:60),INTENT(INOUT) :: GECROS1D - !jref:end REAL, PARAMETER :: MPE = 1.E-6 @@ -2018,8 +2005,7 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in T2MV ,PSNSUN ,PSNSHA , & !out !jref:start QC ,QSFC ,PSFC , & !in - Q2V ,CHV2, CHLEAF, CHUC, & - SH2O,JULIAN, SWDOWN, PRCP, FB, FSR, GECROS1D) ! Gecros + Q2V ,CHV2, CHLEAF, CHUC ) !jref:end END IF @@ -3355,8 +3341,7 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & SHC ,EVG ,EVC ,TR ,GH , & !out T2MV ,PSNSUN ,PSNSHA , & !out QC ,QSFC ,PSFC , & !in - Q2V ,CAH2 ,CHLEAF ,CHUC, & !inout - SH2O,JULIAN, SWDOWN, PRCP, FB, FSR, GECROS1D) ! Gecros + Q2V ,CAH2 ,CHLEAF ,CHUC ) !inout ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3378,8 +3363,8 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & INTEGER, INTENT(IN) :: ISNOW !actual no. of snow layers INTEGER, INTENT(IN) :: VEGTYP !vegetation physiology type REAL, INTENT(IN) :: FVEG !greeness vegetation fraction (-) - REAL, INTENT(INOUT) :: SAV !solar rad absorbed by veg (w/m2) - REAL, INTENT(INOUT) :: SAG !solar rad absorbed by ground (w/m2) + REAL, INTENT(IN) :: SAV !solar rad absorbed by veg (w/m2) + REAL, INTENT(IN) :: SAG !solar rad absorbed by ground (w/m2) REAL, INTENT(IN) :: LWDN !atmospheric longwave radiation (w/m2) REAL, INTENT(IN) :: UR !wind speed at height zlvl (m/s) REAL, INTENT(IN) :: UU !wind speed in eastward dir (m/s) @@ -3560,15 +3545,6 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & INTEGER :: LITER !Last iteration - REAL, INTENT(IN) :: JULIAN, SWDOWN, PRCP, FB - REAL, INTENT(INOUT) :: FSR - REAL,DIMENSION(1:60), INTENT(INOUT) :: GECROS1D - REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2O !soil liquid water - - REAL :: ROOTD, WUL, WLL, Thickness, TLAIE, GLAIE, TLAI, GLAI, FRSU - INTEGER :: NROOT, J - - REAL :: T, TDC !Kelvin to degree Celsius with limit -50 to +50 character(len=80) :: message @@ -3716,56 +3692,6 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & RSSHA ,PSNSHA,ILOC ,JLOC ) !out END IF - ! Call Gecros - IF (opt_crop == 2) then - IF ((GECROS1D(41).GT.0).and.(GECROS1D(42).LT.0.)) then !Gecros - Thickness = 0. - NROOT = 0 - ROOTD = GECROS1D(33) - WUL = 0. - WLL = 0. - - DO J = 1,NSOIL - Thickness = Thickness + DZSNSO (J) - if (Thickness.lt.ROOTD/100.) then - NROOT = NROOT + 1 - endif - ENDDO - - NROOT = NROOT + 1 - NROOT = MAX(1,NROOT) - - Thickness = 0. - - DO J = 1,NROOT - Thickness = Thickness + DZSNSO (J) - if (Thickness.gt.ROOTD/100.) then - WUL = WUL + ((ROOTD/100.-Thickness+DZSNSO(J))*1000.*(SH2O(J)-parameters%SMCWLT(J))) - else - WUL = WUL + (DZSNSO(J)*1000.*(SH2O(J)-parameters%SMCWLT(J))) - endif - ENDDO - - DO J = 1,NSOIL - WLL = WLL + (DZSNSO(J)*1000.*(SH2O(J)-parameters%SMCWLT(J))) - ENDDO - WLL = WLL - WUL - - CALL gecros (JULIAN, DT, 1, RB, RAHC, RAHG+RSURF, FB, SNOWH , & !I - UR, SFCTMP, EAIR, SWDOWN, LWDN, PRCP, WUL, WLL , & !I - parameters%SMCWLT(1), parameters%DLEAF , & !I - GECROS1D , & !H - SAV, SAG, FSR, FRSU, RSSUN, RSSHA) !O - - GLAI = GECROS1D(49) - TLAI = GECROS1D(50) - - ! effective LAIs - TLAIE = MIN(6.,TLAI / FVEG) - GLAIE = MIN(6.,GLAI / FVEG) - ENDIF - ENDIF - END IF ! prepare for sensible heat flux above veg. @@ -3783,12 +3709,7 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & CAW = 1./RAWC CEW = FWET*VAIE/RB - IF (OPT_CROP /= 2) THEN - CTW = (1.-FWET)*(LAISUNE/(RB+RSSUN) + LAISHAE/(RB+RSSHA)) - ELSE - !RSSUN and RSSHA are in resistance per unit LAI in the Jarvis and Ball-Berry!. RSSUN and RSSHA of Gecros are in s/m - CTW = (1.-FWET)*(1./(RB/(FRSU*GLAIE)+RSSUN) + 1./(RB/((1.-FRSU)*GLAIE)+RSSHA)) !transpiration conductance leaf to canopy air - ENDIF + CTW = (1.-FWET)*(LAISUNE/(RB+RSSUN) + LAISHAE/(RB+RSSHA)) CGW = 1./(RAWG+RSURF) COND = CAW + CEW + CTW + CGW AEA = (EAIR*CAW + ESTG*CGW) / COND @@ -9044,53 +8965,6 @@ END SUBROUTINE PSN_CROP ! end subroutine bvocflux ! ================================================================================================== -!***************** SUBROUTINES FOR GECROS CROP SIMULATION *************** -!*----------------------------------------------------------------------* -!* SUBROUTINE EMERG * -!* Purpose: This subroutine calculates germination and emergence of * -!* the crop * -!* * -!* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) * -!* * -!* name type meaning units class * -!* ---- ---- ------- ----- ----- * -!* nowdate C12 Actual date dd.mm.yy I * -!* DT R4 Time step of integration s I * -!* DD R4 Drilling depth cm I * -!* TSOIL R4 Soil temperature in first layer K I * -!* TBEM R4 Temperature threshold oC I * -!* EMA R4 Intercept of function for emergence oC I * -!* EMB R4 Slope of function for emergence oC I * -!* TTEM R4 Cumulative temperature sum for emergence oC I/O * -!* EMERGENCE LOG Flag for emergence - I/O * -!* STATE_GECROS(41) = emerged yes/no * -!* STATE_GECROS(43) = TTEM * -!*----------------------------------------------------------------------* - -SUBROUTINE EMERG(DT, TSOIL, DD, TBEM, EMA, EMB, STATE_GECROS) - IMPLICIT NONE - REAL, INTENT(IN) :: DT, TSOIL, DD, TBEM, EMA, EMB - REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS - REAL :: EMTH, TINT - SAVE - - IF ((TSOIL-273.15).LT.TBEM) THEN - ELSE - STATE_GECROS(43) = STATE_GECROS(43) + (TSOIL-273.15-TBEM)/(86400./DT) - ENDIF - - EMTH = EMA + EMB*DD - - IF (STATE_GECROS(43).GT.EMTH) THEN - STATE_GECROS(41)=1. -! write(*,*) 'Crop emerged on ', nowdate -! read(*,*) - ELSE - STATE_GECROS(41)=-1. - ENDIF - - RETURN -END SUBROUTINE EMERG ! ********************************* end of carbon subroutines ***************************** ! ================================================================================================== @@ -9120,7 +8994,7 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc INTEGER, INTENT(IN) :: iopt_rsf !surface resistance (1->Sakaguchi/Zeng; 2->Seller; 3->mod Sellers; 4->1+snow) INTEGER, INTENT(IN) :: iopt_soil !soil parameters set-up option INTEGER, INTENT(IN) :: iopt_pedo !pedo-transfer function - INTEGER, INTENT(IN) :: iopt_crop !crop model option (0->none; 1->Liu et al.; 2->Gecros) + INTEGER, INTENT(IN) :: iopt_crop !crop model option (0->none; 1->Liu et al. 2016) ! ------------------------------------------------------------------------------------------------- diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index 1f3f74da07..9ac8a7d8bc 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -78,7 +78,6 @@ SUBROUTINE surface_driver( & & ,smcwtdxy ,rechxy ,deeprechxy,fdepthxy,areaxy ,rivercondxy, riverbedxy & & ,eqzwt ,pexpxy ,qrfxy ,qspringxy,qslatxy ,qrfsxy ,qspringsxy & & ,smoiseq ,wtddt ,stepwtd & - ,gecros_state & ! Optional gecros crop & ,opt_thcnd & ! Noah UA changes & ,ua_phys,flx4,fvb,fbur,fgsn & @@ -913,8 +912,6 @@ SUBROUTINE surface_driver( & SMCWTDXY, RECHXY, DEEPRECHXY, FDEPTHXY, AREAXY, RIVERCONDXY, RIVERBEDXY, & EQZWT, PEXPXY, QRFXY, QSPRINGXY, QSLATXY, QRFSXY, QSPRINGSXY - REAL, OPTIONAL, DIMENSION(ims:ime,60,jms:jme) :: gecros_state ! Optional gecros crop - INTEGER, INTENT(IN ):: OPT_THCND ! Noah UA changes LOGICAL, INTENT(IN) :: ua_phys @@ -3100,7 +3097,6 @@ SUBROUTINE surface_driver( & ZSNSOXY, SNICEXY, SNLIQXY, LFMASSXY, RTMASSXY, STMASSXY, & WOODXY, STBLCPXY, FASTCPXY, LAI, XSAIXY, TAUSSXY, & SMOISEQ, SMCWTDXY,DEEPRECHXY, RECHXY, GRAINXY, GDDXY,PGSXY, & ! IN/OUT Noah MP only - GECROS_STATE, & ! IN/OUT gecros model T2MVXY, T2MBXY, Q2MVXY, Q2MBXY, & TRADXY, NEEXY, GPPXY, NPPXY, FVEGXY, RUNSFXY, & RUNSBXY, ECANXY, EDIRXY, ETRANXY, FSAXY, FIRAXY, & diff --git a/run/README.namelist b/run/README.namelist index 397245df1e..0e9e194e63 100644 --- a/run/README.namelist +++ b/run/README.namelist @@ -1256,7 +1256,6 @@ Options for use with the Noah-MP Land Surface Model: geogrid must have been run with GEOGRID.TBL.ARW.noahmp, use with caution 0 -> No crop model, will run default dynamic vegetation 1 -> Liu, et al. 2016 - 2 -> Gecros (Genotype-by-Environment interaction on CROp growth Simulator) Yin and van Laar, 2005 / &fdda