diff --git a/BRAMS/src/ed2/edcp_driver.F90 b/BRAMS/src/ed2/edcp_driver.F90 index 6272d0b93..3bfbb1e29 100644 --- a/BRAMS/src/ed2/edcp_driver.F90 +++ b/BRAMS/src/ed2/edcp_driver.F90 @@ -66,7 +66,8 @@ subroutine ed_coup_driver() integer :: jd2 integer :: ierr integer :: igr - integer :: ping + integer :: ping + logical :: new_day real :: wtime1 real :: wtime2 real :: wtime_start ! wall time diff --git a/ED/src/driver/ed_driver.F90 b/ED/src/driver/ed_driver.F90 index c0bd69053..4533e9d81 100644 --- a/ED/src/driver/ed_driver.F90 +++ b/ED/src/driver/ed_driver.F90 @@ -182,7 +182,8 @@ subroutine ed_driver() !---------------------------------------------------------------------------------------! - if (trim(runtype) == 'HISTORY' ) then + select case (trim(runtype)) + case ('HISTORY') !------------------------------------------------------------------------------------! ! Initialize the model state as a replicate image of a previous state. ! !------------------------------------------------------------------------------------! @@ -207,7 +208,7 @@ subroutine ed_driver() if (nnodetot /= 1 ) call MPI_Barrier(MPI_COMM_WORLD,ierr) #endif !------------------------------------------------------------------------------------! - else + case default !------------------------------------------------------------------------------------! ! Initialize state properties of polygons/sites/patches/cohorts. ! @@ -215,7 +216,7 @@ subroutine ed_driver() if (mynum == nnodetot) write (unit=*,fmt='(a)') ' [+] Load_Ecosystem_State...' call load_ecosystem_state() !------------------------------------------------------------------------------------! - end if + end select !---------------------------------------------------------------------------------------! ! In case the runs is going to produce detailed output, we eliminate all patches ! @@ -291,11 +292,14 @@ subroutine ed_driver() ! Initialise some derived variables. Skip this in case the simulation is resuming ! ! from HISTORY. ! !---------------------------------------------------------------------------------------! - if (trim(runtype) /= 'HISTORY' ) then + select case (trim(runtype)) + case ('HISTORY') + continue + case default do ifm=1,ngrids call update_derived_props(edgrid_g(ifm)) end do - end if + end select !---------------------------------------------------------------------------------------! @@ -304,11 +308,14 @@ subroutine ed_driver() ! Initialise drought phenology. This should be done after the soil moisture has ! ! been set up. ! !---------------------------------------------------------------------------------------! - if (runtype /= 'HISTORY') then + select case (trim(runtype)) + case ('HISTORY') + continue + case default do ifm=1,ngrids call first_phenology(edgrid_g(ifm)) end do - end if + end select !---------------------------------------------------------------------------------------! diff --git a/ED/src/driver/ed_met_driver.f90 b/ED/src/driver/ed_met_driver.f90 index e6fb7449a..a5e6d6628 100644 --- a/ED/src/driver/ed_met_driver.f90 +++ b/ED/src/driver/ed_met_driver.f90 @@ -345,6 +345,8 @@ subroutine read_met_drivers_init type(edtype) , pointer :: cgrid character(len=str_len) :: infile integer :: igr + integer :: year_cyc + integer :: year_cyc_2 integer :: year_use integer :: iformat integer :: iv @@ -526,9 +528,18 @@ subroutine read_met_drivers_init !------------------------------------------------------------------------------------! ! We now retrieve the met driver year based on the stored sequence. ! !------------------------------------------------------------------------------------! - iyear = current_time%year-iyeara+1 - year_use = metyears(iyear) - + !----- If we need to recycle over years, find the appropriate year to apply. --------! + year_cyc = current_time%year + ncyc = metcycf - metcyc1 + 1 + !----- If we are after the last year... ---------------------------------------------! + do while(year_cyc > metcycf) + year_cyc = year_cyc - ncyc + end do + !----- If we are before the first year... -------------------------------------------! + do while(year_cyc < metcyc1) + year_cyc = year_cyc + ncyc + end do + !------------------------------------------------------------------------------------! gridloop: do igr = 1,ngrids @@ -536,7 +547,7 @@ subroutine read_met_drivers_init !----- Loop over the different file formats --------------------------------------! formloop: do iformat = 1, nformats - + !------------------------------------------------------------------------------! ! SPECIAL CASE FOR CO2: ! ! Usually we do not want to cycle CO2 but only the other meteorology. ! @@ -545,8 +556,30 @@ subroutine read_met_drivers_init !------------------------------------------------------------------------------! not_cycle_co2 = (met_nv(iformat) == 1 .and. trim(met_vars(iformat,1)) == 'co2') if (not_cycle_co2) then - year_use = current_time%year - endif + !---------------------------------------------------------------------------! + ! Make sure that the special CO2 file exists. If not, then use the ! + ! default met cycle years. ! + !---------------------------------------------------------------------------! + write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)),current_time%year & + ,mname(current_time%month),'.h5' + inquire(file=trim(infile),exist=exans) + if (exans) then + !----- Allow CO2 outside met cycle. -------------------------------------! + year_use = current_time%year + !------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. -----------! + year_use = year_cyc + !------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. --------------! + year_use = year_cyc + !---------------------------------------------------------------------------! + end if + !------------------------------------------------------------------------------! + !----- Create the file name and check whether it exists. ----------------------! write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)), year_use & ,mname(current_time%month),'.h5' @@ -555,10 +588,17 @@ subroutine read_met_drivers_init call shdf5_open_f(trim(infile),'R') else write (unit=*,fmt='(a)' ) '------------------------------' - write (unit=*,fmt='(a,1x,i12)') ' - METCYC1 =',metcyc1 - write (unit=*,fmt='(a,1x,i12)') ' - METCYCF =',metcycf - write (unit=*,fmt='(a,1x,i12)') ' - IYEAR =',iyear - write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE =',year_use + write (unit=*,fmt='(a,1x,a)' ) ' - MET_VARS =',trim(met_vars(iformat,1)) + write (unit=*,fmt='(a,1x,l1)' ) ' - CYCLE_CO2 =',.not. not_cycle_co2 + write (unit=*,fmt='(a,1x,i12)') ' - METCYC1 =',metcyc1 + write (unit=*,fmt='(a,1x,i12)') ' - METCYCF =',metcycf + write (unit=*,fmt='(a,1x,i12)') ' - NCYC =',ncyc + write (unit=*,fmt='(a,1x,i12)') ' - NYEARS =',nyears + write (unit=*,fmt='(a,1x,i12)') ' - IYEAR =',iyear + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CURR =',current_time%month + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CURR =',current_time%year + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC =',year_cyc + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE =',year_use write (unit=*,fmt='(a)' ) '------------------------------' call fatal_error('Cannot open met driver input file '//trim(infile)//'!' & ,'read_met_drivers_init','ed_met_driver.f90') @@ -573,7 +613,7 @@ subroutine read_met_drivers_init !----- Loop over variables. and read the data. --------------------------------! do iv = 1, met_nv(iformat) offset = 0 - call read_ol_file(infile,iformat, iv, mname(current_time%month) & + call read_ol_file(infile,iformat, iv, mname(current_time%month) & ,current_time%year, offset, cgrid) end do @@ -584,27 +624,61 @@ subroutine read_met_drivers_init !------------------------------------------------------------------------------! ! For all interpolated variables, we also need the next time. ! !------------------------------------------------------------------------------! - !------ Find next month and year ----------------------------------------------! - m2 = current_time%month + 1 - + !------------------------------------------------------------------------------! - ! If this takes us into the next year, take the next year in sequence and ! - ! reset month to January. ! + ! Find next month and year. If this takes us into the next year, increment ! + ! year and reset month to January. ! !------------------------------------------------------------------------------! - if (m2 == 13) then - m2 = 1 - y2 = current_time%year + 1 - else - !----- Otherwise, use the same year. ---------------------------------------! - y2 = current_time%year - end if - iyear = y2 - iyeara + 1 - year_use_2 = metyears(iyear) - - ! Again consider the special case of not cycling co2 + m2 = current_time%month + 1 + select case (m2) + case (13) + m2 = 1 + y2 = current_time%year + 1 + year_cyc_2 = y2 + + !----- If we are now after the last year... --------------------------------! + do while(year_cyc_2 > metcycf) + year_cyc_2 = year_cyc_2 - ncyc + end do + !---------------------------------------------------------------------------! + + !----- If we are now before the first year... ------------------------------! + do while(year_cyc_2 < metcyc1) + year_cyc_2 = year_cyc_2 + ncyc + end do + !---------------------------------------------------------------------------! + case default + !---- Same year as the previous month. -------------------------------------! + y2 = current_time%year + year_cyc_2 = year_cyc + !---------------------------------------------------------------------------! + end select + + !------------------------------------------------------------------------------# + ! Again consider the special case of not cycling co2 + !------------------------------------------------------------------------------# if (not_cycle_co2) then - year_use_2 = y2 - endif + !---------------------------------------------------------------------------! + ! Make sure that the special CO2 file exists. If not, then use the ! + ! default met cycle years. ! + !---------------------------------------------------------------------------! + write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)),y2,mname(m2),'.h5' + inquire(file=trim(infile),exist=exans) + if (exans) then + !----- Allow CO2 outside met cycle. -------------------------------------! + year_use_2 = y2 + !------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. -----------! + year_use_2 = year_cyc_2 + !------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. --------------! + year_use_2 = year_cyc_2 + !---------------------------------------------------------------------------! + end if !----- Now, open the file once. -----------------------------------------------! write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)), year_use_2 & ,mname(m2),'.h5' @@ -612,10 +686,26 @@ subroutine read_met_drivers_init if (exans) then call shdf5_open_f(trim(infile),'R') else + write (unit=*,fmt='(a)' ) '------------------------------' + write (unit=*,fmt='(a,1x,a)' ) ' - MET_VARS =',trim(met_vars(iformat,1)) + write (unit=*,fmt='(a,1x,l1)' ) ' - CYCLE_CO2 =',.not. not_cycle_co2 + write (unit=*,fmt='(a,1x,i12)') ' - METCYC1 =',metcyc1 + write (unit=*,fmt='(a,1x,i12)') ' - METCYCF =',metcycf + write (unit=*,fmt='(a,1x,i12)') ' - NCYC =',ncyc + write (unit=*,fmt='(a,1x,i12)') ' - NYEARS =',nyears + write (unit=*,fmt='(a,1x,i12)') ' - IYEAR =',iyear + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CURR =',current_time%month + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CURR =',current_time%year + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC =',year_cyc + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE =',year_use + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CYC_2 =',m2 + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC_2 =',year_cyc_2 + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE_2 =',year_use_2 + write (unit=*,fmt='(a)' ) '------------------------------' call fatal_error('Cannot open met driver input file '//trim(infile)//'!' & ,'read_met_drivers_init','ed_met_driver.f90') end if - + !----- Loop over variables. ---------------------------------------------------! varloop: do iv = 1, met_nv(iformat) @@ -674,6 +764,8 @@ subroutine read_met_drivers type(edtype) , pointer :: cgrid character(len=str_len) :: infile integer :: igr + integer :: year_cyc + integer :: year_cyc_2 integer :: year_use integer :: ncyc integer :: iformat @@ -694,17 +786,17 @@ subroutine read_met_drivers !----- If we need to recycle over years, find the appropriate year to apply. --------! - year_use = current_time%year - ncyc = metcycf - metcyc1 + 1 + year_cyc = current_time%year + ncyc = metcycf - metcyc1 + 1 !----- If we are after the last year... ---------------------------------------------! - do while(year_use > metcycf) - year_use = year_use - ncyc + do while(year_cyc > metcycf) + year_cyc = year_cyc - ncyc end do !----- If we are before the first year... -------------------------------------------! - do while(year_use < metcyc1) - year_use = year_use + ncyc + do while(year_cyc < metcyc1) + year_cyc = year_cyc + ncyc end do gridloop: do igr=1,ngrids @@ -713,7 +805,7 @@ subroutine read_met_drivers !----- Loop over the different file formats --------------------------------------! formloop: do iformat = 1, nformats - + !------------------------------------------------------------------------------! ! SPECIAL CASE FOR CO2: ! ! Usually we do not want to cycle CO2 but only the other meteorology. ! @@ -722,8 +814,32 @@ subroutine read_met_drivers !------------------------------------------------------------------------------! not_cycle_co2 = (met_nv(iformat) == 1 .and. trim(met_vars(iformat,1)) == 'co2') if (not_cycle_co2) then - year_use = current_time%year - endif + !---------------------------------------------------------------------------! + ! Make sure that the special CO2 file exists. If not, then use the ! + ! default met cycle years. ! + !---------------------------------------------------------------------------! + write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)),current_time%year & + ,mname(current_time%month),'.h5' + inquire(file=trim(infile),exist=exans) + if (exans) then + !----- Allow CO2 outside met cycle. -------------------------------------! + year_use = current_time%year + !------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. -----------! + year_use = year_cyc + !------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. --------------! + year_use = year_cyc + !---------------------------------------------------------------------------! + end if + !------------------------------------------------------------------------------! + + + !----- Create the file name and check whether it exists. ----------------------! write(infile,'(a,i4.4,a,a)')trim(met_names(iformat)), year_use, & mname(current_time%month),'.h5' @@ -732,6 +848,17 @@ subroutine read_met_drivers if(exans)then call shdf5_open_f(trim(infile),'R') else + write (unit=*,fmt='(a)' ) '------------------------------' + write (unit=*,fmt='(a,1x,a)' ) ' - MET_VARS =',trim(met_vars(iformat,1)) + write (unit=*,fmt='(a,1x,l1)' ) ' - CYCLE_CO2 =',.not. not_cycle_co2 + write (unit=*,fmt='(a,1x,i12)') ' - METCYC1 =',metcyc1 + write (unit=*,fmt='(a,1x,i12)') ' - METCYCF =',metcycf + write (unit=*,fmt='(a,1x,i12)') ' - NCYC =',ncyc + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CURR =',current_time%month + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CURR =',current_time%year + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC =',year_cyc + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE =',year_use + write (unit=*,fmt='(a)' ) '------------------------------' call fatal_error('Cannot open met driver input file '//trim(infile)//'!' & ,'read_met_drivers','ed_met_driver.f90') end if @@ -766,37 +893,63 @@ subroutine read_met_drivers !------------------------------------------------------------------------------! ! For all interpolated variables, we also need the next time. ! !------------------------------------------------------------------------------! - !------ Find next month and year ----------------------------------------------! - m2 = current_time%month + 1 - y2 = current_time%year - year_use_2 = year_use - - + !------------------------------------------------------------------------------! - ! If this takes us into the next year, increment year and reset month to ! - ! January. ! + ! Find next month and year. If this takes us into the next year, increment ! + ! year and reset month to January. ! !------------------------------------------------------------------------------! - if(m2 == 13)then + m2 = current_time%month + 1 + select case (m2) + case (13) m2 = 1 y2 = current_time%year + 1 - year_use_2 = y2 + year_cyc_2 = y2 !----- If we are now after the last year... --------------------------------! - do while(year_use_2 > metcycf) - year_use_2 = year_use_2 - ncyc + do while(year_cyc_2 > metcycf) + year_cyc_2 = year_cyc_2 - ncyc end do - + !---------------------------------------------------------------------------! + !----- If we are now before the first year... ------------------------------! - do while(year_use_2 < metcyc1) - year_use_2 = year_use_2 + ncyc + do while(year_cyc_2 < metcyc1) + year_cyc_2 = year_cyc_2 + ncyc end do - end if + !---------------------------------------------------------------------------! + case default + !---- Same year as the previous month. -------------------------------------! + y2 = current_time%year + year_cyc_2 = year_cyc + !---------------------------------------------------------------------------! + end select - ! Again, consider the special case for not_cycle_co2 + !---- Again, consider the special case for not_cycle_co2. ---------------------! if (not_cycle_co2) then - year_use_2 = y2 - endif - + !---------------------------------------------------------------------------! + ! Make sure that the special CO2 file exists. If not, then use the ! + ! default met cycle years. ! + !---------------------------------------------------------------------------! + write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)),y2 & + ,mname(current_time%month),'.h5' + inquire(file=trim(infile),exist=exans) + if (exans) then + !----- Allow CO2 outside met cycle. -------------------------------------! + year_use_2 = y2 + !------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. -----------! + year_use_2 = year_cyc_2 + !------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------! + else + !----- Typical case. Pool data from the meteorological cycle. --------------! + year_use_2 = year_cyc_2 + !---------------------------------------------------------------------------! + end if + !------------------------------------------------------------------------------! + + !----- Now, open the file once. -----------------------------------------------! write(infile,fmt='(a,i4.4,a,a)') trim(met_names(iformat)), year_use_2 & ,mname(m2),'.h5' @@ -804,6 +957,20 @@ subroutine read_met_drivers if (exans) then call shdf5_open_f(trim(infile),'R') else + write (unit=*,fmt='(a)' ) '------------------------------' + write (unit=*,fmt='(a,1x,a)' ) ' - MET_VARS =',trim(met_vars(iformat,1)) + write (unit=*,fmt='(a,1x,l1)' ) ' - CYCLE_CO2 =',.not. not_cycle_co2 + write (unit=*,fmt='(a,1x,i12)') ' - METCYC1 =',metcyc1 + write (unit=*,fmt='(a,1x,i12)') ' - METCYCF =',metcycf + write (unit=*,fmt='(a,1x,i12)') ' - NCYC =',ncyc + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CURR =',current_time%month + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CURR =',current_time%year + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC =',year_cyc + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE =',year_use + write (unit=*,fmt='(a,1x,i12)') ' - MONTH_CYC_2 =',m2 + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_CYC_2 =',year_cyc_2 + write (unit=*,fmt='(a,1x,i12)') ' - YEAR_USE_2 =',year_use_2 + write (unit=*,fmt='(a)' ) '------------------------------' call fatal_error ('Cannot open met driver input file '//trim(infile)//'!' & ,'read_met_drivers','ed_met_driver.f90') end if diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index f1dc10627..ea34fe56f 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -1563,7 +1563,7 @@ subroutine site_disturbance_rates(year, cgrid) !---------------------------------------------------------------------------! case default !------ Read anthropogenic disturbance from external data set. -------------! - if (clutime%landuse(12) < 0 .or. clutime%landuse(14) < 0) then + if (clutime%landuse(12) < 0. .or. clutime%landuse(14) < 0.) then find_target = .true. cpoly%primary_harvest_target (isi) = 0. cpoly%secondary_harvest_target(isi) = 0. @@ -2265,6 +2265,9 @@ subroutine increment_patch_vars(csite,np,cp,area_fac,cb_enthalpy,can_exner,cb_ma csite%fmean_sfcw_mass (np) = csite%fmean_sfcw_mass (np) & + csite%fmean_sfcw_mass (cp) & * area_fac + csite%fmean_snowfac (np) = csite%fmean_snowfac (np) & + + csite%fmean_snowfac (cp) & + * area_fac csite%fmean_rshort_gnd (np) = csite%fmean_rshort_gnd (np) & + csite%fmean_rshort_gnd (cp) & * area_fac @@ -2487,6 +2490,9 @@ subroutine increment_patch_vars(csite,np,cp,area_fac,cb_enthalpy,can_exner,cb_ma csite%dmean_sfcw_fliq ( np) = csite%dmean_sfcw_fliq ( np) & + csite%dmean_sfcw_fliq ( cp) & * area_fac + csite%dmean_snowfac ( np) = csite%dmean_snowfac ( np) & + + csite%dmean_snowfac ( cp) & + * area_fac csite%dmean_rshort_gnd ( np) = csite%dmean_rshort_gnd ( np) & + csite%dmean_rshort_gnd ( cp) & * area_fac @@ -2770,6 +2776,9 @@ subroutine increment_patch_vars(csite,np,cp,area_fac,cb_enthalpy,can_exner,cb_ma csite%mmean_sfcw_fliq ( np) = csite%mmean_sfcw_fliq ( np) & + csite%mmean_sfcw_fliq ( cp) & * area_fac + csite%mmean_snowfac ( np) = csite%mmean_snowfac ( np) & + + csite%mmean_snowfac ( cp) & + * area_fac csite%mmean_rshort_gnd ( np) = csite%mmean_rshort_gnd ( np) & + csite%mmean_rshort_gnd ( cp) & * area_fac @@ -3065,6 +3074,9 @@ subroutine increment_patch_vars(csite,np,cp,area_fac,cb_enthalpy,can_exner,cb_ma csite%qmean_sfcw_fliq ( :,np) = csite%qmean_sfcw_fliq ( :,np) & + csite%qmean_sfcw_fliq ( :,cp) & * area_fac + csite%qmean_snowfac ( :,np) = csite%qmean_snowfac ( :,np) & + + csite%qmean_snowfac ( :,cp) & + * area_fac csite%qmean_soil_energy (:,:,np) = csite%qmean_soil_energy (:,:,np) & + csite%qmean_soil_energy (:,:,cp) & * area_fac diff --git a/ED/src/dynamics/euler_driver.f90 b/ED/src/dynamics/euler_driver.f90 index 416104399..a1bcbb0b5 100644 --- a/ED/src/dynamics/euler_driver.f90 +++ b/ED/src/dynamics/euler_driver.f90 @@ -97,7 +97,7 @@ subroutine euler_timestep(cgrid) ! Update the monthly rainfall. ! !------------------------------------------------------------------------------! imon = current_time%month - cpoly%avg_monthly_pcpg(imon,isi) = cpoly%avg_monthly_pcpg(imon,isi) & + cpoly%avg_monthly_accp(imon,isi) = cpoly%avg_monthly_accp(imon,isi) & + cmet%pcpg * dtlsm !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/farq_katul.f90 b/ED/src/dynamics/farq_katul.f90 index 529d34a04..0e184f0b1 100644 --- a/ED/src/dynamics/farq_katul.f90 +++ b/ED/src/dynamics/farq_katul.f90 @@ -247,9 +247,6 @@ subroutine katul_lphys(ib,can_prss,can_rhos,can_shv,can_co2,ipft,leaf_par,leaf_t !update photosynthetic parameters with water stress select case (h2o_plant_lim) - case (0,1,2,3) - ! use fsw to account for water stress in photosyn_driv - water_stress_factor = 1. case (4) ! leaf water potential will influence stomata optimization ! at two different scales @@ -268,6 +265,9 @@ subroutine katul_lphys(ib,can_prss,can_rhos,can_shv,can_co2,ipft,leaf_par,leaf_t 1. / (1. + & 0.1 * (leaf_psi / leaf_psi_tlp(ipft)) ** 6.0))) lambda8 = lambda8 * dble(exp(stoma_beta(ipft) * dmax_leaf_psi)) + case default + ! use fsw to account for water stress in photosyn_driv + water_stress_factor = 1. end select !thispft(ib)%vm0 = thispft(ib)%vm0 * dble(water_stress_factor) @@ -848,8 +848,8 @@ subroutine photosynthesis_stomata_solver8(ib,gsc,limit_case, real(kind=8) :: k1,k2 !! Variable used in photosynthesis equation real(kind=8) :: a,b,c !! Coefficients of the quadratic equation to solve ci real(kind=8) :: rad !! sqrt(b2-4ac) - real(kind=8) :: dbdg,dcdg !! derivatives of b,c wrt. gsc - + real(kind=8) :: dbdg !! derivatives of b wrt. gsc + real(kind=8) :: dcdg !! derivatives of c wrt. gsc !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/fire.f90 b/ED/src/dynamics/fire.f90 index 0df3bcbf3..b93c9703e 100644 --- a/ED/src/dynamics/fire.f90 +++ b/ED/src/dynamics/fire.f90 @@ -62,7 +62,7 @@ subroutine fire_frequency(cgrid) real :: fuel real :: ignition_rate real :: mean_fire_intensity - real :: sum_pcpg + real :: sum_accp logical :: people_around !------------------------------------------------------------------------------------! @@ -97,8 +97,8 @@ subroutine fire_frequency(cgrid) ! Find the total rainfall of the past year and reset the counter for this ! ! month. ! !------------------------------------------------------------------------------! - sum_pcpg = sum(cpoly%avg_monthly_pcpg(:,isi)) - cpoly%avg_monthly_pcpg(imon,isi) = 0. + sum_accp = sum(cpoly%avg_monthly_accp(:,isi)) + cpoly%avg_monthly_accp(imon,isi) = 0. !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 320513cbf..39a27c12a 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -142,6 +142,10 @@ subroutine find_lambda_harvest(cpoly,isi,onsp,lambda_harvest) ! croplands and pastures. ! !---------------------------------------------------------------------------------! select case(ilu) + case (1,8) + !---- Pasture or cropland. Do nothing. ----------------------------------------! + continue + !------------------------------------------------------------------------------! case (2:7) hcoh_loop: do ico=1,cpatch%ncohorts ipft = cpatch%pft(ico) diff --git a/ED/src/dynamics/heun_driver.f90 b/ED/src/dynamics/heun_driver.f90 index c3c9631b6..d5aae0d2d 100644 --- a/ED/src/dynamics/heun_driver.f90 +++ b/ED/src/dynamics/heun_driver.f90 @@ -98,7 +98,7 @@ subroutine heun_timestep(cgrid) ! Update the monthly rainfall. ! !------------------------------------------------------------------------------! imon = current_time%month - cpoly%avg_monthly_pcpg(imon,isi) = cpoly%avg_monthly_pcpg(imon,isi) & + cpoly%avg_monthly_accp(imon,isi) = cpoly%avg_monthly_accp(imon,isi) & + cmet%pcpg * dtlsm !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/hybrid_driver.f90 b/ED/src/dynamics/hybrid_driver.f90 index 7aaf1d59f..54689f0d7 100644 --- a/ED/src/dynamics/hybrid_driver.f90 +++ b/ED/src/dynamics/hybrid_driver.f90 @@ -107,7 +107,7 @@ subroutine hybrid_timestep(cgrid) ! Update the monthly rainfall. ! !---------------------------------------------------------------------! imon = current_time%month - cpoly%avg_monthly_pcpg(imon,isi) = cpoly%avg_monthly_pcpg(imon,isi) & + cpoly%avg_monthly_accp(imon,isi) = cpoly%avg_monthly_accp(imon,isi) & + cmet%pcpg * dtlsm !---------------------------------------------------------------------! diff --git a/ED/src/dynamics/reproduction.f90 b/ED/src/dynamics/reproduction.f90 index 170fec5a8..52ab6a4d6 100644 --- a/ED/src/dynamics/reproduction.f90 +++ b/ED/src/dynamics/reproduction.f90 @@ -179,6 +179,7 @@ subroutine reproduction_driver(cgrid,month,veget_dyn_on) !----- The big loops start here. -------------------------------------------------! polyloop: do ipy = 1,cgrid%npolygons + cpoly => cgrid%polygon(ipy) !------------------------------------------------------------------------------! ! Check whether this is late spring/early summer. This is needed for ! @@ -187,8 +188,8 @@ subroutine reproduction_driver(cgrid,month,veget_dyn_on) !------------------------------------------------------------------------------! late_spring = (cgrid%lat(ipy) >= 0.0 .and. month == 6) .or. & (cgrid%lat(ipy) < 0.0 .and. month == 12) + !------------------------------------------------------------------------------! - cpoly => cgrid%polygon(ipy) siteloop_sort: do isi = 1,cpoly%nsites csite => cpoly%site(isi) @@ -1286,8 +1287,8 @@ subroutine seed_dispersal(cpoly,late_spring) ! of them will land again in this patch, and we correct for this further ! ! down. ! !------------------------------------------------------------------------! - csite%cbudget_seedrain(donpa) = csite%cbudget_seedrain(donpa) & - - bseed_maygo * frqsumi + donsite%cbudget_seedrain(donpa) = donsite%cbudget_seedrain(donpa) & + - bseed_maygo * frqsumi !------------------------------------------------------------------------! @@ -1314,7 +1315,7 @@ subroutine seed_dispersal(cpoly,late_spring) ! (4) RPY = DPA * AD * AR (1->3) ! ! (5) RPA = DPA * AD (4->2, regardless of the patch) ! !------------------------------------------------------------------------! - bseed_xpatch = bseed_maygo * csite%area(donpa) * cpoly%area(donsi) + bseed_xpatch = bseed_maygo * donsite%area(donpa) * cpoly%area(donsi) !------------------------------------------------------------------------! @@ -1343,8 +1344,8 @@ subroutine seed_dispersal(cpoly,late_spring) ! subtracted all the non-local dispersal outside the receptor site ! ! loop. ! !------------------------------------------------------------------! - csite%cbudget_seedrain(recpa) = csite%cbudget_seedrain(recpa) & - + bseed_xpatch * frqsumi + recsite%cbudget_seedrain(recpa) = recsite%cbudget_seedrain(recpa) & + + bseed_xpatch * frqsumi !------------------------------------------------------------------! diff --git a/ED/src/dynamics/rk4_copy_patch.f90 b/ED/src/dynamics/rk4_copy_patch.f90 index 9fbeff87e..fc7556a37 100644 --- a/ED/src/dynamics/rk4_copy_patch.f90 +++ b/ED/src/dynamics/rk4_copy_patch.f90 @@ -2092,6 +2092,8 @@ subroutine initp2modelp(hdid,initp,csite,ipa,nighttime,wbudget_loss2atm,ebudget_ + csite%ground_shv (ipa) * dtlsm_o_frqsum csite%fmean_can_ggnd (ipa) = csite%fmean_can_ggnd (ipa) & + csite%ggnet (ipa) * dtlsm_o_frqsum + csite%fmean_snowfac (ipa) = csite%fmean_snowfac (ipa) & + + csite%snowfac (ipa) * dtlsm_o_frqsum !------------------------------------------------------------------------------------! ! Snow/pounding layers. We keep track of the total, not individual layers. ! ! Energy will be integrated as an extensive variable, we will convert it by the ! diff --git a/ED/src/dynamics/rk4_driver.F90 b/ED/src/dynamics/rk4_driver.F90 index ae758614c..38ba31935 100644 --- a/ED/src/dynamics/rk4_driver.F90 +++ b/ED/src/dynamics/rk4_driver.F90 @@ -112,7 +112,7 @@ subroutine rk4_timestep(cgrid) ! Update the monthly rainfall. ! !------------------------------------------------------------------------------! imon = current_time%month - cpoly%avg_monthly_pcpg(imon,isi) = cpoly%avg_monthly_pcpg(imon,isi) & + cpoly%avg_monthly_accp(imon,isi) = cpoly%avg_monthly_accp(imon,isi) & + cmet%pcpg * dtlsm !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/vegetation_dynamics.f90 b/ED/src/dynamics/vegetation_dynamics.f90 index 645ae24c9..c30e6dcd7 100644 --- a/ED/src/dynamics/vegetation_dynamics.f90 +++ b/ED/src/dynamics/vegetation_dynamics.f90 @@ -34,7 +34,7 @@ subroutine veg_dynamics_driver(new_month,new_year,gr_tfact0,veget_dyn_on) , yr_day ! ! intent(in) use mem_polygons , only : maxpatch ! ! intent(in) use average_utils , only : normalize_ed_today_vars & ! sub-routine - , normalize_ed_todaynpp_vars & ! sub-routine + , copy_today_to_dmean_vars & ! sub-routine , zero_ed_today_vars ! ! sub-routine use canopy_radiation_coms, only : ihrzrad ! ! intent(in) use hrzshade_utils , only : split_hrzshade & ! sub-routine @@ -129,7 +129,7 @@ subroutine veg_dynamics_driver(new_month,new_year,gr_tfact0,veget_dyn_on) !------ update dmean and mmean values for NPP allocation terms ------------------! - call normalize_ed_todayNPP_vars(cgrid) + call copy_today_to_dmean_vars(cgrid) !---------------------------------------------------------------------------------! diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index 98b14d90b..f991239ae 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -525,7 +525,7 @@ subroutine init_decomp_params() rh0 = 0.700 ! 0.701 ! 0.425 rh_q10 = 1.500 ! 1.500 ! 1.893 rh_p_smoist = 1.600 ! 0.836 ! 0.606 - rh_p_oxygen = 0.600 ! 0.404 ! 0.164 + rh_p_oxygen = 0.450 ! 0.404 ! 0.164 !---------------------------------------------------------------------------------------! diff --git a/ED/src/init/ed_type_init.f90 b/ED/src/init/ed_type_init.f90 index 456c93fd4..ad93868f8 100644 --- a/ED/src/init/ed_type_init.f90 +++ b/ED/src/init/ed_type_init.f90 @@ -1128,6 +1128,7 @@ subroutine init_ed_patch_vars(csite,ipaa,ipaz,lsl) csite%fmean_sfcw_mass (ipaa:ipaz) = 0.0 csite%fmean_sfcw_temp (ipaa:ipaz) = 0.0 csite%fmean_sfcw_fliq (ipaa:ipaz) = 0.0 + csite%fmean_snowfac (ipaa:ipaz) = 0.0 csite%fmean_rshort_gnd (ipaa:ipaz) = 0.0 csite%fmean_par_gnd (ipaa:ipaz) = 0.0 csite%fmean_rlong_gnd (ipaa:ipaz) = 0.0 @@ -1210,6 +1211,7 @@ subroutine init_ed_patch_vars(csite,ipaa,ipaz,lsl) csite%dmean_sfcw_mass (ipaa:ipaz) = 0.0 csite%dmean_sfcw_temp (ipaa:ipaz) = 0.0 csite%dmean_sfcw_fliq (ipaa:ipaz) = 0.0 + csite%dmean_snowfac (ipaa:ipaz) = 0.0 csite%dmean_rshort_gnd (ipaa:ipaz) = 0.0 csite%dmean_par_gnd (ipaa:ipaz) = 0.0 csite%dmean_rlong_gnd (ipaa:ipaz) = 0.0 @@ -1286,6 +1288,7 @@ subroutine init_ed_patch_vars(csite,ipaa,ipaz,lsl) csite%mmean_sfcw_mass (ipaa:ipaz) = 0.0 csite%mmean_sfcw_temp (ipaa:ipaz) = 0.0 csite%mmean_sfcw_fliq (ipaa:ipaz) = 0.0 + csite%mmean_snowfac (ipaa:ipaz) = 0.0 csite%mmean_rshort_gnd (ipaa:ipaz) = 0.0 csite%mmean_par_gnd (ipaa:ipaz) = 0.0 csite%mmean_rlong_gnd (ipaa:ipaz) = 0.0 @@ -1409,6 +1412,7 @@ subroutine init_ed_patch_vars(csite,ipaa,ipaz,lsl) csite%qmean_sfcw_mass (:,ipaa:ipaz) = 0.0 csite%qmean_sfcw_temp (:,ipaa:ipaz) = 0.0 csite%qmean_sfcw_fliq (:,ipaa:ipaz) = 0.0 + csite%qmean_snowfac (:,ipaa:ipaz) = 0.0 csite%qmean_rshort_gnd (:,ipaa:ipaz) = 0.0 csite%qmean_par_gnd (:,ipaa:ipaz) = 0.0 csite%qmean_rlong_gnd (:,ipaa:ipaz) = 0.0 @@ -1522,6 +1526,7 @@ subroutine init_ed_site_vars(cpoly) , writing_eorq & ! intent(in) , writing_dcyc & ! intent(in) , economics_scheme ! ! intent(in) + use consts_coms , only : huge_num ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! type(polygontype), target :: cpoly @@ -1674,7 +1679,7 @@ subroutine init_ed_site_vars(cpoly) ! Initialise the minimum monthly temperature with a very large value, this is ! ! going to be reduced as the canopy temperature is updated. ! !------------------------------------------------------------------------------------! - cpoly%min_monthly_temp(:) = huge(1.) + cpoly%min_monthly_temp(:) = huge_num !------------------------------------------------------------------------------------! @@ -1684,7 +1689,7 @@ subroutine init_ed_site_vars(cpoly) ! by actual rainfall after 12 months. In the future we may initialise with climato- ! ! logical rainfall. ! !------------------------------------------------------------------------------------! - cpoly%avg_monthly_pcpg(:,:) = 500. + cpoly%avg_monthly_accp(:,:) = 500. !------------------------------------------------------------------------------------! @@ -2084,6 +2089,7 @@ subroutine init_ed_poly_vars(cgrid) cgrid%fmean_sfcw_mass (ipy) = 0.0 cgrid%fmean_sfcw_temp (ipy) = 0.0 cgrid%fmean_sfcw_fliq (ipy) = 0.0 + cgrid%fmean_snowfac (ipy) = 0.0 cgrid%fmean_rshort_gnd (ipy) = 0.0 cgrid%fmean_par_gnd (ipy) = 0.0 cgrid%fmean_rlong_gnd (ipy) = 0.0 @@ -2257,6 +2263,7 @@ subroutine init_ed_poly_vars(cgrid) cgrid%dmean_sfcw_mass (ipy) = 0.0 cgrid%dmean_sfcw_temp (ipy) = 0.0 cgrid%dmean_sfcw_fliq (ipy) = 0.0 + cgrid%dmean_snowfac (ipy) = 0.0 cgrid%dmean_rshort_gnd (ipy) = 0.0 cgrid%dmean_par_gnd (ipy) = 0.0 cgrid%dmean_rlong_gnd (ipy) = 0.0 @@ -2414,6 +2421,7 @@ subroutine init_ed_poly_vars(cgrid) cgrid%mmean_sfcw_mass (ipy) = 0.0 cgrid%mmean_sfcw_temp (ipy) = 0.0 cgrid%mmean_sfcw_fliq (ipy) = 0.0 + cgrid%mmean_snowfac (ipy) = 0.0 cgrid%mmean_rshort_gnd (ipy) = 0.0 cgrid%mmean_par_gnd (ipy) = 0.0 cgrid%mmean_rlong_gnd (ipy) = 0.0 @@ -2656,6 +2664,7 @@ subroutine init_ed_poly_vars(cgrid) cgrid%qmean_sfcw_mass (:,ipy) = 0.0 cgrid%qmean_sfcw_temp (:,ipy) = 0.0 cgrid%qmean_sfcw_fliq (:,ipy) = 0.0 + cgrid%qmean_snowfac (:,ipy) = 0.0 cgrid%qmean_rshort_gnd (:,ipy) = 0.0 cgrid%qmean_par_gnd (:,ipy) = 0.0 cgrid%qmean_rlong_gnd (:,ipy) = 0.0 diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index 77de83090..b2839475d 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -112,7 +112,7 @@ subroutine read_landuse_matrix write (unit=*,fmt='(a)') ' so MAX_LU_YEARS >= IYEARZ-IYEARA+1, then recompile' write (unit=*,fmt='(a)') ' your model.' call fatal_error ('Simulation is too long for anthropogenic disturbance.' & - ,'landuse_init','landuse_init.f90') + ,'read_landuse_matrix','landuse_init.f90') end if !------------------------------------------------------------------------------------! @@ -301,7 +301,7 @@ subroutine read_landuse_matrix write (unit=*,fmt='(a,1x,i6)') ' - Nharvest:',nharvest write (unit=*,fmt='(a)') '------------------------------------' call fatal_error('Land use area is zero, it doesn''t make any sense!' & - ,'landuse_init','landuse_init.f90') + ,'read_landuse_matrix','landuse_init.f90') else lu_area_i = 1. / lu_area end if diff --git a/ED/src/io/average_utils.f90 b/ED/src/io/average_utils.f90 index 152c09294..722b61051 100644 --- a/ED/src/io/average_utils.f90 +++ b/ED/src/io/average_utils.f90 @@ -493,6 +493,9 @@ subroutine aggregate_polygon_fmean(cgrid) cgrid%fmean_sfcw_mass (ipy) = cgrid%fmean_sfcw_mass (ipy) & + csite%fmean_sfcw_mass (ipa) & * patch_wgt + cgrid%fmean_snowfac (ipy) = cgrid%fmean_snowfac (ipy) & + + csite%fmean_snowfac (ipa) & + * patch_wgt cgrid%fmean_rshort_gnd (ipy) = cgrid%fmean_rshort_gnd (ipy) & + csite%fmean_rshort_gnd (ipa) & * patch_wgt @@ -1425,11 +1428,13 @@ subroutine zero_ed_fmean_vars(cgrid) cgrid%fmean_sfcw_mass ( ipy) = 0.0 cgrid%fmean_sfcw_temp ( ipy) = 0.0 cgrid%fmean_sfcw_fliq ( ipy) = 0.0 + cgrid%fmean_snowfac ( ipy) = 0.0 cgrid%fmean_soil_energy (:,ipy) = 0.0 cgrid%fmean_soil_mstpot (:,ipy) = 0.0 cgrid%fmean_soil_water (:,ipy) = 0.0 cgrid%fmean_soil_temp (:,ipy) = 0.0 cgrid%fmean_soil_fliq (:,ipy) = 0.0 + cgrid%fmean_snowfac ( ipy) = 0.0 cgrid%fmean_rshort_gnd ( ipy) = 0.0 cgrid%fmean_par_gnd ( ipy) = 0.0 cgrid%fmean_rlong_gnd ( ipy) = 0.0 @@ -1574,6 +1579,7 @@ subroutine zero_ed_fmean_vars(cgrid) csite%fmean_sfcw_mass ( ipa) = 0.0 csite%fmean_sfcw_temp ( ipa) = 0.0 csite%fmean_sfcw_fliq ( ipa) = 0.0 + csite%fmean_snowfac ( ipa) = 0.0 csite%fmean_soil_energy (:,ipa) = 0.0 csite%fmean_soil_mstpot (:,ipa) = 0.0 csite%fmean_soil_water (:,ipa) = 0.0 @@ -2059,6 +2065,9 @@ subroutine integrate_ed_dmean_vars(cgrid) cgrid%dmean_sfcw_mass (ipy) = cgrid%dmean_sfcw_mass (ipy) & + cgrid%fmean_sfcw_mass (ipy) & * frqsum_o_daysec + cgrid%dmean_snowfac (ipy) = cgrid%dmean_snowfac (ipy) & + + cgrid%fmean_snowfac (ipy) & + * frqsum_o_daysec cgrid%dmean_rshort_gnd (ipy) = cgrid%dmean_rshort_gnd (ipy) & + cgrid%fmean_rshort_gnd (ipy) & * frqsum_o_daysec @@ -2363,6 +2372,9 @@ subroutine integrate_ed_dmean_vars(cgrid) csite%dmean_sfcw_mass (ipa) = csite%dmean_sfcw_mass (ipa) & + csite%fmean_sfcw_mass (ipa) & * frqsum_o_daysec + csite%dmean_snowfac (ipa) = csite%dmean_snowfac (ipa) & + + csite%fmean_snowfac (ipa) & + * frqsum_o_daysec csite%dmean_rshort_gnd (ipa) = csite%dmean_rshort_gnd (ipa) & + csite%fmean_rshort_gnd (ipa) & * frqsum_o_daysec @@ -2812,10 +2824,11 @@ end subroutine normalize_ed_today_vars !=======================================================================================! !=======================================================================================! - ! SUBROUTINE: NORMALIZE_ED_TODAYNPP_VARS - !> \brief This subroutine will scale the daily NPP allocation terms + ! SUBROUTINE: COPY_TODAY_TO_DMEAN_VARS + !> \brief This subroutine scales the daily NPP allocation terms and transfer today + !! variables to dmean. !---------------------------------------------------------------------------------------! - subroutine normalize_ed_todayNPP_vars(cgrid) + subroutine copy_today_to_dmean_vars(cgrid) use ed_state_vars , only : edtype & ! structure , polygontype & ! structure , sitetype & ! structure @@ -2833,47 +2846,70 @@ subroutine normalize_ed_todayNPP_vars(cgrid) integer :: isi integer :: ipa integer :: ico + !------------------------------------------------------------------------------------! + + !------ Skip sub-routine and do nothing in case we are not writing daily averages. --! + if (.not. writing_long) return + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Loop through polygons. ! + !------------------------------------------------------------------------------------! polyloop: do ipy=1,cgrid%npolygons cpoly => cgrid%polygon(ipy) + + !---------------------------------------------------------------------------------! + ! Loop through sites. ! + !---------------------------------------------------------------------------------! siteloop: do isi=1,cpoly%nsites csite => cpoly%site(isi) - patchloop: do ipa=1,csite%npatches + !------------------------------------------------------------------------------! + ! Loop through patches. ! + !------------------------------------------------------------------------------! + patchloop: do ipa=1,csite%npatches cpatch => csite%patch(ipa) - - !----- Included a loop so it won't crash with empty cohorts... -------------! + + !---------------------------------------------------------------------------! + ! Loop through cohorts. This must be a loop, so it won't crash with ! + ! empty patches. ! + !---------------------------------------------------------------------------! cohortloop: do ico=1,cpatch%ncohorts !------------------------------------------------------------------------! ! We now update the daily means of NPP allocation terms ! ! and we convert them to kgC/plant/yr ! !------------------------------------------------------------------------! - if (writing_long) then - cpatch%dmean_nppleaf (ico) = cpatch%today_nppleaf (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppfroot (ico) = cpatch%today_nppfroot (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppsapwood(ico) = cpatch%today_nppsapwood(ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppbark (ico) = cpatch%today_nppbark (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppcroot (ico) = cpatch%today_nppcroot (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppseeds (ico) = cpatch%today_nppseeds (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppwood (ico) = cpatch%today_nppwood (ico) & - * yr_day / cpatch%nplant (ico) - cpatch%dmean_nppdaily (ico) = cpatch%today_nppdaily (ico) & - * yr_day / cpatch%nplant (ico) - end if + cpatch%dmean_nppleaf (ico) = cpatch%today_nppleaf (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppfroot (ico) = cpatch%today_nppfroot (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppsapwood(ico) = cpatch%today_nppsapwood(ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppbark (ico) = cpatch%today_nppbark (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppcroot (ico) = cpatch%today_nppcroot (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppseeds (ico) = cpatch%today_nppseeds (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppwood (ico) = cpatch%today_nppwood (ico) & + * yr_day / cpatch%nplant (ico) + cpatch%dmean_nppdaily (ico) = cpatch%today_nppdaily (ico) & + * yr_day / cpatch%nplant (ico) + !------------------------------------------------------------------------! end do cohortloop + !---------------------------------------------------------------------------! end do patchloop + !------------------------------------------------------------------------------! end do siteloop + !---------------------------------------------------------------------------------! end do polyloop + !------------------------------------------------------------------------------------! return - end subroutine normalize_ed_todayNPP_vars + end subroutine copy_today_to_dmean_vars !=======================================================================================! !=======================================================================================! @@ -3702,6 +3738,7 @@ subroutine zero_ed_dmean_vars(cgrid) cgrid%dmean_sfcw_mass (ipy) = 0.0 cgrid%dmean_sfcw_temp (ipy) = 0.0 cgrid%dmean_sfcw_fliq (ipy) = 0.0 + cgrid%dmean_snowfac (ipy) = 0.0 cgrid%dmean_soil_energy (:,ipy) = 0.0 cgrid%dmean_soil_mstpot (:,ipy) = 0.0 cgrid%dmean_soil_water (:,ipy) = 0.0 @@ -3824,6 +3861,7 @@ subroutine zero_ed_dmean_vars(cgrid) csite%dmean_sfcw_mass (ipa) = 0.0 csite%dmean_sfcw_temp (ipa) = 0.0 csite%dmean_sfcw_fliq (ipa) = 0.0 + csite%dmean_snowfac (ipa) = 0.0 csite%dmean_soil_energy (:,ipa) = 0.0 csite%dmean_soil_mstpot (:,ipa) = 0.0 csite%dmean_soil_water (:,ipa) = 0.0 @@ -4426,6 +4464,9 @@ subroutine integrate_ed_mmean_vars(cgrid) cgrid%mmean_sfcw_mass (ipy) = cgrid%mmean_sfcw_mass (ipy) & + cgrid%dmean_sfcw_mass (ipy) & * ndaysi + cgrid%mmean_snowfac (ipy) = cgrid%mmean_snowfac (ipy) & + + cgrid%dmean_snowfac (ipy) & + * ndaysi cgrid%mmean_soil_energy (:,ipy) = cgrid%mmean_soil_energy (:,ipy) & + cgrid%dmean_soil_energy (:,ipy) & * ndaysi @@ -4951,6 +4992,9 @@ subroutine integrate_ed_mmean_vars(cgrid) csite%mmean_sfcw_mass (ipa) = csite%mmean_sfcw_mass (ipa) & + csite%dmean_sfcw_mass (ipa) & * ndaysi + csite%mmean_snowfac (ipa) = csite%mmean_snowfac (ipa) & + + csite%dmean_snowfac (ipa) & + * ndaysi csite%mmean_soil_energy (:,ipa) = csite%mmean_soil_energy (:,ipa) & + csite%dmean_soil_energy (:,ipa) & * ndaysi @@ -6086,6 +6130,7 @@ subroutine zero_ed_mmean_vars(cgrid) cgrid%mmean_sfcw_mass (ipy) = 0.0 cgrid%mmean_sfcw_temp (ipy) = 0.0 cgrid%mmean_sfcw_fliq (ipy) = 0.0 + cgrid%mmean_snowfac (ipy) = 0.0 cgrid%mmean_soil_energy (:,ipy) = 0.0 cgrid%mmean_soil_mstpot (:,ipy) = 0.0 cgrid%mmean_soil_water (:,ipy) = 0.0 @@ -6274,6 +6319,7 @@ subroutine zero_ed_mmean_vars(cgrid) csite%mmean_sfcw_mass (ipa) = 0.0 csite%mmean_sfcw_temp (ipa) = 0.0 csite%mmean_sfcw_fliq (ipa) = 0.0 + csite%mmean_snowfac (ipa) = 0.0 csite%mmean_soil_energy (:,ipa) = 0.0 csite%mmean_soil_mstpot (:,ipa) = 0.0 csite%mmean_soil_water (:,ipa) = 0.0 @@ -6838,6 +6884,9 @@ subroutine integrate_ed_qmean_vars(cgrid) cgrid%qmean_sfcw_mass (t,ipy) = cgrid%qmean_sfcw_mass (t,ipy) & + cgrid%fmean_sfcw_mass (ipy) & * ndaysi + cgrid%qmean_snowfac (t,ipy) = cgrid%qmean_snowfac (t,ipy) & + + cgrid%fmean_snowfac (ipy) & + * ndaysi cgrid%qmean_soil_energy (:,t,ipy) = cgrid%qmean_soil_energy (:,t,ipy) & + cgrid%fmean_soil_energy (:,ipy) & * ndaysi @@ -7230,6 +7279,9 @@ subroutine integrate_ed_qmean_vars(cgrid) csite%qmean_sfcw_mass (t,ipa) = csite%qmean_sfcw_mass (t,ipa) & + csite%fmean_sfcw_mass (ipa) & * ndaysi + csite%qmean_snowfac (t,ipa) = csite%qmean_snowfac (t,ipa) & + + csite%fmean_snowfac (ipa) & + * ndaysi csite%qmean_soil_energy (:,t,ipa) = csite%qmean_soil_energy (:,t,ipa) & + csite%fmean_soil_energy (:,ipa) & * ndaysi @@ -7796,6 +7848,25 @@ subroutine normalize_ed_qmean_vars(cgrid) + !---------------------------------------------------------------------------! + ! Find the derived properties for the canopy air space. ! + !---------------------------------------------------------------------------! + do t=1,ndcycle + can_exner = press2exner ( csite%qmean_can_prss (t,ipa)) + csite%qmean_can_temp(t,ipa) = extheta2temp( can_exner & + , csite%qmean_can_theta(t,ipa)) + csite%qmean_can_rhos(t,ipa) = idealdenssh ( csite%qmean_can_prss (t,ipa) & + , csite%qmean_can_temp (t,ipa) & + , csite%qmean_can_shv (t,ipa)) + csite%qmean_can_dmol(t,ipa) = idealdmolsh ( csite%qmean_can_prss (t,ipa) & + , csite%qmean_can_temp (t,ipa) & + , csite%qmean_can_shv (t,ipa)) + end do + !---------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------! ! Soil matric potential, temperature, and liquid water. ! !---------------------------------------------------------------------------! @@ -8163,6 +8234,7 @@ subroutine zero_ed_qmean_vars(cgrid) cgrid%qmean_sfcw_mass (:,ipy) = 0.0 cgrid%qmean_sfcw_temp (:,ipy) = 0.0 cgrid%qmean_sfcw_fliq (:,ipy) = 0.0 + cgrid%qmean_snowfac (:,ipy) = 0.0 cgrid%qmean_soil_energy (:,:,ipy) = 0.0 cgrid%qmean_soil_mstpot (:,:,ipy) = 0.0 cgrid%qmean_soil_water (:,:,ipy) = 0.0 @@ -8309,6 +8381,7 @@ subroutine zero_ed_qmean_vars(cgrid) csite%qmean_sfcw_mass (:,ipa) = 0.0 csite%qmean_sfcw_temp (:,ipa) = 0.0 csite%qmean_sfcw_fliq (:,ipa) = 0.0 + csite%qmean_snowfac (:,ipa) = 0.0 csite%qmean_soil_energy (:,:,ipa) = 0.0 csite%qmean_soil_mstpot (:,:,ipa) = 0.0 csite%qmean_soil_water (:,:,ipa) = 0.0 diff --git a/ED/src/io/ed_init_history.f90 b/ED/src/io/ed_init_history.f90 index f4d91cc49..d9c77c35a 100644 --- a/ED/src/io/ed_init_history.f90 +++ b/ED/src/io/ed_init_history.f90 @@ -1004,6 +1004,8 @@ subroutine fill_history_grid_p11dmean(cgrid,ipy,py_index) ,'DMEAN_SFCW_TEMP_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%dmean_sfcw_fliq (ipy:ipy) & ,'DMEAN_SFCW_FLIQ_PY ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(cgrid%dmean_snowfac (ipy:ipy) & + ,'DMEAN_SNOWFAC_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%dmean_rshort_gnd (ipy:ipy) & ,'DMEAN_RSHORT_GND_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%dmean_par_gnd (ipy:ipy) & @@ -1432,6 +1434,8 @@ subroutine fill_history_grid_p11mmean(cgrid,ipy,py_index) ,'MMEAN_SFCW_TEMP_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%mmean_sfcw_fliq (ipy:ipy) & ,'MMEAN_SFCW_FLIQ_PY ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(cgrid%mmean_snowfac (ipy:ipy) & + ,'MMEAN_SNOWFAC_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%mmean_rshort_gnd (ipy:ipy) & ,'MMEAN_RSHORT_GND_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%mmean_par_gnd (ipy:ipy) & @@ -2055,6 +2059,8 @@ subroutine fill_history_grid_m11(cgrid,ipy,py_index) ,'QMEAN_SFCW_TEMP_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%qmean_sfcw_fliq (:,ipy) & ,'QMEAN_SFCW_FLIQ_PY ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(cgrid%qmean_snowfac (:,ipy) & + ,'QMEAN_SNOWFAC_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%qmean_rshort_gnd (:,ipy) & ,'QMEAN_RSHORT_GND_PY ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%qmean_par_gnd (:,ipy) & @@ -3257,8 +3263,8 @@ subroutine fill_history_polygon(cpoly,pysi_index,nsites_global,nsites_now,is_bur memoffs (2) = 0_8 call hdf_getslab_r(cpoly%lambda_fire & ,'LAMBDA_FIRE ',dsetrank,iparallel,.true.,foundvar) - call hdf_getslab_r(cpoly%avg_monthly_pcpg & - ,'AVG_MONTHLY_PCPG ',dsetrank,iparallel,.true.,foundvar) + call hdf_getslab_r(cpoly%avg_monthly_accp & + ,'AVG_MONTHLY_ACCP ',dsetrank,iparallel,.true.,foundvar) call hdf_getslab_r(cpoly%crop_yield & ,'CROP_YIELD_SI ',dsetrank,iparallel,.true.,foundvar) !------------------------------------------------------------------------------------! @@ -4069,6 +4075,8 @@ subroutine fill_history_site(csite,sipa_index,npatches_global,is_burnt) ,'DMEAN_SFCW_TEMP_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%dmean_sfcw_fliq & ,'DMEAN_SFCW_FLIQ_PA ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(csite%dmean_snowfac & + ,'DMEAN_SNOWFAC_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%dmean_rshort_gnd & ,'DMEAN_RSHORT_GND_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%dmean_par_gnd & @@ -4238,6 +4246,8 @@ subroutine fill_history_site(csite,sipa_index,npatches_global,is_burnt) ,'MMEAN_SFCW_TEMP_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%mmean_sfcw_fliq & ,'MMEAN_SFCW_FLIQ_PA ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(csite%mmean_snowfac & + ,'MMEAN_SNOWFAC_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%mmean_rshort_gnd & ,'MMEAN_RSHORT_GND_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%mmean_par_gnd & @@ -4443,6 +4453,8 @@ subroutine fill_history_site(csite,sipa_index,npatches_global,is_burnt) ,'QMEAN_SFCW_TEMP_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%qmean_sfcw_fliq & ,'QMEAN_SFCW_FLIQ_PA ',dsetrank,iparallel,.false.,foundvar) + call hdf_getslab_r(csite%qmean_snowfac & + ,'QMEAN_SNOWFAC_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%qmean_rshort_gnd & ,'QMEAN_RSHORT_GND_PA ',dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(csite%qmean_par_gnd & diff --git a/ED/src/io/ed_read_ed10_20_history.f90 b/ED/src/io/ed_read_ed10_20_history.f90 index 083c6e93e..23d16cf17 100644 --- a/ED/src/io/ed_read_ed10_20_history.f90 +++ b/ED/src/io/ed_read_ed10_20_history.f90 @@ -882,11 +882,16 @@ subroutine read_ed10_ed20_history_file !------------------------------------------------------------------! !------------------------------------------------------------------! - ! Initialise SLA with the look-up table value, this may be ! - ! updated during phenology initialisation, but an initial assign- ! - ! ment is needed to obtain area indices. ! + ! Initialise SLA, Vm0, Rd0, and with the look-up table value. ! + ! These variables may be updated during phenology initialisation, ! + ! or trait plasticity, but they must have an initial assignment so ! + ! we can even calculate the initial area indices and inicial trait ! + ! values needed for the trait update. ! !------------------------------------------------------------------! - cpatch%sla(ic2) = SLA(ipft(ic)) + cpatch%sla (ic2) = SLA (ipft(ic)) + cpatch%vm_bar(ic2) = Vm0 (ipft(ic)) + cpatch%rd_bar(ic2) = Rd0 (ipft(ic)) + cpatch%llspan(ic2) = leaf_lifespan(ipft(ic)) !------------------------------------------------------------------! diff --git a/ED/src/memory/ed_state_vars.F90 b/ED/src/memory/ed_state_vars.F90 index 0c111ae29..d779778fa 100644 --- a/ED/src/memory/ed_state_vars.F90 +++ b/ED/src/memory/ed_state_vars.F90 @@ -1869,6 +1869,7 @@ module ed_state_vars real,pointer,dimension(:) :: fmean_sfcw_mass ! tiny_sfcwater_mass) then csite%fmean_sfcw_energy (recp) = csite%fmean_sfcw_energy (recp) & @@ -7775,6 +7777,8 @@ subroutine fuse_2_patches(csite,donp,recp,mzg,mzs,cmet,lsl,ntext_soil,green_leaf * csite%dmean_sfcw_mass (recp) * rawgt csite%dmean_sfcw_mass (recp) = csite%dmean_sfcw_mass (donp) * dawgt & + csite%dmean_sfcw_mass (recp) * rawgt + csite%dmean_snowfac (recp) = csite%dmean_snowfac (donp) * dawgt & + + csite%dmean_snowfac (recp) * rawgt !----- Check whether there is enough surface water. ---------------------------! if (csite%dmean_sfcw_mass(recp) > tiny_sfcwater_mass) then csite%dmean_sfcw_energy (recp) = csite%dmean_sfcw_energy (recp) & @@ -8169,6 +8173,8 @@ subroutine fuse_2_patches(csite,donp,recp,mzg,mzs,cmet,lsl,ntext_soil,green_leaf * csite%mmean_sfcw_mass (recp) * rawgt csite%mmean_sfcw_mass (recp) = csite%mmean_sfcw_mass (donp) * dawgt & + csite%mmean_sfcw_mass (recp) * rawgt + csite%mmean_snowfac (recp) = csite%mmean_snowfac (donp) * dawgt & + + csite%mmean_snowfac (recp) * rawgt !----- Check whether there is enough surface water. ---------------------------! if (csite%mmean_sfcw_mass(recp) > tiny_sfcwater_mass) then csite%mmean_sfcw_energy (recp) = csite%mmean_sfcw_energy (recp) & @@ -8520,6 +8526,8 @@ subroutine fuse_2_patches(csite,donp,recp,mzg,mzs,cmet,lsl,ntext_soil,green_leaf * csite%qmean_sfcw_mass (t,recp) * rawgt csite%qmean_sfcw_mass (t,recp) = csite%qmean_sfcw_mass (t,donp) * dawgt & + csite%qmean_sfcw_mass (t,recp) * rawgt + csite%qmean_snowfac (t,recp) = csite%qmean_snowfac (t,donp) * dawgt & + + csite%qmean_snowfac (t,recp) * rawgt !----- Check whether there is enough surface water. ------------------------! if (csite%qmean_sfcw_mass(t,recp) > tiny_sfcwater_mass) then csite%qmean_sfcw_energy (t,recp) = csite%qmean_sfcw_energy(t,recp) &