From 3c99e1067091354e2cce56dd35f14749612f87dc Mon Sep 17 00:00:00 2001 From: dabail10 Date: Thu, 27 Sep 2018 09:38:38 -0600 Subject: [PATCH] Update CICE with CMIP changes. (#191) * Remove some dummy arguments and unused variables. * Dummy variable cleanup * Cleanup * Remove icepack from commit * Reset icepack master * Fix * Fix * Update icepack * fix * Fix flush ifdef logic * Uninitialized variables. * Additional initialization * Update dummy * more dummy fixes * More hobart changes * Update Hobart Intel settings * Additional Hobart Intel changes * Update icepack * Some tweaks to dummy variables * Remove some initialization code. * Initialization needed for CICE * Add gravit to ice_dyn_shared.F90 * CMIP history changes * CMIP changes * Update icepack * Redo istat changes * Update icepack to CMIP --- cicecore/cicedynB/analysis/ice_history.F90 | 2057 +++++++++++++++-- .../cicedynB/analysis/ice_history_shared.F90 | 155 ++ cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 11 + cicecore/cicedynB/general/ice_flux.F90 | 31 +- cicecore/cicedynB/general/ice_step_mod.F90 | 27 +- cicecore/cicedynB/infrastructure/ice_grid.F90 | 4 +- .../io/io_binary/ice_history_write.F90 | 1 + .../io/io_netcdf/ice_history_write.F90 | 3 + .../io/io_pio/ice_history_write.F90 | 3 + cicecore/drivers/cice/CICE_RunMod.F90 | 16 +- icepack | 2 +- 11 files changed, 2123 insertions(+), 187 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index f06559608..8ef859b65 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/cicecore/cicedynB/analysis/ice_history.F90 @@ -1,5 +1,4 @@ !======================================================================= - ! Driver for core history output ! ! The following variables are currently hard-wired as snapshots @@ -31,7 +30,7 @@ module ice_history use ice_kinds_mod - use ice_constants, only: c0, c1, c2, c100, p25, & + use ice_constants, only: c0, c1, c2, c100, p001, p25, p5, & mps_to_cmpdy, kg_to_g, spval use ice_fileunits, only: nu_nml, nml_filename, nu_diag, & get_fileunit, release_fileunit @@ -184,6 +183,84 @@ subroutine init_hist (dt) ! to prevent array-out-of-bounds when aggregating if (f_fmeltt_ai(1:1) /= 'x') f_fmelttn_ai = f_fmeltt_ai + ! Turn on all CMIP fields in one go. + + if (f_CMIP(1:1) /= 'x') then + f_sithick = 'mxxxx' + f_sisnthick = 'mxxxx' + f_siage = 'mxxxx' + f_sitemptop = 'mxxxx' + f_sitempsnic = 'mxxxx' + f_sitempbot = 'mxxxx' + f_sispeed = 'mxxxx' + f_siu = 'mxxxx' + f_siv = 'mxxxx' + f_sidmasstranx = 'mxxxx' + f_sidmasstrany = 'mxxxx' + f_sistrxdtop = 'mxxxx' + f_sistrydtop = 'mxxxx' + f_sistrxubot = 'mxxxx' + f_sistryubot = 'mxxxx' + f_sicompstren = 'mxxxx' + f_sialb = 'mxxxx' + f_sihc = 'mxxxx' + f_sisnhc = 'mxxxx' + f_sidconcth = 'mxxxx' + f_sidconcdyn = 'mxxxx' + f_sidmassth = 'mxxxx' + f_sidmassdyn = 'mxxxx' + f_sidmassgrowthwat = 'mxxxx' + f_sidmassgrowthbot = 'mxxxx' + f_sidmasssi = 'mxxxx' + f_sidmassevapsubl = 'mxxxx' + f_sndmasssubl = 'mxxxx' + f_sidmassmelttop = 'mxxxx' + f_sidmassmeltbot = 'mxxxx' + f_sidmasslat = 'mxxxx' + f_sndmasssnf = 'mxxxx' + f_sndmassmelt = 'mxxxx' + f_siflswdtop = 'mxxxx' + f_siflswutop = 'mxxxx' + f_siflswdbot = 'mxxxx' + f_sifllwdtop = 'mxxxx' + f_sifllwutop = 'mxxxx' + f_siflsenstop = 'mxxxx' + f_siflsensupbot = 'mxxxx' + f_sifllatstop = 'mxxxx' + f_siflcondtop = 'mxxxx' + f_siflcondbot = 'mxxxx' + f_sipr = 'mxxxx' + f_sifb = 'mxxxx' + f_siflsaltbot = 'mxxxx' + f_siflfwbot = 'mxxxx' + f_siflfwdrain = 'mxxxx' + f_siforcetiltx = 'mxxxx' + f_siforcetilty = 'mxxxx' + f_siforcecoriolx = 'mxxxx' + f_siforcecorioly = 'mxxxx' + f_siforceintstrx = 'mxxxx' + f_siforceintstry = 'mxxxx' + f_sidragtop = 'mxxxx' + f_sistreave = 'mxxxx' + f_sistremax = 'mxxxx' + f_sirdgthick = 'mxxxx' + f_siitdconc = 'mxxxx' + f_siitdthick = 'mxxxx' + f_siitdsnthick = 'mxxxx' + f_aicen = 'mxxxx' + endif + + if (f_CMIP(2:2) == 'd') then + f_icepresent = f_CMIP + f_aice = f_CMIP + f_sithick = f_CMIP + f_sisnthick = f_CMIP + f_sitemptop = f_CMIP + f_siu = f_CMIP + f_siv = f_CMIP + f_sispeed = f_CMIP + endif + #ifndef ncdf f_bounds = .false. #endif @@ -311,6 +388,68 @@ subroutine init_hist (dt) call broadcast_scalar (f_frz_onset, master_task) call broadcast_scalar (f_aisnap, master_task) call broadcast_scalar (f_hisnap, master_task) + call broadcast_scalar (f_sithick, master_task) + call broadcast_scalar (f_siage, master_task) + call broadcast_scalar (f_sisnthick, master_task) + call broadcast_scalar (f_sitemptop, master_task) + call broadcast_scalar (f_sitempsnic, master_task) + call broadcast_scalar (f_sitempbot, master_task) + call broadcast_scalar (f_siu, master_task) + call broadcast_scalar (f_siv, master_task) + call broadcast_scalar (f_sidmasstranx, master_task) + call broadcast_scalar (f_sidmasstrany, master_task) + call broadcast_scalar (f_sistrxdtop, master_task) + call broadcast_scalar (f_sistrydtop, master_task) + call broadcast_scalar (f_sistrxubot, master_task) + call broadcast_scalar (f_sistryubot, master_task) + call broadcast_scalar (f_sicompstren, master_task) + call broadcast_scalar (f_sispeed, master_task) + call broadcast_scalar (f_sialb, master_task) + call broadcast_scalar (f_sihc, master_task) + call broadcast_scalar (f_sisnhc, master_task) + call broadcast_scalar (f_sidconcth, master_task) + call broadcast_scalar (f_sidconcdyn, master_task) + call broadcast_scalar (f_sidmassth, master_task) + call broadcast_scalar (f_sidmassdyn, master_task) + call broadcast_scalar (f_sidmassgrowthwat, master_task) + call broadcast_scalar (f_sidmassgrowthbot, master_task) + call broadcast_scalar (f_sidmasssi, master_task) + call broadcast_scalar (f_sidmassevapsubl, master_task) + call broadcast_scalar (f_sndmasssubl, master_task) + call broadcast_scalar (f_sidmassmelttop, master_task) + call broadcast_scalar (f_sidmassmeltbot, master_task) + call broadcast_scalar (f_sidmasslat, master_task) + call broadcast_scalar (f_sndmasssnf, master_task) + call broadcast_scalar (f_sndmassmelt, master_task) + call broadcast_scalar (f_siflswdtop, master_task) + call broadcast_scalar (f_siflswutop, master_task) + call broadcast_scalar (f_siflswdbot, master_task) + call broadcast_scalar (f_sifllwdtop, master_task) + call broadcast_scalar (f_sifllwutop, master_task) + call broadcast_scalar (f_siflsenstop, master_task) + call broadcast_scalar (f_siflsensupbot, master_task) + call broadcast_scalar (f_sifllatstop, master_task) + call broadcast_scalar (f_siflcondtop, master_task) + call broadcast_scalar (f_siflcondbot, master_task) + call broadcast_scalar (f_sipr, master_task) + call broadcast_scalar (f_sifb, master_task) + call broadcast_scalar (f_siflsaltbot, master_task) + call broadcast_scalar (f_siflfwbot, master_task) + call broadcast_scalar (f_siflfwdrain, master_task) + call broadcast_scalar (f_siforcetiltx, master_task) + call broadcast_scalar (f_siforcetilty, master_task) + call broadcast_scalar (f_siforcecoriolx, master_task) + call broadcast_scalar (f_siforcecorioly, master_task) + call broadcast_scalar (f_siforceintstrx, master_task) + call broadcast_scalar (f_siforceintstry, master_task) + call broadcast_scalar (f_siitdconc, master_task) + call broadcast_scalar (f_siitdthick, master_task) + call broadcast_scalar (f_siitdsnthick, master_task) + call broadcast_scalar (f_sidragtop, master_task) + call broadcast_scalar (f_sistreave, master_task) + call broadcast_scalar (f_sistremax, master_task) + call broadcast_scalar (f_sirdgthick, master_task) + call broadcast_scalar (f_aicen, master_task) call broadcast_scalar (f_vicen, master_task) call broadcast_scalar (f_vsnon, master_task) @@ -467,7 +606,7 @@ subroutine init_hist (dt) "if >0, new ice forms; if <0, ice melts", c1, c0, & ns1, f_frzmlt) - call define_hist_field(n_fswfac,"fswfac","1",tstr2D, tcstr, & + call define_hist_field(n_fswfac,"scale_factor","1",tstr2D, tcstr, & "shortwave scaling factor", & "ratio of netsw new:old", c1, c0, & ns1, f_fswfac) @@ -928,6 +1067,298 @@ subroutine init_hist (dt) "weighted by ice area", c1, c0, & ns1, f_FY) + ! CMIP 2D variables + + call define_hist_field(n_sithick,"sithick","m",tstr2D, tcstr, & + "sea ice thickness", & + "volume divided by area", c1, c0, & + ns1, f_sithick) + + call define_hist_field(n_siage,"siage","s",tstr2D, tcstr, & + "sea ice age", & + "none", c1, c0, & + ns1, f_siage) + + call define_hist_field(n_sisnthick,"sisnthick","m",tstr2D, tcstr, & + "sea ice snow thickness", & + "snow volume divided by area", c1, c0, & + ns1, f_sisnthick) + + call define_hist_field(n_sitemptop,"sitemptop","K",tstr2D, tcstr, & + "sea ice surface temperature", & + "none", c1, c0, & + ns1, f_sitemptop) + + call define_hist_field(n_sitempsnic,"sitempsnic","K",tstr2D, tcstr, & + "snow ice interface temperature", & + "surface temperature when no snow present", c1, c0, & + ns1, f_sitempsnic) + + call define_hist_field(n_sitempbot,"sitempbot","K",tstr2D, tcstr, & + "sea ice bottom temperature", & + "none", c1, c0, & + ns1, f_sitempbot) + + call define_hist_field(n_siu,"siu","m/s",ustr2D, ucstr, & + "ice x velocity component", & + "none", c1, c0, & + ns1, f_siu) + + call define_hist_field(n_siv,"siv","m/s",ustr2D, ucstr, & + "ice y velocity component", & + "none", c1, c0, & + ns1, f_siv) + + call define_hist_field(n_sidmasstranx,"sidmasstranx","kg/s",ustr2D, ucstr, & + "x component of snow and sea ice mass transport", & + "none", c1, c0, & + ns1, f_sidmasstranx) + + call define_hist_field(n_sidmasstrany,"sidmasstrany","kg/s",ustr2D, ucstr, & + "y component of snow and sea ice mass transport", & + "none", c1, c0, & + ns1, f_sidmasstrany) + + call define_hist_field(n_sistrxdtop,"sistrxdtop","N m-2",ustr2D, ucstr, & + "x component of atmospheric stress on sea ice", & + "none", c1, c0, & + ns1, f_sistrxdtop) + + call define_hist_field(n_sistrydtop,"sistrydtop","N m-2",ustr2D, ucstr, & + "y component of atmospheric stress on sea ice", & + "none", c1, c0, & + ns1, f_sistrydtop) + + call define_hist_field(n_sistrxubot,"sistrxubot","N m-2",ustr2D, ucstr, & + "x component of ocean stress on sea ice", & + "none", c1, c0, & + ns1, f_sistrxubot) + + call define_hist_field(n_sistryubot,"sistryubot","N m-2",ustr2D, ucstr, & + "y component of ocean stress on sea ice", & + "none", c1, c0, & + ns1, f_sistryubot) + + call define_hist_field(n_sicompstren,"sicompstren","N m-1",tstr2D, tcstr, & + "compressive sea ice strength", & + "none", c1, c0, & + ns1, f_sicompstren) + + call define_hist_field(n_sispeed,"sispeed","m/s",ustr2D, ucstr, & + "ice speed", & + "none", c1, c0, & + ns1, f_sispeed) + + call define_hist_field(n_sialb,"sialb","1",tstr2D, tcstr, & + "sea ice albedo", & + "none", c1, c0, & + ns1, f_sialb) + + call define_hist_field(n_sihc,"sihc","J m-2",tstr2D, tcstr, & + "sea ice heat content", & + "none", c1, c0, & + ns1, f_sihc) + + call define_hist_field(n_sisnhc,"sisnhc","J m-2",tstr2D, tcstr, & + "snow heat content", & + "none", c1, c0, & + ns1, f_sisnhc) + + call define_hist_field(n_sidconcth,"sidconcth","1/s",tstr2D, tcstr, & + "sea ice area change from thermodynamics", & + "none", c1, c0, & + ns1, f_sidconcth) + + call define_hist_field(n_sidconcdyn,"sidconcdyn","1/s",tstr2D, tcstr, & + "sea ice area change from dynamics", & + "none", c1, c0, & + ns1, f_sidconcdyn) + + call define_hist_field(n_sidmassth,"sidmassth","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from thermodynamics", & + "none", c1, c0, & + ns1, f_sidmassth) + + call define_hist_field(n_sidmassdyn,"sidmassdyn","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from dynamics", & + "none", c1, c0, & + ns1, f_sidmassdyn) + + call define_hist_field(n_sidmassgrowthwat,"sidmassgrowthwat","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from frazil", & + "none", c1, c0, & + ns1, f_sidmassgrowthwat) + + call define_hist_field(n_sidmassgrowthbot,"sidmassgrowthbot","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from basal growth", & + "none", c1, c0, & + ns1, f_sidmassgrowthbot) + + call define_hist_field(n_sidmasssi,"sidmasssi","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from snow-ice formation", & + "none", c1, c0, & + ns1, f_sidmasssi) + + call define_hist_field(n_sidmassevapsubl,"sidmassevapsubl","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change from evaporation and sublimation", & + "none", c1, c0, & + ns1, f_sidmassevapsubl) + + call define_hist_field(n_sndmasssubl,"sndmassubl","kg m-2 s-1",tstr2D, tcstr, & + "snow mass change from evaporation and sublimation", & + "none", c1, c0, & + ns1, f_sndmasssubl) + + call define_hist_field(n_sidmassmelttop,"sidmassmelttop","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change top melt", & + "none", c1, c0, & + ns1, f_sidmassmelttop) + + call define_hist_field(n_sidmassmeltbot,"sidmassmeltbot","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change bottom melt", & + "none", c1, c0, & + ns1, f_sidmassmeltbot) + + call define_hist_field(n_sidmasslat,"sidmasslat","kg m-2 s-1",tstr2D, tcstr, & + "sea ice mass change lateral melt", & + "none", c1, c0, & + ns1, f_sidmasslat) + + call define_hist_field(n_sndmasssnf,"sndmasssnf","kg m-2 s-1",tstr2D, tcstr, & + "snow mass change from snow fall", & + "none", c1, c0, & + ns1, f_sndmasssnf) + + call define_hist_field(n_sndmassmelt,"sndmassmelt","kg m-2 s-1",tstr2D, tcstr, & + "snow mass change from snow melt", & + "none", c1, c0, & + ns1, f_sndmassmelt) + + call define_hist_field(n_siflswdtop,"siflswdtop","W/m2",tstr2D, tcstr, & + "down shortwave flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_siflswdtop) + + call define_hist_field(n_siflswutop,"siflswutop","W/m2",tstr2D, tcstr, & + "upward shortwave flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_siflswutop) + + call define_hist_field(n_siflswdbot,"siflswdbot","W/m2",tstr2D, tcstr, & + "down shortwave flux at bottom of ice", & + "positive downward", c1, c0, & + ns1, f_siflswdbot) + + call define_hist_field(n_sifllwdtop,"sifllwdtop","W/m2",tstr2D, tcstr, & + "down longwave flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_sifllwdtop) + + call define_hist_field(n_sifllwutop,"sifllwutop","W/m2",tstr2D, tcstr, & + "upward longwave flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_sifllwutop) + + call define_hist_field(n_siflsenstop,"siflsenstop","W/m2",tstr2D, tcstr, & + "sensible heat flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_siflsenstop) + + call define_hist_field(n_siflsensupbot,"siflsensupbot","W/m2",tstr2D, tcstr, & + "sensible heat flux at bottom of sea ice", & + "positive downward", c1, c0, & + ns1, f_siflsensupbot) + + call define_hist_field(n_sifllatstop,"sifllatstop","W/m2",tstr2D, tcstr, & + "latent heat flux over sea ice", & + "positive downward", c1, c0, & + ns1, f_sifllatstop) + + call define_hist_field(n_siflcondtop,"siflcondtop","W/m2",tstr2D, tcstr, & + "conductive heat flux at top of sea ice", & + "positive downward", c1, c0, & + ns1, f_siflcondtop) + + call define_hist_field(n_siflcondbot,"siflcondbot","W/m2",tstr2D, tcstr, & + "conductive heat flux at bottom of sea ice", & + "positive downward", c1, c0, & + ns1, f_siflcondbot) + + call define_hist_field(n_sipr,"sipr","kg m-2 s-1",tstr2D, tcstr, & + "rainfall over sea ice", & + "none", c1, c0, & + ns1, f_sipr) + + call define_hist_field(n_sifb,"sifb","m",tstr2D, tcstr, & + "sea ice freeboard above sea level", & + "none", c1, c0, & + ns1, f_sifb) + + call define_hist_field(n_siflsaltbot,"siflsaltbot","kg m-2 s-1",tstr2D, tcstr, & + "salt flux from sea ice", & + "positive downward", c1, c0, & + ns1, f_siflsaltbot) + + call define_hist_field(n_siflfwbot,"siflfwbot","kg m-2 s-1",tstr2D, tcstr, & + "fresh water flux from sea ice", & + "positive downward", c1, c0, & + ns1, f_siflfwbot) + + call define_hist_field(n_siflfwdrain,"siflfwdrain","kg m-2 s-1",tstr2D, tcstr, & + "fresh water drainage through sea ice", & + "positive downward", c1, c0, & + ns1, f_siflfwdrain) + + call define_hist_field(n_sidragtop,"sidragtop","1",tstr2D, tcstr, & + "atmospheric drag over sea ice", & + "none", c1, c0, & + ns1, f_sidragtop) + + call define_hist_field(n_sirdgthick,"sirdgthick","m",tstr2D, tcstr, & + "sea ice ridge thickness", & + "vrdg divided by ardg", c1, c0, & + ns1, f_sirdgthick) + + call define_hist_field(n_siforcetiltx,"siforcetiltx","N m-2",tstr2D, tcstr, & + "sea surface tilt term", & + "none", c1, c0, & + ns1, f_siforcetiltx) + + call define_hist_field(n_siforcetilty,"siforcetilty","N m-2",tstr2D, tcstr, & + "sea surface tile term", & + "none", c1, c0, & + ns1, f_siforcetilty) + + call define_hist_field(n_siforcecoriolx,"siforcecoriolx","N m-2",tstr2D, tcstr, & + "coriolis term", & + "none", c1, c0, & + ns1, f_siforcecoriolx) + + call define_hist_field(n_siforcecorioly,"siforcecorioly","N m-2",tstr2D, tcstr, & + "coriolis term", & + "none", c1, c0, & + ns1, f_siforcecorioly) + + call define_hist_field(n_siforceintstrx,"siforceintstrx","N m-2",tstr2D, tcstr, & + "internal stress term", & + "none", c1, c0, & + ns1, f_siforceintstrx) + + call define_hist_field(n_siforceintstry,"siforceintstry","N m-2",tstr2D, tcstr, & + "internal stress term", & + "none", c1, c0, & + ns1, f_siforceintstry) + + call define_hist_field(n_sistreave,"sistreave","N m-1",ustr2D, ucstr, & + "average normal stress", & + "sistreave is instantaneous", c1, c0, & + ns1, f_sistreave) + + call define_hist_field(n_sistremax,"sistremax","N m-1",ustr2D, ucstr, & + "maximum shear stress", & + "sistremax is instantaneous", c1, c0, & + ns1, f_sistremax) + endif ! if (histfreq(ns1) /= 'x') then enddo ! ns1 @@ -993,6 +1424,19 @@ subroutine init_hist (dt) "multilayer scheme", c1, c0, & ns1, f_keffn_top) + ! CMIP 3D + call define_hist_field(n_siitdconc,"siitdconc","1",tstr3Dc, tcstr, & + "ice area, categories","none", c1, c0, & + ns1, f_siitdconc) + + call define_hist_field(n_siitdthick,"siitdthick","m",tstr3Dc, tcstr, & + "ice thickness, categories","none", c1, c0, & + ns1, f_siitdthick) + + call define_hist_field(n_siitdsnthick,"siitdsnthick","m",tstr3Dc, tcstr, & + "snow thickness, categories","none", c1, c0, & + ns1, f_siitdsnthick) + endif ! if (histfreq(ns1) /= 'x') then enddo ! ns1 @@ -1190,7 +1634,7 @@ subroutine accum_hist (dt) use ice_blocks, only: block, get_block, nx_block, ny_block use ice_domain, only: blocks_ice, nblocks - use ice_grid, only: tmask, lmask_n, lmask_s + use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu use ice_calendar, only: new_year, write_history, & write_ic, time, histfreq, nstreams, month, & new_month @@ -1198,27 +1642,28 @@ subroutine accum_hist (dt) yieldstress11, yieldstress12, yieldstress22 use ice_dyn_shared, only: kdyn, principal_stress use ice_flux, only: fsw, flw, fsnow, frain, sst, sss, uocn, vocn, & - frzmlt_init, fswfac, fswabs, fswthru, alvdr, alvdf, alidr, alidf, & - albice, albsno, albpnd, coszen, flat, fsens, flwout, evap, & - Tair, Tref, Qref, congel, frazil, snoice, dsnow, & + frzmlt_init, scale_factor, fswabs, fswthru, alvdr, alvdf, alidr, alidf, & + albice, albsno, albpnd, coszen, flat, fsens, flwout, evap, evaps, evapi, & + Tair, Tref, Qref, congel, frazil, frazil_diag, snoice, dsnow, & melts, meltb, meltt, meltl, fresh, fsalt, fresh_ai, fsalt_ai, & - fhocn, fhocn_ai, uatm, vatm, fbot, & + 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, fsurfn, fcondtopn, flatn, fsensn, albcnt, & + fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, & stressp_1, stressm_1, stress12_1, & stressp_2, & stressp_3, & stressp_4, sig1, sig2, sigP, & mlt_onset, frz_onset, dagedtt, dagedtd, fswint_ai, keffn_top, & - snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai - use ice_arrays_column, only: snowfracn + snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai, update_ocn_f + use ice_arrays_column, only: snowfracn, Cdn_atm use ice_history_shared ! almost everything use ice_history_write, only: ice_write_hist 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_drag, only: accum_hist_drag + use icepack_mushy_physics, only: density_brine, liquid_fraction, temperature_mush use ice_state ! almost everything use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite @@ -1242,12 +1687,19 @@ subroutine accum_hist (dt) sn ! temporary variable for salinity real (kind=dbl_kind), dimension (nx_block,ny_block) :: & - worka, workb + worka, workb, ravgip + + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat_hist) :: & + ravgipn, worka3 real (kind=dbl_kind) :: awtvdr, awtidr, awtvdf, awtidf, puny, secday + 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 - integer (kind=int_kind) :: nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY, nt_Tsfc + 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 type (block) :: & this_block ! block information for current block @@ -1255,11 +1707,14 @@ subroutine accum_hist (dt) call icepack_query_parameters(awtvdr_out=awtvdr, awtidr_out=awtidr, & awtvdf_out=awtvdf, awtidf_out=awtidf, puny_out=puny, secday_out=secday) - call icepack_query_parameters(formdrag_out=formdrag, skl_bgc_out=skl_bgc) + call icepack_query_parameters(Tffresh_out=Tffresh, rhoi_out=rhoi, rhos_out=rhos, & + 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) 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_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) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -1328,7 +1783,8 @@ subroutine accum_hist (dt) !--------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, & - !$OMP k,n,qn,ns,sn,worka,workb,Tinz4d,Sinz4d,Tsnz4d) + !$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob, & + !$OMP worka,workb,worka3,Tinz4d,Sinz4d,Tsnz4d) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -1400,7 +1856,7 @@ subroutine accum_hist (dt) call accum_hist_field(n_frzmlt, iblk, frzmlt_init(:,:,iblk), a2D) if (f_fswfac (1:1) /= 'x') & - call accum_hist_field(n_fswfac, iblk, fswfac(:,:,iblk), a2D) + call accum_hist_field(n_fswfac, iblk, scale_factor(:,:,iblk), a2D) if (f_fswabs (1:1) /= 'x') & call accum_hist_field(n_fswabs, iblk, fswabs(:,:,iblk), a2D) @@ -1573,156 +2029,898 @@ subroutine accum_hist (dt) call accum_hist_field(n_icepresent, iblk, worka(:,:), a2D) endif - ! 3D category fields - if (f_aicen (1:1) /= 'x') & - call accum_hist_field(n_aicen-n2D, iblk, ncat_hist, & - aicen(:,:,1:ncat_hist,iblk), a3Dc) - if (f_vicen (1:1) /= 'x') & - call accum_hist_field(n_vicen-n2D, iblk, ncat_hist, & - vicen(:,:,1:ncat_hist,iblk), a3Dc) - if (f_vsnon (1:1) /= 'x') & - call accum_hist_field(n_vsnon-n2D, iblk, ncat_hist, & - vsnon(:,:,1:ncat_hist,iblk), a3Dc) - if (f_snowfracn(1:1) /= 'x') & - call accum_hist_field(n_snowfracn-n2D, iblk, ncat_hist, & - snowfracn(:,:,1:ncat_hist,iblk), a3Dc) - if (f_snowfracn(1:1) /= 'x') & - call accum_hist_field(n_snowfracn-n2D, iblk, ncat_hist, & - snowfracn(:,:,1:ncat_hist,iblk), a3Dc) - if (f_keffn_top (1:1) /= 'x') & - call accum_hist_field(n_keffn_top-n2D, iblk, ncat_hist, & - keffn_top(:,:,1:ncat_hist,iblk), a3Dc) - if (f_fsurfn_ai (1:1) /= 'x') & - call accum_hist_field(n_fsurfn_ai-n2D, iblk, ncat_hist, & - fsurfn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) - if (f_fcondtopn_ai (1:1) /= 'x') & - call accum_hist_field(n_fcondtopn_ai-n2D, iblk, ncat_hist, & - fcondtopn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) - if (f_flatn_ai (1:1) /= 'x') & - call accum_hist_field(n_flatn_ai-n2D, iblk, ncat_hist, & - flatn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) - if (f_fsensn_ai (1:1) /= 'x') & - call accum_hist_field(n_fsensn_ai-n2D, iblk, ncat_hist, & - fsensn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) - ! Calculate surface heat flux that causes melt (calculated by the - ! atmos in HadGEM3 so needed for checking purposes) - if (f_fmelttn_ai (1:1) /= 'x') & - call accum_hist_field(n_fmelttn_ai-n2D, iblk, ncat_hist, & - max(fsurfn(:,:,1:ncat_hist,iblk) - fcondtopn(:,:,1:ncat_hist,iblk),c0) & - *aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + ! 2D CMIP fields -! example for 3D field (x,y,z) -! if (f_field3dz (1:1) /= 'x') & -! call accum_hist_field(n_field3dz-n3Dccum, iblk, nzilyr, & -! field3dz(:,:,1:nzilyr,iblk), a3Dz) + if (f_sithick(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) worka(i,j) = vice(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sithick, iblk, worka(:,:), a2D) + endif - ! 4D category fields - if (f_Tinz (1:1) /= 'x') then - Tinz4d(:,:,:,:) = c0 - do n = 1, ncat_hist - do j = jlo, jhi - do i = ilo, ihi - do k = 1, nzilyr - qn = trcrn(i,j,nt_qice+k-1,n,iblk) - sn = trcrn(i,j,nt_sice+k-1,n,iblk) - Tinz4d(i,j,k,n) = icepack_ice_temperature(qn,sn) - enddo - enddo - enddo - enddo - call accum_hist_field(n_Tinz-n3Dacum, iblk, nzilyr, ncat_hist, & - Tinz4d(:,:,1:nzilyr,1:ncat_hist), a4Di) + if (f_siage(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk)*trcr(i,j,nt_iage,iblk) + enddo + enddo + call accum_hist_field(n_siage, iblk, worka(:,:), a2D) endif - if (f_Sinz (1:1) /= 'x') then - Sinz4d(:,:,:,:) = c0 - do n = 1, ncat_hist - do j = jlo, jhi - do i = ilo, ihi - if (vicen(i,j,n,iblk) > puny) then - Sinz4d(i,j,1:nzilyr,n) = trcrn(i,j,nt_sice:nt_sice+nzilyr-1,n,iblk) - endif - enddo - enddo - enddo - call accum_hist_field(n_Sinz-n3Dacum, iblk, nzilyr, ncat_hist, & - Sinz4d(:,:,1:nzilyr,1:ncat_hist), a4Di) + + if (f_sisnthick(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (vsno(i,j,iblk) > puny) & + worka(i,j) = vsno(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sisnthick, iblk, worka(:,:), a2D) endif - - if (f_Tsnz (1:1) /= 'x') then - Tsnz4d(:,:,:,:) = c0 - do n = 1, ncat_hist - do j = jlo, jhi - do i = ilo, ihi - do k = 1, nzslyr - qn = trcrn(i,j,nt_qsno+k-1,n,iblk) - Tsnz4d(i,j,k,n) = icepack_snow_temperature(qn) - enddo - enddo - enddo - enddo - call accum_hist_field(n_Tsnz-n4Dicum, iblk, nzslyr, ncat_hist, & - Tsnz4d(:,:,1:nzslyr,1:ncat_hist), a4Ds) + + if (f_sitemptop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*(trcr(i,j,nt_Tsfc,iblk)+Tffresh) + enddo + enddo + call accum_hist_field(n_sitemptop, iblk, worka(:,:), a2D) endif - - ! Calculate aggregate surface melt flux by summing category values - if (f_fmeltt_ai(1:1) /= 'x') then - do ns = 1, nstreams - if (n_fmeltt_ai(ns) /= 0) then - worka(:,:) = c0 - do j = jlo, jhi - do i = ilo, ihi - if (tmask(i,j,iblk)) then - do n=1,ncat_hist - worka(i,j) = worka(i,j) + a3Dc(i,j,n,n_fmelttn_ai(ns)-n2D,iblk) - enddo ! n - endif ! tmask - enddo ! i - enddo ! j - a2D(:,:,n_fmeltt_ai(ns),iblk) = worka(:,:) - endif - enddo - endif - !--------------------------------------------------------------- - ! accumulate other history output - !--------------------------------------------------------------- + if (f_sitempsnic(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (vsno(i,j,iblk) > puny .and. aice_init(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*(Tsnice(i,j,iblk)/aice_init(i,j,iblk)+Tffresh) + else + worka(i,j) = aice(i,j,iblk)*(trcr(i,j,nt_Tsfc,iblk)+Tffresh) + endif + enddo + enddo + call accum_hist_field(n_sitempsnic, iblk, worka(:,:), a2D) + endif - ! mechanical redistribution - call accum_hist_mechred (iblk) + if (f_sitempbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice_init(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*(Tbot(i,j,iblk)/aice_init(i,j,iblk)+Tffresh) + enddo + enddo + call accum_hist_field(n_sitempbot, iblk, worka(:,:), a2D) + endif - ! melt ponds - if (tr_pond) call accum_hist_pond (iblk) + if (f_siu(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk)*uvel(i,j,iblk) + enddo + enddo + call accum_hist_field(n_siu, iblk, worka(:,:), a2D) + endif - ! biogeochemistry - if (tr_aero .or. tr_brine .or. skl_bgc) call accum_hist_bgc (iblk) + if (f_siv(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk)*vvel(i,j,iblk) + enddo + enddo + call accum_hist_field(n_siv, iblk, worka(:,:), a2D) + endif - ! form drag - if (formdrag) call accum_hist_drag (iblk) + if (f_sispeed(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) worka(i,j) = aice(i,j,iblk) & + * sqrt(uvel(i,j,iblk)*uvel(i,j,iblk)+vvel(i,j,iblk)*vvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sispeed, iblk, worka(:,:), a2D) + endif - enddo ! iblk - !$OMP END PARALLEL DO + if (f_sidmasstranx(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))*dyu(i,j,iblk) & + + rhos*p5*(vsno(i+1,j,iblk)+vsno(i,j,iblk))*dyu(i,j,iblk)) & + * p5*(uvel(i,j-1,iblk)+uvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sidmasstranx, iblk, worka(:,:), a2D) + endif - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) + if (f_sidmasstrany(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))*dxu(i,j,iblk) & + + rhos*p5*(vsno(i,j+1,iblk)+vsno(i,j,iblk))*dxu(i,j,iblk)) & + * p5*(vvel(i-1,j,iblk)+vvel(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sidmasstrany, iblk, worka(:,:), a2D) + endif - !--------------------------------------------------------------- - ! Write output files at prescribed intervals - !--------------------------------------------------------------- + if (f_sistrxdtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice_init(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*(aice(i,j,iblk)*strairx(i,j,iblk)/aice_init(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sistrxdtop, iblk, worka(:,:), a2D) + endif - nstrm = nstreams - if (write_ic) nstrm = 1 + if (f_sistrydtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice_init(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*(aice(i,j,iblk)*strairy(i,j,iblk)/aice_init(i,j,iblk)) + enddo + enddo + call accum_hist_field(n_sistrydtop, iblk, worka(:,:), a2D) + endif - do ns = 1, nstrm - if (write_history(ns) .or. write_ic) then + if (f_sistrxubot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*strocnx(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sistrxubot, iblk, worka(:,:), a2D) + endif - !--------------------------------------------------------------- - ! Mask out land points and convert units - !--------------------------------------------------------------- + if (f_sistryubot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*strocny(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sistryubot, iblk, worka(:,:), a2D) + endif + + if (f_sicompstren(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) & + worka(i,j) = aice(i,j,iblk)*strength(i,j,iblk) + enddo + enddo + call accum_hist_field(n_sicompstren, iblk, worka(:,:), a2D) + endif + + if (f_sialb(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (fsw(i,j,iblk) > puny .and. aice_init(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*(fsw(i,j,iblk)-fswabs(i,j,iblk) & + * aice(i,j,iblk)/aice_init(i,j,iblk)) / fsw(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sialb, iblk, worka(:,:), a2D) + endif + + if (f_sihc(1:1) /= 'x') then + worka(:,:) = c0 + do k = 1,nzilyr + do j = jlo, jhi + do i = ilo, ihi + worka(i,j) = worka(i,j) + trcr(i,j,nt_qice+k-1,iblk)*vice(i,j,iblk)/real(nzilyr,kind=dbl_kind) + enddo + enddo + enddo + call accum_hist_field(n_sihc, iblk, worka(:,:), a2D) + endif + + if (f_sisnhc(1:1) /= 'x') then + worka(:,:) = c0 + do k = 1,nzslyr + do j = jlo, jhi + do i = ilo, ihi + worka(i,j) = worka(i,j) + trcr(i,j,nt_qsno+k-1,iblk)*vsno(i,j,iblk)/real(nzslyr,kind=dbl_kind) + enddo + enddo + enddo + call accum_hist_field(n_sisnhc, iblk, worka(:,:), a2D) + endif + + if (f_sidconcth(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = daidtt(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sidconcth, iblk, worka(:,:), a2D) + endif + + if (f_sidconcdyn(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = daidtd(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sidconcdyn, iblk, worka(:,:), a2D) + endif + + if (f_sidmassth(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = dvidtt(i,j,iblk) * rhoi + endif + enddo + enddo + call accum_hist_field(n_sidmassth, iblk, worka(:,:), a2D) + endif + + if (f_sidmassdyn(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = dvidtd(i,j,iblk) * rhoi + endif + enddo + enddo + call accum_hist_field(n_sidmassdyn, iblk, worka(:,:), a2D) + endif + + if (f_sidmassgrowthwat(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice_init(i,j,iblk) > puny) then + worka(i,j) = frazil(i,j,iblk)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmassgrowthwat, iblk, worka(:,:), a2D) + endif + + if (f_sidmassgrowthbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = congel(i,j,iblk)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmassgrowthbot, iblk, worka(:,:), a2D) + endif + + if (f_sidmasssi(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)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmasssi, iblk, worka(:,:), a2D) + endif + + if (f_sidmassevapsubl(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = evapi(i,j,iblk)*rhoi + endif + enddo + enddo + call accum_hist_field(n_sidmassevapsubl, iblk, worka(:,:), a2D) + endif + + if (f_sidmassmelttop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = meltt(i,j,iblk)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmassmelttop, iblk, worka(:,:), a2D) + endif + + if (f_sidmassmeltbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = meltb(i,j,iblk)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmassmeltbot, iblk, worka(:,:), a2D) + endif + + if (f_sidmasslat(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = meltl(i,j,iblk)*rhoi/dt + endif + enddo + enddo + call accum_hist_field(n_sidmasslat, iblk, worka(:,:), a2D) + endif + + if (f_sndmasssubl(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = evaps(i,j,iblk)*rhos + endif + enddo + enddo + call accum_hist_field(n_sndmasssubl, iblk, worka(:,:), a2D) + endif + + if (f_sndmasssnf(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fsnow(i,j,iblk)*rhos + endif + enddo + enddo + call accum_hist_field(n_sndmasssnf, iblk, worka(:,:), a2D) + endif + + 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 + endif + enddo + enddo + call accum_hist_field(n_sndmassmelt, iblk, worka(:,:), a2D) + endif + + if (f_siflswdtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (fsw(i,j,iblk) > puny .and. aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fsw(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflswdtop, iblk, worka(:,:), a2D) + endif + + if (f_siflswutop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (fsw(i,j,iblk) > puny .and. aice_init(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*(fsw(i,j,iblk)-fswabs(i,j,iblk) & + * aice(i,j,iblk)/aice_init(i,j,iblk)) + endif + enddo + enddo + call accum_hist_field(n_siflswutop, iblk, worka(:,:), a2D) + endif + + if (f_siflswdbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fswthru(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflswdbot, iblk, worka(:,:), a2D) + endif + + if (f_sifllwdtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*flw(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sifllwdtop, iblk, worka(:,:), a2D) + endif + + if (f_sifllwutop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*flwout(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sifllwutop, iblk, worka(:,:), a2D) + endif + + if (f_siflsenstop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fsens(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflsenstop, iblk, worka(:,:), a2D) + endif + + if (f_siflsensupbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fhocn(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflsensupbot, iblk, worka(:,:), a2D) + endif + + if (f_sifllatstop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*flat(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sifllatstop, iblk, worka(:,:), a2D) + endif + + if (f_siflcondtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fcondtop(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflcondtop, iblk, worka(:,:), a2D) + endif + + if (f_siflcondbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice_init(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fcondbot(i,j,iblk)/aice_init(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siflcondbot, iblk, worka(:,:), a2D) + endif + + if (f_sipr(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*frain(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sipr, iblk, worka(:,:), a2D) + endif + + if (f_sifb(1:1) /= 'x') then + worka(:,:) = c0 + rho_ice = rhoi + rho_ocn = rhow + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + if (ktherm == 2) then + rho_ocn = density_brine(sss(i,j,iblk)) + rho_ice = c0 + do k = 1, nzilyr + Tice = temperature_mush(trcr(i,j,nt_qice+k-1,iblk),trcr(i,j,nt_sice+k-1,iblk)) + Sbr = trcr(i,j,nt_sice+k-1,iblk) + phi = liquid_fraction(Tice,Sbr) + rhob = density_brine(Sbr) + rho_ice = rho_ice + min(phi*rhob+(c1-phi)*rhoi,rho_ocn) + enddo + rho_ice = rho_ice / real(nzilyr,kind=dbl_kind) + endif + worka(i,j) = ((rho_ocn-rho_ice)*vice(i,j,iblk) - rhos*vsno(i,j,iblk))/rho_ocn +! if (worka(i,j) < c0) then +! write(nu_diag,*) 'negative fb',rho_ocn,rho_ice,rhos +! write(nu_diag,*) vice(i,j,iblk),vsno(i,j,iblk) +! endif + endif + enddo + enddo + call accum_hist_field(n_sifb, iblk, worka(:,:), a2D) + endif + + if (f_siflsaltbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then +! Add in frazil flux + if (.not. update_ocn_f) then + if ( ktherm == 2) then + dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt + else + dfresh = -rhoi*frazil(i,j,iblk)/dt + endif + endif + dfsalt = ice_ref_salinity*p001*dfresh + worka(i,j) = aice(i,j,iblk)*(fsalt(i,j,iblk)+dfsalt) + endif + enddo + enddo + call accum_hist_field(n_siflsaltbot, iblk, worka(:,:), a2D) + endif + + if (f_siflfwbot(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then +! Add in frazil flux +! Add in frazil flux + if (.not. update_ocn_f) then + if ( ktherm == 2) then + dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt + else + dfresh = -rhoi*frazil(i,j,iblk)/dt + endif + endif + worka(i,j) = aice(i,j,iblk)*(fresh(i,j,iblk)+dfresh) + endif + enddo + enddo + call accum_hist_field(n_siflfwbot, iblk, worka(:,:), a2D) + endif + + if (f_siflfwdrain(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*(frain(i,j,iblk)+melts(i,j,iblk)+meltt(i,j,iblk)) + endif + enddo + enddo + call accum_hist_field(n_siflfwdrain, iblk, worka(:,:), a2D) + endif + + if (f_sidragtop(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*Cdn_atm(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_sidragtop, iblk, worka(:,:), a2D) + endif + + if (f_sirdgthick(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk)*(c1 - trcr(i,j,nt_alvl,iblk)) > puny) then + worka(i,j) = vice(i,j,iblk) * (c1 - trcr(i,j,nt_vlvl,iblk)) & + / (aice(i,j,iblk) * (c1 - trcr(i,j,nt_alvl,iblk))) + endif + enddo + enddo + call accum_hist_field(n_sirdgthick, iblk, worka(:,:), a2D) + endif + + if (f_siforcetiltx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*strtltx(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforcetiltx, iblk, worka(:,:), a2D) + endif + + if (f_siforcetilty(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*strtlty(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforcetilty, iblk, worka(:,:), a2D) + endif + + if (f_siforcecoriolx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*fm(i,j,iblk)*vvel(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforcecoriolx, iblk, worka(:,:), a2D) + endif + + if (f_siforcecorioly(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = -aice(i,j,iblk)*fm(i,j,iblk)*uvel(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforcecorioly, iblk, worka(:,:), a2D) + endif + + if (f_siforceintstrx(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*strintx(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforceintstrx, iblk, worka(:,:), a2D) + endif + + if (f_siforceintstry(1:1) /= 'x') then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (aice(i,j,iblk) > puny) then + worka(i,j) = aice(i,j,iblk)*strinty(i,j,iblk) + endif + enddo + enddo + call accum_hist_field(n_siforceintstry, iblk, worka(:,:), a2D) + endif + + ! 3D category fields + if (f_aicen (1:1) /= 'x') & + call accum_hist_field(n_aicen-n2D, iblk, ncat_hist, & + aicen(:,:,1:ncat_hist,iblk), a3Dc) + if (f_vicen (1:1) /= 'x') & + call accum_hist_field(n_vicen-n2D, iblk, ncat_hist, & + vicen(:,:,1:ncat_hist,iblk), a3Dc) + if (f_vsnon (1:1) /= 'x') & + call accum_hist_field(n_vsnon-n2D, iblk, ncat_hist, & + vsnon(:,:,1:ncat_hist,iblk), a3Dc) + if (f_snowfracn(1:1) /= 'x') & + call accum_hist_field(n_snowfracn-n2D, iblk, ncat_hist, & + snowfracn(:,:,1:ncat_hist,iblk), a3Dc) + if (f_snowfracn(1:1) /= 'x') & + call accum_hist_field(n_snowfracn-n2D, iblk, ncat_hist, & + snowfracn(:,:,1:ncat_hist,iblk), a3Dc) + if (f_keffn_top (1:1) /= 'x') & + call accum_hist_field(n_keffn_top-n2D, iblk, ncat_hist, & + keffn_top(:,:,1:ncat_hist,iblk), a3Dc) + if (f_fsurfn_ai (1:1) /= 'x') & + call accum_hist_field(n_fsurfn_ai-n2D, iblk, ncat_hist, & + fsurfn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + if (f_fcondtopn_ai (1:1) /= 'x') & + call accum_hist_field(n_fcondtopn_ai-n2D, iblk, ncat_hist, & + fcondtopn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + if (f_flatn_ai (1:1) /= 'x') & + call accum_hist_field(n_flatn_ai-n2D, iblk, ncat_hist, & + flatn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + if (f_fsensn_ai (1:1) /= 'x') & + call accum_hist_field(n_fsensn_ai-n2D, iblk, ncat_hist, & + fsensn(:,:,1:ncat_hist,iblk)*aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + ! Calculate surface heat flux that causes melt (calculated by the + ! atmos in HadGEM3 so needed for checking purposes) + if (f_fmelttn_ai (1:1) /= 'x') & + call accum_hist_field(n_fmelttn_ai-n2D, iblk, ncat_hist, & + max(fsurfn(:,:,1:ncat_hist,iblk) - fcondtopn(:,:,1:ncat_hist,iblk),c0) & + *aicen_init(:,:,1:ncat_hist,iblk), a3Dc) + + if (f_siitdconc (1:1) /= 'x') then + worka3(:,:,:) = c0 + do n = 1,ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (aicen(i,j,n,iblk) > puny) then + worka3(i,j,n) = aicen(i,j,n,iblk) + endif + enddo + enddo + enddo + call accum_hist_field(n_siitdconc-n2D, iblk, ncat_hist, worka3(:,:,:), a3Dc) + endif + + if (f_siitdthick (1:1) /= 'x') then + worka3(:,:,:) = c0 + do n = 1,ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (aicen(i,j,n,iblk) > puny) then + worka3(i,j,n) = vicen(i,j,n,iblk) + endif + enddo + enddo + enddo + call accum_hist_field(n_siitdthick-n2D, iblk, ncat_hist, worka3(:,:,:), a3Dc) + endif + + if (f_siitdsnthick (1:1) /= 'x') then + worka3(:,:,:) = c0 + do n = 1,ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (aicen(i,j,n,iblk) > puny) then + worka3(i,j,n) = vsnon(i,j,n,iblk) + endif + enddo + enddo + enddo + call accum_hist_field(n_siitdsnthick-n2D, iblk, ncat_hist, worka3(:,:,:), a3Dc) + endif + +! example for 3D field (x,y,z) +! if (f_field3dz (1:1) /= 'x') & +! call accum_hist_field(n_field3dz-n3Dccum, iblk, nzilyr, & +! field3dz(:,:,1:nzilyr,iblk), a3Dz) + + ! 4D category fields + if (f_Tinz (1:1) /= 'x') then + Tinz4d(:,:,:,:) = c0 + do n = 1, ncat_hist + do j = jlo, jhi + do i = ilo, ihi + do k = 1, nzilyr + qn = trcrn(i,j,nt_qice+k-1,n,iblk) + sn = trcrn(i,j,nt_sice+k-1,n,iblk) + Tinz4d(i,j,k,n) = icepack_ice_temperature(qn,sn) + enddo + enddo + enddo + enddo + call accum_hist_field(n_Tinz-n3Dacum, iblk, nzilyr, ncat_hist, & + Tinz4d(:,:,1:nzilyr,1:ncat_hist), a4Di) + endif + if (f_Sinz (1:1) /= 'x') then + Sinz4d(:,:,:,:) = c0 + do n = 1, ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (vicen(i,j,n,iblk) > puny) then + Sinz4d(i,j,1:nzilyr,n) = trcrn(i,j,nt_sice:nt_sice+nzilyr-1,n,iblk) + endif + enddo + enddo + enddo + call accum_hist_field(n_Sinz-n3Dacum, iblk, nzilyr, ncat_hist, & + Sinz4d(:,:,1:nzilyr,1:ncat_hist), a4Di) + endif + + if (f_Tsnz (1:1) /= 'x') then + Tsnz4d(:,:,:,:) = c0 + do n = 1, ncat_hist + do j = jlo, jhi + do i = ilo, ihi + do k = 1, nzslyr + qn = trcrn(i,j,nt_qsno+k-1,n,iblk) + Tsnz4d(i,j,k,n) = icepack_snow_temperature(qn) + enddo + enddo + enddo + enddo + call accum_hist_field(n_Tsnz-n4Dicum, iblk, nzslyr, ncat_hist, & + Tsnz4d(:,:,1:nzslyr,1:ncat_hist), a4Ds) + endif + + ! Calculate aggregate surface melt flux by summing category values + if (f_fmeltt_ai(1:1) /= 'x') then + do ns = 1, nstreams + if (n_fmeltt_ai(ns) /= 0) then + worka(:,:) = c0 + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + do n=1,ncat_hist + worka(i,j) = worka(i,j) + a3Dc(i,j,n,n_fmelttn_ai(ns)-n2D,iblk) + enddo ! n + endif ! tmask + enddo ! i + enddo ! j + a2D(:,:,n_fmeltt_ai(ns),iblk) = worka(:,:) + endif + enddo + endif + + !--------------------------------------------------------------- + ! accumulate other history output + !--------------------------------------------------------------- + + ! mechanical redistribution + call accum_hist_mechred (iblk) + + ! melt ponds + if (tr_pond) call accum_hist_pond (iblk) + + ! biogeochemistry + if (tr_aero .or. tr_brine .or. skl_bgc) call accum_hist_bgc (iblk) + + ! form drag + if (formdrag) call accum_hist_drag (iblk) + + enddo ! iblk + !$OMP END PARALLEL DO + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + !--------------------------------------------------------------- + ! Write output files at prescribed intervals + !--------------------------------------------------------------- + + nstrm = nstreams + if (write_ic) nstrm = 1 + + do ns = 1, nstrm + if (write_history(ns) .or. write_ic) then + + !--------------------------------------------------------------- + ! Mask out land points and convert units + !--------------------------------------------------------------- ravgct = c1/avgct(ns) !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, & - !$OMP n,nn,ravgctz) + !$OMP n,nn,ravgctz,ravgip,ravgipn) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -1730,6 +2928,34 @@ subroutine accum_hist (dt) jlo = this_block%jlo jhi = this_block%jhi + ! Ice fraction really needs to be on one of the history + ! streams, but in case it is not. + + if (n_aice(ns) > 0) then + do j = jlo, jhi + do i = ilo, ihi + if (a2D(i,j,n_aice(ns),iblk) > puny) then + ravgip(i,j) = c1/(a2D(i,j,n_aice(ns),iblk)) + else + ravgip(i,j) = c0 + endif + enddo ! i + enddo ! j + endif + if (n_aicen(ns) > n2D) then + do k=1,ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (a3Dc(i,j,k,n_aicen(ns)-n2D,iblk) > puny) then + ravgipn(i,j,k) = c1/(a3Dc(i,j,k,n_aicen(ns)-n2D,iblk)) + else + ravgipn(i,j,k) = c0 + endif + enddo ! i + enddo ! j + enddo ! k + endif + do n = 1, num_avail_hist_fields_2D if (avail_hist_fields(n)%vhistfreq == histfreq(ns)) then @@ -1744,6 +2970,503 @@ subroutine accum_hist (dt) enddo ! i enddo ! j + ! Only average for timesteps when ice present + if (index(avail_hist_fields(n)%vname,'sithick') /= 0) then + if (f_sithick(1:1) /= 'x' .and. n_sithick(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sithick(ns),iblk) = & + a2D(i,j,n_sithick(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sithick(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siage') /= 0) then + if (f_siage(1:1) /= 'x' .and. n_siage(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siage(ns),iblk) = & + a2D(i,j,n_siage(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siage(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sisnthick') /= 0) then + if (f_sisnthick(1:1) /= 'x' .and. n_sisnthick(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sisnthick(ns),iblk) = & + a2D(i,j,n_sisnthick(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sisnthick(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sitemptop') /= 0) then + if (f_sitemptop(1:1) /= 'x' .and. n_sitemptop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sitemptop(ns),iblk) = & + a2D(i,j,n_sitemptop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sitemptop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sitempsnic') /= 0) then + if (f_sitempsnic(1:1) /= 'x' .and. n_sitempsnic(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sitempsnic(ns),iblk) = & + a2D(i,j,n_sitempsnic(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sitempsnic(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sitempbot') /= 0) then + if (f_sitempbot(1:1) /= 'x' .and. n_sitempbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sitempbot(ns),iblk) = & + a2D(i,j,n_sitempbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sitempbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siu') /= 0) then + if (f_siu(1:1) /= 'x' .and. n_siu(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siu(ns),iblk) = & + a2D(i,j,n_siu(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siu(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siv') /= 0) then + if (f_siv(1:1) /= 'x' .and. n_siv(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siv(ns),iblk) = & + a2D(i,j,n_siv(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siv(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sistrxdtop') /= 0) then + if (f_sistrxdtop(1:1) /= 'x' .and. n_sistrxdtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sistrxdtop(ns),iblk) = & + a2D(i,j,n_sistrxdtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sistrxdtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sistrydtop') /= 0) then + if (f_sistrydtop(1:1) /= 'x' .and. n_sistrydtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sistrydtop(ns),iblk) = & + a2D(i,j,n_sistrydtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sistrydtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sistrxubot') /= 0) then + if (f_sistrxubot(1:1) /= 'x' .and. n_sistrxubot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sistrxubot(ns),iblk) = & + a2D(i,j,n_sistrxubot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sistrxubot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sistryubot') /= 0) then + if (f_sistryubot(1:1) /= 'x' .and. n_sistryubot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sistryubot(ns),iblk) = & + a2D(i,j,n_sistryubot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sistryubot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sicompstren') /= 0) then + if (f_sicompstren(1:1) /= 'x' .and. n_sicompstren(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sicompstren(ns),iblk) = & + a2D(i,j,n_sicompstren(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sicompstren(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sispeed') /= 0) then + if (f_sispeed(1:1) /= 'x' .and. n_sispeed(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sispeed(ns),iblk) = & + a2D(i,j,n_sispeed(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sispeed(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sialb') /= 0) then + if (f_sialb(1:1) /= 'x' .and. n_sialb(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sialb(ns),iblk) = & + a2D(i,j,n_sialb(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sialb(ns),iblk) = spval + if (albcnt(i,j,iblk,ns) <= puny) a2D(i,j,n_sialb(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflswdtop') /= 0) then + if (f_siflswdtop(1:1) /= 'x' .and. n_siflswdtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflswdtop(ns),iblk) = & + a2D(i,j,n_siflswdtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflswdtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflswutop') /= 0) then + if (f_siflswutop(1:1) /= 'x' .and. n_siflswutop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflswutop(ns),iblk) = & + a2D(i,j,n_siflswutop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflswutop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflswdbot') /= 0) then + if (f_siflswdbot(1:1) /= 'x' .and. n_siflswdbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflswdbot(ns),iblk) = & + a2D(i,j,n_siflswdbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflswdbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sifllwdtop') /= 0) then + if (f_sifllwdtop(1:1) /= 'x' .and. n_sifllwdtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sifllwdtop(ns),iblk) = & + a2D(i,j,n_sifllwdtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sifllwdtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sifllwutop') /= 0) then + if (f_sifllwutop(1:1) /= 'x' .and. n_sifllwutop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sifllwutop(ns),iblk) = & + a2D(i,j,n_sifllwutop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sifllwutop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflsenstop') /= 0) then + if (f_siflsenstop(1:1) /= 'x' .and. n_siflsenstop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflsenstop(ns),iblk) = & + a2D(i,j,n_siflsenstop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflsenstop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflsensupbot') /= 0) then + if (f_siflsensupbot(1:1) /= 'x' .and. n_siflsensupbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflsensupbot(ns),iblk) = & + a2D(i,j,n_siflsensupbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflsensupbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sifllatstop') /= 0) then + if (f_sifllatstop(1:1) /= 'x' .and. n_sifllatstop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sifllatstop(ns),iblk) = & + a2D(i,j,n_sifllatstop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sifllatstop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sipr') /= 0) then + if (f_sipr(1:1) /= 'x' .and. n_sipr(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sipr(ns),iblk) = & + a2D(i,j,n_sipr(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sipr(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sifb') /= 0) then + if (f_sifb(1:1) /= 'x' .and. n_sifb(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sifb(ns),iblk) = & + a2D(i,j,n_sifb(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sifb(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflcondtop') /= 0) then + if (f_siflcondtop(1:1) /= 'x' .and. n_siflcondtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflcondtop(ns),iblk) = & + a2D(i,j,n_siflcondtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflcondtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflcondbot') /= 0) then + if (f_siflcondbot(1:1) /= 'x' .and. n_siflcondbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflcondbot(ns),iblk) = & + a2D(i,j,n_siflcondbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflcondbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflsaltbot') /= 0) then + if (f_siflsaltbot(1:1) /= 'x' .and. n_siflsaltbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflsaltbot(ns),iblk) = & + a2D(i,j,n_siflsaltbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflsaltbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflfwbot') /= 0) then + if (f_siflfwbot(1:1) /= 'x' .and. n_siflfwbot(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflfwbot(ns),iblk) = & + a2D(i,j,n_siflfwbot(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflfwbot(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siflfwdrain') /= 0) then + if (f_siflfwdrain(1:1) /= 'x' .and. n_siflfwdrain(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siflfwdrain(ns),iblk) = & + a2D(i,j,n_siflfwdrain(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siflfwdrain(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sidragtop') /= 0) then + if (f_sidragtop(1:1) /= 'x' .and. n_sidragtop(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sidragtop(ns),iblk) = & + a2D(i,j,n_sidragtop(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sidragtop(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'sirdgthick') /= 0) then + if (f_sirdgthick(1:1) /= 'x' .and. n_sirdgthick(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_sirdgthick(ns),iblk) = & + a2D(i,j,n_sirdgthick(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_sirdgthick(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforcetiltx') /= 0) then + if (f_siforcetiltx(1:1) /= 'x' .and. n_siforcetiltx(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforcetiltx(ns),iblk) = & + a2D(i,j,n_siforcetiltx(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforcetiltx(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforcetilty') /= 0) then + if (f_siforcetilty(1:1) /= 'x' .and. n_siforcetilty(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforcetilty(ns),iblk) = & + a2D(i,j,n_siforcetilty(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforcetilty(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforcecoriolx') /= 0) then + if (f_siforcecoriolx(1:1) /= 'x' .and. n_siforcecoriolx(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforcecoriolx(ns),iblk) = & + a2D(i,j,n_siforcecoriolx(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforcecoriolx(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforcecorioly') /= 0) then + if (f_siforcecorioly(1:1) /= 'x' .and. n_siforcecorioly(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforcecorioly(ns),iblk) = & + a2D(i,j,n_siforcecorioly(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforcecorioly(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforceintstrx') /= 0) then + if (f_siforceintstrx(1:1) /= 'x' .and. n_siforceintstrx(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforceintstrx(ns),iblk) = & + a2D(i,j,n_siforceintstrx(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforceintstrx(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + if (index(avail_hist_fields(n)%vname,'siforceintstry') /= 0) then + if (f_siforceintstry(1:1) /= 'x' .and. n_siforceintstry(ns) /= 0) then + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + a2D(i,j,n_siforceintstry(ns),iblk) = & + a2D(i,j,n_siforceintstry(ns),iblk)*avgct(ns)*ravgip(i,j) + endif + if (ravgip(i,j) == c0) a2D(i,j,n_siforceintstry(ns),iblk) = spval + enddo ! i + enddo ! j + endif + endif + ! back out albedo/zenith angle dependence if (avail_hist_fields(n)%vname(1:6) == 'albice') then do j = jlo, jhi @@ -1809,6 +3532,7 @@ subroutine accum_hist (dt) do n = 1, num_avail_hist_fields_3Dc nn = n2D + n if (avail_hist_fields(nn)%vhistfreq == histfreq(ns)) then + do k = 1, ncat_hist do j = jlo, jhi do i = ilo, ihi @@ -1821,12 +3545,43 @@ subroutine accum_hist (dt) enddo ! i enddo ! j enddo ! k + if (index(avail_hist_fields(nn)%vname,'siitdthick') /= 0) then + if (f_siitdthick(1:1) /= 'x' .and. n_siitdthick(ns)-n2D /= 0) then + do k = 1, ncat_hist + 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) + endif + enddo ! i + enddo ! j + enddo ! k + endif + endif + if (index(avail_hist_fields(nn)%vname,'siitdsnthick') /= 0) then + if (f_siitdsnthick(1:1) /= 'x' .and. n_siitdsnthick(ns)-n2D /= 0) then + do k = 1, ncat_hist + do j = jlo, jhi + do i = ilo, ihi + if (tmask(i,j,iblk)) then + 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 + enddo ! i + enddo ! j + enddo ! k + endif endif + + endif + enddo ! n do n = 1, num_avail_hist_fields_3Dz nn = n3Dccum + n if (avail_hist_fields(nn)%vhistfreq == histfreq(ns)) then + do k = 1, nzilyr do j = jlo, jhi do i = ilo, ihi @@ -1940,6 +3695,8 @@ subroutine accum_hist (dt) if (n_sig1 (ns) /= 0) a2D(i,j,n_sig1(ns), iblk) = spval if (n_sig2 (ns) /= 0) a2D(i,j,n_sig2(ns), iblk) = spval if (n_sigP (ns) /= 0) a2D(i,j,n_sigP(ns), iblk) = spval + if (n_sistreave(ns) /= 0) a2D(i,j,n_sistreave(ns),iblk) = spval + if (n_sistremax(ns) /= 0) a2D(i,j,n_sistremax(ns),iblk) = spval if (n_mlt_onset(ns) /= 0) a2D(i,j,n_mlt_onset(ns),iblk) = spval if (n_frz_onset(ns) /= 0) a2D(i,j,n_frz_onset(ns),iblk) = spval if (n_hisnap (ns) /= 0) a2D(i,j,n_hisnap(ns), iblk) = spval @@ -1970,6 +3727,10 @@ subroutine accum_hist (dt) sig2 (i,j,iblk)*avail_hist_fields(n_sig2(ns))%cona if (n_sigP (ns) /= 0) a2D(i,j,n_sigP(ns),iblk) = & sigP (i,j,iblk)*avail_hist_fields(n_sigP(ns))%cona + if (n_sistreave(ns) /= 0) a2D(i,j,n_sistreave(ns),iblk) = & + p5*(sig1(i,j,iblk)+sig2(i,j,iblk))*avail_hist_fields(n_sistreave(ns))%cona + if (n_sistremax(ns) /= 0) a2D(i,j,n_sistremax(ns),iblk) = & + p5*(sig1(i,j,iblk)-sig2(i,j,iblk))*avail_hist_fields(n_sistremax(ns))%cona if (n_mlt_onset(ns) /= 0) a2D(i,j,n_mlt_onset(ns),iblk) = & mlt_onset(i,j,iblk) if (n_frz_onset(ns) /= 0) a2D(i,j,n_frz_onset(ns),iblk) = & diff --git a/cicecore/cicedynB/analysis/ice_history_shared.F90 b/cicecore/cicedynB/analysis/ice_history_shared.F90 index 97ad37a1b..8cf3ef0e5 100644 --- a/cicecore/cicedynB/analysis/ice_history_shared.F90 +++ b/cicecore/cicedynB/analysis/ice_history_shared.F90 @@ -235,6 +235,58 @@ module ice_history_shared f_mlt_onset = 'm', f_frz_onset = 'm', & f_iage = 'm', f_FY = 'm', & f_hisnap = 'm', f_aisnap = 'm', & + f_CMIP = 'x', & + f_sithick = 'x', f_sisnthick = 'x', & + f_siage = 'x', & + f_sitemptop = 'x', f_sitempsnic = 'x', & + f_sitempbot = 'x', f_sispeed = 'x', & + f_siu = 'x', f_siv = 'x', & + f_sidmasstranx = 'x', f_sidmasstrany = 'x', & + f_sistrxdtop = 'x', f_sistrydtop = 'x', & + f_sistrxubot = 'x', f_sistryubot = 'x', & + f_sicompstren = 'x', & + f_sialb = 'x', & + f_sihc = 'x', f_sisnhc = 'x', & + f_sidconcth = 'x', f_sidconcdyn = 'x', & + f_sidmassth = 'x', f_sidmassdyn = 'x', & + f_sidmassgrowthwat = 'x', & + f_sidmassgrowthbot = 'x', & + f_sidmasssi = 'x', & + f_sidmassevapsubl = 'x', & + f_sndmasssubl = 'x', & + f_sidmassmelttop = 'x', & + f_sidmassmeltbot = 'x', & + f_sidmasslat = 'x', & + f_sndmasssnf = 'x', & + f_sndmassmelt = 'x', & + f_siflswdtop = 'x', & + f_siflswutop = 'x', & + f_siflswdbot = 'x', & + f_sifllwdtop = 'x', & + f_sifllwutop = 'x', & + f_siflsenstop = 'x', & + f_siflsensupbot = 'x', & + f_sifllatstop = 'x', & + f_siflcondtop = 'x', & + f_siflcondbot = 'x', & + f_sipr = 'x', & + f_sifb = 'x', & + f_siflsaltbot = 'x', & + f_siflfwbot = 'x', & + f_siflfwdrain = 'x', & + f_siforcetiltx = 'x', & + f_siforcetilty = 'x', & + f_siforcecoriolx = 'x', & + f_siforcecorioly = 'x', & + f_siforceintstrx = 'x', & + f_siforceintstry = 'x', & + f_siitdconc = 'x', & + f_siitdthick = 'x', & + f_siitdsnthick = 'x', & + f_sidragtop = 'x', & + f_sirdgthick = 'x', & + f_sistreave = 'x', & + f_sistremax = 'x', & f_aicen = 'x', f_vicen = 'x', & f_vsnon = 'x', & f_trsig = 'm', f_icepresent = 'm', & @@ -323,6 +375,58 @@ module ice_history_shared f_mlt_onset, f_frz_onset, & f_iage, f_FY , & f_hisnap, f_aisnap , & + f_CMIP, & + f_sithick, f_sisnthick, & + f_siage, & + f_sitemptop, f_sitempsnic,& + f_sitempbot, f_sispeed, & + f_siu, f_siv, & + f_sidmasstranx, f_sidmasstrany, & + f_sistrxdtop, f_sistrydtop, & + f_sistrxubot, f_sistryubot, & + f_sicompstren, & + f_sialb, & + f_sihc, f_sisnhc, & + f_sidconcth, f_sidconcdyn,& + f_sidmassth, f_sidmassdyn,& + f_sidmassgrowthwat, & + f_sidmassgrowthbot, & + f_sidmasssi, & + f_sidmassevapsubl, & + f_sndmasssubl, & + f_sidmassmelttop, & + f_sidmassmeltbot, & + f_sidmasslat, & + f_sndmasssnf, & + f_sndmassmelt, & + f_siflswdtop, & + f_siflswutop, & + f_siflswdbot, & + f_sifllwdtop, & + f_sifllwutop, & + f_siflsenstop, & + f_siflsensupbot, & + f_sifllatstop, & + f_siflcondtop, & + f_siflcondbot, & + f_sipr, & + f_sifb, & + f_siflsaltbot, & + f_siflfwbot, & + f_siflfwdrain, & + f_siforcetiltx, & + f_siforcetilty, & + f_siforcecoriolx, & + f_siforcecorioly, & + f_siforceintstrx, & + f_siforceintstry, & + f_siitdconc, & + f_siitdthick, & + f_siitdsnthick, & + f_sidragtop, & + f_sirdgthick, & + f_sistreave, & + f_sistremax, & f_aicen, f_vicen , & f_vsnon, & f_trsig, f_icepresent,& @@ -428,6 +532,57 @@ module ice_history_shared n_dagedtt , n_dagedtd , & n_mlt_onset , n_frz_onset , & n_hisnap , n_aisnap , & + n_sithick , n_sisnthick , & + n_siage, & + n_sitemptop , n_sitempsnic , & + n_sitempbot , n_sispeed, & + n_siu, n_siv, & + n_sidmasstranx, n_sidmasstrany, & + n_sistrxdtop, n_sistrydtop, & + n_sistrxubot, n_sistryubot, & + n_sicompstren, & + n_sialb, & + n_sihc , n_sisnhc, & + n_sidconcth , n_sidconcdyn, & + n_sidmassth , n_sidmassdyn, & + n_sidmassgrowthwat, & + n_sidmassgrowthbot, & + n_sidmasssi, & + n_sidmassevapsubl, & + n_sndmasssubl, & + n_sidmassmelttop, & + n_sidmassmeltbot, & + n_sidmasslat, & + n_sndmasssnf, & + n_sndmassmelt, & + n_siflswdtop, & + n_siflswutop, & + n_siflswdbot, & + n_sifllwdtop, & + n_sifllwutop, & + n_siflsenstop, & + n_siflsensupbot, & + n_sifllatstop, & + n_siflcondtop, & + n_siflcondbot, & + n_sipr, & + n_sifb, & + n_siflsaltbot, & + n_siflfwbot, & + n_siflfwdrain, & + n_siforcetiltx, & + n_siforcetilty, & + n_siforcecoriolx, & + n_siforcecorioly, & + n_siforceintstrx, & + n_siforceintstry, & + n_siitdconc, & + n_siitdthick, & + n_siitdsnthick, & + n_sidragtop, & + n_sirdgthick, & + n_sistreave, & + n_sistremax, & n_trsig , n_icepresent , & n_iage , n_FY , & n_fsurf_ai , & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index f4ea9d7c2..17ffa388d 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -475,6 +475,10 @@ subroutine dyn_prep2 (nx_block, ny_block, & integer (kind=int_kind) :: & i, j, ij +#ifdef coupled + real (kind=dbl_kind) :: gravit +#endif + logical (kind=log_kind), dimension(nx_block,ny_block) :: & iceumask_old ! old-time iceumask @@ -586,6 +590,13 @@ subroutine dyn_prep2 (nx_block, ny_block, & ! Define variables for momentum equation !----------------------------------------------------------------- +#ifdef coupled + call icepack_query_parameters(gravit_out=gravit) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) +#endif + do ij = 1, icellu i = indxui(ij) j = indxuj(ij) diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index 51e453d43..dbdab39b3 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/cicecore/cicedynB/general/ice_flux.F90 @@ -175,7 +175,9 @@ module ice_flux Tref , & ! 2m atm reference temperature (K) Qref , & ! 2m atm reference spec humidity (kg/kg) Uref , & ! 10m atm reference wind speed (m/s) - evap ! evaporative water flux (kg/m^2/s) + evap , & ! evaporative water flux (kg/m^2/s) + evaps , & ! evaporative water flux over snow (kg/m^2/s) + evapi ! evaporative water flux over ice (kg/m^2/s) ! albedos aggregated over categories (if calc_Tsfc) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks), public :: & @@ -218,7 +220,6 @@ module ice_flux real (kind=dbl_kind), & dimension (nx_block,ny_block,max_blocks), public :: & - fswfac , & ! for history scale_factor! scaling factor for shortwave components logical (kind=log_kind), public :: & @@ -259,7 +260,10 @@ module ice_flux real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: & fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2) fcondtop,&! top surface conductive flux (W/m^2) + fcondbot,&! bottom surface conductive flux (W/m^2) fbot, & ! heat flux at bottom surface of ice (excluding excess) (W/m^2) + Tbot, & ! temperature at bottom surface of ice (deg C) + Tsnice, & ! temperature at snow ice interface (deg C) congel, & ! basal ice growth (m/step-->cm/day) frazil, & ! frazil ice growth (m/step-->cm/day) snoice, & ! snow-ice formation (m/step-->cm/day) @@ -279,6 +283,7 @@ module ice_flux dimension (nx_block,ny_block,ncat,max_blocks), public :: & fsurfn, & ! category fsurf fcondtopn,& ! category fcondtop + fcondbotn,& ! category fcondbot fsensn, & ! category sensible heat flux flatn ! category latent heat flux @@ -443,7 +448,7 @@ subroutine init_coupler_flux uocn (:,:,:) = c0 ! surface ocean currents (m/s) vocn (:,:,:) = c0 frzmlt(:,:,:) = c0 ! freezing/melting potential (W/m^2) -! frzmlt_init(:,:,:) = c0 ! freezing/melting potential (W/m^2) + frzmlt_init(:,:,:) = c0 ! freezing/melting potential (W/m^2) sss (:,:,:) = 34.0_dbl_kind ! sea surface salinity (ppt) do iblk = 1, size(Tf,3) @@ -472,10 +477,12 @@ subroutine init_coupler_flux fsens (:,:,:) = c0 flat (:,:,:) = c0 fswabs (:,:,:) = c0 -! fswint_ai(:,:,:) = c0 + fswint_ai(:,:,:) = c0 flwout (:,:,:) = -stefan_boltzmann*Tffresh**4 ! in case atm model diagnoses Tsfc from flwout evap (:,:,:) = c0 + evaps (:,:,:) = c0 + evapi (:,:,:) = c0 Tref (:,:,:) = c0 Qref (:,:,:) = c0 Uref (:,:,:) = c0 @@ -492,7 +499,7 @@ subroutine init_coupler_flux strocnyT(:,:,:) = c0 ! ice-ocean stress, y-direction (T-cell) fresh (:,:,:) = c0 fsalt (:,:,:) = c0 -! fpiond (:,:,:) = c0 + fpond (:,:,:) = c0 fhocn (:,:,:) = c0 fswthru (:,:,:) = c0 fresh_da(:,:,:) = c0 ! data assimilation @@ -581,7 +588,7 @@ subroutine init_flux_ocn fresh (:,:,:) = c0 fsalt (:,:,:) = c0 -! fpond (:,:,:) = c0 + fpond (:,:,:) = c0 fhocn (:,:,:) = c0 fswthru (:,:,:) = c0 faero_ocn(:,:,:,:) = c0 @@ -630,8 +637,11 @@ subroutine init_history_therm fsurf (:,:,:) = c0 fcondtop(:,:,:)= c0 + fcondbot(:,:,:)= c0 congel (:,:,:) = c0 -! fbot (:,:,:) = c0 + fbot (:,:,:) = c0 + Tbot (:,:,:) = c0 + Tsnice (:,:,:) = c0 frazil (:,:,:) = c0 snoice (:,:,:) = c0 dsnow (:,:,:) = c0 @@ -648,9 +658,9 @@ subroutine init_history_therm endif fsurfn (:,:,:,:) = c0 fcondtopn (:,:,:,:) = c0 + fcondbotn (:,:,:,:) = c0 flatn (:,:,:,:) = c0 fsensn (:,:,:,:) = c0 - fpond (:,:,:) = c0 fresh_ai (:,:,:) = c0 fsalt_ai (:,:,:) = c0 fhocn_ai (:,:,:) = c0 @@ -667,11 +677,10 @@ subroutine init_history_therm Cdn_ocn(:,:,:) = dragio Cdn_atm(:,:,:) = (vonkar/log(zref/iceruf)) & * (vonkar/log(zref/iceruf)) ! atmo drag for RASM -! Cdn_atm_ratio(:,:,:)= c0 + Cdn_atm_ratio(:,:,:)= c0 if (formdrag) then Cdn_atm_rdg (:,:,:) = c0 - Cdn_atm_ratio(:,:,:)= c0 Cdn_atm_floe(:,:,:) = c0 Cdn_atm_pond(:,:,:) = c0 Cdn_atm_skin(:,:,:) = c0 @@ -719,7 +728,7 @@ subroutine init_history_dyn sig2 (:,:,:) = c0 taubx (:,:,:) = c0 tauby (:,:,:) = c0 -! strength (:,:,:) = c0 + strength (:,:,:) = c0 strocnx (:,:,:) = c0 strocny (:,:,:) = c0 strairx (:,:,:) = c0 diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90 index 083cbd5eb..cd45c7aca 100644 --- a/cicecore/cicedynB/general/ice_step_mod.F90 +++ b/cicecore/cicedynB/general/ice_step_mod.F90 @@ -77,10 +77,10 @@ subroutine prep_radiation (iblk) call ice_timer_start(timer_sw) ! shortwave -! alvdr_init(:,:,:) = c0 -! alvdf_init(:,:,:) = c0 -! alidr_init(:,:,:) = c0 -! alidf_init(:,:,:) = c0 + alvdr_init(:,:,:) = c0 + alvdf_init(:,:,:) = c0 + alidr_init(:,:,:) = c0 + alidf_init(:,:,:) = c0 this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -140,12 +140,12 @@ subroutine step_therm1 (dt, iblk) use ice_calendar, only: yday use ice_domain, only: blocks_ice use ice_domain_size, only: ncat, nilyr, nslyr, n_aero - use ice_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fbot, & + use ice_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fbot, Tbot, Tsnice, & meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm, & wind, rhoa, potT, Qa, zlvl, strax, stray, flatn, fsensn, fsurfn, fcondtopn, & - flw, fsnow, fpond, sss, mlt_onset, frz_onset, & + flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, & frain, Tair, strairxT, strairyT, fsurf, fcondtop, fsens, & - flat, fswabs, flwout, evap, Tref, Qref, Uref, fresh, fsalt, fhocn, & + flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, Uref, fresh, fsalt, fhocn, & fswthru, meltt, melts, meltb, congel, snoice, & flatn_f, fsensn_f, fsurfn_f, fcondtopn_f use ice_flux_bgc, only: dsnown, faero_atm, faero_ocn @@ -317,11 +317,13 @@ subroutine step_therm1 (dt, iblk) sss (i,j, iblk), Tf (i,j, iblk), & strocnxT (i,j, iblk), strocnyT (i,j, iblk), & fbot (i,j, iblk), & + Tbot (i,j, iblk), Tsnice (i,j, iblk), & frzmlt (i,j, iblk), rside (i,j, iblk), & fsnow (i,j, iblk), frain (i,j, iblk), & fpond (i,j, iblk), & fsurf (i,j, iblk), fsurfn (i,j,:,iblk), & fcondtop (i,j, iblk), fcondtopn (i,j,:,iblk), & + fcondbot (i,j, iblk), fcondbotn (i,j,:,iblk), & fswsfcn (i,j,:,iblk), fswintn (i,j,:,iblk), & fswthrun (i,j,:,iblk), fswabs (i,j, iblk), & flwout (i,j, iblk), & @@ -330,6 +332,7 @@ subroutine step_therm1 (dt, iblk) fsens (i,j, iblk), fsensn (i,j,:,iblk), & flat (i,j, iblk), flatn (i,j,:,iblk), & evap (i,j, iblk), & + evaps (i,j, iblk), evapi (i,j, iblk), & fresh (i,j, iblk), fsalt (i,j, iblk), & fhocn (i,j, iblk), fswthru (i,j, iblk), & flatn_f (i,j,:,iblk), fsensn_f (i,j,:,iblk), & @@ -478,7 +481,7 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset) use ice_blocks, only: nx_block, ny_block use ice_domain, only: nblocks use ice_domain_size, only: ncat - use ice_grid, only: tmask +! use ice_grid, only: tmask use ice_state, only: aicen, trcrn, vicen, vsnon, & aice, trcr, vice, vsno, aice0, trcr_depend, & bound_state, trcr_base, nt_strata, n_trcr_strata @@ -984,14 +987,6 @@ subroutine ocean_mixed_layer (dt, iblk) indxi(:) = 0 indxj(:) = 0 -! this_block = get_block(blocks_ice(iblk),iblk) -! ilo = this_block%ilo -! ihi = this_block%ihi -! jlo = this_block%jlo -! jhi = this_block%jhi - -! do j = jlo, jhi -! do i = ilo, ihi do j = 1, ny_block do i = 1, nx_block diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 8f41d8360..dba4a9c00 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -1945,8 +1945,8 @@ subroutine gridbox_corners ! (1) SW corner, (2) SE corner, (3) NE corner, (4) NW corner !------------------------------------------------------------- -! latu_bounds(:,:,:,:) = c0 -! lonu_bounds(:,:,:,:) = c0 + latu_bounds(:,:,:,:) = c0 + lonu_bounds(:,:,:,:) = c0 !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) do iblk = 1, nblocks diff --git a/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 index e72f050ce..73540f5a9 100644 --- a/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 @@ -149,6 +149,7 @@ subroutine ice_write_hist(ns) .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & + .or. n==n_sistreave(ns) .or. n==n_sistremax(ns) & .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) & .or. n==n_hisnap(ns) .or. n==n_aisnap(ns)) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & diff --git a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 index 855e607a3..f0b1579f0 100644 --- a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 @@ -505,6 +505,8 @@ subroutine ice_write_hist (ns) if (hist_avg) then if (TRIM(avail_hist_fields(n)%vname)/='sig1' & .or.TRIM(avail_hist_fields(n)%vname)/='sig2' & + .or.TRIM(avail_hist_fields(n)%vname)/='sistreave' & + .or.TRIM(avail_hist_fields(n)%vname)/='sistremax' & .or.TRIM(avail_hist_fields(n)%vname)/='sigP') then status = nf90_put_att(ncid,varid,'cell_methods','time: mean') if (status /= nf90_noerr) call abort_ice(subname// & @@ -516,6 +518,7 @@ subroutine ice_write_hist (ns) .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & + .or. n==n_sistreave(ns) .or. n==n_sistremax(ns) & .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) & .or. n==n_hisnap(ns) .or. n==n_aisnap(ns)) then status = nf90_put_att(ncid,varid,'time_rep','instantaneous') diff --git a/cicecore/cicedynB/infrastructure/io/io_pio/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_pio/ice_history_write.F90 index 4372ab213..ff3d2d6fd 100644 --- a/cicecore/cicedynB/infrastructure/io/io_pio/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_pio/ice_history_write.F90 @@ -436,6 +436,8 @@ subroutine ice_write_hist (ns) if (hist_avg .and. histfreq(ns) /= '1') then if (TRIM(avail_hist_fields(n)%vname)/='sig1' & .or.TRIM(avail_hist_fields(n)%vname)/='sig2' & + .or.TRIM(avail_hist_fields(n)%vname)/='sistreave' & + .or.TRIM(avail_hist_fields(n)%vname)/='sistremax' & .or.TRIM(avail_hist_fields(n)%vname)/='sigP') then status = pio_put_att(File,varid,'cell_methods','time: mean') endif @@ -445,6 +447,7 @@ subroutine ice_write_hist (ns) .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & + .or. n==n_sistreave(ns) .or. n==n_sistremax(ns) & .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) & .or. n==n_hisnap(ns) .or. n==n_aisnap(ns)) then status = pio_put_att(File,varid,'time_rep','instantaneous') diff --git a/cicecore/drivers/cice/CICE_RunMod.F90 b/cicecore/drivers/cice/CICE_RunMod.F90 index 05cb998da..fe1aa14c4 100644 --- a/cicecore/drivers/cice/CICE_RunMod.F90 +++ b/cicecore/drivers/cice/CICE_RunMod.F90 @@ -420,17 +420,15 @@ subroutine coupling_prep (iblk) enddo enddo -! this_block = get_block(blocks_ice(iblk),iblk) -! ilo = this_block%ilo -! ihi = this_block%ihi -! jlo = this_block%jlo -! jhi = this_block%jhi + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi do n = 1, ncat -! do j = jlo, jhi -! do i = ilo, ihi - do j = 1, ny_block - do i = 1, nx_block + do j = jlo, jhi + do i = ilo, ihi if (aicen(i,j,n,iblk) > puny) then alvdf(i,j,iblk) = alvdf(i,j,iblk) & diff --git a/icepack b/icepack index 49a1cc9b7..e922ad03f 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 49a1cc9b76b05034cd12b404f93447bfdf8c355b +Subproject commit e922ad03ff32018ae47bfd7e4d73d1554c29ecfd