Skip to content

Commit

Permalink
Merge pull request #2022 from billsacks/runoff_convert_sq_degrees
Browse files Browse the repository at this point in the history
If area units on source grid file are in sq deg, convert to sq radians

This just applies to source grid files in the scrip format.

Without this fix, areas were assumed to be in sq radians, which led to
mapping weights that were off by about a factor of 3000 if areas were
actually in sq degrees.

Test suite: Tested runoff_map tool with a source file with units in sq
            degrees; confirmed that area_a and S values now have the
            right order of magnitude. Also tested with a source file
            with units in sq radians, and confirmed that the resulting
            mapping file is the same as before.
Test baseline: n/a
Test namelist changes: none
Test status: bit for bit

Fixes #2018 

User interface changes?: none

Update gh-pages html (Y/N)?: N

Code review:
  • Loading branch information
alperaltuntas authored Nov 3, 2017
2 parents 4793fd8 + 2a8a76d commit 28f74aa
Showing 1 changed file with 19 additions and 3 deletions.
22 changes: 19 additions & 3 deletions tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -644,13 +644,29 @@ SUBROUTINE map_gridRead(map, rfilename, ofilename, gridtype, lmake_rSCRIP)

rcode = nf_inq_varid (fid,'grid_imask',vid )
rcode = nf_get_var_int (fid,vid ,map%mask_a)

rcode = nf_inq_varid (fid,'grid_area',vid )
if (rcode.eq.0) then
rcode = nf_get_var_double(fid,vid ,map%area_a)
else
if (rcode.ne.0) then
write(6,*) "ERROR: could not find variable grid_area in source grid input file!"
stop
end if
rcode = nf_get_var_double(fid,vid ,map%area_a)
units = "" ! units needs to be emptied before reading from netCDF file
rcode = nf_get_att_text(fid, vid, "units", units)
if (rcode.ne.0) then
write(6,*) "ERROR: No units attribute found for source grid_area variable"
write(6,*) "Please add a units attribute with value 'square radians' or 'square degrees'"
stop
end if
if (trim(units).eq."square radians") then
! Nothing to do
else if (trim(units).eq."square degrees") then
map%area_a = map%area_a * DEGtoRAD * DEGtoRAD
else
write(6,*) "ERROR: Unrecognized units for source grid_area variable: ", trim(units)
write(6,*) "Recognized units are 'square radians' or 'square degrees'"
stop
end if

map%frac_a = map%mask_a * 1.0_r8

Expand Down

0 comments on commit 28f74aa

Please sign in to comment.