Merge CICE-Consortium/Icepack main to E3SM-Project/Icepack main

.github/PULL_REQUEST_TEMPLATE.md

@@ -26,4 +26,4 @@ please refer to: <https://github.com/CICE-Consortium/About-Us/wiki/Resource-Inde
     - [ ] No, does the documentation need to be updated at a later time?
         - [ ] Yes
         - [ ] No 
-- [ ] Please provide any additional information or relevant details below:
+- [ ] Please document the changes in detail, including _why_ the changes are made.  This will become part of the PR commit log.

columnphysics/icepack_itd.F90

@@ -135,7 +135,7 @@ subroutine rebin (ntrcr,    trcr_depend,     &
 
       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
          hin_max   ! category limits (m)
- 
+
       real (kind=dbl_kind), intent(in) :: &
          Tf                ! freezing temperature
 
@@ -785,11 +785,11 @@ subroutine cleanup_itd (dt,          ntrcr,      &
          ntrcr , & ! number of tracers in use
          nbtrcr, & ! number of bio tracers in use
          n_aero    ! number of aerosol tracers
- 
-      real (kind=dbl_kind), intent(in) :: & 
-         dt        ! time step 
- 
-      real (kind=dbl_kind), intent(in) :: & 
+
+      real (kind=dbl_kind), intent(in) :: &
+         dt        ! time step
+
+      real (kind=dbl_kind), intent(in) :: &
          Tf        ! Freezing temperature
 
       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
@@ -1923,7 +1923,7 @@ subroutine icepack_aggregate (ncat,               &
                                    atrcr,     aice,          &
                                    vice ,     vsno,          &
                                    trcr_base, n_trcr_strata, &
-                                   nt_strata, trcr, Tf)   
+                                   nt_strata, trcr, Tf)
       if (icepack_warnings_aborted(subname)) return
 
       deallocate (atrcr)

columnphysics/icepack_mechred.F90

@@ -68,7 +68,6 @@ module icepack_mechred
          Gstar  = p15     , & ! max value of G(h) that participates
                               ! (krdg_partic = 0)
          astar  = p05     , & ! e-folding scale for G(h) participation
-!echmod         astar  = p1        , & ! e-folding scale for G(h) participation
                               ! (krdg_partic = 1)
          maxraft= c1      , & ! max value of hrmin - hi = max thickness
                               ! of ice that rafts (m)
@@ -401,7 +400,7 @@ subroutine ridge_ice (dt,          ndtd,       &
                            msnow_mlt,   esnow_mlt,   &
                            maero,       miso,        &
                            mpond,       Tf,          &
-                           aredistn,    vredistn)    
+                           aredistn,    vredistn)
          if (icepack_warnings_aborted(subname)) return
 
       !-----------------------------------------------------------------

columnphysics/icepack_parameters.F90

@@ -1799,9 +1799,6 @@ subroutine icepack_query_parameters(   &
       if (present(sw_frac_out)           ) sw_frac_out      = sw_frac
       if (present(sw_dtemp_out)          ) sw_dtemp_out     = sw_dtemp
 
-      call icepack_recompute_constants()
-      if (icepack_warnings_aborted(subname)) return
-
       end subroutine icepack_query_parameters
 
 !=======================================================================

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
@@ -1845,7 +1844,7 @@ subroutine compute_dEdd_3bd(                           &
          delr    , & ! snow grain radius interpolation parameter
       ! inherent optical properties (iop) for snow
          Qs      , & ! Snow extinction efficiency
-         ks      , & ! Snow mass extinction coefficient (m^2/kg)  ! echmod - CHECK units and usage
+         ks      , & ! Snow mass extinction coefficient (1/m)
          ws      , & ! Snow single scattering albedo
          gs          ! Snow asymmetry parameter
 
@@ -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))
@@ -4513,7 +4508,7 @@ subroutine compute_dEdd_5bd(                          &
       ! copy/point to table data for local names
       asm_prm_ice_drc => ssp_sasymmdr
       asm_prm_ice_dfs => ssp_sasymmdf
-      ss_alb_ice_drc => ssp_snwalbdr 
+      ss_alb_ice_drc => ssp_snwalbdr
       ss_alb_ice_dfs => ssp_snwalbdf
       ext_cff_mss_ice_drc => ssp_snwextdr
       ext_cff_mss_ice_dfs => ssp_snwextdf
@@ -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