diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index f6c1a7223..dfa2d6422 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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)) @@ -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 @@ -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) @@ -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) & @@ -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) @@ -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) & @@ -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 diff --git a/columnphysics/icepack_shortwave_data.F90 b/columnphysics/icepack_shortwave_data.F90 index 50bdf5eaa..b9ebe614d 100644 --- a/columnphysics/icepack_shortwave_data.F90 +++ b/columnphysics/icepack_shortwave_data.F90 @@ -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) @@ -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 @@ -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, & @@ -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 @@ -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 diff --git a/configuration/scripts/options/set_env.snicar b/configuration/scripts/options/set_env.snicar index 13cd94861..91c70cb4b 100644 --- a/configuration/scripts/options/set_env.snicar +++ b/configuration/scripts/options/set_env.snicar @@ -1,3 +1 @@ -setenv ICE_IOTYPE netcdf -setenv ICE_NXGLOB 1 setenv ICE_SNICARHC true diff --git a/configuration/scripts/options/set_env.snicartest b/configuration/scripts/options/set_env.snicartest deleted file mode 100644 index 7ac74351d..000000000 --- a/configuration/scripts/options/set_env.snicartest +++ /dev/null @@ -1,2 +0,0 @@ -setenv ICE_IOTYPE netcdf -setenv ICE_NXGLOB 1 diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index ce98270fd..43e405dab 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -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 @@ -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 diff --git a/configuration/scripts/tests/snicar_suite.ts b/configuration/scripts/tests/snicar_suite.ts deleted file mode 100644 index 2eb4c9b32..000000000 --- a/configuration/scripts/tests/snicar_suite.ts +++ /dev/null @@ -1,7 +0,0 @@ -# Test Grid PEs Sets BFB-compare -smoke col 1x1 diag1,run1year -smoke col 1x1 debug,run1year -smoke col 1x1 diag1,run1year,snicartest -smoke col 1x1 debug,run1year,snicartest -smoke col 1x1 diag1,run1year,snicar -smoke col 1x1 debug,run1year,snicar