Skip to content

Commit

Permalink
Update 5-band dEdd to support test data (#472)
Browse files Browse the repository at this point in the history 
* Update 5-band dEdd to support test data

Add an interpolation method in icepack_shortwave.F90 for rsnw

Refactor portions of icepack_shortwave to improve robustness

Add snicar and snicartest test cases

* Update the 5-band snicar shortwave rsnw_snicar_tab interpolation

This changes answers but QC tests pass.  The new implementation checks
whether the rsnw_snicar_tab is an array of integer values with deltas
of value 1 (and sets the rsnw_snicar_chk variable).  If so, it uses a
shortcut to find the rsnw_snicar_tab bounds for interpolation, otherwise
it calls the shortwave_search routine.  In testing, the shortcut did not
change performance much, but that could change in the future.

* Clean up debug output

* Refactor the logic to control rsnw table interpolation in the shortwave.
Set a character string, rsnw_datatype, at initialization.
  • Loading branch information
@apcraig
apcraig authored Oct 27, 2023
1 parent 7a0d4e9 commit 6ad0f44
Show file tree
Hide file tree
Showing 6 changed files with 162 additions and 59 deletions.
189 changes: 146 additions & 43 deletions columnphysics/icepack_shortwave.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ module icepack_shortwave
use icepack_tracers, only: tr_zaero, nlt_chl_sw, nlt_zaero_sw
use icepack_tracers, only: n_algae, n_aero, n_zaero
use icepack_tracers, only: nmodal1, nmodal2, max_aero
use icepack_shortwave_data, only: nspint_3bd, nspint_5bd
use icepack_shortwave_data, only: nspint_3bd, nspint_5bd, rsnw_datatype
use icepack_zbgc_shared,only: R_chl2N, F_abs_chl
use icepack_zbgc_shared,only: remap_zbgc
use icepack_orbital, only: compute_coszen
Expand Down Expand Up @@ -101,9 +101,8 @@ module icepack_shortwave
gaer_5bd, kaer_5bd, waer_5bd
use icepack_shortwave_data, only: &
nmbrad_snicar , & ! number of snow grain radii in SNICAR SSP tables
rsnw_snicar_min, & ! minimum snow radius - integer value used for indexing
rsnw_snicar_max ! maximum snow radius - integer value used for indexing
use icepack_shortwave_data, only: &
rsnw_snicar_min, & ! minimum snow radius
rsnw_snicar_max, & ! maximum snow radius
ssp_snwextdr, ssp_snwalbdr, ssp_sasymmdr, &
ssp_snwextdf, ssp_snwalbdf, ssp_sasymmdf, &
rsnw_snicar_tab
Expand Down Expand Up @@ -2135,7 +2134,7 @@ subroutine compute_dEdd_3bd( &
! Cheng: note that aerosol IOPs are related to snow grain radius.
! CICE adjusted snow grain radius rsnw to frsnw in the original 3-band
! scheme, while for SNICAR the snow grain radius is used directly.
ksnow = k - min(k-1,0)
ksnow = max(k,1)
tmp_gs = frsnw(ksnow)

! grain size index
Expand Down Expand Up @@ -2268,7 +2267,7 @@ subroutine compute_dEdd_3bd( &

do k = 0, nslyr
! use top rsnw, rhosnw for snow ssl and rest of top layer
ksnow = k - min(k-1,0)
ksnow = max(k,1)
! find snow iops using input snow density and snow grain radius:
if (frsnw(ksnow) < rsnw_tab(1)) then
Qs = Qs_tab(ns,1)
Expand All @@ -2279,20 +2278,16 @@ subroutine compute_dEdd_3bd( &
ws = ws_tab(ns,nmbrad_snw)
gs = gs_tab(ns,nmbrad_snw)
else
! linear interpolation in rsnw
do nr = 2, nmbrad_snw
if (rsnw_tab(nr-1) <= frsnw(ksnow) .and. &
frsnw(ksnow) < rsnw_tab(nr)) then
delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / &
(rsnw_tab(nr) - rsnw_tab(nr-1))
Qs = Qs_tab(ns,nr-1)*(c1-delr) + &
Qs_tab(ns,nr )* delr
ws = ws_tab(ns,nr-1)*(c1-delr) + &
ws_tab(ns,nr )* delr
gs = gs_tab(ns,nr-1)*(c1-delr) + &
gs_tab(ns,nr )* delr
endif
enddo ! nr
call shortwave_search(frsnw(ksnow),rsnw_tab,nr)
if (icepack_warnings_aborted(subname)) return
delr = (frsnw(ksnow) - rsnw_tab(nr-1)) / &
(rsnw_tab(nr) - rsnw_tab(nr-1))
Qs = Qs_tab(ns,nr-1)*(c1-delr) + &
Qs_tab(ns,nr )* delr
ws = ws_tab(ns,nr-1)*(c1-delr) + &
ws_tab(ns,nr )* delr
gs = gs_tab(ns,nr-1)*(c1-delr) + &
gs_tab(ns,nr )* delr
endif
ks = Qs*((rhosnw(ksnow)/rhoi)*3._dbl_kind / &
(4._dbl_kind*frsnw(ksnow)*1.0e-6_dbl_kind))
Expand Down Expand Up @@ -4677,7 +4672,7 @@ subroutine compute_dEdd_5bd( &
! CICE adjusted snow grain radius rsnw to frsnw, while for
! SNICAR there is no need, the tmp_gs is therefore calculated
! differently from code in subroutine compute_dEdd
ksnow = k - min(k-1,0)
ksnow = max(k,1)
tmp_gs = rsnw(ksnow) ! use rsnw not frsnw

! grain size index
Expand Down Expand Up @@ -4806,7 +4801,7 @@ subroutine compute_dEdd_5bd( &
if (nsky == 1) then ! direct incident
do k = 0, nslyr
! use top rsnw, rhosnw for snow ssl and rest of top layer
ksnow = k - min(k-1,0)
ksnow = max(k,1)
if (rsnw(ksnow) <= rsnw_snicar_min) then
ks = ext_cff_mss_ice_drc(ns,1)
ws = ss_alb_ice_drc (ns,1)
Expand All @@ -4815,16 +4810,18 @@ subroutine compute_dEdd_5bd( &
ks = ext_cff_mss_ice_drc(ns,nmbrad_snicar)
ws = ss_alb_ice_drc (ns,nmbrad_snicar)
gs = asm_prm_ice_drc (ns,nmbrad_snicar)
elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then
! radius = 30 --> nr = 1 in SNICAR table ! NOTE
nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30
ks = ext_cff_mss_ice_drc(ns,nr)
ws = ss_alb_ice_drc (ns,nr)
gs = asm_prm_ice_drc (ns,nr)
else ! linear interpolation in rsnw
nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30
delr = (rsnw(ksnow) - floor(rsnw(ksnow))) & ! hardwired delta radius = 1 in table
/ (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow))) ! denom always = 1
else
! linear interpolation
if (trim(rsnw_datatype) == 'sorted_idelta1') then
! NOTE: Assumes delta rsnw_snicar_tab is 1 and rsnw_snicar_tab are integers
! This is just for performance, could call shortwave_search
nr = ceiling(rsnw(ksnow)) - nint(rsnw_snicar_min) + 1
else
call shortwave_search(rsnw(ksnow),rsnw_snicar_tab,nr)
if (icepack_warnings_aborted(subname)) return
endif
delr = (rsnw(ksnow) - rsnw_snicar_tab(nr-1)) &
/ (rsnw_snicar_tab(nr) - rsnw_snicar_tab(nr-1))
ks = ext_cff_mss_ice_drc(ns,nr-1)*(c1-delr) &
+ ext_cff_mss_ice_drc(ns,nr )* delr
ws = ss_alb_ice_drc (ns,nr-1)*(c1-delr) &
Expand All @@ -4839,7 +4836,7 @@ subroutine compute_dEdd_5bd( &
elseif (nsky == 2) then ! diffuse incident
do k = 0, nslyr
! use top rsnw, rhosnw for snow ssl and rest of top layer
ksnow = k - min(k-1,0)
ksnow = max(k,1)
if (rsnw(ksnow) < rsnw_snicar_min) then
ks = ext_cff_mss_ice_dfs(ns,1)
ws = ss_alb_ice_dfs (ns,1)
Expand All @@ -4848,16 +4845,18 @@ subroutine compute_dEdd_5bd( &
ks = ext_cff_mss_ice_dfs(ns,nmbrad_snicar)
ws = ss_alb_ice_dfs (ns,nmbrad_snicar)
gs = asm_prm_ice_dfs (ns,nmbrad_snicar)
elseif (ceiling(rsnw(ksnow)) - rsnw(ksnow) < 1.0e-3_dbl_kind) then
! radius = 30 --> nr = 1 in SNICAR table ! NOTE
nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30
ks = ext_cff_mss_ice_dfs(ns,nr)
ws = ss_alb_ice_dfs (ns,nr)
gs = asm_prm_ice_dfs (ns,nr)
else ! linear interpolation in rsnw
nr = ceiling(rsnw(ksnow)) - 30 + 1 ! hardwired min radius = 30
delr = (rsnw(ksnow) - floor(rsnw(ksnow))) & ! hardwired delta radius = 1 in table
/ (ceiling(rsnw(ksnow)) - floor(rsnw(ksnow)))
else
! linear interpolation
if (trim(rsnw_datatype) == 'sorted_idelta1') then
! NOTE: delta rsnw_snicar_tab is 1 and rsnw_snicar_tab are integers
! This is just for performance, could call shortwave_search
nr = ceiling(rsnw(ksnow)) - nint(rsnw_snicar_min) + 1
else
call shortwave_search(rsnw(ksnow),rsnw_snicar_tab,nr)
if (icepack_warnings_aborted(subname)) return
endif
delr = (rsnw(ksnow) - rsnw_snicar_tab(nr-1)) &
/ (rsnw_snicar_tab(nr) - rsnw_snicar_tab(nr-1))
ks = ext_cff_mss_ice_dfs(ns,nr-1)*(c1-delr) &
+ ext_cff_mss_ice_dfs(ns,nr )* delr
ws = ss_alb_ice_dfs (ns,nr-1)*(c1-delr) &
Expand Down Expand Up @@ -5400,6 +5399,110 @@ subroutine compute_dEdd_5bd( &

end subroutine compute_dEdd_5bd

!=======================================================================
! This subroutine searches array for val and returns nr such that
! array(nr-1) < val <= array(nr)
! If nr cannot be found, an error is thrown
! This does NOT check that array is sorted because it would be too expensive,
! but it must be sorted to work properly.

subroutine shortwave_search(val,array,nr)

real (kind=dbl_kind), intent(in) :: &
val ! search value

real (kind=dbl_kind), dimension (:), intent(in) :: &
array ! sorted array

integer (kind=int_kind), intent(out) :: &
nr ! index in array >= val

! local variables

integer (kind=int_kind) :: &
nrcnt, & ! counter
nrp, & ! prior nr
nrl, nru, & ! lower and upper search indices
nrsize ! size of array

logical (kind=log_kind) :: &
found ! search flag

character (len=512) :: &
tmpstr ! temporary string

character(len=*),parameter :: subname='(shortwave_search)'


if (rsnw_datatype(1:6) /= 'sorted') then
call icepack_warnings_add(subname//' rsnw_datatype not valid: '//trim(rsnw_datatype))
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
endif

nrsize = size(array)

!debug write(tmpstr,*) "val = ",val
! call icepack_warnings_add(subname//trim(tmpstr))
! write(tmpstr,*) "nrsize = ",nrsize
! call icepack_warnings_add(subname//trim(tmpstr))
! write(tmpstr,*) "array1 = ",array(1)
! call icepack_warnings_add(subname//trim(tmpstr))
! write(tmpstr,*) "arrayn = ",array(nrsize)
! call icepack_warnings_add(subname//trim(tmpstr))

if (nrsize > 10) then
! binary search
nrl = 1
nru = nrsize
nr = (nrl + nru) / 2
found = .false.
nrcnt = 0
do while (.not.found .and. nrcnt < nrsize)
nrcnt = nrcnt + 1
nrp = nr
if (val > array(nr)) then
if (val < array(nr+1)) then
found = .true.
nr = nr + 1
else
nrl = nr + 1
nr = (nrl + nru) / 2
endif
else
if (val > array(nr-1)) then
found = .true.
else
nru = nr - 1
nr = (nrl + nru) / 2
endif
endif
!debug write(tmpstr,*) "iter = ",nrcnt,nrp,nr
! call icepack_warnings_add(subname//trim(tmpstr))
! call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
enddo
if (.not. found) then
call icepack_warnings_add(subname//' ERROR: binary search failed')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
else
! linear search
nr = -1
do nrcnt = 2,nrsize
if (val > array(nrcnt-1) .and. val < array(nrcnt)) then
nr = nrcnt
exit
endif
enddo
if (nr < 1) then
call icepack_warnings_add(subname//' ERROR: linear search failed')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
endif

end subroutine shortwave_search

!=======================================================================

end module icepack_shortwave
Expand Down
17 changes: 12 additions & 5 deletions columnphysics/icepack_shortwave_data.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,9 @@ module icepack_shortwave_data
nspint_3bd = 3, & ! number of solar spectral bands
nspint_5bd = 5 ! number of solar spectral bands (snicar snow)

character (len=char_len), public :: &
rsnw_datatype = 'unknown'

! dEdd 3-band data
real (kind=dbl_kind), dimension (nspint_3bd), public :: &
! inherent optical properties (iop)
Expand Down Expand Up @@ -63,7 +66,9 @@ module icepack_shortwave_data

! dEdd 5-band data SSP SNICAR
integer (kind=int_kind), public :: &
nmbrad_snicar , & ! number of snow grain radii in SNICAR SSP tables
nmbrad_snicar ! number of snow grain radii in SNICAR SSP tables

real (kind=dbl_kind), public :: &
rsnw_snicar_min, & ! minimum snow radius - integer value used for indexing
rsnw_snicar_max ! maximum snow radius - integer value used for indexing

Expand Down Expand Up @@ -107,6 +112,7 @@ subroutine icepack_shortwave_init_dEdd3band
allocate(gs_tab (nspint_3bd,nmbrad_snw))

! snow grain radii (micro-meters) for table
rsnw_datatype = 'sorted'
rsnw_tab = (/ & ! snow grain radius for each table entry (micro-meters)
5._dbl_kind, 7._dbl_kind, 10._dbl_kind, 15._dbl_kind, &
20._dbl_kind, 30._dbl_kind, 40._dbl_kind, 50._dbl_kind, &
Expand Down Expand Up @@ -646,6 +652,7 @@ subroutine icepack_shortwave_init_snicartest()
6._dbl_kind, 37._dbl_kind, 221._dbl_kind, 600._dbl_kind, 1340._dbl_kind/)
rsnw_snicar_min = rsnw_snicar_tab(1) ! minimum snow radius - integer value used for indexing
rsnw_snicar_max = rsnw_snicar_tab(nmbrad_snicar) ! maximum snow radius - integer value used for indexing
rsnw_datatype = 'sorted'

allocate(ssp_snwextdr(nspint_5bd,nmbrad_snicar)) ! extinction coefficient, direct
allocate(ssp_snwextdf(nspint_5bd,nmbrad_snicar)) ! extinction coefficient, diffuse
Expand Down Expand Up @@ -9586,11 +9593,11 @@ subroutine icepack_shortwave_init_snicar()

! Copy to local variables

!echmod - this might not be needed
rsnw_snicar_min = 30 ! minimum snow grain radius
rsnw_snicar_max = 1500 ! maximum snow grain radius
rsnw_snicar_min = 30._dbl_kind ! minimum snow grain radius
rsnw_snicar_max = 1500._dbl_kind ! maximum snow grain radius
rsnw_datatype = 'sorted_idelta1' ! sorted "integers" with constant delta 1
allocate(rsnw_snicar_tab(nmbrad_snicar))
rsnw_snicar_tab(1) = real(rsnw_snicar_min,dbl_kind)
rsnw_snicar_tab(1) = rsnw_snicar_min
do n = 1, nmbrad_snicar-1
rsnw_snicar_tab(n+1) = rsnw_snicar_tab(n) + 1.0_dbl_kind
enddo
Expand Down
2 changes: 0 additions & 2 deletions configuration/scripts/options/set_env.snicar
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1 @@
setenv ICE_IOTYPE netcdf
setenv ICE_NXGLOB 1
setenv ICE_SNICARHC true
2 changes: 0 additions & 2 deletions configuration/scripts/options/set_env.snicartest
View file

This file was deleted.

4 changes: 4 additions & 0 deletions configuration/scripts/tests/base_suite.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ smoke col 1x1 debug,run1year,calcdragio
smoke col 1x1 debug,run1year,modal
smoke col 1x1 debug,run1year,fluxopenw
smoke col 1x1 debug,run1year,dyn,fluxopenw
smoke col 1x1 debug,run1year,debug,snicartest
smoke col 1x1 debug,run1year,debug,snicar
restart col 1x1 debug
restart col 1x1 diag1
restart col 1x1 pondlvl
Expand All @@ -39,4 +41,6 @@ restart col 1x1 fsd12,short
restart col 1x1 snwitdrdg,snwgrain
restart col 1x1 modal
restart col 1x1 fluxopenw
restart col 1x1 snicartest
restart col 1x1 snicar

7 changes: 0 additions & 7 deletions configuration/scripts/tests/snicar_suite.ts
View file

This file was deleted.

0 comments on commit 6ad0f44

Please sign in to comment.