From 2a692af211a0577927d98d3546726c34caeaccf0 Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Thu, 12 Aug 2021 15:45:27 -0600 Subject: [PATCH] Advanced snow physics driver (#621) * initial snow redist/metamorph tracer implementation (no physics) * add rhosmin, snwlvlfac to namelist, change init_snow to init_snowtracers * updating format labels due to changes in Consortium master * add calls to new snow physics routines * update_state after step_snow * fix optional in update_state * wrap step_snow and update_state with tr_snow conditional * snow interactions with thermo (except dEdd) * conflict warning * adding snwgrain namelist flag * use snwredist=bulk instead of 30percent, fix array errors, init snow tracer values, add namelist flags to ice_in * aging lookup table for dry metamorphism (version for testing) * add/clean up diagnostics, begin debugging restart implementation * fixing restarts * add snow tests to base_suite * add history; minor cleanup * fix step_therm1 interface * io fixes for snow1 * bug fix for snow diagnostics * fix line continuation * refactor bgc diagnostics as in Icepack * add snowtable implementation * cleaning up snow test options * send rsnw tracer to dEdd; initialize snow tracers; update -s options for snow * broadcast new snow diagnostics * new snow tests require 5 snow layers and therefore can not start from standard initial condition files with nslyr=1 * update icepack * use local variable for rsnow to prevent array out-of-bound error; cleanup * updating documentation and some code comments * documentation: index * fix history output; include but comment out code for averaging based on snow presence; set rsnw_fall in ice_in for default dEdd and include relevant snow parameters in options files; update doc; clean up * minor cleanup, especially documentation * update icepack Co-authored-by: Andrew Roberts Co-authored-by: Elizabeth Hunke Co-authored-by: apcraig --- .../cicedynB/analysis/ice_diagnostics.F90 | 84 +++- .../cicedynB/analysis/ice_diagnostics_bgc.F90 | 15 +- cicecore/cicedynB/analysis/ice_history.F90 | 60 ++- .../cicedynB/analysis/ice_history_fsd.F90 | 2 +- .../cicedynB/analysis/ice_history_pond.F90 | 8 +- .../cicedynB/analysis/ice_history_shared.F90 | 2 +- .../cicedynB/analysis/ice_history_snow.F90 | 430 ++++++++++++++++++ .../dynamics/ice_transport_driver.F90 | 25 +- cicecore/cicedynB/general/ice_flux.F90 | 7 + cicecore/cicedynB/general/ice_forcing.F90 | 209 ++++++++- cicecore/cicedynB/general/ice_init.F90 | 254 ++++++++++- cicecore/cicedynB/general/ice_step_mod.F90 | 220 +++++++-- .../infrastructure/ice_read_write.F90 | 308 ++++++++++++- .../io/io_binary/ice_restart.F90 | 61 ++- .../io/io_netcdf/ice_restart.F90 | 15 +- .../infrastructure/io/io_pio2/ice_restart.F90 | 15 +- .../drivers/standalone/cice/CICE_InitMod.F90 | 52 ++- .../drivers/standalone/cice/CICE_RunMod.F90 | 27 +- cicecore/shared/ice_arrays_column.F90 | 11 + cicecore/shared/ice_fileunits.F90 | 6 + cicecore/shared/ice_init_column.F90 | 77 +++- cicecore/shared/ice_restart_column.F90 | 91 +++- configuration/scripts/ice_in | 39 ++ .../scripts/options/set_nml.snw30percent | 5 + .../scripts/options/set_nml.snwITDrdg | 10 + .../scripts/options/set_nml.snwgrain | 15 + configuration/scripts/tests/base_suite.ts | 3 + doc/source/cice_index.rst | 34 +- doc/source/developer_guide/dg_driver.rst | 11 +- doc/source/science_guide/sg_tracers.rst | 6 +- doc/source/user_guide/ug_case_settings.rst | 35 ++ icepack | 2 +- 32 files changed, 1993 insertions(+), 146 deletions(-) create mode 100644 cicecore/cicedynB/analysis/ice_history_snow.F90 create mode 100644 configuration/scripts/options/set_nml.snw30percent create mode 100644 configuration/scripts/options/set_nml.snwITDrdg create mode 100644 configuration/scripts/options/set_nml.snwgrain diff --git a/cicecore/cicedynB/analysis/ice_diagnostics.F90 b/cicecore/cicedynB/analysis/ice_diagnostics.F90 index 760952fc5..d4e7066fb 100644 --- a/cicecore/cicedynB/analysis/ice_diagnostics.F90 +++ b/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 diff --git a/cicecore/cicedynB/analysis/ice_diagnostics_bgc.F90 b/cicecore/cicedynB/analysis/ice_diagnostics_bgc.F90 index fa965dfe0..74485a5e2 100644 --- a/cicecore/cicedynB/analysis/ice_diagnostics_bgc.F90 +++ b/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 diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index f91562449..87bcba3f9 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/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 diff --git a/cicecore/cicedynB/analysis/ice_history_fsd.F90 b/cicecore/cicedynB/analysis/ice_history_fsd.F90 index 43afc1e27..7ad81e7d2 100644 --- a/cicecore/cicedynB/analysis/ice_history_fsd.F90 +++ b/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 diff --git a/cicecore/cicedynB/analysis/ice_history_pond.F90 b/cicecore/cicedynB/analysis/ice_history_pond.F90 index de10eb9fb..182865fec 100644 --- a/cicecore/cicedynB/analysis/ice_history_pond.F90 +++ b/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) diff --git a/cicecore/cicedynB/analysis/ice_history_shared.F90 b/cicecore/cicedynB/analysis/ice_history_shared.F90 index 52d268990..9297d78c9 100644 --- a/cicecore/cicedynB/analysis/ice_history_shared.F90 +++ b/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_) ! (2) Add to namelist (here and also in ice_in) ! (3) Add to index list diff --git a/cicecore/cicedynB/analysis/ice_history_snow.F90 b/cicecore/cicedynB/analysis/ice_history_snow.F90 new file mode 100644 index 000000000..5a590af2b --- /dev/null +++ b/cicecore/cicedynB/analysis/ice_history_snow.F90 @@ -0,0 +1,430 @@ +!======================================================================= + +! Snow tracer history output + + module ice_history_snow + + use ice_kinds_mod + use ice_constants, only: c0, c1, mps_to_cmpdy + use ice_domain_size, only: max_nstrm, nslyr + use ice_fileunits, only: nu_nml, nml_filename, & + get_fileunit, release_fileunit + use ice_fileunits, only: nu_diag + use ice_exit, only: abort_ice + use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted + use icepack_intfc, only: icepack_query_parameters, & + icepack_query_tracer_flags, icepack_query_tracer_indices + + implicit none + private + public :: accum_hist_snow, init_hist_snow_2D, init_hist_snow_3Dc + + !--------------------------------------------------------------- + ! flags: write to output file if true or histfreq value + !--------------------------------------------------------------- + + character (len=max_nstrm), public :: & + f_smassice = 'm', f_smassicen = 'x', & + f_smassliq = 'm', f_smassliqn = 'x', & + f_rhos_cmp = 'm', f_rhos_cmpn = 'x', & + f_rhos_cnt = 'm', f_rhos_cntn = 'x', & + f_rsnw = 'm', f_rsnwn = 'x', & + f_meltsliq = 'm', f_fsloss = 'x' + + !--------------------------------------------------------------- + ! namelist variables + !--------------------------------------------------------------- + + namelist / icefields_snow_nml / & + f_smassice, f_smassicen, & + f_smassliq, f_smassliqn, & + f_rhos_cmp, f_rhos_cmpn, & + f_rhos_cnt, f_rhos_cntn, & + f_rsnw, f_rsnwn, & + f_meltsliq, f_fsloss + + !--------------------------------------------------------------- + ! field indices + !--------------------------------------------------------------- + + integer (kind=int_kind), dimension(max_nstrm), public :: & + n_smassice, n_smassicen, & + n_smassliq, n_smassliqn, & + n_rhos_cmp, n_rhos_cmpn, & + n_rhos_cnt, n_rhos_cntn, & + n_rsnw, n_rsnwn, & + n_meltsliq, n_fsloss + +!======================================================================= + + contains + +!======================================================================= + + subroutine init_hist_snow_2D (dt) + + use ice_broadcast, only: broadcast_scalar + use ice_calendar, only: nstreams, histfreq + use ice_communicate, only: my_task, master_task + use ice_history_shared, only: tstr2D, tcstr, define_hist_field + use ice_fileunits, only: nu_nml, nml_filename, & + get_fileunit, release_fileunit + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + integer (kind=int_kind) :: ns + integer (kind=int_kind) :: nml_error ! namelist i/o error flag + real (kind=dbl_kind) :: rhofresh, secday + logical (kind=log_kind) :: tr_snow + character(len=*), parameter :: subname = '(init_hist_snow_2D)' + + call icepack_query_tracer_flags(tr_snow_out=tr_snow) + call icepack_query_parameters(rhofresh_out=rhofresh,secday_out=secday) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (tr_snow) then + + !----------------------------------------------------------------- + ! read namelist + !----------------------------------------------------------------- + + call get_fileunit(nu_nml) + if (my_task == master_task) then + open (nu_nml, file=nml_filename, status='old',iostat=nml_error) + if (nml_error /= 0) then + nml_error = -1 + else + nml_error = 1 + endif + do while (nml_error > 0) + read(nu_nml, nml=icefields_snow_nml,iostat=nml_error) + end do + if (nml_error == 0) close(nu_nml) + endif + call release_fileunit(nu_nml) + + call broadcast_scalar(nml_error, master_task) + if (nml_error /= 0) then + close (nu_nml) + call abort_ice('ice: error reading icefields_snow_nml') + endif + + else ! .not. tr_snow + f_smassice = 'x' + f_smassliq = 'x' + f_rhos_cmp = 'x' + f_rhos_cnt = 'x' + f_rsnw = 'x' + f_smassicen= 'x' + f_smassliqn= 'x' + f_rhos_cmpn= 'x' + f_rhos_cntn= 'x' + f_rsnwn = 'x' + f_meltsliq = 'x' + f_fsloss = 'x' + endif + + call broadcast_scalar (f_smassice, master_task) + call broadcast_scalar (f_smassliq, master_task) + call broadcast_scalar (f_rhos_cmp, master_task) + call broadcast_scalar (f_rhos_cnt, master_task) + call broadcast_scalar (f_rsnw, master_task) + call broadcast_scalar (f_smassicen,master_task) + call broadcast_scalar (f_smassliqn,master_task) + call broadcast_scalar (f_rhos_cmpn,master_task) + call broadcast_scalar (f_rhos_cntn,master_task) + call broadcast_scalar (f_rsnwn, master_task) + call broadcast_scalar (f_meltsliq, master_task) + call broadcast_scalar (f_fsloss, master_task) + + if (tr_snow) then + + ! 2D variables + do ns = 1, nstreams + if (histfreq(ns) /= 'x') then + + if (f_smassice(1:1) /= 'x') & + call define_hist_field(n_smassice,"smassice","kg/m^2",tstr2D, tcstr, & + "ice mass per unit area in snow", & + "none", c1, c0, & + ns, f_smassice) + + if (f_smassliq(1:1) /= 'x') & + call define_hist_field(n_smassliq,"smassliq","kg/m^2",tstr2D, tcstr, & + "liquid mass per unit area in snow", & + "none", c1, c0, & + ns, f_smassliq) + + if (f_rhos_cmp(1:1) /= 'x') & + call define_hist_field(n_rhos_cmp,"rhos_cmp","kg/m^3",tstr2D, tcstr, & + "snow density: compaction", & + "none", c1, c0, & + ns, f_rhos_cmp) + + if (f_rhos_cnt(1:1) /= 'x') & + call define_hist_field(n_rhos_cnt,"rhos_cnt","kg/m^3",tstr2D, tcstr, & + "snow density: content", & + "none", c1, c0, & + ns, f_rhos_cnt) + + if (f_rsnw(1:1) /= 'x') & + call define_hist_field(n_rsnw,"rsnw","10^-6 m",tstr2D, tcstr, & + "average snow grain radius", & + "none", c1, c0, & + ns, f_rsnw) + + if (f_meltsliq(1:1) /= 'x') & + call define_hist_field(n_meltsliq,"meltsliq","kg/m^2/s",tstr2D, tcstr, & + "snow liquid contribution to meltponds", & + "none", c1/dt, c0, & + ns, f_meltsliq) + + if (f_fsloss(1:1) /= 'x') & + call define_hist_field(n_fsloss,"fsloss","kg/m^2/s",tstr2D, tcstr, & + "rate of snow loss to leads (liquid)", & + "none", c1, c0, & + ns, f_fsloss) + + endif ! histfreq(ns) /= 'x' + enddo ! nstreams + endif ! tr_snow + + end subroutine init_hist_snow_2D + +!======================================================================= + + subroutine init_hist_snow_3Dc + + use ice_calendar, only: nstreams, histfreq + use ice_history_shared, only: tstr3Dc, tcstr, define_hist_field + + integer (kind=int_kind) :: ns + logical (kind=log_kind) :: tr_snow + character(len=*), parameter :: subname = '(init_hist_pond_3Dc)' + + call icepack_query_tracer_flags(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__) + + if (tr_snow) then + + ! 3D (category) variables must be looped separately + do ns = 1, nstreams + if (histfreq(ns) /= 'x') then + + if (f_smassicen(1:1) /= 'x') & + call define_hist_field(n_smassicen,"smassicen","kg/m^2",tstr3Dc, tcstr, & + "ice mass per unit area in snow, category", & + "none", c1, c0, & + ns, f_smassicen) + + if (f_smassliqn(1:1) /= 'x') & + call define_hist_field(n_smassliqn,"smassliqn","kg/m^2",tstr3Dc, tcstr, & + "liquid mass per unit area in snow, category", & + "none", c1, c0, & + ns, f_smassliqn) + + if (f_rhos_cmpn(1:1) /= 'x') & + call define_hist_field(n_rhos_cmpn,"rhos_cmpn","kg/m^3",tstr3Dc, tcstr, & + "snow density: compaction, category", & + "none", c1, c0, & + ns, f_rhos_cmpn) + + if (f_rhos_cntn(1:1) /= 'x') & + call define_hist_field(n_rhos_cntn,"rhos_cntn","kg/m^3",tstr3Dc, tcstr, & + "snow density: content, category", & + "none", c1, c0, & + ns, f_rhos_cntn) + + if (f_rsnwn(1:1) /= 'x') & + call define_hist_field(n_rsnwn,"rsnwn","10^-6 m",tstr3Dc, tcstr, & + "average snow grain radius, category", & + "none", c1, c0, & + ns, f_rsnwn) + + endif ! histfreq(ns) /= 'x' + enddo ! ns + + endif ! tr_snow + + end subroutine init_hist_snow_3Dc + +!======================================================================= + +! accumulate average ice quantities or snapshots + + subroutine accum_hist_snow (iblk) + + use ice_arrays_column, only: meltsliq + use ice_blocks, only: block, nx_block, ny_block + use ice_domain, only: blocks_ice + use ice_flux, only: fsloss + use ice_history_shared, only: n2D, a2D, a3Dc, ncat_hist, & + accum_hist_field, nzslyr + use ice_state, only: vsno, vsnon, trcr, trcrn + + integer (kind=int_kind), intent(in) :: & + iblk ! block index + + ! local variables + + integer (kind=int_kind) :: & + i, j, k, n + + integer (kind=int_kind) :: & + nt_smice, nt_smliq, nt_rhos, nt_rsnw + + logical (kind=log_kind) :: tr_snow + + real (kind=dbl_kind), dimension (nx_block,ny_block) :: & + worka + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat_hist) :: & + workb + + character(len=*), parameter :: subname = '(accum_hist_snow)' + + !--------------------------------------------------------------- + ! increment field + !--------------------------------------------------------------- + + call icepack_query_tracer_flags(tr_snow_out=tr_snow) + call icepack_query_tracer_indices(nt_smice_out=nt_smice, & + nt_smliq_out=nt_smliq, nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (allocated(a2D)) then + if (tr_snow) then + + if (f_smassice(1:1)/= 'x') then + worka(:,:) = c0 + do k = 1, nzslyr + worka(:,:) = worka(:,:) & + + trcr(:,:,nt_smice+k-1,iblk) + enddo + worka(:,:) = worka(:,:) * vsno(:,:,iblk) / real(nslyr,kind=dbl_kind) + call accum_hist_field(n_smassice, iblk, worka, a2D) + endif + if (f_smassliq(1:1)/= 'x') then + worka(:,:) = c0 + do k = 1, nzslyr + worka(:,:) = worka(:,:) & + + trcr(:,:,nt_smliq+k-1,iblk) + enddo + worka(:,:) = worka(:,:) * vsno(:,:,iblk) / real(nslyr,kind=dbl_kind) + call accum_hist_field(n_smassliq, iblk, worka, a2D) + endif + if (f_rhos_cmp(1:1)/= 'x') then + worka(:,:) = c0 + do k = 1, nzslyr + worka(:,:) = worka(:,:) & + + trcr(:,:,nt_rhos+k-1,iblk) + enddo + worka(:,:) = worka(:,:) / real(nslyr,kind=dbl_kind) + call accum_hist_field(n_rhos_cmp, iblk, worka, a2D) + endif + if (f_rhos_cnt(1:1)/= 'x') then + worka(:,:) = c0 + do k = 1, nzslyr + worka(:,:) = worka(:,:) & + + trcr(:,:,nt_smice+k-1,iblk) & + + trcr(:,:,nt_smliq+k-1,iblk) + enddo + worka(:,:) = worka(:,:) / real(nslyr,kind=dbl_kind) + call accum_hist_field(n_rhos_cnt, iblk, worka, a2D) + endif + if (f_rsnw(1:1)/= 'x') then + worka(:,:) = c0 + do k = 1, nzslyr + worka(:,:) = worka(:,:) & + + trcr(:,:,nt_rsnw+k-1,iblk) + enddo + worka(:,:) = worka(:,:) / real(nslyr,kind=dbl_kind) + call accum_hist_field(n_rsnw, iblk, worka, a2D) + endif + if (f_meltsliq(1:1)/= 'x') & + call accum_hist_field(n_meltsliq, iblk, & + meltsliq(:,:,iblk), a2D) + if (f_fsloss(1:1)/= 'x') & + call accum_hist_field(n_fsloss, iblk, & + fsloss(:,:,iblk), a2D) + + endif ! allocated(a2D) + + ! 3D category fields + if (allocated(a3Dc)) then + if (f_smassicen(1:1)/= 'x') then + workb(:,:,:) = c0 + do n = 1, ncat_hist + do k = 1, nzslyr + workb(:,:,n) = workb(:,:,n) & + + trcrn(:,:,nt_smice+k-1,n,iblk) + enddo + workb(:,:,n) = workb(:,:,n) & + * vsnon(:,:,n,iblk) / real(nslyr,kind=dbl_kind) + enddo + call accum_hist_field(n_smassicen-n2D, iblk, ncat_hist, workb, a3Dc) + endif + if (f_smassliqn(1:1)/= 'x') then + workb(:,:,:) = c0 + do n = 1, ncat_hist + do k = 1, nzslyr + workb(:,:,n) = workb(:,:,n) & + + trcrn(:,:,nt_smliq+k-1,n,iblk) + enddo + workb(:,:,n) = workb(:,:,n) & + * vsnon(:,:,n,iblk) / real(nslyr,kind=dbl_kind) + enddo + call accum_hist_field(n_smassliqn-n2D, iblk, ncat_hist, workb, a3Dc) + endif + if (f_rhos_cmpn(1:1)/= 'x') then + workb(:,:,:) = c0 + do n = 1, ncat_hist + do k = 1, nzslyr + workb(:,:,n) = workb(:,:,n) & + + trcrn(:,:,nt_rhos+k-1,n,iblk) + enddo + workb(:,:,n) = workb(:,:,n) / real(nslyr,kind=dbl_kind) + enddo + call accum_hist_field(n_rhos_cmpn-n2D, iblk, ncat_hist, workb, a3Dc) + endif + if (f_rhos_cntn(1:1)/= 'x') then + workb(:,:,:) = c0 + do n = 1, ncat_hist + do k = 1, nzslyr + workb(:,:,n) = workb(:,:,n) & + + trcrn(:,:,nt_smice+k-1,n,iblk) & + + trcrn(:,:,nt_smliq+k-1,n,iblk) + enddo + workb(:,:,n) = workb(:,:,n) / real(nslyr,kind=dbl_kind) + enddo + call accum_hist_field(n_rhos_cntn-n2D, iblk, ncat_hist, workb, a3Dc) + endif + if (f_rsnwn(1:1)/= 'x') then + workb(:,:,:) = c0 + do n = 1, ncat_hist + do k = 1, nzslyr + workb(:,:,n) = workb(:,:,n) & + + trcrn(:,:,nt_rsnw+k-1,n,iblk) + enddo + workb(:,:,n) = workb(:,:,n) / real(nslyr,kind=dbl_kind) + enddo + call accum_hist_field(n_rsnwn-n2D, iblk, ncat_hist, workb, a3Dc) + endif + endif ! allocated(a3Dc) + + endif ! tr_snow + + end subroutine accum_hist_snow + +!======================================================================= + + end module ice_history_snow + +!======================================================================= diff --git a/cicecore/cicedynB/dynamics/ice_transport_driver.F90 b/cicecore/cicedynB/dynamics/ice_transport_driver.F90 index e3da6390b..f2dff2367 100644 --- a/cicecore/cicedynB/dynamics/ice_transport_driver.F90 +++ b/cicecore/cicedynB/dynamics/ice_transport_driver.F90 @@ -85,7 +85,9 @@ subroutine init_transport integer (kind=int_kind) :: ntrcr, nt_Tsfc, nt_qice, nt_qsno, & nt_sice, nt_fbri, nt_iage, nt_FY, nt_alvl, nt_vlvl, & - nt_apnd, nt_hpnd, nt_ipnd, nt_fsd, nt_isosno, nt_isoice, nt_bgc_Nit, nt_bgc_S + nt_apnd, nt_hpnd, nt_ipnd, nt_fsd, & + nt_smice, nt_smliq, nt_rhos, nt_rsnw, & + nt_isosno, nt_isoice, nt_bgc_Nit, nt_bgc_S character(len=*), parameter :: subname = '(init_transport)' @@ -94,9 +96,12 @@ subroutine init_transport call icepack_query_tracer_sizes(ntrcr_out=ntrcr) 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_fbri_out=nt_fbri, & - nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_alvl_out=nt_alvl, nt_fsd_out=nt_fsd, & - nt_vlvl_out=nt_vlvl, nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, & - nt_ipnd_out=nt_ipnd, nt_bgc_Nit_out=nt_bgc_Nit, nt_bgc_S_out=nt_bgc_S, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_fsd_out=nt_fsd, & + 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_smice_out=nt_smice, nt_smliq_out=nt_smliq, nt_rhos_out=nt_rhos, & + nt_rsnw_out=nt_rsnw, & + nt_bgc_Nit_out=nt_bgc_Nit, nt_bgc_S_out=nt_bgc_S, & 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, & @@ -195,6 +200,18 @@ subroutine init_transport if (nt-k==nt_ipnd) & write(nu_diag,1000) 'nt_ipnd ',nt,depend(nt),tracer_type(nt),& has_dependents(nt) + if (nt-k==nt_smice) & + write(nu_diag,1000) 'nt_smice ',nt,depend(nt),tracer_type(nt),& + has_dependents(nt) + if (nt-k==nt_smliq) & + write(nu_diag,1000) 'nt_smliq ',nt,depend(nt),tracer_type(nt),& + has_dependents(nt) + if (nt-k==nt_rhos) & + write(nu_diag,1000) 'nt_rhos ',nt,depend(nt),tracer_type(nt),& + has_dependents(nt) + if (nt-k==nt_rsnw) & + write(nu_diag,1000) 'nt_rsnw ',nt,depend(nt),tracer_type(nt),& + has_dependents(nt) if (nt-k==nt_fsd) & write(nu_diag,1000) 'nt_fsd ',nt,depend(nt),tracer_type(nt),& has_dependents(nt) diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index bcc7305ff..23fb9df63 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -218,6 +218,7 @@ module ice_flux fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) + fsloss , & ! rate of snow loss to leads (kg/m^2/s) fswthru , & ! shortwave penetrating to ocean (W/m^2) fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2) @@ -294,6 +295,10 @@ module ice_flux fsensn, & ! category sensible heat flux flatn ! category latent heat flux + real (kind=dbl_kind), & + dimension (:,:,:,:), allocatable, public :: & + snwcnt ! counter for presence of snow + ! As above but these remain grid box mean values i.e. they are not ! divided by aice at end of ice_dynamics. These are used in ! CICE_IN_NEMO for coupling and also for generating @@ -448,6 +453,7 @@ subroutine alloc_flux fresh (nx_block,ny_block,max_blocks), & ! fresh water flux to ocean (kg/m^2/s) fsalt (nx_block,ny_block,max_blocks), & ! salt flux to ocean (kg/m^2/s) fhocn (nx_block,ny_block,max_blocks), & ! net heat flux to ocean (W/m^2) + fsloss (nx_block,ny_block,max_blocks), & ! rate of snow loss to leads (kg/m^2/s) fswthru (nx_block,ny_block,max_blocks), & ! shortwave penetrating to ocean (W/m^2) fswthru_vdr (nx_block,ny_block,max_blocks), & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_vdf (nx_block,ny_block,max_blocks), & ! vis dif shortwave penetrating to ocean (W/m^2) @@ -525,6 +531,7 @@ subroutine alloc_flux fsensn (nx_block,ny_block,ncat,max_blocks), & ! category sensible heat flux flatn (nx_block,ny_block,ncat,max_blocks), & ! category latent heat flux albcnt (nx_block,ny_block,max_blocks,max_nstrm), & ! counter for zenith angle + snwcnt (nx_block,ny_block,max_blocks,max_nstrm), & ! counter for snow salinz (nx_block,ny_block,nilyr+1,max_blocks), & ! initial salinity profile (ppt) Tmltz (nx_block,ny_block,nilyr+1,max_blocks), & ! initial melting temperature (^oC) stat=ierr) diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index a71e6dd17..84bf1d461 100755 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -41,7 +41,7 @@ module ice_forcing field_type_vector, field_loc_NEcorner use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_sea_freezing_temperature - use icepack_intfc, only: icepack_init_wave + use icepack_intfc, only: icepack_init_wave, icepack_init_parameters use icepack_intfc, only: icepack_query_tracer_indices, icepack_query_parameters implicit none @@ -50,7 +50,8 @@ module ice_forcing get_forcing_atmo, get_forcing_ocn, get_wave_spec, & read_clim_data, read_clim_data_nc, & interpolate_data, interp_coeff_monthly, & - read_data_nc_point, interp_coeff + read_data_nc_point, interp_coeff, & + init_snowtable integer (kind=int_kind), public :: & ycycle , & ! number of years in forcing cycle, set by namelist @@ -166,6 +167,16 @@ module ice_forcing integer (kind=int_kind), public :: & Njday_atm ! Number of atm forcing timesteps + character (len=char_len_long), public :: & + snw_filename ! filename for snow lookup table + + character (char_len), public :: & + snw_rhos_fname , & ! snow table 1d rhos field name + snw_Tgrd_fname , & ! snow table 1d Tgrd field name + snw_T_fname , & ! snow table 1d T field name + snw_tau_fname , & ! snow table 3d tau field name + snw_kappa_fname, & ! snow table 3d kappa field name + snw_drdt0_fname ! snow table 3d drdt0 field name ! PRIVATE: @@ -5398,7 +5409,199 @@ end subroutine get_wave_spec !======================================================================= - end module ice_forcing +! initial snow aging lookup table +! +! Dry snow metamorphism table +! snicar_drdt_bst_fit_60_c070416.nc +! Flanner (file metadata units mislabelled) +! drdsdt0 (10^-6 m/hr) tau (10^-6 m) +! + subroutine init_snowtable + + use ice_broadcast, only: broadcast_array, broadcast_scalar + integer (kind=int_kind) :: & + idx_T_max , & ! Table dimensions + idx_rhos_max, & + idx_Tgrd_max + real (kind=dbl_kind), allocatable :: & + snowage_rhos (:), & + snowage_Tgrd (:), & + snowage_T (:), & + snowage_tau (:,:,:), & + snowage_kappa(:,:,:), & + snowage_drdt0(:,:,:) + + ! local variables + + logical (kind=log_kind) :: diag = .false. + + integer (kind=int_kind) :: & + fid ! file id for netCDF file + + character (char_len) :: & + snw_aging_table, & ! aging table setting + fieldname ! field name in netcdf file + + integer (kind=int_kind) :: & + j, k ! indices + + character(len=*), parameter :: subname = '(init_snowtable)' + + !----------------------------------------------------------------- + ! read table of snow aging parameters + !----------------------------------------------------------------- + + call icepack_query_parameters(snw_aging_table_out=snw_aging_table, & + isnw_rhos_out=idx_rhos_max, isnw_Tgrd_out=idx_Tgrd_max, isnw_T_out=idx_T_max) + + if (my_task == master_task) then + write (nu_diag,*) ' ' + write (nu_diag,*) 'Snow aging file:', trim(snw_filename) + endif + + if (snw_aging_table == 'snicar') then + ! just read the 3d data and pass it in + + call ice_open_nc(snw_filename,fid) + + allocate(snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max)) + allocate(snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max)) + allocate(snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max)) + + fieldname = trim(snw_tau_fname) + call ice_read_nc(fid,fieldname,snowage_tau, diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + fieldname = trim(snw_kappa_fname) + call ice_read_nc(fid,fieldname,snowage_kappa,diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + fieldname = trim(snw_drdt0_fname) + call ice_read_nc(fid,fieldname,snowage_drdt0,diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + + call ice_close_nc(fid) + + call broadcast_array(snowage_tau , master_task) + call broadcast_array(snowage_kappa, master_task) + call broadcast_array(snowage_drdt0, master_task) + + if (my_task == master_task) then + write(nu_diag,*) subname,' ' + write(nu_diag,*) subname,' Successfully read snow aging properties:' + write(nu_diag,*) subname,' snw_aging_table = ',trim(snw_aging_table) + write(nu_diag,*) subname,' idx_rhos_max = ',idx_rhos_max + write(nu_diag,*) subname,' idx_Tgrd_max = ',idx_Tgrd_max + write(nu_diag,*) subname,' idx_T_max = ',idx_T_max + write(nu_diag,*) subname,' Data at rhos, Tgrd, T at first index ' + write(nu_diag,*) subname,' snoage_tau (1,1,1) = ',snowage_tau (1,1,1) + write(nu_diag,*) subname,' snoage_kappa (1,1,1) = ',snowage_kappa(1,1,1) + write(nu_diag,*) subname,' snoage_drdt0 (1,1,1) = ',snowage_drdt0(1,1,1) + write(nu_diag,*) subname,' Data at rhos, Tgrd, T at max index' + write(nu_diag,*) subname,' snoage_tau (max,max,max) = ',snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max) + write(nu_diag,*) subname,' snoage_kappa (max,max,max) = ',snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max) + write(nu_diag,*) subname,' snoage_drdt0 (max,max,max) = ',snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max) + endif + + call icepack_init_parameters( & + snowage_tau_in = snowage_tau, & + snowage_kappa_in = snowage_kappa, & + snowage_drdt0_in = snowage_drdt0 ) + + deallocate(snowage_tau) + deallocate(snowage_kappa) + deallocate(snowage_drdt0) + + else + ! read everything and pass it in + + call ice_open_nc(snw_filename,fid) + + fieldname = trim(snw_rhos_fname) + call ice_get_ncvarsize(fid,fieldname,idx_rhos_max) + fieldname = trim(snw_Tgrd_fname) + call ice_get_ncvarsize(fid,fieldname,idx_Tgrd_max) + fieldname = trim(snw_T_fname) + call ice_get_ncvarsize(fid,fieldname,idx_T_max) + + call broadcast_scalar(idx_rhos_max, master_task) + call broadcast_scalar(idx_Tgrd_max, master_task) + call broadcast_scalar(idx_T_max , master_task) + + allocate(snowage_rhos (idx_rhos_max)) + allocate(snowage_Tgrd (idx_Tgrd_max)) + allocate(snowage_T (idx_T_max)) + allocate(snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max)) + allocate(snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max)) + allocate(snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max)) + + fieldname = trim(snw_rhos_fname) + call ice_read_nc(fid,fieldname,snowage_rhos, diag, & + idx_rhos_max) + fieldname = trim(snw_Tgrd_fname) + call ice_read_nc(fid,fieldname,snowage_Tgrd, diag, & + idx_Tgrd_max) + fieldname = trim(snw_T_fname) + call ice_read_nc(fid,fieldname,snowage_T, diag, & + idx_T_max) + + fieldname = trim(snw_tau_fname) + call ice_read_nc(fid,fieldname,snowage_tau, diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + fieldname = trim(snw_kappa_fname) + call ice_read_nc(fid,fieldname,snowage_kappa,diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + fieldname = trim(snw_drdt0_fname) + call ice_read_nc(fid,fieldname,snowage_drdt0,diag, & + idx_rhos_max,idx_Tgrd_max,idx_T_max) + + call ice_close_nc(fid) + + call broadcast_array(snowage_rhos , master_task) + call broadcast_array(snowage_Tgrd , master_task) + call broadcast_array(snowage_T , master_task) + call broadcast_array(snowage_tau , master_task) + call broadcast_array(snowage_kappa, master_task) + call broadcast_array(snowage_drdt0, master_task) + + if (my_task == master_task) then + write(nu_diag,*) subname,' ' + write(nu_diag,*) subname,' Successfully read snow aging properties:' + write(nu_diag,*) subname,' idx_rhos_max = ',idx_rhos_max + write(nu_diag,*) subname,' idx_Tgrd_max = ',idx_Tgrd_max + write(nu_diag,*) subname,' idx_T_max = ',idx_T_max + write(nu_diag,*) subname,' Data at rhos, Tgrd, T = ',snowage_rhos(1),snowage_Tgrd(1),snowage_T(1) + write(nu_diag,*) subname,' snoage_tau (1,1,1) = ',snowage_tau (1,1,1) + write(nu_diag,*) subname,' snoage_kappa (1,1,1) = ',snowage_kappa(1,1,1) + write(nu_diag,*) subname,' snoage_drdt0 (1,1,1) = ',snowage_drdt0(1,1,1) + write(nu_diag,*) subname,' Data at rhos, Tgrd, T = ',snowage_rhos(idx_rhos_max),snowage_Tgrd(idx_Tgrd_max),snowage_T(idx_T_max) + write(nu_diag,*) subname,' snoage_tau (max,max,max) = ',snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max) + write(nu_diag,*) subname,' snoage_kappa (max,max,max) = ',snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max) + write(nu_diag,*) subname,' snoage_drdt0 (max,max,max) = ',snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max) + endif + + call icepack_init_parameters( & + isnw_t_in = idx_T_max, & + isnw_Tgrd_in = idx_Tgrd_max, & + isnw_rhos_in = idx_rhos_max, & + snowage_rhos_in = snowage_rhos, & + snowage_Tgrd_in = snowage_Tgrd, & + snowage_T_in = snowage_T, & + snowage_tau_in = snowage_tau, & + snowage_kappa_in = snowage_kappa, & + snowage_drdt0_in = snowage_drdt0 ) + + deallocate(snowage_rhos) + deallocate(snowage_Tgrd) + deallocate(snowage_T) + deallocate(snowage_tau) + deallocate(snowage_kappa) + deallocate(snowage_drdt0) + + endif + + end subroutine init_snowtable !======================================================================= + end module ice_forcing + +!======================================================================= diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index c0c089ac8..a96b30d98 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -74,7 +74,7 @@ subroutine input_data use ice_arrays_column, only: oceanmixed_ice use ice_restart_column, only: restart_age, restart_FY, restart_lvl, & restart_pond_cesm, restart_pond_lvl, restart_pond_topo, restart_aero, & - restart_fsd, restart_iso + restart_fsd, restart_iso, restart_snow use ice_restart_shared, only: & restart, restart_ext, restart_coszen, restart_dir, restart_file, pointer_file, & runid, runtype, use_restart_time, restart_format, lcdf64 @@ -91,7 +91,10 @@ subroutine input_data bgc_data_type, & ocn_data_type, ocn_data_dir, wave_spec_file, & oceanmixed_file, restore_ocn, trestore, & - ice_data_type + ice_data_type, & + snw_filename, & + snw_tau_fname, snw_kappa_fname, snw_drdt0_fname, & + snw_rhos_fname, snw_Tgrd_fname, snw_T_fname use ice_arrays_column, only: bgc_data_dir, fe_data_type use ice_grid, only: grid_file, gridcpl_file, kmt_file, & bathymetry_file, use_bathymetry, & @@ -129,19 +132,21 @@ subroutine input_data mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, & a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, & phi_c_slow_mode, phi_i_mushy, kalg, atmiter_conv, Pstar, Cstar, & - sw_frac, sw_dtemp, floediam, hfrazilmin, iceruf, iceruf_ocn + sw_frac, sw_dtemp, floediam, hfrazilmin, iceruf, iceruf_ocn, & + rsnw_fall, rsnw_tmax, rhosnew, rhosmin, rhosmax, & + windmin, drhosdwind, snwlvlfac integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, natmiter, & kitd, kcatbound, ktransport character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & - tfrz_option, frzpnd, atmbndy, wave_spec_type + tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, & - sw_redist, calc_dragio + sw_redist, calc_dragio, use_smliq_pnd, snwgrain logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond - logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd, tr_snow logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo integer (kind=int_kind) :: numin, numax ! unit number limits @@ -188,6 +193,7 @@ subroutine input_data tr_pond_cesm, restart_pond_cesm, & tr_pond_lvl, restart_pond_lvl, & tr_pond_topo, restart_pond_topo, & + tr_snow, restart_snow, & tr_iso, restart_iso, & tr_aero, restart_aero, & tr_fsd, restart_fsd, & @@ -228,6 +234,13 @@ subroutine input_data rfracmin, rfracmax, pndaspect, hs1, & hp1 + namelist /snow_nml/ & + snwredist, snwgrain, rsnw_fall, rsnw_tmax, & + rhosnew, rhosmin, rhosmax, snwlvlfac, & + windmin, drhosdwind, use_smliq_pnd, snw_aging_table,& + snw_filename, snw_rhos_fname, snw_Tgrd_fname,snw_T_fname, & + snw_tau_fname, snw_kappa_fname, snw_drdt0_fname + namelist /forcing_nml/ & formdrag, atmbndy, calc_strair, calc_Tsfc, & highfreq, natmiter, atmiter_conv, calc_dragio, & @@ -411,6 +424,25 @@ subroutine input_data rfracmin = 0.15_dbl_kind ! minimum retained fraction of meltwater rfracmax = 0.85_dbl_kind ! maximum retained fraction of meltwater pndaspect = 0.8_dbl_kind ! ratio of pond depth to area fraction + snwredist = 'none' ! type of snow redistribution + snw_aging_table = 'test' ! snow aging lookup table + snw_filename = 'unknown' ! snowtable filename + snw_tau_fname = 'unknown' ! snowtable file tau fieldname + snw_kappa_fname = 'unknown' ! snowtable file kappa fieldname + snw_drdt0_fname = 'unknown' ! snowtable file drdt0 fieldname + snw_rhos_fname = 'unknown' ! snowtable file rhos fieldname + snw_Tgrd_fname = 'unknown' ! snowtable file Tgrd fieldname + snw_T_fname = 'unknown' ! snowtable file T fieldname + snwgrain = .false. ! snow metamorphosis + use_smliq_pnd = .false. ! use liquid in snow for ponds + rsnw_fall = 100.0_dbl_kind ! radius of new snow (10^-6 m) ! advanced snow physics: 54.526 x 10^-6 m + rsnw_tmax = 1500.0_dbl_kind ! maximum snow radius (10^-6 m) + rhosnew = 100.0_dbl_kind ! new snow density (kg/m^3) + rhosmin = 100.0_dbl_kind ! minimum snow density (kg/m^3) + rhosmax = 450.0_dbl_kind ! maximum snow density (kg/m^3) + windmin = 10.0_dbl_kind ! minimum wind speed to compact snow (m/s) + drhosdwind= 27.3_dbl_kind ! wind compaction factor for snow (kg s/m^4) + snwlvlfac = 0.3_dbl_kind ! fractional increase in snow depth for bulk redistribution albicev = 0.78_dbl_kind ! visible ice albedo for h > ahmax albicei = 0.36_dbl_kind ! near-ir ice albedo for h > ahmax albsnowv = 0.98_dbl_kind ! cold snow albedo, visible @@ -474,6 +506,8 @@ subroutine input_data restart_pond_lvl = .false. ! melt ponds restart tr_pond_topo = .false. ! explicit melt ponds (topographic) restart_pond_topo = .false. ! melt ponds restart + tr_snow = .false. ! advanced snow physics + restart_snow = .false. ! advanced snow physics restart tr_iso = .false. ! isotopes restart_iso = .false. ! isotopes restart tr_aero = .false. ! aerosols @@ -547,6 +581,9 @@ subroutine input_data print*,'Reading ponds_nml' read(nu_nml, nml=ponds_nml,iostat=nml_error) if (nml_error /= 0) exit + print*,'Reading snow_nml' + read(nu_nml, nml=snow_nml,iostat=nml_error) + if (nml_error /= 0) exit print*,'Reading forcing_nml' read(nu_nml, nml=forcing_nml,iostat=nml_error) if (nml_error /= 0) exit @@ -737,6 +774,25 @@ subroutine input_data call broadcast_scalar(rfracmin, master_task) call broadcast_scalar(rfracmax, master_task) call broadcast_scalar(pndaspect, master_task) + call broadcast_scalar(snwredist, master_task) + call broadcast_scalar(snw_aging_table, master_task) + call broadcast_scalar(snw_filename, master_task) + call broadcast_scalar(snw_tau_fname, master_task) + call broadcast_scalar(snw_kappa_fname, master_task) + call broadcast_scalar(snw_drdt0_fname, master_task) + call broadcast_scalar(snw_rhos_fname, master_task) + call broadcast_scalar(snw_Tgrd_fname, master_task) + call broadcast_scalar(snw_T_fname, master_task) + call broadcast_scalar(snwgrain, master_task) + call broadcast_scalar(use_smliq_pnd, master_task) + call broadcast_scalar(rsnw_fall, master_task) + call broadcast_scalar(rsnw_tmax, master_task) + call broadcast_scalar(rhosnew, master_task) + call broadcast_scalar(rhosmin, master_task) + call broadcast_scalar(rhosmax, master_task) + call broadcast_scalar(windmin, master_task) + call broadcast_scalar(drhosdwind, master_task) + call broadcast_scalar(snwlvlfac, master_task) call broadcast_scalar(albicev, master_task) call broadcast_scalar(albicei, master_task) call broadcast_scalar(albsnowv, master_task) @@ -800,6 +856,8 @@ subroutine input_data call broadcast_scalar(restart_pond_lvl, master_task) call broadcast_scalar(tr_pond_topo, master_task) call broadcast_scalar(restart_pond_topo, master_task) + call broadcast_scalar(tr_snow, master_task) + call broadcast_scalar(restart_snow, master_task) call broadcast_scalar(tr_iso, master_task) call broadcast_scalar(restart_iso, master_task) call broadcast_scalar(tr_aero, master_task) @@ -880,6 +938,7 @@ subroutine input_data restart_pond_cesm = .false. restart_pond_lvl = .false. restart_pond_topo = .false. + restart_snow = .false. ! tcraig, OK to leave as true, needed for boxrestore case ! restart_ext = .false. else @@ -988,6 +1047,59 @@ subroutine input_data abort_list = trim(abort_list)//":8" endif + if (snwredist(1:4) /= 'none' .and. .not. tr_snow) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: snwredist on but tr_snow=F' + write (nu_diag,*) 'ERROR: Use tr_snow=T for snow redistribution' + endif + abort_list = trim(abort_list)//":37" + endif + if (snwredist(1:4) == 'bulk' .and. .not. tr_lvl) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: snwredist=bulk but tr_lvl=F' + write (nu_diag,*) 'ERROR: Use tr_lvl=T for snow redistribution' + endif + abort_list = trim(abort_list)//":38" + endif + if (snwredist(1:6) == 'ITDrdg' .and. .not. tr_lvl) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: snwredist=ITDrdg but tr_lvl=F' + write (nu_diag,*) 'ERROR: Use tr_lvl=T for snow redistribution' + endif + abort_list = trim(abort_list)//":39" + endif + if (use_smliq_pnd .and. .not. snwgrain) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: use_smliq_pnd = T but' + write (nu_diag,*) 'ERROR: snow metamorphosis not used' + write (nu_diag,*) 'ERROR: Use snwgrain=T with smliq for ponds' + endif + abort_list = trim(abort_list)//":40" + endif + if (use_smliq_pnd .and. .not. tr_snow) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: use_smliq_pnd = T but' + write (nu_diag,*) 'ERROR: snow tracers are not active' + write (nu_diag,*) 'ERROR: Use tr_snow=T with smliq for ponds' + endif + abort_list = trim(abort_list)//":41" + endif + if (snwgrain .and. .not. tr_snow) then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: snwgrain=T but tr_snow=F' + write (nu_diag,*) 'ERROR: Use tr_snow=T for snow metamorphosis' + endif + abort_list = trim(abort_list)//":42" + endif + if (trim(snw_aging_table) /= 'test' .and. & + trim(snw_aging_table) /= 'snicar' .and. & + trim(snw_aging_table) /= 'file') then + if (my_task == master_task) then + write (nu_diag,*) 'ERROR: unknown snw_aging_table = '//trim(snw_aging_table) + endif + abort_list = trim(abort_list)//":43" + endif + if (tr_iso .and. n_iso==0) then if (my_task == master_task) then write(nu_diag,*) subname//' ERROR: isotopes activated but' @@ -1017,7 +1129,7 @@ subroutine input_data if (my_task == master_task) then write(nu_diag,*) subname//' ERROR: nilyr < 1' endif - abort_list = trim(abort_list)//":33" + abort_list = trim(abort_list)//":2" endif if (nslyr < 1) then @@ -1051,6 +1163,13 @@ subroutine input_data abort_list = trim(abort_list)//":10" endif + if (trim(shortwave) /= 'dEdd' .and. snwgrain) then + if (my_task == master_task) then + write (nu_diag,*) 'WARNING: snow grain radius activated but' + write (nu_diag,*) 'WARNING: dEdd shortwave is not.' + endif + endif + if ((rfracmin < -puny .or. rfracmin > c1+puny) .or. & (rfracmax < -puny .or. rfracmax > c1+puny) .or. & (rfracmin > rfracmax)) then @@ -1655,6 +1774,78 @@ subroutine input_data write(nu_diag,1002) ' rfracmin = ', rfracmin,' : minimum fraction of melt water added to ponds' write(nu_diag,1002) ' rfracmax = ', rfracmax,' : maximum fraction of melt water added to ponds' + write(nu_diag,*) ' ' + write(nu_diag,*) ' Snow redistribution/metamorphism tracers' + write(nu_diag,*) '-----------------------------------------' + if (tr_snow) then + write(nu_diag,1010) ' tr_snow = ', tr_snow, & + ' : advanced snow physics' + if (snwredist(1:4) == 'none') then + write(nu_diag,1030) ' snwredist = ', trim(snwredist), & + ' : Snow redistribution scheme turned off' + else + if (snwredist(1:4) == 'bulk') then + write(nu_diag,1030) ' snwredist = ', trim(snwredist), & + ' : Using bulk snow redistribution scheme' + elseif (snwredist(1:6) == 'ITDrdg') then + write(nu_diag,1030) ' snwredist = ', trim(snwredist), & + ' : Using ridging based snow redistribution scheme' + write(nu_diag,1002) ' rhosnew = ', rhosnew, & + ' : new snow density (kg/m^3)' + write(nu_diag,1002) ' rhosmin = ', rhosmin, & + ' : minimum snow density (kg/m^3)' + write(nu_diag,1002) ' rhosmax = ', rhosmax, & + ' : maximum snow density (kg/m^3)' + write(nu_diag,1002) ' windmin = ', windmin, & + ' : minimum wind speed to compact snow (m/s)' + write(nu_diag,1002) ' drhosdwind = ', drhosdwind, & + ' : wind compaction factor (kg s/m^4)' + endif + write(nu_diag,1002) ' snwlvlfac = ', snwlvlfac, & + ' : fractional increase in snow depth for redistribution on ridges' + endif + if (.not. snwgrain) then + write(nu_diag,1010) ' snwgrain = ', snwgrain, & + ' : Snow metamorphosis turned off' + else + write(nu_diag,1010) ' snwgrain = ', snwgrain, & + ' : Using snow metamorphosis scheme' + write(nu_diag,1002) ' rsnw_tmax = ', rsnw_tmax, & + ' : maximum snow radius (10^-6 m)' + endif + write(nu_diag,1002) ' rsnw_fall = ', rsnw_fall, & + ' : radius of new snow (10^-6 m)' + if (snwgrain) then + if (use_smliq_pnd) then + tmpstr2 = ' : Using liquid water in snow for melt ponds' + else + tmpstr2 = ' : NOT using liquid water in snow for melt ponds' + endif + write(nu_diag,1010) ' use_smliq_pnd = ', use_smliq_pnd, trim(tmpstr2) + if (snw_aging_table == 'test') then + tmpstr2 = ' : Using 5x5x1 test matrix of internallly defined snow aging parameters' + write(nu_diag,1030) ' snw_aging_table = ', trim(snw_aging_table),trim(tmpstr2) + elseif (snw_aging_table == 'snicar') then + tmpstr2 = ' : Reading 3D snow aging parameters from SNICAR file' + write(nu_diag,1030) ' snw_aging_table = ', trim(snw_aging_table),trim(tmpstr2) + write(nu_diag,1031) ' snw_filename = ',trim(snw_filename) + write(nu_diag,1031) ' snw_tau_fname = ',trim(snw_tau_fname) + write(nu_diag,1031) ' snw_kappa_fname = ',trim(snw_kappa_fname) + write(nu_diag,1031) ' snw_drdt0_fname = ',trim(snw_drdt0_fname) + elseif (snw_aging_table == 'file') then + tmpstr2 = ' : Reading 1D and 3D snow aging dimensions and parameters from external file' + write(nu_diag,1030) ' snw_aging_table = ', trim(snw_aging_table),trim(tmpstr2) + write(nu_diag,1031) ' snw_filename = ',trim(snw_filename) + write(nu_diag,1031) ' snw_rhos_fname = ',trim(snw_rhos_fname) + write(nu_diag,1031) ' snw_Tgrd_fname = ',trim(snw_Tgrd_fname) + write(nu_diag,1031) ' snw_T_fname = ',trim(snw_T_fname) + write(nu_diag,1031) ' snw_tau_fname = ',trim(snw_tau_fname) + write(nu_diag,1031) ' snw_kappa_fname = ',trim(snw_kappa_fname) + write(nu_diag,1031) ' snw_drdt0_fname = ',trim(snw_drdt0_fname) + endif + endif + endif + write(nu_diag,*) ' ' write(nu_diag,*) ' Primary state variables, tracers' write(nu_diag,*) ' (excluding biogeochemistry)' @@ -1668,6 +1859,7 @@ subroutine input_data if (tr_pond_lvl) write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl,' : level-ice pond formulation' if (tr_pond_topo) write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo,' : topo pond formulation' if (tr_pond_cesm) write(nu_diag,1010) ' tr_pond_cesm = ', tr_pond_cesm,' : CESM pond formulation' + if (tr_snow) write(nu_diag,1010) ' tr_snow = ', tr_snow,' : advanced snow physics' if (tr_iage) write(nu_diag,1010) ' tr_iage = ', tr_iage,' : chronological ice age' if (tr_FY) write(nu_diag,1010) ' tr_FY = ', tr_FY,' : first-year ice area' if (tr_iso) write(nu_diag,1010) ' tr_iso = ', tr_iso,' : diagnostic isotope tracers' @@ -1789,6 +1981,7 @@ subroutine input_data write(nu_diag,1011) ' restart_pond_cesm= ', restart_pond_cesm write(nu_diag,1011) ' restart_pond_lvl = ', restart_pond_lvl write(nu_diag,1011) ' restart_pond_topo= ', restart_pond_topo + write(nu_diag,1011) ' restart_snow = ', restart_snow write(nu_diag,1011) ' restart_iso = ', restart_iso write(nu_diag,1011) ' restart_aero = ', restart_aero write(nu_diag,1011) ' restart_fsd = ', restart_fsd @@ -1853,10 +2046,14 @@ subroutine input_data wave_spec_in=wave_spec, nfreq_in=nfreq, & tfrz_option_in=tfrz_option, kalg_in=kalg, fbot_xfer_type_in=fbot_xfer_type, & Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, iceruf_ocn_in=iceruf_ocn, calc_dragio_in=calc_dragio, & + windmin_in=windmin, drhosdwind_in=drhosdwind, & + rsnw_fall_in=rsnw_fall, rsnw_tmax_in=rsnw_tmax, rhosnew_in=rhosnew, & + snwlvlfac_in=snwlvlfac, rhosmin_in=rhosmin, rhosmax_in=rhosmax, & + snwredist_in=snwredist, snwgrain_in=snwgrain, snw_aging_table_in=trim(snw_aging_table), & sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp) call icepack_init_tracer_flags(tr_iage_in=tr_iage, tr_FY_in=tr_FY, & tr_lvl_in=tr_lvl, tr_iso_in=tr_iso, tr_aero_in=tr_aero, & - tr_fsd_in=tr_fsd, tr_pond_in=tr_pond, & + tr_fsd_in=tr_fsd, tr_snow_in=tr_snow, 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_sizes(ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & nfsd_in=nfsd, n_algae_in=n_algae, n_iso_in=n_iso, n_aero_in=n_aero, & @@ -1913,10 +2110,12 @@ subroutine init_state heat_capacity ! from icepack integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo + logical (kind=log_kind) :: tr_snow, tr_fsd integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw integer (kind=int_kind) :: nt_isosno, nt_isoice, nt_aero, nt_fsd type (block) :: & @@ -1929,12 +2128,15 @@ subroutine init_state call icepack_query_parameters(heat_capacity_out=heat_capacity) call icepack_query_tracer_sizes(ntrcr_out=ntrcr) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_fsd_out=tr_fsd, & - tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo) + tr_lvl_out=tr_lvl, tr_iso_out=tr_iso, tr_aero_out=tr_aero, & + tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo, & + tr_snow_out=tr_snow, tr_fsd_out=tr_fsd) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, 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_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_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) @@ -2011,6 +2213,14 @@ subroutine init_state trcr_depend(nt_hpnd) = 2+nt_apnd ! melt pond depth trcr_depend(nt_ipnd) = 2+nt_apnd ! refrozen pond lid endif + if (tr_snow) then ! snow-volume-weighted snow tracers + do k = 1, nslyr + trcr_depend(nt_smice + k - 1) = 2 ! ice mass in snow + trcr_depend(nt_smliq + k - 1) = 2 ! liquid mass in snow + trcr_depend(nt_rhos + k - 1) = 2 ! effective snow density + trcr_depend(nt_rsnw + k - 1) = 2 ! snow radius + enddo + endif if (tr_fsd) then do it = 1, nfsd trcr_depend(nt_fsd + it - 1) = 0 ! area-weighted floe size distribution @@ -2241,7 +2451,7 @@ subroutine set_state_var (nx_block, ny_block, & indxi, indxj ! compressed indices for cells with aicen > puny real (kind=dbl_kind) :: & - Tsfc, sum, hbar, puny, rhos, Lfresh, rad_to_deg + Tsfc, sum, hbar, puny, rhos, Lfresh, rad_to_deg, rsnw_fall real (kind=dbl_kind), dimension(ncat) :: & ainit, hinit ! initial area, thickness @@ -2257,22 +2467,26 @@ subroutine set_state_var (nx_block, ny_block, & edge_init_nh = 70._dbl_kind, & ! initial ice edge, N.Hem. (deg) edge_init_sh = -60._dbl_kind ! initial ice edge, S.Hem. (deg) - logical (kind=log_kind) :: tr_brine, tr_lvl + logical (kind=log_kind) :: tr_brine, tr_lvl, tr_snow integer (kind=int_kind) :: ntrcr integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw character(len=*), parameter :: subname='(set_state_var)' !----------------------------------------------------------------- call icepack_query_tracer_sizes(ntrcr_out=ntrcr) - call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl) + call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl, & + 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_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl) + nt_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & + nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & + nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw) call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny, & - rad_to_deg_out=rad_to_deg) + rad_to_deg_out=rad_to_deg, rsnw_fall_out=rsnw_fall) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -2304,6 +2518,14 @@ subroutine set_state_var (nx_block, ny_block, & do k = 1, nslyr trcrn(i,j,nt_qsno+k-1,n) = -rhos * Lfresh enddo + if (tr_snow) then + do k = 1, nslyr + trcrn(i,j,nt_rsnw +k-1,n) = rsnw_fall + trcrn(i,j,nt_rhos +k-1,n) = rhos + trcrn(i,j,nt_smice+k-1,n) = rhos + trcrn(i,j,nt_smliq+k-1,n) = c0 + enddo ! nslyr + endif enddo enddo enddo diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90 index d65cf52d3..976e95361 100644 --- a/cicecore/cicedynB/general/ice_step_mod.F90 +++ b/cicecore/cicedynB/general/ice_step_mod.F90 @@ -36,7 +36,7 @@ module ice_step_mod private public :: step_therm1, step_therm2, step_dyn_horiz, step_dyn_ridge, & - prep_radiation, step_radiation, ocean_mixed_layer, & + step_snow, prep_radiation, step_radiation, ocean_mixed_layer, & update_state, biogeochemistry, save_init, step_dyn_wave !======================================================================= @@ -163,7 +163,7 @@ subroutine step_therm1 (dt, iblk) Cdn_ocn, Cdn_ocn_skin, Cdn_ocn_floe, Cdn_ocn_keel, Cdn_atm_ratio, & Cdn_atm, Cdn_atm_skin, Cdn_atm_floe, Cdn_atm_rdg, Cdn_atm_pond, & hfreebd, hdraft, hridge, distrdg, hkeel, dkeel, lfloe, dfloe, & - fswsfcn, fswintn, Sswabsn, Iswabsn, & + fswsfcn, fswintn, Sswabsn, Iswabsn, meltsliqn, meltsliq, & fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf use ice_blocks, only: block, get_block, nx_block, ny_block use ice_calendar, only: yday @@ -172,13 +172,13 @@ subroutine step_therm1 (dt, iblk) use ice_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fbot, Tbot, Tsnice, & meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm, fside, & wind, rhoa, potT, Qa, zlvl, zlvs, strax, stray, flatn, fsensn, fsurfn, fcondtopn, & - flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, & + flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, fsloss, & frain, Tair, strairxT, strairyT, fsurf, fcondtop, fsens, & flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, Uref, fresh, fsalt, fhocn, & fswthru, fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf, & meltt, melts, meltb, congel, snoice, & flatn_f, fsensn_f, fsurfn_f, fcondtopn_f, & - send_i2x_per_cat, fswthrun_ai + send_i2x_per_cat, fswthrun_ai, dsnow use ice_flux_bgc, only: dsnown, faero_atm, faero_ocn, fiso_atm, fiso_ocn, & Qa_iso, Qref_iso, fiso_evap, HDO_ocn, H2_16O_ocn, H2_18O_ocn use ice_grid, only: lmask_n, lmask_s, tmask @@ -211,11 +211,11 @@ subroutine step_therm1 (dt, iblk) integer (kind=int_kind) :: & ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, & nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_qsno, & - nt_isosno, nt_isoice + nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq logical (kind=log_kind) :: & tr_iage, tr_FY, tr_iso, tr_aero, tr_pond, tr_pond_cesm, & - tr_pond_lvl, tr_pond_topo, calc_Tsfc, highfreq + tr_pond_lvl, tr_pond_topo, calc_Tsfc, highfreq, tr_snow real (kind=dbl_kind) :: & uvel_center, & ! cell-centered velocity, x component (m/s) @@ -228,6 +228,9 @@ subroutine step_therm1 (dt, iblk) real (kind=dbl_kind), dimension(n_iso,ncat) :: & isosno, isoice ! kg/m^2 + real (kind=dbl_kind), dimension(nslyr,ncat) :: & + rsnwn, smicen, smliqn + type (block) :: & this_block ! block information for current block @@ -240,13 +243,15 @@ subroutine step_therm1 (dt, iblk) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_iso_out=tr_iso, & tr_aero_out=tr_aero, tr_pond_out=tr_pond, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo) + tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo, & + tr_snow_out=tr_snow) call icepack_query_tracer_indices( & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & nt_iage_out=nt_iage, nt_FY_out=nt_FY, & nt_qice_out=nt_qice, nt_sice_out=nt_sice, & nt_aero_out=nt_aero, nt_qsno_out=nt_qsno, & + nt_rsnw_out=nt_rsnw, nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & 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, & @@ -256,7 +261,9 @@ subroutine step_therm1 (dt, iblk) prescribed_ice = .false. #endif - isosno (:,:) = c0 + rsnwn (:,:) = c0 + smicen (:,:) = c0 + smliqn (:,:) = c0 isoice (:,:) = c0 aerosno(:,:,:) = c0 aeroice(:,:,:) = c0 @@ -302,6 +309,16 @@ subroutine step_therm1 (dt, iblk) vvel_center = c0 endif ! highfreq + if (tr_snow) then + do n = 1, ncat + do k = 1, nslyr + rsnwn (k,n) = trcrn(i,j,nt_rsnw +k-1,n,iblk) + smicen(k,n) = trcrn(i,j,nt_smice+k-1,n,iblk) + smliqn(k,n) = trcrn(i,j,nt_smliq+k-1,n,iblk) + enddo + enddo + endif ! tr_snow + if (tr_iso) then ! trcrn(nt_iso*) has units kg/m^3 do n=1,ncat do k=1,n_iso @@ -350,6 +367,9 @@ subroutine step_therm1 (dt, iblk) ipnd = trcrn (i,j,nt_ipnd,:,iblk), & iage = trcrn (i,j,nt_iage,:,iblk), & FY = trcrn (i,j,nt_FY ,:,iblk), & + rsnwn = rsnwn (:,:), & + smicen = smicen (:,:), & + smliqn = smliqn (:,:), & aerosno = aerosno (:,:,:), & aeroice = aeroice (:,:,:), & isosno = isosno (:,:), & @@ -397,13 +417,14 @@ subroutine step_therm1 (dt, iblk) strocnyT = strocnyT (i,j, iblk), & fbot = fbot (i,j, iblk), & Tbot = Tbot (i,j, iblk), & - Tsnice = Tsnice (i,j, iblk), & + Tsnice = Tsnice (i,j, iblk), & frzmlt = frzmlt (i,j, iblk), & rside = rside (i,j, iblk), & fside = fside (i,j, iblk), & fsnow = fsnow (i,j, iblk), & frain = frain (i,j, iblk), & fpond = fpond (i,j, iblk), & + fsloss = fsloss (i,j, iblk), & fsurf = fsurf (i,j, iblk), & fsurfn = fsurfn (i,j,:,iblk), & fcondtop = fcondtop (i,j, iblk), & @@ -433,10 +454,10 @@ subroutine step_therm1 (dt, iblk) fsalt = fsalt (i,j, iblk), & fhocn = fhocn (i,j, iblk), & fswthru = fswthru (i,j, iblk), & - fswthru_vdr = fswthru_vdr (i,j, iblk),& - fswthru_vdf = fswthru_vdf (i,j, iblk),& - fswthru_idr = fswthru_idr (i,j, iblk),& - fswthru_idf = fswthru_idf (i,j, iblk),& + fswthru_vdr = fswthru_vdr (i,j, iblk), & + fswthru_vdf = fswthru_vdf (i,j, iblk), & + fswthru_idr = fswthru_idr (i,j, iblk), & + fswthru_idf = fswthru_idf (i,j, iblk), & flatn_f = flatn_f (i,j,:,iblk), & fsensn_f = fsensn_f (i,j,:,iblk), & fsurfn_f = fsurfn_f (i,j,:,iblk), & @@ -461,7 +482,10 @@ subroutine step_therm1 (dt, iblk) congeln = congeln (i,j,:,iblk), & snoice = snoice (i,j, iblk), & snoicen = snoicen (i,j,:,iblk), & + dsnow = dsnow (i,j, iblk), & dsnown = dsnown (i,j,:,iblk), & + meltsliq = meltsliq (i,j, iblk), & + meltsliqn = meltsliqn (i,j,:,iblk), & lmask_n = lmask_n (i,j, iblk), & lmask_s = lmask_s (i,j, iblk), & mlt_onset = mlt_onset (i,j, iblk), & @@ -483,6 +507,16 @@ subroutine step_therm1 (dt, iblk) endif + if (tr_snow) then + do n = 1, ncat + do k = 1, nslyr + trcrn(i,j,nt_rsnw +k-1,n,iblk) = rsnwn (k,n) + trcrn(i,j,nt_smice+k-1,n,iblk) = smicen(k,n) + trcrn(i,j,nt_smliq+k-1,n,iblk) = smliqn(k,n) + enddo + enddo + endif ! tr_snow + if (tr_iso) then do n = 1, ncat if (vicen(i,j,n,iblk) > puny) & @@ -685,13 +719,15 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset) use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound real (kind=dbl_kind), intent(in) :: & - dt , & ! time step - offset ! d(age)/dt time offset = dt for thermo, 0 for dyn + dt ! time step - real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & - daidt, & ! change in ice area per time step - dvidt, & ! change in ice volume per time step - dagedt ! change in ice age per time step + real (kind=dbl_kind), dimension(:,:,:), intent(inout), optional :: & + daidt, & ! change in ice area per time step + dvidt, & ! change in ice volume per time step + dagedt ! change in ice age per time step + + real (kind=dbl_kind), intent(in), optional :: & + offset ! d(age)/dt time offset = dt for thermo, 0 for dyn integer (kind=int_kind) :: & iblk, & ! block index @@ -747,6 +783,8 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset) n_trcr_strata = n_trcr_strata(:), & nt_strata = nt_strata(:,:)) + if (present(offset)) then + !----------------------------------------------------------------- ! Compute thermodynamic area and volume tendencies. !----------------------------------------------------------------- @@ -762,7 +800,8 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset) dagedt(i,j,iblk) = (trcr(i,j,nt_iage,iblk) & - dagedt(i,j,iblk)) / dt endif - endif + endif ! tr_iage + endif ! present(offset) enddo ! i enddo ! j @@ -1022,6 +1061,118 @@ subroutine step_dyn_ridge (dt, ndtd, iblk) end subroutine step_dyn_ridge +!======================================================================= +! +! Updates snow tracers +! +! authors: Elizabeth C. Hunke, LANL +! Nicole Jeffery, LANL + + subroutine step_snow (dt, iblk) + + use ice_blocks, only: block, get_block + use ice_calendar, only: nstreams + use ice_domain, only: blocks_ice + use ice_domain_size, only: ncat, nslyr, nilyr + use ice_flux, only: snwcnt, wind, fresh, fhocn, fsloss, fsnow + use ice_state, only: trcrn, vsno, vsnon, vicen, aicen, aice + use icepack_intfc, only: icepack_step_snow + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + integer (kind=int_kind), intent(in) :: & + iblk ! block index + + ! local variables + + integer (kind=int_kind) :: & + nt_smice, nt_smliq, nt_rsnw, & + nt_Tsfc, nt_qice, nt_sice, nt_qsno, & + nt_alvl, nt_vlvl, nt_rhos + + integer (kind=int_kind) :: & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + i, j, & ! horizontal indices + n, & ! category index + ns, & ! history streams index + ipoint ! index for print diagnostic + + real (kind=dbl_kind) :: & + puny + + real (kind=dbl_kind) :: & + fhs ! flag for presence of snow + + character(len=*), parameter :: subname = '(step_snow)' + + type (block) :: & + this_block ! block information for current block + + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + !----------------------------------------------------------------- + ! query icepack values + !----------------------------------------------------------------- + + call icepack_query_parameters(puny_out=puny) + call icepack_query_tracer_indices( & + nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & + nt_rsnw_out=nt_rsnw, nt_Tsfc_out=nt_Tsfc, & + nt_qice_out=nt_qice, nt_sice_out=nt_sice, nt_qsno_out=nt_qsno, & + nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_rhos_out=nt_rhos) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + !----------------------------------------------------------------- + ! Snow redistribution and metamorphosis + !----------------------------------------------------------------- + + do j = jlo, jhi + do i = ilo, ihi + + call icepack_step_snow (dt, nilyr, & + nslyr, ncat, & + wind (i,j, iblk), & + aice (i,j, iblk), & + aicen(i,j,:,iblk), & + vicen(i,j,:,iblk), & + vsnon(i,j,:,iblk), & + trcrn(i,j,nt_Tsfc,:,iblk), & + trcrn(i,j,nt_qice,:,iblk), & ! top layer only + trcrn(i,j,nt_sice,:,iblk), & ! top layer only + trcrn(i,j,nt_qsno:nt_qsno+nslyr-1,:,iblk), & + trcrn(i,j,nt_alvl,:,iblk), & + trcrn(i,j,nt_vlvl,:,iblk), & + trcrn(i,j,nt_smice:nt_smice+nslyr-1,:,iblk), & + trcrn(i,j,nt_smliq:nt_smliq+nslyr-1,:,iblk), & + trcrn(i,j,nt_rsnw:nt_rsnw+nslyr-1,:,iblk), & + trcrn(i,j,nt_rhos:nt_rhos+nslyr-1,:,iblk), & + fresh (i,j,iblk), & + fhocn (i,j,iblk), & + fsloss (i,j,iblk), & + fsnow (i,j,iblk)) + enddo + enddo + + ! increment counter for history averaging + do j = jlo, jhi + do i = ilo, ihi + fhs = c0 + if (vsno(i,j,iblk) > puny) fhs = c1 + do ns = 1, nstreams + snwcnt(i,j,iblk,ns) = snwcnt(i,j,iblk,ns) + fhs + enddo + enddo + enddo + + end subroutine step_snow + !======================================================================= ! ! Computes radiation fields @@ -1067,7 +1218,7 @@ subroutine step_radiation (dt, iblk) this_block ! block information for current block integer (kind=int_kind) :: & - nt_Tsfc, nt_alvl, & + nt_Tsfc, nt_alvl, nt_rsnw, & nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nlt_chl_sw, & ntrcr, nbtrcr, nbtrcr_sw, nt_fbri @@ -1078,13 +1229,14 @@ subroutine step_radiation (dt, iblk) nlt_zaero_sw, nt_zaero logical (kind=log_kind) :: & - tr_bgc_N, tr_zaero, tr_brine, dEdd_algae, modal_aero + tr_bgc_N, tr_zaero, tr_brine, dEdd_algae, modal_aero, snwgrain real (kind=dbl_kind), dimension(ncat) :: & - fbri ! brine height to ice thickness + fbri ! brine height to ice thickness real(kind= dbl_kind), dimension(:,:), allocatable :: & - ztrcr_sw + ztrcr_sw, & ! zaerosols (kg/m^3) and chla (mg/m^3) + rsnow ! snow grain radius tracer (10^-6 m) logical (kind=log_kind) :: & debug, & ! flag for printing debugging information @@ -1099,16 +1251,18 @@ subroutine step_radiation (dt, iblk) call icepack_query_tracer_flags( & tr_brine_out=tr_brine, tr_bgc_N_out=tr_bgc_N, tr_zaero_out=tr_zaero) call icepack_query_tracer_indices( & - nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, & + nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, nt_rsnw_out=nt_rsnw, & 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_bgc_N_out=nt_bgc_N) - call icepack_query_parameters(dEdd_algae_out=dEdd_algae, modal_aero_out=modal_aero) + call icepack_query_parameters(dEdd_algae_out=dEdd_algae, modal_aero_out=modal_aero, & + snwgrain_out=snwgrain) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) allocate(ztrcr_sw(nbtrcr_sw,ncat)) + allocate(rsnow(nslyr,ncat)) this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -1130,10 +1284,16 @@ subroutine step_radiation (dt, iblk) write (nu_diag, *) 'my_task = ',my_task enddo ! ipoint endif - fbri(:) = c0 + fbri (:) = c0 ztrcr_sw(:,:) = c0 + rsnow (:,:) = c0 do n = 1, ncat - if (tr_brine) fbri(n) = trcrn(i,j,nt_fbri,n,iblk) + if (tr_brine) fbri(n) = trcrn(i,j,nt_fbri,n,iblk) + if (snwgrain) then + do k = 1, nslyr + rsnow(k,n) = trcrn(i,j,nt_rsnw+k-1,n,iblk) + enddo + endif enddo if (tmask(i,j,iblk)) then @@ -1182,8 +1342,7 @@ subroutine step_radiation (dt, iblk) albpndn =albpndn (i,j,: ,iblk), apeffn =apeffn (i,j,: ,iblk), & snowfracn=snowfracn(i,j,: ,iblk), & dhsn =dhsn (i,j,: ,iblk), ffracn =ffracn(i,j,:,iblk), & - l_print_point=l_print_point) - + rsnow =rsnow (:,:), l_print_point=l_print_point) endif if (dEdd_algae .and. (tr_zaero .or. tr_bgc_N)) then @@ -1202,6 +1361,7 @@ subroutine step_radiation (dt, iblk) file=__FILE__, line=__LINE__) deallocate(ztrcr_sw) + deallocate(rsnow) call ice_timer_stop(timer_sw) ! shortwave diff --git a/cicecore/cicedynB/infrastructure/ice_read_write.F90 b/cicecore/cicedynB/infrastructure/ice_read_write.F90 index d902c62f8..d488b5693 100644 --- a/cicecore/cicedynB/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedynB/infrastructure/ice_read_write.F90 @@ -68,6 +68,9 @@ module ice_read_write ice_read_nc_xyz, & !ice_read_nc_xyf, & ice_read_nc_point, & + ice_read_nc_1D, & + ice_read_nc_2D, & + ice_read_nc_3D, & ice_read_nc_z end interface @@ -1120,7 +1123,7 @@ subroutine ice_read_nc_xy(fid, nrec, varname, work, diag, & amin, amax, asum ! min, max values and sum of input array ! character (char_len) :: & -! dimname ! dimension name +! dimname ! dimension name real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g1 @@ -1293,7 +1296,7 @@ subroutine ice_read_nc_xyz(fid, nrec, varname, work, diag, & amin, amax, asum ! min, max values and sum of input array ! character (char_len) :: & -! dimname ! dimension name +! dimname ! dimension name real (kind=dbl_kind), dimension(:,:,:), allocatable :: & work_g1 @@ -1472,7 +1475,7 @@ subroutine ice_read_nc_xyf(fid, nrec, varname, work, diag, & amin, amax, asum ! min, max values and sum of input array character (char_len) :: & - dimname ! dimension name + dimname ! dimension name real (kind=dbl_kind), dimension(:,:,:), allocatable :: & work_g1 @@ -1646,7 +1649,7 @@ subroutine ice_read_nc_point(fid, nrec, varname, work, diag, & workg ! temporary work variable character (char_len) :: & - dimname ! dimension name + dimname ! dimension name if (my_task == master_task) then @@ -1700,6 +1703,277 @@ end subroutine ice_read_nc_point !======================================================================= +! Written by T. Craig + + subroutine ice_read_nc_1D(fid, varname, work, diag, & + xdim) + + use ice_fileunits, only: nu_diag + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + xdim ! field dimensions + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (char_len), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), dimension(:), intent(out) :: & + work ! output array + + ! local variables + + character(len=*), parameter :: subname = '(ice_read_nc_1D)' + +#ifdef USE_NETCDF +! netCDF file diagnostics: + integer (kind=int_kind) :: & + varid, & ! netcdf id for field + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + dimlen ! size of dimension + + character (char_len) :: & + dimname ! dimension name + + real (kind=dbl_kind), dimension(xdim) :: & + workg ! output array (real, 8-byte) + + !-------------------------------------------------------------- + + if (my_task == master_task) then + + if (size(work,dim=1) < xdim) then + write(nu_diag,*) subname,' work, dim=1 ',size(work,dim=1),xdim + call abort_ice (subname//' ERROR: work array wrong size '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + !------------------------------------------------------------- + ! Find out ID of required variable + !------------------------------------------------------------- + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call abort_ice (subname//' ERROR: Cannot find variable '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + + !-------------------------------------------------------------- + ! Read array + !-------------------------------------------------------------- + status = nf90_get_var( fid, varid, workg, & + start=(/1/), & + count=(/xdim/) ) + work(1:xdim) = workg(1:xdim) + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (diag) then + write(nu_diag,*) subname, & + ' fid= ',fid, ', xdim = ',xdim, & + ' varname = ',trim(varname) + status = nf90_inquire(fid, nDimensions=ndim, nVariables=nvar) + write(nu_diag,*) 'ndim= ',ndim,', nvar= ',nvar + endif + endif +#else + call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & + file=__FILE__, line=__LINE__) + work = c0 ! to satisfy intent(out) attribute +#endif + + end subroutine ice_read_nc_1D + +!======================================================================= + +! Written by T. Craig + + subroutine ice_read_nc_2D(fid, varname, work, diag, & + xdim, ydim) + + use ice_fileunits, only: nu_diag + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + xdim, ydim ! field dimensions + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (char_len), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), dimension(:,:), intent(out) :: & + work ! output array + + ! local variables + + character(len=*), parameter :: subname = '(ice_read_nc_2D)' + +#ifdef USE_NETCDF +! netCDF file diagnostics: + integer (kind=int_kind) :: & + varid, & ! netcdf id for field + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + dimlen ! size of dimension + + character (char_len) :: & + dimname ! dimension name + + real (kind=dbl_kind), dimension(xdim,ydim) :: & + workg ! output array (real, 8-byte) + + !-------------------------------------------------------------- + + if (my_task == master_task) then + + if (size(work,dim=1) < xdim .or. & + size(work,dim=2) < ydim) then + write(nu_diag,*) subname,' work, dim=1 ',size(work,dim=1),xdim + write(nu_diag,*) subname,' work, dim=2 ',size(work,dim=2),ydim + call abort_ice (subname//' ERROR: work array wrong size '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + !------------------------------------------------------------- + ! Find out ID of required variable + !------------------------------------------------------------- + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call abort_ice (subname//' ERROR: Cannot find variable '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + + !-------------------------------------------------------------- + ! Read array + !-------------------------------------------------------------- + status = nf90_get_var( fid, varid, workg, & + start=(/1,1/), & + count=(/xdim,ydim/) ) + work(1:xdim,1:ydim) = workg(1:xdim, 1:ydim) + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (diag) then + write(nu_diag,*) subname, & + ' fid= ',fid, ', xdim = ',xdim, & + ' ydim= ', ydim, ' varname = ',trim(varname) + status = nf90_inquire(fid, nDimensions=ndim, nVariables=nvar) + write(nu_diag,*) 'ndim= ',ndim,', nvar= ',nvar + endif + endif +#else + call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & + file=__FILE__, line=__LINE__) + work = c0 ! to satisfy intent(out) attribute +#endif + + end subroutine ice_read_nc_2D + +!======================================================================= +!======================================================================= + +! Written by T. Craig + + subroutine ice_read_nc_3D(fid, varname, work, diag, & + xdim, ydim, zdim) + + use ice_fileunits, only: nu_diag + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + xdim, ydim,zdim ! field dimensions + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (char_len), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), dimension(:,:,:), intent(out) :: & + work ! output array + + ! local variables + + character(len=*), parameter :: subname = '(ice_read_nc_3D)' + +#ifdef USE_NETCDF +! netCDF file diagnostics: + integer (kind=int_kind) :: & + varid, & ! netcdf id for field + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + dimlen ! size of dimension + + character (char_len) :: & + dimname ! dimension name + + real (kind=dbl_kind), dimension(xdim,ydim,zdim) :: & + workg ! output array (real, 8-byte) + + !-------------------------------------------------------------- + + if (my_task == master_task) then + + if (size(work,dim=1) < xdim .or. & + size(work,dim=2) < ydim .or. & + size(work,dim=3) < zdim ) then + write(nu_diag,*) subname,' work, dim=1 ',size(work,dim=1),xdim + write(nu_diag,*) subname,' work, dim=2 ',size(work,dim=2),ydim + write(nu_diag,*) subname,' work, dim=3 ',size(work,dim=3),zdim + call abort_ice (subname//' ERROR: work array wrong size '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + !------------------------------------------------------------- + ! Find out ID of required variable + !------------------------------------------------------------- + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call abort_ice (subname//' ERROR: Cannot find variable '//trim(varname), & + file=__FILE__, line=__LINE__ ) + endif + + !-------------------------------------------------------------- + ! Read array + !-------------------------------------------------------------- + status = nf90_get_var( fid, varid, workg, & + start=(/1,1,1/), & + count=(/xdim,ydim,zdim/) ) + work(1:xdim,1:ydim,1:zdim) = workg(1:xdim, 1:ydim, 1:zdim) + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (diag) then + write(nu_diag,*) subname, & + ' fid= ',fid, ', xdim = ',xdim, & + ' ydim= ', ydim,' zdim = ',zdim, ' varname = ',trim(varname) + status = nf90_inquire(fid, nDimensions=ndim, nVariables=nvar) + write(nu_diag,*) 'ndim= ',ndim,', nvar= ',nvar + endif + endif +#else + call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & + file=__FILE__, line=__LINE__) + work = c0 ! to satisfy intent(out) attribute +#endif + + end subroutine ice_read_nc_3D + +!======================================================================= + ! Adapted by Nicole Jeffery, LANL subroutine ice_read_nc_z(fid, nrec, varname, work, diag, & @@ -1739,7 +2013,7 @@ subroutine ice_read_nc_z(fid, nrec, varname, work, diag, & dimlen ! size of dimension character (char_len) :: & - dimname ! dimension name + dimname ! dimension name #endif character(len=*), parameter :: subname = '(ice_read_nc_z)' @@ -1791,7 +2065,7 @@ subroutine ice_read_nc_z(fid, nrec, varname, work, diag, & #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) work = c0 ! to satisfy intent(out) attribute #endif end subroutine ice_read_nc_z @@ -1841,7 +2115,7 @@ subroutine ice_write_nc_xy(fid, nrec, varid, work, diag, & character (char_len) :: & lvarname ! variable name -! dimname ! dimension name +! dimname ! dimension name real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g1 @@ -1914,7 +2188,7 @@ subroutine ice_write_nc_xy(fid, nrec, varid, work, diag, & #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) #endif end subroutine ice_write_nc_xy @@ -1965,7 +2239,7 @@ subroutine ice_write_nc_xyz(fid, nrec, varid, work, diag, & character (char_len) :: & lvarname ! variable name -! dimname ! dimension name +! dimname ! dimension name real (kind=dbl_kind), dimension(:,:,:), allocatable :: & work_g1 @@ -2048,7 +2322,7 @@ subroutine ice_write_nc_xyz(fid, nrec, varid, work, diag, & #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) #endif end subroutine ice_write_nc_xyz @@ -2094,7 +2368,7 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) amin, amax, asum ! min, max values and sum of input array ! character (char_len) :: & -! dimname ! dimension name +! dimname ! dimension name ! real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g3 @@ -2162,7 +2436,7 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) work_g = c0 ! to satisfy intent(out) attribute #endif @@ -2191,7 +2465,7 @@ subroutine ice_close_nc(fid) endif #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) #endif end subroutine ice_close_nc @@ -2249,7 +2523,7 @@ subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, & amin, amax, asum ! min, max values and sum of input array ! character (char_len) :: & -! dimname ! dimension name +! dimname ! dimension name real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g1 @@ -2328,7 +2602,7 @@ subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, & #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) work = c0 ! to satisfy intent(out) attribute #endif @@ -2406,7 +2680,7 @@ subroutine ice_read_vec_nc (fid, nrec, varname, work_g, diag) #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) work_g = c0 ! to satisfy intent(out) attribute #endif @@ -2452,7 +2726,7 @@ subroutine ice_get_ncvarsize(fid,varname,recsize) endif #else call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) recsize = 0 ! to satisfy intent(out) attribute #endif diff --git a/cicecore/cicedynB/infrastructure/io/io_binary/ice_restart.F90 b/cicecore/cicedynB/infrastructure/io/io_binary/ice_restart.F90 index 91d57ea48..a6f42a6a5 100644 --- a/cicecore/cicedynB/infrastructure/io/io_binary/ice_restart.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_binary/ice_restart.F90 @@ -15,11 +15,12 @@ module ice_restart use ice_fileunits, only: nu_diag, nu_rst_pointer use ice_fileunits, only: nu_dump, nu_dump_eap, nu_dump_FY, nu_dump_age use ice_fileunits, only: nu_dump_lvl, nu_dump_pond, nu_dump_hbrine - use ice_fileunits, only: nu_dump_bgc, nu_dump_aero, nu_dump_fsd, nu_dump_iso + use ice_fileunits, only: nu_dump_iso, nu_dump_snow + use ice_fileunits, only: nu_dump_bgc, nu_dump_aero, nu_dump_fsd use ice_fileunits, only: nu_restart, nu_restart_eap, nu_restart_FY, nu_restart_age use ice_fileunits, only: nu_restart_lvl, nu_restart_pond, nu_restart_hbrine use ice_fileunits, only: nu_restart_bgc, nu_restart_aero, nu_restart_fsd - use ice_fileunits, only: nu_restart_iso + use ice_fileunits, only: nu_restart_iso, nu_restart_snow use ice_exit, only: abort_ice use icepack_intfc, only: icepack_query_parameters use icepack_intfc, only: icepack_query_tracer_sizes @@ -57,7 +58,7 @@ subroutine init_restart_read(ice_ic) logical (kind=log_kind) :: & solve_zsal, tr_fsd, & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_pond_cesm, & - tr_pond_topo, tr_pond_lvl, tr_brine + tr_pond_topo, tr_pond_lvl, tr_brine, tr_snow character(len=char_len_long) :: & filename, filename0 @@ -82,7 +83,8 @@ subroutine init_restart_read(ice_ic) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_fsd_out=tr_fsd, & tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, tr_brine_out=tr_brine) + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_snow_out=tr_snow, tr_brine_out=tr_brine) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -285,6 +287,26 @@ subroutine init_restart_read(ice_ic) endif endif + if (tr_snow) then + if (my_task == master_task) then + n = index(filename0,trim(restart_file)) + if (n == 0) call abort_ice(subname//'ERROR: snow restart: filename discrepancy') + string1 = trim(filename0(1:n-1)) + string2 = trim(filename0(n+lenstr(restart_file):lenstr(filename0))) + write(filename,'(a,a,a,a)') & + string1(1:lenstr(string1)), & + restart_file(1:lenstr(restart_file)),'.snow', & + string2(1:lenstr(string2)) + if (restart_ext) then + call ice_open_ext(nu_restart_snow,filename,0) + else + call ice_open(nu_restart_snow,filename,0) + endif + read (nu_restart_snow) iignore,rignore,rignore + write(nu_diag,*) 'Reading ',filename(1:lenstr(filename)) + endif + endif + if (tr_brine) then if (my_task == master_task) then n = index(filename0,trim(restart_file)) @@ -392,7 +414,7 @@ subroutine init_restart_write(filename_spec) logical (kind=log_kind) :: & solve_zsal, tr_fsd, & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_pond_cesm, & - tr_pond_topo, tr_pond_lvl, tr_brine + tr_pond_topo, tr_pond_lvl, tr_brine, tr_snow integer (kind=int_kind) :: & nbtrcr ! number of bgc tracers @@ -408,7 +430,8 @@ subroutine init_restart_write(filename_spec) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_fsd_out=tr_fsd, & tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, tr_brine_out=tr_brine) + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_snow_out=tr_snow, tr_brine_out=tr_brine) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -599,6 +622,26 @@ subroutine init_restart_write(filename_spec) endif + if (tr_snow) then + + write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & + restart_dir(1:lenstr(restart_dir)), & + restart_file(1:lenstr(restart_file)),'.snow.', & + myear,'-',mmonth,'-',mday,'-',msec + + if (restart_ext) then + call ice_open_ext(nu_dump_snow,filename,0) + else + call ice_open(nu_dump_snow,filename,0) + endif + + if (my_task == master_task) then + write(nu_dump_snow) istep1,timesecs,time_forc + write(nu_diag,*) 'Writing ',filename(1:lenstr(filename)) + endif + + endif + if (tr_brine) then write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & @@ -808,7 +851,7 @@ subroutine final_restart() logical (kind=log_kind) :: & solve_zsal, & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_pond_cesm, & - tr_pond_topo, tr_pond_lvl, tr_brine + tr_pond_topo, tr_pond_lvl, tr_brine, tr_snow integer (kind=int_kind) :: & nbtrcr ! number of bgc tracers @@ -822,7 +865,8 @@ subroutine final_restart() call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, & tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, tr_brine_out=tr_brine) + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_snow_out=tr_snow, tr_brine_out=tr_brine) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -838,6 +882,7 @@ subroutine final_restart() if (tr_pond_cesm) close(nu_dump_pond) if (tr_pond_lvl) close(nu_dump_pond) if (tr_pond_topo) close(nu_dump_pond) + if (tr_snow) close(nu_dump_snow) if (tr_brine) close(nu_dump_hbrine) if (solve_zsal .or. nbtrcr > 0) & close(nu_dump_bgc) diff --git a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 index e744caf09..f6002ff40 100644 --- a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 @@ -145,7 +145,7 @@ subroutine init_restart_write(filename_spec) logical (kind=log_kind) :: & solve_zsal, skl_bgc, z_tracers, tr_fsd, & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_pond_cesm, & - tr_pond_topo, tr_pond_lvl, tr_brine, & + tr_pond_topo, tr_pond_lvl, tr_brine, tr_snow, & tr_bgc_N, tr_bgc_C, tr_bgc_Nit, & tr_bgc_Sil, tr_bgc_DMS, & tr_bgc_chl, tr_bgc_Am, & @@ -181,7 +181,8 @@ subroutine init_restart_write(filename_spec) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_fsd_out=tr_fsd, & tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, tr_brine_out=tr_brine, & + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_snow_out=tr_snow, tr_brine_out=tr_brine, & tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, tr_bgc_Nit_out=tr_bgc_Nit, & tr_bgc_Sil_out=tr_bgc_Sil, tr_bgc_DMS_out=tr_bgc_DMS, & tr_bgc_chl_out=tr_bgc_chl, tr_bgc_Am_out=tr_bgc_Am, & @@ -480,6 +481,16 @@ subroutine init_restart_write(filename_spec) call define_rest_field(ncid,'qsno'//trim(nchar),dims) enddo + if (tr_snow) then + do k=1,nslyr + write(nchar,'(i3.3)') k + call define_rest_field(ncid,'smice'//trim(nchar),dims) + call define_rest_field(ncid,'smliq'//trim(nchar),dims) + call define_rest_field(ncid, 'rhos'//trim(nchar),dims) + call define_rest_field(ncid, 'rsnw'//trim(nchar),dims) + enddo + endif + if (tr_fsd) then do k=1,nfsd write(nchar,'(i3.3)') k diff --git a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 index 12d5d8e71..c6d6a02af 100644 --- a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 @@ -151,7 +151,7 @@ subroutine init_restart_write(filename_spec) logical (kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_pond_cesm, & - tr_pond_topo, tr_pond_lvl, tr_brine, & + tr_pond_topo, tr_pond_lvl, tr_brine, tr_snow, & tr_bgc_N, tr_bgc_C, tr_bgc_Nit, & tr_bgc_Sil, tr_bgc_DMS, & tr_bgc_chl, tr_bgc_Am, & @@ -187,7 +187,8 @@ subroutine init_restart_write(filename_spec) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, & tr_iso_out=tr_iso, tr_aero_out=tr_aero, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, tr_brine_out=tr_brine, & + tr_pond_topo_out=tr_pond_topo, tr_pond_lvl_out=tr_pond_lvl, & + tr_snow_out=tr_snow, tr_brine_out=tr_brine, & tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, tr_bgc_Nit_out=tr_bgc_Nit, & tr_bgc_Sil_out=tr_bgc_Sil, tr_bgc_DMS_out=tr_bgc_DMS, & tr_bgc_chl_out=tr_bgc_chl, tr_bgc_Am_out=tr_bgc_Am, & @@ -483,6 +484,16 @@ subroutine init_restart_write(filename_spec) call define_rest_field(File,'qsno'//trim(nchar),dims) enddo + if (tr_snow) then + do k=1,nslyr + write(nchar,'(i3.3)') k + call define_rest_field(File,'smice'//trim(nchar),dims) + call define_rest_field(File,'smliq'//trim(nchar),dims) + call define_rest_field(File, 'rhos'//trim(nchar),dims) + call define_rest_field(File, 'rsnw'//trim(nchar),dims) + enddo + endif + if (tr_fsd) then do k=1,nfsd write(nchar,'(i3.3)') k diff --git a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 index 60f71fa8a..363025b9b 100644 --- a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 +++ b/cicecore/drivers/standalone/cice/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, & @@ -76,7 +77,7 @@ subroutine cice_init 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 @@ -90,7 +91,8 @@ subroutine cice_init use ice_transport_driver, only: init_transport 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)' call init_communicate ! initial setup for message passing @@ -162,7 +164,7 @@ subroutine cice_init 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__) @@ -176,7 +178,7 @@ subroutine cice_init 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__) @@ -207,8 +209,20 @@ subroutine cice_init call get_forcing_atmo ! atmospheric forcing from data call get_forcing_ocn(dt) ! ocean forcing from data + ! 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) @@ -235,12 +249,12 @@ 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, & @@ -248,6 +262,7 @@ subroutine init_restart 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, & @@ -262,12 +277,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)' @@ -282,10 +298,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, & @@ -382,6 +400,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. @@ -398,7 +432,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 diff --git a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 index 08059435f..0fde18e04 100644 --- a/cicecore/drivers/standalone/cice/CICE_RunMod.F90 +++ b/cicecore/drivers/standalone/cice/CICE_RunMod.F90 @@ -151,12 +151,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 @@ -170,7 +171,7 @@ 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 @@ -191,7 +192,7 @@ subroutine ice_step 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__) @@ -317,17 +318,28 @@ subroutine ice_step enddo endif + call ice_timer_start(timer_column) ! column physics + call ice_timer_start(timer_thermo) ! thermodynamics + !----------------------------------------------------------------- - ! albedo, shortwave radiation + ! snow redistribution and metamorphosis !----------------------------------------------------------------- - call ice_timer_start(timer_column) ! column physics - call ice_timer_start(timer_thermo) ! thermodynamics + 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 @@ -383,6 +395,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 diff --git a/cicecore/shared/ice_arrays_column.F90 b/cicecore/shared/ice_arrays_column.F90 index 46ea6f62e..dbad4292c 100644 --- a/cicecore/shared/ice_arrays_column.F90 +++ b/cicecore/shared/ice_arrays_column.F90 @@ -67,6 +67,15 @@ module ice_arrays_column character (len=35), public, allocatable :: c_hi_range(:) + ! icepack_snow.F90 + real (kind=dbl_kind), public, & + dimension (:,:,:), allocatable :: & + meltsliq ! snow melt mass (kg/m^2/step-->kg/m^2/day) + + real (kind=dbl_kind), public, & + dimension (:,:,:,:), allocatable :: & + meltsliqn ! snow melt mass in category n (kg/m^2) + ! icepack_meltpond_lvl.F90 real (kind=dbl_kind), public, & dimension (:,:,:,:), allocatable :: & @@ -354,6 +363,8 @@ subroutine alloc_arrays_column fzsal_g (nx_block,ny_block,max_blocks), & ! Total gravity drainage flux upNO (nx_block,ny_block,max_blocks), & ! nitrate uptake rate (mmol/m^2/d) times aice upNH (nx_block,ny_block,max_blocks), & ! ammonium uptake rate (mmol/m^2/d) times aice + meltsliq (nx_block,ny_block,max_blocks), & ! snow melt mass (kg/m^2) + meltsliqn (nx_block,ny_block,ncat,max_blocks), & ! snow melt mass in category n (kg/m^2) dhsn (nx_block,ny_block,ncat,max_blocks), & ! depth difference for snow on sea ice and pond ice ffracn (nx_block,ny_block,ncat,max_blocks), & ! fraction of fsurfn used to melt ipond alvdrn (nx_block,ny_block,ncat,max_blocks), & ! visible direct albedo (fraction) diff --git a/cicecore/shared/ice_fileunits.F90 b/cicecore/shared/ice_fileunits.F90 index b6b30d47a..ccb518807 100644 --- a/cicecore/shared/ice_fileunits.F90 +++ b/cicecore/shared/ice_fileunits.F90 @@ -51,6 +51,8 @@ module ice_fileunits nu_restart_lvl, & ! restart input file for level ice tracers nu_dump_pond , & ! dump file for restarting melt pond tracer nu_restart_pond,& ! restart input file for melt pond tracer + nu_dump_snow , & ! dump file for restarting snow redist/metamorph tracers + nu_restart_snow,& ! restart input file for snow redist/metamorph tracers nu_dump_fsd , & ! dump file for restarting floe size distribution nu_restart_fsd, & ! restart input file for floe size distribution nu_dump_iso , & ! dump file for restarting isotope tracers @@ -129,6 +131,8 @@ subroutine init_fileunits call get_fileunit(nu_restart_lvl) call get_fileunit(nu_dump_pond) call get_fileunit(nu_restart_pond) + call get_fileunit(nu_dump_snow) + call get_fileunit(nu_restart_snow) call get_fileunit(nu_dump_fsd) call get_fileunit(nu_restart_fsd) call get_fileunit(nu_dump_iso) @@ -218,6 +222,8 @@ subroutine release_all_fileunits call release_fileunit(nu_restart_lvl) call release_fileunit(nu_dump_pond) call release_fileunit(nu_restart_pond) + call release_fileunit(nu_dump_snow) + call release_fileunit(nu_restart_snow) call release_fileunit(nu_dump_fsd) call release_fileunit(nu_restart_fsd) call release_fileunit(nu_dump_iso) diff --git a/cicecore/shared/ice_init_column.F90 b/cicecore/shared/ice_init_column.F90 index 4f4641467..eff39a464 100644 --- a/cicecore/shared/ice_init_column.F90 +++ b/cicecore/shared/ice_init_column.F90 @@ -46,7 +46,7 @@ module ice_init_column init_age, init_FY, init_lvl, init_fsd, & init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & init_aerosol, init_bgc, init_hbrine, init_zbgc, input_zbgc, & - count_tracers, init_isotope + count_tracers, init_isotope, init_snowtracers ! namelist parameters needed locally @@ -214,8 +214,9 @@ subroutine init_shortwave logical (kind=log_kind) :: & l_print_point, & ! flag to print designated grid point diagnostics debug, & ! if true, print diagnostics - dEdd_algae, & ! from icepack - modal_aero ! from icepack + dEdd_algae, & ! use prognostic chla in dEdd radiation + modal_aero, & ! use modal aerosol optical treatment + snwgrain ! use variable snow radius character (char_len) :: shortwave @@ -225,12 +226,13 @@ subroutine init_shortwave real (kind=dbl_kind), dimension(ncat) :: & fbri ! brine height to ice thickness - real(kind=dbl_kind), allocatable :: & - ztrcr_sw(:,:) ! + real(kind= dbl_kind), dimension(:,:), allocatable :: & + ztrcr_sw, & ! zaerosols (kg/m^3) and chla (mg/m^3) + rsnow ! snow grain radius tracer (10^-6 m) logical (kind=log_kind) :: tr_brine, tr_zaero, tr_bgc_n integer (kind=int_kind) :: nt_alvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero, & - nt_fbri, nt_tsfc, ntrcr, nbtrcr, nbtrcr_sw + nt_fbri, nt_tsfc, ntrcr, nbtrcr, nbtrcr_sw, nt_rsnw integer (kind=int_kind), dimension(icepack_max_algae) :: & nt_bgc_N integer (kind=int_kind), dimension(icepack_max_aero) :: & @@ -243,17 +245,19 @@ subroutine init_shortwave call icepack_query_parameters(shortwave_out=shortwave) call icepack_query_parameters(dEdd_algae_out=dEdd_algae) call icepack_query_parameters(modal_aero_out=modal_aero) + call icepack_query_parameters(snwgrain_out=snwgrain) call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr, nbtrcr_sw_out=nbtrcr_sw) call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_zaero_out=tr_zaero, & tr_bgc_n_out=tr_bgc_n) call icepack_query_tracer_indices(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, nt_fbri_out=nt_fbri, nt_tsfc_out=nt_tsfc, & - nt_bgc_N_out=nt_bgc_N, nt_zaero_out=nt_zaero) + nt_bgc_N_out=nt_bgc_N, nt_zaero_out=nt_zaero, nt_rsnw_out=nt_rsnw) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__,line= __LINE__) allocate(ztrcr_sw(nbtrcr_sw, ncat)) + allocate(rsnow(nslyr,ncat)) do iblk=1,nblocks @@ -330,8 +334,14 @@ subroutine init_shortwave fbri(:) = c0 ztrcr_sw(:,:) = c0 + rsnow (:,:) = c0 do n = 1, ncat - if (tr_brine) fbri(n) = trcrn(i,j,nt_fbri,n,iblk) + if (tr_brine) fbri(n) = trcrn(i,j,nt_fbri,n,iblk) + if (snwgrain) then + do k = 1, nslyr + rsnow(k,n) = trcrn(i,j,nt_rsnw+k-1,n,iblk) + enddo + endif enddo if (tmask(i,j,iblk)) then @@ -379,6 +389,7 @@ subroutine init_shortwave albpndn=albpndn(i,j,:,iblk), apeffn=apeffn(i,j,:,iblk), & snowfracn=snowfracn(i,j,:,iblk), & dhsn=dhsn(i,j,:,iblk), ffracn=ffracn(i,j,:,iblk), & + rsnow=rsnow(:,:), & l_print_point=l_print_point, & initonly = .true.) endif @@ -475,6 +486,7 @@ subroutine init_shortwave enddo ! iblk deallocate(ztrcr_sw) + deallocate(rsnow) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & @@ -587,6 +599,29 @@ end subroutine init_meltponds_topo !======================================================================= +! Initialize snow redistribution/metamorphosis tracers (call prior to reading restart data) + + subroutine init_snowtracers(smice, smliq, rhos_cmp, rsnw) + + real(kind=dbl_kind), dimension(:,:,:,:), intent(out) :: & + smice, smliq, rhos_cmp, rsnw + character(len=*),parameter :: subname='(init_snowtracers)' + + real (kind=dbl_kind) :: & + rsnw_fall, & ! snow grain radius of new fallen snow (10^-6 m) + rhos ! snow density (kg/m^3) + + call icepack_query_parameters(rsnw_fall_out=rsnw_fall, rhos_out=rhos) + + rsnw (:,:,:,:) = rsnw_fall + rhos_cmp(:,:,:,:) = rhos + smice (:,:,:,:) = rhos + smliq (:,:,:,:) = c0 + + end subroutine init_snowtracers + +!======================================================================= + ! Initialize floe size distribution tracer (call prior to reading restart data) subroutine init_fsd(floesize) @@ -1776,10 +1811,12 @@ subroutine count_tracers integer (kind=int_kind) :: ntrcr logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_snow logical (kind=log_kind) :: tr_iso, tr_pond_cesm, tr_pond_lvl, tr_pond_topo integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero integer (kind=int_kind) :: nt_fsd, nt_isosno, nt_isoice + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw integer (kind=int_kind) :: & nbtrcr, nbtrcr_sw, & @@ -1862,7 +1899,7 @@ subroutine count_tracers tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_pond_out=tr_pond, & 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_fsd_out=tr_fsd, & - tr_iso_out=tr_iso, & + tr_snow_out=tr_snow, tr_iso_out=tr_iso, & tr_bgc_Nit_out=tr_bgc_Nit, tr_bgc_Am_out =tr_bgc_Am, tr_bgc_Sil_out=tr_bgc_Sil, & tr_bgc_DMS_out=tr_bgc_DMS, tr_bgc_PON_out=tr_bgc_PON, & tr_bgc_N_out =tr_bgc_N, tr_bgc_C_out =tr_bgc_C, tr_bgc_chl_out=tr_bgc_chl, & @@ -1925,6 +1962,21 @@ subroutine count_tracers endif endif + nt_smice = 0 + nt_smliq = 0 + nt_rhos = 0 + nt_rsnw = 0 + if (tr_snow) then + nt_smice = ntrcr + 1 + ntrcr = ntrcr + nslyr ! mass of ice in nslyr snow layers + nt_smliq = ntrcr + 1 + ntrcr = ntrcr + nslyr ! mass of liquid in nslyr snow layers + nt_rhos = ntrcr + 1 + ntrcr = ntrcr + nslyr ! snow density in nslyr layers + nt_rsnw = ntrcr + 1 + ntrcr = ntrcr + nslyr ! snow grain radius in nslyr layers + endif + nt_fsd = 0 if (tr_fsd) then nt_fsd = ntrcr + 1 ! floe size distribution @@ -2212,7 +2264,7 @@ subroutine count_tracers !tcx, +1 here is the unused tracer, want to get rid of it ntrcr = ntrcr + 1 -!tcx, reset unusaed tracer index, eventually get rid of it. +!tcx, reset unused tracer index, eventually get rid of it. if (nt_iage <= 0) nt_iage = ntrcr if (nt_FY <= 0) nt_FY = ntrcr if (nt_alvl <= 0) nt_alvl = ntrcr @@ -2220,6 +2272,10 @@ subroutine count_tracers if (nt_apnd <= 0) nt_apnd = ntrcr if (nt_hpnd <= 0) nt_hpnd = ntrcr if (nt_ipnd <= 0) nt_ipnd = ntrcr + if (nt_smice <= 0) nt_smice = ntrcr + if (nt_smliq <= 0) nt_smliq = ntrcr + if (nt_rhos <= 0) nt_rhos = ntrcr + if (nt_rsnw <= 0) nt_rsnw = ntrcr if (nt_fsd <= 0) nt_fsd = ntrcr if (nt_isosno<= 0) nt_isosno= ntrcr if (nt_isoice<= 0) nt_isoice= ntrcr @@ -2246,6 +2302,7 @@ subroutine count_tracers nt_qice_in=nt_qice, nt_qsno_in=nt_qsno, nt_iage_in=nt_iage, nt_fy_in=nt_fy, & nt_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl, nt_apnd_in=nt_apnd, nt_hpnd_in=nt_hpnd, & nt_ipnd_in=nt_ipnd, nt_fsd_in=nt_fsd, nt_aero_in=nt_aero, & + nt_smice_in=nt_smice, nt_smliq_in=nt_smliq, nt_rhos_in=nt_rhos, nt_rsnw_in=nt_rsnw, & nt_isosno_in=nt_isosno, nt_isoice_in=nt_isoice, nt_fbri_in=nt_fbri, & nt_bgc_Nit_in=nt_bgc_Nit, nt_bgc_Am_in=nt_bgc_Am, nt_bgc_Sil_in=nt_bgc_Sil, & nt_bgc_DMS_in=nt_bgc_DMS, nt_bgc_PON_in=nt_bgc_PON, nt_bgc_S_in=nt_bgc_S, & diff --git a/cicecore/shared/ice_restart_column.F90 b/cicecore/shared/ice_restart_column.F90 index e819b1098..074b37dbe 100644 --- a/cicecore/shared/ice_restart_column.F90 +++ b/cicecore/shared/ice_restart_column.F90 @@ -12,7 +12,7 @@ module ice_restart_column use ice_communicate, only: my_task, master_task use ice_constants, only: c0, c1, p5 use ice_constants, only: field_loc_center, field_type_scalar - use ice_domain_size, only: ncat, nfsd, nblyr + use ice_domain_size, only: ncat, nslyr, nfsd, nblyr use ice_restart,only: read_restart_field, write_restart_field use ice_exit, only: abort_ice use ice_fileunits, only: nu_diag @@ -32,6 +32,7 @@ module ice_restart_column write_restart_pond_cesm, read_restart_pond_cesm, & write_restart_pond_lvl, read_restart_pond_lvl, & write_restart_pond_topo, read_restart_pond_topo, & + write_restart_snow, read_restart_snow, & write_restart_fsd, read_restart_fsd, & write_restart_iso, read_restart_iso, & write_restart_aero, read_restart_aero, & @@ -45,6 +46,7 @@ module ice_restart_column restart_pond_cesm, & ! if .true., read meltponds restart file restart_pond_lvl , & ! if .true., read meltponds restart file restart_pond_topo, & ! if .true., read meltponds restart file + restart_snow , & ! if .true., read snow tracer restart file restart_fsd , & ! if .true., read floe size restart file restart_iso , & ! if .true., read isotope tracer restart file restart_aero , & ! if .true., read aerosol tracer restart file @@ -483,6 +485,93 @@ end subroutine read_restart_pond_topo !======================================================================= +! Dumps all values needed for restarting snow redistribution/metamorphism +! author Elizabeth C. Hunke, LANL + + subroutine write_restart_snow() + + use ice_fileunits, only: nu_dump_snow + use ice_state, only: trcrn + + ! local variables + + logical (kind=log_kind) :: diag + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw, k + character*3 ck + character(len=*),parameter :: subname='(write_restart_snow)' + + call icepack_query_tracer_indices(nt_smice_out=nt_smice, & + nt_smliq_out=nt_smliq, nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + diag = .true. + + !----------------------------------------------------------------- + + do k = 1,nslyr + write(ck,'(i3.3)') k + call write_restart_field(nu_dump_snow,0, trcrn(:,:,nt_smice+k-1,:,:), & + 'ruf8','smice'//trim(ck),ncat,diag) + call write_restart_field(nu_dump_snow,0, trcrn(:,:,nt_smliq+k-1,:,:), & + 'ruf8','smliq'//trim(ck),ncat,diag) + call write_restart_field(nu_dump_snow,0, trcrn(:,:,nt_rhos+k-1,:,:), & + 'ruf8','rhos'//trim(ck),ncat,diag) + call write_restart_field(nu_dump_snow,0, trcrn(:,:,nt_rsnw+k-1,:,:), & + 'ruf8','rsnw'//trim(ck),ncat,diag) + enddo + + end subroutine write_restart_snow + +!======================================================================= + +! Reads all values needed for a restart with snow redistribution/metamorphism +! author Elizabeth C. Hunke, LANL + + subroutine read_restart_snow() + + use ice_fileunits, only: nu_restart_snow + use ice_state, only: trcrn + + ! local variables + + logical (kind=log_kind) :: & + diag + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw, k + character*3 ck + character(len=*),parameter :: subname='(read_restart_snow)' + + call icepack_query_tracer_indices(nt_smice_out=nt_smice, & + nt_smliq_out=nt_smliq, nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + diag = .true. + + if (my_task == master_task) write(nu_diag,*) subname,'min/max snow tracers' + + do k=1,nslyr + write(ck,'(i3.3)') k + call read_restart_field(nu_restart_snow,0,trcrn(:,:,nt_smice+k-1,:,:), & + 'ruf8','smice'//trim(ck),ncat,diag, & + field_type=field_type_scalar,field_loc=field_loc_center) + call read_restart_field(nu_restart_snow,0,trcrn(:,:,nt_smliq+k-1,:,:), & + 'ruf8','smliq'//trim(ck),ncat,diag, & + field_type=field_type_scalar,field_loc=field_loc_center) + call read_restart_field(nu_restart_snow,0,trcrn(:,:,nt_rhos+k-1,:,:), & + 'ruf8','rhos'//trim(ck),ncat,diag, & + field_type=field_type_scalar,field_loc=field_loc_center) + call read_restart_field(nu_restart_snow,0,trcrn(:,:,nt_rsnw+k-1,:,:), & + 'ruf8','rsnw'//trim(ck),ncat,diag, & + field_type=field_type_scalar,field_loc=field_loc_center) + enddo + + end subroutine read_restart_snow + +!======================================================================= + ! Dumps all values needed for restarting ! author Elizabeth C. Hunke, LANL diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index f9631a91d..3dec72963 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -100,6 +100,8 @@ restart_pond_topo = .false. tr_pond_lvl = .true. restart_pond_lvl = .false. + tr_snow = .false. + restart_snow = .false. tr_iso = .false. restart_iso = .false. tr_aero = .false. @@ -197,6 +199,28 @@ pndaspect = 0.8 / +&snow_nml + snwredist = 'none' + snwgrain = .false. + use_smliq_pnd = .false. + rsnw_fall = 100.0 + rsnw_tmax = 1500.0 + rhosnew = 100.0 + rhosmin = 100.0 + rhosmax = 450.0 + windmin = 10.0 + drhosdwind = 27.3 + snwlvlfac = 0.3 + snw_aging_table = 'test' + snw_filename = 'unknown' + snw_rhos_fname = 'unknown' + snw_Tgrd_fname = 'unknown' + snw_T_fname = 'unknown' + snw_tau_fname = 'unknown' + snw_kappa_fname = 'unknown' + snw_drdt0_fname = 'unknown' +/ + &forcing_nml formdrag = .false. atmbndy = 'default' @@ -584,6 +608,21 @@ f_apeff_ai = 'm' / +&icefields_snow_nml + f_smassicen = 'x' + f_smassliqn = 'x' + f_rhos_cmpn = 'x' + f_rhos_cntn = 'x' + f_rsnwn = 'x' + f_smassice = 'm' + f_smassliq = 'm' + f_rhos_cmp = 'm' + f_rhos_cnt = 'm' + f_rsnw = 'm' + f_meltsliq = 'm' + f_fsloss = 'm' +/ + &icefields_bgc_nml f_fiso_atm = 'x' f_fiso_ocn = 'x' diff --git a/configuration/scripts/options/set_nml.snw30percent b/configuration/scripts/options/set_nml.snw30percent new file mode 100644 index 000000000..ecf88ad4e --- /dev/null +++ b/configuration/scripts/options/set_nml.snw30percent @@ -0,0 +1,5 @@ +tr_snow = .true. +snwredist = 'bulk' +snwlvlfac = 0.3 +nslyr = 5 + diff --git a/configuration/scripts/options/set_nml.snwITDrdg b/configuration/scripts/options/set_nml.snwITDrdg new file mode 100644 index 000000000..cddeedec3 --- /dev/null +++ b/configuration/scripts/options/set_nml.snwITDrdg @@ -0,0 +1,10 @@ +tr_snow = .true. +snwredist = 'ITDrdg' +nslyr = 5 +rhosnew = 100.0 +rhosmin = 100.0 +rhosmax = 450.0 +windmin = 10.0 +drhosdwind = 27.3 +snwlvlfac = 0.3 + diff --git a/configuration/scripts/options/set_nml.snwgrain b/configuration/scripts/options/set_nml.snwgrain new file mode 100644 index 000000000..653030385 --- /dev/null +++ b/configuration/scripts/options/set_nml.snwgrain @@ -0,0 +1,15 @@ +tr_snow = .true. +snwgrain = .true. +use_smliq_pnd = .true. +rsnw_fall = 54.526 +rsnw_tmax = 1500.0 +snw_aging_table = 'file' +snw_filename = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/snicar_drdt_bst_fit_60_c04262019.nc' +snw_tau_fname = 'snowEmpiricalGrowthParameterTau' +snw_kappa_fname = 'snowEmpiricalGrowthParameterKappa' +snw_drdt0_fname = 'snowPropertyRate' +snw_rhos_fname = 'nGrainAgingSnowDensity' +snw_Tgrd_fname = 'nGrainAgingTempGradient' +snw_T_fname = 'nGrainAgingTemperature' +nslyr = 5 + diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 78df9ac81..7aac29450 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -58,6 +58,9 @@ restart gx3 4x2 fsd12,debug,short smoke gx3 8x2 fsd12ww3,diag24,run1day smoke gx3 4x1 isotope,debug restart gx3 8x2 isotope +smoke gx3 4x1 snwITDrdg,snwgrain,icdefault,debug +smoke gx3 4x1 snw30percent,icdefault,debug +restart gx3 8x2 snwITDrdg,icdefault,snwgrain restart gx3 4x4 gx3ncarbulk,iobinary restart gx3 4x4 histall,precision8,cdf64 smoke gx3 30x1 bgcz,histall diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 2efcd0335..0a04b5e26 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -168,6 +168,7 @@ either Celsius or Kelvin units). "dms", "dimethyl sulfide concentration", "mmol/m\ :math:`^3`" "dmsp", "dimethyl sulfoniopropionate concentration", "mmol/m\ :math:`^3`" "dpscale", "time scale for flushing in permeable ice", ":math:`1\times 10^{-3}`" + "drhosdwind", "wind compaction factor for snow", "27.3 kg s/m\ :math:`^{4}`" "dragio", "drag coefficient for water on ice", "0.00536" "dSdt_slow_mode", "drainage strength parameter", "" "dsnow", "change in snow thickness", "m" @@ -256,6 +257,7 @@ either Celsius or Kelvin units). "fsnow", "snowfall rate", "kg/m\ :math:`^2`/s" "fsnowrdg", "snow fraction that survives in ridging", "0.5" "fsurf(n)(_f)", "net surface heat flux excluding fcondtop", "W/m\ :math:`^2`" + "fsloss", "rate of snow loss to leads", "kg/m\ :math:`^{2}` s" "fsw", "incoming shortwave radiation", "W/m\ :math:`^2`" "fswabs", "total absorbed shortwave radiation", "W/m\ :math:`^2`" "fswfac", "scaling factor to adjust ice quantities for updated data", "" @@ -393,6 +395,8 @@ either Celsius or Kelvin units). "meltb", "basal ice melt", "m" "meltl", "lateral ice melt", "m" "melts", "snow melt", "m" + "meltsliq", "snow melt mass", "kg/m\ :math:`^{2}`" + "meltsliqn", "snow melt mass in category n", "kg/m\ :math:`^{2}`" "meltt", "top ice melt", "m" "min_salin", "threshold for brine pockets", "0.1 ppt" "mlt_onset", "day of year that surface melt begins", "" @@ -556,14 +560,21 @@ either Celsius or Kelvin units). "rhofresh", "density of fresh water", "1000.0 kg/m\ :math:`^3`" "rhoi", "density of ice", "917. kg/m\ :math:`^3`" "rhos", "density of snow", "330. kg/m\ :math:`^3`" + "rhos_cmp", "density of snow due to wind compaction", "kg/m\ :math:`^3`" + "rhos_cnt", "density of ice and liquid content of snow", "kg/m\ :math:`^3`" "rhosi", "average sea ice density (for hbrine tracer)", "940. kg/m\ :math:`^3`" + "rhosmax", "maximum snow density", "450 kg/m\ :math:`^{3}`" + "rhosmin", "minimum snow density", "100 kg/m\ :math:`^{3}`" + "rhosnew", "new snow density", "100 kg/m\ :math:`^{3}`" "rhow", "density of seawater", "1026. kg/m\ :math:`^3`" "rnilyr", "real(nlyr)", "" "rside", "fraction of ice that melts laterally", "" - "rsnw_fresh", "freshly fallen snow grain radius", "100. :math:`\times` 10\ :math:`^{-6}` m" + "rsnw", "snow grain radius", "10\ :math:`^{-6}` m" + "rsnw_fall", "freshly fallen snow grain radius", "100. :math:`\times` 10\ :math:`^{-6}` m" "rsnw_melt", "melting snow grain radius", "1000. :math:`\times` 10\ :math:`^{-6}` m" "rsnw_nonmelt", "nonmelting snow grain radius", "500. :math:`\times` 10\ :math:`^{-6}` m" "rsnw_sig", "standard deviation of snow grain radius", "250. :math:`\times` 10\ :math:`^{-6}` m" + "rsnw_tmax", "maximum snow radius", "1500. :math:`\times` 10\ :math:`^{-6}` m" "runid", "identifier for run", "" "runtype", "type of initialization used", "" "**S**", "", "" @@ -586,6 +597,25 @@ either Celsius or Kelvin units). "snoice", "snow–ice formation", "m" "snowpatch", "length scale for parameterizing nonuniform snow coverage", "0.02 m" "skl_bgc", "biogeochemistry on/off", "" + "smassice", "mass of ice in snow from smice tracer", "kg/m\ :math:`^2`" + "smassliq", "mass of liquid in snow from smliq tracer", "kg/m\ :math:`^2`" + "snowage_drdt0", "initial rate of change of effective snow radius", " " + "snowage_rhos", "snow aging parameter (density)", " " + "snowage_kappa", "snow aging best-fit parameter", " " + "snowage_tau", "snow aging best-fit parameter", " " + "snowage_T", "snow aging parameter (temperature)", " " + "snowage_Tgrd", "snow aging parameter (temperature gradient)", " " + "snw_aging_table", "snow aging lookup table", " " + "snw_filename", "snowtable filename", " " + "snw_tau_fname", "snowtable file tau fieldname", " " + "snw_kappa_fname", "snowtable file kappa fieldname", " " + "snw_drdt0_fname", "snowtable file drdt0 fieldname", " " + "snw_rhos_fname", "snowtable file rhos fieldname", " " + "snw_Tgrd_fname", "snowtable file Tgrd fieldname", " " + "snw_T_fname", "snowtable file T fieldname", " " + "snwgrain", "activate snow metamorphosis", " " + "snwlvlfac", "fractional increase in snow depth for redistribution on ridges", "0.3" + "snwredist", "type of snow redistribution", " " "spval", "special value (single precision)", ":math:`10^{30}`", "" "spval_dbl", "special value (double precision)", ":math:`10^{30}`", "" "ss_tltx(y)", "sea surface in the x(y) direction", "m/m" @@ -666,6 +696,7 @@ either Celsius or Kelvin units). "update_ocn_f", "if true, include frazil ice fluxes in ocean flux fields", "" "use_leap_years", "if true, include leap days", "" "use_restart_time", "if true, use date from restart file", "" + "use_smliq_pnd", "use liquid in snow for ponds", " " "ustar_min", "minimum friction velocity under ice", "" "ucstr", "string identifying U grid for history variables", "" "uvel", "x-component of ice velocity", "m/s" @@ -691,6 +722,7 @@ either Celsius or Kelvin units). "wave_spectrum", "wave spectrum", "m\ :math:`^2`/s" "wavefreq", "wave frequencies", "1/s" "wind", "wind speed", "m/s" + "windmin", "minimum wind speed to compact snow", "10 m/s" "write_history", "if true, write history now", "" "write_ic", "if true, write initial conditions", "" "write_restart", "if 1, write restart now", "" diff --git a/doc/source/developer_guide/dg_driver.rst b/doc/source/developer_guide/dg_driver.rst index a10cb319a..637e91b68 100644 --- a/doc/source/developer_guide/dg_driver.rst +++ b/doc/source/developer_guide/dg_driver.rst @@ -65,7 +65,6 @@ The initialize calling sequence looks something like:: call init_thermo_vertical ! initialize vertical thermodynamics call icepack_init_itd(ncat, hin_max) ! ice thickness distribution if (tr_fsd) call icepack_init_fsd_bounds ! floe size distribution - call calendar(time) ! determine the initial date call init_forcing_ocn(dt) ! initialize sss and sst from data call init_state ! initialize the ice state call init_transport ! initialize horizontal transport @@ -74,10 +73,13 @@ The initialize calling sequence looks something like:: call init_diags ! initialize diagnostic output points call init_history_therm ! initialize thermo history variables call init_history_dyn ! initialize dynamic history variables + call calc_timesteps ! update timestep counter if not using npt_unit="1" call init_shortwave ! initialize radiative transfer + call advance_timestep ! advance the time step call init_forcing_atmo ! initialize atmospheric forcing (standalone) if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice call get_forcing* ! read forcing data (standalone) + if (tr_snow) call icepack_init_snow ! advanced snow physics See a **CICE_InitMod.F90** file for the latest. @@ -105,6 +107,13 @@ The run sequence within a time loop looks something like:: call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset) enddo + if (tr_snow) then ! advanced snow physics + do iblk = 1, nblocks + call step_snow (dt, iblk) + enddo + call update_state (dt) ! clean up + endif + do iblk = 1, nblocks call step_radiation (dt, iblk) call coupling_prep (iblk) diff --git a/doc/source/science_guide/sg_tracers.rst b/doc/source/science_guide/sg_tracers.rst index bbd18eb1f..215c13d08 100644 --- a/doc/source/science_guide/sg_tracers.rst +++ b/doc/source/science_guide/sg_tracers.rst @@ -90,6 +90,10 @@ is not in use. "tr_iso", "n_iso", "vice, vsno", "nt_iso"," " "tr_brine", " ", "vice", "nt_fbri", " " "tr_fsd","nfsd","aice","nt_fsd"," " + "tr_snow","nslyr","vsno","nt_rsnw"," " + " ","nslyr","vsno","nt_rhos"," " + " ","nslyr","vsno","nt_smice"," " + " ","nslyr","vsno","nt_smliq"," " "solve_zsal", "n_trzs", "fbri or (a,v)ice", "nt_bgc_S", " " "tr_bgc_N", "n_algae", "fbri or (a,v)ice", "nt_bgc_N", "nlt_bgc_N" "tr_bgc_Nit", " ", "fbri or (a,v)ice", "nt_bgc_Nit", "nlt_bgc_Nit" @@ -115,4 +119,4 @@ Users may add any number of additional tracers that are transported conservative provided that the dependency ``trcr_depend`` is defined appropriately. See Section :ref:`addtrcr` for guidance on adding tracers. -Please see the `Icepack documentation `_ for additional information about tracers that depend on other tracers, the floe size distribution, age of the ice, aerosols, water isotopes, brine height, and the sea ice ecosystem. +Please see the `Icepack documentation `_ for additional information about tracers that depend on other tracers, the floe size distribution, advanced snow physics, age of the ice, aerosols, water isotopes, brine height, and the sea ice ecosystem. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index b1f6b4614..65a50120b 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -327,6 +327,7 @@ tracer_nml "``tr_pond_cesm``", "logical", "CESM melt ponds", "``.false.``" "``tr_pond_lvl``", "logical", "level-ice melt ponds", "``.false.``" "``tr_pond_topo``", "logical", "topo melt ponds", "``.false.``" + "``tr_snow``", "logical", "advanced snow physics", "``.false.``" "``restart_aero``", "logical", "restart tracer values from file", "``.false.``" "``restart_age``", "logical", "restart tracer values from file", "``.false.``" "``restart_fsd``", "logical", "restart floe size distribution values from file", "``.false.``" @@ -485,6 +486,40 @@ ponds_nml "``rfracmin``", ":math:`0 \le r_{min} \le 1`", "minimum melt water added to ponds", "0.15" "", "", "", "" +snow_nml +~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. csv-table:: **snow_nml namelist options** + :header: "variable", "options/format", "description", "default value" + :widths: 15, 15, 30, 15 + + "", "", "", "" + "``drhosdwind``", "real", "wind compactions factor for now in kg-s/m^4", "27.3" + "``rhosmax``", "real", "maximum snow density in kg/m^3", "450." + "``rhosmin``", "real", "minimum snow density in kg/m^3", "100." + "``rhosnew``", "real", "new snow density in kg/m^3", "100." + "``rsnw_fall``", "real", "radius of new snow in 1.0e-6 m", "100." + "``rsnw_tmax``", "real", "maximum snow radius in 1.0e-6 m", "1500." + "``snwgrain``", "logical", "snow metamorophsis flag", "``.false.``" + "``snwlvlfac``", "real", "fractional increase in snow", "0.3" + "``snwredist``", "``bulk``", "bulk snow redistribution scheme", "``none``" + "", "``ITD``", "ITD snow redistribution scheme", "" + "", "``ITDrdg``", "ITDrdg snow redistribution scheme", "" + "", "``none``", "snow redistribution scheme off", "" + "``snw_aging_table``", "file", "read 1D and 3D fields for dry metamorophsis lookup table", "test" + "", "snicar", "read 3D fields for dry metamorophsis lookup table", "" + "", "test", "internally generated dry metamorophsis lookup table for testing", "" + "``snw_drdt0_fname``", "string", "snow aging file drdt0 fieldname", "unknown" + "``snw_filename``", "string", "snow aging table data filename", "unknown" + "``snw_kappa_fname``", "string", "snow aging file kappa fieldname", "unknown" + "``snw_rhos_fname``", "string", "snow aging file rhos fieldname", "unknown" + "``snw_T_fname``", "string", "snow aging file T fieldname", "unknown" + "``snw_tau_fname``", "string", "snow aging file tau fieldname", "unknown" + "``snw_Tgrd_fname``", "string", "snow aging file Tgrd fieldname", "unknown" + "``use_smliq_pnd``", "logical", "use liquid in snow for ponds", "``.false.``" + "``windmin``", "real", "minimum wind speed to compact snow in m/s", "10." + "", "", "", "" + forcing_nml ~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/icepack b/icepack index a80472b54..53ffce04d 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit a80472b547aa6d7a85f8ae5e1449273a323e0371 +Subproject commit 53ffce04d729a3e802b397340c4410126b5b0ac9