program mksurfdat !----------------------------------------------------------------------- !BOP ! ! !PROGRAM: mksurfdat ! ! !DESCRIPTION: ! Creates land model surface dataset from original "raw" data files. ! Surface dataset contains model grid, pfts, inland water, glacier, ! soil texture, soil color, soil order,LAI and SAI, urban fraction, and urban ! parameters. ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8, r4 => shr_kind_r4 use fileutils , only : opnfil, getavu use mklaiMod , only : mklai use mkpftMod , only : pft_idx, pft_frc, mkpft, mkpftInit, mkpft_parse_oride, & natpft_lb, natpft_ub, cft_lb, cft_ub, num_cft use mksoilMod , only : soil_sand, soil_clay, mksoiltex, mksoilInit, & soil_color, mksoilcol, mkorganic, & soil_fmax, mkfmax,& soil_order, mksoilord use mkvocefMod , only : mkvocef use mklanwatMod , only : mklakwat, mkwetlnd, mklakparams use mkglcmecMod , only : nglcec, mkglcmec, mkglcmecInit, mkglacier use mkharvestMod , only : mkharvest, mkharvest_init, mkharvest_fieldname, & mkharvest_numtypes, mkharvest_parse_oride use mkurbanparCommonMod, only : mkelev use mkurbanparMod , only : mkurbanInit, mkurban, mkurbanpar, numurbl use mkutilsMod , only : normalize_classes_by_gcell use mkfileMod , only : mkfile use mkvarpar , only : nlevsoi, nlevslp, elev_thresh use mkvarctl use nanMod , only : nan, bigint use mkncdio , only : check_ret use mkdomainMod , only : domain_type, domain_read_map, domain_read, & domain_write use mkgdpMod , only : mkgdp use mkpeatMod , only : mkpeat use mkagfirepkmonthMod , only : mkagfirepkmon use mktopostatsMod , only : mktopostats use mkVICparamsMod , only : mkVICparams use mkCH4inversionMod , only : mkCH4inversion use mksoilphosphorusMod, only : mksoilphosphorus use mkSedMod , only : mkgrvl, mkslp10, mkEROparams ! ! !ARGUMENTS: implicit none include 'netcdf.inc' ! ! !REVISION HISTORY: ! Authors: Gordon Bonan, Sam Levis and Mariana Vertenstein ! Revised: Nan Rosenbloom to add fmax processing. ! 3/18/08: David Lawrence added organic matter processing ! 1/22/09: Keith Oleson added urban parameter processing ! 2/11/13: Sam Levis added abm, peat, and gdp processing for new fire model ! 3/27/15: Xiaoying Shi added soil order for new Phosphorus models ! ! ! !LOCAL VARIABLES: !EOP integer :: nsoicol ! number of model color classes integer :: nsoiord ! number of model order classes integer :: k,m,n ! indices integer :: ni,nj,ns_o ! indices integer :: ier ! error status integer :: ndiag,nfdyn ! unit numbers integer :: ncid ! netCDF id integer :: omode ! netCDF output mode integer :: varid ! netCDF variable id integer :: ndims ! netCDF number of dimensions integer :: beg(4),len(4),dimids(4) ! netCDF dimension sizes integer :: ret ! netCDF return status integer :: ntim ! time sample for dynamic land use integer :: year ! year for dynamic land use logical :: all_veg ! if gridcell will be 100% vegetated land-cover real(r8) :: suma ! sum for error check character(len=256) :: fgrddat ! grid data file character(len=256) :: fsurdat ! output surface data file name character(len=256) :: fsurlog ! output surface log file name character(len=256) :: fdyndat ! dynamic landuse data file name character(len=256) :: fname ! generic filename character(len=256) :: string ! string read in integer :: t1 ! timer real(r8),parameter :: p5 = 0.5_r8 ! constant real(r8),parameter :: p25 = 0.25_r8 ! constant real(r8), allocatable :: landfrac_pft(:) ! PFT data: % land per gridcell real(r8), allocatable :: pctlnd_pft(:) ! PFT data: % of gridcell for PFTs real(r8), allocatable :: pctlnd_pft_dyn(:) ! PFT data: % of gridcell for dyn landuse PFTs integer , allocatable :: pftdata_mask(:) ! mask indicating real or fake land type real(r8), pointer :: pctpft_full(:,:) ! PFT data: % cover of each pft and cft on the vegetated landunits ! ('full' denotes inclusion of CFTs as well as natural PFTs in this array) real(r8), allocatable :: pctnatveg(:) ! percent of grid cell that is natural veg landunit real(r8), allocatable :: pctcrop(:) ! percent of grid cell that is crop landunit real(r8), allocatable :: pctnatpft(:,:) ! % of each pft on the natural veg landunit (adds to 100%) real(r8), allocatable :: pctcft(:,:) ! % of each cft on the crop landunit (adds to 100%) real(r8), pointer :: harvest(:,:) ! harvest data: normalized harvesting real(r8), allocatable :: pctgla(:) ! percent of grid cell that is glacier real(r8), allocatable :: pctglc_gic(:) ! percent of grid cell that is gic (% of glc landunit) real(r8), allocatable :: pctglc_icesheet(:) ! percent of grid cell that is ice sheet (% of glc landunit) real(r8), allocatable :: pctglcmec(:,:) ! glacier_mec pct coverage in each class (% of landunit) real(r8), allocatable :: topoglcmec(:,:) ! glacier_mec sfc elevation in each gridcell and class real(r8), allocatable :: pctglcmec_gic(:,:) ! GIC pct coverage in each class (% of landunit) real(r8), allocatable :: pctglcmec_icesheet(:,:) ! icesheet pct coverage in each class (% of landunit) real(r8), allocatable :: elevclass(:) ! glacier_mec elevation classes real(r8), allocatable :: pctlak(:) ! percent of grid cell that is lake real(r8), allocatable :: pctwet(:) ! percent of grid cell that is wetland real(r8), allocatable :: pcturb(:) ! percent of grid cell that is urbanized (total across all urban classes) real(r8), allocatable :: urbn_classes(:,:) ! percent cover of each urban class, as % of total urban area real(r8), allocatable :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell real(r8), allocatable :: elev(:) ! glc elevation (m) real(r8), allocatable :: topo(:) ! land elevation (m) real(r8), allocatable :: fmax(:) ! fractional saturated area integer , allocatable :: soicol(:) ! soil color integer , allocatable :: soiord(:) ! soil order real(r8), allocatable :: pctsand(:,:) ! soil texture: percent sand real(r8), allocatable :: pctclay(:,:) ! soil texture: percent clay real(r8), allocatable :: ef1_btr(:) ! Isoprene emission factor for broadleaf real(r8), allocatable :: ef1_fet(:) ! Isoprene emission factor for fine/everg real(r8), allocatable :: ef1_fdt(:) ! Isoprene emission factor for fine/dec real(r8), allocatable :: ef1_shr(:) ! Isoprene emission factor for shrubs real(r8), allocatable :: ef1_grs(:) ! Isoprene emission factor for grasses real(r8), allocatable :: ef1_crp(:) ! Isoprene emission factor for crops real(r8), allocatable :: organic(:,:) ! organic matter density (kg/m3) real(r8), allocatable :: gdp(:) ! GDP (x1000 1995 US$/capita) real(r8), allocatable :: fpeat(:) ! peatland fraction of gridcell integer , allocatable :: agfirepkmon(:) ! agricultural fire peak month integer , allocatable :: urban_region(:) ! urban region ID real(r8), allocatable :: topo_stddev(:) ! standard deviation of elevation (m) real(r8), allocatable :: slope(:) ! topographic slope (degrees) real(r8), allocatable :: vic_binfl(:) ! VIC b parameter (unitless) real(r8), allocatable :: vic_ws(:) ! VIC Ws parameter (unitless) real(r8), allocatable :: vic_dsmax(:) ! VIC Dsmax parameter (mm/day) real(r8), allocatable :: vic_ds(:) ! VIC Ds parameter (unitless) real(r8), allocatable :: lakedepth(:) ! lake depth (m) real(r8), allocatable :: f0(:) ! max fractional inundated area (unitless) real(r8), allocatable :: p3(:) ! coefficient for qflx_surf_lag for finundated (s/mm) real(r8), allocatable :: zwt0(:) ! decay factor for finundated (m) real(r8), allocatable :: apatiteP(:) ! apptite phosphorus real(r8), allocatable :: labileP(:) ! labile phosphorus real(r8), allocatable :: occludedP(:) ! occluded phosphorus real(r8), allocatable :: secondaryP(:) ! secondaryP phosphorus real(r8), allocatable :: grvl(:,:) ! soil gravel content (percent) real(r8), allocatable :: slp10(:,:) ! slope percentile (km/km) real(r8), allocatable :: ero_c1(:) ! ELM-Erosion c1 parameter (unitless) real(r8), allocatable :: ero_c2(:) ! ELM-Erosion c2 parameter (unitless) real(r8), allocatable :: ero_c3(:) ! ELM-Erosion c3 parameter (unitless) type(domain_type) :: ldomain character(len=32) :: subname = 'mksrfdat' ! program name namelist /elmexp/ & mksrf_fgrid, & mksrf_gridtype, & mksrf_fvegtyp, & mksrf_fsoitex, & mksrf_forganic, & mksrf_fsoicol, & mksrf_fsoiord, & mksrf_fvocef, & mksrf_flakwat, & mksrf_fwetlnd, & mksrf_fglacier, & mksrf_furbtopo, & mksrf_flndtopo, & mksrf_fmax, & mksrf_furban, & mksrf_flai, & mksrf_fdynuse, & mksrf_fgdp, & mksrf_fpeat, & mksrf_fabm, & mksrf_ftopostats, & mksrf_fvic, & mksrf_fch4, & mksrf_fphosphorus, & mksrf_fgrvl, & mksrf_fslp10, & mksrf_fero, & nglcec, & numpft, & soil_color, & soil_order, & soil_sand, & soil_fmax, & soil_clay, & pft_idx, & pft_frc, & all_urban, & no_inlandwet, & map_fpft, & map_flakwat, & map_fwetlnd, & map_fglacier, & map_fsoitex, & map_fsoicol, & map_fsoiord, & map_furban, & map_furbtopo, & map_flndtopo, & map_fmax, & map_forganic, & map_fvocef, & map_flai, & map_fharvest, & map_fgdp, & map_fpeat, & map_fabm, & map_ftopostats, & map_fvic, & map_fch4, & map_fphosphorus, & map_fgrvl, & map_fslp10, & map_fero, & outnc_large_files, & outnc_double, & outnc_dims, & fsurdat, & fdyndat, & fsurlog !----------------------------------------------------------------------- ! ====================================================================== ! Read input namelist ! ====================================== ! Must specify settings for the output grid: ! ====================================== ! mksrf_fgrid -- Grid dataset ! ====================================== ! Must specify settings for input high resolution datafiles ! ====================================== ! mksrf_fglacier - Glacier dataset ! mksrf_flai ----- Leaf Area Index dataset ! mksrf_flakwat -- Lake water dataset ! mksrf_fwetlnd -- Wetland water dataset ! mksrf_forganic - Organic soil carbon dataset ! mksrf_fmax ----- Max fractional saturated area dataset ! mksrf_fsoicol -- Soil color dataset ! mksrf_fsoiord -- Soil order dataset ! mksrf_fsoitex -- Soil texture dataset ! mksrf_furbtopo-- Topography dataset (for limiting urban areas) ! mksrf_furban --- Urban dataset ! mksrf_fvegtyp -- PFT vegetation type dataset ! mksrf_fvocef -- Volatile Organic Compund Emission Factor dataset ! mksrf_fgdp ----- GDP dataset ! mksrf_fpeat ---- Peatland dataset ! mksrf_fabm ----- Agricultural fire peak month dataset ! mksrf_ftopostats Topography statistics dataset ! mksrf_fvic ----- VIC parameters dataset ! mksrf_fch4 ----- inversion-derived CH4 parameters dataset ! mksrf_fphosphorus Soil phosphorus dataset ! mksrf_fgrvl ---- Soil gravel content dataset ! mksrf_fslp10 --- Slope percentile dataset ! mksrf_fero ----- ELM-Erosion parameters dataset ! ====================================== ! Must specify mapping file for the different datafiles above ! ====================================== ! map_fpft -------- Mapping for mksrf_fvegtyp ! map_flakwat ----- Mapping for mksrf_flakwat ! map_fwetlnd ----- Mapping for mksrf_fwetlnd ! map_fglacier ---- Mapping for mksrf_fglacier ! map_fsoitex ----- Mapping for mksrf_fsoitex ! map_fsoicol ----- Mapping for mksrf_fsoicol ! map_fsoiord ----- Mapping for mksrf_fsoiord ! map_furban ------ Mapping for mksrf_furban ! map_furbtopo ---- Mapping for mksrf_furbtopo ! map_flndtopo ---- Mapping for mksrf_flndtopo ! map_fmax -------- Mapping for mksrf_fmax ! map_forganic ---- Mapping for mksrf_forganic ! map_fvocef ------ Mapping for mksrf_fvocef ! map_flai -------- Mapping for mksrf_flai ! map_fharvest ---- Mapping for mksrf_flai harvesting ! map_fgdp -------- Mapping for mksrf_fgdp ! map_fpeat ------- Mapping for mksrf_fpeat ! map_fabm -------- Mapping for mksrf_fabm ! map_ftopostats -- Mapping for mksrf_ftopostats ! map_fvic -------- Mapping for mksrf_fvic ! map_fch4 -------- Mapping for mksrf_fch4 ! map_fphosphorus - Mapping for mksrf_fphosphorus ! map_fgrvl ------- Mapping for mksrf_fgrvl ! map_fslp10 ------ Mapping for mksrf_fslp10 ! map_fero -------- Mapping for mksrf_fero ! ====================================== ! Optionally specify setting for: ! ====================================== ! mksrf_fdynuse ----- ASCII text file that lists each year of pft files to use ! mksrf_gridtype ---- Type of grid (default is 'global') ! outnc_double ------ If output should be in double precision ! outnc_large_files - If output should be in NetCDF large file format ! nglcec ------------ If you want to change the number of Glacier elevation classes ! ====================================== ! Optional settings to change values for entire area ! ====================================== ! all_urban --------- If entire area is urban ! no_inlandwet ------ If wetland should be set to 0% over land ! soil_color -------- If you want to change the soil_color to this value everywhere ! soil_order -------- If you want to change the soil_order to this value everywhere ! soil_clay --------- If you want to change the soil_clay % to this value everywhere ! soil_fmax --------- If you want to change the soil_fmax to this value everywhere ! soil_sand --------- If you want to change the soil_sand % to this value everywhere ! pft_idx ----------- If you want to change to 100% veg covered with given PFT indices ! pft_frc ----------- Fractions that correspond to the pft_idx above ! ================== ! numpft (if different than default of 16) ! ====================================================================== write(6,*) 'Attempting to initialize control settings .....' mksrf_gridtype = 'global' outnc_large_files = .false. outnc_double = .true. all_urban = .false. no_inlandwet = .true. read(5, elmexp, iostat=ier) if (ier /= 0) then write(6,*)'error: namelist input resulted in error code ',ier call abort() endif write (6,*) 'Attempting to create surface boundary data .....' write (6,'(72a1)') ("-",n=1,60) ! ---------------------------------------------------------------------- ! Error check namelist input ! ---------------------------------------------------------------------- if (mksrf_fgrid /= ' ')then fgrddat = mksrf_fgrid write(6,*)'mksrf_fgrid = ',mksrf_fgrid else write (6,*)'must specify mksrf_fgrid' call abort() endif if (trim(mksrf_gridtype) == 'global' .or. & trim(mksrf_gridtype) == 'regional') then write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype) else write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype) write (6,*)'illegal mksrf_gridtype, must be global or regional ' call abort() endif if ( outnc_large_files )then write(6,*)'Output file in NetCDF 64-bit large_files format' end if if ( outnc_double )then write(6,*)'Output ALL data in file as 64-bit' end if if ( all_urban )then write(6,*) 'Output ALL data in file as 100% urban' end if if ( no_inlandwet )then write(6,*) 'Set wetland to 0% over land' end if ! ! Call module initialization routines ! call mksoilInit( ) call mkpftInit( all_urban, all_veg ) allocate ( elevclass(nglcec+1) ) call mkglcmecInit (elevclass) call mkurbanInit (mksrf_furban) if ( all_veg )then write(6,*) 'Output ALL data in file as 100% vegetated' end if ! ---------------------------------------------------------------------- ! Determine land model grid, fractional land and land mask ! ---------------------------------------------------------------------- write(6,*)'calling domain_read' if ( .not. domain_read_map(ldomain, fgrddat) )then call domain_read(ldomain, fgrddat) end if write(6,*)'finished domain_read' ! Invalidate mask and frac for ldomain !ldomain%mask = bigint !ldomain%frac = nan ! Determine if will have 1d output if (ldomain%ni /= -9999 .and. ldomain%nj /= -9999) then write(6,*)'fsurdat is 2d lat/lon grid' write(6,*)'nlon= ',ldomain%ni,' nlat= ',ldomain%nj if (outnc_dims == 1) then write(6,*)' writing output file in 1d gridcell format' end if else write(6,*)'fsurdat is 1d gridcell grid' outnc_dims = 1 end if outnc_1d = .false. if ((ldomain%ni == -9999 .and. ldomain%nj == -9999) .or. outnc_dims==1) then outnc_1d = .true. write(6,*)'output file will be 1d' end if ! ---------------------------------------------------------------------- ! Allocate and initialize dynamic memory ! ---------------------------------------------------------------------- ns_o = ldomain%ns allocate ( landfrac_pft(ns_o) , & pctlnd_pft(ns_o) , & pftdata_mask(ns_o) , & pctpft_full(ns_o,0:numpft) , & pctnatveg(ns_o) , & pctcrop(ns_o) , & pctnatpft(ns_o,natpft_lb:natpft_ub), & pctcft(ns_o,cft_lb:cft_ub) , & pctgla(ns_o) , & pctlak(ns_o) , & pctwet(ns_o) , & pcturb(ns_o) , & urban_region(ns_o) , & urbn_classes(ns_o,numurbl) , & urbn_classes_g(ns_o,numurbl) , & pctsand(ns_o,nlevsoi) , & pctclay(ns_o,nlevsoi) , & soicol(ns_o) , & soiord(ns_o) , & gdp(ns_o) , & fpeat(ns_o) , & agfirepkmon(ns_o) , & topo_stddev(ns_o) , & slope(ns_o) , & vic_binfl(ns_o) , & vic_ws(ns_o) , & vic_dsmax(ns_o) , & vic_ds(ns_o) , & lakedepth(ns_o) , & f0(ns_o) , & p3(ns_o) , & zwt0(ns_o) , & apatiteP(ns_o) , & labileP(ns_o) , & occludedP(ns_o) , & secondaryP(ns_o) , & grvl(ns_o,nlevsoi) , & slp10(ns_o,nlevslp) , & ero_c1(ns_o) , & ero_c2(ns_o) , & ero_c3(ns_o) ) landfrac_pft(:) = spval pctlnd_pft(:) = spval pftdata_mask(:) = -999 pctpft_full(:,:) = spval pctnatveg(:) = spval pctcrop(:) = spval pctnatpft(:,:) = spval pctcft(:,:) = spval pctgla(:) = spval pctlak(:) = spval pctwet(:) = spval pcturb(:) = spval urban_region(:) = -999 urbn_classes(:,:) = spval urbn_classes_g(:,:) = spval pctsand(:,:) = spval pctclay(:,:) = spval soicol(:) = -999 soiord(:) = -999 gdp(:) = spval fpeat(:) = spval agfirepkmon(:) = -999 topo_stddev(:) = spval slope(:) = spval vic_binfl(:) = spval vic_ws(:) = spval vic_dsmax(:) = spval vic_ds(:) = spval lakedepth(:) = spval f0(:) = spval p3(:) = spval zwt0(:) = spval apatiteP(:) = spval labileP(:) = spval occludedP(:) = spval secondaryP(:) = spval grvl(:,:) = spval slp10(:,:) = spval ero_c1(:) = spval ero_c2(:) = spval ero_c3(:) = spval ! ---------------------------------------------------------------------- ! Open diagnostic output log file ! ---------------------------------------------------------------------- if (fsurlog == ' ') then write(6,*)' must specify fsurlog in namelist' stop else ndiag = getavu(); call opnfil (fsurlog, ndiag, 'f') end if if (mksrf_fgrid /= ' ')then write (ndiag,*)'using fractional land data from file= ', & trim(mksrf_fgrid),' to create the surface dataset' endif if (trim(mksrf_gridtype) == 'global' .or. & trim(mksrf_gridtype) == 'regional') then write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype) endif write(ndiag,*) 'PFTs from: ',trim(mksrf_fvegtyp) write(ndiag,*) 'fmax from: ',trim(mksrf_fmax) write(ndiag,*) 'glaciers from: ',trim(mksrf_fglacier) write(ndiag,*) ' with: ', nglcec, ' glacier elevation classes' write(ndiag,*) 'urban topography from: ',trim(mksrf_furbtopo) write(ndiag,*) 'land topography from: ',trim(mksrf_flndtopo) write(ndiag,*) 'urban from: ',trim(mksrf_furban) write(ndiag,*) 'inland lake from: ',trim(mksrf_flakwat) write(ndiag,*) 'inland wetland from: ',trim(mksrf_fwetlnd) write(ndiag,*) 'soil texture from: ',trim(mksrf_fsoitex) write(ndiag,*) 'soil organic from: ',trim(mksrf_forganic) write(ndiag,*) 'soil color from: ',trim(mksrf_fsoicol) write(ndiag,*) 'soil order from: ',trim(mksrf_fsoiord) write(ndiag,*) 'VOC emission factors from: ',trim(mksrf_fvocef) write(ndiag,*) 'gdp from: ',trim(mksrf_fgdp) write(ndiag,*) 'peat from: ',trim(mksrf_fpeat) write(ndiag,*) 'abm from: ',trim(mksrf_fabm) write(ndiag,*) 'topography statistics from: ',trim(mksrf_ftopostats) write(ndiag,*) 'VIC parameters from: ',trim(mksrf_fvic) write(ndiag,*) 'CH4 parameters from: ',trim(mksrf_fch4) write(ndiag,*) 'Soil phosphorus from: ',trim(mksrf_fphosphorus) write(ndiag,*)' mapping for pft ',trim(map_fpft) write(ndiag,*)' mapping for lake water ',trim(map_flakwat) write(ndiag,*)' mapping for wetland ',trim(map_fwetlnd) write(ndiag,*)' mapping for glacier ',trim(map_fglacier) write(ndiag,*)' mapping for soil texture ',trim(map_fsoitex) write(ndiag,*)' mapping for soil color ',trim(map_fsoicol) write(ndiag,*)' mapping for soil order ',trim(map_fsoiord) write(ndiag,*)' mapping for soil organic ',trim(map_forganic) write(ndiag,*)' mapping for urban ',trim(map_furban) write(ndiag,*)' mapping for fmax ',trim(map_fmax) write(ndiag,*)' mapping for VOC pct emis ',trim(map_fvocef) write(ndiag,*)' mapping for harvest ',trim(map_fharvest) write(ndiag,*)' mapping for lai/sai ',trim(map_flai) write(ndiag,*)' mapping for urb topography ',trim(map_furbtopo) write(ndiag,*)' mapping for land topography ',trim(map_flndtopo) write(ndiag,*)' mapping for GDP ',trim(map_fgdp) write(ndiag,*)' mapping for peatlands ',trim(map_fpeat) write(ndiag,*)' mapping for ag fire pk month ',trim(map_fabm) write(ndiag,*)' mapping for topography stats ',trim(map_ftopostats) write(ndiag,*)' mapping for VIC parameters ',trim(map_fvic) write(ndiag,*)' mapping for CH4 parameters ',trim(map_fch4) write(ndiag,*)' mapping for soil phosphorus ',trim(map_fphosphorus) write(ndiag,*)' mapping for soil gravel ',trim(map_fgrvl) write(ndiag,*)' mapping for slope percentile ',trim(map_fslp10) write(ndiag,*)' mapping for erosion params ',trim(map_fero) if (mksrf_fdynuse /= ' ') then write(6,*)'mksrf_fdynuse = ',trim(mksrf_fdynuse) end if ! ---------------------------------------------------------------------- ! Make surface dataset fields ! ---------------------------------------------------------------------- ! Make PFTs [pctpft_full] from dataset [fvegtyp] call mkpft(ldomain, mapfname=map_fpft, fpft=mksrf_fvegtyp, & ndiag=ndiag, pctlnd_o=pctlnd_pft, pctpft_o=pctpft_full ) ! Make inland water [pctlak, pctwet] [flakwat] [fwetlnd] call mklakwat (ldomain, mapfname=map_flakwat, datfname=mksrf_flakwat, & ndiag=ndiag, zero_out=all_urban.or.all_veg, lake_o=pctlak) call mkwetlnd (ldomain, mapfname=map_fwetlnd, datfname=mksrf_fwetlnd, & ndiag=ndiag, zero_out=all_urban.or.all_veg.or.no_inlandwet, swmp_o=pctwet) ! Make glacier fraction [pctgla] from [fglacier] dataset call mkglacier (ldomain, mapfname=map_fglacier, datfname=mksrf_fglacier, & ndiag=ndiag, zero_out=all_urban.or.all_veg, glac_o=pctgla) ! Make soil texture [pctsand, pctclay] [fsoitex] call mksoiltex (ldomain, mapfname=map_fsoitex, datfname=mksrf_fsoitex, & ndiag=ndiag, sand_o=pctsand, clay_o=pctclay) ! Make soil color classes [soicol] [fsoicol] call mksoilcol (ldomain, mapfname=map_fsoicol, datfname=mksrf_fsoicol, & ndiag=ndiag, soil_color_o=soicol, nsoicol=nsoicol) ! Make soil order classes [soiord] [fsoiord] call mksoilord (ldomain, mapfname=map_fsoiord, datfname=mksrf_fsoiord, & ndiag=ndiag, pctglac_o=pctgla, soil_order_o=soiord, nsoiord=nsoiord) ! Make fmax [fmax] from [fmax] dataset allocate(fmax(ns_o)) fmax(:) = spval call mkfmax (ldomain, mapfname=map_fmax, datfname=mksrf_fmax, & ndiag=ndiag, fmax_o=fmax) ! Make GDP data [gdp] from [gdp] call mkgdp (ldomain, mapfname=map_fgdp, datfname=mksrf_fgdp, & ndiag=ndiag, gdp_o=gdp) ! Make peat data [fpeat] from [peatf] call mkpeat (ldomain, mapfname=map_fpeat, datfname=mksrf_fpeat, & ndiag=ndiag, peat_o=fpeat) ! Make agricultural fire peak month data [abm] from [abm] call mkagfirepkmon (ldomain, mapfname=map_fabm, datfname=mksrf_fabm, & ndiag=ndiag, agfirepkmon_o=agfirepkmon) ! Make urban fraction [pcturb] from [furban] dataset call mkurban (ldomain, mapfname=map_furban, datfname=mksrf_furban, & ndiag=ndiag, zero_out=all_veg, urbn_o=pcturb, urbn_classes_o=urbn_classes, & region_o=urban_region) ! WJS (9-25-12): Note about topo datasets: Until now, there have been two topography ! datasets: flndtopo & fglctopo. flndtopo is used to create the TOPO variable, which I ! believe is used to downscale grid cell-level climate to glc_mec columns (10-26-12: ! Now I'm not surue about this: I think TOPO might actually come from a different file ! in CLM, and TOPO on the surface dataset may be unused). Until now, fglctopo was used ! for dividing pct_glacier data into multiple elevation classes in ! mkglcmecMod. However, it is no longer needed for this purpose, since elevation data ! is now folded into fglacier. fglctopo has also been used to screen urban points (I'm ! not sure why fglctopo rather than flndtopo was chosen for that purpose). ! ! For now, I am keeping fglctopo around simply for the urban screening purpose. To ! make its purpose clear, I am renaming it to furbtopo. I had planned to switch to a ! new topo file that is consistent with the topo data that are implicitly included in ! fglacier (i.e., a file that gives the topo that's used for glc purposes, even though ! fglctopo itself isn't used for glc purposes any more). However, this caused problems ! in coming up with a new elev_thresh. Thus, for now I am continuing to use the old ! fglctopo file, which no longer has any meaning with respect to glc (and again, I am ! renaming it to furbtopo to make it clear that it is not connected with glc). ! ! In the longer term, a better solution for this urban screening would probably be to ! modify the raw urban data. In that case, I believe we could remove furbtopo. ! ! Why was TOPO created from flndtopo rather than fglctopo? It seems like, for the ! purpose of downscaling, this TOPO variable should ideally represent CAM's ! topographic height. For that purpose, flndtopo is more appropriate, because it seems ! to have come from CAM's topo dataset. However, I believe that many (all??) CAM ! resolutions use some sort of smoothed topography. So the ideal thing to do would be ! for CLM to get its grid cell-level topography from CAM at initialization. If that ! were done, then I think flndtopo and the TOPO variable on the surface dataset could ! go away. (Update 10-26-12: it actually looks to me like CLM's TOPO comes from a ! different source entirely (flndtopo in CLM), so it may be that TOPO on the surface ! dataset isn't currently used for anything!) ! Make elevation [elev] from [ftopo, ffrac] dataset ! Used only to screen pcturb ! Screen pcturb by elevation threshold from elev dataset if ( .not. all_urban .and. .not. all_veg )then allocate(elev(ns_o)) elev(:) = spval call mkelev (ldomain, mapfname=map_furbtopo, datfname=mksrf_furbtopo, & varname='TOPO_ICE', ndiag=ndiag, elev_o=elev) where (elev .gt. elev_thresh) pcturb = 0._r8 end where deallocate(elev) end if ! Determine topography allocate(topo(ns_o)) call mkelev (ldomain, mapfname=map_flndtopo, datfname=mksrf_flndtopo, & varname='TOPO', ndiag=ndiag, elev_o=topo) ! Compute topography statistics [topo_stddev, slope] from [ftopostats] call mktopostats (ldomain, mapfname=map_ftopostats, datfname=mksrf_ftopostats, & ndiag=ndiag, topo_stddev_o=topo_stddev, slope_o=slope) ! Make VIC parameters [binfl, ws, dsmax, ds] from [fvic] call mkVICparams (ldomain, mapfname=map_fvic, datfname=mksrf_fvic, ndiag=ndiag, & binfl_o=vic_binfl, ws_o=vic_ws, dsmax_o=vic_dsmax, ds_o=vic_ds) ! Make lake depth [lakedepth] from [flakwat] call mklakparams (ldomain, mapfname=map_flakwat, datfname=mksrf_flakwat, ndiag=ndiag, & lakedepth_o=lakedepth) ! Make inversion-derived CH4 parameters [f0, p3, zwt0] from [fch4] call mkCH4inversion (ldomain, mapfname=map_fch4, datfname=mksrf_fch4, ndiag=ndiag, & f0_o=f0, p3_o=p3, zwt0_o=zwt0) ! Make organic matter density [organic] [forganic] allocate (organic(ns_o,nlevsoi)) organic(:,:) = spval call mkorganic (ldomain, mapfname=map_forganic, datfname=mksrf_forganic, & ndiag=ndiag, organic_o=organic) ! Make VOC emission factors for isoprene & ! [ef1_btr,ef1_fet,ef1_fdt,ef1_shr,ef1_grs,ef1_crp] allocate ( ef1_btr(ns_o) , & ef1_fet(ns_o) , & ef1_fdt(ns_o) , & ef1_shr(ns_o) , & ef1_grs(ns_o) , & ef1_crp(ns_o) ) ef1_btr(:) = 0._r8 ef1_fet(:) = 0._r8 ef1_fdt(:) = 0._r8 ef1_shr(:) = 0._r8 ef1_grs(:) = 0._r8 ef1_crp(:) = 0._r8 call mkvocef (ldomain, mapfname=map_fvocef, datfname=mksrf_fvocef, ndiag=ndiag, & ef_btr_o=ef1_btr, ef_fet_o=ef1_fet, ef_fdt_o=ef1_fdt, & ef_shr_o=ef1_shr, ef_grs_o=ef1_grs, ef_crp_o=ef1_crp) ! Do landuse changes such as for the poles, etc. call change_landuse( ldomain, dynpft=.false. ) call mksoilphosphorus (ldomain, mapfname=map_fphosphorus, datfname=mksrf_fphosphorus, & ndiag=ndiag, apatiteP_o=apatiteP, labileP_o=labileP, occludedP_o=occludedP, & secondaryP_o=secondaryP) call mkgrvl(ldomain, mapfname=map_fgrvl, datfname=mksrf_fgrvl, ndiag=ndiag, grvl_o=grvl) call mkslp10(ldomain, mapfname=map_fslp10, datfname=mksrf_fslp10, ndiag=ndiag, slp10_o=slp10) call mkEROparams(ldomain, mapfname=map_fero, datfname=mksrf_fero, ndiag=ndiag, & ero_c1_o=ero_c1, ero_c2_o=ero_c2, ero_c3_o=ero_c3) do n = 1,ns_o ! Assume wetland and/or lake when dataset landmask implies ocean ! (assume medium soil color (15), soil order(15) and loamy texture). ! Also set pftdata_mask here ! Note that pctpft_full is NOT adjusted here, so that we still have information ! about the landunit breakdown into PFTs (pctnatveg and pctcrop will later become 0 ! due to wetland and lake adding to 100%). if (pctlnd_pft(n) < 1.e-6_r8) then pftdata_mask(n) = 0 soicol(n) = 15 soiord(n) = 15 pctwet(n) = 100._r8 - pctlak(n) pcturb(n) = 0._r8 pctgla(n) = 0._r8 pctsand(n,:) = 43._r8 pctclay(n,:) = 18._r8 organic(n,:) = 0._r8 grvl(n,:) = 0._r8 slp10(n,:) = 0._r8 ero_c1(n) = 0._r8 ero_c2(n) = 0._r8 ero_c3(n) = 0._r8 else pftdata_mask(n) = 1 end if ! Truncate all percentage fields on output grid. This is needed to ! insure that wt is zero (not a very small number such as ! 1e-16) where it really should be zero do k = 1,nlevsoi pctsand(n,k) = float(nint(pctsand(n,k))) pctclay(n,k) = float(nint(pctclay(n,k))) grvl(n,k) = float(nint(grvl(n,k))) end do pctlak(n) = float(nint(pctlak(n))) pctwet(n) = float(nint(pctwet(n))) pctgla(n) = float(nint(pctgla(n))) ! Make sure sum of land cover types does not exceed 100. If it does, ! subtract excess from most dominant land cover. suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) if (suma > 250._r4) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb and pctgla is greater than 250%' write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n) call abort() else if (suma > 100._r4) then pctlak(n) = pctlak(n) * 100._r8/suma pctwet(n) = pctwet(n) * 100._r8/suma pcturb(n) = pcturb(n) * 100._r8/suma pctgla(n) = pctgla(n) * 100._r8/suma end if end do call normalizencheck_landuse(ldomain) ! Write out sum of PFT's do k = natpft_lb,natpft_ub suma = 0._r8 do n = 1,ns_o suma = suma + pctnatpft(n,k)*pctnatveg(n)/100._r8 enddo write(6,*) 'sum over domain of pft ',k,suma enddo write(6,*) do k = cft_lb,cft_ub suma = 0._r8 do n = 1,ns_o suma = suma + pctcft(n,k)*pctcrop(n)/100._r8 enddo write(6,*) 'sum over domain of cft ',k,suma enddo write(6,*) ! Make final values of percent urban by class ! This call needs to occur after all corrections are made to pcturb call normalize_classes_by_gcell(urbn_classes, pcturb, urbn_classes_g) ! Make glacier multiple elevation classes [pctglcmec,topoglcmec] from [fglacier] dataset ! This call needs to occur after pctgla has been adjusted for the final time if ( nglcec > 0 )then allocate (pctglcmec(ns_o,nglcec), & topoglcmec(ns_o,nglcec), & pctglcmec_gic(ns_o,nglcec), & pctglcmec_icesheet(ns_o,nglcec)) allocate (pctglc_gic(ns_o)) allocate (pctglc_icesheet(ns_o)) pctglcmec(:,:) = spval topoglcmec(:,:) = spval pctglcmec_gic(:,:) = spval pctglcmec_icesheet(:,:) = spval pctglc_gic(:) = spval pctglc_icesheet(:) = spval call mkglcmec (ldomain, mapfname=map_fglacier, & datfname_fglacier=mksrf_fglacier, ndiag=ndiag, & pctglcmec_o=pctglcmec, topoglcmec_o=topoglcmec, & pctglcmec_gic_o=pctglcmec_gic, pctglcmec_icesheet_o=pctglcmec_icesheet, & pctglc_gic_o=pctglc_gic, pctglc_icesheet_o=pctglc_icesheet) end if ! Determine fractional land from pft dataset do n = 1,ns_o landfrac_pft(n) = pctlnd_pft(n)/100._r8 end do ! ---------------------------------------------------------------------- ! Create surface dataset ! ---------------------------------------------------------------------- ! Create netCDF surface dataset. if (fsurdat == ' ') then write(6,*)' must specify fsurdat in namelist' stop end if call mkfile(ldomain, trim(fsurdat), dynlanduse = .false.) call domain_write(ldomain, fsurdat) call check_ret(nf_open(trim(fsurdat), nf_write, ncid), subname) call check_ret(nf_set_fill (ncid, nf_nofill, omode), subname) ! Write fields OTHER THAN lai, sai, heights, and urban parameters to netcdf surface dataset call check_ret(nf_inq_varid(ncid, 'natpft', varid), subname) call check_ret(nf_put_var_int(ncid, varid, (/(n,n=natpft_lb,natpft_ub)/)), subname) if (num_cft > 0) then call check_ret(nf_inq_varid(ncid, 'cft', varid), subname) call check_ret(nf_put_var_int(ncid, varid, (/(n,n=cft_lb,cft_ub)/)), subname) end if call check_ret(nf_inq_varid(ncid, 'PFTDATA_MASK', varid), subname) call check_ret(nf_put_var_int(ncid, varid, pftdata_mask), subname) call check_ret(nf_inq_varid(ncid, 'LANDFRAC_PFT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, landfrac_pft), subname) call check_ret(nf_inq_varid(ncid, 'mxsoil_color', varid), subname) call check_ret(nf_put_var_int(ncid, varid, nsoicol), subname) call check_ret(nf_inq_varid(ncid, 'SOIL_COLOR', varid), subname) call check_ret(nf_put_var_int(ncid, varid, soicol), subname) call check_ret(nf_inq_varid(ncid, 'mxsoil_order', varid), subname) call check_ret(nf_put_var_int(ncid, varid, nsoiord), subname) call check_ret(nf_inq_varid(ncid, 'SOIL_ORDER', varid), subname) call check_ret(nf_put_var_int(ncid, varid, soiord), subname) call check_ret(nf_inq_varid(ncid, 'PCT_SAND', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctsand), subname) call check_ret(nf_inq_varid(ncid, 'PCT_CLAY', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctclay), subname) call check_ret(nf_inq_varid(ncid, 'PCT_WETLAND', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctwet), subname) call check_ret(nf_inq_varid(ncid, 'PCT_LAKE', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctlak), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLACIER', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctgla), subname) if ( nglcec > 0 )then call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctglcmec), subname) call check_ret(nf_inq_varid(ncid, 'GLC_MEC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, elevclass), subname) call check_ret(nf_inq_varid(ncid, 'TOPO_GLC_MEC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, topoglcmec), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_GIC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctglcmec_gic), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_ICESHEET', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctglcmec_icesheet), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLC_GIC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctglc_gic), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLC_ICESHEET', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctglc_icesheet), subname) end if call check_ret(nf_inq_varid(ncid, 'TOPO', varid), subname) call check_ret(nf_put_var_double(ncid, varid, topo), subname) call check_ret(nf_inq_varid(ncid, 'PCT_URBAN', varid), subname) call check_ret(nf_put_var_double(ncid, varid, urbn_classes_g), subname) call check_ret(nf_inq_varid(ncid, 'PCT_NATVEG', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctnatveg), subname) call check_ret(nf_inq_varid(ncid, 'PCT_CROP', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctcrop), subname) call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctnatpft), subname) if (num_cft > 0) then call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctcft), subname) end if call check_ret(nf_inq_varid(ncid, 'FMAX', varid), subname) call check_ret(nf_put_var_double(ncid, varid, fmax), subname) call check_ret(nf_inq_varid(ncid, 'gdp', varid), subname) call check_ret(nf_put_var_double(ncid, varid, gdp), subname) call check_ret(nf_inq_varid(ncid, 'peatf', varid), subname) call check_ret(nf_put_var_double(ncid, varid, fpeat), subname) call check_ret(nf_inq_varid(ncid, 'abm', varid), subname) call check_ret(nf_put_var_int(ncid, varid, agfirepkmon), subname) call check_ret(nf_inq_varid(ncid, 'SLOPE', varid), subname) call check_ret(nf_put_var_double(ncid, varid, slope), subname) call check_ret(nf_inq_varid(ncid, 'STD_ELEV', varid), subname) call check_ret(nf_put_var_double(ncid, varid, topo_stddev), subname) call check_ret(nf_inq_varid(ncid, 'binfl', varid), subname) call check_ret(nf_put_var_double(ncid, varid, vic_binfl), subname) call check_ret(nf_inq_varid(ncid, 'Ws', varid), subname) call check_ret(nf_put_var_double(ncid, varid, vic_ws), subname) call check_ret(nf_inq_varid(ncid, 'Dsmax', varid), subname) call check_ret(nf_put_var_double(ncid, varid, vic_dsmax), subname) call check_ret(nf_inq_varid(ncid, 'Ds', varid), subname) call check_ret(nf_put_var_double(ncid, varid, vic_ds), subname) call check_ret(nf_inq_varid(ncid, 'LAKEDEPTH', varid), subname) call check_ret(nf_put_var_double(ncid, varid, lakedepth), subname) call check_ret(nf_inq_varid(ncid, 'F0', varid), subname) call check_ret(nf_put_var_double(ncid, varid, f0), subname) call check_ret(nf_inq_varid(ncid, 'P3', varid), subname) call check_ret(nf_put_var_double(ncid, varid, p3), subname) call check_ret(nf_inq_varid(ncid, 'ZWT0', varid), subname) call check_ret(nf_put_var_double(ncid, varid, zwt0), subname) call check_ret(nf_inq_varid(ncid, 'EF1_BTR', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_btr), subname) call check_ret(nf_inq_varid(ncid, 'EF1_FET', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_fet), subname) call check_ret(nf_inq_varid(ncid, 'EF1_FDT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_fdt), subname) call check_ret(nf_inq_varid(ncid, 'EF1_SHR', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_shr), subname) call check_ret(nf_inq_varid(ncid, 'EF1_GRS', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_grs), subname) call check_ret(nf_inq_varid(ncid, 'EF1_CRP', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ef1_crp), subname) call check_ret(nf_inq_varid(ncid, 'ORGANIC', varid), subname) call check_ret(nf_put_var_double(ncid, varid, organic), subname) call check_ret(nf_inq_varid(ncid, 'URBAN_REGION_ID', varid), subname) call check_ret(nf_put_var_int(ncid, varid, urban_region), subname) call check_ret(nf_inq_varid(ncid, 'APATITE_P', varid), subname) call check_ret(nf_put_var_double(ncid, varid, apatiteP), subname) call check_ret(nf_inq_varid(ncid, 'LABILE_P', varid), subname) call check_ret(nf_put_var_double(ncid, varid, labileP), subname) call check_ret(nf_inq_varid(ncid, 'OCCLUDED_P', varid), subname) call check_ret(nf_put_var_double(ncid, varid, occludedP), subname) call check_ret(nf_inq_varid(ncid, 'SECONDARY_P', varid), subname) call check_ret(nf_put_var_double(ncid, varid, secondaryP), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GRVL', varid), subname) call check_ret(nf_put_var_double(ncid, varid, grvl), subname) call check_ret(nf_inq_varid(ncid, 'SLP_P10', varid), subname) call check_ret(nf_put_var_double(ncid, varid, slp10), subname) call check_ret(nf_inq_varid(ncid, 'parEro_c1', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ero_c1), subname) call check_ret(nf_inq_varid(ncid, 'parEro_c2', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ero_c2), subname) call check_ret(nf_inq_varid(ncid, 'parEro_c3', varid), subname) call check_ret(nf_put_var_double(ncid, varid, ero_c3), subname) ! Deallocate arrays NOT needed for dynamic-pft section of code deallocate ( organic ) deallocate ( ef1_btr, ef1_fet, ef1_fdt, ef1_shr, ef1_grs, ef1_crp ) if ( nglcec > 0 ) deallocate ( pctglcmec, topoglcmec) if ( nglcec > 0 ) deallocate ( pctglc_gic, pctglc_icesheet) deallocate ( elevclass ) deallocate ( fmax ) deallocate ( pctsand, pctclay ) deallocate ( soicol ) deallocate ( soiord ) deallocate ( gdp, fpeat, agfirepkmon ) deallocate ( topo_stddev, slope ) deallocate ( vic_binfl, vic_ws, vic_dsmax, vic_ds ) deallocate ( lakedepth ) deallocate ( f0, p3, zwt0 ) deallocate ( apatiteP, labileP, occludedP, secondaryP ) deallocate ( grvl, slp10, ero_c1, ero_c2, ero_c3 ) ! Synchronize the disk copy of a netCDF dataset with in-memory buffers call check_ret(nf_sync(ncid), subname) ! ---------------------------------------------------------------------- ! Make Urban Parameters from raw input data and write to surface dataset ! Write to netcdf file is done inside mkurbanpar routine ! ---------------------------------------------------------------------- write(6,*)'calling mkurbanpar' call mkurbanpar(datfname=mksrf_furban, ncido=ncid, region_o=urban_region, & urbn_classes_gcell_o=urbn_classes_g) ! ---------------------------------------------------------------------- ! Make LAI and SAI from 1/2 degree data and write to surface dataset ! Write to netcdf file is done inside mklai routine ! ---------------------------------------------------------------------- write(6,*)'calling mklai' call mklai(ldomain, mapfname=map_flai, datfname=mksrf_flai, & ndiag=ndiag, ncido=ncid ) ! Close surface dataset call check_ret(nf_close(ncid), subname) write (6,'(72a1)') ("-",n=1,60) write (6,*)' land model surface data set successfully created for ', & 'grid of size ',ns_o ! ---------------------------------------------------------------------- ! Create dynamic land use dataset if appropriate ! ---------------------------------------------------------------------- if (mksrf_fdynuse /= ' ') then write(6,*)'creating dynamic land use dataset' allocate(pctlnd_pft_dyn(ns_o)) call mkharvest_init( ns_o, spval, harvest, mksrf_fvegtyp ) if (fdyndat == ' ') then write(6,*)' must specify fdyndat in namelist if mksrf_fdynuse is not blank' stop end if ! Define dimensions and global attributes call mkfile(ldomain, fdyndat, dynlanduse=.true.) ! Write fields other pft to dynamic land use dataset call domain_write(ldomain, fdyndat) call check_ret(nf_open(trim(fdyndat), nf_write, ncid), subname) call check_ret(nf_set_fill (ncid, nf_nofill, omode), subname) call check_ret(nf_inq_varid(ncid, 'natpft', varid), subname) call check_ret(nf_put_var_int(ncid, varid, (/(n,n=natpft_lb,natpft_ub)/)), subname) if (num_cft > 0) then call check_ret(nf_inq_varid(ncid, 'cft', varid), subname) call check_ret(nf_put_var_int(ncid, varid, (/(n,n=cft_lb,cft_ub)/)), subname) end if call check_ret(nf_inq_varid(ncid, 'PFTDATA_MASK', varid), subname) call check_ret(nf_put_var_int(ncid, varid, pftdata_mask), subname) call check_ret(nf_inq_varid(ncid, 'LANDFRAC_PFT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, landfrac_pft), subname) call check_ret(nf_inq_varid(ncid, 'PCT_WETLAND', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctwet), subname) call check_ret(nf_inq_varid(ncid, 'PCT_LAKE', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctlak), subname) call check_ret(nf_inq_varid(ncid, 'PCT_GLACIER', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctgla), subname) call check_ret(nf_inq_varid(ncid, 'PCT_URBAN', varid), subname) call check_ret(nf_put_var_double(ncid, varid, urbn_classes_g), subname) call check_ret(nf_inq_varid(ncid, 'PCT_NATVEG', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctnatveg), subname) call check_ret(nf_inq_varid(ncid, 'PCT_CROP', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctcrop), subname) if (num_cft > 0) then call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname) call check_ret(nf_put_var_double(ncid, varid, pctcft), subname) end if ! Synchronize the disk copy of a netCDF dataset with in-memory buffers call check_ret(nf_sync(ncid), subname) ! Read in each dynamic pft landuse dataset nfdyn = getavu(); call opnfil (mksrf_fdynuse, nfdyn, 'f') ntim = 0 do ! Read input pft data read(nfdyn, '(A195,1x,I4)', iostat=ier) string, year if (ier /= 0) exit ! ! If pft fraction override is set, than intrepret string as PFT and harvesting override values ! if ( any(pft_frc > 0.0_r8 ) )then fname = ' ' call mkpft_parse_oride(string) call mkharvest_parse_oride(string) write(6,*)'PFT and harvesting values are ',trim(string),' year is ',year ! ! Otherwise intrepret string as a filename with PFT and harvesting values in it ! else fname = string write(6,*)'input pft dynamic dataset is ',trim(fname),' year is ',year end if ntim = ntim + 1 ! Re-allocate pctpft_full, because it is deallocated in normalizencheck_landuse ! for safety allocate(pctpft_full(ns_o,0:numpft)) ! Create pctpft data at model resolution call mkpft(ldomain, mapfname=map_fpft, fpft=fname, & ndiag=ndiag, pctlnd_o=pctlnd_pft_dyn, pctpft_o=pctpft_full ) ! Create harvesting data at model resolution call mkharvest( ldomain, mapfname=map_fharvest, datfname=fname, & ndiag=ndiag, harv_o=harvest ) ! Consistency check on input land fraction do n = 1,ns_o if (pctlnd_pft_dyn(n) /= pctlnd_pft(n)) then write(6,*) subname,' error: pctlnd_pft for dynamics data = ',& pctlnd_pft_dyn(n), ' not equal to pctlnd_pft for surface data = ',& pctlnd_pft(n),' at n= ',n if ( trim(fname) == ' ' )then write(6,*) ' PFT string = ', string else write(6,*) ' PFT file = ', fname end if call abort() end if end do call change_landuse(ldomain, dynpft=.true.) call normalizencheck_landuse(ldomain) ! Output pctpft data for current year call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname) call check_ret(nf_inq_varndims(ncid, varid, ndims), subname) call check_ret(nf_inq_vardimid(ncid, varid, dimids), subname) beg(1:ndims-1) = 1 do n = 1,ndims-1 call check_ret(nf_inq_dimlen(ncid, dimids(n), len(n)), subname) end do len(ndims) = 1 beg(ndims) = ntim call check_ret(nf_put_vara_double(ncid, varid, beg, len, pctnatpft), subname) do k = 1, mkharvest_numtypes() call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(k)), varid), subname) call check_ret(nf_inq_varndims(ncid, varid, ndims), subname) call check_ret(nf_inq_vardimid(ncid, varid, dimids), subname) beg(1:ndims-1) = 1 do n = 1,ndims-1 call check_ret(nf_inq_dimlen(ncid, dimids(n), len(n)), subname) end do len(ndims) = 1 beg(ndims) = ntim call check_ret(nf_put_vara_double(ncid, varid, beg, len, harvest(:,k)), subname) end do call check_ret(nf_inq_varid(ncid, 'YEAR', varid), subname) call check_ret(nf_put_vara_int(ncid, varid, ntim, 1, year), subname) call check_ret(nf_inq_varid(ncid, 'time', varid), subname) call check_ret(nf_put_vara_int(ncid, varid, ntim, 1, year), subname) call check_ret(nf_inq_varid(ncid, 'input_pftdata_filename', varid), subname) call check_ret(nf_put_vara_text(ncid, varid, (/ 1, ntim /), (/ len_trim(string), 1 /), trim(string) ), subname) ! Synchronize the disk copy of a netCDF dataset with in-memory buffers call check_ret(nf_sync(ncid), subname) end do ! end of read loop call check_ret(nf_close(ncid), subname) end if ! end of if-create dynamic landust dataset ! ---------------------------------------------------------------------- ! Close diagnostic dataset ! ---------------------------------------------------------------------- close (ndiag) write (6,*) write (6,*) 'Surface data output file = ',trim(fsurdat) write (6,*) ' This file contains the land model surface data' write (6,*) 'Diagnostic log file = ',trim(fsurlog) write (6,*) ' See this file for a summary of the dataset' write (6,*) write (6,*) 'Successfully created surface dataset' !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: change_landuse ! ! !INTERFACE: subroutine change_landuse( ldomain, dynpft ) ! ! !DESCRIPTION: ! ! Do landuse changes such as for the poles, etc. ! ! !USES: implicit none ! ! !ARGUMENTS: type(domain_type) :: ldomain logical, intent(in) :: dynpft ! if part of the dynpft section of code ! ! !REVISION HISTORY: ! 9/10/09: Erik Kluzek spin off subroutine from original embedded code ! !EOP ! ! !LOCAL VARIABLES: logical :: first_time = .true. ! flag if this is the first pass through or not integer ,parameter :: bdtroptree = 6 ! Index for broadleaf decidious tropical tree integer ,parameter :: bdtemptree = 7 ! Index for broadleaf decidious temperate tree integer ,parameter :: bdtempshrub = 10 ! Index for broadleaf decidious temperate shrub real(r8),parameter :: troplat = 23.5_r8 ! Latitude to define as tropical integer :: n,ns_o ! indices character(len=32) :: subname = 'change_landuse' ! subroutine name !----------------------------------------------------------------------- ns_o = ldomain%ns do n = 1,ns_o ! Set pfts 7 and 10 to 6 in the tropics to avoid lais > 1000 ! Using P. Thornton's method found in surfrdMod.F90 in clm3.5 if (abs(ldomain%latc(n))0._r8) then pctpft_full(n,bdtroptree) = pctpft_full(n,bdtroptree) + pctpft_full(n,bdtemptree) pctpft_full(n,bdtemptree) = 0._r8 if ( first_time ) write (6,*) subname, ' Warning: all wgt of pft ', & bdtemptree, ' now added to pft ', bdtroptree end if if (abs(ldomain%latc(n))0._r8) then pctpft_full(n,bdtroptree) = pctpft_full(n,bdtroptree) + pctpft_full(n,bdtempshrub) pctpft_full(n,bdtempshrub) = 0._r8 if ( first_time ) write (6,*) subname, ' Warning: all wgt of pft ', & bdtempshrub, ' now added to pft ', bdtroptree end if first_time = .false. ! If have pole points on grid - set south pole to glacier ! north pole is assumed as non-land ! Note that pctpft_full is NOT adjusted here, so that we still have information ! about the landunit breakdown into PFTs (pctnatveg and pctcrop will later become 0 ! due to glacier adding to 100%). if (abs((ldomain%latc(n) - 90._r8)) < 1.e-6_r8) then pctlak(n) = 0._r8 pctwet(n) = 0._r8 pcturb(n) = 0._r8 pctgla(n) = 100._r8 if ( .not. dynpft )then organic(n,:) = 0._r8 ef1_btr(n) = 0._r8 ef1_fet(n) = 0._r8 ef1_fdt(n) = 0._r8 ef1_shr(n) = 0._r8 ef1_grs(n) = 0._r8 ef1_crp(n) = 0._r8 end if end if end do end subroutine change_landuse !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: normalizencheck_landuse ! ! !INTERFACE: subroutine normalizencheck_landuse(ldomain) ! ! !DESCRIPTION: ! ! Normalize land use and make sure things add up to 100% as well as ! checking that things are as they should be. ! ! Note: this deallocates pctpft_full ! ! !USES: use mkutilsMod, only : remove_small_cover use mkpftMod , only : mkpft_normalize implicit none ! !ARGUMENTS: type(domain_type) :: ldomain ! ! !REVISION HISTORY: ! 9/10/09: Erik Kluzek spin off subroutine from original embedded code ! !EOP ! ! !LOCAL VARIABLES: integer :: m,k,n,ns_o ! indices integer :: nsmall ! number of small PFT values for a single check integer :: nsmall_tot ! total number of small PFT values in all grid cells real(r8) :: suma ! sum for error check real(r8) :: suma2 ! another sum for error check real(r8) :: pcturb_rescaled ! % urban normalized as if it were on the veg landunit real(r8) :: bare_urb_diff ! difference between bare soil and urban % real(r8) :: pcturb_excess ! excess urban % not accounted for by bare soil real(r8) :: sumpft ! sum of non-baresoil pfts real(r8) :: sum8, sum8a ! sum for error check real(r4) :: sum4a ! sum for error check real(r8), parameter :: toosmallPFT = 1.e-10_r8 ! tolerance for PFT's to ignore character(len=32) :: subname = 'normalizencheck_landuse' ! subroutine name !----------------------------------------------------------------------- ns_o = ldomain%ns do n = 1,ns_o ! The following corrects pctpft_full to account for urban area. Urban needs to be ! handled specially because we replace bare soil preferentially with urban, rather ! than rescaling all PFTs equally. To accomplish this, we first rescale urban to ! the value it would have if it were in the pctpft_full array (pcturb_rescaled), which ! gives the fraction of each pft on the (vegetated + crop) landunits. Before ! (conceptually) adding urban to this array, pctpft_full added up to 100% of the ! (vegetated + crop) landunits, so we preserve this concept by adjusting pctpft_full so ! that sum(pctpft_full) + pcturb_rescaled = 100%. After doing this, we can rescale ! pctpft_full so it accounts for the other special landunits - but since we have already ! accounted for urban, note that that final rescaling only considers the sum of the ! non-urban special landunits. if (pcturb(n) .gt. 0._r8) then suma = pctlak(n)+pctwet(n)+pctgla(n) if (suma >= 100._r8) then write(6,*) 'ERROR: pcturb > 0 and other special landunits add to >= 100' write(6,*) 'n, pcturb, suma = ', n, pcturb(n), suma call abort() end if ! Rescale pcturb as if it were going on the (vegetated + crop) landunits pcturb_rescaled = pcturb(n)/(0.01_r8*(100._r8 - suma)) ! Replace bare soil preferentially with urban bare_urb_diff = pctpft_full(n,0) - pcturb_rescaled pctpft_full(n,0) = max(0._r8,bare_urb_diff) pcturb_excess = abs(min(0._r8,bare_urb_diff)) ! For any urban not accounted for by bare soil, replace other PFTs ! proportionally sumpft = sum(pctpft_full(n,1:numpft)) if (pcturb_excess > 0._r8 .and. sumpft > 0._r8) then do m = 1, numpft pctpft_full(n,m) = pctpft_full(n,m) - pcturb_excess*pctpft_full(n,m)/sumpft if ( pctpft_full(n,m) < 0.0_r8 )then write (6,*)'pctpft_full < 0.0 = ', pctpft_full(n,m), & ' n, m, suma, pcturb_excess, sumpft = ', n, m, suma, pcturb_excess, sumpft ! Note that the tolerance here (0.00001_r8) matches the ! tolerance for the error check on the sum in mkpftMod: mkpft if ( abs(pctpft_full(n,m)) > 1.e-5_r8 )then call abort() end if pctpft_full(n,m) = 0.0_r8 end if end do end if ! Confirm that we have done the rescaling correctly: ! make sure sum(pctpft_full) + pcturb_rescaled = 100% ! Note that the tolerance here (0.00001_r8) matches the tolerance for the error ! check on the sum in mkpftMod: mkpft if (abs( (sum(pctpft_full(n,0:numpft)) + pcturb_rescaled) - 100._r8 ) > 0.00001_r8) then write(6,*) 'ERROR: pctpft_full + pcturb_rescaled not equal to 100' write(6,*) 'n, sum(pctpft_full), pcturb_rescaled, diff = ' write(6,*) n, sum(pctpft_full(n,0:numpft)), pcturb_rescaled, abs( (sum(pctpft_full(n,0:numpft)) + pcturb_rescaled) - 100._r8 ) call abort() end if end if ! Normalize pctpft_full to be the remainder of [100 - (special landunits)] ! Note that we have already accounted for pcturb, above, so here we just account ! for the other special landunits suma = pctlak(n) + pctwet(n) + pctgla(n) call mkpft_normalize (pctpft_full(n,:), suma, pctnatveg(n), pctcrop(n), pctnatpft(n,:), pctcft(n,:)) end do ! Deallocate no-longer-needed pctpft_full so it won't accidentally be updated further deallocate(pctpft_full) nsmall_tot = 0 do n = 1,ns_o suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctnatveg(n) + pctcrop(n) if (suma < 90._r8) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla, pctnatveg and pctcrop is less than 90' write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n) call abort() else if (suma > 100._r8 + 1.e-4_r8) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla, pctnatveg and pctcrop is greater than 100' write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,sum= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma call abort() else ! If the coverage of any PFT or CFT is too small at the gridcell level, set its ! % cover to 0, then renormalize everything else as needed call remove_small_cover(pctnatveg(n), pctnatpft(n,:), nsmall, toosmallPFT, suma) nsmall_tot = nsmall_tot + nsmall call remove_small_cover(pctcrop(n), pctcft(n,:), nsmall, toosmallPFT, suma) nsmall_tot = nsmall_tot + nsmall suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctnatveg(n) + pctcrop(n) if ( abs(suma - 100.0_r8) > 2.0*epsilon(suma) )then pctlak(n) = pctlak(n) * 100._r8/suma pctwet(n) = pctwet(n) * 100._r8/suma pcturb(n) = pcturb(n) * 100._r8/suma pctgla(n) = pctgla(n) * 100._r8/suma pctnatveg(n) = pctnatveg(n) * 100._r8/suma pctcrop(n) = pctcrop(n) * 100._r8/suma end if end if ! Roundoff error fix suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) suma2 = pctnatveg(n) + pctcrop(n) if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. & (suma2 > 0.0_r8 .and. suma2 < 1.e-6_r8) ) then write (6,*) 'Special land units near 100%, but not quite for n,suma =',n,suma write (6,*) 'Adjusting special land units to 100%' if (pctlak(n) >= 25._r8) then pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n)) else if (pctwet(n) >= 25._r8) then pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n)) else if (pcturb(n) >= 25._r8) then pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n)) else if (pctgla(n) >= 25._r8) then pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n)) else write (6,*) subname, 'Error: sum of special land units nearly 100% but none is >= 25% at ', & 'n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma = ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma call abort() end if pctnatveg(n) = 0._r8 pctcrop(n) = 0._r8 end if if ( any(pctnatpft(n,:)*pctnatveg(n) > 0.0_r8 .and. pctnatpft(n,:)*pctnatveg(n) < toosmallPFT ) .or. & any(pctcft(n,:)*pctcrop(n) > 0.0_r8 .and. pctcft(n,:)*pctcrop(n) < toosmallPFT )) then write (6,*) 'pctnatpft or pctcft is small at n=', n write (6,*) 'pctnatpft = ', pctnatpft(n,:) write (6,*) 'pctcft = ', pctcft(n,:) write (6,*) 'pctnatveg = ', pctnatveg write (6,*) 'pctcrop = ', pctcrop call abort() end if suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) if (suma < 100._r8-epsilon(suma) .and. suma > (100._r8 - 4._r8*epsilon(suma))) then write (6,*) subname, 'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n) call abort() end if suma = suma + pctnatveg(n) + pctcrop(n) if ( abs(suma-100._r8) > 1.e-10_r8) then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla, pctnatveg and pctcrop is NOT equal to 100' write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,sum= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),& pctnatveg(n),pctcrop(n), suma call abort() end if end do ! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100% if ( .not. outnc_double )then do n = 1,ns_o sum8 = real(pctlak(n),r4) sum8 = sum8 + real(pctwet(n),r4) sum8 = sum8 + real(pcturb(n),r4) sum8 = sum8 + real(pctgla(n),r4) sum4a = real(pctnatveg(n),r4) sum4a = sum4a + real(pctcrop(n),r4) if ( sum4a==0.0_r4 .and. sum8 < 100._r4-2._r4*epsilon(sum4a) )then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8 write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n), pctnatveg(n),pctcrop(n) call abort() end if end do else do n = 1,ns_o sum8 = pctlak(n) sum8 = sum8 + pctwet(n) sum8 = sum8 + pcturb(n) sum8 = sum8 + pctgla(n) sum8a = pctnatveg(n) sum8a = sum8a + pctcrop(n) if ( sum8a==0._r8 .and. sum8 < (100._r8-4._r8*epsilon(sum8)) )then write (6,*) subname, ' error: sum of pctlak, pctwet,', & 'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8 write (6,*) 'Total error, error/epsilon = ',100._r8-sum8, ((100._r8-sum8)/epsilon(sum8)) write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,epsilon= ', & n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n), epsilon(sum8) call abort() end if end do end if do n = 1,ns_o do k = natpft_lb,natpft_ub if ( pctnatpft(n,k)*pctnatveg(n) < 0.0_r8 )then write (6,*)'pctnatpft*pctnatveg < 0.0 = ', pctnatpft(n,k)*pctnatveg(n) call abort() end if end do do k = cft_lb,cft_ub if ( pctcft(n,k)*pctcrop(n) < 0.0_r8 )then write (6,*)'pctcft*pctcrop < 0.0 = ', pctcft(n,k)*pctcrop(n) call abort() end if end do end do ! Make sure that there is no vegetation outside the pft mask do n = 1,ns_o if (pftdata_mask(n) == 0 .and. (pctnatveg(n) > 0 .or. pctcrop(n) > 0)) then write (6,*)'vegetation found outside the pft mask at n=',n write (6,*)'pctnatveg,pctcrop=', pctnatveg, pctcrop call abort() end if end do ! Make sure that sums at the landunit level all add to 100% ! (Note that we don't check pctglcmec here, because it isn't computed at the point ! that this subroutine is called -- but the check of sum(pctglcmec) is done in ! mkglcmecMod) do n = 1,ns_o if (abs(sum(pctnatpft(n,:)) - 100._r8) > 1.e-12_r8) then write(6,*) 'sum(pctnatpft(n,:)) != 100: ', n, sum(pctnatpft(n,:)) call abort() end if if (size(pctcft,2) > 0) then if (abs(sum(pctcft(n,:)) - 100._r8) > 1.e-12_r8) then write(6,*) 'sum(pctcft(n,:)) != 100: ', n, sum(pctcft(n,:)) call abort() end if end if if (abs(sum(urbn_classes(n,:)) - 100._r8) > 1.e-12_r8) then write(6,*) 'sum(urbn_classes(n,:)) != 100: ', n, sum(urbn_classes(n,:)) call abort() end if end do if ( nsmall_tot > 0 )then write (6,*)'number of small pft = ', nsmall_tot end if end subroutine normalizencheck_landuse end program mksurfdat