Dynamical lakes and new lake cover map #1073
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -64,7 +64,7 @@ program mksurfdat | |
| integer :: k,m,n ! indices | ||
| integer :: ni,nj,ns_o ! indices | ||
| integer :: ier ! error status | ||
| integer :: ndiag,nfdyn,nfdynlak ! unit numbers | ||
| integer :: ncid ! netCDF id | ||
| integer :: omode ! netCDF output mode | ||
| integer :: varid ! netCDF variable id | ||
|
|
@@ -80,6 +80,7 @@ program mksurfdat | |
| character(len=256) :: fdyndat ! dynamic landuse data file name | ||
| character(len=256) :: fname ! generic filename | ||
| character(len=256) :: fhrvname ! generic harvest filename | ||
| character(len=256) :: flakname ! generic lake filename | ||
| character(len=256) :: string ! string read in | ||
| integer :: t1 ! timer | ||
| real(r8),parameter :: p5 = 0.5_r8 ! constant | ||
|
|
@@ -134,7 +135,10 @@ program mksurfdat | |
| 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 :: pctlak_orig(:) ! percent lake of gridcell before dynamic land use adjustments | ||
| real(r8), allocatable :: pctwet_orig(:) ! percent wetland of gridcell before dynamic land use adjustments | ||
| real(r8), allocatable :: pcturb_orig(:) ! percent urban of gridcell before dynamic land use adjustments | ||
| real(r8), allocatable :: pctgla_orig(:) ! percent glacier of gridcell before dynamic land use adjustments | ||
| real(r8) :: std_elev = -999.99_r8 ! Standard deviation of elevation (m) to use for entire grid | ||
|
|
||
| integer, allocatable :: harvind1D(:) ! Indices of 1D harvest fields | ||
|
|
@@ -150,11 +154,11 @@ program mksurfdat | |
| type(harvestDataType) :: harvdata | ||
|
|
||
| namelist /clmexp/ & | ||
| mksrf_fgrid, & | ||
| mksrf_gridtype, & | ||
| mksrf_fvegtyp, & | ||
| mksrf_fhrvtyp, & | ||
| mksrf_fsoitex, & | ||
| mksrf_forganic, & | ||
| mksrf_fsoicol, & | ||
| mksrf_fvocef, & | ||
|
|
@@ -167,6 +171,7 @@ program mksurfdat | |
| mksrf_furban, & | ||
| mksrf_flai, & | ||
| mksrf_fdynuse, & | ||
| mksrf_fdynlak, & | ||
| mksrf_fgdp, & | ||
| mksrf_fpeat, & | ||
| mksrf_fsoildepth, & | ||
|
|
@@ -278,6 +283,7 @@ program mksurfdat | |
| ! Optionally specify setting for: | ||
| ! ====================================== | ||
| ! mksrf_fdynuse ----- ASCII text file that lists each year of pft files to use | ||
| ! mksrf_fdynlak ----- ASCII text file that list each year of dynlake 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 | ||
|
|
@@ -459,7 +465,11 @@ program mksurfdat | |
| vic_dsmax(ns_o) , & | ||
| vic_ds(ns_o) , & | ||
| lakedepth(ns_o) , & | ||
| glacier_region(ns_o) , & | ||
| pctlak_orig(ns_o) , & | ||
| pctwet_orig(ns_o) , & | ||
| pcturb_orig(ns_o) , & | ||
| pctgla_orig(ns_o) ) | ||
| landfrac_pft(:) = spval | ||
| pctlnd_pft(:) = spval | ||
| pftdata_mask(:) = -999 | ||
|
|
@@ -709,62 +719,12 @@ program mksurfdat | |
|
|
||
| call change_landuse( ldomain, dynpft=.false. ) | ||
|
|
||
| ! Save special land unit areas of surface dataset | ||
| pctlak_orig(:) = pctlak(:) | ||
| pctwet_orig(:) = pctwet(:) | ||
| pcturb_orig(:) = pcturb(:) | ||
| pctgla_orig(:) = pctgla(:) | ||
|
|
||
| call normalizencheck_landuse(ldomain) | ||
|
|
||
| ! Write out sum of PFT's | ||
|
|
@@ -1062,14 +1022,11 @@ program mksurfdat | |
|
|
||
| ! Deallocate arrays NOT needed for dynamic-pft section of code | ||
|
|
||
| deallocate ( ef1_btr, ef1_fet, ef1_fdt, ef1_shr, ef1_grs, ef1_crp ) | ||
| deallocate ( pctglcmec, topoglcmec) | ||
| if ( outnc_3dglc ) deallocate ( pctglc_gic, pctglc_icesheet) | ||
| deallocate ( elevclass ) | ||
| deallocate ( fmax ) | ||
| deallocate ( gdp, fpeat, agfirepkmon ) | ||
| deallocate ( soildepth ) | ||
| deallocate ( topo_stddev, slope ) | ||
|
|
@@ -1118,7 +1075,7 @@ program mksurfdat | |
| 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) | ||
|
|
||
| ! Synchronize the disk copy of a netCDF dataset with in-memory buffers | ||
|
|
||
|
|
@@ -1127,6 +1084,9 @@ program mksurfdat | |
| ! Read in each dynamic pft landuse dataset | ||
|
|
||
| nfdyn = getavu(); call opnfil (mksrf_fdynuse, nfdyn, 'f') | ||
|
|
||
| ! Read in dynamic lake dataset | ||
| nfdynlak = getavu(); call opnfil (mksrf_fdynlak, nfdynlak, 'f') | ||
|
|
||
| pctnatpft_max = pctnatpft | ||
| pctcft_max = pctcft | ||
|
|
@@ -1159,16 +1119,26 @@ program mksurfdat | |
| call abort() | ||
| end if | ||
| end if | ||
|
|
||
|
|
||
| ! Read input lake pct data | ||
| read(nfdynlak, '(A195,1x,I4)', iostat=ier) string, year | ||
|
There was a problem hiding this comment.
There was a problem hiding this comment.
How exactly we handle this will depend on whether we're going to put transient lake areas on all of our landuse timeseries files moving forward. @ekluzek do you feel we should just go ahead and do that, assuming we have the data? I'm thinking yes, but do you have thoughts on it? |
||
| if (ier /= 0) exit | ||
|
|
||
| flakname = string | ||
| write(6,*)'input lake dynamic dataset for year ', year, ' is : ', trim(flakname) | ||
|
|
||
|
|
||
| ntim = ntim + 1 | ||
|
|
||
| ! Create pctpft data at model resolution | ||
|
|
||
| call mkpft(ldomain, mapfname=map_fpft, fpft=fname, & | ||
| ndiag=ndiag, pctlnd_o=pctlnd_pft_dyn, pctnatpft_o=pctnatpft, pctcft_o=pctcft ) | ||
| ! Create harvesting data at model resolution | ||
|
|
||
| call mkharvest(ldomain, mapfname=map_fharvest, datfname=fhrvname, & | ||
| ndiag=ndiag, harvdata=harvdata ) | ||
|
|
||
| ! Consistency check on input land fraction | ||
|
|
@@ -1186,9 +1156,17 @@ program mksurfdat | |
| call abort() | ||
| end if | ||
| end do | ||
|
|
||
|
|
||
|
There was a problem hiding this comment.
|
||
| ! Create pctlak data at model resolution (use original mapping file from lake data) | ||
| call mklakwat (ldomain, mapfname=map_flakwat, datfname=flakname, & | ||
| ndiag=ndiag, zero_out=all_urban.or.all_veg, lake_o=pctlak) | ||
|
|
||
| ! For landunits NOT read each year: reset to their pre-adjustment values in preparation for redoing landunit area normalization | ||
| pctwet(:) = pctwet_orig(:) | ||
| pcturb(:) = pcturb_orig(:) | ||
| pctgla(:) = pctgla_orig(:) | ||
|
|
||
| call normalizencheck_landuse(ldomain) | ||
|
|
||
| call update_max_array(pctnatpft_max,pctnatpft) | ||
|
|
@@ -1207,6 +1185,9 @@ program mksurfdat | |
| call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctcft)) | ||
| end if | ||
|
|
||
| call check_ret(nf_inq_varid(ncid, 'PCT_LAKE', varid), subname) | ||
| call ncd_put_time_slice(ncid, varid, ntim, pctlak) | ||
|
|
||
| call harvdata%getFieldsIdx( harvind1D, harvind2D ) | ||
| do k = 1, harvdata%num1Dfields() | ||
| call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(harvind1D(k),constant=.false.)), varid), subname) | ||
|
|
@@ -1370,13 +1351,66 @@ subroutine normalizencheck_landuse(ldomain) | |
| character(len=32) :: subname = 'normalizencheck_landuse' ! subroutine name | ||
| !----------------------------------------------------------------------- | ||
|
|
||
| ! ------------------------------------------------------------------------------- | ||
| ! Normalize vegetated and lake area so that vegetated + other special area is 100% | ||
| ! ------------------------------------------------------------------------------- | ||
|
|
||
| ns_o = ldomain%ns | ||
| do n = 1,ns_o | ||
|
|
||
| ! 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))) | ||
| end do | ||
| pctlak(n) = float(nint(pctlak(n))) | ||
| pctwet(n) = float(nint(pctwet(n))) | ||
| pctgla(n) = float(nint(pctgla(n))) | ||
|
|
||
| ! Assume wetland, glacier and/or lake when dataset landmask implies ocean | ||
| ! (assume medium soil color (15) and loamy texture). | ||
| ! Also set pftdata_mask here | ||
|
|
||
| if (pctlnd_pft(n) < 1.e-6_r8) then | ||
| pftdata_mask(n) = 0 | ||
| soicol(n) = 15 | ||
| if (pctgla(n) < 1.e-6_r8) then | ||
| pctwet(n) = 100._r8 - pctlak(n) | ||
| pctgla(n) = 0._r8 | ||
| else | ||
| pctwet(n) = max(100._r8 - pctgla(n) - pctlak(n), 0.0_r8) | ||
| end if | ||
| pcturb(n) = 0._r8 | ||
| call pctnatpft(n)%set_pct_l2g(0._r8) | ||
| call pctcft(n)%set_pct_l2g(0._r8) | ||
| pctsand(n,:) = 43._r8 | ||
| pctclay(n,:) = 18._r8 | ||
| organic(n,:) = 0._r8 | ||
| else | ||
| pftdata_mask(n) = 1 | ||
| end if | ||
|
|
||
| ! 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 | ||
|
|
||
|
|
||
| ! Check preconditions | ||
| if ( pctlak(n) < 0.0_r8 )then | ||
| write(6,*) subname, ' ERROR: pctlak is negative!' | ||
|
|
@@ -1398,7 +1432,7 @@ subroutine normalizencheck_landuse(ldomain) | |
| write(6,*) 'n, pctgla = ', n, pctgla(n) | ||
| call abort() | ||
| end if | ||
| suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) | ||
| if (suma > (100._r8 + tol_loose)) then | ||
| write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be' | ||
|
|
@@ -1408,6 +1442,7 @@ subroutine normalizencheck_landuse(ldomain) | |
| call abort() | ||
| end if | ||
|
|
||
|
|
||
| ! First normalize vegetated (natural veg + crop) cover so that the total of | ||
| ! (vegetated + (special excluding urban)) is 100%. We'll deal with urban later. | ||
| ! | ||
|
|
@@ -1416,7 +1451,7 @@ subroutine normalizencheck_landuse(ldomain) | |
| ! will work properly regardless of the initial area of natural veg + crop (even if | ||
| ! that initial area is 0%). | ||
|
|
||
| suma = pctlak(n) + pctwet(n)+ pctgla(n) | ||
| new_total_veg_pct = 100._r8 - suma | ||
| ! correct for rounding error: | ||
| new_total_veg_pct = max(new_total_veg_pct, 0._r8) | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What happens if the dynamic lake file isn't provided? We need a way for it to abort. We also should decide on if lakes should always be dynamic for transient conditions, or if we you pick and choose? It would probably be best to always have it on the datasets, but then decide in CTSM whether to activate it or not. We have something similar for harvest.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As I mention in a comment below, I'm thinking we should put the lake file names on the same file as the pft & harvest file names. @ekluzek does that make sense to you, or do you see a reason to keep them on separate files? My assumption (like yours, @ekluzek ) is that we'll always put the transient lake data on the landuse timeseries files, but this will only be activated if the flag
do_transient_lakesis true. By default I think this will be false for now, while this is still an experimental setting.