Update Icepack interfaces consistent with recent changes in Icepack (#…

…379)

* update icepack interfaces consistent with icepack 31734de3b376

* update icepack

* remove n_aero from icepack therm1 and therm2 arguments
migrate use of non icepack_intfc to icepack_intfc
rename some icepack mushy routines to meet naming convention
works with icepack #96970594235d81b

* update interface call

* update icepack and update CICE version number to 6.1

cicecore/cicedynB/analysis/ice_history.F90

@@ -1742,8 +1742,9 @@ subroutine accum_hist (dt)
       use ice_history_mechred, only: accum_hist_mechred
       use ice_history_pond, only: accum_hist_pond
       use ice_history_drag, only: accum_hist_drag
+      use icepack_intfc, only: icepack_mushy_density_brine, icepack_mushy_liquid_fraction
+      use icepack_intfc, only: icepack_mushy_temperature_mush
       use ice_history_fsd, only: accum_hist_fsd
-      use icepack_mushy_physics, only: density_brine, liquid_fraction, temperature_mush
       use ice_state ! almost everything
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite
 
@@ -2697,13 +2698,13 @@ subroutine accum_hist (dt)
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
                  if (ktherm == 2) then
-                    rho_ocn = density_brine(sss(i,j,iblk))
+                    rho_ocn = icepack_mushy_density_brine(sss(i,j,iblk))
                     rho_ice = c0
                     do k = 1, nzilyr
-                       Tice = temperature_mush(trcr(i,j,nt_qice+k-1,iblk),trcr(i,j,nt_sice+k-1,iblk))
+                       Tice = icepack_mushy_temperature_mush(trcr(i,j,nt_qice+k-1,iblk),trcr(i,j,nt_sice+k-1,iblk))
                        Sbr = trcr(i,j,nt_sice+k-1,iblk)
-                       phi = liquid_fraction(Tice,Sbr)
-                       rhob = density_brine(Sbr)
+                       phi = icepack_mushy_liquid_fraction(Tice,Sbr)
+                       rhob = icepack_mushy_density_brine(Sbr)
                        rho_ice = rho_ice + min(phi*rhob+(c1-phi)*rhoi,rho_ocn)
                     enddo
                     rho_ice = rho_ice / real(nzilyr,kind=dbl_kind)

cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90

@@ -36,6 +36,8 @@ module ice_dyn_evp_1d
   use ice_kinds_mod
   use ice_fileunits, only: nu_diag
   use ice_exit, only: abort_ice
+  use icepack_intfc, only: icepack_query_parameters
+  use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
 
   implicit none
   private
@@ -244,7 +246,6 @@ subroutine stress_i(NA_len, &
 
     use ice_kinds_mod
     use ice_constants, only: p027, p055, p111, p166, p222, p25, p333, p5, c1p5, c1
-    use icepack_parameters, only: puny
     use ice_dyn_shared, only: ecci, denom1, arlx1i, Ktens, revp
 
     implicit none
@@ -264,6 +265,8 @@ subroutine stress_i(NA_len, &
     !-- local variables
 
     integer (kind=int_kind) :: iw,il,iu
+    real    (kind=dbl_kind) ::                                                   &
+      puny
     real    (kind=DBL_KIND) ::                                                   &
       divune, divunw, divuse, divusw,tensionne, tensionnw, tensionse, tensionsw, & 
       shearne, shearnw, shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw   , &
@@ -279,6 +282,11 @@ subroutine stress_i(NA_len, &
     character(len=*), parameter :: subname = '(stress_i)'
     !---------------------------------------
 
+    call icepack_query_parameters(puny_out=puny)
+    call icepack_warnings_flush(nu_diag)
+    if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+        file=__FILE__, line=__LINE__)
+
 #ifdef _OPENACC
     !$acc parallel                                                                 &
     !$acc present(ee,ne,se,strength,uvel,vvel,dxt,dyt,                             &
@@ -509,7 +517,6 @@ subroutine stress_l(NA_len, tarear, &
 
     use ice_kinds_mod
     use ice_constants, only: p027, p055, p111, p166, p222, p25, p333, p5, c1p5, c0, c1
-    use icepack_parameters, only: puny
     use ice_dyn_shared, only: ecci, denom1, arlx1i, Ktens, revp
 
     implicit none
@@ -531,6 +538,8 @@ subroutine stress_l(NA_len, tarear, &
     !-- local variables
 
     integer (kind=int_kind) :: iw,il,iu
+    real    (kind=dbl_kind) ::                                                   &
+      puny
     real    (kind=DBL_KIND) ::                                                   &
       divune, divunw, divuse, divusw,tensionne, tensionnw, tensionse, tensionsw, & 
       shearne, shearnw, shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw   , &
@@ -546,6 +555,11 @@ subroutine stress_l(NA_len, tarear, &
     character(len=*), parameter :: subname = '(stress_l)'
     !---------------------------------------
 
+    call icepack_query_parameters(puny_out=puny)
+    call icepack_warnings_flush(nu_diag)
+    if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+        file=__FILE__, line=__LINE__)
+
 #ifdef _OPENACC
     !$acc parallel                                                                 &
     !$acc present(ee,ne,se,strength,uvel,vvel,dxt,dyt,tarear,                      &
@@ -867,7 +881,6 @@ subroutine stepu_last(NA_len, rhow, &
 
     use ice_kinds_mod
     use ice_constants, only: c0, c1
-    use icepack_parameters, only: puny
     use ice_dyn_shared, only: brlx, revp, basalstress
 
     implicit none
@@ -1397,7 +1410,6 @@ subroutine evp_kernel_v2
 
     use ice_constants, only : c0
     use ice_dyn_shared, only: ndte
-    use icepack_intfc, only: icepack_query_parameters
     use ice_communicate, only: my_task, master_task
     implicit none
 
@@ -1413,6 +1425,9 @@ subroutine evp_kernel_v2
 
       !- Read constants...
       call icepack_query_parameters(rhow_out=rhow)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+          file=__FILE__, line=__LINE__)
       na=NA_len
       nb=NA_len
       navel=NAVEL_len

cicecore/cicedynB/general/ice_init.F90

@@ -27,6 +27,7 @@ module ice_init
       use icepack_intfc, only: icepack_init_trcr
       use icepack_intfc, only: icepack_init_parameters
       use icepack_intfc, only: icepack_init_tracer_flags
+      use icepack_intfc, only: icepack_init_tracer_indices
       use icepack_intfc, only: icepack_query_tracer_flags
       use icepack_intfc, only: icepack_query_tracer_numbers
       use icepack_intfc, only: icepack_query_tracer_indices
@@ -1260,6 +1261,10 @@ subroutine input_data
       call icepack_init_tracer_flags(tr_iage_in=tr_iage, tr_FY_in=tr_FY, &
          tr_lvl_in=tr_lvl, tr_aero_in=tr_aero, tr_fsd_in=tr_fsd, tr_pond_in=tr_pond, &
          tr_pond_cesm_in=tr_pond_cesm, tr_pond_lvl_in=tr_pond_lvl, tr_pond_topo_in=tr_pond_topo)
+      call icepack_init_tracer_indices(ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
 
  1000    format (a30,2x,f9.2)  ! a30 to align formatted, unformatted statements
  1005    format (a30,2x,f12.6) ! float
@@ -1968,7 +1973,6 @@ subroutine boxslotcyl_data_vel(i,        j,       &
       use ice_constants, only: c2, c12, p5, cm_to_m
       use ice_domain_size, only: nx_global, ny_global
       use ice_grid, only: dxrect
-      use icepack_parameters, only: secday, pi
 
       integer (kind=int_kind), intent(in) :: &
          i, j,               & ! local indices
@@ -1980,14 +1984,21 @@ subroutine boxslotcyl_data_vel(i,        j,       &
          uvel, vvel            ! ice velocity
 
       ! local variables
-         
+
       real (kind=dbl_kind) :: &
+         pi             , & ! pi
+         secday         , & ! seconds per day
          max_vel        , & ! max velocity
          domain_length  , & ! physical domain length
          period             ! rotational period
 
       character(len=*), parameter :: subname = '(boxslotcyl_data_vel)'
 
+      call icepack_query_parameters(secday_out=secday, pi_out=pi)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
       domain_length = dxrect*cm_to_m*nx_global
       period        = c12*secday               ! 12 days rotational period
       max_vel       = pi*domain_length/period

cicecore/cicedynB/general/ice_step_mod.F90

@@ -24,7 +24,7 @@ module ice_step_mod
       use icepack_intfc, only: icepack_step_wavefracture
       use icepack_intfc, only: icepack_step_radiation
       use icepack_intfc, only: icepack_ocn_mixed_layer, icepack_atm_boundary
-      use icepack_intfc, only: icepack_biogeochemistry, icepack_init_OceanConcArray
+      use icepack_intfc, only: icepack_biogeochemistry, icepack_load_ocean_bio_array
       use icepack_intfc, only: icepack_max_algae, icepack_max_nbtrcr, icepack_max_don
       use icepack_intfc, only: icepack_max_doc, icepack_max_dic, icepack_max_aero
       use icepack_intfc, only: icepack_max_fe
@@ -286,7 +286,7 @@ subroutine step_therm1 (dt, iblk)
 
          if (tmask(i,j,iblk)) &
          call icepack_step_therm1(dt=dt, ncat=ncat,            &
-                      nilyr=nilyr, nslyr=nslyr, n_aero=n_aero, &
+                      nilyr=nilyr, nslyr=nslyr,                &
                       aicen_init   = aicen_init  (i,j,:,iblk), &
                       vicen_init   = vicen_init  (i,j,:,iblk), &
                       vsnon_init   = vsnon_init  (i,j,:,iblk), &
@@ -511,7 +511,7 @@ subroutine step_therm2 (dt, iblk)
          if (tr_fsd) &
          wave_sig_ht(i,j,iblk) = c4*SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:)))
 
-         call icepack_step_therm2(dt=dt, ncat=ncat, n_aero=n_aero, &
+         call icepack_step_therm2(dt=dt, ncat=ncat, &
                       nltrcr=nltrcr, nilyr=nilyr, nslyr=nslyr, nblyr=nblyr, &
                       hin_max    = hin_max   (:),          &   
                       aicen      = aicen     (i,j,:,iblk), &
@@ -974,6 +974,9 @@ subroutine step_radiation (dt, iblk)
          nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nlt_chl_sw, &
          ntrcr, nbtrcr, nbtrcr_sw, nt_fbri
 
+      integer (kind=int_kind), dimension(icepack_max_algae) :: &
+         nt_bgc_N
+
       integer (kind=int_kind), dimension(icepack_max_aero) :: &
          nlt_zaero_sw, nt_zaero
 
@@ -984,7 +987,6 @@ subroutine step_radiation (dt, iblk)
          fbri                 ! brine height to ice thickness
 
       real(kind= dbl_kind), dimension(:,:), allocatable :: &
-         ztrcr    , &
          ztrcr_sw
 
       logical (kind=log_kind) :: &
@@ -1003,14 +1005,13 @@ subroutine step_radiation (dt, iblk)
          nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, &
          nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, &
          nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, &
-         nt_fbri_out=nt_fbri, nt_zaero_out=nt_zaero)
+         nt_fbri_out=nt_fbri, nt_zaero_out=nt_zaero, nt_bgc_N_out=nt_bgc_N)
       call icepack_query_parameters(dEdd_algae_out=dEdd_algae, modal_aero_out=modal_aero)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
-      allocate(ztrcr(ntrcr,ncat))
-      allocate(ztrcr_sw(ntrcr,ncat))
+      allocate(ztrcr_sw(nbtrcr_sw,ncat))
 
       this_block = get_block(blocks_ice(iblk),iblk)         
       ilo = this_block%ilo
@@ -1035,20 +1036,13 @@ subroutine step_radiation (dt, iblk)
          fbri(:) = c0
          ztrcr_sw(:,:) = c0
          do n = 1, ncat
-           do k = 1, ntrcr
-             ztrcr(k,n) = trcrn(i,j,k,n,iblk)
-           enddo
            if (tr_brine)  fbri(n) = trcrn(i,j,nt_fbri,n,iblk)
          enddo
 
          if (tmask(i,j,iblk)) then
 
             call icepack_step_radiation (dt=dt,   ncat=ncat,                  &
-                         n_algae=n_algae, tr_zaero=tr_zaero, nblyr=nblyr,     &
-                         ntrcr=ntrcr, nbtrcr_sw=nbtrcr_sw,                    &
-                         nilyr=nilyr, nslyr=nslyr, n_aero=n_aero,             &
-                         n_zaero=n_zaero, nlt_chl_sw=nlt_chl_sw,              &
-                         nlt_zaero_sw=nlt_zaero_sw(:),                        &
+                         nblyr=nblyr, nilyr=nilyr, nslyr=nslyr,               &
                          dEdd_algae=dEdd_algae,                               &
                          swgrid=swgrid(:),        igrid=igrid(:),             &
                          fbri=fbri(:),                                        &
@@ -1061,8 +1055,9 @@ subroutine step_radiation (dt, iblk)
                          hpndn=trcrn(i,j,nt_hpnd,:,iblk),                     &
                          ipndn=trcrn(i,j,nt_ipnd,:,iblk),                     &
                          aeron=trcrn(i,j,nt_aero:nt_aero+4*n_aero-1,:,iblk),  &
-                         zbion=ztrcr_sw,                                      &
-                         trcrn=ztrcr,                                         &
+                         bgcNn=trcrn(i,j,nt_bgc_N(1):nt_bgc_N(1)+n_algae*(nblyr+3)-1,:,iblk), &
+                         zaeron=trcrn(i,j,nt_zaero(1):nt_zaero(1)+n_zaero*(nblyr+3)-1,:,iblk), &
+                         trcrn_bgcsw=ztrcr_sw,                                &
                          TLAT=TLAT(i,j,iblk),     TLON=TLON(i,j,iblk),        &
                          calendar_type=calendar_type,                         &
                          days_per_year=days_per_year,                         &
@@ -1104,7 +1099,6 @@ subroutine step_radiation (dt, iblk)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
-      deallocate(ztrcr)
       deallocate(ztrcr_sw)
 
       call ice_timer_stop(timer_sw)     ! shortwave
@@ -1350,7 +1344,7 @@ subroutine biogeochemistry (dt, iblk)
       do j = jlo, jhi
       do i = ilo, ihi    
 
-         call icepack_init_OceanConcArray(max_nbtrcr = icepack_max_nbtrcr, &
+         call icepack_load_ocean_bio_array(max_nbtrcr = icepack_max_nbtrcr, &
                 max_algae = icepack_max_algae, max_don = icepack_max_don, &
                 max_doc   = icepack_max_doc,   max_dic = icepack_max_dic, &
                 max_aero  = icepack_max_aero,  max_fe  = icepack_max_fe,  &

cicecore/shared/ice_arrays_column.F90

@@ -298,17 +298,15 @@ module ice_arrays_column
 !=======================================================================
 
       subroutine alloc_arrays_column
-
-      ! Allocate column arrays
-
-      use ice_exit, only: abort_ice
-      integer (int_kind) :: max_nbtrcr, max_algae, max_aero, &
-         nmodal1, nmodal2, max_don
-      integer (int_kind) :: ierr, ntrcr
+        ! Allocate column arrays
+        use ice_exit, only: abort_ice
+        integer (int_kind) :: max_nbtrcr, max_algae, max_aero, &
+           nmodal1, nmodal2, max_don, nbtrcr_sw
+        integer (int_kind) :: ierr, ntrcr
 
       character(len=*),parameter :: subname='(alloc_arrays_column)'
 
-      call icepack_query_tracer_numbers(ntrcr_out=ntrcr)
+      call icepack_query_tracer_numbers(ntrcr_out=ntrcr, nbtrcr_sw_out=nbtrcr_sw)
       call icepack_query_tracer_sizes( max_nbtrcr_out=max_nbtrcr, &
          max_algae_out=max_algae, max_aero_out=max_aero, &
          nmodal1_out=nmodal1, nmodal2_out=nmodal2, max_don_out=max_don)
@@ -390,7 +388,7 @@ subroutine alloc_arrays_column
          ocean_bio_all(nx_block,ny_block,max_nbtrcr,max_blocks), & ! fixed order, all values even for tracers false
          ice_bio_net  (nx_block,ny_block,max_nbtrcr,max_blocks), & ! depth integrated tracer (mmol/m^2) 
          snow_bio_net (nx_block,ny_block,max_nbtrcr,max_blocks), & ! depth integrated snow tracer (mmol/m^2)
-         trcrn_sw     (nx_block,ny_block,ntrcr,ncat,max_blocks), & ! bgc tracers active in the delta-Eddington shortwave 
+         trcrn_sw     (nx_block,ny_block,nbtrcr_sw,ncat,max_blocks), & ! bgc tracers active in the delta-Eddington shortwave 
          algal_peak   (nx_block,ny_block,max_algae,max_blocks), & ! vertical location of algal maximum, 0 if no maximum 
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of Memory2')