From 616f437b1012654d40558b928a42126c1f00aab2 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 2 Nov 2017 13:30:36 -0600 Subject: [PATCH 1/2] 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. --- .../runoff_to_ocn/src/map_mod.F90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 b/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 index 7cf87a7d846..4509f21cd04 100644 --- a/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 +++ b/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 @@ -644,13 +644,23 @@ 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 (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: ", trim(units) + stop + end if map%frac_a = map%mask_a * 1.0_r8 From 2a8a76d02ea5643166619c1fd253ad7a30713c76 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 2 Nov 2017 17:05:16 -0600 Subject: [PATCH 2/2] Abort if no units are found for grid_area Also give a more helpful message if grid_area has unrecognized units --- .../gen_mapping_files/runoff_to_ocn/src/map_mod.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 b/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 index 4509f21cd04..1a7dc51ac99 100644 --- a/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 +++ b/tools/mapping/gen_mapping_files/runoff_to_ocn/src/map_mod.F90 @@ -653,12 +653,18 @@ SUBROUTINE map_gridRead(map, rfilename, ofilename, gridtype, lmake_rSCRIP) 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: ", trim(units) + 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