diff --git a/bld/check_decomp.csh b/bld/check_decomp.csh index 1331efc9..149c62f7 100755 --- a/bld/check_decomp.csh +++ b/bld/check_decomp.csh @@ -1,6 +1,6 @@ #!/bin/csh -set res = (gx3v7 gx1v6 tx0.1v2 ne240np4) +set res = (gx3v7 gx1v6 tx0.1v2 ne240np4 tnx1v1 tnx0.25v1 tnx1v3 tnx1v4) set pe1 = (1 1 1 1) set pex = (500 4000 60000 60000) diff --git a/cime_config/buildnml b/cime_config/buildnml index 8c66a08c..cedc755d 100755 --- a/cime_config/buildnml +++ b/cime_config/buildnml @@ -47,8 +47,10 @@ def _create_namelists(case, confdir, infile, nmlgen): cice_mode = case.get_value("CICE_MODE") hgrid = case.get_value('ICE_GRID') + ocn_model = case.get_value('COMP_OCN') config['hgrid'] = hgrid config['cice_mode'] = cice_mode + config['ocn_model'] = ocn_model cice_config_opts = case.get_value('CICE_CONFIG_OPTS') if "cice5" in cice_config_opts: @@ -149,6 +151,7 @@ def _create_namelists(case, confdir, infile, nmlgen): 'dynamics_nml', 'shortwave_nml', 'ponds_nml', + 'snowphys_nml', 'forcing_nml', 'domain_nml', 'zbgc_nml', diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index f6e62ecb..ae693f32 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -5,9 +5,10 @@ - Sea ICE (cice) model version 5 + Sea ICE (cice) model version 5 :prescribed cice :with modifications appropriate for CMIP6 experiments + :with NORESM modifications appropriate for CMIP6 experiments @@ -50,17 +51,6 @@ - - char - - - cam4=.true. - - run_component_cice - env_run.xml - CICE specific namelist settings for -namelist option - - logical true,false @@ -156,6 +146,8 @@ $SRCROOT/components/cice/cime_config/usermods_dirs/cmip6 + $SRCROOT/components/cice/cime_config/usermods_dirs/noresm-cmip6 run_component_cice env_case.xml diff --git a/cime_config/namelist_definition_cice.xml b/cime_config/namelist_definition_cice.xml index c068cc87..47fa9bcb 100644 --- a/cime_config/namelist_definition_cice.xml +++ b/cime_config/namelist_definition_cice.xml @@ -377,6 +377,13 @@ bin nc + bin + bin + bin + bin + bin + bin + bin @@ -396,6 +403,13 @@ tripole regional regional + tripole + tripole + tripole + tripole + tripole + tripole + tripole @@ -418,6 +432,13 @@ $DIN_LOC_ROOT/ice/cice/ar9v3_20101118.grid $DIN_LOC_ROOT/ice/cice/ar9v3_20101118.grid $DIN_LOC_ROOT/ice/cice/domain.ocn.us20.20110426.nc + $DIN_LOC_ROOT/ocn/pop/tnx2v1/grid/horiz_grid_20130206.ieeer8 + $DIN_LOC_ROOT/ocn/pop/tnx1v1/grid/horiz_grid_20120120.ieeer8 + $DIN_LOC_ROOT/ocn/micom/tnx1v3/20170612/horiz_grid_20170612.ieeer8 + $DIN_LOC_ROOT/ocn/micom/tnx1v4/20170601/horiz_grid_20170601.ieeer8 + $DIN_LOC_ROOT/ocn/pop/tnx0.25v1/grid/horiz_grid_20130930.ieeer8 + $DIN_LOC_ROOT/ocn/micom/tnx0.25v3/20170623/horiz_grid_20170623.ieeer8 + $DIN_LOC_ROOT/ocn/micom/tnx0.25v4/20170619/horiz_grid_20170619.ieeer8 @@ -440,6 +461,13 @@ $DIN_LOC_ROOT/ocn/pop/ar9v3/grid/kmt.ar9v3.ocn.ChanBarrier.20100622.ieeei4 $DIN_LOC_ROOT/ice/cice/ar9v3_20110106.kmt $DIN_LOC_ROOT/ice/cice/domain.ocn.us20.20110426.nc + $DIN_LOC_ROOT/ocn/pop/tnx2v1/grid/topography_20130206.ieeei4 + $DIN_LOC_ROOT/ocn/pop/tnx1v1/grid/topography_20120120.ieeei4 + $DIN_LOC_ROOT/ocn/micom/tnx1v3/20170612/topography_20170612.ieeei4 + $DIN_LOC_ROOT/ocn/micom/tnx1v4/20170601/topography_20170601.ieeei4 + $DIN_LOC_ROOT/ocn/pop/tnx0.25v1/grid/topography_20130930.ieeei4 + $DIN_LOC_ROOT/ocn/micom/tnx0.25v3/20170623/topography_20170623.ieeei4 + $DIN_LOC_ROOT/ocn/micom/tnx0.25v4/20170619/topography_20170619.ieeei4 @@ -539,6 +567,13 @@ tripole tripole tripole + tripole + tripole + tripole + tripole + tripole + tripole + tripole @@ -1337,6 +1372,47 @@ + + + + + + char + snowphys + snowphys_nml + + Blowing snow parameterization (none, lecomte2013) + + + lecomte2013 + + + + + real + snowphys + snowphys_nml + + Constant snow density (real) + + + 330.0 + + + + + real + snowphys + snowphys_nml + + Snow thermal conductivity (real) + + + 0.3 + + + + @@ -1451,17 +1527,18 @@ - - - - - - - - - - - + + logical + forcing + forcing_nml + + include frazil water/salt fluxes in ocean fluxes + + + .false. + .true. + + @@ -2935,7 +3012,24 @@ - + + + + + char + history + icefields_nml + + + + mxxxx + mdhxx + xxxxx + + + + + char history icefields_nml @@ -2948,6 +3042,8 @@ + + char history diff --git a/cime_config/usermods_dirs/noresm-cmip6/user_nl_cice b/cime_config/usermods_dirs/noresm-cmip6/user_nl_cice new file mode 100644 index 00000000..43a2019c --- /dev/null +++ b/cime_config/usermods_dirs/noresm-cmip6/user_nl_cice @@ -0,0 +1,30 @@ +!---------------------------------------------------------------------------------- +! Users should add all user specific namelist changes below in the form of +! namelist_var = new_namelist_value +! Note - that it does not matter what namelist group the namelist_var belongs to +!---------------------------------------------------------------------------------- + +histfreq = 'm','d','x','x','x' +histfreq_n = 1,1,1,1,1 +f_CMIP = 'mdxxx' +f_hi ="mxxxx" +f_hs="mxxxx" +f_fswdn="mxxxx" +f_fswabs="mxxxx" +f_congel="mxxxx" +f_frazil="mxxxx" +f_meltt="mxxxx" +f_melts="mxxxx" +f_meltb="mxxxx" +f_meltl="mxxxx" +f_fswthru="mxxxx" +f_dvidtt="mxxxx" +f_dvidtd="mxxxx" +f_daidtt="mxxxx" +f_daidtd="mxxxx" +f_apond_ai="mxxxx" +f_hpond_ai="mxxxx" +f_apeff_ai="mxxxx" +f_snowfrac="mxxxx" +f_aicen="mxxxx" +f_snowfracn="mxxxx" diff --git a/src/drivers/cesm/ice_constants.F90 b/src/drivers/cesm/ice_constants.F90 index 06584c53..0ca2d319 100644 --- a/src/drivers/cesm/ice_constants.F90 +++ b/src/drivers/cesm/ice_constants.F90 @@ -31,12 +31,20 @@ module ice_constants public +!jd + !----------------------------------------------------------------- + ! physical variables, should not be declared as constants + !----------------------------------------------------------------- + real (kind=dbl_kind) :: & + rhos = 330.0_dbl_kind ,&! density of snow (kg/m^3) + ksno = 0.30_dbl_kind ! thermal conductivity of snow (W/m/deg) +!jd !----------------------------------------------------------------- ! physical constants !----------------------------------------------------------------- real (kind=dbl_kind), parameter :: & - rhos = 330.0_dbl_kind ,&! density of snow (kg/m^3) +!jd rhos = 330.0_dbl_kind ,&! density of snow (kg/m^3) rhoi = SHR_CONST_RHOICE ,&! density of ice (kg/m^3) rhow = SHR_CONST_RHOSW ,&! density of seawater (kg/m^3) cp_air = SHR_CONST_CPDAIR ,&! specific heat of air (J/kg/K) @@ -87,7 +95,7 @@ module ice_constants kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg) kseaice= 2.00_dbl_kind ,&! thermal conductivity of sea ice (W/m/deg) ! (used in zero layer thermodynamics option) - ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) +!jd ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg) zref = 10._dbl_kind ,&! reference height for stability (m) ! hs0 = 0.03_dbl_kind, &! parameter for delta-Eddington snow frac ! hsmin = 0.0001_dbl_kind, &! minimum snow thickness for dEdd diff --git a/src/drivers/cesm/ice_prescribed_mod.F90 b/src/drivers/cesm/ice_prescribed_mod.F90 index 97b43ae6..f81bc0be 100644 --- a/src/drivers/cesm/ice_prescribed_mod.F90 +++ b/src/drivers/cesm/ice_prescribed_mod.F90 @@ -92,11 +92,13 @@ module ice_prescribed_mod real (kind=dbl_kind), parameter :: & cp_sno = 0.0_dbl_kind & ! specific heat of snow (J/kg/K) , rLfi = Lfresh*rhoi & ! latent heat of fusion ice (J/m^3) - , rLfs = Lfresh*rhos & ! latent heat of fusion snow (J/m^3) +!jd These are not used in the code and gives compilation error when changing r +!jd rhos to a namelist variable (not parameter) +!jd , rLfs = Lfresh*rhos & ! latent heat of fusion snow (J/m^3) , rLvi = Lvap*rhoi & ! latent heat of vapor*rhoice (J/m^3) - , rLvs = Lvap*rhos & ! latent heat of vapor*rhosno (J/m^3) +!jd , rLvs = Lvap*rhos & ! latent heat of vapor*rhosno (J/m^3) , rcpi = cp_ice*rhoi & ! heat capacity of fresh ice (J/m^3) - , rcps = cp_sno*rhos & ! heat capacity of snow (J/m^3) +!jd , rcps = cp_sno*rhos & ! heat capacity of snow (J/m^3) , rcpidepressT = rcpi*depressT & ! param for finding T(z) from q (J/m^3) , rLfidepressT = rLfi*depressT ! param for heat capacity (J deg/m^3) ! heat capacity of sea ice, rhoi*C=rcpi+rLfidepressT*salinity/T^2 diff --git a/src/source/ice_flux.F90 b/src/source/ice_flux.F90 index 8414e628..36811258 100644 --- a/src/source/ice_flux.F90 +++ b/src/source/ice_flux.F90 @@ -66,6 +66,9 @@ module ice_flux strinty , & ! divergence of internal ice stress, y (N/m^2) daidtd , & ! ice area tendency due to transport (1/s) dvidtd , & ! ice volume tendency due to transport (m/s) +!jd + dvsdtd , & ! snow volume tendency due to transport (m/s) +!jd dagedtd , & ! ice age tendency due to transport (s/s) dardg1dt, & ! rate of area loss by ridging ice (1/s) dardg2dt, & ! rate of area gain by new ridges (1/s) @@ -630,6 +633,9 @@ subroutine init_history_therm Cdn_ocn_floe, Cdn_ocn_skin, formdrag use ice_state, only: aice, vice, trcr, tr_iage, nt_iage use ice_constants, only: vonkar,zref,iceruf +!jd + use ice_snowphys, only: snow2ocn, snowfonice +!jd fsurf (:,:,:) = c0 fcondtop(:,:,:)= c0 @@ -638,6 +644,10 @@ subroutine init_history_therm frazil (:,:,:) = c0 frazil_diag (:,:,:) = c0 snoice (:,:,:) = c0 +!jd + snow2ocn (:,:,:) = c0 + snowfonice(:,:,:) = c0 +!jd dsnow (:,:,:) = c0 meltt (:,:,:) = c0 melts (:,:,:) = c0 @@ -704,7 +714,8 @@ end subroutine init_history_therm subroutine init_history_dyn - use ice_state, only: aice, vice, trcr, tr_iage, nt_iage +!jd use ice_state, only: aice, vice, trcr, tr_iage, nt_iage + use ice_state, only: aice, vice, vsno, trcr, tr_iage, nt_iage sig1 (:,:,:) = c0 sig2 (:,:,:) = c0 @@ -722,6 +733,9 @@ subroutine init_history_dyn opening (:,:,:) = c0 daidtd (:,:,:) = aice(:,:,:) ! temporary initial area dvidtd (:,:,:) = vice(:,:,:) ! temporary initial volume +!jd + dvsdtd (:,:,:) = vsno(:,:,:) ! temporary initial volume +!jd if (tr_iage) & dagedtd (:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age fm (:,:,:) = c0 @@ -777,6 +791,10 @@ subroutine merge_fluxes (nx_block, ny_block, & meltt, melts, & meltb, & congel, snoice, & +!jd + snow2ocnn, snow2ocn, & + snowfonicen,snowfonice, & +!jd Uref, Urefn, & Qref_iso, Qrefn_iso, & fiso_evap,fiso_evapn, & @@ -821,6 +839,10 @@ subroutine merge_fluxes (nx_block, ny_block, & meltbn , & ! bottom ice melt (m) meltsn , & ! snow melt (m) congeln , & ! congelation ice growth (m) +!jd + snow2ocnn, &! snow blowing of ice + snowfonicen,&! fraction of snow kept on ice +!jd snoicen ! snow-ice growth (m) real (kind=dbl_kind), dimension(nx_block,ny_block), optional, intent(in):: & @@ -859,6 +881,10 @@ subroutine merge_fluxes (nx_block, ny_block, & meltb , & ! bottom ice melt (m) melts , & ! snow melt (m) congel , & ! congelation ice growth (m) +!jd + snow2ocn, &! snow blowing of ice + snowfonice,&! fraction of snow kept on ice +!jd snoice ! snow-ice growth (m) real (kind=dbl_kind), dimension(nx_block,ny_block), optional, & @@ -937,6 +963,9 @@ subroutine merge_fluxes (nx_block, ny_block, & congel (i,j) = congel (i,j) + congeln (i,j)*aicen(i,j) snoice (i,j) = snoice (i,j) + snoicen (i,j)*aicen(i,j) + ! Blowing snow + snow2ocn (i,j) = snow2ocn (i,j) + snow2ocnn(i,j)*aicen(i,j) + snowfonice(i,j)= snowfonice(i,j)+snowfonicen(i,j)*aicen(i,j) enddo ! ij end subroutine merge_fluxes diff --git a/src/source/ice_history.F90 b/src/source/ice_history.F90 index b0f98537..b55e58cd 100644 --- a/src/source/ice_history.F90 +++ b/src/source/ice_history.F90 @@ -184,6 +184,17 @@ subroutine init_hist (dt) f_siv = 'mxxxx' f_sidmasstranx = 'mxxxx' f_sidmasstrany = 'mxxxx' +!jd noresm-start + f_simasstranx = 'mxxxx' + f_simasstrany = 'mxxxx' + f_sndmasstranx = 'mxxxx' + f_sndmasstrany = 'mxxxx' + f_aicetranx = 'mxxxx' + f_aicetrany = 'mxxxx' + f_simass = 'mxxxx' + f_sisnmass = 'mxxxx' + f_sisal = 'mxxxx' +!jd end f_sistrxdtop = 'mxxxx' f_sistrydtop = 'mxxxx' f_sistrxubot = 'mxxxx' @@ -205,7 +216,10 @@ subroutine init_hist (dt) f_sidmassmeltbot = 'mxxxx' f_sidmasslat = 'mxxxx' f_sndmasssnf = 'mxxxx' + f_sndmasssi = 'mxxxx' f_sndmassmelt = 'mxxxx' + f_sndmasswindrif = 'mxxxx'; f_sn2oc='xxxxx' + f_sndmassdyn = 'mxxxx' f_siflswdtop = 'mxxxx' f_siflswutop = 'mxxxx' f_siflswdbot = 'mxxxx' @@ -231,7 +245,9 @@ subroutine init_hist (dt) f_sistreave = 'mxxxx' f_sistremax = 'mxxxx' f_sirdgthick = 'mxxxx' - f_siitdconc = 'mxxxx' + f_sisal = 'mxxxx' +!jd f_siitdconc = 'mxxxx' + f_siitdconc = 'xxxxx' ! same as aicen f_siitdthick = 'mxxxx' f_siitdsnthick = 'mxxxx' @@ -239,7 +255,7 @@ subroutine init_hist (dt) f_iage = 'xxxxx' f_snoice = 'xxxxx' f_strength = 'xxxxx' - f_divu = 'xxxxx' +!jd f_divu = 'xxxxx' f_flat = 'xxxxx' f_flat_ai = 'xxxxx' f_flwdn = 'xxxxx' @@ -261,7 +277,7 @@ subroutine init_hist (dt) f_strocny = 'xxxxx' f_rain = 'xxxxx' f_snow = 'xxxxx' - f_shear = 'xxxxx' +!jd f_shear = 'xxxxx' f_Tsfc = 'xxxxx' f_uvel = 'xxxxx' f_vvel = 'xxxxx' @@ -371,6 +387,10 @@ subroutine init_hist (dt) call broadcast_scalar (f_frazil, master_task) call broadcast_scalar (f_snoice, master_task) call broadcast_scalar (f_dsnow, master_task) +!jd + call broadcast_scalar (f_sn2oc, master_task) + call broadcast_scalar (f_snfonice, master_task) +!jd call broadcast_scalar (f_meltt, master_task) call broadcast_scalar (f_melts, master_task) call broadcast_scalar (f_meltb, master_task) @@ -421,6 +441,17 @@ subroutine init_hist (dt) call broadcast_scalar (f_siv, master_task) call broadcast_scalar (f_sidmasstranx, master_task) call broadcast_scalar (f_sidmasstrany, master_task) +!jd + call broadcast_scalar (f_simasstranx, master_task) + call broadcast_scalar (f_simasstrany, master_task) + call broadcast_scalar (f_sndmasstranx, master_task) + call broadcast_scalar (f_sndmasstrany, master_task) + call broadcast_scalar (f_aicetranx, master_task) + call broadcast_scalar (f_aicetrany, master_task) + call broadcast_scalar (f_simass, master_task) + call broadcast_scalar (f_sisnmass, master_task) + call broadcast_scalar (f_sisal, master_task) +!jd call broadcast_scalar (f_sistrxdtop, master_task) call broadcast_scalar (f_sistrydtop, master_task) call broadcast_scalar (f_sistrxubot, master_task) @@ -443,6 +474,11 @@ subroutine init_hist (dt) call broadcast_scalar (f_sidmassmeltbot, master_task) call broadcast_scalar (f_sidmasslat, master_task) call broadcast_scalar (f_sndmasssnf, master_task) +!jd noresm + call broadcast_scalar (f_sndmasssi, master_task) + call broadcast_scalar (f_sndmasswindrif, master_task) + call broadcast_scalar (f_sndmassdyn, master_task) +!jd end call broadcast_scalar (f_sndmassmelt, master_task) call broadcast_scalar (f_siflswdtop, master_task) call broadcast_scalar (f_siflswutop, master_task) @@ -788,6 +824,18 @@ subroutine init_hist (dt) "snow formation", & "none", mps_to_cmpdy/dt, c0, & ns1, f_dsnow) + +!jd + call define_hist_field(n_sn2oc,"snow2ocn","kg/m2/s",tstr2D, tcstr, & + "Snow blowing into ocean", & + "none", c1, c0, & + ns1, f_sn2oc) + + call define_hist_field(n_snfonice,"snfonice","1",tstr2D, tcstr, & + "Fraction of snow remaining on ice during preciptation and wind", & + "none", c1, c0, & + ns1, f_snfonice) +!jd call define_hist_field(n_meltt,"meltt","cm/day",tstr2D, tcstr, & "top ice melt", & @@ -1121,6 +1169,65 @@ subroutine init_hist (dt) "y component of snow and sea ice mass transport", & "none", c1, c0, & ns1, f_sidmasstrany) + + +!jd noresm-start +!jd f_simasstranx = 'x', + call define_hist_field(n_simasstranx,"simasstranx","kg/s",ustr2D, ucstr, & + "x component of sea ice mass transport", & + "none", c1, c0, & + ns1, f_simasstranx) + +!jd f_simasstrany = 'x', & + call define_hist_field(n_simasstrany,"simasstrany","kg/s",ustr2D, ucstr, & + "y component of sea ice mass transport", & + "none", c1, c0, & + ns1, f_simasstrany) + + +!jd f_sndmasstranx = 'x', + call define_hist_field(n_sndmasstranx,"sndmasstranx","kg/s",ustr2D, ucstr, & + "x component of snow mass transport", & + "none", c1, c0, & + ns1, f_sndmasstranx) + +!jd f_sndmasstrany = 'x', & + call define_hist_field(n_sndmasstrany,"sndmasstrany","kg/s",ustr2D, ucstr, & + "y component of snow mass transport", & + "none", c1, c0, & + ns1, f_sndmasstrany) + +!jd f_aicetranx = 'x', + call define_hist_field(n_aicetranx,"aicetranx","m^2/s",ustr2D, ucstr, & + "x component of ice area transport", & + "none", c1, c0, & + ns1, f_aicetranx) + +!jd f_aicetrany = 'x', & + call define_hist_field(n_aicetrany,"aicetrany","m^2/s",ustr2D, ucstr, & + "y component of ice area transport", & + "none", c1, c0, & + ns1, f_aicetrany) + + + call define_hist_field(n_simass,"simass","kg/m^2",tstr2D, tcstr, & + "Sea ice mass per Area", & + "none", c1, c0, & + ns1, f_simass) + + call define_hist_field(n_sisnmass,"sisnmass","kg/m^2",tstr2D, tcstr, & + "Snow Mass per Area", & + "none", c1, c0, & + ns1, f_sisnmass) + + call define_hist_field(n_sisal,"sisal","0.001",tstr2D, tcstr, & + "Mean Sea Ice Salinity of all ice (used in thermodynamics)", & + "none", c1, c0, & + ns1, f_sisal) + + +!jd end + call define_hist_field(n_sistrxdtop,"sistrxdtop","N m-2",ustr2D, ucstr, & "x component of atmospheric stress on sea ice", & @@ -1231,7 +1338,25 @@ subroutine init_hist (dt) "snow mass change from snow fall", & "none", c1, c0, & ns1, f_sndmasssnf) - + +!jd start + call define_hist_field(n_sndmasssi,"sndmasssi","kg m-2 s-1",tstr2D, tcstr, & + "Snow Mass Rate of Change Through Snow-to-Ice Conversion", & + "none", c1, c0, & + ns1, f_sndmasssi) + + call define_hist_field(n_sndmasswindrif,"sndmasswindrif","kg m-2 s-1",tstr2D, tcstr, & + "Snow mass change through wind drift of snow", & + "none", c1, c0, & + ns1, f_sndmasswindrif) + + call define_hist_field(n_sndmassdyn,"sndmassdyn","kg m-2 s-1",tstr2D, tcstr, & + "Snow mass rate of change due to dynamics", & + "none", c1, c0, & + ns1, f_sndmassdyn) + + +!jd end call define_hist_field(n_sndmassmelt,"sndmassmelt","kg m-2 s-1",tstr2D, tcstr, & "snow mass change from snow melt", & "none", c1, c0, & @@ -1653,7 +1778,7 @@ subroutine accum_hist (dt) awtvdr, awtidr, awtvdf, awtidf, Lfresh, rhoi, rhos, rhow, cp_ice, spval_dbl, hs_min, & p001, ice_ref_salinity use ice_domain, only: blocks_ice, nblocks - use ice_grid, only: tmask, lmask_n, lmask_s, tarea, dxu, dyu + use ice_grid, only: tmask, lmask_n, lmask_s, tarea, dxu, dyu,HTE,HTN use ice_calendar, only: new_year, write_history, & write_ic, time, histfreq, nstreams, month, & new_month @@ -1668,6 +1793,9 @@ subroutine accum_hist (dt) fhocn, fhocn_ai, uatm, vatm, & fswthru_ai, strairx, strairy, strtltx, strtlty, strintx, strinty, & strocnx, strocny, fm, daidtt, dvidtt, daidtd, dvidtd, fsurf, & +!jd + dvsdtd, & +!jd fcondtop, fsurfn, fcondtopn, flatn, fsensn, albcnt, prs_sig, & fcondbot, fcondbotn, update_ocn_f, & stressp_1, stressm_1, stress12_1, & @@ -1676,6 +1804,7 @@ subroutine accum_hist (dt) stressp_4, stressm_4, stress12_4, sig1, sig2, & mlt_onset, frz_onset, dagedtt, dagedtd, fswint_ai, keffn_top, & snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai, Tbot, Tsnice + use ice_snowphys, only: snowfonice, snow2ocn use ice_atmo, only: formdrag, Cd_atm use ice_meltpond_cesm, only: hs0 use ice_history_shared ! almost everything @@ -1973,6 +2102,12 @@ subroutine accum_hist (dt) call accum_hist_field(n_snoice, iblk, snoice(:,:,iblk), a2D) if (f_dsnow (1:1) /= 'x') & call accum_hist_field(n_dsnow, iblk, dsnow(:,:,iblk), a2D) +!jd + if (f_sn2oc (1:1) /= 'x') & + call accum_hist_field(n_sn2oc, iblk, snow2ocn(:,:,iblk), a2D) + if (f_snfonice (1:1) /= 'x') & + call accum_hist_field(n_snfonice, iblk, snowfonice(:,:,iblk), a2D) +!jd if (f_meltt (1:1) /= 'x') & call accum_hist_field(n_meltt, iblk, meltt(:,:,iblk), a2D) if (f_melts (1:1) /= 'x') & @@ -2193,6 +2328,124 @@ subroutine accum_hist (dt) call accum_hist_field(n_sidmasstrany, iblk, worka(:,:), a2D) endif +!jd noresm-start + + if (f_simasstranx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhoi*p5*(vice(i+1,j,iblk)+vice(i,j,iblk))*HTE(i,j,iblk) & + * p5*(uvel(i,j-1,iblk)+uvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_simasstranx, iblk, worka(:,:), a2D) + endif + + if (f_simasstrany(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhoi*p5*(vice(i,j+1,iblk)+vice(i,j,iblk))*HTN(i,j,iblk) & + * p5*(vvel(i-1,j,iblk)+vvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_simasstrany, iblk, worka(:,:), a2D) + endif + + + if (f_sndmasstranx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhos*p5*(vsno(i+1,j,iblk)+vsno(i,j,iblk))*HTE(i,j,iblk) & + * p5*(uvel(i,j-1,iblk)+uvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sndmasstranx, iblk, worka(:,:), a2D) + endif + + if (f_sndmasstrany(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhos*p5*(vsno(i,j+1,iblk)+vsno(i,j,iblk))*HTN(i,j,iblk) & + * p5*(vvel(i-1,j,iblk)+vvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sndmasstrany, iblk, worka(:,:), a2D) + endif + + + if (f_aicetranx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = p5*(aice(i+1,j,iblk)+aice(i,j,iblk))*HTE(i,j,iblk) & + * p5*(uvel(i,j-1,iblk)+uvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_aicetranx, iblk, worka(:,:), a2D) + endif + + if (f_aicetrany(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = p5*(aice(i,j+1,iblk)+aice(i,j,iblk))*HTN(i,j,iblk) & + * p5*(vvel(i-1,j,iblk)+vvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_aicetrany, iblk, worka(:,:), a2D) + endif + + + if (f_simass(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhoi*vice(i,j,iblk) + enddo + enddo + call accum_hist_field(n_simass, iblk, worka(:,:), a2D) + endif + + if (f_sisnmass(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = rhos*vsno(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sisnmass, iblk, worka(:,:), a2D) + endif + + if (f_sisal(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + do k = 1, nzilyr + worka(i,j) = worka(i,j) + trcr(i,j,nt_sice+k-1,iblk) + enddo + worka(i,j) = aice(i,j,iblk) * worka(i,j) / nzilyr + endif + enddo + enddo + call accum_hist_field(n_sisal, iblk, worka(:,:), a2D) + endif + + + +!jd end + if (f_sistrxdtop(1:1) /= 'x') then worka(:,:) = c0 do j = jlo, jhi @@ -2374,7 +2627,9 @@ subroutine accum_hist (dt) do j = jlo, jhi do i = ilo, ihi if (aice(i,j,iblk) > puny) then - worka(i,j) = evapi(i,j,iblk)*rhoi +!jd worka(i,j) = evapi(i,j,iblk)*rhoi +!jd Evaps has unit kg/m^2, convert to kg/m^2/s + worka(i,j) = evapi(i,j,iblk)/dt endif enddo enddo @@ -2422,7 +2677,9 @@ subroutine accum_hist (dt) do j = jlo, jhi do i = ilo, ihi if (aice(i,j,iblk) > puny) then - worka(i,j) = evaps(i,j,iblk)*rhos +!jd worka(i,j) = evaps(i,j,iblk)*rhos +!jd Evaps has unit kg/m^2, convert to kg/m^2/s + worka(i,j) = evaps(i,j,iblk)/dt endif enddo enddo @@ -2441,12 +2698,52 @@ subroutine accum_hist (dt) call accum_hist_field(n_sndmasssnf, iblk, worka(:,:), a2D) endif +!jd noresm start + if (f_sndmasssi(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = -snoice(i,j,iblk)*rhos/dt + endif + enddo + enddo + call accum_hist_field(n_sndmasssi, iblk, worka(:,:), a2D) + endif + + if (f_sndmasswindrif(1:1) /= 'x') then + worka(:,:)=-snow2ocn(:,:,iblk) +!jd worka(:,:) = c0 +!jd do j = jlo, jhi +!jd do i = ilo, ihi +!jd if (aice(i,j,iblk) > puny) then +!jd worka(i,j) = snow2ocn(i,j,iblk) +!jd endif +!jd enddo +!jd enddo + call accum_hist_field(n_sndmasswindrif, iblk, worka(:,:), a2D) + endif + + if (f_sndmassdyn(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = dvsdtd(i,j,iblk) * rhos + endif + enddo + enddo + call accum_hist_field(n_sndmassdyn, iblk, worka(:,:), a2D) + endif + + +!jd end if (f_sndmassmelt(1:1) /= 'x') then worka(:,:) = c0 do j = jlo, jhi do i = ilo, ihi if (aice(i,j,iblk) > puny) then - worka(i,j) = melts(i,j,iblk)*rhos/dt + worka(i,j) = -melts(i,j,iblk)*rhos/dt endif enddo enddo @@ -3560,7 +3857,22 @@ subroutine accum_hist (dt) enddo ! j endif endif +!jd noresm-start + if (index(avail_hist_fields(n)%vname,'sisal') /= 0) then + if (f_sisal(1:1) /= 'x' .and. n_sisal(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sisal(ns),iblk) = & + a2D(i,j,n_sisal(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sisal(ns),iblk) = spval_dbl + enddo ! i + enddo ! j + endif + endif +!jd end ! back out albedo/zenith angle dependence if (avail_hist_fields(n)%vname(1:6) == 'albice') then do j = jlo, jhi @@ -3594,7 +3906,8 @@ subroutine accum_hist (dt) do k = 1, ncat_hist do j = jlo, jhi do i = ilo, ihi - if (.not. tmask(i,j,iblk) .or. ravgipn(i,j,k) == c0) then ! mask out land points +!jd if (.not. tmask(i,j,iblk) .or. ravgipn(i,j,k) == c0) then ! mask out land points + if (.not. tmask(i,j,iblk) ) then ! mask out land points a3Dc(i,j,k,n,iblk) = spval_dbl else ! convert units a3Dc(i,j,k,n,iblk) = avail_hist_fields(nn)%cona*a3Dc(i,j,k,n,iblk) & @@ -3609,9 +3922,11 @@ subroutine accum_hist (dt) do j = jlo, jhi do i = ilo, ihi if (tmask(i,j,iblk)) then - a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk) = & - a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk)*avgct(ns)*ravgipn(i,j,k) + a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk) = & + a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk)*avgct(ns)*ravgipn(i,j,k) endif + if (ravgipn(i,j,k) == c0) a3Dc(i,j,k,n_siitdthick(ns)-n2D,iblk) = spval_dbl + enddo ! i enddo ! j enddo ! k @@ -3626,6 +3941,7 @@ subroutine accum_hist (dt) a3Dc(i,j,k,n_siitdsnthick(ns)-n2D,iblk) = & a3Dc(i,j,k,n_siitdsnthick(ns)-n2D,iblk)*avgct(ns)*ravgipn(i,j,k) endif + if (ravgipn(i,j,k) == c0) a3Dc(i,j,k,n_siitdsnthick(ns)-n2D,iblk) = spval_dbl enddo ! i enddo ! j enddo ! k diff --git a/src/source/ice_history_shared.F90 b/src/source/ice_history_shared.F90 index 5bc3b420..033d6d13 100644 --- a/src/source/ice_history_shared.F90 +++ b/src/source/ice_history_shared.F90 @@ -215,6 +215,9 @@ module ice_history_shared f_Tref = 'm', f_Qref = 'm', & f_congel = 'm', f_frazil = 'm', & f_snoice = 'm', f_dsnow = 'm', & +!jd + f_sn2oc = 'm', f_snfonice = 'm', & +!jd f_meltt = 'm', f_melts = 'm', & f_meltb = 'm', f_meltl = 'm', & f_fresh = 'm', f_fresh_ai = 'm', & @@ -242,6 +245,13 @@ module ice_history_shared f_sitempbot = 'x', f_sispeed = 'x', & f_siu = 'x', f_siv = 'x', & f_sidmasstranx = 'x', f_sidmasstrany = 'x', & +!jd noresm-start + f_simasstranx = 'x', f_simasstrany = 'x', & + f_sndmasstranx = 'x', f_sndmasstrany = 'x', & + f_aicetranx = 'x', f_aicetrany = 'x', & + f_simass = 'x', f_sisnmass = 'x', & + f_sisal = 'x', & +!jd end f_sistrxdtop = 'x', f_sistrydtop = 'x', & f_sistrxubot = 'x', f_sistryubot = 'x', & f_sicompstren = 'x', & @@ -258,6 +268,11 @@ module ice_history_shared f_sidmassmeltbot = 'x', & f_sidmasslat = 'x', & f_sndmasssnf = 'x', & +!jd noresm start + f_sndmasssi = 'x', & + f_sndmasswindrif = 'x', & + f_sndmassdyn = 'x', & +!jd end f_sndmassmelt = 'x', & f_siflswdtop = 'x', & f_siflswutop = 'x', & @@ -355,6 +370,9 @@ module ice_history_shared f_Tref, f_Qref , & f_congel, f_frazil , & f_snoice, f_dsnow , & +!jd + f_sn2oc, f_snfonice , & +!jd f_meltt, f_melts , & f_meltb, f_meltl , & f_fresh, f_fresh_ai , & @@ -382,6 +400,13 @@ module ice_history_shared f_sitempbot, f_sispeed, & f_siu, f_siv, & f_sidmasstranx, f_sidmasstrany, & +!jd noresm-start + f_simasstranx, f_simasstrany, & + f_sndmasstranx, f_sndmasstrany, & + f_aicetranx, f_aicetrany, & + f_simass , f_sisnmass , & + f_sisal , & +!jd end f_sistrxdtop, f_sistrydtop, & f_sistrxubot, f_sistryubot, & f_sicompstren, & @@ -398,6 +423,11 @@ module ice_history_shared f_sidmassmeltbot, & f_sidmasslat, & f_sndmasssnf, & +!jd noresm start + f_sndmasssi, & + f_sndmasswindrif, & + f_sndmassdyn, & +!jd end f_sndmassmelt, & f_siflswdtop, & f_siflswutop, & @@ -511,6 +541,9 @@ module ice_history_shared n_Tref , n_Qref , & n_congel , n_frazil , & n_snoice , n_dsnow , & +!jd + n_sn2oc , n_snfonice , & +!jd n_meltt , n_melts , & n_meltb , n_meltl , & n_fresh , n_fresh_ai , & @@ -521,6 +554,13 @@ module ice_history_shared n_sitempbot , n_sispeed, & n_siu, n_siv, & n_sidmasstranx, n_sidmasstrany, & +!jd noresm-start + n_simasstranx, n_simasstrany, & + n_sndmasstranx,n_sndmasstrany, & + n_aicetranx, n_aicetrany, & + n_simass , n_sisnmass , & + n_sisal, & +!jd end n_sistrxdtop, n_sistrydtop, & n_sistrxubot, n_sistryubot, & n_sicompstren, & @@ -537,6 +577,11 @@ module ice_history_shared n_sidmassmeltbot, & n_sidmasslat, & n_sndmasssnf, & +!jd noresm start + n_sndmasssi, & + n_sndmasswindrif, & + n_sndmassdyn, & +!jd end n_sndmassmelt, & n_siflswdtop, & n_siflswutop, & diff --git a/src/source/ice_init.F90 b/src/source/ice_init.F90 index 9e954047..6e08f688 100644 --- a/src/source/ice_init.F90 +++ b/src/source/ice_init.F90 @@ -39,7 +39,7 @@ subroutine input_data use ice_age, only: restart_age use ice_broadcast, only: broadcast_scalar, broadcast_array - use ice_constants, only: c0, c1, puny + use ice_constants, only: c0, c1, puny, rhos, ksno use ice_diagnostics, only: diag_file, print_global, print_points, latpnt, lonpnt use ice_domain_size, only: max_nstrm, nilyr, nslyr, max_ntrcr, ncat, n_aero, n_iso use ice_fileunits, only: nu_nml, nu_diag, nml_filename, diag_type, & @@ -83,6 +83,7 @@ subroutine input_data use ice_meltpond_topo, only: hp1, restart_pond_topo use ice_meltpond_lvl, only: restart_pond_lvl, dpscale, frzpnd, & rfracmin, rfracmax, pndaspect, hs1 + use ice_snowphys, only: blowingsnow use ice_aerosol, only: restart_aero use ice_isotope, only: restart_iso use ice_therm_shared, only: ktherm, calc_Tsfc, conduct @@ -152,6 +153,9 @@ subroutine input_data rfracmin, rfracmax, pndaspect, hs1, & hp1 + namelist /snowphys_nml/ & + blowingsnow, rhos, ksno + namelist /forcing_nml/ & atmbndy, fyear_init, ycycle, atm_data_format,& atm_data_type, atm_data_dir, calc_strair, calc_Tsfc, & @@ -266,6 +270,7 @@ subroutine input_data albsnowv = 0.98_dbl_kind ! cold snow albedo, visible albsnowi = 0.70_dbl_kind ! cold snow albedo, near IR ahmax = 0.3_dbl_kind ! thickness above which ice albedo is constant (m) + blowingsnow='none' ! No blowing snow as default, alternative 'lecomte2013' atmbndy = 'default' ! or 'constant' fyear_init = 1900 ! first year of forcing cycle @@ -367,6 +372,9 @@ subroutine input_data print*,'Reading ponds_nml' read(nu_nml, nml=ponds_nml,iostat=nml_error) if (nml_error /= 0) exit + print*,'Reading snowphys_nml' + read(nu_nml, nml=snowphys_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 @@ -746,6 +754,9 @@ subroutine input_data call broadcast_scalar(rfracmin, master_task) call broadcast_scalar(rfracmax, master_task) call broadcast_scalar(pndaspect, master_task) + call broadcast_scalar(blowingsnow, master_task) + call broadcast_scalar(rhos, master_task) + call broadcast_scalar(ksno, master_task) call broadcast_scalar(albicev, master_task) call broadcast_scalar(albicei, master_task) call broadcast_scalar(albsnowv, master_task) @@ -948,6 +959,10 @@ subroutine input_data if (tr_pond .and. .not. tr_pond_lvl) & write(nu_diag,1000) ' pndaspect = ', pndaspect + write(nu_diag,1030) ' blowingsnow = ', & + trim(blowingsnow) + write(nu_diag,1000) ' rhos = ', rhos + write(nu_diag,1000) ' ksno = ', ksno write(nu_diag,1020) ' ktherm = ', ktherm if (ktherm == 1) & write(nu_diag,1030) ' conduct = ', conduct diff --git a/src/source/ice_itd.F90 b/src/source/ice_itd.F90 index e75060d9..756773cf 100644 --- a/src/source/ice_itd.F90 +++ b/src/source/ice_itd.F90 @@ -2563,9 +2563,11 @@ subroutine zerolayer_check (nx_block, ny_block, & n , & ! category index ij ! combined horizontal index - real (kind=dbl_kind), parameter :: & - max_error = puny*Lfresh*rhos ! max error in zero layer energy check - ! (so max volume error = puny) +!jd real (kind=dbl_kind), parameter :: & +!jd max_error = puny*Lfresh*rhos ! max error in zero layer energy check +!jd ! (so max volume error = puny) + real (kind=dbl_kind) :: & + max_error ! max error in zero layer energy check real (kind=dbl_kind), dimension (nx_block,ny_block,ncat) :: & eicen ! energy of melting for each ice layer (J/m^2) @@ -2583,7 +2585,10 @@ subroutine zerolayer_check (nx_block, ny_block, & !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- - +!jd + max_error = puny*Lfresh*rhos ! max error in zero layer energy check + ! (so max volume error = puny) +!jd l_stop = .false. istop = 0 jstop = 0 diff --git a/src/source/ice_orbital.F90 b/src/source/ice_orbital.F90 index 0f36287c..bf8b472a 100644 --- a/src/source/ice_orbital.F90 +++ b/src/source/ice_orbital.F90 @@ -68,6 +68,11 @@ subroutine compute_coszen (nx_block, ny_block, & use ice_calendar, only: yday, sec, calendar_type, nextsw_cday, days_per_year use ice_constants, only: c0, c2, p5, pi, secday use shr_orb_mod, only: shr_orb_decl +!+tht +#ifdef CESMCOUPLED + use shr_orb_mod, only: shr_orb_cosz +#endif +!-tht integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -116,6 +121,21 @@ subroutine compute_coszen (nx_block, ny_block, & coszen(:,:) = c0 ! sun at horizon +!+tht +#ifdef CESMCOUPLED +!DIR$ CONCURRENT !Cray +!cdir nodep !NEC +!ocl novrec !Fujitsu + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) +!+tht outlined again + coszen(i,j) = shr_orb_cosz(ydayp1, & + tlat(i,j),tlon(i,j),delta,dt) + enddo + + endif +#else !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu @@ -125,15 +145,12 @@ subroutine compute_coszen (nx_block, ny_block, & !lipscomb - function inlined to improve vector efficiency ! coszen(i,j) = shr_orb_cosz(ydayp1, & ! tlat(i,j),tlon(i,j),delta) - coszen(i,j) = sin(tlat(i,j))*sin(delta) & + cos(tlat(i,j))*cos(delta) & *cos((sec/secday-p5)*c2*pi + tlon(i,j)) !cos(hour angle) enddo - -#ifdef CESMCOUPLED - endif #endif +!-tht end subroutine compute_coszen diff --git a/src/source/ice_snowphys_mod.F90 b/src/source/ice_snowphys_mod.F90 new file mode 100644 index 00000000..e7a14055 --- /dev/null +++ b/src/source/ice_snowphys_mod.F90 @@ -0,0 +1,62 @@ + + +! Modifications of amount and physical properties of snow on ice. + +! author: Jens Boldingh Debernard, MET-Norway + + + module ice_snowphys + + use ice_kinds_mod + use ice_blocks, only: nx_block, ny_block + use ice_domain_size, only: ncat + use ice_constants + use ice_domain_size, only: max_blocks + + implicit none + save + + private + public :: snowphys_snowfonice + + character (len=char_len), public :: & + blowingsnow ! pond refreezing parameterization + + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat,max_blocks), public :: & + snowfonicen, & ! Fraction of snow kept on ice during snowfall + snow2ocnn ! Mass fluks of snow blowing indirectly into ocean + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: & + snowfonice, & ! Fraction of snow kept on ice during snowfall (agregate) + snow2ocn ! Mass fluks of snow blowing indirectly into ocean (agreate) + + contains + + subroutine snowphys_snowfonice(icells, indxi, indxj, snwfonice,ai) + integer (kind=int_kind), intent(in) :: & + icells ! number of cells with aicen > puny + + integer (kind=int_kind), dimension(nx_block*ny_block), & + intent(in) :: & + indxi, indxj ! compressed indices for cells with aicen > puny + + + real (kind=dbl_kind), dimension (nx_block,ny_block),intent(in) :: ai + real (kind=dbl_kind), dimension (nx_block,ny_block),intent(inout) :: & + snwfonice + + integer (kind=int_kind) :: ij, i,j + + do ij=1,icells + i=indxi(ij) + j=indxj(ij) + +!jd snwfonice(i,j)=0.9 +! if (aice(i,j,iblk) > p01 ) & + snwfonice(i,j)=(c1 - max(c0,c1-ai(i,j))**p6 )/ai(i,j) + enddo + + end subroutine snowphys_snowfonice + + end module ice_snowphys diff --git a/src/source/ice_step_mod.F90 b/src/source/ice_step_mod.F90 index 2e311921..c0230af6 100644 --- a/src/source/ice_step_mod.F90 +++ b/src/source/ice_step_mod.F90 @@ -201,6 +201,10 @@ subroutine step_therm1 (dt, iblk) use ice_therm_shared, only: calc_Tsfc use ice_therm_vertical, only: frzmlt_bottom_lateral, thermo_vertical use ice_timers, only: ice_timer_start, ice_timer_stop, timer_ponds +!jd - Blowing snow modification + use ice_snowphys, only: blowingsnow,snowphys_snowfonice, & + snowfonicen, snow2ocnn, snowfonice, snow2ocn +!jd real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -269,7 +273,7 @@ subroutine step_therm1 (dt, iblk) istop, jstop ! indices of grid cell where model aborts real (kind=dbl_kind), dimension (nx_block,ny_block) :: & - worka, workb + worka, workb, fsnown l_stop = .false. @@ -370,6 +374,7 @@ subroutine step_therm1 (dt, iblk) dfloe (:,:,iblk), ncat) endif + do n = 1, ncat meltsn(:,:,n,iblk) = c0 @@ -379,6 +384,15 @@ subroutine step_therm1 (dt, iblk) snoicen(:,:,n,iblk) = c0 dsnown(:,:,n,iblk) = c0 ! Tsf_icen(:,:,n,iblk) = c0 +!jd + if (trim(blowingsnow)=='lecomte2013') then + snowfonicen(:,:,n,iblk) = c0 + else + snowfonicen(:,:,n,iblk) = c1 + end if + fsnown(:,:) = fsnow(:,:,iblk) + snow2ocnn(:,:,n,iblk) = c0 +!jd !----------------------------------------------------------------- ! Identify cells with nonzero ice area @@ -492,6 +506,23 @@ subroutine step_therm1 (dt, iblk) trcrn(:,:,nt_FY,n,iblk)) endif + !----------------------------------------------------------------- + ! Calculate amount of snow blowing directly into the ocean from this + ! cathegory. (jd) First part, split + !----------------------------------------------------------------- + if (trim(blowingsnow)=='lecomte2013') then + call snowphys_snowfonice(icells, & + indxi, indxj, & + snowfonicen(:,:,n,iblk), & + aice(:,:,iblk)) + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + fsnown(i,j)=fsnow(i,j,iblk)*snowfonicen(i,j,n,iblk) + snow2ocnn(i,j,n,iblk)=fsnow(i,j,iblk) - fsnown(i,j) + end do + + end if !----------------------------------------------------------------- ! Vertical thermodynamics: Heat conduction, growth and melting. !----------------------------------------------------------------- @@ -524,7 +555,8 @@ subroutine step_therm1 (dt, iblk) vicen(:,:,n,iblk), vsnon(:,:,n,iblk), & flw (:,:,iblk), potT (:,:,iblk), & Qa (:,:,iblk), rhoa (:,:,iblk), & - fsnow (:,:,iblk), fpond (:,:,iblk), & +!jd fsnow (:,:,iblk), fpond (:,:,iblk), & + fsnown (:,:), fpond (:,:,iblk), & fbot, Tbotn, & Tsnicen, & sss (:,:,iblk), & @@ -580,6 +612,18 @@ subroutine step_therm1 (dt, iblk) enddo enddo + !----------------------------------------------------------------- + ! Snowfall and energy to the ocean + !----------------------------------------------------------------- + + do ij = 1,icells + i = indxi(ij) + j = indxj(ij) + fhocnn(i,j) = fhocnn(i,j) & + - Lfresh*snow2ocnn(i,j,n,iblk) + freshn(i,j) = freshn(i,j) + snow2ocnn(i,j,n,iblk) + enddo + !----------------------------------------------------------------- ! Aerosol update !----------------------------------------------------------------- @@ -593,7 +637,8 @@ subroutine step_therm1 (dt, iblk) meltbn(:,:,n,iblk), & congeln(:,:,n,iblk), & snoicen(:,:,n,iblk), & - fsnow(:,:,iblk), & +!jd fsnow(:,:,iblk), & + fsnown(:,:), & trcrn(:,:,:,n,iblk), & aicen_init(:,:,n,iblk), & vicen_init(:,:,n,iblk), & @@ -619,7 +664,8 @@ subroutine step_therm1 (dt, iblk) congeln(:,:,n,iblk), & snoicen(:,:,n,iblk), & evapn, & - fsnow(:,:,iblk), & +!jd fsnow(:,:,iblk), & + fsnown(:,:), & Qrefn_iso(:,:,:), & trcrn(:,:,:,n,iblk), & aicen_init(:,:,n,iblk), & @@ -755,6 +801,10 @@ subroutine step_therm1 (dt, iblk) meltt (:,:,iblk), melts (:,:,iblk), & meltb (:,:,iblk), & congel (:,:,iblk), snoice (:,:,iblk), & +!jd + snow2ocnn(:,:,n,iblk), snow2ocn(:,:,iblk), & + snowfonicen(:,:,n,iblk), snowfonice(:,:,iblk), & +!jd Uref=Uref(:,:,iblk), Urefn=Urefn, & Qref_iso=Qref_iso(:,:,:,iblk), & Qrefn_iso=Qrefn_iso(:,:,:),& @@ -1164,7 +1214,8 @@ subroutine step_dynamics (dt, ndtd) use ice_dyn_evp, only: evp use ice_dyn_eap, only: eap use ice_dyn_shared, only: kdyn - use ice_flux, only: daidtd, dvidtd, init_history_dyn, dagedtd +!jd use ice_flux, only: daidtd, dvidtd, init_history_dyn, dagedtd + use ice_flux, only: daidtd, dvidtd, dvsdtd, init_history_dyn, dagedtd use ice_grid, only: tmask use ice_itd, only: aggregate use ice_state, only: nt_qsno, trcrn, vsnon, aicen, vicen, ntrcr, & @@ -1262,6 +1313,9 @@ subroutine step_dynamics (dt, ndtd) do j = jlo,jhi do i = ilo,ihi dvidtd(i,j,iblk) = (vice(i,j,iblk) - dvidtd(i,j,iblk)) /dt +!jd + dvsdtd(i,j,iblk) = (vsno(i,j,iblk) - dvsdtd(i,j,iblk)) /dt +!jd daidtd(i,j,iblk) = (aice(i,j,iblk) - daidtd(i,j,iblk)) /dt if (tr_iage) & dagedtd(i,j,iblk)= (trcr(i,j,nt_iage,iblk)-dagedtd(i,j,iblk))/dt diff --git a/src/source/ice_therm_mushy.F90 b/src/source/ice_therm_mushy.F90 index 13fd5e66..a35d7e07 100644 --- a/src/source/ice_therm_mushy.F90 +++ b/src/source/ice_therm_mushy.F90 @@ -3708,9 +3708,12 @@ function temperature_snow(zqsn) result(zTsn) zTsn ! snow layer temperature (C) real(kind=dbl_kind), parameter :: & - A = c1 / (rhos * cp_ice) , & +!jd A = c1 / (rhos * cp_ice) , & B = Lfresh / cp_ice - +!jd + real(kind=dbl_kind) :: A + A = c1 / (rhos * cp_ice) +!jd zTsn = A * zqsn + B end function temperature_snow diff --git a/src/source/ice_therm_vertical.F90 b/src/source/ice_therm_vertical.F90 index 65fd0c95..c2e3028d 100644 --- a/src/source/ice_therm_vertical.F90 +++ b/src/source/ice_therm_vertical.F90 @@ -241,6 +241,9 @@ subroutine thermo_vertical (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block) :: & fadvocn ! advective heat flux to ocean +!jd Approximation of Tsnice + real (kind=dbl_kind) :: ksh, kih + !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- @@ -422,8 +425,14 @@ subroutine thermo_vertical (nx_block, ny_block, & i = indxi(ij) j = indxj(ij) - Tsnice(i,j) = (hslyr(ij)*zTsn(ij,nslyr) + hilyr(ij)*zTin(ij,1)) & - / (hslyr(ij)+hilyr(ij)) +!jd Tsnice(i,j) = (hslyr(ij)*zTsn(ij,nslyr) + hilyr(ij)*zTin(ij,1)) & +!jd / (hslyr(ij)+hilyr(ij)) + + ksh = ksno*hilyr(ij) + kih = kice*hslyr(ij) + Tsnice(i,j) = (ksh*zTsn(ij,nslyr) + kih*zTin(ij,1)) / (ksh + kih ) + + enddo !----------------------------------------------------------------- @@ -2387,7 +2396,7 @@ subroutine conservation_check_vthermo(nx_block, ny_block, & - fsnow(i,j)*Lfresh - fadvocn(i,j)) * dt ferr = abs(efinal(ij)-einit(ij)-einp) / dt - if (ferr > ferrmax) then + if (ferr > c2*ferrmax) then l_stop = .true. istop = i jstop = j