diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index ba1b1de6..6320986e 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -42,7 +42,7 @@ subroutine icedrv_initialize ! use icedrv_diagnostics, only: icedrv_diagnostics_debug use icedrv_flux, only: init_coupler_flux, init_history_therm, & init_flux_atm_ocn - use icedrv_forcing, only: init_forcing, get_forcing, get_wave_spec + use icedrv_forcing, only: init_forcing, get_forcing, get_wave_spec, precalc_forc use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, fiso_default, init_forcing_bgc use icedrv_restart_shared, only: restart use icedrv_init, only: input_data, init_state, init_grid2, init_fsd @@ -137,7 +137,11 @@ subroutine icedrv_initialize call init_forcing ! initialize forcing (standalone) if (skl_bgc .or. z_tracers) call init_forcing_bgc !cn if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice - call get_forcing(istep1) ! get forcing from data arrays + if (precalc_forc) then + call get_forcing(istep) ! precalculated arrays are indexed by istep + else + call get_forcing(istep1) ! get forcing from data arrays + endif if (tr_snow) then call icepack_init_snow ! snow aging table diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index a7bc82bd..fb784cc5 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -32,7 +32,7 @@ module icedrv_RunMod subroutine icedrv_run use icedrv_calendar, only: istep, istep1, time, dt, stop_now, calendar - use icedrv_forcing, only: get_forcing, get_wave_spec + use icedrv_forcing, only: get_forcing, get_wave_spec, precalc_forc use icedrv_forcing_bgc, only: faero_default, fiso_default, get_forcing_bgc use icedrv_flux, only: init_flux_atm_ocn use icedrv_history, only: history_format, history_close @@ -74,7 +74,11 @@ subroutine icedrv_run file=__FILE__,line= __LINE__) if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice - call get_forcing(istep1) ! get forcing from data arrays + if (precalc_forc) then + call get_forcing(istep) ! precalculated arrays are indexed by istep + else + call get_forcing(istep1) ! get forcing from data arrays + endif ! biogeochemistry forcing if (tr_iso) call fiso_default ! default values diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 321a8491..80d5eec4 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -129,6 +129,9 @@ module icedrv_flux qdp , & ! deep ocean heat flux (W/m^2), negative upward hmix ! mixed layer depth (m) + real (kind=dbl_kind), public :: & + sst_init ! initial sea surface temperature (C) + ! water isotopes real (kind=dbl_kind), dimension (nx), public :: & HDO_ocn , & ! seawater concentration of HDO (kg/kg) @@ -184,6 +187,12 @@ module icedrv_flux fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) fswthru_idf ! nir dif shortwave penetrating to ocean (W/m^2) + ! fixed ocean mixed layer properties (are overwritten by forcing data) + real (kind=dbl_kind), public :: & + sss_fixed , & ! Sea surface salinity (PSU) + qdp_fixed , & ! Deep ocean heat flux (negative upward, W/m^2) + hmix_fixed ! Mixed layer depth (m) + ! internal real (kind=dbl_kind), & @@ -485,8 +494,8 @@ subroutine init_coupler_flux uocn (:) = c0 ! surface ocean currents (m/s) vocn (:) = c0 frzmlt (:) = c0 ! freezing/melting potential (W/m^2) - sss (:) = 34.0_dbl_kind ! sea surface salinity (ppt) - sst (:) = -1.8_dbl_kind ! sea surface temperature (C) + sss (:) = sss_fixed ! sea surface salinity (ppt) + sst (:) = sst_init ! sea surface temperature (C) sstdat (:) = sst(:) ! sea surface temperature (C) ! water isotopes from ocean @@ -501,8 +510,8 @@ subroutine init_coupler_flux if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - qdp (:) = c0 ! deep ocean heat flux (W/m^2) - hmix (:) = c20 ! ocean mixed layer depth + qdp (:) = qdp_fixed ! deep ocean heat flux (W/m^2) + hmix (:) = hmix_fixed ! ocean mixed layer depth !----------------------------------------------------------------- ! fluxes sent to atmosphere diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90 index 98e390e2..fe1eabb4 100644 --- a/configuration/driver/icedrv_forcing.F90 +++ b/configuration/driver/icedrv_forcing.F90 @@ -10,6 +10,7 @@ module icedrv_forcing use icedrv_domain_size, only: nx use icedrv_calendar, only: time, nyr, dayyr, mday, month, secday use icedrv_calendar, only: daymo, daycal, dt, yday, sec + use icedrv_calendar, only: npt, use_leap_years, time0, year_init use icedrv_constants, only: nu_diag, nu_forcing, nu_open_clos use icedrv_constants, only: c0, c1, c2, c10, c100, p5, c4, c24 use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -20,14 +21,18 @@ module icedrv_forcing use icedrv_flux, only: zlvl, Tair, potT, rhoa, uatm, vatm, wind, & strax, stray, fsw, swvdr, swvdf, swidr, swidf, Qa, flw, frain, & fsnow, sst, sss, uocn, vocn, qdp, hmix, Tf, opening, closing, sstdat +#ifdef USE_NETCDF + use netcdf +#endif implicit none private + public :: init_forcing, get_forcing, interp_coeff, & interp_coeff_monthly, get_wave_spec - integer (kind=int_kind), parameter :: & - ntime = 8760 ! number of data points in time + integer (kind=int_kind) :: & + ntime ! number of data points in time integer (kind=int_kind), public :: & ycycle , & ! number of years in forcing cycle @@ -35,47 +40,45 @@ module icedrv_forcing fyear , & ! current year in forcing cycle fyear_final ! last year in cycle - real (kind=dbl_kind), dimension(ntime) :: & - fsw_data, & ! field values at temporal data points - cldf_data, & - fsnow_data, & - Tair_data, & - uatm_data, & - vatm_data, & - wind_data, & - strax_data, & - stray_data, & - rhum_data, & - Qa_data, & - rhoa_data, & - potT_data, & - flw_data, & - qdp_data, & - sst_data, & - sss_data, & - uocn_data, & - vocn_data, & - frain_data, & - swvdr_data, & - swvdf_data, & - swidr_data, & - swidf_data, & - zlvl_data, & - hmix_data + real (kind=dbl_kind), allocatable :: & + fsw_data(:), & ! field values at temporal data points + cldf_data(:), & + fsnow_data(:), & + Tair_data(:), & + uatm_data(:), & + vatm_data(:), & + wind_data(:), & + strax_data(:), & + stray_data(:), & + rhum_data(:), & + Qa_data(:), & + rhoa_data(:), & + potT_data(:), & + flw_data(:), & + qdp_data(:), & + sst_data(:), & + sss_data(:), & + uocn_data(:), & + vocn_data(:), & + frain_data(:), & + swvdr_data(:), & + swvdf_data(:), & + swidr_data(:), & + swidf_data(:), & + zlvl_data(:), & + hmix_data(:), & + open_data(:), & + clos_data(:) real (kind=dbl_kind), dimension(nx) :: & sst_temp - real (kind=dbl_kind), dimension(ntime) :: & - open_data, & - clos_data - character(char_len), public :: & atm_data_format, & ! 'bin'=binary or 'nc'=netcdf ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf bgc_data_format, & ! 'bin'=binary or 'nc'=netcdf - atm_data_type, & ! 'default', 'clim', 'CFS' - ocn_data_type, & ! 'default', 'SHEBA' + atm_data_type, & ! 'default', 'clim', 'CFS', 'MDF' + ocn_data_type, & ! 'default', 'SHEBA' 'MDF' bgc_data_type, & ! 'default', 'ISPOL', 'NICE' lateral_flux_type, & ! 'uniform_ice', 'open_water' atm_data_file, & ! atmospheric forcing data file @@ -94,8 +97,9 @@ module icedrv_forcing frcidf = 0.17_dbl_kind ! frac of incoming sw in near IR diffuse band logical (kind=log_kind), public :: & - oceanmixed_ice , & ! if true, use internal ocean mixed layer - restore_ocn ! restore sst if true + oceanmixed_ice , & ! if true, use internal ocean mixed layer + restore_ocn , & ! restore sst if true + precalc_forc ! whether to precalculate forcing real (kind=dbl_kind), public :: & trest, & ! restoring time scale (sec) @@ -120,6 +124,32 @@ subroutine init_forcing character(len=*), parameter :: subname='(init_forcing)' + ! Initialize ntime and allocate data arrays + if (precalc_forc) then + if (trim(atm_data_type(1:3)) /= 'MDF') & + call icedrv_system_abort(string=subname//& + 'precalc_forc should only be used with MDF atmosphere', & + file=__FILE__,line=__LINE__) + if (.not. ((trim(ocn_data_type(1:3)) == 'MDF') & + .or. (trim(ocn_data_type(1:7)) == 'default'))) & + call icedrv_system_abort(string=subname//& + 'precalc_forc should only be used with MDF ocean or'//& + ' default ocean', file=__FILE__,line=__LINE__) + ntime = npt + else + ntime = 8760 + endif + allocate(fsw_data(ntime), cldf_data(ntime), fsnow_data(ntime), & + Tair_data(ntime), uatm_data(ntime), vatm_data(ntime), & + wind_data(ntime), strax_data(ntime), stray_data(ntime), & + rhum_data(ntime), Qa_data(ntime), rhoa_data(ntime), & + potT_data(ntime), flw_data(ntime), qdp_data(ntime), & + sst_data(ntime), sss_data(ntime), uocn_data(ntime), & + vocn_data(ntime), frain_data(ntime), swvdr_data(ntime), & + swvdf_data(ntime), swidr_data(ntime), swidf_data(ntime), & + zlvl_data(ntime), hmix_data(ntime), open_data(ntime), & + clos_data(ntime)) + fyear = fyear_init + mod(nyr-1,ycycle) ! current year fyear_final = fyear_init + ycycle - 1 ! last year in forcing cycle @@ -153,6 +183,7 @@ subroutine init_forcing fsnow_data(:) = fsnow(i) ! snowfall rate (kg/m^2 s) qdp_data(:) = qdp (i) ! deep ocean heat flux (W/m^2) sss_data(:) = sss (i) ! sea surface salinity + hmix_data(:)= hmix (i) ! ocean mixed layer depth (m) uocn_data(:) = uocn (i) ! ocean current components (m/s) vocn_data(:) = vocn (i) cldf_data(:) = c0 ! cloud fraction @@ -161,6 +192,7 @@ subroutine init_forcing if (trim(atm_data_type(1:4)) == 'clim') call atm_climatological if (trim(atm_data_type(1:5)) == 'ISPOL') call atm_ISPOL if (trim(atm_data_type(1:4)) == 'NICE') call atm_NICE + if (trim(atm_data_type(1:6)) == 'MDF') call atm_MDF if (trim(ocn_data_type(1:5)) == 'SHEBA') call ice_open_clos if (restore_ocn) then @@ -176,6 +208,7 @@ subroutine init_forcing if (trim(ocn_data_type(1:5)) == 'ISPOL') call ocn_ISPOL if (trim(ocn_data_type(1:4)) == 'NICE') call ocn_NICE + if (trim(ocn_data_type(1:6)) == 'MDF') call ocn_MDF call prepare_forcing (Tair_data, fsw_data, & cldf_data, & @@ -227,164 +260,194 @@ subroutine get_forcing(timestep) character(len=*), parameter :: subname='(get_forcing)' - if (trim(atm_data_type) == 'CFS') then - ! calculate data index corresponding to current timestep - i = mod(timestep-1,ntime)+1 ! repeat forcing cycle - mlast = i - mnext = mlast - c1intp = c1 - c2intp = c0 - - ! fill all grid boxes with the same forcing data - Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) - Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) - uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) - vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) - fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) - flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) - fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) + if (precalc_forc) then + ! Fill all grid boxes with same forcing data + Tair (:) = Tair_data(timestep) + Qa (:) = Qa_data(timestep) + uatm (:) = uatm_data(timestep) + vatm (:) = vatm_data(timestep) + fsnow(:) = fsnow_data(timestep) + flw (:) = flw_data(timestep) + fsw (:) = fsw_data(timestep) ! derived (or not otherwise set) - potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) - wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) - strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) - stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) - rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) - frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) - swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) - swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) - swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) - swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) - - elseif (trim(atm_data_type) == 'clim') then - midmonth = 15 ! assume data is given on 15th of every month - recslot = 1 ! latter half of month - if (mday < midmonth) recslot = 2 ! first half of month - if (recslot == 1) then - mlast = month - mnext = mod(month ,12) + 1 - else ! recslot = 2 - mlast = mod(month+10,12) + 1 - mnext = month - endif - call interp_coeff_monthly(recslot, c1intp, c2intp) + potT (:) = potT_data(timestep) + wind (:) = wind_data(timestep) + strax(:) = strax_data(timestep) + stray(:) = stray_data(timestep) + rhoa (:) = rhoa_data(timestep) + frain(:) = frain_data(timestep) + swvdr(:) = swvdr_data(timestep) + swvdf(:) = swvdf_data(timestep) + swidr(:) = swidr_data(timestep) + swidf(:) = swidf_data(timestep) + + ! Ocean forcing + sst_temp(:) = sst_data(timestep) + sss (:) = sss_data(timestep) + uocn (:) = uocn_data(timestep) + vocn (:) = vocn_data(timestep) + qdp (:) = qdp_data(timestep) + + else + if (trim(atm_data_type) == 'CFS') then + ! calculate data index corresponding to current timestep + i = mod(timestep-1,ntime)+1 ! repeat forcing cycle + mlast = i + mnext = mlast + c1intp = c1 + c2intp = c0 + + ! fill all grid boxes with the same forcing data + Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) + Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) + uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) + vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) + fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) + flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) + fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) + + ! derived (or not otherwise set) + potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) + wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) + strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) + stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) + rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) + frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) + swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) + swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) + swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) + swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) + + elseif (trim(atm_data_type) == 'clim') then + midmonth = 15 ! assume data is given on 15th of every month + recslot = 1 ! latter half of month + if (mday < midmonth) recslot = 2 ! first half of month + if (recslot == 1) then + mlast = month + mnext = mod(month ,12) + 1 + else ! recslot = 2 + mlast = mod(month+10,12) + 1 + mnext = month + endif + call interp_coeff_monthly(recslot, c1intp, c2intp) + + ! fill all grid boxes with the same forcing data + Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) + Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) + uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) + vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) + fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) + flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) + fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) + + ! derived (or not otherwise set) + potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) + wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) + strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) + stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) + rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) + frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) + swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) + swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) + swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) + swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) + + elseif (trim(atm_data_type) == 'ISPOL') then + + offndy = 0 ! first data record (Julian day) + offset = real(offndy,dbl_kind)*secday + dataloc = 1 ! data located at middle of interval + maxrec = 365 + recslot = 2 + recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1 + mlast = mod(recnum+maxrec-2,maxrec) + 1 + mnext = mod(recnum-1, maxrec) + 1 + call interp_coeff (recnum, recslot, secday, dataloc, & + c1intp, c2intp, offset) - ! fill all grid boxes with the same forcing data Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) - flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) - fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) - - ! derived (or not otherwise set) - potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) - wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) - strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) - stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) - rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) - frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) - swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) - swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) - swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) - swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) - - elseif (trim(atm_data_type) == 'ISPOL') then - - offndy = 0 ! first data record (Julian day) - offset = real(offndy,dbl_kind)*secday - dataloc = 1 ! data located at middle of interval - maxrec = 365 - recslot = 2 - recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1 - mlast = mod(recnum+maxrec-2,maxrec) + 1 - mnext = mod(recnum-1, maxrec) + 1 - call interp_coeff (recnum, recslot, secday, dataloc, & - c1intp, c2intp, offset) - - Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) - Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) - uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) - vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) - fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) - ! derived (or not otherwise set) - potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) - wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) - strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) - stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) - rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) - frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) - - sec6hr = secday/c4; ! seconds in 6 hours - offndy = 0 - maxrec = 1460 - recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr) - recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1 ! data begins on 16 June 2004 - recslot = 2 - mlast = mod(recnum+maxrec-2,maxrec) + 1 - mnext = mod(recnum-1, maxrec) + 1 - call interp_coeff (recnum, recslot, sec6hr, dataloc, & - c1intp, c2intp, offset) + ! derived (or not otherwise set) + potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) + wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) + strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) + stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) + rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) + frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) + + sec6hr = secday/c4; ! seconds in 6 hours + offndy = 0 + maxrec = 1460 + recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr) + recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1 ! data begins on 16 June 2004 + recslot = 2 + mlast = mod(recnum+maxrec-2,maxrec) + 1 + mnext = mod(recnum-1, maxrec) + 1 + call interp_coeff (recnum, recslot, sec6hr, dataloc, & + c1intp, c2intp, offset) - fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) - flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) + fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) + flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) - ! derived - swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) - swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) - swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) - swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) + ! derived + swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) + swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) + swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) + swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) - elseif (trim(atm_data_type) == 'NICE') then + elseif (trim(atm_data_type) == 'NICE') then - offndy = 0 ! first data record (Julian day) - offset = real(offndy,dbl_kind)*secday - dataloc = 1 ! data located in middle of interval - maxrec = 365 - recslot = 2 - recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1 - mlast = mod(recnum+maxrec-2,maxrec) + 1 - mnext = mod(recnum-1, maxrec) + 1 - call interp_coeff (recnum, recslot, secday, dataloc, & - c1intp, c2intp, offset) + offndy = 0 ! first data record (Julian day) + offset = real(offndy,dbl_kind)*secday + dataloc = 1 ! data located in middle of interval + maxrec = 365 + recslot = 2 + recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1 + mlast = mod(recnum+maxrec-2,maxrec) + 1 + mnext = mod(recnum-1, maxrec) + 1 + call interp_coeff (recnum, recslot, secday, dataloc, & + c1intp, c2intp, offset) - Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) - Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) - uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) - vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) - fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) + Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext) + Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext) + uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext) + vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext) + fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext) - ! derived (or not otherwise set) - potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) - wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) - strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) - stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) - rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) - frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) - - sec6hr = secday/c4; ! seconds in 6 hours - maxrec = 1460 - dataloc = 2 ! data located at end of interval - recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr) - recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1 - recslot = 2 - mlast = mod(recnum+maxrec-2,maxrec) + 1 - mnext = mod(recnum-1, maxrec) + 1 - call interp_coeff (recnum, recslot, sec6hr, dataloc, & - c1intp, c2intp, offset) + ! derived (or not otherwise set) + potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext) + wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext) + strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext) + stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext) + rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext) + frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext) + + sec6hr = secday/c4; ! seconds in 6 hours + maxrec = 1460 + dataloc = 2 ! data located at end of interval + recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr) + recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1 + recslot = 2 + mlast = mod(recnum+maxrec-2,maxrec) + 1 + mnext = mod(recnum-1, maxrec) + 1 + call interp_coeff (recnum, recslot, sec6hr, dataloc, & + c1intp, c2intp, offset) - fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) - flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) + fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext) + flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext) - ! derived - swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) - swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) - swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) - swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) + ! derived + swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext) + swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext) + swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext) + swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext) - endif + endif ! possible bug: is the ocean data also offset to the beginning of the field campaigns? @@ -441,6 +504,7 @@ subroutine get_forcing(timestep) qdp (:) = c1intp * qdp_data(mlast) + c2intp * qdp_data(mnext) endif + endif call finish_ocn_forcing(sst_temp) @@ -534,6 +598,10 @@ subroutine atm_climatological ! 6 W/m2 warming of mixed layer from deep ocean qdp_data(:) = -6.0 ! 2 W/m2 from deep + 4 W/m2 counteracting larger ! SH+LH with bulk transfer than in MU 71 + ! Warn that this overwrites default and namelist value + write(nu_diag,*) subname + write(nu_diag,*) 'WARNING: atm_data_type = clim overwrites '//& + 'oceanic heat flux convergence from default or namelist' end subroutine atm_climatological @@ -972,6 +1040,395 @@ subroutine atm_NICE end subroutine atm_NICE +!======================================================================= + + subroutine atm_MDF + + integer (kind=int_kind) :: & + nt, & ! timestep index for Icepack arrays + i, & ! index for forcing data arrays + bound, & ! bound for subsetting data + dimlen, & ! length of the data arrays + ncid, & ! NetCDF file id + dimid, & ! NetCDF dimension id + status, & ! NetCDF status flag + varid ! NetCDF variable id + + integer (kind=8), allocatable :: & + data_time(:) ! array for time array in forcing data + + integer (kind=8), dimension(ntime) :: & + model_time ! array for Icepack minutely time + + real (kind=dbl_kind) :: & + model_time0 ! start time for model in seconds since 1970 + + real (kind=dbl_kind), allocatable :: & + data(:) ! data array from file + + character (char_len) :: & + calendar_type, & ! data calendar type + varname + + character (char_len_long) :: & + filename, & + time_basis ! time basis for data + + integer (kind=int_kind), dimension(ntime, 2) :: & + data_sections ! 2D array for indices corresponding + ! to which data values should be averaged to + ! create the model forcing values + + real (kind=dbl_kind), parameter :: & + Gregorian_year = 365.2425, & ! days in Gregorian year per cf standard + model_miss_val = -9999.00, & ! missing value for internal use + leg4_end_time = 1596034800, &! end of leg 4 in seconds since 1970 + leg5_start_time= 1598451600 ! start of leg 5 in seconds since 1970 + + character(len=*), parameter :: subname='(atm_MDF)' + + filename = trim(data_dir)//'/MDF/'//trim(atm_data_file) + + if (atm_data_format /= 'nc') then + call icedrv_system_abort(string=subname//& + ' ERROR: only NetCDF input implemented for atm_MDF', & + file=__FILE__,line=__LINE__) + else +#ifdef USE_NETCDF + ! Open forcing file + status = nf90_open(trim(filename), nf90_nowrite, ncid) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt open netcdf file', & + file=__FILE__,line=__LINE__) + + ! Create array for the time step values in seconds since 1970 + ! CF standard calendar is Gregorian + ! May have strange behavior if dt is not an integer + model_time0 = (year_init - 1970) * Gregorian_year * 24 * 3600 + time0 + do nt = 1, ntime + model_time(nt) = int(model_time0 + dt * nt, kind=8) + enddo + + ! Warn if simulation includes leg 4-5 transition + if ((model_time(1) < leg5_start_time) .and. & + (model_time(ntime) > leg4_end_time) .and. & + (index(atm_data_file, 'MOSAiC') > 0)) then + write(nu_diag,*) subname + write(nu_diag,*) 'WARNING: Forcing may be from MOSAiC '// & + 'and time includes MOSAIC leg 4-5 repositioning. ' // & + 'If so, forcing interpolation is not valid.' + endif + + ! Read, average, and interpolate forcing data from each variable + ! Moving average forcing values into model arrays + call load_var_MDF("tas", Tair_data, ncid, model_time) + call load_var_MDF("hus", Qa_data, ncid, model_time) + call load_var_MDF("uas", uatm_data, ncid, model_time) + call load_var_MDF("vas", vatm_data, ncid, model_time) + call load_var_MDF("rlds", flw_data, ncid, model_time) + call load_var_MDF("rsds", fsw_data, ncid, model_time) + call load_var_MDF("prsn", fsnow_data, ncid, model_time) + + ! Currently no rainfall data, to do + frain_data(:) = c0 + +#else + call icedrv_system_abort(string=subname//& + ' ERROR: atm_data_format = "nc" requires USE_NETCDF', & + file=__FILE__,line=__LINE__) +#endif + endif + + end subroutine atm_MDF + +!======================================================================= + +#ifdef USE_NETCDF + subroutine MDF_average(data_var_name, model_var_arr, & + data_var_len, ncid, data_sections, model_miss_val) + + character(len=*), intent(in) :: & + data_var_name ! Name of the variable in the MDF forcing file + + real (kind=dbl_kind), dimension(ntime), intent(out) :: & + model_var_arr ! array to place averaged forcing data in + + integer (kind=int_kind), intent(in) :: & + data_var_len, & ! Size of data array in MDF forcing file + ncid ! NetCDF file id + + integer (kind=int_kind), dimension(ntime, 2) :: & + data_sections ! indices for which data values to average + + real (kind=dbl_kind), intent(in) :: & + model_miss_val ! for when there is no data in a time step + + ! Local variables + real (kind=dbl_kind), dimension(data_var_len) :: & + data_var_arr ! array for data from forcing file + + real (kind=dbl_kind) :: & + work, & ! variable for averaging + data_miss_val, & ! value of missing data + count ! counter for data to average + + integer (kind=int_kind) :: & + status, & ! NetCDF status flag + nt, & ! timestep index for Icepack arrays + i, & ! index for forcing data arrays + varid ! NetCDF variable id + + character(len=*), parameter :: subname='(MDF_average)' + + ! Allocate get data and missing value from file + status = nf90_inq_varid(ncid, trim(data_var_name), varid) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//data_var_name//' var id', & + file=__FILE__,line=__LINE__) + status = nf90_get_att(ncid, varid, "missing_value", data_miss_val) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//data_var_name//' missing value', & + file=__FILE__,line=__LINE__) + status = nf90_get_var(ncid, varid, data_var_arr) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//data_var_name//' values', & + file=__FILE__,line=__LINE__) + + ! For each model time point average non-missing data values + do nt = 1, ntime + count = 0 + work = c0 + do i = data_sections(nt, 1), data_sections(nt, 2) + if (data_var_arr(i) /= data_miss_val) then + work = work + data_var_arr(i) + count = count + 1 + endif + end do + if (count > 0) then + model_var_arr(nt) = work / count + else + model_var_arr(nt) = model_miss_val + endif + end do + + end subroutine MDF_average +#endif + +!======================================================================= + + subroutine MDF_interpolate(model_var_arr, model_miss_val) + + real (kind=dbl_kind), dimension(ntime), intent(inout) :: & + model_var_arr ! array to place averaged forcing data in + + real (kind=dbl_kind), intent(in) :: & + model_miss_val ! for when there is no data in a time step + + integer (kind=int_kind) :: & + mlast, & ! index of last present data + nt, m, & ! model timestep indices + count ! counter for missing values + + character(len=*), parameter :: subname='(MDF_interpolate)' + + ! Interpolate, extrapolate for first and last values + if (model_var_arr(1) == model_miss_val) then + mlast = 0 + else + mlast = 1 + endif + do nt = 2, ntime + if (model_var_arr(nt) == model_miss_val) then + ! Do nothing (i.e., allow nt to increment) unless we're at end + if (nt == ntime) then + do m = mlast + 1, nt + model_var_arr(m) = model_var_arr(mlast) + end do + endif + else if ((nt - mlast) == 1) then + ! No missing data, increment mlast + mlast = nt + else + ! If we're at the start extrapolate to fill + if (mlast==0) then + do m = mlast + 1, nt - 1 + model_var_arr(m) = model_var_arr(nt) + end do + else + ! Interpolate missing data + do m = mlast + 1, nt - 1 + model_var_arr(m) = model_var_arr(mlast) & + + (model_var_arr(nt) - model_var_arr(mlast)) & + * (m - mlast) / (nt - mlast) + end do + endif + mlast = nt + endif + end do + + end subroutine MDF_interpolate +!======================================================================= + + subroutine load_var_MDF(data_var_name, model_var_arr, ncid, & + model_time) + + character(len=*), intent(in) :: & + data_var_name ! Name of the variable in the MDF forcing file + + real (kind=dbl_kind), dimension(ntime), intent(out) :: & + model_var_arr ! array to place forcing data in + + integer (kind=int_kind), intent(in) :: & + ncid ! NetCDF file id + + integer (kind=8), dimension(ntime), intent(in) :: & + model_time ! model time array + + ! Local variables + + integer (kind=int_kind) :: & + nt, & ! timestep index for Icepack arrays + i, & ! index for forcing data arrays + bound, & ! bound for subsetting data + dimlen, & ! length of the data arrays + dimid, & ! NetCDF dimension id + status, & ! NetCDF status flag + nvardims,& ! number of dimensions for variable + varid ! NetCDF variable id + + integer (kind=8), allocatable :: & + data_time(:) ! array for time array in forcing data + + integer, dimension(1) :: & + vardimids ! dimension id for variable + + character (char_len) :: & + calendar_type, & ! data calendar type + dimname ! name for variables dimension + + character (char_len_long) :: & + time_basis ! time basis for data + + integer (kind=int_kind), dimension(ntime, 2) :: & + data_sections ! 2D array for indices corresponding + ! to which data values should be averaged to + ! create the model forcing values + + real (kind=dbl_kind), parameter :: & + Gregorian_year = 365.2425, & ! days in Gregorian year per cf standard + model_miss_val = -9999.00 ! missing value for internal use + + character(len=*), parameter :: subname='(load_var_MDF)' + +#ifdef USE_NETCDF + + ! Get varid and missing value from file + status = nf90_inq_varid(ncid, trim(data_var_name), varid) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//data_var_name//' var id', & + file=__FILE__,line=__LINE__) + ! Get information about the dimension for this variable + status = nf90_inquire_variable(ncid, varid, ndims=nvardims) + if (nvardims /= 1) call icedrv_system_abort(& + string=subname//data_var_name//' has more than 1 dimension', & + file=__FILE__,line=__LINE__) + status = nf90_inquire_variable(ncid, varid, dimids=vardimids) + status = nf90_inquire_dimension(ncid, vardimids(1), & + name=dimname, len=dimlen) + ! Check that dimname matches pattern + if (dimname(1:4) /= 'time') call icedrv_system_abort(& + string=subname//data_var_name//' dimension name is not timeXXXX', & + file=__FILE__,line=__LINE__) + if (verify(trim(dimname(5:)), "0123456789") /= 0) call icedrv_system_abort(& + string=subname//data_var_name//' dimension name is not timeXXXX', & + file=__FILE__,line=__LINE__) + allocate (data_time(dimlen)) + ! Check that cadence variable exists and calendars match + status = nf90_inq_varid(ncid, trim(dimname), varid) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//trim(dimname)//' var id', & + file=__FILE__,line=__LINE__) + status = nf90_get_att(ncid, varid, "calendar", calendar_type) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get calendar attribute', & + file=__FILE__,line=__LINE__) + ! In future this check could be replaced with calendar matching + if (calendar_type /= "standard" .or. .not. use_leap_years) then + call icedrv_system_abort(& + string=subname//'Forcing calendar not standard or not using leap years',& + file=__FILE__,line=__LINE__) + endif + ! Get the time array + !! Note, in the file the value is actually unsigned, need to make sure this + ! doesn't cause issues since Fortran 90 doesn't support unsigned ints. + status = nf90_get_var(ncid, varid, data_time) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt read '//trim(dimname)//' values', & + file=__FILE__,line=__LINE__) + ! Convert the data time from minutes into seconds for compatability w/ icepack + data_time = data_time * 60 + + ! Check that the time basis in forcing file is correct + status = nf90_get_att(ncid, varid, "units", time_basis) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt get '//trim(dimname)//' units', & + file=__FILE__,line=__LINE__) + if (time_basis /= "minutes since 1970-01-01 00:00:00") then + call icedrv_system_abort(& + string=subname//'Time basis is not minutes since 1970',& + file=__FILE__,line=__LINE__) + endif + + ! Check that we are not extrapolating forcing outside of time bounds + if (model_time(1) < data_time(1)) call icedrv_system_abort(& + string=subname//'Simulation starts before forcing',& + file=__FILE__,line=__LINE__) + if (model_time(ntime) > data_time(dimlen)) call icedrv_system_abort(& + string=subname//'Simulation ends after forcing',& + file=__FILE__,line=__LINE__) + + ! data_sections is a 2D array where the first dimension + ! is the same length as model_time. The 1D array at each index + ! contains the start and stop indices of the data to be averaged + ! into each model timestep + ! Get the first start index + bound = model_time(1) - (model_time(2) - model_time(1))/2 + i = 1 + do while (data_time(i) < bound) + i = i + 1 + end do + data_sections(1, 1) = i + do nt = 1, ntime - 1 + ! Bound is halfway between this time step and the next + bound = (model_time(nt + 1) + model_time(nt))/2 + do while (data_time(i) < bound) + i = i + 1 + end do ! i - 1 is now the last element in timestep nt + data_sections(nt, 2) = i - 1 + data_sections(nt + 1, 1) = i + end do + ! Get the last index + bound = model_time(ntime) + (model_time(ntime) - model_time(ntime - 1))/2 + i = dimlen + do while (data_time(i) > bound) + i = i - 1 + end do + data_sections(ntime, 2) = i + + ! Moving average forcing values into model arrays + call MDF_average(data_var_name, model_var_arr, dimlen, ncid, & + data_sections, model_miss_val) + ! Linearly interpolate missing values + call MDF_interpolate(model_var_arr, model_miss_val) + +#else + call icedrv_system_abort(string=subname//& + ' load_var_MDF requires USE_NETCDF', & + file=__FILE__,line=__LINE__) +#endif + + end subroutine load_var_MDF + !======================================================================= subroutine ocn_NICE @@ -1070,7 +1527,102 @@ subroutine ocn_ISPOL qdp_data (i) = qdp (i) end do - end subroutine ocn_ISPOL + end subroutine ocn_ISPOL + +!======================================================================= + + subroutine ocn_MDF + + integer (kind=int_kind) :: & + nt, & ! timestep index for Icepack arrays + i, & ! index for forcing data arrays + bound, & ! bound for subsetting data + dimlen, & ! length of the data arrays + ncid, & ! NetCDF file id + dimid, & ! NetCDF dimension id + status, & ! NetCDF status flag + varid ! NetCDF variable id + + integer (kind=8), allocatable :: & + data_time(:) ! array for time array in forcing data + + integer (kind=8), dimension(ntime) :: & + model_time ! array for Icepack minutely time + + real (kind=dbl_kind) :: & + model_time0 ! start time for model in seconds since 1970 + + real (kind=dbl_kind), allocatable :: & + data(:) ! data array from file + + character (char_len) :: & + calendar_type, & ! data calendar type + varname + + character (char_len_long) :: & + filename, & + time_basis ! time basis for data + + integer (kind=int_kind), dimension(ntime, 2) :: & + data_sections ! 2D array for indices corresponding + ! to which data values should be averaged to + ! create the model forcing values + + real (kind=dbl_kind), parameter :: & + Gregorian_year = 365.2425, & ! days in Gregorian year per cf standard + model_miss_val = -9999.00, & ! missing value for internal use + leg4_end_time = 1596034800, &! end of leg 4 in seconds since 1970 + leg5_start_time= 1598451600 ! start of leg 5 in seconds since 1970 + + character(len=*), parameter :: subname='(ocn_MDF)' + + filename = trim(data_dir)//'/MDF/'//trim(ocn_data_file) + + if (ocn_data_format /= 'nc') then + call icedrv_system_abort(string=subname//& + ' ERROR: only NetCDF input implemented for ocn_MDF', & + file=__FILE__,line=__LINE__) + else +#ifdef USE_NETCDF + ! Open forcing file + status = nf90_open(trim(filename), nf90_nowrite, ncid) + if (status /= nf90_noerr) call icedrv_system_abort(& + string=subname//'Couldnt open netcdf file', & + file=__FILE__,line=__LINE__) + + ! Create array for the time step values in seconds since 1970 + ! CF standard calendar is Gregorian + ! May have strange behavior if dt is not an integer + model_time0 = (year_init - 1970) * Gregorian_year * 24 * 3600 + time0 + do nt = 1, ntime + model_time(nt) = int(model_time0 + dt * nt, kind=8) + enddo + + ! Warn if simulation includes leg 4-5 transition + if ((model_time(1) < leg5_start_time) .and. & + (model_time(ntime) > leg4_end_time) .and. & + (index(ocn_data_file, 'MOSAiC') > 0)) then + write(nu_diag,*) subname + write(nu_diag,*) 'WARNING: Forcing may be from MOSAiC '// & + 'and time includes MOSAIC leg 4-5 repositioning. ' // & + 'If so, forcing interpolation is not valid.' + endif + + ! Read, average, and interpolate forcing data from each variable + ! Moving average forcing values into model arrays + call load_var_MDF("so", sss_data, ncid, model_time) + call load_var_MDF("mlotst", hmix_data, ncid, model_time) + call load_var_MDF("hfsot", qdp_data, ncid, model_time) + call load_var_MDF("tos", sst_data, ncid, model_time) + +#else + call icedrv_system_abort(string=subname//& + ' ERROR: ocn_data_format = "nc" requires USE_NETCDF', & + file=__FILE__,line=__LINE__) +#endif + endif + + end subroutine ocn_MDF !======================================================================= diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index f342d74f..9a3d4748 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -10,6 +10,7 @@ module icedrv_init use icedrv_constants, only: nu_diag, ice_stdout, nu_diag_out, nu_nml use icedrv_constants, only: c0, c1, c2, c3, p2, p5 use icedrv_domain_size, only: nx + use icedrv_flux, only: sst_init use icepack_intfc, only: icepack_init_parameters use icepack_intfc, only: icepack_init_fsd use icepack_intfc, only: icepack_init_tracer_flags @@ -43,6 +44,12 @@ module icedrv_init lmask_n, & ! northern hemisphere mask lmask_s ! southern hemisphere mask + real (kind=dbl_kind), public :: & + hi_init_slab, & ! initial ice thickness for slab cell (nx=2) + hsno_init_slab, & ! initial snow thickness for slab cell (nx=2) + hbar_init_itd, & ! hbar for ice thickness for itd cell (nx=3) + hsno_init_itd ! initial snow thickness for itd cell (nx=3) + !======================================================================= contains @@ -66,6 +73,7 @@ subroutine input_data use icedrv_restart_shared, only: restart, restart_dir, restart_file, restart_format use icedrv_flux, only: l_mpond_fresh, cpl_bgc use icedrv_flux, only: default_season + use icedrv_flux, only: sss_fixed, qdp_fixed, hmix_fixed use icedrv_forcing, only: precip_units, fyear_init, ycycle use icedrv_forcing, only: atm_data_type, ocn_data_type, bgc_data_type use icedrv_forcing, only: atm_data_file, ocn_data_file, bgc_data_file @@ -74,6 +82,7 @@ subroutine input_data use icedrv_forcing, only: data_dir use icedrv_forcing, only: oceanmixed_ice, restore_ocn, trestore use icedrv_forcing, only: snw_ssp_table, lateral_flux_type + use icedrv_forcing, only: precalc_forc ! local variables @@ -138,7 +147,9 @@ subroutine input_data ice_ic, restart, restart_dir, restart_file, & restart_format, & dumpfreq, diagfreq, diag_file, cpl_bgc, & - conserv_check, history_format + conserv_check, history_format, & + hi_init_slab, hsno_init_slab, hbar_init_itd, hsno_init_itd, & + sst_init namelist /grid_nml/ & kcatbound @@ -183,7 +194,9 @@ subroutine input_data atm_data_file, ocn_data_file, bgc_data_file, & ice_data_file, & atm_data_format, ocn_data_format, bgc_data_format, & - data_dir, trestore, restore_ocn + data_dir, trestore, restore_ocn, & + sss_fixed, qdp_fixed, hmix_fixed, & + precalc_forc namelist /tracer_nml/ & tr_iage, & @@ -264,6 +277,11 @@ subroutine input_data history_format = 'none' ! if 'nc', write history files. Otherwise do nothing ice_ic = 'default' ! initial conditions are specified in the code ! otherwise, the filename for reading restarts + hi_init_slab = c2 ! initial ice thickness for slab cell (nx=2) + hsno_init_slab = c0 ! initial snow thickness for slab cell (nx=2) + hbar_init_itd = c3 ! hbar for ice thickness for itd cell (nx=3) + hsno_init_itd = 0.25_dbl_kind ! initial snow thickness for itd cell (nx=3) + sst_init = -1.8_dbl_kind ! initial mixed layer temperature (all cells) ndtd = 1 ! dynamic time steps per thermodynamic time step l_mpond_fresh = .false. ! logical switch for including meltpond freshwater ! flux feedback to ocean model @@ -290,6 +308,7 @@ subroutine input_data restore_ocn = .false. ! restore sst if true trestore = 90 ! restoring timescale, days (0 instantaneous) snw_ssp_table = 'test' ! snow table type, test or snicar + precalc_forc = .false. ! whether to precalculate forcing ! extra tracers tr_iage = .false. ! ice age @@ -673,6 +692,11 @@ subroutine input_data write(nu_diag,1030) ' restart_format = ', trim(restart_format) write(nu_diag,1030) ' history_format = ', trim(history_format) write(nu_diag,1030) ' ice_ic = ', trim(ice_ic) + write(nu_diag,1005) ' hi_init_slab = ', hi_init_slab + write(nu_diag,1005) ' hsno_init_slab = ', hsno_init_slab + write(nu_diag,1005) ' hbar_init_itd = ', hbar_init_itd + write(nu_diag,1005) ' hsno_init_itd = ', hsno_init_itd + write(nu_diag,1005) ' sst_init = ', sst_init write(nu_diag,1010) ' conserv_check = ', conserv_check write(nu_diag,1020) ' kitd = ', kitd write(nu_diag,1020) ' kcatbound = ', kcatbound @@ -778,6 +802,7 @@ subroutine input_data write(nu_diag,1030) ' ocn_data_file = ', trim(ocn_data_file) write(nu_diag,1030) ' bgc_data_file = ', trim(bgc_data_file) write(nu_diag,1030) ' ice_data_file = ', trim(ice_data_file) + write(nu_diag,1010) ' precalc_forc = ', precalc_forc if (trim(atm_data_type)=='default') & write(nu_diag,1030) ' default_season = ', trim(default_season) @@ -1342,9 +1367,6 @@ subroutine set_state_var (nx, & real (kind=dbl_kind), dimension(nslyr) :: & qsn ! snow enthalpy (J/m3) - real (kind=dbl_kind), parameter :: & - hsno_init = 0.25_dbl_kind ! initial snow thickness (m) - logical (kind=log_kind) :: tr_brine, tr_lvl, tr_fsd, tr_snow integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice, nt_fsd integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl @@ -1408,24 +1430,31 @@ subroutine set_state_var (nx, & i = 1 ! ice-free ! already initialized above + if (i <= nx) then ! need to set ocean parameters + sst(i) = sst_init + endif !----------------------------------------------------------------- - i = 2 ! 2-m slab, no snow + i = 2 ! 100% ice concentration slab, thickness and snow from namelist if (i <= nx) then - if (3 <= ncat) then - n = 3 - ainit(n) = c1 ! assumes we are using the default ITD boundaries - hinit(n) = c2 - else - ainit(ncat) = c1 - hinit(ncat) = c2 - endif + sst(i) = sst_init ! initial ocean temperature + do n = 1, ncat + if (hi_init_slab <= hin_max(n)) then + ainit(n) = c1 + hinit(n) = hi_init_slab + exit + endif + enddo + if (hi_init_slab > hin_max(ncat)) then + ainit(ncat) = c1 + hinit(ncat) = hi_init_slab + endif do n = 1, ncat ! ice volume, snow volume aicen(i,n) = ainit(n) vicen(i,n) = hinit(n) * ainit(n) ! m - vsnon(i,n) = c0 + vsnon(i,n) = hsno_init_slab * ainit(n) ! tracers call icepack_init_enthalpy(Tair = Tair(i), & Tf = Tf(i), & @@ -1469,8 +1498,9 @@ subroutine set_state_var (nx, & i = 3 ! full thickness distribution if (i <= nx) then + sst(i) = sst_init ! initial category areas in cells with ice - hbar = c3 ! initial ice thickness with greatest area + hbar = hbar_init_itd ! initial ice thickness with greatest area ! Note: the resulting average ice thickness ! tends to be less than hbar due to the ! nonlinear distribution of ice thicknesses @@ -1494,7 +1524,7 @@ subroutine set_state_var (nx, & ! ice volume, snow volume aicen(i,n) = ainit(n) vicen(i,n) = hinit(n) * ainit(n) ! m - vsnon(i,n) = min(aicen(i,n)*hsno_init,p2*vicen(i,n)) + vsnon(i,n) = min(aicen(i,n)*hsno_init_itd,p2*vicen(i,n)) ! tracers call icepack_init_enthalpy(Tair = Tair(i), & Tf = Tf(i), & diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index d831c2d5..a3fc6d25 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -7,6 +7,11 @@ npt = 8760 ndtd = 1 ice_ic = 'default' + hi_init_slab = 2.0 + hsno_init_slab = 0.0 + hbar_init_itd = 3.0 + hsno_init_itd = 0.25 + sst_init = -1.8 restart = .false. restart_dir = './restart/' restart_format = 'bin' @@ -58,7 +63,7 @@ albicev = 0.78 albicei = 0.36 albsnowv = 0.98 - albsnowi = 0.70 + albsnowi = 0.70 ahmax = 0.3 R_ice = 0. R_pnd = 0. @@ -121,6 +126,9 @@ wave_spec_type = 'none' restore_ocn = .false. trestore = 90 + sss_fixed = 34.0 + qdp_fixed = 0.0 + hmix_fixed = 20.0 precip_units = 'mks' default_season = 'spring' atm_data_type = 'clim' @@ -137,6 +145,7 @@ atm_data_format = 'bin' ocn_data_format = 'bin' bgc_data_format = 'bin' + precalc_forc = .false. / &dynamics_nml @@ -167,7 +176,7 @@ tr_bgc_PON = .false. tr_bgc_hum = .false. tr_bgc_DON = .false. - tr_bgc_Fe = .false. + tr_bgc_Fe = .false. grid_o = 0.006 grid_o_t = 0.006 l_sk = 2.0 @@ -175,17 +184,17 @@ l_skS = 0.028 phi_snow = -1.0 initbio_frac = 1.0 - frazil_scav = 0.8 - ratio_Si2N_diatoms = 1.8 + frazil_scav = 0.8 + ratio_Si2N_diatoms = 1.8 ratio_Si2N_sp = 0.0 ratio_Si2N_phaeo = 0.0 - ratio_S2N_diatoms = 0.03 - ratio_S2N_sp = 0.03 + ratio_S2N_diatoms = 0.03 + ratio_S2N_sp = 0.03 ratio_S2N_phaeo = 0.03 ratio_Fe2C_diatoms = 0.0033 ratio_Fe2C_sp = 0.0033 ratio_Fe2C_phaeo = 0.1 - ratio_Fe2N_diatoms = 0.023 + ratio_Fe2N_diatoms = 0.023 ratio_Fe2N_sp = 0.023 ratio_Fe2N_phaeo = 0.7 ratio_Fe2DON = 0.023 diff --git a/configuration/scripts/options/set_nml.atmmosaic b/configuration/scripts/options/set_nml.atmmosaic new file mode 100644 index 00000000..4e747203 --- /dev/null +++ b/configuration/scripts/options/set_nml.atmmosaic @@ -0,0 +1,14 @@ +atm_data_type = 'MDF' +atm_data_file = 'MOSAiC_kazr_snow_MDF_20191005_20201001.nc' +atm_data_format = 'nc' +use_leap_years = .true. +fyear_init = 2019 +year_init = 2019 +istep0 = 7968 +npt = 6120 +precalc_forc = .true. +hi_init_slab = 0.80 +hsno_init_slab = 0.09 +qdp_fixed = -1.0 +sss_fixed = 32.0 +hmix_fixed = 45.0 \ No newline at end of file diff --git a/configuration/scripts/options/set_nml.ocnmosaic b/configuration/scripts/options/set_nml.ocnmosaic new file mode 100644 index 00000000..d75acbe7 --- /dev/null +++ b/configuration/scripts/options/set_nml.ocnmosaic @@ -0,0 +1,6 @@ +ocn_data_type = 'MDF' +ocn_data_file = 'MOSAiC_ocn_MDF_20191006_20200919.nc' +ocn_data_format = 'nc' +ustar_min = 0.005 +precalc_forc = .true. +ice_data_file = 'unknown_ice_data_file' \ No newline at end of file diff --git a/doc/source/developer_guide/dg_driver.rst b/doc/source/developer_guide/dg_driver.rst index b8389b90..14dae9d9 100755 --- a/doc/source/developer_guide/dg_driver.rst +++ b/doc/source/developer_guide/dg_driver.rst @@ -40,7 +40,7 @@ Overview The icepack driver exists to test the column physics. At the present time, it is hardwired to run 4 different gridcells on one processor with the same forcing used for all gridcells. -There is no MPI and no threading built into the icepack driver. There is limited IO capabilities, -no history files, and no netcdf restart files. The model generally runs very quickly. +There is no MPI and no threading built into the icepack driver. There is limited IO capabilities. +The model generally runs very quickly. Forcing data and details on these data are available in :ref:`force`. diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index 712d1d89..1e0cbeaa 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -206,7 +206,9 @@ section :ref:`tabnamelist`. "H2_16O_ocn", "concentration of H2_16O isotope in ocean", "kg/kg" "H2_18O_ocn", "concentration of H2_18O isotope in ocean", "kg/kg" "HDO_ocn", "concentration of HDO isotope in ocean", "kg/kg" + "hbar_init_itd", ":math:`\bullet` initial modal ice thickness for itd-initialized grid cell", "3. m" "hfrazilmin", "minimum thickness of new frazil ice", "0.05 m" + "hi_init_slab", ":math:`\bullet` initial ice thickness for slab-initialized grid cell", "2. m" "hi_min", "minimum ice thickness for thinnest ice category", "m" "hi_ssl", "ice surface scattering layer thickness", "0.05 m" "hicen", "ice thickness in category n", "m" @@ -215,6 +217,7 @@ section :ref:`tabnamelist`. "hin_max", "category thickness limits", "m" "history_format", "turns on netcdf history output if set to 'nc'", "" "hmix", "ocean mixed layer depth", "20. m" + "hmix_fixed", ":math:`\bullet` constant ocean mixed layer depth", "20. m" "hour", "hour of the year", "" "hp0", "pond depth at which shortwave transition to bare ice occurs", "0.2 m" "hp1", ":math:`\bullet` critical ice lid thickness for topo ponds (dEdd)", "0.01 m" @@ -224,6 +227,8 @@ section :ref:`tabnamelist`. "hs0", ":math:`\bullet` snow depth at which transition to ice occurs (dEdd)", "" "hs1", ":math:`\bullet` snow depth of transition to pond ice", "0.03 m" "hs_ssl", "snow surface scattering layer thickness", "0.04 m" + "hsno_init_itd", ":math:`\bullet` initial snow depth for itd-initialized grid cell", "0.25 m" + "hsno_init_slab", ":math:`\bullet` initial snow depth for slab-initialized grid cell", "0. m" "Hstar", "determines mean thickness of ridged ice", "25. m" "**I**", "", "" "i0vis","fraction of penetrating visible solar radiation", "0.70" @@ -357,13 +362,15 @@ section :ref:`tabnamelist`. "potT", "atmospheric potential temperature", "K" "PP_net", "total primary productivity per grid cell", "mg C/m\ :math:`^2`/s" "precip_units", ":math:`\bullet` liquid precipitation data units", "" + "precalc_forc", ":math:`\bullet` if true, average/interpolate forcing data on initialization", "" "print_points", ":math:`\bullet` if true, print point data", "F" "Pstar", "ice strength parameter", "2.75\ :math:`\times`\ 10\ :math:`^4`\ N/m\ :math:`^2`" "puny", "a small positive number", "1\ :math:`\times`\ 10\ :math:`^{-11}`" "**Q**", "", "" "Qa", "specific humidity at 10 m", "kg/kg" "Qa_iso", "water isotope specific humidity at 10 m", "kg/kg" - "qdp", "deep ocean heat flux", "W/m\ :math:`^2`" + "qdp", "oceanic heat flux convergence", "W/m\ :math:`^2`" + "qdp_fixed", ":math:`\bullet` constant oceanic heat flux convergence", "0. W/m\ :math:`^2`" "qqqice", "for saturated specific humidity over ice", "1.16378\ :math:`\times`\ 10\ :math:`^7`\ kg/m\ :math:`^3`" "qqqocn", "for saturated specific humidity over ocean", "6.275724\ :math:`\times`\ 10\ :math:`^6`\ kg/m\ :math:`^3`" "Qref", "2m atmospheric reference specific humidity", "kg/kg" @@ -432,7 +439,9 @@ section :ref:`tabnamelist`. "spval_dbl", "special value (double precision)", ":math:`10^{30}`", "" "ss_tltx(y)", "sea surface in the x(y) direction", "m/m" "sss", "sea surface salinity", "ppt" + "sss_fixed", ":math:`\bullet` constant mixed layer salinity", "34. ppt" "sst", "sea surface temperature", "C" + "sst_init", ":math:`\bullet` initial mixed layer temperature", "-1.8\ :math:`^\circ`\ C" "Sswabs", "shortwave radiation absorbed in snow layers", "W/m\ :math:`^2`" "stefan-boltzmann", "Stefan-Boltzmann constant", "5.67\ :math:`\times`\ 10\ :math:`^{-8}` W/m\ :math:`^2`\ K\ :math:`^4`" "stop_now", "if 1, end program execution", "" diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 2296cbd9..1ebb511f 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -843,6 +843,13 @@ @article{Edwards2012 doi = {https://doi.org/10.4319/lo.2012.57.2.0554}, year = {2012} } +@misc{Clemens-Sewall24 + author = "Clemens-Sewall, D. and Cox, C. and Schulz, K. and Raphael, I.", + title = "{Merged Datasets for the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) Central Observatory in the Arctic Ocean (2019-2020)}", + year = {2024}, + doi = {http://dx.doi.org/10.18739/A2GX44W6J}, + howpublished = {Arctic Data Center} +} % ********************************************** % For new entries, see example entry in BIB_TEMPLATE.txt % ********************************************** diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 3dea32fb..87c200f5 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -142,8 +142,12 @@ setup_nml "", "``y``", "write restart every ``dumpfreq_n`` years", "" "``dump_last``", "true/false", "write restart at end of run", "false" "``dt``", "seconds", "thermodynamics time step length", "3600." + "``hbar_init_itd``", "real", "initial modal ice thickness for itd-initialized grid cell in m", "3.0" + "``hi_init_slab``", "real", "initial ice thickness for slab-initialized grid cell in m", "2.0" "``history_format``", "``cdf``", "history file output in netcdf format", "``none``" "","``none``","no history output","" + "``hsno_init_itd``", "real", "initial snow depth for itd-initialized grid cell in m", "0.25" + "``hsno_init_slab``", "real", "initial snow depth for slab-initialized grid cell in m", "0.0" "``ice_ic``", "``default``", "latitude and sst dependent initial condition", "``default``" "", "``none``", "no ice", "" "", "'path/file'", "restart file name", "" @@ -155,6 +159,7 @@ setup_nml "``restart_file``", "string", "output file prefix for restart dump", "'iced'" "``restart_format``", "``bin``", "restart file output in binary format", "``bin``" "","``cdf``","restart file output in netcdf format","" + "``sst_init``", "real", "initial ocean mixed layer temperature in C", "-1.8" "``use_leap_years``", "logical", "include leap days", "``.false.``" "``year_init``", "integer", "the initial year if not using restart", "0" "", "", "", "" @@ -358,6 +363,7 @@ forcing_nml "``formdrag``", "logical", "calculate form drag", "``.false.``" "``fyear_init``", "integer", "first year of atmospheric forcing data", "1998" "``highfreq``", "logical", "high-frequency atmo coupling", "``.false.``" + "``hmix_fixed``", "real", "constant ocean mixed layer depth in m", "20.0" "``lateral_flux_type``", "``uniform_ice``", "flux ice with identical properties into the cell when closing (Icepack only)", "" "", "``none``", "advect open water into the cell when closing (Icepack only)", "" "``ice_data_file``", "string", "file containing ice opening, closing data", "' '" @@ -371,13 +377,16 @@ forcing_nml "", "``ISPOL``", "ISPOL experiment data (see :ref:`force`)", "" "", "``NICE``", "N-ICE experiment data (see :ref:`force`)", "" "", "``SHEBA``", "Opening/closing dataset from SHEBA", "" + "``precalc_forc``", "logical", "average/interpolate forcing data on initialization", "``.false.``" "``precip_units``", "``mks``", "liquid precipitation data units", "``mks``" "", "``mm_per_month``", "", "" "", "``mm_per_sec``", "(same as MKS units)", "" "", "``m_per_sec``", "", "" + "``qdp_fixed``", "real", "constant oceanic heat flux convergence W/m^2", "0.0" "``restore_ocn``", "logical", "restore sst to data", "``.false.``" "``saltflux_option``", "``constant``","salt flux is referenced to a constant salinity","``constant``" "","``prognostic``","use actual sea ice bulk salinity in flux" + "``sss_fixed``", "real", "constant ocean mixed layer salinity in ppt", "34.0" "``tfrz_option``","``constant``", "constant ocean freezing temperature (Tocnfrz)","``mushy``" "", "``linear_salt``", "linear function of salinity (ktherm=1)", "" "", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8^{\circ} C`)", "" diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 76a0db32..bd1e3945 100755 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -118,17 +118,74 @@ will be described in more detail in the :ref:`tabnamelist`. Two namelist variables control model initialization, ``ice_ic`` and ``restart``. Setting ``ice_ic`` = 'default' causes the model to run using -initial values set in the code. To start +initial values set in the code and the namelist. To start from a file **filename**, set ``restart`` = .true. and ``ice_ic`` = **filename**. When restarting using the Icepack driver, for simplicity the tracers are assumed to be set the same way (on/off) as in the run that created the restart file; i.e. that the restart file contains exactly the information needed for the new run. CICE is more flexible in this regard. +When the model is not running from a restart file (i.e., ``ice_ic`` = 'default' +and ``restart`` = .false.), there are namelist options that control the initial +snow depth , ice thicknesses and mixed layer temperature (``sst_init``). +For the slab-initialized grid cell (``nx`` = 2), the run starts with a single +ice thickness category having 100% ice cover . ``hsno_init_slab`` and +``hi_init_slab`` define the initial snow depth and ice thickness for that +ice thickness category. For the itd-initialized grid cell (``nx`` = 3), the ice +thickness in each category is set to the midpoint of that category's ice +thickness range (excluding the last category, which is set to 1 m thicker than +the lower bound). The area fraction of each category is set according to a +normalized, downward-facing parabolic function of ice thickness, where the +maximum of the parabola is ``hbar_init_itd`` and the area fraction of open +water is zero. All thickness categories are initialized with a snow depth of +``hsno_init_itd``. + For stand-alone runs, -routines in **icedrv\_forcing.F90** read and interpolate data from files, -and are intended merely for testing, although they can also provide guidance for -the user to write his or her own routines. +routines in **icedrv\_forcing.F90** read and interpolate data from files. +The namelist variables ``precalc_forc``, ``atm_data_type``, +``atm_data_format``, ``ocn_data_type``, and ``ocn_data_format`` control +how the forcing data is handled. If ``precalc_forc`` = .false. and the +``atm/ocn_data_type`` = 'bin', when ``init_forcing`` is called, a +subroutine for the specific dataset (e.g., ``atm_CFS``) stores the +forcing data in the ``*_data`` arrays in essentially the same format +that the raw data is present in, without timestamp information. Then, +at each timestep the ``get_forcing`` subroutine has a code block for +each forcing dataset that contains the forcing's time basis and +interpolates the forcing data to the given timestep. The forcing data +that is available in 'bin' format are intended merely for code testing, +not scientific results. + +If ``precalc_forc`` = .true., ``atm/ocn_data_type`` = 'MDF', and +``atm/ocn_data_format`` = 'nc', then ``init_forcing`` reads data from +netCDF files formatted according to the Merged Data File (MDF), +conventions which includes timestamp information. During initialization, +the forcing data is averaged/interpolated to the Icepack timestep and +stored the the ``*_data`` arrays. The ``get_forcing`` subroutine then +simply queries the ``*_data`` arrays at each timestep. MDF forcing data +are expected to come from observations, and hence may contain missing +data and may be present at a shorter or longer sampling interval than +the Icepack timestep. To handle variable frequencies and missing data, +forcing data are first temporally-averaged for each timestep and then +interpolated. For a give data variable (``var_data``), The +``MDF_average`` subroutine takes the average of all forcing datapoints +within each ``timestep`` +- 0.5 ``dt`` excluding missing values and +stores the results in ``var_data``. If there is no valid data within +0.5 ``dt`` of a given ``timestep`` (e.g., most timesteps if ``dt`` is +much smaller than the sampling interval) then a missing value is placed +in ``var_data(timestep)``. Then, the ``MDF_interpolate`` subroutine +linearly interpolates missing values in ``var_data``. The MDF +conventions were developed by the Year of Polar Prediction supersite +Model Intercomparison Project (`Uttal et al., 2024 +`_) and a `python toolbox +`_ is available to build MDF +files from raw data. The ``ocn_MDF`` subroutine currently assumes that +the oceanic heat flux convergence (``qdp``) is equal to the turbulent +heat flux over the thermocline. + +If no ocean forcing is provided, namelist variables provide constant +values of the ocean mixed layer salinity (``sss_fixed``), thickness +(``hmix_fixed``), and oceanic heat flux convergence (``qdp_fixed``). If +forcing data is provided then these variables are ignored. .. _parameters: diff --git a/doc/source/user_guide/ug_running.rst b/doc/source/user_guide/ug_running.rst index 782ef56b..337901ec 100755 --- a/doc/source/user_guide/ug_running.rst +++ b/doc/source/user_guide/ug_running.rst @@ -736,6 +736,32 @@ cases. Current filenames can be found in the options scripts in At present, only the opening and closing rates (1/s) are used from the forcing data. In the namelist, set ``ocn_data_type = SHEBA`` to use this forcing in Icepack. + d) **Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC)** + + Atmospheric and oceanic forcing are available from the 2019-2020 Multidisciplinary Drifting Observatory + for the Study of Arctic Climate (MOSAiC) expedition :cite:`Clemens-Sewall24`. The atmospheric forcing is + available minutely and the oceanic forcing is available daily. These data are based on observations from + a collection of drifting ice floes in the Arctic Ocean. MOSAiC consisted of two drift experiments, + in which the R/V Polarstern was moored to a drifting ice floe and continuous observations were made. The + first drift started in October 2019 North of the Laptev Sea (85.4N, 129.2E) and ended in July 2020 in + the Fram Strait (79.2N, 2.6W). The second drift was from August to September, 2020 in the vicinity of + the North Pole (approximately 88.6N, 110.6E). + + Atmospheric forcing fields from :cite:`Clemens-Sewall24` consist of 2-m air temperature (K), 2-m + specific humidity (kg/kg), 2-m wind velocity in the u and v directions (m/s), downward solar radiation + (:math:`W/m^2`), and downward longwave radiation (:math:`W/m^2`). For the first drift, snowfall rate + (:math:`kg/m^2/s`) is also available. Oceanic forcing fields consist of mixed layer salinity (PSU), + mixed layer depth (m), turbulent heat fluxes over the halocline and thermocline (:math:`W/m^2`), + mixed layer temperature (K), ocean-ice friction velocity (m/s), and sea ice drift speed (m/s). Note, + icepack's driver currently lacks a mechanism to use the ocean-ice friction velocity and/or sea ice + drift speed in the thermodynamics calculations. Additionally, caution should be used when using the + ocean mixed layer temperatures directly because for much of the experiment icepack's parameterization + of the ocean mixed layer freezing point differs from the true freezing point by ~0.5 C. + + MOSAiC data are available exclusively as MDF-formatted netCDF files (see :ref:`init`). + + + 3) **Climatological** - Maykut and Untersteiner 1971 :cite:`Maykut71` The climatological forcing consists of a monthly climatology of downward radiative fluxes, air temperature,