diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 166f4104c..736732a96 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -7,7 +7,8 @@ module icepack_parameters use icepack_kinds - use icepack_warnings, only: icepack_warnings_aborted + use icepack_warnings, only: icepack_warnings_aborted, & + icepack_warnings_add, icepack_warnings_setabort implicit none private @@ -317,7 +318,7 @@ module icepack_parameters character (len=char_len), public :: & snwredist = 'none', & ! type of snow redistribution - snw_aging_table = 'test' ! lookup table: 'snicar' or 'test' + snw_aging_table = 'test' ! lookup table: 'snicar' or 'test' or 'file' logical (kind=log_kind), public :: & use_smliq_pnd = .false. , & ! use liquid in snow for ponds @@ -340,10 +341,14 @@ module icepack_parameters isnw_rhos ! maximum snow density index ! dry snow aging parameters + real (kind=dbl_kind), dimension(:), allocatable, public :: & + snowage_rhos, & ! snowage table dimension data for rhos (kg/m^3) + snowage_Tgrd, & ! snowage table dimension data for temp gradient (deg K/m) + snowage_T ! snowage table dimension data for temperature (deg K) real (kind=dbl_kind), dimension(:,:,:), allocatable, public :: & - snowage_tau, & ! (10^-6 m) - snowage_kappa, & ! - snowage_drdt0 ! (10^-6 m/hr) + snowage_tau, & ! snowage table 3D data for tau (10^-6 m) + snowage_kappa, & ! snowage table 3D data for kappa (10^-6 m) + snowage_drdt0 ! snowage table 3D data for drdt0 (10^-6 m/hr) !----------------------------------------------------------------------- ! Parameters for biogeochemistry @@ -454,6 +459,7 @@ subroutine icepack_init_parameters( & snwredist_in, use_smliq_pnd_in, rsnw_fall_in, rsnw_tmax_in, & rhosnew_in, rhosmin_in, rhosmax_in, windmin_in, drhosdwind_in, & snwlvlfac_in, isnw_T_in, isnw_Tgrd_in, isnw_rhos_in, & + snowage_rhos_in, snowage_Tgrd_in, snowage_T_in, & snowage_tau_in, snowage_kappa_in, snowage_drdt0_in, & snw_aging_table_in) @@ -513,7 +519,7 @@ subroutine icepack_init_parameters( & ! 1 = Bitz and Lipscomb 1999 ! 2 = mushy layer theory - character (char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & conduct_in, & ! 'MU71' or 'bubbly' fbot_xfer_type_in ! transfer coefficient type for ice-ocean heat flux @@ -538,7 +544,7 @@ subroutine icepack_init_parameters( & phi_c_slow_mode_in , & ! liquid fraction porosity cutoff for slow mode phi_i_mushy_in ! liquid fraction of congelation ice - character(len=char_len), intent(in), optional :: & + character(len=*), intent(in), optional :: & tfrz_option_in ! form of ocean freezing temperature ! 'minus1p8' = -1.8 C ! 'linear_salt' = -depressT * sss @@ -561,7 +567,7 @@ subroutine icepack_init_parameters( & awtvdf_in, & ! visible, diffuse awtidf_in ! near IR, diffuse - character (len=char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & shortwave_in, & ! shortwave method, 'ccsm3' or 'dEdd' albedo_type_in ! albedo parameterization, 'ccsm3' or 'constant' ! shortwave='dEdd' overrides this parameter @@ -629,7 +635,7 @@ subroutine icepack_init_parameters( & qqqocn_in, & ! for qsat over ocn TTTocn_in ! for qsat over ocn - character (len=char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & atmbndy_in ! atmo boundary method, 'default' ('ccsm3') or 'constant' logical (kind=log_kind), intent(in), optional :: & @@ -669,14 +675,14 @@ subroutine icepack_init_parameters( & logical (kind=log_kind), intent(in), optional :: & wave_spec_in ! if true, use wave forcing - character (len=char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & wave_spec_type_in ! type of wave spectrum forcing !----------------------------------------------------------------------- ! Parameters for biogeochemistry !----------------------------------------------------------------------- - character(char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & bgc_flux_type_in ! type of ocean-ice piston velocity ! 'constant', 'Jin2006' @@ -737,7 +743,7 @@ subroutine icepack_init_parameters( & hs0_in ! snow depth for transition to bare sea ice (m) ! level-ice ponds - character (len=char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & frzpnd_in ! pond refreezing parameterization real (kind=dbl_kind), intent(in), optional :: & @@ -755,7 +761,7 @@ subroutine icepack_init_parameters( & ! Parameters for snow redistribution, metamorphosis !----------------------------------------------------------------------- - character (len=char_len), intent(in), optional :: & + character (len=*), intent(in), optional :: & snwredist_in, & ! type of snow redistribution snw_aging_table_in ! snow aging lookup table @@ -771,11 +777,18 @@ subroutine icepack_init_parameters( & rhosmax_in, & ! maximum snow density (kg/m^3) windmin_in, & ! minimum wind speed to compact snow (m/s) drhosdwind_in, & ! wind compaction factor (kg s/m^4) - snwlvlfac_in, & ! fractional increase in snow depth + snwlvlfac_in ! fractional increase in snow depth + + integer (kind=int_kind), intent(in), optional :: & isnw_T_in, & ! maxiumum temperature index isnw_Tgrd_in, & ! maxiumum temperature gradient index isnw_rhos_in ! maxiumum snow density index + real (kind=dbl_kind), dimension(:), intent(in), optional :: & + snowage_rhos_in, & ! snowage dimension data + snowage_Tgrd_in, & ! + snowage_T_in ! + real (kind=dbl_kind), dimension(:,:,:), intent(in), optional :: & snowage_tau_in, & ! (10^-6 m) snowage_kappa_in, &! @@ -911,9 +924,109 @@ subroutine icepack_init_parameters( & if (present(isnw_T_in) ) isnw_T = isnw_T_in if (present(isnw_Tgrd_in) ) isnw_Tgrd = isnw_Tgrd_in if (present(isnw_rhos_in) ) isnw_rhos = isnw_rhos_in - if (present(snowage_tau_in) ) snowage_tau = snowage_tau_in - if (present(snowage_kappa_in) ) snowage_kappa = snowage_kappa_in - if (present(snowage_drdt0_in) ) snowage_drdt0 = snowage_drdt0_in + + ! check array sizes and re/allocate if necessary + if (present(snowage_rhos_in) ) then + if (size(snowage_rhos_in) /= isnw_rhos) then + call icepack_warnings_add(subname//' incorrect size of snowage_rhos_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_rhos)) then + allocate(snowage_rhos(isnw_rhos)) + snowage_rhos = snowage_rhos_in + elseif (size(snowage_rhos) /= isnw_rhos) then + deallocate(snowage_rhos) + allocate(snowage_rhos(isnw_rhos)) + snowage_rhos = snowage_rhos_in + else + snowage_rhos = snowage_rhos_in + endif + endif + + ! check array sizes and re/allocate if necessary + if (present(snowage_Tgrd_in) ) then + if (size(snowage_Tgrd_in) /= isnw_Tgrd) then + call icepack_warnings_add(subname//' incorrect size of snowage_Tgrd_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_Tgrd)) then + allocate(snowage_Tgrd(isnw_Tgrd)) + snowage_Tgrd = snowage_Tgrd_in + elseif (size(snowage_Tgrd) /= isnw_Tgrd) then + deallocate(snowage_Tgrd) + allocate(snowage_Tgrd(isnw_Tgrd)) + snowage_Tgrd = snowage_Tgrd_in + else + snowage_Tgrd = snowage_Tgrd_in + endif + endif + + ! check array sizes and re/allocate if necessary + if (present(snowage_T_in) ) then + if (size(snowage_T_in) /= isnw_T) then + call icepack_warnings_add(subname//' incorrect size of snowage_T_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_T)) then + allocate(snowage_T(isnw_T)) + snowage_T = snowage_T_in + elseif (size(snowage_T) /= isnw_T) then + deallocate(snowage_T) + allocate(snowage_T(isnw_T)) + snowage_T = snowage_T_in + else + snowage_T = snowage_T_in + endif + endif + + ! check array sizes and re/allocate if necessary + if (present(snowage_tau_in) ) then + if (size(snowage_tau_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then + call icepack_warnings_add(subname//' incorrect size of snowage_tau_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_tau)) then + allocate(snowage_tau(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_tau = snowage_tau_in + elseif (size(snowage_tau) /= isnw_T*isnw_Tgrd*isnw_rhos) then + deallocate(snowage_tau) + allocate(snowage_tau(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_tau = snowage_tau_in + else + snowage_tau = snowage_tau_in + endif + endif + + ! check array sizes and re/allocate if necessary + if (present(snowage_kappa_in) ) then + if (size(snowage_kappa_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then + call icepack_warnings_add(subname//' incorrect size of snowage_kappa_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_kappa)) then + allocate(snowage_kappa(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_kappa = snowage_kappa_in + elseif (size(snowage_kappa) /= isnw_T*isnw_Tgrd*isnw_rhos) then + deallocate(snowage_kappa) + allocate(snowage_kappa(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_kappa = snowage_kappa_in + else + snowage_kappa = snowage_kappa_in + endif + endif + + ! check array sizes and re/allocate if necessary + if (present(snowage_drdt0_in) ) then + if (size(snowage_drdt0_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then + call icepack_warnings_add(subname//' incorrect size of snowage_drdt0_in') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + elseif (.not.allocated(snowage_drdt0)) then + allocate(snowage_drdt0(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_drdt0 = snowage_drdt0_in + elseif (size(snowage_drdt0) /= isnw_T*isnw_Tgrd*isnw_rhos) then + deallocate(snowage_drdt0) + allocate(snowage_drdt0(isnw_T,isnw_Tgrd,isnw_rhos)) + snowage_drdt0 = snowage_drdt0_in + else + snowage_drdt0 = snowage_drdt0_in + endif + endif + if (present(bgc_flux_type_in) ) bgc_flux_type = bgc_flux_type_in if (present(z_tracers_in) ) z_tracers = z_tracers_in if (present(scale_bgc_in) ) scale_bgc = scale_bgc_in @@ -1009,6 +1122,7 @@ subroutine icepack_query_parameters( & snwredist_out, use_smliq_pnd_out, rsnw_fall_out, rsnw_tmax_out, & rhosnew_out, rhosmin_out, rhosmax_out, windmin_out, drhosdwind_out, & snwlvlfac_out, isnw_T_out, isnw_Tgrd_out, isnw_rhos_out, & + snowage_rhos_out, snowage_Tgrd_out, snowage_T_out, & snowage_tau_out, snowage_kappa_out, snowage_drdt0_out, & snw_aging_table_out) @@ -1077,7 +1191,7 @@ subroutine icepack_query_parameters( & ! 1 = Bitz and Lipscomb 1999 ! 2 = mushy layer theory - character (char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & conduct_out, & ! 'MU71' or 'bubbly' fbot_xfer_type_out ! transfer coefficient type for ice-ocean heat flux @@ -1102,7 +1216,7 @@ subroutine icepack_query_parameters( & phi_c_slow_mode_out , & ! liquid fraction porosity cutoff for slow mode phi_i_mushy_out ! liquid fraction of congelation ice - character(len=char_len), intent(out), optional :: & + character(len=*), intent(out), optional :: & tfrz_option_out ! form of ocean freezing temperature ! 'minus1p8' = -1.8 C ! 'linear_salt' = -depressT * sss @@ -1125,7 +1239,7 @@ subroutine icepack_query_parameters( & awtvdf_out, & ! visible, diffuse awtidf_out ! near IR, diffuse - character (len=char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & shortwave_out, & ! shortwave method, 'ccsm3' or 'dEdd' albedo_type_out ! albedo parameterization, 'ccsm3' or 'constant' ! shortwave='dEdd' overrides this parameter @@ -1193,7 +1307,7 @@ subroutine icepack_query_parameters( & qqqocn_out, & ! for qsat over ocn TTTocn_out ! for qsat over ocn - character (len=char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & atmbndy_out ! atmo boundary method, 'default' ('ccsm3') or 'constant' logical (kind=log_kind), intent(out), optional :: & @@ -1233,15 +1347,15 @@ subroutine icepack_query_parameters( & logical (kind=log_kind), intent(out), optional :: & wave_spec_out ! if true, use wave forcing - character (len=char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & wave_spec_type_out ! type of wave spectrum forcing !----------------------------------------------------------------------- ! Parameters for biogeochemistry !----------------------------------------------------------------------- - character(char_len), intent(out), optional :: & - bgc_flux_type_out ! type of ocean-ice piston velocity + character (len=*), intent(out), optional :: & + bgc_flux_type_out ! type of ocean-ice piston velocity ! 'constant', 'Jin2006' logical (kind=log_kind), intent(out), optional :: & @@ -1301,7 +1415,7 @@ subroutine icepack_query_parameters( & hs0_out ! snow depth for transition to bare sea ice (m) ! level-ice ponds - character (len=char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & frzpnd_out ! pond refreezing parameterization real (kind=dbl_kind), intent(out), optional :: & @@ -1319,7 +1433,7 @@ subroutine icepack_query_parameters( & ! Parameters for snow redistribution, metamorphosis !----------------------------------------------------------------------- - character (len=char_len), intent(out), optional :: & + character (len=*), intent(out), optional :: & snwredist_out, & ! type of snow redistribution snw_aging_table_out ! snow aging lookup table @@ -1335,11 +1449,18 @@ subroutine icepack_query_parameters( & rhosmax_out, & ! maximum snow density (kg/m^3) windmin_out, & ! minimum wind speed to compact snow (m/s) drhosdwind_out, & ! wind compaction factor (kg s/m^4) - snwlvlfac_out, & ! fractional increase in snow depth + snwlvlfac_out ! fractional increase in snow depth + + integer (kind=int_kind), intent(out), optional :: & isnw_T_out, & ! maxiumum temperature index isnw_Tgrd_out, & ! maxiumum temperature gradient index isnw_rhos_out ! maxiumum snow density index + real (kind=dbl_kind), dimension(:), intent(out), optional :: & + snowage_rhos_out, & ! snowage dimension data + snowage_Tgrd_out, & ! + snowage_T_out ! + real (kind=dbl_kind), dimension(:,:,:), intent(out), optional :: & snowage_tau_out, & ! (10^-6 m) snowage_kappa_out, &! @@ -1515,6 +1636,9 @@ subroutine icepack_query_parameters( & if (present(isnw_T_out) ) isnw_T_out = isnw_T if (present(isnw_Tgrd_out) ) isnw_Tgrd_out = isnw_Tgrd if (present(isnw_rhos_out) ) isnw_rhos_out = isnw_rhos + if (present(snowage_rhos_out) ) snowage_rhos_out = snowage_rhos + if (present(snowage_Tgrd_out) ) snowage_Tgrd_out = snowage_Tgrd + if (present(snowage_T_out) ) snowage_T_out = snowage_T if (present(snowage_tau_out) ) snowage_tau_out = snowage_tau if (present(snowage_kappa_out) ) snowage_kappa_out= snowage_kappa if (present(snowage_drdt0_out) ) snowage_drdt0_out= snowage_drdt0 @@ -1713,9 +1837,12 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " isnw_T = ", isnw_T write(iounit,*) " isnw_Tgrd = ", isnw_Tgrd write(iounit,*) " isnw_rhos = ", isnw_rhos - write(iounit,*) " snowage_tau = ", snowage_tau - write(iounit,*) " snowage_kappa = ", snowage_kappa - write(iounit,*) " snowage_drdt0 = ", snowage_drdt0 + write(iounit,*) " snowage_rhos = ", snowage_rhos(1) + write(iounit,*) " snowage_Tgrd = ", snowage_Tgrd(1) + write(iounit,*) " snowage_T = ", snowage_T(1) + write(iounit,*) " snowage_tau = ", snowage_tau(1,1,1) + write(iounit,*) " snowage_kappa = ", snowage_kappa(1,1,1) + write(iounit,*) " snowage_drdt0 = ", snowage_drdt0(1,1,1) write(iounit,*) " bgc_flux_type = ", bgc_flux_type write(iounit,*) " z_tracers = ", z_tracers write(iounit,*) " scale_bgc = ", scale_bgc diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90 index 5eeecd6e5..7f3088e7e 100644 --- a/columnphysics/icepack_shortwave.F90 +++ b/columnphysics/icepack_shortwave.F90 @@ -1050,6 +1050,7 @@ subroutine run_dEdd(dt, ncat, & hs0, hsnlvl, & rhosnwn(:), rsnwn(:), & rsnow(:,n)) + if (icepack_warnings_aborted(subname)) return endif ! snwredist fpn = c0 ! fraction of ice covered in pond diff --git a/columnphysics/icepack_snow.F90 b/columnphysics/icepack_snow.F90 index a714d4fa3..2139dbef4 100644 --- a/columnphysics/icepack_snow.F90 +++ b/columnphysics/icepack_snow.F90 @@ -14,10 +14,12 @@ module icepack_snow use icepack_parameters, only: snwredist, rsnw_fall, rsnw_tmax, rhosnew use icepack_parameters, only: rhosmin, rhosmax, windmin, drhosdwind use icepack_parameters, only: isnw_T, isnw_Tgrd, isnw_rhos + use icepack_parameters, only: snowage_rhos, snowage_Tgrd, snowage_T use icepack_parameters, only: snowage_tau, snowage_kappa, snowage_drdt0 use icepack_parameters, only: snw_aging_table use icepack_warnings, only: icepack_warnings_add, icepack_warnings_setabort + use icepack_warnings, only: icepack_warnings_aborted implicit none private @@ -28,11 +30,18 @@ module icepack_snow S_r = 0.033_dbl_kind, & ! irreducible saturation (Anderson 1976) S_wet= 4.22e-5_dbl_kind ! (um^3/s) wet metamorphism parameters - ! shift in indices for aging test table - integer (kind=int_kind) :: & - iT_shift, & ! shift in maximum temperature index - iTgrd_shift, & ! shift in maximum temperature gradient index - irhos_shift ! shift in maximum snow density index + real (kind=dbl_kind) :: & + min_rhos, & ! snowtable axis data, assumes linear data + del_rhos, & + min_Tgrd, & + del_Tgrd, & + min_T , & + del_T + + logical (kind=log_kind) :: & + lin_rhos = .false., & ! flag to specify whether the snowtable dimensions are linear + lin_Tgrd = .false., & ! this will allow quick lookup + lin_T = .false. !======================================================================= @@ -51,33 +60,96 @@ subroutine icepack_init_snow ! local variables + integer (kind=int_kind) :: n + character (len=*),parameter :: subname='(icepack_init_snow)' !----------------------------------------------------------------- ! Snow metamorphism lookup table !----------------------------------------------------------------- + ! if snw_aging_table = 'snicar' + ! best-fit parameters are read from a table (netcdf) + ! snowage_tau, snowage_kappa, snowage_drdt0 + ! 11 temperatures from 223.15 to 273.15 K by 5.0 + ! 31 temperature gradients from 0 to 300 K/m by 10 + ! 8 snow densities from 50 to 400 kg/m3 by 50 + + ! if snw_aging_table = 'test' + ! for testing Icepack without netcdf, + ! use a subsampled, hard-coded table + ! 5 temperatures from 243.15 by 5.0 K + ! 5 temperature gradients from 0 to 40 K/m by 10 + ! 1 snow density at 300 kg/m3 + + ! if snw_aging_table = 'file' + ! all data has to be passed into icepack_parameters + + ! tcraig, should be 223.15 and 243.15 below + if (snwgrain) then if (trim(snw_aging_table) == 'snicar') then ! read netcdf file - isnw_T = 11 ! maxiumum temperature index - isnw_Tgrd = 31 ! maxiumum temperature gradient index isnw_rhos = 8 ! maxiumum snow density index + isnw_Tgrd = 31 ! maxiumum temperature gradient index + isnw_T = 11 ! maxiumum temperature index + min_rhos = 50.0_dbl_kind ! snowtable dimension data + del_rhos = 50.0_dbl_kind + lin_rhos = .true. + min_Tgrd = 0.0_dbl_kind + del_Tgrd = 10.0_dbl_kind + lin_Tgrd = .true. + min_T = 223.0_dbl_kind + del_T = 5.0_dbl_kind + lin_T = .true. allocate (snowage_tau (isnw_rhos,isnw_Tgrd,isnw_T)) allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T)) allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T)) - iT_shift = 0 ! use entire table - iTgrd_shift = 0 ! - irhos_shift = 0 ! + allocate (snowage_rhos (isnw_rhos)) + allocate (snowage_Tgrd (isnw_Tgrd)) + allocate (snowage_T (isnw_T)) + do n = 1, isnw_rhos + snowage_rhos(n) = min_rhos + real((n-1),dbl_kind)*del_rhos + enddo + do n = 1, isnw_Tgrd + snowage_Tgrd(n) = min_Tgrd + real((n-1),dbl_kind)*del_Tgrd + enddo + do n = 1, isnw_T + snowage_T(n) = min_T + real((n-1),dbl_kind)*del_T + enddo + + elseif (trim(snw_aging_table) == 'file') then + isnw_rhos = -1 ! maxiumum snow density index + isnw_Tgrd = -1 ! maxiumum temperature gradient index + isnw_T = -1 ! maxiumum temperature index + elseif (trim(snw_aging_table) == 'test') then - isnw_T = 5 ! maxiumum temperature index - isnw_Tgrd = 5 ! maxiumum temperature gradient index isnw_rhos = 1 ! maxiumum snow density index + isnw_Tgrd = 5 ! maxiumum temperature gradient index + isnw_T = 5 ! maxiumum temperature index + min_rhos = 300.0_dbl_kind ! snowtable dimension data + del_rhos = 50.0_dbl_kind + lin_rhos = .true. + min_Tgrd = 0.0_dbl_kind + del_Tgrd = 10.0_dbl_kind + lin_Tgrd = .true. + min_T = 243.0_dbl_kind + del_T = 5.0_dbl_kind + lin_T = .true. allocate (snowage_tau (isnw_rhos,isnw_Tgrd,isnw_T)) allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T)) allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T)) - iT_shift = 4 ! index adjustments for subset of table - iTgrd_shift = 0 ! - irhos_shift = 5 ! + allocate (snowage_rhos (isnw_rhos)) + allocate (snowage_Tgrd (isnw_Tgrd)) + allocate (snowage_T (isnw_T)) + do n = 1, isnw_rhos + snowage_rhos(n) = min_rhos + real((n-1),dbl_kind)*del_rhos + enddo + do n = 1, isnw_Tgrd + snowage_Tgrd(n) = min_Tgrd + real((n-1),dbl_kind)*del_Tgrd + enddo + do n = 1, isnw_T + snowage_T(n) = min_T + real((n-1),dbl_kind)*del_T + enddo ! Subset of dry snow aging parameters snowage_tau = reshape((/ & @@ -107,6 +179,10 @@ subroutine icepack_init_snow 3.85605846, 3.85668001, 3.85844559, 3.86073682, 3.863199, & 5.0861907, 5.08765668, 5.09200195, 5.09665276, 5.10079895/), & (/isnw_rhos,isnw_Tgrd,isnw_T/)) + else + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//'ERROR: snw_aging_table value') + return endif endif @@ -207,6 +283,7 @@ subroutine icepack_step_snow(dt, nilyr, & smice, smliq, & rhos_effn, rhos_eff, & rhos_cmpn, rhos_cmp) + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Redistribute snow based on wind @@ -224,6 +301,7 @@ subroutine icepack_step_snow(dt, nilyr, & fresh, fhocn, & fsloss, rhos_cmpn, & fsnow) + if (icepack_warnings_aborted(subname)) return endif vsno = c0 @@ -260,6 +338,7 @@ subroutine icepack_step_snow(dt, nilyr, & Tsfc, zTin1, & hsn, zqsn, & smice, smliq) + if (icepack_warnings_aborted(subname)) return endif end subroutine icepack_step_snow @@ -655,6 +734,7 @@ subroutine snow_redist(dt, nslyr, ncat, wind, ain, vin, vsn, zqsn, & zs1(:), zs2(:), & hslyr, hsn_new(n), & zqsn(:,n)) + if (icepack_warnings_aborted(subname)) return endif ! nslyr > 1 endif ! |dhsn| > puny endif ! ain > puny @@ -841,6 +921,7 @@ subroutine update_snow_radius (dt, ncat, nslyr, nilyr, rsnw, hin, & drsnw_dry, zqsn(:,n), Tsfc(n), & zTin(n), hsn(n), hin(n), & smice(:,n), smliq(:,n)) + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! wet metamorphism @@ -848,6 +929,7 @@ subroutine update_snow_radius (dt, ncat, nslyr, nilyr, rsnw, hin, & do k = 1,nslyr call snow_wet_metamorph (dt, drsnw_wet(k), rsnw(k,n), & smice(k,n), smliq(k,n)) + if (icepack_warnings_aborted(subname)) return rsnw(k,n) = min(rsnw_tmax, rsnw(k,n) + drsnw_dry(k) + drsnw_wet(k)) enddo @@ -870,18 +952,18 @@ end subroutine update_snow_radius subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, & Tsfc, zTin1, hsn, hin, smice, smliq) - ! Vapor redistribution: Method is to retrieve 3 best-fit parameters that - ! depend on snow temperature, temperature gradient, and density, - ! that are derived from the microphysical model described in: - ! Flanner and Zender (2006), Linking snowpack microphysics and albedo - ! evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834. - ! The parametric equation has the form: - ! dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), where: - ! r is the effective radius, - ! tau and kappa are best-fit parameters, - ! drdt_0 is the initial rate of change of effective radius, and - ! dr_fresh is the difference between the current and fresh snow states - ! (r_current - r_fresh). + ! Vapor redistribution: Method is to retrieve 3 best-fit parameters that + ! depend on snow temperature, temperature gradient, and density, + ! that are derived from the microphysical model described in: + ! Flanner and Zender (2006), Linking snowpack microphysics and albedo + ! evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834. + ! The parametric equation has the form: + ! dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), where: + ! r is the effective radius, + ! tau and kappa are best-fit parameters, + ! drdt_0 is the initial rate of change of effective radius, and + ! dr_fresh is the difference between the current and fresh snow states + ! (r_current - r_fresh). integer (kind=int_kind), intent(in) :: & nslyr, & ! number of snow layers @@ -930,10 +1012,51 @@ subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, & dzi, & ! ice layer thickness (m) dz ! dzs + dzi (m) + logical (kind=log_kind) :: & + first_call = .true. ! first call flag + + character (char_len) :: & + string ! generic string for writing messages + character (len=*),parameter :: subname='(snow_dry_metamorph)' -! Needed for variable snow density not currently modeled -! calculate density based on liquid and ice content of snow + !----------------------------------------------------------------- + ! On the first call, check that the table is setup properly + ! Check sizes of 1D and 3D data + ! Check that the 1D data is regularly spaced and set min, del, and lin values + ! for each 1D variable. This info will be used later for the table lookup. + !----------------------------------------------------------------- + + if (first_call) then + if (isnw_rhos < 1 .or. isnw_Tgrd < 1 .or. isnw_T < 1 .or. & + size(snowage_rhos) /= isnw_rhos .or. & + size(snowage_Tgrd) /= isnw_Tgrd .or. & + size(snowage_T) /= isnw_T .or. & + size(snowage_tau) /= isnw_rhos*isnw_Tgrd*isnw_T .or. & + size(snowage_kappa) /= isnw_rhos*isnw_Tgrd*isnw_T .or. & + size(snowage_drdt0) /= isnw_rhos*isnw_Tgrd*isnw_T) then + write(string,'(a,3i4)') subname//' snowtable size1 = ',isnw_rhos, isnw_Tgrd, isnw_T + call icepack_warnings_add(string) + write(string,'(a,3i4)') subname//' snowtable size2 = ',size(snowage_rhos),size(snowage_Tgrd),size(snowage_T) + call icepack_warnings_add(string) + write(string,'(a,3i9)') subname//' snowtable size3 = ',size(snowage_tau),size(snowage_kappa),size(snowage_drdt0) + call icepack_warnings_add(string) + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//'ERROR: arrays sizes error') + return + endif + call snowtable_check_dimension(snowage_rhos, min_rhos, del_rhos, lin_rhos) + if (icepack_warnings_aborted(subname)) return + call snowtable_check_dimension(snowage_Tgrd, min_Tgrd, del_Tgrd, lin_Tgrd) + if (icepack_warnings_aborted(subname)) return + call snowtable_check_dimension(snowage_T , min_T , del_T , lin_T ) + if (icepack_warnings_aborted(subname)) return + endif + + !----------------------------------------------------------------- + ! Needed for variable snow density not currently modeled + ! calculate density based on liquid and ice content of snow + !----------------------------------------------------------------- drsnw_dry(:) = c0 zTsn(:) = c0 @@ -966,43 +1089,56 @@ subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, & zdTdz(nslyr) = min(c10*isnw_Tgrd, zdTdz(nslyr)) endif - ! if snw_aging_table = 'snicar' - ! best-fit parameters are read from a table (netcdf) - ! snowage_tau, snowage_kappa, snowage_drdt0 - ! 11 temperatures from 225 to 273 K - ! 31 temperature gradients from 0 to 300 K/m - ! 8 snow densities from 0 to 350 kg/m3 + do k = 1, nslyr + zrhos(k) = smice(k) + smliq(k) - ! if snw_aging_table = 'test' - ! for testing Icepack without netcdf, - ! use a subsampled, hard-coded table - ! 5 temperatures - ! 5 temperature gradients - ! 1 snow density + !----------------------------------------------------------------- + ! best-fit table indices: + ! Will abort if 1D data is not regularly spaced (lin_* flag must be true) + ! Leave option for doing a search into non regularly spaced 1D data in the future + ! via a binary search or similar. + !----------------------------------------------------------------- - do k = 1, nslyr - zrhos(k) = smice(k) + smliq(k) - - ! best-fit table indices: - T_idx = nint(abs(zTsn(k)+ Tffresh - 223.0_dbl_kind) / 5.0_dbl_kind, kind=int_kind) - iT_shift - Tgrd_idx = nint(zdTdz(k) / 10.0_dbl_kind, kind=int_kind) - iTgrd_shift - !rhos_idx = nint(zrhos(k)-50.0_dbl_kind) / 50.0_dbl_kind, kind=int_kind) - irhos_shift ! variable density - rhos_idx = nint((rhos-50.0_dbl_kind) / 50.0_dbl_kind, kind=int_kind) - irhos_shift ! fixed density - - ! boundary check: - T_idx = min(isnw_T, max(1,T_idx+1)) - Tgrd_idx = min(isnw_Tgrd, max(1,Tgrd_idx+1)) - rhos_idx = min(isnw_rhos, max(1,rhos_idx+1)) - - bst_tau = snowage_tau (rhos_idx,Tgrd_idx,T_idx) - bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx) - bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx) - - ! change in snow effective radius, using best-fit parameters - dr_fresh = max(c0,rsnw(k)-rsnw_fall) - drsnw_dry(k) = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1/bst_kappa))& - * (dt/3600.0_dbl_kind) - enddo + if (lin_rhos) then + rhos_idx = nint( (rhos - min_rhos) / del_rhos, kind=int_kind) + 1 + else + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for rhos not supported yet') + return + endif + + if (lin_Tgrd) then + Tgrd_idx = nint( (zdTdz(k) - min_Tgrd) / del_Tgrd, kind=int_kind) + 1 + else + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for Tgrd not supported yet') + return + endif + + if (lin_T) then + T_idx = nint(abs(zTsn(k)+ Tffresh - min_T ) / del_T , kind=int_kind) + 1 + else + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + call icepack_warnings_add(subname//'ERROR: nonlinear lookup table for T not supported yet') + return + endif + + ! boundary check: + rhos_idx = min(isnw_rhos, max(1,rhos_idx)) + Tgrd_idx = min(isnw_Tgrd, max(1,Tgrd_idx)) + T_idx = min(isnw_T , max(1,T_idx )) + + bst_tau = snowage_tau (rhos_idx,Tgrd_idx,T_idx) + bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx) + bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx) + + ! change in snow effective radius, using best-fit parameters + dr_fresh = max(c0,rsnw(k)-rsnw_fall) + drsnw_dry(k) = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1/bst_kappa))& + * (dt/3600.0_dbl_kind) + enddo + + first_call = .false. end subroutine snow_dry_metamorph @@ -1047,6 +1183,50 @@ end subroutine snow_wet_metamorph !======================================================================= +! Analyze 1D array for regular spacing for snow table lookup +! and set the min, del, and lin values. +! Tolerance for regular spacing set at 1.0e-8 * typical array value + + subroutine snowtable_check_dimension(array, min_fld, del_fld, lin_fld) + + real (kind=dbl_kind), intent(in), dimension(:) :: & + array ! array to check + + real (kind=dbl_kind), intent(inout) :: & + min_fld, & ! min value if linear + del_fld ! delta value if linear + + logical (kind=log_kind), intent(inout) :: & + lin_fld ! linear flag + + ! local temporary variables + + integer (kind=int_kind) :: n, asize + + real (kind=dbl_kind) :: tolerance ! tolerance for linear checking + + character (len=*),parameter :: subname='(snowtable_check_dimension)' + + asize = size(array) + + if (asize == 1) then + min_fld = array(1) + del_fld = array(1) ! arbitrary + lin_fld = .true. + else + lin_fld = .true. + min_fld = array(1) + del_fld = array(2) - array(1) + tolerance = 1.0e-08_dbl_kind * max(abs(array(1)),abs(array(2))) ! relative to typical value + do n = 3,asize + if (abs(array(n) - array(n-1) - del_fld) > tolerance) lin_fld = .false. + enddo + endif + + end subroutine snowtable_check_dimension + +!======================================================================= + ! Conversions between ice mass, liquid water mass in snow subroutine drain_snow (dt, nslyr, vsnon, aicen, & diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index ac2d4d1dc..f34947579 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -485,7 +485,6 @@ subroutine two_stage_solver_snow(nilyr, nslyr, & Spond, sss, & q, dSdt, & w ) - if (icepack_warnings_aborted(subname)) return ! halt if solver failed if (icepack_warnings_aborted(subname)) return @@ -531,7 +530,6 @@ subroutine two_stage_solver_snow(nilyr, nslyr, & Spond, sss, & q, dSdt, & w ) - if (icepack_warnings_aborted(subname)) return ! halt if solver failed if (icepack_warnings_aborted(subname)) return @@ -586,7 +584,6 @@ subroutine two_stage_solver_snow(nilyr, nslyr, & q, dSdt, & w ) - if (icepack_warnings_aborted(subname)) return ! halt if solver failed if (icepack_warnings_aborted(subname)) return @@ -635,7 +632,6 @@ subroutine two_stage_solver_snow(nilyr, nslyr, & Spond, sss, & q, dSdt, & w ) - if (icepack_warnings_aborted(subname)) return ! halt if solver failed if (icepack_warnings_aborted(subname)) return @@ -798,7 +794,6 @@ subroutine two_stage_solver_nosnow(nilyr, nslyr, & Spond, sss, & q, dSdt, & w ) - if (icepack_warnings_aborted(subname)) return ! halt if solver failed if (icepack_warnings_aborted(subname)) return diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index e1c1816fd..f3a833a6a 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -411,7 +411,13 @@ either Celsius or Kelvin units). "shortwave", ":math:`\bullet` flag for shortwave parameterization ('default' or 'dEdd')", "" "sil", "silicate concentration", "mmol/m\ :math:`^3`" "sk_l", "skeletal layer thickness", "0.03 m" - "snoice", "snow-ice formation", "m" + "snowage_drdt0", "snowage table 3D data for drdt0 (10^-6 m/hr)", "" + "snowage_kappa", "snowage table 3D data for kappa (10^-6 m)", "" + "snowage_rhos", "snowage table dimension data for rhos (kg/m^3)", "" + "snowage_T", "snowage table dimension data for temperature (deg K)", "" + "snowage_tau", "snowage table 3D data for tau (10^-6 m)", "" + "snowage_Tgrd", "snowage table dimension data for temp gradient (deg K/m)", "" + "snoice"", "snow-ice formation", "m"" "snowpatch", "length scale for parameterizing nonuniform snow coverage", "0.02 m" "skl_bgc", ":math:`\bullet` biogeochemistry on/off", "" "spval", "special value (single precision)", ":math:`10^{30}`", ""