From acfc046f02c3e09f935833748e8853f393e66fc3 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Wed, 1 Mar 2023 07:58:26 -0800 Subject: [PATCH 1/4] Update implementation of some optional arguments in Icepack (#429) * Update implementation of some optional arguments in Icepack - Remove local copies where possible - Check optional arguments Remove public interface declarations where not needed Clean up some intent statements * Add write_diags argument to icepack_init_fsd_bounds Update documentation * Clean up meltsliq implementation in icepack_step_therm1 --- columnphysics/icepack_fsd.F90 | 12 +- columnphysics/icepack_itd.F90 | 10 +- columnphysics/icepack_mechred.F90 | 39 ++--- columnphysics/icepack_orbital.F90 | 15 ++ columnphysics/icepack_shortwave.F90 | 156 ++++++----------- columnphysics/icepack_therm_bl99.F90 | 59 +++---- columnphysics/icepack_therm_itd.F90 | 5 +- columnphysics/icepack_therm_vertical.F90 | 37 ++-- columnphysics/icepack_tracers.F90 | 2 +- configuration/driver/icedrv_InitMod.F90 | 3 +- doc/source/user_guide/interfaces.include | 204 ++++++++++++----------- doc/source/user_guide/lg_sequence.rst | 10 +- 12 files changed, 241 insertions(+), 311 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 4d89ca845..bee46cd35 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -120,18 +120,10 @@ subroutine icepack_init_fsd_bounds(nfsd, & real (kind=dbl_kind), dimension(:), allocatable :: & lims - logical (kind=log_kind) :: & - l_write_diags ! local write diags - character(len=8) :: c_fsd1,c_fsd2 character(len=2) :: c_nf character(len=*), parameter :: subname='(icepack_init_fsd_bounds)' - l_write_diags = .true. - if (present(write_diags)) then - l_write_diags = write_diags - endif - if (nfsd.eq.24) then allocate(lims(24+1)) @@ -230,7 +222,8 @@ subroutine icepack_init_fsd_bounds(nfsd, & c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m' enddo - if (l_write_diags) then + if (present(write_diags)) then + if (write_diags) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname @@ -244,6 +237,7 @@ subroutine icepack_init_fsd_bounds(nfsd, & write(warnstr,*) ' ' call icepack_warnings_add(warnstr) endif + endif end subroutine icepack_init_fsd_bounds diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 00f0f768f..3c18a4e13 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -997,10 +997,12 @@ subroutine cleanup_itd (dt, ntrcr, & faero_ocn(it) = faero_ocn(it) + dfaero_ocn(it) enddo endif - if (tr_iso) then - do it = 1, n_iso - fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it) - enddo + if (present(fiso_ocn)) then + if (tr_iso) then + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it) + enddo + endif endif if (present(flux_bio)) then do it = 1, nbtrcr diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 4dc87cb28..ea049761f 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -58,10 +58,7 @@ module icepack_mechred implicit none private - public :: ridge_ice, & - asum_ridging, & - ridge_itd, & - icepack_ice_strength, & + public :: icepack_ice_strength, & icepack_step_ridge real (kind=dbl_kind), parameter :: & @@ -113,7 +110,7 @@ subroutine ridge_ice (dt, ndtd, & dardg1ndt, dardg2ndt, & dvirdgndt, & araftn, vraftn, & - closing_flag,closing ) + closing ) integer (kind=int_kind), intent(in) :: & ndtd , & ! number of dynamics subcycles @@ -161,7 +158,6 @@ subroutine ridge_ice (dt, ndtd, & krdg_redist ! selects redistribution function logical (kind=log_kind), intent(in) :: & - closing_flag, &! flag if closing is valid tr_brine ! if .true., brine height differs from ice thickness ! optional history fields @@ -296,7 +292,7 @@ subroutine ridge_ice (dt, ndtd, & ! Compute the area opening and closing. !----------------------------------------------------------------- - if (closing_flag) then + if (present(opening) .and. present(closing)) then opning = opening closing_net = closing @@ -595,11 +591,13 @@ subroutine ridge_ice (dt, ndtd, & faero_ocn(it) = faero_ocn(it) + maero(it)*dti enddo endif - if (tr_iso) then - ! check size fiso_ocn vs n_iso ??? - do it = 1, n_iso - fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti - enddo + if (present(fiso_ocn)) then + if (tr_iso) then + ! check size fiso_ocn vs n_iso ??? + do it = 1, n_iso + fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti + enddo + endif endif if (present(fpond)) then fpond = fpond - mpond ! units change later @@ -1826,12 +1824,6 @@ subroutine icepack_step_ridge (dt, ndtd, & real (kind=dbl_kind) :: & dtt ! thermo time step - real (kind=dbl_kind) :: & - l_closing ! local rate of closing due to divergence/shear (1/s) - - logical (kind=log_kind) :: & - l_closing_flag ! flag if closing is passed - logical (kind=log_kind), save :: & first_call = .true. ! first call flag @@ -1859,14 +1851,6 @@ subroutine icepack_step_ridge (dt, ndtd, & ! it may be out of whack, which the ridging helps fix).-ECH !----------------------------------------------------------------- - if (present(closing)) then - l_closing_flag = .true. - l_closing = closing - else - l_closing_flag = .false. - l_closing = c0 - endif - call ridge_ice (dt, ndtd, & ncat, n_aero, & nilyr, nslyr, & @@ -1892,8 +1876,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dardg1ndt, dardg2ndt, & dvirdgndt, & araftn, vraftn, & - l_closing_flag, & - l_closing ) + closing ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- diff --git a/columnphysics/icepack_orbital.F90 b/columnphysics/icepack_orbital.F90 index 4c7c53ccf..13f9d4824 100644 --- a/columnphysics/icepack_orbital.F90 +++ b/columnphysics/icepack_orbital.F90 @@ -177,11 +177,24 @@ subroutine compute_coszen (tlat, tlon, & real (kind=dbl_kind) :: ydayp1 ! day of year plus one time step + logical (kind=log_kind), save :: & + first_call = .true. ! first call flag + character(len=*),parameter :: subname='(compute_coszen)' ! Solar declination for next time step #ifdef CESMCOUPLED + if (icepack_chkoptargflag(first_call)) then + if (.not.(present(days_per_year) .and. & + present(nextsw_cday) .and. & + present(calendar_type))) then + call icepack_warnings_add(subname//' error in CESMCOUPLED args') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif + if (calendar_type == "GREGORIAN") then ydayp1 = min(nextsw_cday, real(days_per_year,kind=dbl_kind)) else @@ -206,6 +219,8 @@ subroutine compute_coszen (tlat, tlon, & endif #endif + first_call = .false. + end subroutine compute_coszen !=============================================================================== diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 1063c4b29..d30214345 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -66,10 +66,7 @@ module icepack_shortwave implicit none private - public :: run_dEdd, & - shortwave_ccsm3, & - compute_shortwave_trcr, & - icepack_prep_radiation, & + public :: icepack_prep_radiation, & icepack_step_radiation real (kind=dbl_kind), parameter :: & @@ -178,11 +175,11 @@ subroutine shortwave_ccsm3 (aicen, vicen, & alvdfns, & ! visible, diffuse, snow (fraction) alidfns ! near-ir, diffuse, snow (fraction) - real (kind=dbl_kind), dimension(:), allocatable :: & - l_fswthru_vdr , & ! vis dir SW through ice to ocean (W m-2) - l_fswthru_vdf , & ! vis dif SW through ice to ocean (W m-2) - l_fswthru_idr , & ! nir dir SW through ice to ocean (W m-2) - l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) + real (kind=dbl_kind) :: & + l_fswthru_vdr, & ! vis dir SW through ice to ocean (W m-2) + l_fswthru_vdf, & ! vis dif SW through ice to ocean (W m-2) + l_fswthru_idr, & ! nir dir SW through ice to ocean (W m-2) + l_fswthru_idf ! nir dif SW through ice to ocean (W m-2) character(len=*),parameter :: subname='(shortwave_ccsm3)' @@ -190,11 +187,6 @@ subroutine shortwave_ccsm3 (aicen, vicen, & ! Solar radiation: albedo and absorbed shortwave !----------------------------------------------------------------- - allocate(l_fswthru_vdr(ncat)) - allocate(l_fswthru_vdf(ncat)) - allocate(l_fswthru_idr(ncat)) - allocate(l_fswthru_idf(ncat)) - ! For basic shortwave, set coszen to a constant between 0 and 1. coszen = p5 ! sun above the horizon @@ -295,29 +287,24 @@ subroutine shortwave_ccsm3 (aicen, vicen, & fswsfc=fswsfc(n), & fswint=fswint(n), & fswthru=fswthru(n), & - fswthru_vdr=l_fswthru_vdr(n),& - fswthru_vdf=l_fswthru_vdf(n),& - fswthru_idr=l_fswthru_idr(n),& - fswthru_idf=l_fswthru_idf(n),& + fswthru_vdr=l_fswthru_vdr,& + fswthru_vdf=l_fswthru_vdf,& + fswthru_idr=l_fswthru_idr,& + fswthru_idf=l_fswthru_idf,& fswpenl=fswpenl(:,n), & Iswabs=Iswabs(:,n)) if (icepack_warnings_aborted(subname)) return + if(present(fswthru_vdr)) fswthru_vdr(n) = l_fswthru_vdr + if(present(fswthru_vdf)) fswthru_vdf(n) = l_fswthru_vdf + if(present(fswthru_idr)) fswthru_idr(n) = l_fswthru_idr + if(present(fswthru_idf)) fswthru_idf(n) = l_fswthru_idf + endif ! aicen > puny enddo ! ncat - if(present(fswthru_vdr)) fswthru_vdr = l_fswthru_vdr - if(present(fswthru_vdf)) fswthru_vdf = l_fswthru_vdf - if(present(fswthru_idr)) fswthru_idr = l_fswthru_idr - if(present(fswthru_idf)) fswthru_idf = l_fswthru_idf - - deallocate(l_fswthru_vdr) - deallocate(l_fswthru_vdf) - deallocate(l_fswthru_idr) - deallocate(l_fswthru_idf) - end subroutine shortwave_ccsm3 !======================================================================= @@ -879,8 +866,10 @@ subroutine run_dEdd(dt, ncat, & fswthrun_idr, & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real(kind=dbl_kind), dimension(:,:), intent(inout), optional :: & + rsnow ! snow grain radius tracer (10^-6 m) + real(kind=dbl_kind), dimension(:,:), intent(inout) :: & - rsnow , & ! snow grain radius tracer (10^-6 m) Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) Iswabsn , & ! SW radiation absorbed in ice layers (W m-2) fswpenln ! visible SW entering ice layers (W m-2) @@ -903,6 +892,7 @@ subroutine run_dEdd(dt, ncat, & alvl ! area fraction of level ice real (kind=dbl_kind), dimension (nslyr) :: & + l_rsnow , & ! local array for snow grain radius tracer (10^-6 m) rhosnwn , & ! snow density (kg/m3) rsnwn ! snow grain radius (micrometers) @@ -929,7 +919,7 @@ subroutine run_dEdd(dt, ncat, & logical (kind=log_kind) :: & linitonly ! local initonly value - real (kind=dbl_kind), dimension(:), allocatable :: & + real (kind=dbl_kind) :: & l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W m-2) l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W m-2) l_fswthrun_idr , & ! nir dir SW through ice to ocean (W m-2) @@ -937,16 +927,13 @@ subroutine run_dEdd(dt, ncat, & character(len=*),parameter :: subname='(run_dEdd)' - allocate(l_fswthrun_vdr(ncat)) - allocate(l_fswthrun_vdf(ncat)) - allocate(l_fswthrun_idr(ncat)) - allocate(l_fswthrun_idf(ncat)) - linitonly = .false. if (present(initonly)) then linitonly = initonly endif + l_rsnow = c0 + ! cosine of the zenith angle #ifdef CESMCOUPLED call compute_coszen (tlat, tlon, & @@ -973,13 +960,17 @@ subroutine run_dEdd(dt, ncat, & if (aicen(n) > puny) then + if (present(rsnow)) then + l_rsnow(:) = rsnow(:,n) + endif + call shortwave_dEdd_set_snow(nslyr, R_snw, & dT_mlt, rsnw_mlt, & aicen(n), vsnon(n), & Tsfcn(n), fsn, & hs0, hsn, & rhosnwn, rsnwn, & - rsnow(:,n)) + l_rsnow(:)) if (icepack_warnings_aborted(subname)) return ! set pond properties @@ -1002,7 +993,7 @@ subroutine run_dEdd(dt, ncat, & Tsfcn(n), fsn, & hs0, hsnlvl, & rhosnwn(:), rsnwn(:), & - rsnow(:,n)) + l_rsnow(:)) if (icepack_warnings_aborted(subname)) return endif ! snwredist @@ -1122,10 +1113,10 @@ subroutine run_dEdd(dt, ncat, & alidrn(n), alidfn(n), & fswsfcn(n), fswintn(n), & fswthru=fswthrun(n), & - fswthru_vdr=l_fswthrun_vdr(n), & - fswthru_vdf=l_fswthrun_vdf(n), & - fswthru_idr=l_fswthrun_idr(n), & - fswthru_idf=l_fswthrun_idf(n), & + fswthru_vdr=l_fswthrun_vdr, & + fswthru_vdf=l_fswthrun_vdf, & + fswthru_idr=l_fswthrun_idr, & + fswthru_idf=l_fswthrun_idf, & Sswabs=Sswabsn(:,n), & Iswabs=Iswabsn(:,n), & albice=albicen(n), & @@ -1137,26 +1128,23 @@ subroutine run_dEdd(dt, ncat, & if (icepack_warnings_aborted(subname)) return - if (.not. snwgrain) then - do k = 1,nslyr - rsnow(k,n) = rsnwn(k) ! for history - enddo + if(present(fswthrun_vdr)) fswthrun_vdr(n) = l_fswthrun_vdr + if(present(fswthrun_vdf)) fswthrun_vdf(n) = l_fswthrun_vdf + if(present(fswthrun_idr)) fswthrun_idr(n) = l_fswthrun_idr + if(present(fswthrun_idf)) fswthrun_idf(n) = l_fswthrun_idf + + if (present(rsnow)) then + if (.not. snwgrain) then + do k = 1,nslyr + rsnow(k,n) = rsnwn(k) ! set rsnow for history + enddo + endif endif endif ! aicen > puny enddo ! ncat - if(present(fswthrun_vdr)) fswthrun_vdr = l_fswthrun_vdr - if(present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf - if(present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr - if(present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf - - deallocate(l_fswthrun_vdr) - deallocate(l_fswthrun_vdf) - deallocate(l_fswthrun_idr) - deallocate(l_fswthrun_idf) - end subroutine run_dEdd !======================================================================= @@ -4106,39 +4094,14 @@ subroutine icepack_step_radiation (dt, ncat, & integer (kind=int_kind) :: & n ! thickness category index - logical (kind=log_kind) :: & - linitonly ! local flag for initonly - real(kind=dbl_kind) :: & hin, & ! Ice thickness (m) hbri ! brine thickness (m) - real (kind=dbl_kind), dimension(:), allocatable :: & - l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) - l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) - l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) - l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) - - real (kind=dbl_kind), dimension(:,:), allocatable :: & - l_rsnow ! snow grain radius tracer (10^-6 m) - character(len=*),parameter :: subname='(icepack_step_radiation)' - allocate(l_fswthrun_vdr(ncat)) - allocate(l_fswthrun_vdf(ncat)) - allocate(l_fswthrun_idr(ncat)) - allocate(l_fswthrun_idf(ncat)) - hin = c0 hbri = c0 - linitonly = .false. - if (present(initonly)) then - linitonly = initonly - endif - - allocate(l_rsnow (nslyr,ncat)) - l_rsnow = c0 - if (present(rsnow)) l_rsnow = rsnow ! Initialize do n = 1, ncat @@ -4175,7 +4138,7 @@ subroutine icepack_step_radiation (dt, ncat, & enddo endif - if (calc_Tsfc) then + if (calc_Tsfc) then if (trim(shortwave) == 'dEdd') then ! delta Eddington call run_dEdd(dt, ncat, & @@ -4209,10 +4172,10 @@ subroutine icepack_step_radiation (dt, ncat, & alidrn, alidfn, & fswsfcn, fswintn, & fswthrun=fswthrun, & - fswthrun_vdr=l_fswthrun_vdr, & - fswthrun_vdf=l_fswthrun_vdf, & - fswthrun_idr=l_fswthrun_idr, & - fswthrun_idf=l_fswthrun_idf, & + fswthrun_vdr=fswthrun_vdr, & + fswthrun_vdf=fswthrun_vdf, & + fswthrun_idr=fswthrun_idr, & + fswthrun_idf=fswthrun_idf, & fswpenln=fswpenln, & Sswabsn=Sswabsn, & Iswabsn=Iswabsn, & @@ -4223,9 +4186,9 @@ subroutine icepack_step_radiation (dt, ncat, & snowfracn=snowfracn, & dhsn=dhsn, & ffracn=ffracn, & - rsnow=l_rsnow, & + rsnow=rsnow, & l_print_point=l_print_point, & - initonly=linitonly) + initonly=initonly) if (icepack_warnings_aborted(subname)) return elseif (trim(shortwave) == 'ccsm3') then @@ -4243,10 +4206,10 @@ subroutine icepack_step_radiation (dt, ncat, & alvdfn, alidfn, & fswsfcn, fswintn, & fswthru=fswthrun, & - fswthru_vdr=l_fswthrun_vdr,& - fswthru_vdf=l_fswthrun_vdf,& - fswthru_idr=l_fswthrun_idr,& - fswthru_idf=l_fswthrun_idf,& + fswthru_vdr=fswthrun_vdr,& + fswthru_vdf=fswthrun_vdf,& + fswthru_idr=fswthrun_idr,& + fswthru_idf=fswthrun_idf,& fswpenl=fswpenln, & Iswabs=Iswabsn, & Sswabs=Sswabsn, & @@ -4300,17 +4263,6 @@ subroutine icepack_step_radiation (dt, ncat, & endif ! calc_Tsfc - if (present(fswthrun_vdr)) fswthrun_vdr = l_fswthrun_vdr - if (present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf - if (present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr - if (present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf - - deallocate(l_fswthrun_vdr) - deallocate(l_fswthrun_vdf) - deallocate(l_fswthrun_idr) - deallocate(l_fswthrun_idf) - deallocate(l_rsnow) - end subroutine icepack_step_radiation !======================================================================= diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 7aa553c95..7a090ef91 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -28,7 +28,7 @@ module icepack_therm_bl99 implicit none private - public :: surface_fluxes, temperature_changes + public :: temperature_changes real (kind=dbl_kind), parameter :: & betak = 0.13_dbl_kind, & ! constant in formula for k (W m-1 ppt-1) @@ -79,8 +79,7 @@ subroutine temperature_changes (dt, & real (kind=dbl_kind), intent(in) :: & dt ! time step - real (kind=dbl_kind), & - intent(in) :: & + real (kind=dbl_kind), intent(in) :: & rhoa , & ! air density (kg/m^3) flw , & ! incoming longwave radiation (W/m^2) potT , & ! air potential temperature (K) @@ -89,8 +88,7 @@ subroutine temperature_changes (dt, & lhcoef , & ! transfer coefficient for latent heat Tbot ! ice bottom surface temperature (deg C) - real (kind=dbl_kind), & - intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) fswint ! SW absorbed in ice interior below surface (W m-2) @@ -99,12 +97,10 @@ subroutine temperature_changes (dt, & hslyr , & ! snow layer thickness (m) einit ! initial energy of melting (J m-2) - real (kind=dbl_kind), dimension (nslyr), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nslyr), intent(inout) :: & Sswabs ! SW radiation absorbed in snow layers (W m-2) - real (kind=dbl_kind), dimension (nilyr), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nilyr), intent(inout) :: & Iswabs ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), intent(inout):: & @@ -117,21 +113,17 @@ subroutine temperature_changes (dt, & real (kind=dbl_kind), intent(out):: & fcondbot ! downward cond flux at bottom surface (W m-2) - real (kind=dbl_kind), & - intent(inout):: & + real (kind=dbl_kind), intent(inout):: & Tsf ! ice/snow surface temperature, Tsfcn - real (kind=dbl_kind), dimension (nilyr), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nilyr), intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) zTin ! internal ice layer temperatures - real (kind=dbl_kind), dimension (nilyr), & - intent(in) :: & + real (kind=dbl_kind), dimension (nilyr), intent(in) :: & zSin ! internal ice layer salinities - real (kind=dbl_kind), dimension (nslyr), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nslyr), intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zTsn ! internal snow layer temperatures @@ -833,8 +825,7 @@ subroutine conductivity (l_snow, & zTin , & ! internal ice layer temperatures zSin ! internal ice layer salinities - real (kind=dbl_kind), dimension (nilyr+nslyr+1), & - intent(out) :: & + real (kind=dbl_kind), dimension (nilyr+nslyr+1), intent(out) :: & kh ! effective conductivity at interfaces (W m-2 deg-1) ! local variables @@ -937,21 +928,18 @@ subroutine surface_fluxes (Tsf, fswsfc, & shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat - real (kind=dbl_kind), & - intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & fsensn , & ! surface downward sensible heat (W m-2) flatn , & ! surface downward latent heat (W m-2) flwoutn , & ! upward LW at surface (W m-2) fsurfn ! net flux to top surface, excluding fcondtopn - real (kind=dbl_kind), & - intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) dflwout_dT ! deriv of flwout wrt Tsf (W m-2 deg-1) - real (kind=dbl_kind), & - intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & dfsurf_dT ! derivative of fsurfn wrt Tsf character(len=*),parameter :: subname='(surface_fluxes)' @@ -1001,8 +989,7 @@ subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, & nilyr , & ! number of ice layers nslyr ! number of snow layers - logical (kind=log_kind), & - intent(in) :: & + logical (kind=log_kind), intent(in) :: & l_snow , & ! true if snow temperatures are computed l_cold ! true if surface temperature is computed @@ -1025,12 +1012,10 @@ subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, & Tsn_init ! snow temp at beginning of time step ! Note: no absorbed SW in snow layers - real (kind=dbl_kind), dimension (nslyr+nilyr+1), & - intent(in) :: & + real (kind=dbl_kind), dimension (nslyr+nilyr+1), intent(in) :: & kh ! effective conductivity at layer interfaces - real (kind=dbl_kind), dimension (nslyr+nilyr+1), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nslyr+nilyr+1), intent(inout) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements @@ -1250,8 +1235,7 @@ subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, & nilyr , & ! number of ice layers nslyr ! number of snow layers - logical (kind=log_kind), & - intent(in) :: & + logical (kind=log_kind), intent(in) :: & l_snow ! true if snow temperatures are computed real (kind=dbl_kind), intent(in) :: & @@ -1266,19 +1250,16 @@ subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, & Tsn_init ! snow temp at beginning of time step ! Note: no absorbed SW in snow layers - real (kind=dbl_kind), dimension (nslyr+nilyr+1), & - intent(in) :: & + real (kind=dbl_kind), dimension (nslyr+nilyr+1), intent(in) :: & kh ! effective conductivity at layer interfaces - real (kind=dbl_kind), dimension (nslyr+nilyr+1), & - intent(inout) :: & + real (kind=dbl_kind), dimension (nslyr+nilyr+1), intent(inout) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements rhs ! rhs of tri-diagonal matrix eqn. - real (kind=dbl_kind), intent(in), & - optional :: & + real (kind=dbl_kind), intent(in) :: & fcondtopn ! conductive flux at top sfc, positive down (W/m^2) ! local variables diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index e4af94abc..93ed8bebd 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -56,10 +56,7 @@ module icepack_therm_itd implicit none private - public :: linear_itd, & - add_new_ice, & - lateral_melt, & - icepack_step_therm2 + public :: icepack_step_therm2 !======================================================================= diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 565e82a59..ac7cf232a 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -61,9 +61,7 @@ module icepack_therm_vertical implicit none private - public :: frzmlt_bottom_lateral, & - thermo_vertical, & - icepack_step_therm1 + public :: icepack_step_therm1 !======================================================================= @@ -2268,15 +2266,19 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2) dsnow , & ! change in snow depth (m/step-->cm/day) - meltsliq , & ! mass of snow melt (kg/m^2) fsloss ! rate of snow loss to leads (kg/m^2/s) + real (kind=dbl_kind), intent(out), optional :: & + meltsliq ! mass of snow melt (kg/m^2) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & Qa_iso , & ! isotope specific humidity (kg/kg) Qref_iso , & ! isotope 2m atm ref spec humidity (kg/kg) fiso_atm , & ! isotope deposition rate (kg/m^2 s) fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) - fiso_evap , & ! isotope evaporation (kg/m^2/s) + fiso_evap ! isotope evaporation (kg/m^2/s) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & meltsliqn ! mass of snow melt (kg/m^2) real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: & @@ -2396,7 +2398,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & smice , & ! tracer for mass of ice in snow (kg/m^3) smliq ! tracer for mass of liquid in snow (kg/m^3) - real (kind=dbl_kind), allocatable, dimension(:) :: & + real (kind=dbl_kind), dimension(ncat) :: & l_meltsliqn ! mass of snow melt local (kg/m^2) real (kind=dbl_kind) :: & @@ -2457,11 +2459,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & smice(:) = c0 smliq(:) = c0 - allocate(l_meltsliqn(ncat)) - l_meltsliqn = c0 - if (present(meltsliqn)) l_meltsliqn = meltsliqn l_meltsliq = c0 - if (present(meltsliq )) l_meltsliq = meltsliq + l_meltsliqn = c0 !----------------------------------------------------------------- ! Initialize rate of snow loss to leads @@ -2555,7 +2554,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & congeln(n) = c0 snoicen(n) = c0 dsnown (n) = c0 - l_meltsliqn(n) = c0 Trefn = c0 Qrefn = c0 @@ -2880,14 +2878,14 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & meltb=meltb, snoicen=snoicen(n),& dsnow=dsnow, dsnown=dsnown(n), & congel=congel, snoice=snoice, & - meltsliq=l_meltsliq, & - meltsliqn=l_meltsliqn(n), & - Uref=Uref, Urefn=Urefn, & - Qref_iso=Qref_iso, & - Qrefn_iso=Qrefn_iso, & - fiso_ocn=fiso_ocn, & - fiso_ocnn=fiso_ocnn, & - fiso_evap=fiso_evap, & + meltsliq=l_meltsliq, & + meltsliqn=l_meltsliqn(n), & + Uref=Uref, Urefn=Urefn, & + Qref_iso=Qref_iso, & + Qrefn_iso=Qrefn_iso, & + fiso_ocn=fiso_ocn, & + fiso_ocnn=fiso_ocnn, & + fiso_evap=fiso_evap, & fiso_evapn=fiso_evapn) if (icepack_warnings_aborted(subname)) return @@ -2923,7 +2921,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (present(meltsliqn )) meltsliqn = l_meltsliqn if (present(meltsliq )) meltsliq = l_meltsliq - deallocate(l_meltsliqn) !----------------------------------------------------------------- ! Calculate ponds from the topographic scheme diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 9f2f2f77c..c8135dea1 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -907,7 +907,7 @@ end subroutine icepack_query_tracer_indices subroutine icepack_write_tracer_indices(iounit) - integer, intent(in), optional :: iounit + integer, intent(in) :: iounit !autodocument_end diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index ec77a9707..26df28a29 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -102,7 +102,8 @@ subroutine icedrv_initialize floe_rad_l=floe_rad_l, & ! fsd size lower bound in m (radius) floe_rad_c=floe_rad_c, & ! fsd size bin centre in m (radius) floe_binwidth=floe_binwidth, & ! fsd size bin width in m (radius) - c_fsd_range=c_fsd_range) ! string for history output + c_fsd_range=c_fsd_range , & ! string for history output + write_diags=.true.) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted(subname)) then call icedrv_system_abort(file=__FILE__,line=__LINE__) diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include index 5919e7742..95db22ba6 100644 --- a/doc/source/user_guide/interfaces.include +++ b/doc/source/user_guide/interfaces.include @@ -55,18 +55,18 @@ icepack_atm_boundary shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat - real (kind=dbl_kind), intent(in), optional, dimension(:) :: & + real (kind=dbl_kind), intent(in), dimension(:), optional :: & Qa_iso ! specific isotopic humidity (kg/kg) - real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & + real (kind=dbl_kind), intent(inout), dimension(:), optional :: & Qref_iso ! reference specific isotopic humidity (kg/kg) - real (kind=dbl_kind), optional, intent(in) :: & + real (kind=dbl_kind), intent(in), optional :: & uvel , & ! x-direction ice speed (m/s) vvel , & ! y-direction ice speed (m/s) zlvs ! atm level height for scalars (if different than zlvl) (m) - real (kind=dbl_kind), optional, intent(out) :: & + real (kind=dbl_kind), intent(out), optional :: & Uref ! reference height wind speed (m/s) @@ -585,7 +585,7 @@ icepack_step_ridge faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) real (kind=dbl_kind), dimension(:,:), intent(inout) :: & @@ -2029,10 +2029,10 @@ icepack_step_therm2 yday ! day of year ! water isotopes - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) - real (kind=dbl_kind), optional, intent(in) :: & + real (kind=dbl_kind), intent(in), optional :: & HDO_ocn , & ! ocean concentration of HDO (kg/kg) H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) @@ -2276,15 +2276,15 @@ icepack_step_therm1 zlvs) integer (kind=int_kind), intent(in) :: & - ncat , & ! number of thickness categories - nilyr , & ! number of ice layers - nslyr ! number of snow layers + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt , & ! time step - uvel , & ! x-component of velocity (m/s) - vvel , & ! y-component of velocity (m/s) - strax , & ! wind stress components (N/m^2) + uvel , & ! x-component of velocity (m/s) + vvel , & ! y-component of velocity (m/s) + strax , & ! wind stress components (N/m^2) stray , & ! yday ! day of year @@ -2297,42 +2297,42 @@ icepack_step_therm1 real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration - vice , & ! volume per unit area of ice (m) - vsno , & ! volume per unit area of snow (m) + vice , & ! volume per unit area of ice (m) + vsno , & ! volume per unit area of snow (m) zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) - uatm , & ! wind velocity components (m/s) - vatm , & - wind , & ! wind speed (m/s) - potT , & ! air potential temperature (K) - Tair , & ! air temperature (K) - Qa , & ! specific humidity (kg/kg) - rhoa , & ! air density (kg/m^3) - frain , & ! rainfall rate (kg/m^2 s) - fsnow , & ! snowfall rate (kg/m^2 s) - fpond , & ! fresh water flux to ponds (kg/m^2/s) - fresh , & ! fresh water flux to ocean (kg/m^2/s) - fsalt , & ! salt flux to ocean (kg/m^2/s) - fhocn , & ! net heat flux to ocean (W/m^2) - fswthru , & ! shortwave penetrating to ocean (W/m^2) + uatm , & ! wind velocity components (m/s) + vatm , & ! (m/s) + wind , & ! wind speed (m/s) + potT , & ! air potential temperature (K) + Tair , & ! air temperature (K) + Qa , & ! specific humidity (kg/kg) + rhoa , & ! air density (kg/m^3) + frain , & ! rainfall rate (kg/m^2 s) + fsnow , & ! snowfall rate (kg/m^2 s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fsalt , & ! salt flux to ocean (kg/m^2/s) + fhocn , & ! net heat flux to ocean (W/m^2) + fswthru , & ! shortwave penetrating to ocean (W/m^2) 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) - fsens , & ! sensible heat flux (W/m^2) - flat , & ! latent heat flux (W/m^2) + fsens , & ! sensible heat flux (W/m^2) + flat , & ! latent heat flux (W/m^2) fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2) - flw , & ! incoming longwave radiation (W/m^2) - flwout , & ! outgoing longwave radiation (W/m^2) - evap , & ! evaporative water flux (kg/m^2/s) - evaps , & ! evaporative water flux over snow (kg/m^2/s) + flw , & ! incoming longwave radiation (W/m^2) + flwout , & ! outgoing longwave radiation (W/m^2) + 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) congel , & ! basal ice growth (m/step-->cm/day) snoice , & ! snow-ice formation (m/step-->cm/day) - Tref , & ! 2m atm reference temperature (K) - Qref , & ! 2m atm reference spec humidity (kg/kg) - Uref , & ! 10m atm reference wind speed (m/s) + Tref , & ! 2m atm reference temperature (K) + Qref , & ! 2m atm reference spec humidity (kg/kg) + Uref , & ! 10m atm reference wind speed (m/s) Cdn_atm , & ! atm drag coefficient Cdn_ocn , & ! ocn drag coefficient - hfreebd , & ! freeboard (m) + hfreebd , & ! freeboard (m) hdraft , & ! draft of ice + snow column (Stoessel1993) hridge , & ! ridge height distrdg , & ! distance between ridges @@ -2353,14 +2353,14 @@ icepack_step_therm1 strocnxT , & ! ice-ocean stress, x-direction strocnyT , & ! ice-ocean stress, y-direction fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) - frzmlt , & ! freezing/melting potential (W/m^2) + frzmlt , & ! freezing/melting potential (W/m^2) rside , & ! fraction of ice that melts laterally - fside , & ! lateral heat flux (W/m^2) - sst , & ! sea surface temperature (C) - Tf , & ! freezing temperature (C) - Tbot , & ! ice bottom surface temperature (deg C) - Tsnice , & ! snow ice interface temperature (deg C) - sss , & ! sea surface salinity (ppt) + fside , & ! lateral heat flux (W/m^2) + sst , & ! sea surface temperature (C) + Tf , & ! freezing temperature (C) + Tbot , & ! ice bottom surface temperature (deg C) + Tsnice , & ! snow ice interface temperature (deg C) + sss , & ! sea surface salinity (ppt) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) meltb , & ! basal ice melt (m/step-->cm/day) @@ -2373,86 +2373,88 @@ icepack_step_therm1 fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2) dsnow , & ! change in snow depth (m/step-->cm/day) - meltsliq , & ! mass of snow melt (kg/m^2) - fsloss ! rate of snow loss to leads (kg/m^2/s) - - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & - Qa_iso , & ! isotope specific humidity (kg/kg) - Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg) - fiso_atm , & ! isotope deposition rate (kg/m^2 s) - fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) - fiso_evap , & ! isotope evaporation (kg/m^2/s) - meltsliqn ! mass of snow melt (kg/m^2) - - real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: & - rsnwn , & ! snow grain radius (10^-6 m) - smicen , & ! tracer for mass of ice in snow (kg/m^3) - smliqn ! tracer for mass of liquid in snow (kg/m^3) - - real (kind=dbl_kind), optional, intent(in) :: & - HDO_ocn , & ! ocean concentration of HDO (kg/kg) - H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) - H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) + meltsliq , & ! mass of snow melt (kg/m^2) + fsloss ! rate of snow loss to leads (kg/m^2/s) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + Qa_iso , & ! isotope specific humidity (kg/kg) + Qref_iso , & ! isotope 2m atm ref spec humidity (kg/kg) + fiso_atm , & ! isotope deposition rate (kg/m^2 s) + fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) + fiso_evap , & ! isotope evaporation (kg/m^2/s) + meltsliqn ! mass of snow melt (kg/m^2) + + real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: & + rsnwn , & ! snow grain radius (10^-6 m) + smicen , & ! tracer for mass of ice in snow (kg/m^3) + smliqn ! tracer for mass of liq in snow (kg/m^3) + + real (kind=dbl_kind), intent(in), optional :: & + HDO_ocn , & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice - vicen_init , & ! volume per unit area of ice (m) - vsnon_init , & ! volume per unit area of snow (m) + vicen_init , & ! volume per unit area of ice (m) + vsnon_init , & ! volume per unit area of snow (m) aicen , & ! concentration of ice - vicen , & ! volume per unit area of ice (m) - vsnon , & ! volume per unit area of snow (m) + vicen , & ! volume per unit area of ice (m) + vsnon , & ! volume per unit area of snow (m) Tsfc , & ! ice/snow surface temperature, Tsfcn alvl , & ! level ice area fraction vlvl , & ! level ice volume fraction apnd , & ! melt pond area fraction - hpnd , & ! melt pond depth (m) - ipnd , & ! melt pond refrozen lid thickness (m) + hpnd , & ! melt pond depth (m) + ipnd , & ! melt pond refrozen lid thickness (m) iage , & ! volume-weighted ice age FY , & ! area-weighted first-year ice area fsurfn , & ! net flux to top surface, excluding fcondtop - fcondtopn , & ! downward cond flux at top surface (W m-2) + fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) - flatn , & ! latent heat flux (W m-2) - fsensn , & ! sensible heat flux (W m-2) + flatn , & ! latent heat flux (W m-2) + fsensn , & ! sensible heat flux (W m-2) fsurfn_f , & ! net flux to top surface, excluding fcondtop - fcondtopn_f , & ! downward cond flux at top surface (W m-2) - flatn_f , & ! latent heat flux (W m-2) - fsensn_f , & ! sensible heat flux (W m-2) - fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) - fswthrun , & ! SW through ice to ocean (W/m^2) + fcondtopn_f , & ! downward cond flux at top surface (W m-2) + flatn_f , & ! latent heat flux (W m-2) + fsensn_f , & ! sensible heat flux (W m-2) + fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) fswintn , & ! SW absorbed in ice interior, below surface (W m-2) - faero_atm , & ! aerosol deposition rate (kg/m^2 s) - faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) + faero_atm , & ! aerosol deposition rate (kg/m^2 s) + faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond - meltsn , & ! snow melt (m) - melttn , & ! top ice melt (m) - meltbn , & ! bottom ice melt (m) - congeln , & ! congelation ice growth (m) - snoicen , & ! snow-ice growth (m) + meltsn , & ! snow melt (m) + melttn , & ! top ice melt (m) + meltbn , & ! bottom ice melt (m) + congeln , & ! congelation ice growth (m) + snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) - real (kind=dbl_kind), optional, dimension(:), intent(inout) :: & - fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) - fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) - fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) - fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), intent(in) :: & + fswthrun ! SW through ice to ocean (W/m^2) + + real (kind=dbl_kind), dimension(:), intent(in), optional :: & + fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) + fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind), dimension(:,:), intent(inout) :: & - zqsn , & ! snow layer enthalpy (J m-3) - zqin , & ! ice layer enthalpy (J m-3) + zqsn , & ! snow layer enthalpy (J m-3) + zqin , & ! ice layer enthalpy (J m-3) zSin , & ! internal ice layer salinities Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) - Iswabsn ! SW radiation absorbed in ice layers (W m-2) + Iswabsn ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & - aerosno , & ! snow aerosol tracer (kg/m^2) - aeroice ! ice aerosol tracer (kg/m^2) + aerosno , & ! snow aerosol tracer (kg/m^2) + aeroice ! ice aerosol tracer (kg/m^2) - real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: & - isosno , & ! snow isotope tracer (kg/m^2) - isoice ! ice isotope tracer (kg/m^2) + real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: & + isosno , & ! snow isotope tracer (kg/m^2) + isoice ! ice isotope tracer (kg/m^2) icepack_tracers.F90 @@ -2776,7 +2778,7 @@ icepack_write_tracer_indices subroutine icepack_write_tracer_indices(iounit) - integer, intent(in), optional :: iounit + integer, intent(in) :: iounit diff --git a/doc/source/user_guide/lg_sequence.rst b/doc/source/user_guide/lg_sequence.rst index 628f96bff..09127b8f2 100755 --- a/doc/source/user_guide/lg_sequence.rst +++ b/doc/source/user_guide/lg_sequence.rst @@ -42,7 +42,6 @@ icepack_kinds module:: icepack_tracers defines a handful of parameters that provide information about maximum array sizes for static dimensioning:: - use icepack_tracers, only: icepack_max_iso => max_iso use icepack_tracers, only: icepack_max_nbtrcr => max_nbtrcr use icepack_tracers, only: icepack_max_algae => max_algae use icepack_tracers, only: icepack_max_dic => max_dic @@ -50,6 +49,7 @@ about maximum array sizes for static dimensioning:: use icepack_tracers, only: icepack_max_don => max_don use icepack_tracers, only: icepack_max_fe => max_fe use icepack_tracers, only: icepack_max_aero => max_aero + use icepack_tracers, only: icepack_max_iso => max_iso use icepack_tracers, only: icepack_nmodal1 => nmodal1 use icepack_tracers, only: icepack_nmodal2 => nmodal2 use icepack_parameters,only: icepack_nspint => nspint @@ -111,8 +111,14 @@ impact of waves on sea ice:: use icepack_wavefracspec, only: icepack_init_wave use icepack_wavefracspec, only: icepack_step_wavefracture +icepack_snow provides a routine to initialize the snow physics +and a routine to update the snow physics:: + + use icepack_snow, only: icepack_init_snow + use icepack_snow, only: icepack_step_snow + icepack_shortwave provides a routine to initialize the radiation -computation and an routine to update the radiation computation:: +computation and a routine to update the radiation computation:: use icepack_shortwave, only: icepack_prep_radiation use icepack_shortwave, only: icepack_step_radiation From 37e215b5329463591d2cce228883fd34aa0ea3be Mon Sep 17 00:00:00 2001 From: davidclemenssewall <51925832+davidclemenssewall@users.noreply.github.com> Date: Fri, 3 Mar 2023 10:13:20 -0700 Subject: [PATCH 2/4] driver: icedrv_restart: Added option for restart files to be read/written in NetCDF format (#427) * Trial fix for netCDF start time problem * driver: icedrv_history: fix start time for NetCDF history output The start time for netCDF output was hardcoded to be at 00:00:00 and there was a bug in icedrv_calendar such that idate0 was fixed at Jan. 1. Both bugs are fixed. NOTE: if dt and istep0 are set such that the start time contains fractional seconds, they will be rounded down in the NetCDF output to the nearest integer. * driver: icedrv_history: fix start time for NetCDF History output Start of netCDF time axis is fixed at Jan. 01 00:00:00 of init_year and time value calculation is now fixed. * Added restart_format variable to icedrv_restart_shared and to namelist. No functionality implemented yet, just checking it's added correctly. * need to also add new namelist to icepack_in * Added error handling for incorrect restart format, should probably update final_restart at some point * NetCDF restart file created, next need to populate fields * NetCDF Restart output is now populated * if block error handling in restartfile * driver: icedrv_restart: enabled restart/ice_ic files to be written/read as NetCDF Added the option for restart (aka initial condition) files to be in NetCDF instead of binary for ease of setting custom inital conditions * modified set_nml.ionetcdf to default to netcdf restart files if we set ionetcdf * Fixed error in which some tracers were not being writting. NetCDF restart files are now identical * fixed istep1 = istep0 but and added test_suite for netcdf * removed bgc from the netcdf restart test suite * put netcdf functions into conditional compilation blocks * Removed unused variables (sec0, h0, m0, s0) from icedrv_history * Removing unused variables in icedrv_history (after merging in from consortium-main) * Changed history_cdf namelist to history_format and updated documentation * Fixed warning message if restart_format is incorrect * halfway done reorganizing restartfile * finished reorganization in restartfile, need to do dumpfile next * Halfway done editing dumpfile * finished editing dumpfile * updated namelist settings and documentation. all tests pass * removed the netcdf_nobgc test suite * Fix minor issues found during testing - icedrv_init write format for history_format - dims allocation in icedrv_restart - cleaned up some intent lines (old) - add return at end of io_suite.ts - update namelist documentation of history_format, restart_format Update to current main trunk * Update history_format default Check history_format and restart_format input values --------- Co-authored-by: apcraig --- .gitignore | 6 + configuration/driver/icedrv_RunMod.F90 | 8 +- configuration/driver/icedrv_history.F90 | 9 +- configuration/driver/icedrv_init.F90 | 28 +- configuration/driver/icedrv_restart.F90 | 944 +++++++++++++++--- .../driver/icedrv_restart_shared.F90 | 3 + configuration/scripts/icepack_in | 3 +- configuration/scripts/options/set_nml.histcdf | 1 + .../scripts/options/set_nml.ionetcdf | 1 - configuration/scripts/options/set_nml.restcdf | 1 + configuration/scripts/tests/io_suite.ts | 7 + doc/source/icepack_index.rst | 3 +- doc/source/user_guide/ug_case_settings.rst | 5 +- doc/source/user_guide/ug_implementation.rst | 11 +- 14 files changed, 853 insertions(+), 177 deletions(-) create mode 100644 configuration/scripts/options/set_nml.histcdf delete mode 100644 configuration/scripts/options/set_nml.ionetcdf create mode 100644 configuration/scripts/options/set_nml.restcdf diff --git a/.gitignore b/.gitignore index f32b5b0ea..b3a5c79ae 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,9 @@ doc/build # Ignore testsuite file/directories testsuite* caselist* + +# Ignore compiled .mod files +*.mod + +# Ignore test case directories (no consistent name so we'll just ignore conda) +conda_* \ No newline at end of file diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index 4b7170663..ff963673f 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -35,7 +35,7 @@ subroutine icedrv_run use icedrv_forcing, only: get_forcing, get_wave_spec 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_cdf, history_close + use icedrv_history, only: history_format, history_close logical (kind=log_kind) :: skl_bgc, z_tracers, tr_aero, tr_zaero, & wave_spec, tr_fsd, tr_iso @@ -63,7 +63,7 @@ subroutine icedrv_run call calendar(time) ! at the end of the timestep if (stop_now >= 1) then - if (history_cdf) call history_close() + if (history_format == 'nc') call history_close() exit timeLoop endif @@ -101,7 +101,7 @@ subroutine ice_step use icedrv_diagnostics_bgc, only: hbrine_diags, zsal_diags, bgc_diags use icedrv_flux, only: init_history_therm, init_history_bgc, & daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd, init_history_dyn - use icedrv_history, only: history_cdf, history_write + use icedrv_history, only: history_format, history_write use icedrv_restart, only: dumpfile, final_restart use icedrv_restart_bgc, only: write_restart_bgc use icedrv_step, only: prep_radiation, step_therm1, step_therm2, & @@ -225,7 +225,7 @@ subroutine ice_step if (tr_brine) call hbrine_diags endif - if (history_cdf) then + if (history_format == 'nc') then call history_write() endif diff --git a/configuration/driver/icedrv_history.F90 b/configuration/driver/icedrv_history.F90 index 1f284d4c7..049566b7b 100644 --- a/configuration/driver/icedrv_history.F90 +++ b/configuration/driver/icedrv_history.F90 @@ -22,7 +22,8 @@ module icedrv_history ! history output file info - logical (kind=log_kind), public :: history_cdf ! flag to turn on cdf history files + character (len=char_len), public :: & + history_format ! format of history files, only supported type is 'nc' character (len=char_len_long) :: hist_file ! hist file name @@ -68,11 +69,7 @@ subroutine history_write() count1(1), count2(2), count3(3), count4(4), & ! cdf start/count arrays varid, & ! cdf varid status, & ! cdf status flag - iflag, & ! history file attributes - sec0, & ! number of seconds into the day at istep0 - h0, & ! start hour - m0, & ! start minute - s0 ! start second + iflag ! history file attributes character (len=8) :: & cdate ! date string diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 9918bcbe9..8e4958447 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -62,8 +62,8 @@ subroutine input_data use icedrv_calendar, only: year_init, istep0 use icedrv_calendar, only: dumpfreq, diagfreq, dump_last use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years - use icedrv_history, only: history_cdf - use icedrv_restart_shared, only: restart, restart_dir, restart_file + use icedrv_history, only: history_format + use icedrv_restart_shared, only: restart, restart_dir, restart_file, restart_format use icedrv_flux, only: update_ocn_f, l_mpond_fresh, cpl_bgc use icedrv_flux, only: default_season use icedrv_forcing, only: precip_units, fyear_init, ycycle @@ -133,8 +133,9 @@ subroutine input_data days_per_year, use_leap_years, year_init, istep0, & dt, npt, ndtd, dump_last, & ice_ic, restart, restart_dir, restart_file, & + restart_format, & dumpfreq, diagfreq, diag_file, cpl_bgc, & - conserv_check, history_cdf + conserv_check, history_format namelist /grid_nml/ & kcatbound @@ -254,9 +255,10 @@ subroutine input_data restart = .false. ! if true, read restart files for initialization restart_dir = './' ! write to executable dir for default restart_file = 'iced' ! restart file name prefix - history_cdf = .false. ! history netcdf file flag - ice_ic = 'default' ! initial conditions are specified in the code - ! otherwise, the filename for reading restarts + restart_format = 'bin' ! default restart format is binary, other option 'nc' + 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 ndtd = 1 ! dynamic time steps per thermodynamic time step l_mpond_fresh = .false. ! logical switch for including meltpond freshwater ! flux feedback to ocean model @@ -544,6 +546,16 @@ subroutine input_data call icedrv_system_abort(file=__FILE__,line=__LINE__) endif + if (restart_format /= 'bin' .and. restart_format /= 'nc') then + write (nu_diag,*) 'WARNING: restart_format value unknown '//trim(restart_format) + call icedrv_system_abort(file=__FILE__,line=__LINE__) + endif + + if (history_format /= 'none' .and. history_format /= 'nc') then + write (nu_diag,*) 'WARNING: history_format value unknown '//trim(history_format) + call icedrv_system_abort(file=__FILE__,line=__LINE__) + endif + if (tr_aero .and. trim(shortwave) /= 'dEdd') then write (nu_diag,*) 'WARNING: aerosols activated but dEdd' write (nu_diag,*) 'WARNING: shortwave is not.' @@ -645,7 +657,9 @@ subroutine input_data trim(restart_dir) write(nu_diag,*) ' restart_file = ', & trim(restart_file) - write(nu_diag,1010) ' history_cdf = ', history_cdf + write(nu_diag,1030) ' restart_format = ', & + trim(restart_format) + write(nu_diag,1030) ' history_format = ', trim(history_format) write(nu_diag,*) ' ice_ic = ', & trim(ice_ic) write(nu_diag,1010) ' conserv_check = ', conserv_check diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90 index 67bd4a033..745ec9210 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -10,10 +10,14 @@ module icedrv_restart use icedrv_constants, only: nu_diag, nu_restart, nu_dump use icedrv_constants, only: c0, c1, p5 use icedrv_restart_shared, only: restart, restart_dir, restart_file, lenstr + use icedrv_restart_shared, only: restart_format use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_tracer_flags, icepack_query_tracer_indices use icepack_intfc, only: icepack_query_parameters use icedrv_system, only: icedrv_system_abort +#ifdef USE_NETCDF + use netcdf +#endif implicit none private :: & @@ -26,6 +30,22 @@ module icedrv_restart write_restart_fsd, read_restart_fsd, & write_restart_iso, read_restart_iso, & write_restart_aero, read_restart_aero + + character (len=3), private :: nchar ! + + logical (kind=log_kind), private :: diag ! netCDF diagnostics flag + +#ifdef USE_NETCDF + private :: & + define_rest_field, write_restart_field_nc2D, & + write_restart_field_nc1D, ice_write_nc2D, & + ice_write_nc1D, read_restart_field_net2D, & + read_restart_field_net1D, ice_read_nc2D, & + ice_read_nc1D + + integer (kind=int_kind), private :: & + ncid ! ID for NetCDF file +#endif public :: dumpfile, restartfile, final_restart, & write_restart_field, read_restart_field @@ -37,21 +57,26 @@ module icedrv_restart !======================================================================= !======================================================================= -!---! these subroutines write/read Fortran unformatted data files .. +!---! these subroutines write/read restart data files .. +! Which can either be in unformatted Fortran format +! (restart_format = 'bin') or NetCDF (restart_format = 'nc') !======================================================================= ! Dumps all values needed for a restart -! author Elizabeth C. Hunke, LANL +! authors Elizabeth C. Hunke, LANL +! David Clemens-Sewall, NCAR subroutine dumpfile use icedrv_calendar, only: sec, month, mday, nyr, istep1 use icedrv_calendar, only: time, time_forc, year_init - use icedrv_domain_size, only: nilyr, nslyr, ncat + use icedrv_domain_size, only: nilyr, nslyr, ncat, n_aero, nfsd, nx, n_iso use icedrv_forcing, only: oceanmixed_ice use icedrv_flux, only: scale_factor, swvdr, swvdf, swidr, swidf - use icedrv_flux, only: sst, frzmlt + use icedrv_flux, only: sst, frzmlt, frz_onset, fsnow use icedrv_state, only: aicen, vicen, vsnon, trcrn + use icedrv_arrays_column, only: dhsn, ffracn + use icedrv_arrays_column, only: first_ice, first_ice_real ! local variables @@ -67,32 +92,83 @@ subroutine dumpfile tr_pond_topo, tr_pond_lvl, tr_snow, tr_fsd ! solve_zsal, skl_bgc, z_tracers + integer (kind=int_kind) :: dims(2) + character(len=char_len_long) :: filename character(len=*), parameter :: subname='(dumpfile)' + +#ifdef USE_NETCDF + ! local variables if we're writing in NetCDF + integer (kind=int_kind) :: & + dimid_ni, & ! netCDF identifiers + dimid_ncat, & ! + iflag, & ! netCDF creation flag + status ! status variable from netCDF routine +#endif - ! construct path/file + ! get year iyear = nyr + year_init - 1 - - write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & - restart_dir(1:lenstr(restart_dir)), & - restart_file(1:lenstr(restart_file)),'.', & - iyear,'-',month,'-',mday,'-',sec - + + ! query which tracers are active and their indices call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & - nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) - call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, & - tr_brine_out=tr_brine, & - tr_pond_topo_out=tr_pond_topo, & - tr_pond_lvl_out=tr_pond_lvl,tr_snow_out=tr_snow,tr_fsd_out=tr_fsd) + nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) + + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, & + tr_brine_out=tr_brine, & + tr_pond_topo_out=tr_pond_topo, & + tr_pond_lvl_out=tr_pond_lvl,tr_snow_out=tr_snow,tr_fsd_out=tr_fsd) ! call icepack_query_parameters(solve_zsal_out=solve_zsal, & ! skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & - file=__FILE__,line= __LINE__) + file=__FILE__,line= __LINE__) + + dims(:) = 0 + if (restart_format == 'bin') then + write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & + restart_dir(1:lenstr(restart_dir)), & + restart_file(1:lenstr(restart_file)),'.', & + iyear,'-',month,'-',mday,'-',sec + + open(nu_dump,file=filename,form='unformatted') + write(nu_dump) istep1,time,time_forc + else if (restart_format == 'nc') then +#ifdef USE_NETCDF + ! Change this if you want diagnostic output + diag = .false. + + write(filename,'(a,a,a,i4.4,a,i2.2,a,i2.2,a,i5.5)') & + restart_dir(1:lenstr(restart_dir)), & + restart_file(1:lenstr(restart_file)),'.', & + iyear,'-',month,'-',mday,'-',sec + filename = trim(filename) // '.nc' + + ! Open NetCDF restart file in define mode + iflag = nf90_clobber + status = nf90_create(trim(filename), iflag, ncid) + if (status /= nf90_noerr) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + ! Set time attributes for restart file + status = nf90_put_att(ncid,nf90_global,'istep1',istep1) + status = nf90_put_att(ncid,nf90_global,'time',time) + status = nf90_put_att(ncid,nf90_global,'time_forc',time_forc) + status = nf90_put_att(ncid,nf90_global,'nyr',iyear) + status = nf90_put_att(ncid,nf90_global,'month',month) + status = nf90_put_att(ncid,nf90_global,'mday',mday) + status = nf90_put_att(ncid,nf90_global,'sec',sec) + ! Set dimensions for restart file + status = nf90_def_dim(ncid,'ni',nx,dimid_ni) + status = nf90_def_dim(ncid,'ncat',ncat,dimid_ncat) + + ! Get dimension IDs + dims(1) = dimid_ni + dims(2) = dimid_ncat +#else + call icedrv_system_abort(string=subname//' ERROR: restart_format = "nc" requires USE_NETCDF',file=__FILE__,line=__LINE__) +#endif + endif - open(nu_dump,file=filename,form='unformatted') - write(nu_dump) istep1,time,time_forc write(nu_diag,*) 'Writing ',filename(1:lenstr(filename)) !----------------------------------------------------------------- @@ -101,51 +177,66 @@ subroutine dumpfile ! tracers are written to their own dump/restart binary files. !----------------------------------------------------------------- - call write_restart_field(nu_dump,aicen(:,:),ncat) - call write_restart_field(nu_dump,vicen(:,:),ncat) - call write_restart_field(nu_dump,vsnon(:,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_Tsfc,:),ncat) + call write_restart_field(nu_dump,aicen(:,:),ncat,'aicen',dims) + call write_restart_field(nu_dump,vicen(:,:),ncat,'vicen',dims) + call write_restart_field(nu_dump,vsnon(:,:),ncat,'vsnon',dims) + call write_restart_field(nu_dump,trcrn(:,nt_Tsfc,:),ncat,'Tsfcn',dims) do k=1,nilyr - call write_restart_field(nu_dump,trcrn(:,nt_sice+k-1,:),ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump,trcrn(:,nt_sice+k-1,:),ncat,& + 'sice'//trim(nchar),dims) enddo do k=1,nilyr - call write_restart_field(nu_dump,trcrn(:,nt_qice+k-1,:),ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump,trcrn(:,nt_qice+k-1,:),ncat,& + 'qice'//trim(nchar),dims) enddo do k=1,nslyr - call write_restart_field(nu_dump,trcrn(:,nt_qsno+k-1,:),ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump,trcrn(:,nt_qsno+k-1,:),ncat,& + 'qsno'//trim(nchar),dims) enddo !----------------------------------------------------------------- ! radiation fields !----------------------------------------------------------------- - call write_restart_field(nu_dump,scale_factor,1) - call write_restart_field(nu_dump,swvdr,1) - call write_restart_field(nu_dump,swvdf,1) - call write_restart_field(nu_dump,swidr,1) - call write_restart_field(nu_dump,swidf,1) + call write_restart_field(nu_dump,scale_factor,1,'scale_factor',dims) + call write_restart_field(nu_dump,swvdr,1,'swvdr',dims) + call write_restart_field(nu_dump,swvdf,1,'swvdf',dims) + call write_restart_field(nu_dump,swidr,1,'swidr',dims) + call write_restart_field(nu_dump,swidf,1,'swidf',dims) !----------------------------------------------------------------- ! for mixed layer model !----------------------------------------------------------------- if (oceanmixed_ice) then - call write_restart_field(nu_dump,sst,1) - call write_restart_field(nu_dump,frzmlt,1) + call write_restart_field(nu_dump,sst,1,'sst',dims) + call write_restart_field(nu_dump,frzmlt,1,'frzmlt',dims) endif ! tracers - if (tr_iage) call write_restart_age() ! ice age tracer - if (tr_FY) call write_restart_FY() ! first-year area tracer - if (tr_lvl) call write_restart_lvl() ! level ice tracer - if (tr_pond_lvl) call write_restart_pond_lvl() ! level-ice melt ponds - if (tr_pond_topo) call write_restart_pond_topo() ! topographic melt ponds - if (tr_snow) call write_restart_snow() ! snow metamorphosis tracers - if (tr_iso) call write_restart_iso() ! ice isotopes - if (tr_aero) call write_restart_aero() ! ice aerosols - if (tr_brine) call write_restart_hbrine() ! brine height - if (tr_fsd) call write_restart_fsd() ! floe size distribution -! called in icedrv_RunMod.F90 to prevent circular dependencies -! if (solve_zsal .or. skl_bgc .or. z_tracers) & -! call write_restart_bgc ! biogeochemistry + if (tr_iage) call write_restart_age(dims) ! ice age tracer + if (tr_FY) call write_restart_FY(dims) ! first-year area tracer + if (tr_lvl) call write_restart_lvl(dims) ! level ice tracer + if (tr_pond_lvl) call write_restart_pond_lvl(dims) ! level-ice melt ponds + if (tr_pond_topo) call write_restart_pond_topo(dims) ! topographic melt ponds + if (tr_snow) call write_restart_snow(dims) ! snow metamorphosis tracers + if (tr_iso) call write_restart_iso(dims) ! ice isotopes + if (tr_aero) call write_restart_aero(dims) ! ice aerosols + if (tr_brine) call write_restart_hbrine(dims) ! brine height + if (tr_fsd) call write_restart_fsd(dims) ! floe size distribution + ! called in icedrv_RunMod.F90 to prevent circular dependencies + ! if (solve_zsal .or. skl_bgc .or. z_tracers) & + ! call write_restart_bgc ! biogeochemistry + + ! Close netcdf file + if (restart_format == 'nc') then +#ifdef USE_NETCDF + status = nf90_close(ncid) + if (status /= nf90_noerr) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) +#endif + endif end subroutine dumpfile @@ -158,15 +249,20 @@ subroutine restartfile (ice_ic) use icedrv_calendar, only: istep0, istep1, time, time_forc use icepack_intfc, only: icepack_aggregate - use icedrv_domain_size, only: nilyr, nslyr, ncat - use icedrv_domain_size, only: max_ntrcr, nx + use icedrv_domain_size, only: nilyr, nslyr, ncat, nfsd, nblyr + use icedrv_domain_size, only: max_ntrcr, nx, n_iso, n_aero use icedrv_flux, only: swvdr, swvdf, swidr, swidf use icedrv_flux, only: sst, frzmlt, scale_factor + use icedrv_flux, only: frz_onset,fsnow use icedrv_forcing, only: oceanmixed_ice use icedrv_init, only: tmask use icedrv_state, only: trcr_depend, aice, vice, vsno, trcr use icedrv_state, only: aice0, aicen, vicen, vsnon, trcrn, aice_init use icedrv_state, only: trcr_base, nt_strata, n_trcr_strata + use icedrv_restart_shared, only: restart_format + use icedrv_arrays_column, only: dhsn, ffracn, hin_max + use icedrv_arrays_column, only: first_ice, first_ice_real + use icepack_tracers, only: ntrcr, nbtrcr character (*), optional :: ice_ic @@ -185,20 +281,14 @@ subroutine restartfile (ice_ic) character(len=char_len_long) :: filename character(len=*), parameter :: subname='(restartfile)' - if (present(ice_ic)) then - filename = trim(ice_ic) - else - call icedrv_system_abort(string=subname//'no ice_ic present', & - file=__FILE__,line=__LINE__) - endif - - write(nu_diag,*) 'Using restart dump=', trim(filename) - open(nu_restart,file=filename,form='unformatted') - read (nu_restart) istep0,time,time_forc - write(nu_diag,*) 'Restart read at istep=',istep0,time,time_forc - + ! local variable for reading from a netcdf file + integer (kind=int_kind) :: & + status + + ! Query tracers call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & - nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) + nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_iso_out=tr_iso, & tr_brine_out=tr_brine, & @@ -208,6 +298,40 @@ subroutine restartfile (ice_ic) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) + if (present(ice_ic)) then + filename = trim(ice_ic) + else + call icedrv_system_abort(string=subname//'no ice_ic present', & + file=__FILE__,line=__LINE__) + endif + + write(nu_diag,*) 'Using restart dump=', trim(filename) + + if (restart_format == 'bin') then + open(nu_restart,file=filename,form='unformatted') + read (nu_restart) istep0,time,time_forc + else if (restart_format == 'nc') then +#ifdef USE_NETCDF + ! set this to .true. for netcdf diagnostic output + diag = .false. + ! Open restart files + 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__) + + status = nf90_get_att(ncid, nf90_global, 'istep1', istep0) + status = nf90_get_att(ncid, nf90_global, 'time', time) + status = nf90_get_att(ncid, nf90_global, 'time_forc', time_forc) +#else + call icedrv_system_abort(string=subname//' ERROR: restart_format = "nc" requires USE_NETCDF',file=__FILE__,line=__LINE__) +#endif + else + call icedrv_system_abort(string=subname//'unrecognized restart format', & + file=__FILE__,line=__LINE__) + endif + + write(nu_diag,*) 'Restart read at istep=',istep0,time,time_forc + istep1 = istep0 !----------------------------------------------------------------- @@ -217,24 +341,30 @@ subroutine restartfile (ice_ic) !----------------------------------------------------------------- write(nu_diag,*) 'min/max area, vol ice, vol snow, Tsfc' - call read_restart_field(nu_restart,aicen,ncat) - call read_restart_field(nu_restart,vicen,ncat) - call read_restart_field(nu_restart,vsnon,ncat) - call read_restart_field(nu_restart,trcrn(:,nt_Tsfc,:),ncat) + call read_restart_field(nu_restart,aicen,ncat,'aicen') + call read_restart_field(nu_restart,vicen,ncat,'vicen') + call read_restart_field(nu_restart,vsnon,ncat,'vsnon') + call read_restart_field(nu_restart,trcrn(:,nt_Tsfc,:),ncat,'Tsfcn') write(nu_diag,*) 'min/max sice for each layer' do k=1,nilyr - call read_restart_field(nu_restart,trcrn(:,nt_sice+k-1,:),ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart,trcrn(:,nt_sice+k-1,:),ncat, & + 'sice'//trim(nchar)) enddo write(nu_diag,*) 'min/max qice for each layer' do k=1,nilyr - call read_restart_field(nu_restart,trcrn(:,nt_qice+k-1,:),ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart,trcrn(:,nt_qice+k-1,:),ncat, & + 'qice'//trim(nchar)) enddo write(nu_diag,*) 'min/max qsno for each layer' do k=1,nslyr - call read_restart_field(nu_restart,trcrn(:,nt_qsno+k-1,:),ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart,trcrn(:,nt_qsno+k-1,:),ncat, & + 'qsno'//trim(nchar)) enddo !----------------------------------------------------------------- @@ -243,11 +373,11 @@ subroutine restartfile (ice_ic) write(nu_diag,*) 'min/max radiation fields' - call read_restart_field(nu_restart,scale_factor,1) - call read_restart_field(nu_restart,swvdr,1) - call read_restart_field(nu_restart,swvdf,1) - call read_restart_field(nu_restart,swidr,1) - call read_restart_field(nu_restart,swidf,1) + call read_restart_field(nu_restart,scale_factor,1,'scale_factor') + call read_restart_field(nu_restart,swvdr,1,'swvdr') + call read_restart_field(nu_restart,swvdf,1,'swvdf') + call read_restart_field(nu_restart,swidr,1,'swidr') + call read_restart_field(nu_restart,swidf,1,'swidf') !----------------------------------------------------------------- ! for mixed layer model @@ -255,8 +385,8 @@ subroutine restartfile (ice_ic) if (oceanmixed_ice) then write(nu_diag,*) 'min/max sst, frzmlt' - call read_restart_field(nu_restart,sst,1) - call read_restart_field(nu_restart,frzmlt,1) + call read_restart_field(nu_restart,sst,1,'sst') + call read_restart_field(nu_restart,frzmlt,1,'frzmlt') endif ! tracers @@ -315,7 +445,7 @@ end subroutine restartfile ! Reads a single restart field ! author Chris Newman, LANL - subroutine read_restart_field(nu,work,ndim) + subroutine read_restart_field(nu,work,ndim,name) use icedrv_domain_size, only: nx @@ -325,6 +455,9 @@ subroutine read_restart_field(nu,work,ndim) real (kind=dbl_kind), dimension(nx,ndim), intent(inout) :: & work ! input array (real, 8-byte) + + character (len=*), intent(in), optional :: & + name ! local variables @@ -339,10 +472,25 @@ subroutine read_restart_field(nu,work,ndim) character(len=*), parameter :: subname='(read_restart_field)' - do n = 1, ndim - read(nu) (work2(i), i=1,nx) - work(:,n) = work2(:) - enddo + if (restart_format == 'bin') then + do n = 1, ndim + read(nu) (work2(i), i=1,nx) + work(:,n) = work2(:) + enddo + else if (restart_format == 'nc') then +#ifdef USE_NETCDF + if (ndim == 1) then + call read_restart_field_net1D(ncid,1,work,name,1,diag) + else + call read_restart_field_net2D(ncid,1,work,name,ndim,diag) + endif +#else + call icedrv_system_abort(string=subname//' ERROR: restart_format = "nc" requires USE_NETCDF',file=__FILE__,line=__LINE__) +#endif + else + call icedrv_system_abort(string=subname//'unrecognized restart format', & + file=__FILE__,line=__LINE__) + endif minw = minval(work) maxw = maxval(work) @@ -356,7 +504,7 @@ end subroutine read_restart_field ! Writes a single restart field. ! author Chris Newman, LANL - subroutine write_restart_field(nu,work,ndim) + subroutine write_restart_field(nu,work,ndim,vname,dims) use icedrv_domain_size, only: nx @@ -366,6 +514,12 @@ subroutine write_restart_field(nu,work,ndim) real (kind=dbl_kind), dimension(nx,ndim), intent(in) :: & work ! input array (real, 8-byte) + + character (len=*), intent(in), optional :: & + vname ! variable name for netcdf output + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs ! local variables @@ -380,10 +534,43 @@ subroutine write_restart_field(nu,work,ndim) character(len=*), parameter :: subname='(write_restart_field)' - do n = 1, ndim - work2(:) = work(:,n) - write(nu) (work2(i), i=1,nx) - enddo +#ifdef USE_NETCDF + ! local variables if we're writing in NetCDF + integer (kind=int_kind) :: & + status, & ! Status variable from netCDF routine + varid ! Variable ID +#endif + + if (restart_format == 'bin') then + do n = 1, ndim + work2(:) = work(:,n) + write(nu) (work2(i), i=1,nx) + enddo + else if (restart_format == 'nc') then +#ifdef USE_NETCDF + ! Check whether the variable already exists in the netcdf file + status = nf90_inq_varid(ncid,trim(vname),varid) + if (status == nf90_enotvar) then ! if not add it to the file + status = nf90_redef(ncid) ! Put file into define mode + if (ndim == 1) then + ! Handle 1D cases + status = nf90_def_var(ncid,trim(vname),nf90_double,dims(1:1),varid) + else + status = nf90_def_var(ncid,trim(vname),nf90_double,dims,varid) + endif + endif + status = nf90_enddef(ncid) ! put file into data mode + if (ndim == 1) then + call write_restart_field_nc1D(ncid,0,work,'ruf8',vname,1,diag) + else + call write_restart_field_nc2D(ncid,0,work,'ruf8',vname,ndim,diag) + endif +#else + call icedrv_system_abort(string=subname//' ERROR: restart_format = "nc" requires USE_NETCDF',file=__FILE__,line=__LINE__) +#endif + else + write (nu_diag,*) 'WARNING: Restart format must be either "bin" or "nc", no restart file written' + endif minw = minval(work) maxw = maxval(work) @@ -415,11 +602,14 @@ end subroutine final_restart ! authors Elizabeth C. Hunke, LANL ! David A. Bailey, NCAR - subroutine write_restart_pond_topo() + subroutine write_restart_pond_topo(dims) use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_apnd, nt_hpnd, nt_ipnd character(len=*), parameter :: subname='(write_restart_pond_topo)' @@ -429,9 +619,9 @@ subroutine write_restart_pond_topo() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - call write_restart_field(nu_dump,trcrn(:,nt_apnd,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_hpnd,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_ipnd,:),ncat) + call write_restart_field(nu_dump,trcrn(:,nt_apnd,:),ncat,'apnd',dims) + call write_restart_field(nu_dump,trcrn(:,nt_hpnd,:),ncat,'hpnd',dims) + call write_restart_field(nu_dump,trcrn(:,nt_ipnd,:),ncat,'ipnd',dims) end subroutine write_restart_pond_topo @@ -457,9 +647,9 @@ subroutine read_restart_pond_topo() write(nu_diag,*) 'min/max topo ponds' - call read_restart_field(nu_restart,trcrn(:,nt_apnd,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_hpnd,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_ipnd,:),ncat) + call read_restart_field(nu_restart,trcrn(:,nt_apnd,:),ncat,'apnd') + call read_restart_field(nu_restart,trcrn(:,nt_hpnd,:),ncat,'hpnd') + call read_restart_field(nu_restart,trcrn(:,nt_ipnd,:),ncat,'ipnd') end subroutine read_restart_pond_topo @@ -469,11 +659,14 @@ end subroutine read_restart_pond_topo ! ! authors Elizabeth C. Hunke, LANL - subroutine write_restart_snow() + subroutine write_restart_snow(dims) use icedrv_state, only: trcrn use icedrv_domain_size, only: nslyr, ncat + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_smice, nt_smliq, nt_rhos, nt_rsnw, k character(len=*), parameter :: subname='(write_restart_snow)' @@ -484,10 +677,11 @@ subroutine write_restart_snow() file=__FILE__,line= __LINE__) do k=1,nslyr - call write_restart_field(nu_dump,trcrn(:,nt_smice+k-1,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_smliq+k-1,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_rhos +k-1,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_rsnw +k-1,:),ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump,trcrn(:,nt_smice+k-1,:),ncat,'smice'//trim(nchar),dims) + call write_restart_field(nu_dump,trcrn(:,nt_smliq+k-1,:),ncat,'smliq'//trim(nchar),dims) + call write_restart_field(nu_dump,trcrn(:,nt_rhos +k-1,:),ncat,'rhos'//trim(nchar),dims) + call write_restart_field(nu_dump,trcrn(:,nt_rsnw +k-1,:),ncat,'rsnw'//trim(nchar),dims) enddo end subroutine write_restart_snow @@ -515,10 +709,11 @@ subroutine read_restart_snow() write(nu_diag,*) 'min/max snow metamorphosis tracers' do k=1,nslyr - call read_restart_field(nu_restart,trcrn(:,nt_smice+k-1,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_smliq+k-1,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_rhos +k-1,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_rsnw +k-1,:),ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart,trcrn(:,nt_smice+k-1,:),ncat,'smice'//trim(nchar)) + call read_restart_field(nu_restart,trcrn(:,nt_smliq+k-1,:),ncat,'smliq'//trim(nchar)) + call read_restart_field(nu_restart,trcrn(:,nt_rhos +k-1,:),ncat,'rhos'//trim(nchar)) + call read_restart_field(nu_restart,trcrn(:,nt_rsnw +k-1,:),ncat,'rsnw'//trim(nchar)) enddo end subroutine read_restart_snow @@ -528,10 +723,14 @@ end subroutine read_restart_snow ! Dumps all values needed for restarting ! author Elizabeth C. Hunke, LANL - subroutine write_restart_age() - + subroutine write_restart_age(dims) + use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_iage character(len=*), parameter :: subname='(write_restart_age)' @@ -540,7 +739,7 @@ subroutine write_restart_age() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - call write_restart_field(nu_dump,trcrn(:,nt_iage,:),ncat) + call write_restart_field(nu_dump,trcrn(:,nt_iage,:),ncat,'iage',dims) end subroutine write_restart_age @@ -562,8 +761,8 @@ subroutine read_restart_age() call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - - call read_restart_field(nu_restart,trcrn(:,nt_iage,:),ncat) + + call read_restart_field(nu_restart,trcrn(:,nt_iage,:),ncat,'iage') end subroutine read_restart_age @@ -572,10 +771,14 @@ end subroutine read_restart_age ! Dumps all values needed for restarting ! author Elizabeth C. Hunke, LANL - subroutine write_restart_fsd() + subroutine write_restart_fsd(dims) use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat, nfsd + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_fsd, k character(len=*), parameter :: subname='(write_restart_fsd)' @@ -585,7 +788,9 @@ subroutine write_restart_fsd() file=__FILE__,line= __LINE__) do k = 1, nfsd - call write_restart_field(nu_dump,trcrn(:,nt_fsd+k-1,:),ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump,trcrn(:,nt_fsd+k-1,:),ncat,& + 'fsd'//trim(nchar),dims) end do end subroutine write_restart_fsd @@ -610,7 +815,8 @@ subroutine read_restart_fsd() file=__FILE__,line= __LINE__) do k =1, nfsd - call read_restart_field(nu_restart,trcrn(:,nt_fsd+k-1,:),ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart,trcrn(:,nt_fsd+k-1,:),ncat,'fsd'//trim(nchar)) end do end subroutine read_restart_fsd @@ -620,11 +826,15 @@ end subroutine read_restart_fsd ! Dumps all values needed for restarting ! author Elizabeth C. Hunke, LANL - subroutine write_restart_FY() + subroutine write_restart_FY(dims) use icedrv_flux, only: frz_onset use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_FY character(len=*), parameter :: subname='(write_restart_FY)' @@ -633,8 +843,8 @@ subroutine write_restart_FY() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - call write_restart_field(nu_dump,trcrn(:,nt_FY,:),ncat) - call write_restart_field(nu_dump,frz_onset,1) + call write_restart_field(nu_dump,trcrn(:,nt_FY,:),ncat,'FY',dims) + call write_restart_field(nu_dump,frz_onset,1,'frz_onset',dims) end subroutine write_restart_FY @@ -658,11 +868,11 @@ subroutine read_restart_FY() write(nu_diag,*) 'min/max first-year ice area' - call read_restart_field(nu_restart,trcrn(:,nt_FY,:),ncat) + call read_restart_field(nu_restart,trcrn(:,nt_FY,:),ncat,'FY') write(nu_diag,*) 'min/max frz_onset' - call read_restart_field(nu_restart,frz_onset,1) + call read_restart_field(nu_restart,frz_onset,1,'frz_onset') end subroutine read_restart_FY @@ -672,10 +882,14 @@ end subroutine read_restart_FY ! ! author Elizabeth C. Hunke, LANL - subroutine write_restart_lvl() + subroutine write_restart_lvl(dims) use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_alvl, nt_vlvl character(len=*), parameter :: subname='(write_restart_lvl)' @@ -684,8 +898,8 @@ subroutine write_restart_lvl() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - call write_restart_field(nu_dump,trcrn(:,nt_alvl,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_vlvl,:),ncat) + call write_restart_field(nu_dump,trcrn(:,nt_alvl,:),ncat,'alvl',dims) + call write_restart_field(nu_dump,trcrn(:,nt_vlvl,:),ncat,'vlvl',dims) end subroutine write_restart_lvl @@ -709,8 +923,8 @@ subroutine read_restart_lvl() write(nu_diag,*) 'min/max level ice area, volume' - call read_restart_field(nu_restart,trcrn(:,nt_alvl,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_vlvl,:),ncat) + call read_restart_field(nu_restart,trcrn(:,nt_alvl,:),ncat,'alvl') + call read_restart_field(nu_restart,trcrn(:,nt_vlvl,:),ncat,'vlvl') end subroutine read_restart_lvl @@ -720,12 +934,16 @@ end subroutine read_restart_lvl ! ! authors Elizabeth C. Hunke, LANL - subroutine write_restart_pond_lvl() + subroutine write_restart_pond_lvl(dims) use icedrv_arrays_column, only: dhsn, ffracn use icedrv_flux, only: fsnow use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + integer (kind=int_kind) :: nt_apnd, nt_hpnd, nt_ipnd character(len=*), parameter :: subname='(write_restart_pond_lvl)' @@ -735,12 +953,12 @@ subroutine write_restart_pond_lvl() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - call write_restart_field(nu_dump,trcrn(:,nt_apnd,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_hpnd,:),ncat) - call write_restart_field(nu_dump,trcrn(:,nt_ipnd,:),ncat) - call write_restart_field(nu_dump,fsnow(:),1) - call write_restart_field(nu_dump,dhsn(:,:),ncat) - call write_restart_field(nu_dump,ffracn(:,:),ncat) + call write_restart_field(nu_dump,trcrn(:,nt_apnd,:),ncat,'apnd',dims) + call write_restart_field(nu_dump,trcrn(:,nt_hpnd,:),ncat,'hpnd',dims) + call write_restart_field(nu_dump,trcrn(:,nt_ipnd,:),ncat,'ipnd',dims) + call write_restart_field(nu_dump,fsnow(:),1,'fsnow',dims) + call write_restart_field(nu_dump,dhsn(:,:),ncat,'dhsn',dims) + call write_restart_field(nu_dump,ffracn(:,:),ncat,'ffracn',dims) end subroutine write_restart_pond_lvl @@ -767,12 +985,12 @@ subroutine read_restart_pond_lvl() write(nu_diag,*) 'min/max level-ice ponds' - call read_restart_field(nu_restart,trcrn(:,nt_apnd,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_hpnd,:),ncat) - call read_restart_field(nu_restart,trcrn(:,nt_ipnd,:),ncat) - call read_restart_field(nu_restart,fsnow(:),1) - call read_restart_field(nu_restart,dhsn(:,:),ncat) - call read_restart_field(nu_restart,ffracn(:,:),ncat) + call read_restart_field(nu_restart,trcrn(:,nt_apnd,:),ncat,'apnd') + call read_restart_field(nu_restart,trcrn(:,nt_hpnd,:),ncat,'hpnd') + call read_restart_field(nu_restart,trcrn(:,nt_ipnd,:),ncat,'ipnd') + call read_restart_field(nu_restart,fsnow(:),1,'fsnow') + call read_restart_field(nu_restart,dhsn(:,:),ncat,'dhsn') + call read_restart_field(nu_restart,ffracn(:,:),ncat,'ffracn') end subroutine read_restart_pond_lvl @@ -784,11 +1002,14 @@ end subroutine read_restart_pond_lvl ! David Bailey, NCAR ! Marika Holland, NCAR - subroutine write_restart_aero() + subroutine write_restart_aero(dims) use icedrv_domain_size, only: n_aero use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs ! local variables @@ -805,10 +1026,15 @@ subroutine write_restart_aero() write(nu_diag,*) 'write_restart_aero (aerosols)' do k = 1, n_aero - call write_restart_field(nu_dump, trcrn(:,nt_aero +(k-1)*4,:), ncat) - call write_restart_field(nu_dump, trcrn(:,nt_aero+1+(k-1)*4,:), ncat) - call write_restart_field(nu_dump, trcrn(:,nt_aero+2+(k-1)*4,:), ncat) - call write_restart_field(nu_dump, trcrn(:,nt_aero+3+(k-1)*4,:), ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump, trcrn(:,nt_aero +(k-1)*4,:),ncat,& + 'aerosnossl'//trim(nchar),dims) + call write_restart_field(nu_dump, trcrn(:,nt_aero+1+(k-1)*4,:),ncat,& + 'aerosnoint'//trim(nchar),dims) + call write_restart_field(nu_dump, trcrn(:,nt_aero+2+(k-1)*4,:),ncat,& + 'aeroicessl'//trim(nchar),dims) + call write_restart_field(nu_dump, trcrn(:,nt_aero+3+(k-1)*4,:),ncat,& + 'aeroiceint'//trim(nchar),dims) enddo end subroutine write_restart_aero @@ -844,10 +1070,11 @@ subroutine read_restart_aero() write(nu_diag,*) 'read_restart_aero (aerosols)' do k = 1, n_aero - call read_restart_field(nu_restart, trcrn(:,nt_aero +(k-1)*4,:), ncat) - call read_restart_field(nu_restart, trcrn(:,nt_aero+1+(k-1)*4,:), ncat) - call read_restart_field(nu_restart, trcrn(:,nt_aero+2+(k-1)*4,:), ncat) - call read_restart_field(nu_restart, trcrn(:,nt_aero+3+(k-1)*4,:), ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart, trcrn(:,nt_aero +(k-1)*4,:),ncat,'aerosnossl'//trim(nchar)) + call read_restart_field(nu_restart, trcrn(:,nt_aero+1+(k-1)*4,:),ncat,'aerosnoint'//trim(nchar)) + call read_restart_field(nu_restart, trcrn(:,nt_aero+2+(k-1)*4,:),ncat,'aeroicessl'//trim(nchar)) + call read_restart_field(nu_restart, trcrn(:,nt_aero+3+(k-1)*4,:),ncat,'aeroiceint'//trim(nchar)) enddo end subroutine read_restart_aero @@ -860,11 +1087,14 @@ end subroutine read_restart_aero ! David Bailey, NCAR ! Marika Holland, NCAR - subroutine write_restart_iso() + subroutine write_restart_iso(dims) use icedrv_domain_size, only: n_iso use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat + + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs ! local variables @@ -882,8 +1112,11 @@ subroutine write_restart_iso() write(nu_diag,*) 'write_restart_iso (isotopes)' do k = 1, n_iso - call write_restart_field(nu_dump, trcrn(:,nt_isosno+(k-1),:), ncat) - call write_restart_field(nu_dump, trcrn(:,nt_isoice+(k-1),:), ncat) + write(nchar,'(i3.3)') k + call write_restart_field(nu_dump, trcrn(:,nt_isosno+(k-1),:),ncat,& + 'isosno'//trim(nchar),dims) + call write_restart_field(nu_dump, trcrn(:,nt_isoice+(k-1),:),ncat,& + 'isoice'//trim(nchar),dims) enddo end subroutine write_restart_iso @@ -920,15 +1153,16 @@ subroutine read_restart_iso() write(nu_diag,*) 'read_restart_iso (isotopes)' do k = 1, n_iso - call read_restart_field(nu_restart, trcrn(:,nt_isosno+(k-1),:), ncat) - call read_restart_field(nu_restart, trcrn(:,nt_isoice+(k-1),:), ncat) + write(nchar,'(i3.3)') k + call read_restart_field(nu_restart, trcrn(:,nt_isosno+(k-1),:),ncat,'isosno'//trim(nchar)) + call read_restart_field(nu_restart, trcrn(:,nt_isoice+(k-1),:),ncat,'isoice'//trim(nchar)) enddo end subroutine read_restart_iso !======================================================================= - subroutine write_restart_hbrine() + subroutine write_restart_hbrine(dims) ! Dumps all values needed for a hbrine restart ! author Elizabeth C. Hunke, LANL @@ -937,6 +1171,9 @@ subroutine write_restart_hbrine() use icedrv_state, only: trcrn use icedrv_domain_size, only: ncat, nx + integer (kind=int_kind), intent(in), optional :: & + dims(:) ! netcdf dimension IDs + ! local variables integer (kind=int_kind) :: & @@ -959,8 +1196,8 @@ subroutine write_restart_hbrine() enddo ! n enddo ! i - call write_restart_field(nu_dump,trcrn(:,nt_fbri,:),ncat) - call write_restart_field(nu_dump,first_ice_real(:,:),ncat) + call write_restart_field(nu_dump,trcrn(:,nt_fbri,:),ncat,'fbri',dims) + call write_restart_field(nu_dump,first_ice_real(:,:),ncat,'first_ice',dims) end subroutine write_restart_hbrine @@ -989,8 +1226,8 @@ subroutine read_restart_hbrine() write(nu_diag,*) 'read brine restart' - call read_restart_field(nu_restart,trcrn(:,nt_fbri,:),ncat) - call read_restart_field(nu_restart,first_ice_real(:,:),ncat) + call read_restart_field(nu_restart,trcrn(:,nt_fbri,:),ncat,'fbri') + call read_restart_field(nu_restart,first_ice_real(:,:),ncat,'first_ice') do i = 1, nx do n = 1, ncat @@ -1005,7 +1242,410 @@ subroutine read_restart_hbrine() end subroutine read_restart_hbrine !======================================================================= +#ifdef USE_NETCDF + subroutine define_rest_field(ncid, vname, dims) + +! Defines a field in NetCDF restart file +! author Chris Riedel, NCAR + + character (len=*) , intent(in) :: vname + integer (kind=int_kind), intent(in) :: dims(:) + integer (kind=int_kind), intent(in) :: ncid + + integer (kind=int_kind) :: varid + + integer (kind=int_kind) :: & + status ! status variable from netCDF routine + + status = nf90_def_var(ncid,trim(vname),nf90_double,dims,varid) + + end subroutine define_rest_field +!======================================================================= + subroutine write_restart_field_nc2D(ncid,nrec,work,atype,vname,ndim3,diag) + +! Write a 2D field (nx, ncat) in NetCDF restart +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + !use ice_fileunits, only: nu_diag + !use ice_read_write, only: ice_write, ice_write_nc + + integer (kind=int_kind), intent(in) :: & + ncid , & ! unit number + ndim3 , & ! third dimension + nrec ! record number (0 for sequential access) + real (kind=dbl_kind), intent(in), dimension(nx,ncat) :: & + work ! input array (real, 8-byte) + character (len=4), intent(in) :: & + atype ! format for output array + ! (real/integer, 4-byte/8-byte) + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: vname + + ! local variables + + integer (kind=int_kind) :: & + n, & ! dimension counter + varid, & ! variable id + status ! status variable from netCDF routine + + + status = nf90_inq_varid(ncid,trim(vname),varid) + write(nu_diag,*) 'Writing out ',trim(vname) + call ice_write_nc2D(ncid, 1, varid, work, diag) + + end subroutine write_restart_field_nc2D +!======================================================================= + subroutine write_restart_field_nc1D(ncid,nrec,work,atype,vname,ndim3,diag) + +! Write a 1D field (nx) in NetCDF restart +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + ncid , & ! unit number + ndim3 , & ! third dimension + nrec ! record number (0 for sequential access) + real (kind=dbl_kind), intent(in), dimension(nx) :: & + work ! input array (real, 8-byte) + character (len=4), intent(in) :: & + atype ! format for output array + ! (real/integer, 4-byte/8-byte) + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: vname + + ! local variables + + integer (kind=int_kind) :: & + n, & ! dimension counter + varid, & ! variable id + status ! status variable from netCDF routine + + + status = nf90_inq_varid(ncid,trim(vname),varid) + if (status /= 0) then + write(nu_diag,*) 'Writing out ',trim(vname) + write(nu_diag,*) 'Erros Status ',status + call icedrv_system_abort(string='Write out restart', & + file=__FILE__,line= __LINE__) + else + write(nu_diag,*) 'Writing out ',trim(vname) + endif + call ice_write_nc1D(ncid, 1, varid, work, diag) + end subroutine write_restart_field_nc1D +!======================================================================= + subroutine ice_write_nc2D(fid, nrec, varid, work, diag) + +! Write a 2D field (nx, ncat) in NetCDF restart +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + varid , & ! variable id + nrec ! record number + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + real (kind=dbl_kind), intent(in), dimension(nx,ncat) :: & + work ! output array (real, 8-byte) + + ! local variables + + integer (kind=int_kind) :: & + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + id, & ! dimension index + dimlen ! size of dimension + integer (kind=int_kind) :: start2(2),count2(2) + + real (kind=dbl_kind) :: & + amin, amax, asum ! min, max values and sum of input array + + character (char_len) :: & + dimname ! dimension name + + integer (kind=int_kind) :: ny,numDims + + ny = ncat + start2(1) = 1 + count2(1) = nx + start2(2) = 1 + count2(2) = ncat + + status = nf90_put_var( fid, varid, work, & + start=start2, & + count=count2) + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (diag) then + amin = minval(work) + amax = maxval(work) + !asum = sum (work) + write(nu_diag,*) ' min, max =', amin, amax + endif + + !deallocate(work) + + end subroutine ice_write_nc2D +!======================================================================= + subroutine ice_write_nc1D(fid, nrec, varid, work, diag) + +! Write a 1D field (nx) in NetCDF restart +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + varid , & ! variable id + nrec ! record number + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + real (kind=dbl_kind), intent(in), dimension(nx) :: & + work ! output array (real, 8-byte) + + ! local variables + + integer (kind=int_kind) :: & + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + id, & ! dimension index + dimlen ! size of dimension + integer (kind=int_kind) :: start2(1),count2(1) + + real (kind=dbl_kind) :: & + amin, amax, asum ! min, max values and sum of input array + + character (char_len) :: & + dimname ! dimension name + + integer (kind=int_kind) :: ny,numDims + + ny = ncat + start2(1) = 1 + count2(1) = nx + + status = nf90_put_var( fid, varid, work, & + start=start2, & + count=count2) + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (diag) then + amin = minval(work) + amax = maxval(work) + !asum = sum (work) + write(nu_diag,*) ' min, max =', amin, amax + endif + + !deallocate(work) + + end subroutine ice_write_nc1D +!======================================================================= + subroutine read_restart_field_net2D(nu,nrec,work,vname,ndim3, & + diag) + +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + nu , & ! unit number (not used for netcdf) + ndim3 , & ! third dimension + nrec ! record number (0 for sequential access) + + real (kind=dbl_kind), intent(inout), dimension(nx,ncat) :: & + work ! input array (real, 8-byte) + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: vname + + ! local variables + + integer (kind=int_kind) :: & + n, & ! number of dimensions for variable + varid, & ! variable id + status ! status variable from netCDF routine + + call ice_read_nc2D(ncid, 1, vname, work, diag) + + end subroutine read_restart_field_net2D +!======================================================================= + subroutine read_restart_field_net1D(nu,nrec,work,vname,ndim3, & + diag) + +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + nu , & ! unit number (not used for netcdf) + ndim3 , & ! third dimension + nrec ! record number (0 for sequential access) + + real (kind=dbl_kind), intent(inout), dimension(nx) :: & + work ! input array (real, 8-byte) + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: vname + + ! local variables + + integer (kind=int_kind) :: & + n, & ! number of dimensions for variable + varid, & ! variable id + status ! status variable from netCDF routine + + call ice_read_nc1D(ncid, 1, vname, work, diag) + + end subroutine read_restart_field_net1D +!======================================================================= + subroutine ice_read_nc2D(fid, nrec, varname, work, diag) + +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + nrec ! record number + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), intent(out), dimension(nx,ncat) :: & + work ! output array (real, 8-byte) + + integer (kind=int_kind) :: & + varid , & ! variable id + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + id, & ! dimension index + dimlen ! size of dimension + + real (kind=dbl_kind) :: & + amin, amax, asum ! min, max values and sum of input array + + character (char_len) :: & + dimname ! dimension name + + integer (kind=int_kind) :: start2(2),count2(2) + + integer (kind=int_kind) :: ny,numDims + + ny = ncat + start2(1) = 1 + count2(1) = nx + start2(2) = 1 + count2(2) = ncat + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call icedrv_system_abort(string='ice_read_nc2D:Var not found '//trim(varname), & + file=__FILE__,line=__LINE__) + endif + + status = nf90_get_var( fid, varid, work, & + start=start2, & + count=count2 ) + + + if (diag) then + amin = minval(work) + amax = maxval(work) + asum = sum (work) + write(nu_diag,*) ' min, max, sum =', amin, amax, asum + endif + + end subroutine ice_read_nc2D +!======================================================================= + subroutine ice_read_nc1D(fid, nrec, varname, work, diag) + +! author Chris Riedel, NCAR + + use icedrv_domain_size, only: ncat,nx + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + nrec ! record number + + logical (kind=log_kind), intent(in) :: & + diag ! if true, write diagnostic output + + character (len=*), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), intent(out), dimension(nx) :: & + work ! output array (real, 8-byte) + + integer (kind=int_kind) :: & + varid , & ! variable id + status, & ! status output from netcdf routines + ndim, nvar, & ! sizes of netcdf file + id, & ! dimension index + dimlen ! size of dimension + + real (kind=dbl_kind) :: & + amin, amax, asum ! min, max values and sum of input array + + character (char_len) :: & + dimname ! dimension name + + integer (kind=int_kind) :: start2(1),count2(1) + + integer (kind=int_kind) :: ny,numDims + + ny = ncat + start2(1) = 1 + count2(1) = nx + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call icedrv_system_abort(string='ice_read_nc2D:Var not found'//trim(varname), & + file=__FILE__,line=__LINE__) + endif + + status = nf90_get_var( fid, varid, work, & + start=start2, & + count=count2 ) + + if (diag) then + amin = minval(work) + amax = maxval(work) + asum = sum (work) + write(nu_diag,*) ' min, max, sum =', amin, amax, asum + endif + + end subroutine ice_read_nc1D +#endif +!======================================================================= end module icedrv_restart !======================================================================= diff --git a/configuration/driver/icedrv_restart_shared.F90 b/configuration/driver/icedrv_restart_shared.F90 index 43b135d11..d4078a93a 100644 --- a/configuration/driver/icedrv_restart_shared.F90 +++ b/configuration/driver/icedrv_restart_shared.F90 @@ -14,6 +14,9 @@ module icedrv_restart_shared restart_file , & ! output file for restart dump restart_dir ! directory name for restart dump + character (len=char_len), public :: & + restart_format ! format of restart files 'nc' + !======================================================================= contains diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index 505dad4c6..bc9e2942c 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -9,11 +9,12 @@ ice_ic = 'default' restart = .false. restart_dir = './restart/' + restart_format = 'bin' dumpfreq = 'y' dump_last = .false. diagfreq = 24 diag_file = 'ice_diag' - history_cdf = .false. + history_format = 'none' cpl_bgc = .false. conserv_check = .false. / diff --git a/configuration/scripts/options/set_nml.histcdf b/configuration/scripts/options/set_nml.histcdf new file mode 100644 index 000000000..a1b6d4289 --- /dev/null +++ b/configuration/scripts/options/set_nml.histcdf @@ -0,0 +1 @@ +history_format = 'nc' diff --git a/configuration/scripts/options/set_nml.ionetcdf b/configuration/scripts/options/set_nml.ionetcdf deleted file mode 100644 index 6aaa620e7..000000000 --- a/configuration/scripts/options/set_nml.ionetcdf +++ /dev/null @@ -1 +0,0 @@ -history_cdf = .true. diff --git a/configuration/scripts/options/set_nml.restcdf b/configuration/scripts/options/set_nml.restcdf new file mode 100644 index 000000000..014a61c4e --- /dev/null +++ b/configuration/scripts/options/set_nml.restcdf @@ -0,0 +1 @@ +restart_format = 'nc' diff --git a/configuration/scripts/tests/io_suite.ts b/configuration/scripts/tests/io_suite.ts index 2af2980d5..cff8d809b 100644 --- a/configuration/scripts/tests/io_suite.ts +++ b/configuration/scripts/tests/io_suite.ts @@ -1,2 +1,9 @@ restart col 1x1 debug,ionetcdf smoke col 1x1 run1year,diag1,ionetcdf +restart col 1x1 debug,ionetcdf,histcdf +smoke col 1x1 run1year,diag1,ionetcdf,histcdf +restart col 1x1 debug,ionetcdf,restcdf +smoke col 1x1 run1year,diag1,ionetcdf,restcdf +restart col 1x1 debug,ionetcdf,histcdf,restcdf +smoke col 1x1 run1year,diag1,ionetcdf,histcdf,restcdf + diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index 439e4c74e..08a3b0d94 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -210,7 +210,7 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "highfreq", ":math:`\bullet` high-frequency atmo coupling", "F" "hin_old", "ice thickness prior to growth/melt", "m" "hin_max", "category thickness limits", "m" - "history_cdf", "flag to turn on netcdf history output", "F" + "history_format", "turns on netcdf history output if set to 'nc'", "" "hmix", "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" @@ -383,6 +383,7 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "restart_bgc", ":math:`\bullet` if true, read bgc restart file", "" "restart_dir", ":math:`\bullet` path to restart/dump files", "" "restart_file", ":math:`\bullet` restart file prefix", "" + "restart_format", "history files are read/written in binary or netcdf format if set to 'bin' or 'nc' respectively", "bin" "restart_[tracer]", ":math:`\bullet` if true, read tracer restart file", "" "restore_bgc", ":math:`\bullet` if true, restore nitrate/silicate to data", "" "restore_ice", ":math:`\bullet` if true, restore ice state along lateral boundaries", "" diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 85e75d4ce..6d4f877ca 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -140,7 +140,8 @@ 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." - "``history_cdf``", "logical", "netcdf history output", "``.false.``" + "``history_format``", "``cdf``", "history file output in netcdf format", "``none``" + "","``none``","no history output","" "``ice_ic``", "``default``", "latitude and sst dependent initial condition", "``default``" "", "``none``", "no ice", "" "", "'path/file'", "restart file name", "" @@ -150,6 +151,8 @@ setup_nml "``restart``", "logical", "initialize using restart file", "``.false.``" "``restart_dir``", "string", "path to restart directory", "'./'" "``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","" "``use_leap_years``", "logical", "include leap days", "``.false.``" "``year_init``", "integer", "the initial year if not using restart", "0" "", "", "", "" diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index b4239b87e..3af5df796 100755 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -165,7 +165,7 @@ Icepack writes diagnostic information for each grid cell as a separate file, Restart files ~~~~~~~~~~~~~ -Icepack provides restart data in binary unformatted format. The restart files +Icepack provides restart data in binary unformatted format or netcdf. The restart files created by the Icepack driver contain all of the variables needed for a full, exact restart. The filename begins with the character string ‘iced.’ and is placed in the directory specified by the namelist variable @@ -174,13 +174,15 @@ variable ``dumpfreq``. The namelist variable ``ice_ic`` contains the pointer to the filename from which the restart data is to be read and the namelist option ``restart`` must be set to ``.true.`` to use the file. ``dump_last`` namelist can also be set to true to trigger restarts automatically -at then end of runs. +at then end of runs. The default restart file format is binary and is set +by setting the ``restart_format`` namelist to 'bin'. Netcdf restart files +are set by setting ``restart_format`` namelist to 'nc'. History files ~~~~~~~~~~~~~ Icepack has a primitive netcdf history capability that is turned on with the -``history_cdf`` namelist. When ``history_cdf`` is set to true, history files +``history_format`` namelist. When ``history_format`` is set to 'nc', history files are created for each run with a naming convention of **icepack.h.yyyymmdd.nc** in the run directory history directory. The yyyymmdd is the start date for each run. @@ -197,7 +199,8 @@ written, the USE_NETCDF C preprocessor directive must be set during compilation. is done by setting ``ICE_IOTYPE`` to ``netcdf`` in **icepack.settings**. In addition, the machine env and Macros files must include support for compilation with NetCDF. The ``icepack.setup -s`` option ``ionetcdf`` will set the ICE_IOTYPE to netcdf, which turns on -the USE_NETCDF C preprocessor. ``ionetcdf`` also sets the ``history_cdf`` flag to true. +the USE_NETCDF C preprocessor. ``ionetcdf`` also sets the ``history_format`` and +``restart_format`` flags to 'nc'. .. _bgc-hist: From a4779cc71125b40a7db3a4da8512247cbf2b0955 Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Tue, 14 Mar 2023 19:27:54 -0600 Subject: [PATCH 3/4] FSD bug fixes (#424) * FSD bug fixes for lateral melt and weld * put back declare sicen * another small welding fix * Change print statements to use icepack_warnings * Fix some comments * Fix some comments * Fix some comments * Check for wlat present when tr_fsd = .true. * Fix some syntax errors --------- Co-authored-by: cmbitz --- columnphysics/icepack_fsd.F90 | 14 ++- columnphysics/icepack_therm_itd.F90 | 125 +++++++++++++---------- columnphysics/icepack_therm_vertical.F90 | 29 +++--- configuration/driver/icedrv_flux.F90 | 1 + configuration/driver/icedrv_step.F90 | 6 +- 5 files changed, 96 insertions(+), 79 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index bee46cd35..1454fdefc 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -699,7 +699,10 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & DO WHILE (elapsed_t.lt.dt) nsubt = nsubt + 1 - if (nsubt.gt.100) print *, 'latg not converging' + if (nsubt.gt.100) then + write(warnstr,*) subname,'latg not converging' + call icepack_warnings_add(warnstr) + endif ! finite differences df_flx(:) = c0 ! NB could stay zero if all in largest FS cat @@ -711,8 +714,6 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & df_flx(k) = f_flx(k+1) - f_flx(k) end do -! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0' - dafsd_tmp(:) = c0 do k = 1, nfsd dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) & @@ -931,7 +932,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, & gain, loss ! welding tendencies real(kind=dbl_kind) :: & - prefac , & ! multiplies kernel kern , & ! kernel subdt , & ! subcycling time step for stability (s) elapsed_t ! elapsed subcycling time @@ -942,7 +942,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, & afsdn (:,:) = c0 afsd_init(:) = c0 stability = c0 - prefac = p5 do n = 1, ncat @@ -986,8 +985,7 @@ subroutine fsd_weld_thermo (ncat, nfsd, & if (k > i) then kern = c_weld * floe_area_c(i) * aicen(n) loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j) - if (i.eq.j) prefac = c1 ! otherwise 0.5 - gain(k) = gain(k) + prefac*kern*afsd_tmp(i)*afsd_tmp(j) + gain(k) = gain(k) + kern*afsd_tmp(i)*afsd_tmp(j) end if end do end do @@ -1011,11 +1009,11 @@ subroutine fsd_weld_thermo (ncat, nfsd, & END DO ! time + afsdn(:,n) = afsd_tmp(:) call icepack_cleanup_fsdn (nfsd, afsdn(:,n)) if (icepack_warnings_aborted(subname)) return do k = 1, nfsd - afsdn(k,n) = afsd_tmp(k) trcrn(nt_fsd+k-1,n) = afsdn(k,n) ! history/diagnostics d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 93ed8bebd..3aac1ac70 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -884,7 +884,7 @@ subroutine lateral_melt (dt, ncat, & fhocn, faero_ocn, & fiso_ocn, & rside, meltl, & - fside, sss, & + fside, wlat, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & @@ -914,6 +914,9 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), intent(in) :: & rside , & ! fraction of ice that melts laterally + wlat ! lateral melt rate (m/s) + + real (kind=dbl_kind), intent(inout) :: & fside ! lateral heat flux (W/m^2) real (kind=dbl_kind), intent(inout) :: & @@ -955,7 +958,7 @@ subroutine lateral_melt (dt, ncat, & dfsalt , & ! change in fsalt dvssl , & ! snow surface layer volume dvint , & ! snow interior layer - cat1_arealoss, tmp ! + bin1_arealoss, tmp ! logical (kind=log_kind) :: & flag ! .true. if there could be lateral melting @@ -965,7 +968,6 @@ subroutine lateral_melt (dt, ncat, & vicen_init, & ! volume per unit area of ice (m) G_radialn , & ! rate of lateral melt (m/s) delta_an , & ! change in the ITD - qin , & ! enthalpy rsiden ! delta_an/aicen real (kind=dbl_kind), dimension (nfsd,ncat) :: & @@ -979,11 +981,9 @@ subroutine lateral_melt (dt, ncat, & real (kind=dbl_kind), dimension(nfsd+1) :: & f_flx ! -!echmod - for average qin - real (kind=dbl_kind), intent(in) :: & - sss real (kind=dbl_kind) :: & - Ti, Si0, qi0, sicen, & + sicen, & + etot, & ! column energy per itd cat, for FSD code elapsed_t, & ! FSD subcycling subdt ! FSD timestep (s) @@ -996,12 +996,11 @@ subroutine lateral_melt (dt, ncat, & dfsalt = c0 dvssl = c0 dvint = c0 - cat1_arealoss = c0 + bin1_arealoss = c0 tmp = c0 vicen_init = c0 G_radialn = c0 delta_an = c0 - qin = c0 rsiden = c0 afsdn = c1 afsdn_init = c0 @@ -1018,59 +1017,57 @@ subroutine lateral_melt (dt, ncat, & d_afsd_latm(:) = c0 end if - if (tr_fsd .and. fside < c0) then + if (tr_fsd .and. wlat > puny) then flag = .true. - - -! echmod - using category values would be preferable to the average value - ! Compute average enthalpy of ice (taken from add_new_ice) - if (sss > c2 * dSin0_frazil) then - Si0 = sss - dSin0_frazil - else - Si0 = sss**2 / (c4*dSin0_frazil) - endif - Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1) - qi0 = enthalpy_mush(Ti, Si0) - + ! for FSD rside and fside not yet computed correctly, redo here + fside = c0 do n = 1, ncat - if (ktherm == 2) then ! mushy - do k = 1, nilyr - !qin(n) = qin(n) & - ! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind) - qin(n) = qi0 - enddo - else - qin(n) = -rhoi*Lfresh - endif - if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative + G_radialn(n) = -wlat ! negative - if (G_radialn(n) < -puny) then + if (any(afsdn(:,n) < c0)) then + write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n + call icepack_warnings_add(warnstr) + endif + bin1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt & + * G_radialn(n) / floe_binwidth(1) - if (any(afsdn(:,n) < c0)) print*,& - 'lateral_melt B afsd < 0',n + delta_an(n) = c0 + do k = 1, nfsd + delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) & + * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0 + end do - cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt & - * G_radialn(n) / floe_binwidth(1) + ! add negative area loss from fsd + delta_an(n) = delta_an(n) - bin1_arealoss - delta_an(n) = c0 - do k = 1, nfsd - delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) & - * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0 - end do + if (delta_an(n) > c0) then + write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n) + call icepack_warnings_add(warnstr) + endif - ! add negative area loss from fsd - delta_an(n) = delta_an(n) - cat1_arealoss + ! following original code, not necessary for fsd + if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1) - if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n) + if (rsiden(n) < c0) then + write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n) + call icepack_warnings_add(warnstr) + endif - ! following original code, not necessary for fsd - if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1) + ! melting energy/unit area in each column, etot < 0 + etot = c0 + do k = 1, nslyr + etot = etot + trcrn(nt_qsno+k-1,n) * vsnon(n)/real(nslyr,kind=dbl_kind) + enddo - if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n) + do k = 1, nilyr + etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind) + enddo ! nilyr + + ! lateral heat flux, fside < 0 + fside = fside + rsiden(n)*etot/dt - end if ! G_radialn enddo ! ncat else if (rside > c0) then ! original, non-fsd implementation @@ -1132,8 +1129,10 @@ subroutine lateral_melt (dt, ncat, & DO WHILE (elapsed_t.lt.dt) nsubt = nsubt + 1 - if (nsubt.gt.100) & - print *, 'latm not converging' + if (nsubt.gt.100) then + write(warnstr,*) subname, 'latm not converging' + call icepack_warnings_add(warnstr) + endif ! finite differences df_flx(:) = c0 @@ -1146,8 +1145,10 @@ subroutine lateral_melt (dt, ncat, & df_flx(k) = f_flx(k+1) - f_flx(k) end do - if (abs(sum(df_flx(:))) > puny) & - print*,'sum(df_flx)/=0' + if (abs(sum(df_flx(:))) > puny) then + write(warnstr,*) subname, 'sum(df_flx) /= 0' + call icepack_warnings_add(warnstr) + endif ! this term ensures area conservation tmp = SUM(afsd_tmp(:)/floe_rad_c(:)) @@ -1243,6 +1244,8 @@ subroutine lateral_melt (dt, ncat, & if (tr_fsd) then + trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn + call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) ) if (icepack_warnings_aborted(subname)) return @@ -1970,7 +1973,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & Tf, sss, & salinz, & rside, meltl, & - fside, & + fside, wlat, & frzmlt, frazil, & frain, fpond, & fresh, fsalt, & @@ -2010,10 +2013,12 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) rside , & ! fraction of ice that melts laterally - fside , & ! lateral heat flux (W/m^2) frzmlt , & ! freezing/melting potential (W/m^2) wave_sig_ht ! significant height of waves in ice (m) + real (kind=dbl_kind), intent(in), optional :: & + wlat ! lateral melt rate (m/s) + real (kind=dbl_kind), dimension(:), intent(in) :: & wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s) @@ -2052,6 +2057,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration aice0 , & ! concentration of open water + fside , & ! lateral heat flux (W/m^2) frain , & ! rainfall rate (kg/m^2 s) fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) @@ -2122,6 +2128,13 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & return endif endif + if (tr_fsd) then + if (.not.(present(wlat))) then + call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif endif !----------------------------------------------------------------- @@ -2221,7 +2234,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & fhocn, faero_ocn, & fiso_ocn, & rside, meltl, & - fside, sss, & + fside, wlat, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index ac7cf232a..32a9e8f4f 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -493,7 +493,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & - fside) + fside, wlat) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -527,6 +527,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind), intent(out) :: & Tbot , & ! ice bottom surface temperature (deg C) fbot , & ! heat flux to ice bottom (W/m^2) + wlat , & ! lateral melt rate (m/s) rside , & ! fraction of ice that melts laterally fside ! lateral heat flux (W/m^2) @@ -543,7 +544,6 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind) :: & deltaT , & ! SST - Tbot >= 0 ustar , & ! skin friction velocity for fbot (m/s) - wlat , & ! lateral melt rate (m/s) xtmp ! temporary variable ! Parameters for bottom melting @@ -616,31 +616,24 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & do n = 1, ncat etot = c0 - qavg = c0 - ! melting energy/unit area in each column, etot < 0 do k = 1, nslyr etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind) - qavg = qavg + qsnon(k,n) enddo do k = 1, nilyr etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind) - qavg = qavg + qicen(k,n) enddo ! nilyr ! lateral heat flux, fside < 0 - if (tr_fsd) then ! floe size distribution - fside = fside + wlat*qavg - else ! default floe size - fside = fside + rside*etot/dt - endif + fside = fside + rside*etot/dt enddo ! n !----------------------------------------------------------------- ! Limit bottom and lateral heat fluxes if necessary. + ! FYI: fside is not yet correct for fsd, may need to adjust fbot further !----------------------------------------------------------------- xtmp = frzmlt/(fbot + fside + puny) @@ -2119,7 +2112,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fbot , & Tbot , Tsnice , & frzmlt , rside , & - fside , & + fside , wlat , & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & @@ -2260,6 +2253,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & mlt_onset , & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) + real (kind=dbl_kind), intent(out), optional :: & + wlat ! lateral melt rate (m/s) + real (kind=dbl_kind), intent(inout), optional :: & fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2) @@ -2449,6 +2445,13 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & call icepack_warnings_setabort(.true.,__FILE__,__LINE__) return endif + if (tr_fsd) then + if (.not.(present(wlat))) then + call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif endif !----------------------------------------------------------------- @@ -2543,7 +2546,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & - fside) + fside, wlat) if (icepack_warnings_aborted(subname)) return diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 7e0e992ee..749e52396 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -276,6 +276,7 @@ module icedrv_flux real (kind=dbl_kind), dimension (nx), public :: & rside , & ! fraction of ice that melts laterally fside , & ! lateral heat flux (W/m^2) + wlat , & ! lateral melt rate (m/s) fsw , & ! incoming shortwave radiation (W/m^2) coszen , & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv, & ! convergence term for ridging (1/s) diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index bc6279d22..0b5e74dda 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -111,7 +111,7 @@ subroutine step_therm1 (dt) use icedrv_arrays_column, only: meltsliqn, meltsliq use icedrv_calendar, only: yday use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx - use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, & + use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, wlat, & fbot, Tbot, Tsnice use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn @@ -324,6 +324,7 @@ subroutine step_therm1 (dt) fbot = fbot(i), frzmlt = frzmlt(i), & Tbot = Tbot(i), Tsnice = Tsnice(i), & rside = rside(i), fside = fside(i), & + wlat = wlat(i), & fsnow = fsnow(i), frain = frain(i), & fpond = fpond(i), fsloss = fsloss(i), & fsurf = fsurf(i), fsurfn = fsurfn(i,:), & @@ -433,7 +434,7 @@ subroutine step_therm2 (dt) use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, & nltrcr, nx, nfsd use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset - use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside + use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn use icedrv_init, only: tmask @@ -496,6 +497,7 @@ subroutine step_therm2 (dt) nt_strata=nt_strata(1:ntrcr,:), & Tf=Tf(i), sss=sss(i), & salinz=salinz(i,:), fside=fside(i), & + wlat=wlat(i), & rside=rside(i), meltl=meltl(i), & frzmlt=frzmlt(i), frazil=frazil(i), & frain=frain(i), fpond=fpond(i), & From 1eae173fafa9b2995ba4b3952e6592f2ade0a237 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Wed, 15 Mar 2023 14:25:31 -0700 Subject: [PATCH 4/4] Fix readthedocs errors, update documentation (#432) * update documentation to fix latex (pdf) errors --- doc/source/science_guide/sg_snow.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/science_guide/sg_snow.rst b/doc/source/science_guide/sg_snow.rst index 7b0ea3dce..3b8db8a77 100755 --- a/doc/source/science_guide/sg_snow.rst +++ b/doc/source/science_guide/sg_snow.rst @@ -73,8 +73,8 @@ where :math:`\rho_s` and :math:`\rho_{max}` are the effective snow density and t :math:`\Phi_E \Delta t` represents the maximum snow mass per unit area that may be suspended from each category, subject to the total mass (per unit area) available on each category. -Erosion begins when the instantaneous wind speed :math:`V` exceeds the seasonal wind speed required to compact the snow to a density :math:`\rho_s`, :math:`V^* = (\rho_s − \beta)/\alpha`. :math:`\sigma_{ITD}` is the standard deviation of the ice thicknesses from the thickness distribution :math:`g` within the grid cell. :math:`\gamma` is a tuning coefficient for -the eroded mass, which :cite:`Lecomte15` set to :math:`10^{-5}` kg m :math:`^{-2}`. From :cite:`Lecomte13`, :math:`\rho_s = 44.6V^* + 174` kg m :math:`^{−3}` for seasonal mean wind speed :math:`V`, i.e. :math:`\alpha=174` kg m :math:`^{-3}` and :math:`\beta=44.6` kg s m :math:`^{-4}`. +Erosion begins when the instantaneous wind speed :math:`V` exceeds the seasonal wind speed required to compact the snow to a density :math:`\rho_s`, :math:`V^* = (\rho_s - \beta)/\alpha`. :math:`\sigma_{ITD}` is the standard deviation of the ice thicknesses from the thickness distribution :math:`g` within the grid cell. :math:`\gamma` is a tuning coefficient for +the eroded mass, which :cite:`Lecomte15` set to :math:`10^{-5}` kg m :math:`^{-2}`. From :cite:`Lecomte13`, :math:`\rho_s = 44.6V^* + 174` kg m :math:`^{-3}` for seasonal mean wind speed :math:`V`, i.e. :math:`\alpha=174` kg m :math:`^{-3}` and :math:`\beta=44.6` kg s m :math:`^{-4}`. In :cite:`Lecomte15`, the fraction of this suspended snow lost in leads is