Update for new snow scheme + introduce stop_now_cpl for coupling (#641)

cicecore/drivers/nuopc/dmi/CICE_InitMod.F90

@@ -18,6 +18,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, &
@@ -81,15 +82,14 @@ subroutine cice_init(mpi_comm)
       use ice_flux, only: init_coupler_flux, init_history_therm, &
           init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
-          get_forcing_atmo, get_forcing_ocn, get_wave_spec
+          get_forcing_atmo, get_forcing_ocn, get_wave_spec, init_snowtable
       use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, &
           faero_default, faero_optics, alloc_forcing_bgc, fiso_default
       use ice_grid, only: init_grid1, init_grid2, alloc_grid
       use ice_history, only: init_hist, accum_hist
       use ice_restart_shared, only: restart, runtype
       use ice_init, only: input_data, init_state
-      use ice_init_column, only: init_thermo_vertical, init_shortwave, &
-          init_zbgc, input_zbgc, count_tracers
+      use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers
       use ice_kinds_mod
       use ice_restoring, only: ice_HaloRestore_init
       use ice_timers, only: timer_total, init_ice_timers, ice_timer_start
@@ -99,7 +99,8 @@ subroutine cice_init(mpi_comm)
          mpi_comm ! communicator for sequential ccsm
 
       logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, &
-         tr_iso, tr_fsd, wave_spec
+         tr_iso, tr_fsd, wave_spec, tr_snow
+      character(len=char_len) :: snw_aging_table
       character(len=*), parameter :: subname = '(cice_init)'
 
       if (present(mpi_comm)) then
@@ -176,7 +177,7 @@ subroutine cice_init(mpi_comm)
       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__)
@@ -190,7 +191,7 @@ subroutine cice_init(mpi_comm)
       call calc_timesteps       ! update timestep counter if not using npt_unit="1"
 
       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__)
@@ -222,8 +223,20 @@ subroutine cice_init(mpi_comm)
       call get_forcing_ocn(dt)  ! ocean forcing from data
 #endif
 
+      ! 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
+
       ! isotopes
       if (tr_iso)     call fiso_default                 ! default values
+
       ! aerosols
       ! if (tr_aero)  call faero_data                   ! data file
       ! if (tr_zaero) call fzaero_data                  ! data file (gx1)
@@ -252,19 +265,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, &
@@ -279,12 +293,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_snow, tr_fsd, tr_iso, tr_aero, tr_brine, &
           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)'
@@ -299,10 +314,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, &
@@ -399,6 +416,22 @@ 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.
@@ -415,7 +448,7 @@ subroutine init_restart
          if (restart_iso) then
             call read_restart_iso
          else
-            do iblk = 1, nblocks 
+            do iblk = 1, nblocks
                call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), &
                                  trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk))
             enddo ! iblk

cicecore/drivers/nuopc/dmi/CICE_RunMod.F90

@@ -41,7 +41,7 @@ module CICE_RunMod
 !         Philip W. Jones, LANL
 !         William H. Lipscomb, LANL
 
-      subroutine CICE_Run
+      subroutine CICE_Run(stop_now_cpl)
 
       use ice_calendar, only: istep, istep1, dt, stop_now, advance_timestep
       use ice_forcing, only: get_forcing_atmo, get_forcing_ocn, &
@@ -54,6 +54,7 @@ subroutine CICE_Run
       logical (kind=log_kind) :: &
           tr_iso, tr_aero, tr_zaero, skl_bgc, z_tracers, wave_spec, tr_fsd
       character(len=*), parameter :: subname = '(CICE_Run)'
+      logical (kind=log_kind), optional, intent(in) :: stop_now_cpl
 
    !--------------------------------------------------------------------
    !  initialize error code and step timer
@@ -84,6 +85,9 @@ subroutine CICE_Run
 
          call advance_timestep()     ! advance time
 
+         if (present(stop_now_cpl)) then
+           if (stop_now_cpl) return
+         endif
 #ifndef CICE_IN_NEMO
 #ifndef CICE_DMI
          if (stop_now >= 1) exit timeLoop
@@ -142,7 +146,7 @@ subroutine ice_step
 
       use ice_boundary, only: ice_HaloUpdate
       use ice_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep
-      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
@@ -155,12 +159,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
@@ -174,19 +179,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__)
@@ -223,14 +237,36 @@ 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
+
                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
 
          enddo ! iblk
@@ -256,38 +292,80 @@ 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
 
-      !-----------------------------------------------------------------
-      ! albedo, shortwave radiation
-      !-----------------------------------------------------------------
+         if (debug_model) then
+            plabeld = 'post dynamics'
+            do iblk = 1, nblocks
+               call debug_ice (iblk, plabeld)
+            enddo
+         endif
 
          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
 
+      !-----------------------------------------------------------------
+      ! albedo, shortwave radiation
+      !-----------------------------------------------------------------
+
             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
 
@@ -325,6 +403,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