Advanced snow physics driver

cicecore/cicedynB/analysis/ice_diagnostics.F90

@@ -14,6 +14,7 @@ module ice_diagnostics
       use ice_communicate, only: my_task, master_task
       use ice_constants, only: c0, c1
       use ice_calendar, only: istep1
+      use ice_domain_size, only: nslyr
       use ice_fileunits, only: nu_diag
       use ice_fileunits, only: flush_fileunit
       use ice_exit, only: abort_ice
@@ -142,15 +143,19 @@ subroutine runtime_diags (dt)
          i, j, k, n, iblk, nc, &
          ktherm, &
          nt_tsfc, nt_aero, nt_fbri, nt_apnd, nt_hpnd, nt_fsd, &
-         nt_isosno, nt_isoice
+         nt_isosno, nt_isoice, nt_rsnw, nt_rhos, nt_smice, nt_smliq
 
       logical (kind=log_kind) :: &
-         tr_pond_topo, tr_brine, tr_iso, tr_aero, calc_Tsfc, tr_fsd
+         tr_pond_topo, tr_brine, tr_iso, tr_aero, calc_Tsfc, tr_fsd, &
+         tr_snow, snwgrain
 
       real (kind=dbl_kind) :: &
          rhow, rhos, rhoi, puny, awtvdr, awtidr, awtvdf, awtidf, &
          rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh
 
+      character (len=char_len) :: &
+         snwredist
+
       ! hemispheric state quantities
       real (kind=dbl_kind) :: &
          umaxn,   hmaxn,   shmaxn,    arean,   snwmxn, extentn, shmaxnt, &
@@ -190,7 +195,8 @@ subroutine runtime_diags (dt)
          pTsfc, pevap, pfswabs, pflwout, pflat, pfsens, &
          pfsurf, pfcondtop, psst, psss, pTf, hiavg, hsavg, hbravg, &
          pfhocn, psalt, fsdavg, &
-         pmeltt, pmeltb, pmeltl, psnoice, pdsnow, pfrazil, pcongel
+         pmeltt, pmeltb, pmeltl, psnoice, pdsnow, pfrazil, pcongel, &
+         prsnwavg, prhosavg, psmicetot, psmliqtot, psmtot
 
       real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
          work1, work2
@@ -199,15 +205,19 @@ subroutine runtime_diags (dt)
 
       call icepack_query_parameters(ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc)
       call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_aero_out=tr_aero, &
-           tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
+           tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, &
+           tr_snow_out=tr_snow)
       call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc, &
            nt_aero_out=nt_aero, nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, &
-           nt_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
+           nt_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
+           nt_rsnw_out=nt_rsnw, nt_rhos_out=nt_rhos, &
+           nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
       call icepack_query_parameters(Tffresh_out=Tffresh, rhos_out=rhos, &
            rhow_out=rhow, rhoi_out=rhoi, puny_out=puny, &
            awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, &
            rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, &
-           ice_ref_salinity_out=ice_ref_salinity)
+           ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, &
+           snwgrain_out=snwgrain)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
@@ -825,6 +835,26 @@ subroutine runtime_diags (dt)
                      enddo
                   endif
                endif
+               if (tr_snow) then      ! snow tracer quantities
+                  prsnwavg (n) = c0   ! avg snow grain radius
+                  prhosavg (n) = c0   ! avg snow density
+                  psmicetot(n) = c0   ! total mass of ice in snow (kg/m2)
+                  psmliqtot(n) = c0   ! total mass of liquid in snow (kg/m2)
+                  psmtot   (n) = c0   ! total mass of snow volume (kg/m2)
+                  if (vsno(i,j,iblk) > c0) then
+                     do k = 1, nslyr
+                        prsnwavg (n) = prsnwavg (n) + trcr(i,j,nt_rsnw +k-1,iblk) ! snow grain radius
+                        prhosavg (n) = prhosavg (n) + trcr(i,j,nt_rhos +k-1,iblk) ! compacted snow density
+                        psmicetot(n) = psmicetot(n) + trcr(i,j,nt_smice+k-1,iblk) * vsno(i,j,iblk)
+                        psmliqtot(n) = psmliqtot(n) + trcr(i,j,nt_smliq+k-1,iblk) * vsno(i,j,iblk)
+                     end do
+                  endif
+                  psmtot   (n) = rhos * vsno(i,j,iblk) ! mass of ice in standard density snow
+                  prsnwavg (n) = prsnwavg (n) / real(nslyr,kind=dbl_kind) ! snow grain radius
+                  prhosavg (n) = prhosavg (n) / real(nslyr,kind=dbl_kind) ! compacted snow density
+                  psmicetot(n) = psmicetot(n) / real(nslyr,kind=dbl_kind) ! mass of ice in snow
+                  psmliqtot(n) = psmliqtot(n) / real(nslyr,kind=dbl_kind) ! mass of liquid in snow
+               end if
                psalt(n) = c0
                if (vice(i,j,iblk) /= c0) psalt(n) = work2(i,j,iblk)/vice(i,j,iblk)
                pTsfc(n) = trcr(i,j,nt_Tsfc,iblk)   ! ice/snow sfc temperature
@@ -877,6 +907,11 @@ subroutine runtime_diags (dt)
             call broadcast_scalar(pmeltl   (n), pmloc(n))
             call broadcast_scalar(psnoice  (n), pmloc(n))
             call broadcast_scalar(pdsnow   (n), pmloc(n))
+            call broadcast_scalar(psmtot   (n), pmloc(n))
+            call broadcast_scalar(prsnwavg (n), pmloc(n))
+            call broadcast_scalar(prhosavg (n), pmloc(n))
+            call broadcast_scalar(psmicetot(n), pmloc(n))
+            call broadcast_scalar(psmliqtot(n), pmloc(n))
             call broadcast_scalar(pfrazil  (n), pmloc(n))
             call broadcast_scalar(pcongel  (n), pmloc(n))
             call broadcast_scalar(pdhi     (n), pmloc(n))
@@ -1060,6 +1095,26 @@ subroutine runtime_diags (dt)
         write(nu_diag,900) 'effective dhi (m)      = ',pdhi(1),pdhi(2)
         write(nu_diag,900) 'effective dhs (m)      = ',pdhs(1),pdhs(2)
         write(nu_diag,900) 'intnl enrgy chng(W/m^2)= ',pde (1),pde (2)
+
+        if (tr_snow) then
+           if (trim(snwredist) /= 'none') then
+              write(nu_diag,900) 'avg snow density(kg/m3)= ',prhosavg(1) &
+                                                            ,prhosavg(2)
+           endif
+           if (snwgrain) then
+              write(nu_diag,900) 'avg snow grain radius  = ',prsnwavg(1) &
+                                                            ,prsnwavg(2)
+              write(nu_diag,900) 'mass ice in snow(kg/m2)= ',psmicetot(1) &
+                                                            ,psmicetot(2)
+              write(nu_diag,900) 'mass liq in snow(kg/m2)= ',psmliqtot(1) &
+                                                            ,psmliqtot(2)
+              write(nu_diag,900) 'mass std snow   (kg/m2)= ',psmtot(1) &
+                                                            ,psmtot(2)
+              write(nu_diag,900) 'max  ice+liq    (kg/m2)= ',rhow * hsavg(1) &
+                                                            ,rhow * hsavg(2)
+           endif
+        endif
+
         write(nu_diag,*) '----------ocn----------'
         write(nu_diag,900) 'sst (C)                = ',psst(1),psst(2)
         write(nu_diag,900) 'sss (ppt)              = ',psss(1),psss(2)
@@ -1597,19 +1652,21 @@ subroutine print_state(plabel,i,j,iblk)
            rad_to_deg, puny, rhoi, lfresh, rhos, cp_ice
 
       integer (kind=int_kind) :: n, k, nt_Tsfc, nt_qice, nt_qsno, nt_fsd, &
-           nt_isosno, nt_isoice, nt_sice
+           nt_isosno, nt_isoice, nt_sice, nt_smice, nt_smliq
 
-      logical (kind=log_kind) :: tr_fsd, tr_iso
+      logical (kind=log_kind) :: tr_fsd, tr_iso, tr_snow
 
       type (block) :: &
          this_block           ! block information for current block
 
       character(len=*), parameter :: subname = '(print_state)'
 
-      call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
+      call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, &
+           tr_snow_out=tr_snow)
       call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
            nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, &
-           nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
+           nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
+           nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
       call icepack_query_parameters( &
            rad_to_deg_out=rad_to_deg, puny_out=puny, rhoi_out=rhoi, lfresh_out=lfresh, &
            rhos_out=rhos, cp_ice_out=cp_ice)
@@ -1639,8 +1696,11 @@ subroutine print_state(plabel,i,j,iblk)
          endif
          write(nu_diag,*) 'Tsfcn',trcrn(i,j,nt_Tsfc,n,iblk)
          if (tr_fsd) write(nu_diag,*) 'afsdn',trcrn(i,j,nt_fsd,n,iblk) ! fsd cat 1
-!         if (tr_iso) write(nu_diag,*) 'isosno',trcrn(i,j,nt_isosno,n,iblk) ! isotopes in snow
-!         if (tr_iso) write(nu_diag,*) 'isoice',trcrn(i,j,nt_isoice,n,iblk) ! isotopes in ice
+! layer 1 diagnostics
+!         if (tr_iso)  write(nu_diag,*) 'isosno',trcrn(i,j,nt_isosno,n,iblk) ! isotopes in snow
+!         if (tr_iso)  write(nu_diag,*) 'isoice',trcrn(i,j,nt_isoice,n,iblk) ! isotopes in ice
+!         if (tr_snow) write(nu_diag,*) 'smice', trcrn(i,j,nt_smice, n,iblk) ! ice mass in snow
+!         if (tr_snow) write(nu_diag,*) 'smliq', trcrn(i,j,nt_smliq, n,iblk) ! liquid mass in snow
          write(nu_diag,*) ' '
 
 ! dynamics (transport and/or ridging) causes the floe size distribution to become non-normal

cicecore/cicedynB/analysis/ice_diagnostics_bgc.F90

@@ -937,19 +937,18 @@ subroutine zsal_diags
                enddo
                if (aice(i,j,iblk) > c0) &
                   psice_rho(n) = psice_rho(n)/aice(i,j,iblk)
-               if (tr_brine .and. aice(i,j,iblk) > c0) &
+               if (tr_brine .and. aice(i,j,iblk) > c0) then
                   phinS(n) = trcr(i,j,nt_fbri,iblk)*vice(i,j,iblk)/aice(i,j,iblk)
-
-               if (aicen(i,j,1,iblk)> c0) then
-                  if (tr_brine) phinS1(n) = trcrn(i,j,nt_fbri,1,iblk) &
-                                          * vicen(i,j,1,iblk)/aicen(i,j,1,iblk)
+                  phbrn(n) = (c1 - rhosi/rhow)*vice(i,j,iblk)/aice(i,j,iblk) &
+                                 - rhos/rhow  *vsno(i,j,iblk)/aice(i,j,iblk)
+               endif
+               if (tr_brine .and. aicen(i,j,1,iblk)> c0) then
+                  phinS1(n) = trcrn(i,j,nt_fbri,1,iblk) &
+                            * vicen(i,j,1,iblk)/aicen(i,j,1,iblk)
                   pdh_top1(n) = dhbr_top(i,j,1,iblk)
                   pdh_bot1(n) = dhbr_bot(i,j,1,iblk)
                   pdarcy_V1(n) = darcy_V(i,j,1,iblk)
                endif  
-               if (tr_brine .AND. aice(i,j,iblk) > c0) &
-                  phbrn(n) = (c1 - rhosi/rhow)*vice(i,j,iblk)/aice(i,j,iblk) &
-                                 - rhos/rhow  *vsno(i,j,iblk)/aice(i,j,iblk)
                do k = 1, nblyr+1
                   pbTiz(n,k) = c0
                   piDin(n,k) = c0

cicecore/cicedynB/analysis/ice_history.F90

@@ -67,10 +67,11 @@ subroutine init_hist (dt)
           histfreq_n, nstreams
       use ice_domain_size, only: max_blocks, max_nstrm, nilyr, nslyr, nblyr, ncat, nfsd
       use ice_dyn_shared, only: kdyn
-      use ice_flux, only: mlt_onset, frz_onset, albcnt
+      use ice_flux, only: mlt_onset, frz_onset, albcnt, snwcnt
       use ice_history_shared ! everything
       use ice_history_mechred, only: init_hist_mechred_2D, init_hist_mechred_3Dc
       use ice_history_pond, only: init_hist_pond_2D, init_hist_pond_3Dc
+      use ice_history_snow, only: init_hist_snow_2D, init_hist_snow_3Dc
       use ice_history_bgc, only:init_hist_bgc_2D, init_hist_bgc_3Dc, &
           init_hist_bgc_3Db, init_hist_bgc_3Da
       use ice_history_drag, only: init_hist_drag_2D
@@ -86,7 +87,7 @@ subroutine init_hist (dt)
       real (kind=dbl_kind) :: rhofresh, Tffresh, secday, rad_to_deg
       logical (kind=log_kind) :: formdrag
       logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_brine
-      logical (kind=log_kind) :: tr_fsd
+      logical (kind=log_kind) :: tr_fsd, tr_snow
       logical (kind=log_kind) :: skl_bgc, solve_zsal, solve_zbgc, z_tracers
       integer (kind=int_kind) :: n, ns, ns1, ns2
       integer (kind=int_kind), dimension(max_nstrm) :: &
@@ -115,7 +116,7 @@ subroutine init_hist (dt)
          solve_zsal_out=solve_zsal, solve_zbgc_out=solve_zbgc, z_tracers_out=z_tracers)
       call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
          tr_lvl_out=tr_lvl, tr_pond_out=tr_pond, tr_aero_out=tr_aero, &
-         tr_brine_out=tr_brine, tr_fsd_out=tr_fsd)
+         tr_brine_out=tr_brine, 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__)
@@ -1426,6 +1427,9 @@ subroutine init_hist (dt)
       ! floe size distribution
       call init_hist_fsd_2D
 
+      ! advanced snow physics
+      call init_hist_snow_2D (dt)
+
       !-----------------------------------------------------------------
       ! 3D (category) variables looped separately for ordering
       !-----------------------------------------------------------------
@@ -1501,6 +1505,9 @@ subroutine init_hist (dt)
       ! biogeochemistry
       call init_hist_bgc_3Dc
 
+      ! advanced snow physics
+      call init_hist_snow_3Dc
+
       !-----------------------------------------------------------------
       ! 3D (vertical) variables must be looped separately
       !-----------------------------------------------------------------
@@ -1688,6 +1695,7 @@ subroutine init_hist (dt)
       if (allocated(a4Df)) a4Df(:,:,:,:,:,:) = c0
       avgct(:) = c0
       albcnt(:,:,:,:) = c0
+      snwcnt(:,:,:,:) = c0
 
       if (restart .and. yday >= c2) then
 ! restarting midyear gives erroneous onset dates
@@ -1726,7 +1734,7 @@ subroutine accum_hist (dt)
           fhocn, fhocn_ai, uatm, vatm, fbot, Tbot, Tsnice, &
           fswthru_ai, strairx, strairy, strtltx, strtlty, strintx, strinty, &
           taubx, tauby, strocnx, strocny, fm, daidtt, dvidtt, daidtd, dvidtd, fsurf, &
-          fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, &
+          fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, snwcnt, &
           stressp_1, stressm_1, stress12_1, &
           stressp_2, &
           stressp_3, &
@@ -1739,6 +1747,8 @@ subroutine accum_hist (dt)
       use ice_history_bgc, only: accum_hist_bgc
       use ice_history_mechred, only: accum_hist_mechred
       use ice_history_pond, only: accum_hist_pond
+      use ice_history_snow, only: accum_hist_snow, &
+          f_rhos_cmp, f_rhos_cnt, n_rhos_cmp, n_rhos_cnt
       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
@@ -1775,7 +1785,7 @@ subroutine accum_hist (dt)
       real (kind=dbl_kind) :: Tffresh, rhoi, rhos, rhow, ice_ref_salinity
       real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt
       logical (kind=log_kind) :: formdrag, skl_bgc
-      logical (kind=log_kind) :: tr_pond, tr_aero, tr_brine
+      logical (kind=log_kind) :: tr_pond, tr_aero, tr_brine, tr_snow
       integer (kind=int_kind) :: ktherm
       integer (kind=int_kind) :: nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY, nt_Tsfc, &
                                  nt_alvl, nt_vlvl
@@ -1791,7 +1801,7 @@ subroutine accum_hist (dt)
            rhow_out=rhow, ice_ref_salinity_out=ice_ref_salinity)
       call icepack_query_parameters(formdrag_out=formdrag, skl_bgc_out=skl_bgc, ktherm_out=ktherm)
       call icepack_query_tracer_flags(tr_pond_out=tr_pond, tr_aero_out=tr_aero, &
-           tr_brine_out=tr_brine)
+           tr_brine_out=tr_brine, tr_snow_out=tr_snow)
       call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_qice_out=nt_qice, &
            nt_qsno_out=nt_qsno, nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_Tsfc_out=nt_Tsfc, &
            nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl)
@@ -3040,6 +3050,9 @@ subroutine accum_hist (dt)
          ! floe size distribution
          call accum_hist_fsd (iblk)
 
+         ! advanced snow physics
+         call accum_hist_snow (iblk)
+
       enddo                     ! iblk
       !$OMP END PARALLEL DO
 
@@ -3669,7 +3682,38 @@ subroutine accum_hist (dt)
               enddo             ! j
               endif
 
-              endif
+! snwcnt averaging is not working correctly
+! for now, these history fields will have zeroes includes in the averages
+!              if (avail_hist_fields(n)%vname(1:8) == 'rhos_cmp') then
+!              do j = jlo, jhi
+!              do i = ilo, ihi
+!                 if (tmask(i,j,iblk)) then
+!                    ravgctz = c0
+!                    if (snwcnt(i,j,iblk,ns) > puny) &
+!                        ravgctz = c1/snwcnt(i,j,iblk,ns)
+!                    if (f_rhos_cmp (1:1) /= 'x' .and. n_rhos_cmp(ns) /= 0) &
+!                       a2D(i,j,n_rhos_cmp(ns),iblk) = &
+!                       a2D(i,j,n_rhos_cmp(ns),iblk)*avgct(ns)*ravgctz
+!                 endif
+!              enddo             ! i
+!              enddo             ! j
+!              endif
+!              if (avail_hist_fields(n)%vname(1:8) == 'rhos_cnt') then
+!              do j = jlo, jhi
+!              do i = ilo, ihi
+!                 if (tmask(i,j,iblk)) then
+!                    ravgctz = c0
+!                    if (snwcnt(i,j,iblk,ns) > puny) &
+!                        ravgctz = c1/snwcnt(i,j,iblk,ns)
+!                    if (f_rhos_cnt (1:1) /= 'x' .and. n_rhos_cnt(ns) /= 0) &
+!                       a2D(i,j,n_rhos_cnt(ns),iblk) = &
+!                       a2D(i,j,n_rhos_cnt(ns),iblk)*avgct(ns)*ravgctz
+!                 endif
+!              enddo             ! i
+!              enddo             ! j
+!              endif
+
+              endif             ! avail_hist_fields(n)%vhistfreq == histfreq(ns)
            enddo                ! n
 
            do n = 1, num_avail_hist_fields_3Dc
@@ -3992,10 +4036,12 @@ subroutine accum_hist (dt)
            if (allocated(a4Df)) a4Df(:,:,:,:,:,:) = c0
            avgct(:) = c0
            albcnt(:,:,:,:) = c0
+           snwcnt(:,:,:,:) = c0
            write_ic = .false.        ! write initial condition once at most
         else
            avgct(ns) = c0
            albcnt(:,:,:,ns) = c0
+           snwcnt(:,:,:,ns) = c0
         endif
 !        if (write_history(ns)) albcnt(:,:,:,ns) = c0
 

cicecore/cicedynB/analysis/ice_history_fsd.F90

@@ -303,7 +303,7 @@ subroutine accum_hist_fsd (iblk)
 
       integer (kind=int_kind) :: &
          i, j, n, k, & ! loop indices
-         nt_fsd        ! ! fsd tracer index
+         nt_fsd        ! fsd tracer index
       logical (kind=log_kind) :: tr_fsd
       real (kind=dbl_kind) :: floeshape, puny
 

cicecore/cicedynB/analysis/ice_history_pond.F90

@@ -75,15 +75,15 @@ subroutine init_hist_pond_2D
       logical (kind=log_kind) :: tr_pond
       character(len=*), parameter :: subname = '(init_hist_pond_2D)'
 
-      !-----------------------------------------------------------------
-      ! read namelist
-      !-----------------------------------------------------------------
-
       call icepack_query_tracer_flags(tr_pond_out=tr_pond)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
+      !-----------------------------------------------------------------
+      ! read namelist
+      !-----------------------------------------------------------------
+
       call get_fileunit(nu_nml)
       if (my_task == master_task) then
          open (nu_nml, file=nml_filename, status='old',iostat=nml_error)

cicecore/cicedynB/analysis/ice_history_shared.F90

@@ -59,7 +59,7 @@ module ice_history_shared
 
       !---------------------------------------------------------------
       ! Instructions for adding a field: (search for 'example')
-      !     Here:
+      !     Here or in ice_history_[process].F90:
       ! (1) Add to frequency flags (f_<field>)
       ! (2) Add to namelist (here and also in ice_in)
       ! (3) Add to index list