Add a new diagnostic variable dvsdtd. (CICE-Consortium#940)

* Add dvsdtd diagnostic for SIMIP sndmassdyn

* Fix the units of SIMIP snow variables

* Fix calls to update_state

* Add present for dagedt and fix indents

* This PR has been updated to include dvsdtt (thermodynamic tendency
for snow). Also, the calls to update_state have been updated in all
of the drivers. Note that I can only test the NUOPC/CMEPS and standalone
drivers, so it would be good if others could test.

* Update the opticep unit test to be consistent with the latest changes.

---------

Co-authored-by: apcraig <[email protected]>

cicecore/cicedyn/analysis/ice_history.F90

@@ -362,6 +362,7 @@ subroutine init_hist (dt)
          f_sidmasslat = 'mxxxx'
          f_sndmasssnf = 'mxxxx'
          f_sndmassmelt = 'mxxxx'
+         f_sndmassdyn = 'mxxxx'
          f_siflswdtop = 'mxxxx'
          f_siflswutop = 'mxxxx'
          f_siflswdbot = 'mxxxx'
@@ -402,6 +403,11 @@ subroutine init_hist (dt)
          f_siu = f_CMIP
          f_siv = f_CMIP
          f_sispeed = f_CMIP
+         f_sndmasssubl = f_CMIP
+         f_sndmasssnf = f_CMIP
+         f_sndmassmelt = f_CMIP
+         f_sndmassdyn = f_CMIP
+         f_sidmasssi = f_CMIP
       endif
 
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
@@ -654,6 +660,7 @@ subroutine init_hist (dt)
       call broadcast_scalar (f_sidmasslat, master_task)
       call broadcast_scalar (f_sndmasssnf, master_task)
       call broadcast_scalar (f_sndmassmelt, master_task)
+      call broadcast_scalar (f_sndmassdyn, master_task)
       call broadcast_scalar (f_siflswdtop, master_task)
       call broadcast_scalar (f_siflswutop, master_task)
       call broadcast_scalar (f_siflswdbot, master_task)
@@ -1648,7 +1655,7 @@ subroutine init_hist (dt)
              "none", c1, c0,         &
              ns1, f_sidmassevapsubl)
 
-         call define_hist_field(n_sndmasssubl,"sndmassubl","kg m-2 s-1",tstr2D, tcstr,  &
+         call define_hist_field(n_sndmasssubl,"sndmasssubl","kg m-2 s-1",tstr2D, tcstr,  &
              "snow mass change from evaporation and sublimation", &
              "none", c1, c0,         &
              ns1, f_sndmasssubl)
@@ -1678,6 +1685,11 @@ subroutine init_hist (dt)
              "none", c1, c0,         &
              ns1, f_sndmassmelt)
 
+         call define_hist_field(n_sndmassdyn,"sndmassdyn","kg m-2 s-1",tstr2D, tcstr,  &
+             "snow mass change from dynamics ridging",                      &
+             "none", c1, c0,         &
+             ns1, f_sndmassdyn)
+
          call define_hist_field(n_siflswdtop,"siflswdtop","W/m2",tstr2D, tcstr, &
              "down shortwave flux over sea ice", &
              "positive downward", c1, c0,                            &
@@ -2160,7 +2172,7 @@ subroutine accum_hist (dt)
           taubxN, taubyN, strocnxN, strocnyN, &
           strairxE, strairyE, strtltxE, strtltyE, strintxE, strintyE, &
           taubxE, taubyE, strocnxE, strocnyE, &
-          fmU, fmN, fmE, daidtt, dvidtt, daidtd, dvidtd, fsurf, &
+          fmU, fmN, fmE, daidtt, dvidtt, daidtd, dvidtd, dvsdtd, fsurf, &
           fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, snwcnt, &
           stressp_1, stressm_1, stress12_1, &
           stresspT, stressmT, stress12T, &
@@ -3068,7 +3080,7 @@ subroutine accum_hist (dt)
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = evaps(i,j,iblk)*rhos
+                 worka(i,j) = evaps(i,j,iblk)
               endif
            enddo
            enddo
@@ -3080,7 +3092,7 @@ subroutine accum_hist (dt)
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*fsnow(i,j,iblk)*rhos
+                 worka(i,j) = aice(i,j,iblk)*fsnow(i,j,iblk)
               endif
            enddo
            enddo
@@ -3099,6 +3111,18 @@ subroutine accum_hist (dt)
            call accum_hist_field(n_sndmassmelt, iblk, worka(:,:), a2D)
          endif
 
+         if (f_sndmassdyn(1:1) /= 'x') then
+           worka(:,:) = c0
+           do j = jlo, jhi
+           do i = ilo, ihi
+              if (aice(i,j,iblk) > puny) then
+                 worka(i,j) = dvsdtd(i,j,iblk)*rhos
+              endif
+           enddo
+           enddo
+           call accum_hist_field(n_sndmassdyn, iblk, worka(:,:), a2D)
+         endif
+
          if (f_siflswdtop(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi

cicecore/cicedyn/analysis/ice_history_shared.F90

@@ -314,6 +314,7 @@ module ice_history_shared
            f_sidmasslat = 'x', &
            f_sndmasssnf = 'x', &
            f_sndmassmelt = 'x', &
+           f_sndmassdyn = 'x', &
            f_siflswdtop = 'x', &
            f_siflswutop = 'x', &
            f_siflswdbot = 'x', &
@@ -485,6 +486,7 @@ module ice_history_shared
            f_sidmasslat, &
            f_sndmasssnf, &
            f_sndmassmelt, &
+           f_sndmassdyn, &
            f_siflswdtop, &
            f_siflswutop, &
            f_siflswdbot, &
@@ -684,6 +686,7 @@ module ice_history_shared
            n_sidmasslat,  &
            n_sndmasssnf,  &
            n_sndmassmelt,  &
+           n_sndmassdyn,  &
            n_siflswdtop,  &
            n_siflswutop,  &
            n_siflswdbot,  &

cicecore/cicedyn/general/ice_flux.F90

@@ -107,6 +107,7 @@ module ice_flux
          strintyE, & ! divergence of internal ice stress, y at E points (N/m^2)
          daidtd  , & ! ice area tendency due to transport   (1/s)
          dvidtd  , & ! ice volume tendency due to transport (m/s)
+         dvsdtd  , & ! snow volume tendency due to transport (m/s)
          dagedtd , & ! ice age tendency due to transport (s/s)
          dardg1dt, & ! rate of area loss by ridging ice (1/s)
          dardg2dt, & ! rate of area gain by new ridges (1/s)
@@ -319,6 +320,7 @@ module ice_flux
          dsnow,  & ! change in snow thickness (m/step-->cm/day)
          daidtt, & ! ice area tendency thermo.   (s^-1)
          dvidtt, & ! ice volume tendency thermo. (m/s)
+         dvsdtt, & ! snow volume tendency thermo. (m/s)
          dagedtt,& ! ice age tendency thermo.    (s/s)
          mlt_onset, &! day of year that sfc melting begins
          frz_onset, &! day of year that freezing begins (congel or frazil)
@@ -419,6 +421,7 @@ subroutine alloc_flux
          strintyU   (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y (N/m^2)
          daidtd     (nx_block,ny_block,max_blocks), & ! ice area tendency due to transport   (1/s)
          dvidtd     (nx_block,ny_block,max_blocks), & ! ice volume tendency due to transport (m/s)
+         dvsdtd     (nx_block,ny_block,max_blocks), & ! snow volume tendency due to transport (m/s)
          dagedtd    (nx_block,ny_block,max_blocks), & ! ice age tendency due to transport (s/s)
          dardg1dt   (nx_block,ny_block,max_blocks), & ! rate of area loss by ridging ice (1/s)
          dardg2dt   (nx_block,ny_block,max_blocks), & ! rate of area gain by new ridges (1/s)
@@ -530,6 +533,7 @@ subroutine alloc_flux
          dsnow      (nx_block,ny_block,max_blocks), & ! change in snow thickness (m/step-->cm/day)
          daidtt     (nx_block,ny_block,max_blocks), & ! ice area tendency thermo.   (s^-1)
          dvidtt     (nx_block,ny_block,max_blocks), & ! ice volume tendency thermo. (m/s)
+         dvsdtt     (nx_block,ny_block,max_blocks), & ! snow volume tendency thermo. (m/s)
          dagedtt    (nx_block,ny_block,max_blocks), & ! ice age tendency thermo.    (s/s)
          mlt_onset  (nx_block,ny_block,max_blocks), & ! day of year that sfc melting begins
          frz_onset  (nx_block,ny_block,max_blocks), & ! day of year that freezing begins (congel or frazil)
@@ -918,7 +922,7 @@ end subroutine init_flux_ocn
 
       subroutine init_history_therm
 
-      use ice_state, only: aice, vice, trcr
+      use ice_state, only: aice, vice, vsno, trcr
       use ice_arrays_column, only: &
           hfreebd, hdraft, hridge, distrdg, hkeel, dkeel, lfloe, dfloe, &
           Cdn_atm_skin, Cdn_atm_floe, Cdn_atm_pond, Cdn_atm_rdg, &
@@ -965,6 +969,7 @@ subroutine init_history_therm
       meltl  (:,:,:) = c0
       daidtt (:,:,:) = aice(:,:,:) ! temporary initial area
       dvidtt (:,:,:) = vice(:,:,:) ! temporary initial volume
+      dvsdtt (:,:,:) = vsno(:,:,:) ! temporary initial volume
       if (tr_iage) then
          dagedtt(:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age
       else
@@ -1022,7 +1027,7 @@ end subroutine init_history_therm
 
       subroutine init_history_dyn
 
-      use ice_state, only: aice, vice, trcr, strength, divu, shear, vort
+      use ice_state, only: aice, vice, vsno, trcr, strength, divu, shear, vort
       use ice_grid,  only: grid_ice
 
       logical (kind=log_kind) :: &
@@ -1061,6 +1066,7 @@ subroutine init_history_dyn
       opening (:,:,:) = c0
       daidtd  (:,:,:) = aice(:,:,:) ! temporary initial area
       dvidtd  (:,:,:) = vice(:,:,:) ! temporary initial volume
+      dvsdtd  (:,:,:) = vsno(:,:,:) ! temporary initial volume
       if (tr_iage) &
          dagedtd (:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age
       fmU     (:,:,:) = c0

cicecore/cicedyn/general/ice_step_mod.F90

@@ -750,7 +750,7 @@ end subroutine step_therm2
 !
 ! authors: Elizabeth Hunke, LANL
 
-      subroutine update_state (dt, daidt, dvidt, dagedt, offset)
+      subroutine update_state (dt, daidt, dvidt, dvsdt, dagedt, offset)
 
       use ice_domain_size, only: ncat
 !     use ice_grid, only: tmask
@@ -766,6 +766,7 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset)
       real (kind=dbl_kind), dimension(:,:,:), intent(inout), optional :: &
          daidt, & ! change in ice area per time step
          dvidt, & ! change in ice volume per time step
+         dvsdt, & ! change in snow volume per time step
          dagedt   ! change in ice age per time step
 
       real (kind=dbl_kind), intent(in), optional :: &
@@ -827,25 +828,26 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset)
                                    nt_strata     = nt_strata(:,:),   &
                                    Tf            = Tf(i,j,iblk))
 
-         if (present(offset)) then
-
-      !-----------------------------------------------------------------
-      ! Compute thermodynamic area and volume tendencies.
-      !-----------------------------------------------------------------
-
-         daidt(i,j,iblk) = (aice(i,j,iblk) - daidt(i,j,iblk)) / dt
-         dvidt(i,j,iblk) = (vice(i,j,iblk) - dvidt(i,j,iblk)) / dt
-         if (tr_iage) then
-            if (offset > c0) then                 ! thermo
-               if (trcr(i,j,nt_iage,iblk) > c0) &
-               dagedt(i,j,iblk) = (trcr(i,j,nt_iage,iblk) &
-                                - dagedt(i,j,iblk) - offset) / dt
-            else                                  ! dynamics
-               dagedt(i,j,iblk) = (trcr(i,j,nt_iage,iblk) &
-                                - dagedt(i,j,iblk)) / dt
-            endif
-         endif ! tr_iage
-         endif ! present(offset)
+            if (present(offset)) then
+
+            !-----------------------------------------------------------------
+            ! Compute thermodynamic area and volume tendencies.
+            !-----------------------------------------------------------------
+
+               if (present(daidt)) daidt(i,j,iblk) = (aice(i,j,iblk) - daidt(i,j,iblk)) / dt
+               if (present(dvidt)) dvidt(i,j,iblk) = (vice(i,j,iblk) - dvidt(i,j,iblk)) / dt
+               if (present(dvsdt)) dvsdt(i,j,iblk) = (vsno(i,j,iblk) - dvsdt(i,j,iblk)) / dt
+               if (present(dagedt) .and. tr_iage) then
+                  if (offset > c0) then                 ! thermo
+                     if (trcr(i,j,nt_iage,iblk) > c0) &
+                     dagedt(i,j,iblk) = (trcr(i,j,nt_iage,iblk) &
+                                      - dagedt(i,j,iblk) - offset) / dt
+                  else                                  ! dynamics
+                     dagedt(i,j,iblk) = (trcr(i,j,nt_iage,iblk) &
+                                      - dagedt(i,j,iblk)) / dt
+                  endif
+               endif ! tr_iage
+            endif ! present(offset)
 
          enddo ! i
          enddo ! j

cicecore/drivers/direct/hadgem3/CICE_RunMod.F90

@@ -140,7 +140,7 @@ subroutine ice_step
       use ice_dyn_eap, only: write_restart_eap
       use ice_dyn_shared, only: kdyn
       use ice_flux, only: scale_factor, init_history_therm, &
-          daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd
+          daidtt, daidtd, dvidtt, dvidtd, dvsdtt, dvsdtd, dagedtt, dagedtd
       use ice_history, only: accum_hist
       use ice_history_bgc, only: init_history_bgc
       use ice_restart, only: final_restart
@@ -238,7 +238,8 @@ subroutine ice_step
 
          ! clean up, update tendency diagnostics
          offset = dt
-         call update_state (dt, daidtt, dvidtt, dagedtt, offset)
+         call update_state (dt=dt, daidt=daidtt, dvidt=dvidtt, dvsdt=dvsdtt, &
+                            dagedt=dagedtt, offset=offset)
 
          call ice_timer_stop(timer_thermo) ! thermodynamics
          call ice_timer_stop(timer_column) ! column physics
@@ -265,7 +266,8 @@ subroutine ice_step
 
             ! clean up, update tendency diagnostics
             offset = c0
-            call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
+            call update_state (dt=dt_dyn, daidt=daidtd, dvidt=dvidtd, dvsdt=dvsdtd, &
+                               dagedt=dagedtd, offset=offset)
 
          enddo
 

cicecore/drivers/direct/nemo_concepts/CICE_RunMod.F90

@@ -140,7 +140,7 @@ subroutine ice_step
       use ice_dyn_eap, only: write_restart_eap
       use ice_dyn_shared, only: kdyn
       use ice_flux, only: scale_factor, init_history_therm, &
-          daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd
+          daidtt, daidtd, dvidtt, dvidtd, dvsdtt, dvsdtd, dagedtt, dagedtd
       use ice_history, only: accum_hist
       use ice_history_bgc, only: init_history_bgc
       use ice_restart, only: final_restart
@@ -238,7 +238,8 @@ subroutine ice_step
 
          ! clean up, update tendency diagnostics
          offset = dt
-         call update_state (dt, daidtt, dvidtt, dagedtt, offset)
+         call update_state (dt=dt, daidt=daidtt, dvidt=dvidtt, dvsdt=dvsdtt, &
+                            dagedt=dagedtt, offset=offset)
 
          call ice_timer_stop(timer_thermo) ! thermodynamics
          call ice_timer_stop(timer_column) ! column physics
@@ -265,7 +266,8 @@ subroutine ice_step
 
             ! clean up, update tendency diagnostics
             offset = c0
-            call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
+            call update_state (dt=dt_dyn, daidt=daidtd, dvidt=dvidtd, dvsdt=dvsdtd, &
+                               dagedt=dagedtd, offset=offset)
 
          enddo
 

cicecore/drivers/mct/cesm1/CICE_RunMod.F90

@@ -132,7 +132,7 @@ subroutine ice_step
       use ice_dyn_eap, only: write_restart_eap
       use ice_dyn_shared, only: kdyn, kridge
       use ice_flux, only: scale_factor, init_history_therm, &
-          daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd
+          daidtt, daidtd, dvidtt, dvidtd, dvsdtt, dvsdtd, dagedtt, dagedtd
       use ice_history, only: accum_hist
       use ice_history_bgc, only: init_history_bgc
       use ice_restart, only: final_restart
@@ -261,7 +261,8 @@ subroutine ice_step
 
          ! clean up, update tendency diagnostics
          offset = dt
-         call update_state (dt, daidtt, dvidtt, dagedtt, offset)
+         call update_state (dt=dt, daidt=daidtt, dvidt=dvidtt, dvsdt=dvsdtt, &
+                            dagedt=dagedtt, offset=offset)
 
          call ice_timer_stop(timer_thermo) ! thermodynamics
          call ice_timer_stop(timer_column) ! column physics
@@ -302,7 +303,8 @@ subroutine ice_step
 
             ! clean up, update tendency diagnostics
             offset = c0
-            call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
+            call update_state (dt=dt_dyn, daidt=daidtd, dvidt=dvidtd, dvsdt=dvsdtd, &
+                               dagedt=dagedtd, offset=offset)
 
          enddo
 
@@ -326,7 +328,7 @@ subroutine ice_step
                call step_snow (dt, iblk)
             enddo
             !$OMP END PARALLEL DO
-            call update_state (dt) ! clean up
+            call update_state (dt=dt) ! clean up
          endif
 
          !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)

cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90

@@ -119,7 +119,7 @@ subroutine ice_step
       use ice_dyn_eap, only: write_restart_eap
       use ice_dyn_shared, only: kdyn, kridge
       use ice_flux, only: scale_factor, init_history_therm, &
-          daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd
+          daidtt, daidtd, dvidtt, dvidtd, dvsdtt, dvsdtd, dagedtt, dagedtd
       use ice_history, only: accum_hist
       use ice_history_bgc, only: init_history_bgc
       use ice_restart, only: final_restart
@@ -250,7 +250,8 @@ subroutine ice_step
 
          ! clean up, update tendency diagnostics
          offset = dt
-         call update_state (dt, daidtt, dvidtt, dagedtt, offset)
+         call update_state (dt=dt, daidt=daidtt, dvidt=dvidtt, dvsdt=dvsdtt, &
+                            dagedt=dagedtt, offset=offset)
 
          call ice_timer_stop(timer_thermo) ! thermodynamics
          call ice_timer_stop(timer_column) ! column physics
@@ -291,7 +292,8 @@ subroutine ice_step
 
             ! clean up, update tendency diagnostics
             offset = c0
-            call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
+            call update_state (dt=dt_dyn, daidt=daidtd, dvidt=dvidtd, dvsdt=dvsdtd, &
+                               dagedt=dagedtd, offset=offset)
 
          enddo
          if (debug_model) then
@@ -318,7 +320,7 @@ subroutine ice_step
             do iblk = 1, nblocks
                call step_snow (dt, iblk)
             enddo
-            call update_state (dt) ! clean up
+            call update_state (dt=dt) ! clean up
          endif
 
 !MHRI: CHECK THIS OMP