recreate cap update branch (CICE-Consortium#42)

* add debug_model feature
* add required variables and calls for tr_snow

cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90

@@ -8,6 +8,7 @@ module CICE_InitMod
   use icepack_intfc, only: icepack_aggregate
   use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
   use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave
+  use icepack_intfc, only: icepack_init_snow
   use icepack_intfc, only: icepack_configure
   use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
   use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags
@@ -83,7 +84,7 @@ subroutine cice_init2()
     use ice_dyn_vp           , only: init_vp
     use ice_flux             , only: init_coupler_flux, init_history_therm
     use ice_flux             , only: init_history_dyn, init_flux_atm, init_flux_ocn
-    use ice_forcing          , only: init_forcing_ocn
+    use ice_forcing          , only: init_forcing_ocn, init_snowtable
     use ice_forcing_bgc      , only: get_forcing_bgc, get_atm_bgc
     use ice_forcing_bgc      , only: faero_default, faero_optics, alloc_forcing_bgc, fiso_default
     use ice_history          , only: init_hist, accum_hist
@@ -95,7 +96,8 @@ subroutine cice_init2()
     use ice_transport_driver , only: init_transport
 
     logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers
-    logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec
+    logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow
+    character(len=char_len) :: snw_aging_table
     character(len=*), parameter :: subname = '(cice_init2)'
     !----------------------------------------------------
 
@@ -145,7 +147,7 @@ subroutine cice_init2()
     call ice_HaloRestore_init ! restored boundary conditions
 
     call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, &
-         wave_spec_out=wave_spec)
+         wave_spec_out=wave_spec, snw_aging_table_out=snw_aging_table)
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(trim(subname), &
          file=__FILE__,line= __LINE__)
@@ -158,7 +160,7 @@ subroutine cice_init2()
     call init_history_dyn     ! initialize dynamic history variables
 
     call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero)
-    call icepack_query_tracer_flags(tr_iso_out=tr_iso)
+    call icepack_query_tracer_flags(tr_iso_out=tr_iso, tr_snow_out=tr_snow)
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(trim(subname), &
          file=__FILE__,line= __LINE__)
@@ -167,6 +169,17 @@ subroutine cice_init2()
        call faero_optics !initialize aerosol optical property tables
     end if
 
+    ! snow aging lookup table initialization
+    if (tr_snow) then         ! advanced snow physics
+       call icepack_init_snow()
+       call icepack_warnings_flush(nu_diag)
+       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+          file=__FILE__, line=__LINE__)
+       if (snw_aging_table(1:4) /= 'test') then
+          call init_snowtable()
+       endif
+    endif
+
     ! Initialize shortwave components using swdn from previous timestep
     ! if restarting. These components will be scaled to current forcing
     ! in prep_radiation.
@@ -199,19 +212,20 @@ subroutine init_restart()
     use ice_calendar, only: calendar
     use ice_constants, only: c0
     use ice_domain, only: nblocks
-    use ice_domain_size, only: ncat, n_iso, n_aero, nfsd
+    use ice_domain_size, only: ncat, n_iso, n_aero, nfsd, nslyr
     use ice_dyn_eap, only: read_restart_eap
     use ice_dyn_shared, only: kdyn
     use ice_grid, only: tmask
     use ice_init, only: ice_ic
-    use ice_init_column, only: init_age, init_FY, init_lvl, &
+    use ice_init_column, only: init_age, init_FY, init_lvl, init_snowtracers, &
          init_meltponds_cesm,  init_meltponds_lvl, init_meltponds_topo, &
          init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd
     use ice_restart_column, only: restart_age, read_restart_age, &
          restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, &
          restart_pond_cesm, read_restart_pond_cesm, &
          restart_pond_lvl, read_restart_pond_lvl, &
          restart_pond_topo, read_restart_pond_topo, &
+         restart_snow, read_restart_snow, &
          restart_fsd, read_restart_fsd, &
          restart_iso, read_restart_iso, &
          restart_aero, read_restart_aero, &
@@ -226,12 +240,13 @@ subroutine init_restart()
          iblk            ! block index
     logical(kind=log_kind) :: &
          tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, &
-         tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, &
+         tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, tr_snow, &
          skl_bgc, z_tracers, solve_zsal
     integer(kind=int_kind) :: &
          ntrcr
     integer(kind=int_kind) :: &
          nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, &
+         nt_smice, nt_smliq, nt_rhos, nt_rsnw, &
          nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice
 
     character(len=*), parameter :: subname = '(init_restart)'
@@ -247,10 +262,12 @@ subroutine init_restart()
     call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
          tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, &
          tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, &
-         tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
+         tr_snow_out=tr_snow, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
     call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, &
          nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, &
          nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, &
+         nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, &
+         nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw, &
          nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
@@ -347,6 +364,21 @@ subroutine init_restart()
           enddo ! iblk
        endif ! .not. restart_pond
     endif
+    ! snow redistribution/metamorphism
+    if (tr_snow) then
+       if (trim(runtype) == 'continue') restart_snow = .true.
+       if (restart_snow) then
+          call read_restart_snow
+       else
+          do iblk = 1, nblocks
+             call init_snowtracers(trcrn(:,:,nt_smice:nt_smice+nslyr-1,:,iblk), &
+                                   trcrn(:,:,nt_smliq:nt_smliq+nslyr-1,:,iblk), &
+                                   trcrn(:,:,nt_rhos :nt_rhos +nslyr-1,:,iblk), &
+                                   trcrn(:,:,nt_rsnw :nt_rsnw +nslyr-1,:,iblk))
+          enddo ! iblk
+       endif
+    endif
+
     ! floe size distribution
     if (tr_fsd) then
        if (trim(runtype) == 'continue') restart_fsd = .true.
@@ -441,7 +473,6 @@ subroutine init_restart()
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
-
   end subroutine init_restart
 
   !=======================================================================

cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90

@@ -110,7 +110,7 @@ subroutine ice_step
       use ice_boundary, only: ice_HaloUpdate
       use ice_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep
       use ice_calendar, only: idate, msec
-      use ice_diagnostics, only: init_mass_diags, runtime_diags
+      use ice_diagnostics, only: init_mass_diags, runtime_diags, debug_model, debug_ice
       use ice_diagnostics_bgc, only: hbrine_diags, zsal_diags, bgc_diags
       use ice_domain, only: halo_info, nblocks
       use ice_dyn_eap, only: write_restart_eap
@@ -123,12 +123,13 @@ subroutine ice_step
       use ice_restart_column, only: write_restart_age, write_restart_FY, &
           write_restart_lvl, write_restart_pond_cesm, write_restart_pond_lvl, &
           write_restart_pond_topo, write_restart_aero, write_restart_fsd, &
-          write_restart_iso, write_restart_bgc, write_restart_hbrine
+          write_restart_iso, write_restart_bgc, write_restart_hbrine, &
+          write_restart_snow
       use ice_restart_driver, only: dumpfile
       use ice_restoring, only: restore_ice, ice_HaloRestore
       use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, &
           update_state, step_dyn_horiz, step_dyn_ridge, step_radiation, &
-          biogeochemistry, save_init, step_dyn_wave
+          biogeochemistry, save_init, step_dyn_wave, step_snow
       use ice_timers, only: ice_timer_start, ice_timer_stop, &
           timer_diags, timer_column, timer_thermo, timer_bound, &
           timer_hist, timer_readwrite
@@ -144,19 +145,28 @@ subroutine ice_step
          offset          ! d(age)/dt time offset
 
       logical (kind=log_kind) :: &
-          tr_iage, tr_FY, tr_lvl, tr_fsd, &
+          tr_iage, tr_FY, tr_lvl, tr_fsd, tr_snow, &
           tr_pond_cesm, tr_pond_lvl, tr_pond_topo, tr_brine, tr_iso, tr_aero, &
           calc_Tsfc, skl_bgc, solve_zsal, z_tracers, wave_spec
 
       character(len=*), parameter :: subname = '(ice_step)'
 
+      character (len=char_len) :: plabeld
+
+      if (debug_model) then
+         plabeld = 'beginning time step'
+         do iblk = 1, nblocks
+            call debug_ice (iblk, plabeld)
+         enddo
+      endif
+
       call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc, skl_bgc_out=skl_bgc, &
            solve_zsal_out=solve_zsal, z_tracers_out=z_tracers, ktherm_out=ktherm, &
            wave_spec_out=wave_spec)
       call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
            tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, &
            tr_pond_topo_out=tr_pond_topo, tr_brine_out=tr_brine, tr_aero_out=tr_aero, &
-           tr_iso_out=tr_iso, tr_fsd_out=tr_fsd)
+           tr_iso_out=tr_iso, tr_fsd_out=tr_fsd, tr_snow_out=tr_snow)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
@@ -201,15 +211,33 @@ subroutine ice_step
       !-----------------------------------------------------------------
 
                if (calc_Tsfc) call prep_radiation (iblk)
+               if (debug_model) then
+                  plabeld = 'post prep_radiation'
+                  call debug_ice (iblk, plabeld)
+               endif
 
       !-----------------------------------------------------------------
       ! thermodynamics and biogeochemistry
       !-----------------------------------------------------------------
 
                call step_therm1     (dt, iblk) ! vertical thermodynamics
+               if (debug_model) then
+                  plabeld = 'post step_therm1'
+                  call debug_ice (iblk, plabeld)
+               endif
+
                call biogeochemistry (dt, iblk) ! biogeochemistry
+               if (debug_model) then
+                  plabeld = 'post biogeochemistry'
+                  call debug_ice (iblk, plabeld)
+               endif
+
                if (.not.prescribed_ice) &
                   call step_therm2  (dt, iblk) ! ice thickness distribution thermo
+                  if (debug_model) then
+                     plabeld = 'post step_therm2'
+                     call debug_ice (iblk, plabeld)
+                  endif
 
             endif ! ktherm > 0
 
@@ -237,19 +265,37 @@ subroutine ice_step
 
             ! momentum, stress, transport
             call step_dyn_horiz (dt_dyn)
+            if (debug_model) then
+               plabeld = 'post step_dyn_horiz'
+               do iblk = 1, nblocks
+                  call debug_ice (iblk, plabeld)
+               enddo ! iblk
+            endif
 
             ! ridging
             !$OMP PARALLEL DO PRIVATE(iblk)
             do iblk = 1, nblocks
                if (kridge > 0) call step_dyn_ridge (dt_dyn, ndtd, iblk)
             enddo
             !$OMP END PARALLEL DO
+            if (debug_model) then
+               plabeld = 'post step_dyn_ridge'
+               do iblk = 1, nblocks
+                  call debug_ice (iblk, plabeld)
+               enddo ! iblk
+            endif
 
             ! clean up, update tendency diagnostics
             offset = c0
             call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
 
          enddo
+         if (debug_model) then
+            plabeld = 'post dynamics'
+            do iblk = 1, nblocks
+               call debug_ice (iblk, plabeld)
+            enddo
+         endif
 
          endif  ! not prescribed ice
 
@@ -260,18 +306,36 @@ subroutine ice_step
          call ice_timer_start(timer_column)  ! column physics
          call ice_timer_start(timer_thermo)  ! thermodynamics
 
+      !-----------------------------------------------------------------
+      ! snow redistribution and metamorphosis
+      !-----------------------------------------------------------------
+
+         if (tr_snow) then         ! advanced snow physics
+            do iblk = 1, nblocks
+               call step_snow (dt, iblk)
+            enddo
+            call update_state (dt) ! clean up
+         endif
+
 !MHRI: CHECK THIS OMP
          !$OMP PARALLEL DO PRIVATE(iblk)
          do iblk = 1, nblocks
 
             if (ktherm >= 0) call step_radiation (dt, iblk)
+            if (debug_model) then
+               plabeld = 'post step_radiation'
+               call debug_ice (iblk, plabeld)
+            endif
 
       !-----------------------------------------------------------------
       ! get ready for coupling and the next time step
       !-----------------------------------------------------------------
 
             call coupling_prep (iblk)
-
+            if (debug_model) then
+               plabeld = 'post coupling_prep'
+               call debug_ice (iblk, plabeld)
+            endif
          enddo ! iblk
          !$OMP END PARALLEL DO
 
@@ -309,6 +373,7 @@ subroutine ice_step
             if (tr_pond_cesm) call write_restart_pond_cesm
             if (tr_pond_lvl)  call write_restart_pond_lvl
             if (tr_pond_topo) call write_restart_pond_topo
+            if (tr_snow)      call write_restart_snow
             if (tr_fsd)       call write_restart_fsd
             if (tr_iso)       call write_restart_iso
             if (tr_aero)      call write_restart_aero
@@ -634,11 +699,12 @@ subroutine sfcflux_to_ocn(nx_block,   ny_block,     &
 
       real (kind=dbl_kind)    :: &
           puny, &          !
+          Lsub, &          !
           rLsub            ! 1/Lsub
 
       character(len=*), parameter :: subname = '(sfcflux_to_ocn)'
 
-      call icepack_query_parameters(puny_out=puny)
+      call icepack_query_parameters(puny_out=puny, Lsub_out=Lsub)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)