Skip to content

bugfix: fix potential divide by zero in stomatal resistance diagnostic #633

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
barlage wants to merge 3 commits into from
Closed

bugfix: fix potential divide by zero in stomatal resistance diagnostic #633

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
84deb5f
Upgrade WRFDA RTTOV interface to v12.1 (#624)
Sep 13, 2018
3ced231
Upgrade WRFDA CRTM code base to v2.3.0 (#627)
Sep 14, 2018
ac98022
Bug fix introduced in #399 where divide by zero can occur if LAI is zero
Sep 24, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,7 +137,7 @@ all_wrfvar :
$(MAKE) MODULE_DIRS="$(DA_WRFVAR_MODULES)" ext
$(MAKE) MODULE_DIRS="$(DA_WRFVAR_MODULES)" toolsdir
if [ $(CRTM) -ne 0 ] ; then \
(cd var/external/crtm_2.2.3; $(MAKE) $(J)) ; \
(cd var/external/crtm_2.3.0; $(MAKE) $(J)) ; \
fi
if [ $(BUFR) ] ; then \
(cd var/external/bufr; \
Expand Down
71 changes: 38 additions & 33 deletions compile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -335,8 +335,8 @@ else
setenv BUFR 1
endif
setenv CRTM_CPP "-DCRTM"
setenv CRTM_LIB "-L../external/crtm_2.2.3/libsrc -lCRTM"
setenv CRTM_SRC "-I../external/crtm_2.2.3/libsrc"
setenv CRTM_LIB "-L../external/crtm_2.3.0/libsrc -lCRTM"
setenv CRTM_SRC "-I../external/crtm_2.3.0/libsrc"
#setenv SFC_CRTM `grep '^SFC' configure.wrf | awk '{print $3}' | sed -e 's/\// /g' | awk '{print $NF}'`
#setenv ABI_CRTM `grep '^SFC' configure.wrf | sed -n 's/.*\(\-m[0-9]\{2\}\).*/\1/p'`
setenv CRTM 1
Expand All @@ -349,37 +349,6 @@ else
setenv CRTM_SRC " "
setenv CRTM 0
endif
set RTTOV = ( `grep "^RTTOVPATH" configure.wrf | cut -d"=" -f2-` )
if ( $RTTOV == "" ) then
setenv RTTOV_LIB " "
setenv RTTOV_SRC " "
unsetenv RTTOV
else
echo " "
echo "Compiling with RTTOV libraries in:"
echo $RTTOV
echo " "
if ( ! $?BUFR ) then
echo " "
echo "BUFR library is needed for radiance data ingest."
echo "setting BUFR=1"
echo " "
setenv BUFR 1
endif
if ( -e ${RTTOV}/lib/librttov11.1.0_main.a ) then
setenv RTTOV_LIB "-L${RTTOV}/lib -lrttov11.1.0_coef_io -lrttov11.1.0_emis_atlas -lrttov11.1.0_main"
else if ( -e ${RTTOV}/lib/librttov11.2.0_main.a ) then
setenv RTTOV_LIB "-L${RTTOV}/lib -lrttov11.2.0_coef_io -lrttov11.2.0_emis_atlas -lrttov11.2.0_main"
else if ( -e ${RTTOV}/lib/librttov11_main.a ) then
setenv RTTOV_LIB "-L${RTTOV}/lib -lrttov11_coef_io -lrttov11_emis_atlas -lrttov11_main"
else
echo "Can not find a compatible RTTOV library! Please ensure that your RTTOV build was successful,"
echo "your 'RTTOV' environment variable is set correctly, and you are using a supported version of RTTOV."
echo "Currently supported versions are 11.1, 11.2, and 11.3"
exit 1
endif
setenv RTTOV_SRC "-I${RTTOV}/include -I${RTTOV}/mod"
endif
set hdf5path = ( `grep "^HDF5PATH" configure.wrf | cut -d"=" -f2-` )
if ( $hdf5path == "" ) then
setenv HDF5_INC ""
Expand All @@ -392,6 +361,42 @@ else
setenv HDF5_INC "-I${hdf5path}/include"
setenv HDF5 1
endif
set RTTOV = ( `grep "^RTTOVPATH" configure.wrf | cut -d"=" -f2-` )
if ( $RTTOV == "" ) then
setenv RTTOV_LIB " "
setenv RTTOV_SRC " "
unsetenv RTTOV
else
if ( $hdf5path == "" ) then
echo "As of version 12.1 of RTTOV, WRFDA requires HDF5 in order utilize the RTTOV library."
echo "RTTOV emissivity atlas files are now provided only in HDF5 format."
echo "Please supply an HDF5 path prior to configure or unset RTTOV."
exit 1
else
echo " "
echo "Compiling with RTTOV libraries in:"
echo $RTTOV
echo " "
if ( ! $?BUFR ) then
echo " "
echo "BUFR library is needed for radiance data ingest."
echo "setting BUFR=1"
echo " "
setenv BUFR 1
endif
if ( -e ${RTTOV}/lib/librttov12_main.a ) then
setenv RTTOV_LIB "-L${hdf5path}/lib -lhdf5hl_fortran -lhdf5_hl -lhdf5_fortran -lhdf5 -L${RTTOV}/lib -lrttov12_coef_io -lrttov12_emis_atlas -lrttov12_main -lrttov12_hdf"
else
echo "Can not find a compatible RTTOV library! Please ensure that your RTTOV build was successful,"
echo "your 'RTTOV' environment variable is set correctly, and you are using a supported version of RTTOV."
echo "Current supported version(s): 12.1"

exit 1
endif
setenv RTTOV_SRC "-I${RTTOV}/include -I${RTTOV}/mod"
endif
endif

if ( $?CLOUD_CV ) then
setenv CLOUD_CV_CPP "-DCLOUD_CV"
else
Expand Down
2 changes: 1 addition & 1 deletion phys/module_sf_noahmpdrv.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -979,7 +979,7 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN
RB = MAX(RB, 0.0)
! New Calculation of total Canopy/Stomatal Conductance Based on Bonan et al. (2011)
! -- Inverse of Canopy Resistance (below)
IF(RSSUN .le. 0.0 .and. RSSHA .le. 0.0) THEN
IF(RSSUN .le. 0.0 .or. RSSHA .le. 0.0 .or. LAISUN .eq. 0.0 .or. LAISHA .eq. 0.0) THEN
RS (I,J) = 0.0
ELSE
RS (I,J) = ((1.0/(RSSUN+RB)*LAISUN) + ((1.0/(RSSHA+RB))*LAISHA))
Expand Down
27 changes: 13 additions & 14 deletions var/da/da_radiance/da_get_innov_vector_rttov.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ subroutine da_get_innov_vector_rttov (it, grid, ob, iv)
integer :: i, j, k ! Index dimension.
integer :: nlevels ! Number of obs levels.
integer :: nchanprof, errorstatus
integer :: ir_atlas_version, mw_atlas_version

character(len=256) :: atlas_path
real*8 :: seap, icep, lndp, snop
Expand All @@ -42,7 +41,7 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
type(aux_vars_type), allocatable :: aux_vars(:)

type(rttov_chanprof), allocatable :: chanprof(:)
type(profile_type), allocatable :: profiles(:)
type(rttov_profile), allocatable :: profiles(:)

! variables for computing clwp
real, allocatable :: dpf(:,:), clw(:,:), pf(:,:)
Expand Down Expand Up @@ -203,20 +202,19 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
if ( rttov_emis_atlas_ir > 0 .or. rttov_emis_atlas_mw > 0 ) then
! set up emissivity atlas
atlas_path = 'emis_data/'
ir_atlas_version = 100
mw_atlas_version = 100 ! TELSEM
if ( rttov_emis_atlas_mw == 2 ) mw_atlas_version = 200 ! CNRW
write(unit=message(1),fmt='(A,A)') &
'Setting up emissivity atlas for instrument ', trim(iv%instid(inst)%rttovid_string)
call da_message(message(1:1))
call rttov_setup_emis_atlas( &
errorstatus, & ! out
opts(inst), & ! in
grid%start_month, & ! in
coefs(inst), & ! in
atlas_type(inst), & ! in
atlas, & ! inout
atlas_id(inst), & ! in, optional
path = trim(atlas_path), & ! in, optional
ir_atlas_ver = ir_atlas_version, & ! in, optional
mw_atlas_ver = mw_atlas_version) ! in, optional
coefs = coefs(inst)) ! in

if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
(/"failure in setting up emissivity atlas"/))
Expand All @@ -234,11 +232,11 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
! latitude, longitude, surftype are used for retreiving emis from atlas
! zenangle is used by MW emmisivity atlas
! snow_frac is used only by IR emmisivity atlas
profiles(n-n1+1)%latitude = iv%instid(inst)%info%lat(1,n)
profiles(n-n1+1)%longitude = iv%instid(inst)%info%lon(1,n)
profiles(n-n1+1)%zenangle = iv%instid(inst)%satzen(n)
profiles(n-n1+1)%skin%surftype = iv%instid(inst)%surftype(n)
profiles(n-n1+1)%snow_frac = iv%instid(inst)%snow_frac(n)
profiles(n-n1+1)%latitude = iv%instid(inst)%info%lat(1,n)
profiles(n-n1+1)%longitude = iv%instid(inst)%info%lon(1,n)
profiles(n-n1+1)%zenangle = iv%instid(inst)%satzen(n)
profiles(n-n1+1)%skin%surftype = iv%instid(inst)%surftype(n)
profiles(n-n1+1)%skin%snow_fraction = iv%instid(inst)%snow_frac(n)
end do

! Retrieve values from atlas
Expand All @@ -248,6 +246,7 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
chanprof, & ! in
profiles, & ! in
coefs(inst), & ! in
atlas, &! in
emissivity=emissivity(:)%emis_in ) ! out
if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
Expand Down Expand Up @@ -377,7 +376,7 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
deallocate (emissivity)

if ( rttov_emis_atlas_ir > 0 .or. rttov_emis_atlas_mw > 0 ) then
call rttov_deallocate_emis_atlas(coefs(inst))
call rttov_deallocate_emis_atlas(atlas)
end if

end do ! end loop for sensor
Expand Down
4 changes: 2 additions & 2 deletions var/da/da_radiance/da_radiance.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@ module da_radiance
init_constants_derived, &
rttov_platform_name, rttov_inst_name, crtm_sensor_name ! names used by both RTTOV and CRTM
#ifdef RTTOV
use module_radiance, only : coefs, rttov_coefs, profile_type, radiance_type, &
transmission_type,errorstatus_success,gas_id_watervapour,rttov_emissivity
use module_radiance, only : coefs, rttov_coefs, rttov_profile, rttov_radiance, &
rttov_transmission,errorstatus_success,gas_id_watervapour,rttov_emissivity
#endif
#ifdef CRTM
use module_radiance, only : crtm_channelinfo_type, crtm_platform_name, crtm_init, &
Expand Down
7 changes: 4 additions & 3 deletions var/da/da_radiance/da_rttov.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,11 @@ module da_rttov
i_kind,r_kind, r_double, &
one, zero, three,deg2rad, q2ppmv, &
coefs, opts,opts_rt_ir, rttov_inst_name
use module_radiance, only : rttov_options, rttov_opts_rt_ir, rttov_coefs, profile_type, &
transmission_type, radiance_type, rttov_chanprof, &
use module_radiance, only : rttov_options, rttov_opts_rt_ir, rttov_coefs, rttov_profile, &
rttov_transmission, rttov_radiance, rttov_chanprof, &
jpim, jprb, errorstatus_success, errorstatus_fatal, gas_id_watervapour, &
sensor_id_ir, sensor_id_mw, sensor_id_hi,rttov_emissivity
atlas, atlas_type, atlas_id, atlas_type_ir, atlas_type_mw, &
sensor_id_ir, sensor_id_mw, sensor_id_hi, sensor_id_po, rttov_emissivity

use da_control, only : max_ob_levels,missing_r, &
v_interp_p, v_interp_h, tovs_batch, gravity, &
Expand Down
33 changes: 18 additions & 15 deletions var/da/da_radiance/da_rttov_ad.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,16 +20,16 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &

! RTTOV input parameters
type (rttov_chanprof), allocatable :: chanprof(:)
type (profile_type), allocatable :: profiles(:), profiles_ad(:)
type (rttov_profile), allocatable :: profiles(:), profiles_ad(:)
type (rttov_emissivity), allocatable :: emissivity(:), emissivity_ad(:)
logical, allocatable :: calcemis(:)

! RTTOV out parameters
integer :: errorstatus

! RTTOV inout parameters
type (radiance_type) :: radiance, radiance_ad
type (transmission_type) :: transmission, transmission_ad
type (rttov_radiance) :: radiance, radiance_ad
type (rttov_transmission) :: transmission, transmission_ad

call da_trace_entry("da_rttov_ad")

Expand Down Expand Up @@ -94,7 +94,7 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
profiles(n) % cloud(:,:) = 0.0
profiles(n) % cfrac(:) = 0.0
profiles(n) % idg = 1
profiles(n) % ish = 1
profiles(n) % ice_scheme = 1
end if

profiles(n) % skin % surftype = aux_vars (n) % surftype
Expand All @@ -106,7 +106,8 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
end if
end if

if ( coefs(inst)%coef%id_sensor == sensor_id_mw ) then
if ( coefs(inst)%coef%id_sensor == sensor_id_mw .or. &
coefs(inst)%coef%id_sensor == sensor_id_po ) then
if ( profiles(n) % skin % surftype == 2 ) then
profiles(n) % skin % fastem (1) = 2.2
profiles(n) % skin % fastem (2) = 3.7
Expand Down Expand Up @@ -164,7 +165,7 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
call rttov_alloc_transmission( &
& errorstatus, &
& transmission, &
& nlevels-1, &
& nlevels, &
& nchanprof, &
& asw, &
& init = .true. )
Expand All @@ -177,7 +178,7 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
call rttov_alloc_transmission( &
& errorstatus, &
& transmission_ad, &
& nlevels-1, &
& nlevels, &
& nchanprof, &
& asw, &
& init = .true. )
Expand All @@ -192,7 +193,7 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
& errorstatus, &
& nchanprof, &
& radiance, &
& nlevels-1, &
& nlevels, &
& asw, &
& init=.true. )
if ( errorstatus /= errorstatus_success ) then
Expand All @@ -205,7 +206,7 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
& errorstatus, &
& nchanprof, &
& radiance_ad, &
& nlevels-1, &
& nlevels, &
& asw, &
& init=.true. )
if ( errorstatus /= errorstatus_success ) then
Expand All @@ -218,11 +219,13 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
radiance_ad % clear ((n-1)*nchanl+1:n*nchanl) = 0.0
end do

if (coefs(inst)%coef%id_sensor == 1 .or. coefs(inst)%coef%id_sensor == 3) then ! infrared sensor
if ( coefs(inst)%coef%id_sensor == sensor_id_ir .or. &
coefs(inst)%coef%id_sensor == sensor_id_hi ) then ! infrared sensor
calcemis(1:nchanprof) = .true.
emissivity(1:nchanprof)%emis_in = 0.0
emissivity_ad(1:nchanprof)%emis_in = 0.0
else if (coefs(inst)%coef%id_sensor == 2) then ! microwave sensor
else if ( coefs(inst)%coef%id_sensor == sensor_id_mw .or. &
coefs(inst)%coef%id_sensor == sensor_id_po ) then ! microwave sensor
do n = 1, nprofiles
if ( profiles(n) % skin % surftype == 1) then ! sea
calcemis((n-1)*nchanl+1:n*nchanl) = .true.
Expand Down Expand Up @@ -291,24 +294,24 @@ subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &

asw = 0 ! deallocation
! deallocate radiance arrays
call rttov_alloc_rad (errorstatus,nchanprof,radiance,nlevels-1,asw)
call rttov_alloc_rad (errorstatus,nchanprof,radiance,nlevels,asw)
if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
(/"radiance deallocation error"/))
end if
call rttov_alloc_rad (errorstatus,nchanprof,radiance_ad,nlevels-1,asw)
call rttov_alloc_rad (errorstatus,nchanprof,radiance_ad,nlevels,asw)
if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
(/"radiance AD deallocation error"/))
end if

! deallocate transmission arrays
call rttov_alloc_transmission (errorstatus,transmission,nlevels-1,nchanprof,asw)
call rttov_alloc_transmission (errorstatus,transmission,nlevels,nchanprof,asw)
if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
(/"transmission deallocation error"/))
end if
call rttov_alloc_transmission (errorstatus,transmission_ad,nlevels-1,nchanprof,asw)
call rttov_alloc_transmission (errorstatus,transmission_ad,nlevels,nchanprof,asw)
if ( errorstatus /= errorstatus_success ) then
call da_error(__FILE__,__LINE__, &
(/"transmission AD deallocation error"/))
Expand Down
Loading