diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index fec4c6823..df50dd99e 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -85,8 +85,14 @@ module ice_dyn_shared logical (kind=log_kind), public :: & basalstress ! if true, basal stress for landfast on + ! basal stress parameters real (kind=dbl_kind), public :: & - k1 ! 1st free parameter for landfast parameterization + k1, & ! 1st free parameter for landfast parameterization + k2, & ! second free parameter (N/m^3) for landfast parametrization + alphab, & ! alphab=Cb factor in Lemieux et al 2015 + threshold_hw ! max water depth for grounding + ! see keel data from Amundrud et al. 2004 (JGR) + !======================================================================= @@ -887,11 +893,7 @@ subroutine basal_stress_coeff (nx_block, ny_block, & au, & ! concentration of ice at u location hu, & ! volume per unit area of ice at u location (mean thickness) hwu, & ! water depth at u location - hcu, & ! critical thickness at u location - k2 = 15.0_dbl_kind , & ! second free parameter (N/m^3) for landfast parametrization - alphab = 20.0_dbl_kind, & ! alphab=Cb factor in Lemieux et al 2015 - threshold_hw = 30.0_dbl_kind ! max water depth for grounding - ! see keel data from Amundrud et al. 2004 (JGR) + hcu ! critical thickness at u location integer (kind=int_kind) :: & i, j, ij diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index 8e3e1476a..ee7d7301c 100755 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -117,7 +117,7 @@ module ice_forcing ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf atm_data_type, & ! 'default', 'monthly', 'ncar', ! 'LYq' or 'hadgem' or 'oned' or - ! 'JRA55' + ! 'JRA55_gx1' or 'JRA55_gx3' bgc_data_type, & ! 'default', 'clim' ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned', ! 'hadgem_sst' or 'hadgem_sst_uvocn' @@ -241,8 +241,8 @@ subroutine init_forcing_atmo file=__FILE__, line=__LINE__) endif - if (use_leap_years .and. (trim(atm_data_type) /= 'JRA55' .and. & - trim(atm_data_type) /= 'default' .and. & + if (use_leap_years .and. (trim(atm_data_type) /= 'JRA55_gx1' .and. & + trim(atm_data_type) /= 'JRA55_gx3' .and. & trim(atm_data_type) /= 'hycom' .and. & trim(atm_data_type) /= 'box2001')) then write(nu_diag,*) 'use_leap_years option is currently only supported for' @@ -259,8 +259,10 @@ subroutine init_forcing_atmo call NCAR_files(fyear) elseif (trim(atm_data_type) == 'LYq') then call LY_files(fyear) - elseif (trim(atm_data_type) == 'JRA55') then - call JRA55_files(fyear) + elseif (trim(atm_data_type) == 'JRA55_gx1') then + call JRA55_gx1_files(fyear) + elseif (trim(atm_data_type) == 'JRA55_gx3') then + call JRA55_gx3_files(fyear) elseif (trim(atm_data_type) == 'hadgem') then call hadgem_files(fyear) elseif (trim(atm_data_type) == 'monthly') then @@ -556,7 +558,9 @@ subroutine get_forcing_atmo call ncar_data elseif (trim(atm_data_type) == 'LYq') then call LY_data - elseif (trim(atm_data_type) == 'JRA55') then + elseif (trim(atm_data_type) == 'JRA55_gx1') then + call JRA55_data(fyear) + elseif (trim(atm_data_type) == 'JRA55_gx3') then call JRA55_data(fyear) elseif (trim(atm_data_type) == 'hadgem') then call hadgem_data @@ -1395,7 +1399,11 @@ subroutine file_year (data_file, yr) i = index(data_file,'.nc') - 5 tmpname = data_file write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' - elseif (trim(atm_data_type) == 'JRA55') then ! netcdf + elseif (trim(atm_data_type) == 'JRA55_gx1') then ! netcdf + i = index(data_file,'.nc') - 5 + tmpname = data_file + write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' + elseif (trim(atm_data_type) == 'JRA55_gx3') then ! netcdf i = index(data_file,'.nc') - 5 tmpname = data_file write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' @@ -2025,12 +2033,12 @@ subroutine LY_files (yr) endif ! master_task end subroutine LY_files - subroutine JRA55_files(yr) + subroutine JRA55_gx1_files(yr) ! integer (kind=int_kind), intent(in) :: & yr ! current forcing year - character(len=*), parameter :: subname = '(JRA55_files)' + character(len=*), parameter :: subname = '(JRA55_gx1_files)' uwind_file = & trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc' @@ -2040,8 +2048,23 @@ subroutine JRA55_files(yr) write (nu_diag,*) 'Atmospheric data files:' write (nu_diag,*) trim(uwind_file) endif - end subroutine JRA55_files + end subroutine JRA55_gx1_files + subroutine JRA55_gx3_files(yr) +! + integer (kind=int_kind), intent(in) :: & + yr ! current forcing year + + character(len=*), parameter :: subname = '(JRA55_gx3_files)' + uwind_file = & + trim(atm_data_dir)//'/8XDAILY/JRA55_gx3_03hr_forcing_2005.nc' + call file_year(uwind_file,yr) + if (my_task == master_task) then + write (nu_diag,*) ' ' + write (nu_diag,*) 'Atmospheric data files:' + write (nu_diag,*) trim(uwind_file) + endif + end subroutine JRA55_gx3_files !======================================================================= ! ! read Large and Yeager atmospheric data diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index c0cc1d6fe..41ff70aec 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -95,7 +95,8 @@ subroutine input_data dxrect, dyrect use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, & kevp_kernel, & - basalstress, k1, Ktens, e_ratio, coriolis, & + basalstress, k1, k2, alphab, threshold_hw, & + Ktens, e_ratio, coriolis, & kridge, ktransport, brlx, arlx use ice_transport_driver, only: advection use ice_restoring, only: restore_ice @@ -184,7 +185,7 @@ subroutine input_data advection, coriolis, kridge, ktransport, & kstrength, krdg_partic, krdg_redist, mu_rdg, & e_ratio, Ktens, Cf, basalstress, & - k1 + k1, k2, alphab, threshold_hw namelist /shortwave_nml/ & shortwave, albedo_type, & @@ -296,6 +297,9 @@ subroutine input_data close_boundaries = .false. ! true = set land on edges of grid basalstress= .false. ! if true, basal stress for landfast is on k1 = 8.0_dbl_kind ! 1st free parameter for landfast parameterization + k2 = 15.0_dbl_kind ! dah: second free parameter (N/m^3) for landfast parametrization + alphab = 20.0_dbl_kind ! alphab=Cb factor in Lemieux et al 2015 + threshold_hw = 30.0_dbl_kind ! max water depth for grounding Ktens = 0.0_dbl_kind ! T=Ktens*P (tensile strength: see Konig and Holland, 2010) e_ratio = 2.0_dbl_kind ! EVP ellipse aspect ratio advection = 'remap' ! incremental remapping transport scheme @@ -569,6 +573,9 @@ subroutine input_data call broadcast_scalar(Cf, master_task) call broadcast_scalar(basalstress, master_task) call broadcast_scalar(k1, master_task) + call broadcast_scalar(k2, master_task) + call broadcast_scalar(alphab, master_task) + call broadcast_scalar(threshold_hw, master_task) call broadcast_scalar(Ktens, master_task) call broadcast_scalar(e_ratio, master_task) call broadcast_scalar(advection, master_task) @@ -1036,8 +1043,12 @@ subroutine input_data write(nu_diag,1000) ' mu_rdg = ', mu_rdg if (kstrength == 1) & write(nu_diag,1000) ' Cf = ', Cf + write(nu_diag,1010) ' basalstress = ', basalstress write(nu_diag,1005) ' k1 = ', k1 + write(nu_diag,1005) ' k2 = ', k2 + write(nu_diag,1005) ' alphab = ', alphab + write(nu_diag,1005) ' threshold_hw = ', threshold_hw write(nu_diag,1005) ' Ktens = ', Ktens write(nu_diag,1005) ' e_ratio = ', e_ratio write(nu_diag,1030) ' advection = ', & diff --git a/cicecore/cicedynB/infrastructure/ice_read_write.F90 b/cicecore/cicedynB/infrastructure/ice_read_write.F90 index 8eef90431..f497db49b 100644 --- a/cicecore/cicedynB/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedynB/infrastructure/ice_read_write.F90 @@ -28,6 +28,11 @@ module ice_read_write implicit none private + + integer (kind=int_kind), parameter, private :: & + bits_per_byte = 8 ! number of bits per byte. + ! used to determine RecSize in ice_open + public :: ice_open, & ice_open_ext, & ice_open_nc, & @@ -86,7 +91,7 @@ subroutine ice_open(nu, filename, nbits, algn) nbits ! no. of bits per variable (0 for sequential access) integer (kind=int_kind), intent(in), optional :: algn - integer (kind=int_kind) :: RecSize, Remnant + integer (kind=int_kind) :: RecSize, Remnant, nbytes character (*) :: filename @@ -99,7 +104,13 @@ subroutine ice_open(nu, filename, nbits, algn) open(nu,file=filename,form='unformatted') else ! direct access - RecSize = nx_global*ny_global*nbits/8 + + ! use nbytes to compute RecSize. + ! this prevents integer overflow with large global grids using nbits + ! nx*ny*nbits > 2^31 -1 (i.e., global grid 9000x7054x64) + nbytes = nbits/bits_per_byte + RecSize = nx_global*ny_global*nbytes + if (present(algn)) then ! If data is keept in blocks using given sizes (=algn) ! Used in eg. HYCOM binary files, which are stored as "blocks" dividable by 16384 bit (=algn) @@ -131,6 +142,8 @@ subroutine ice_open_ext(nu, filename, nbits) integer (kind=int_kind), intent(in) :: & nu , & ! unit number nbits ! no. of bits per variable (0 for sequential access) + + integer (kind=int_kind) :: RecSize, nbytes character (*) :: filename @@ -150,7 +163,12 @@ subroutine ice_open_ext(nu, filename, nbits) nx = nx_global + 2*nghost ny = ny_global + 2*nghost - open(nu,file=filename,recl=nx*ny*nbits/8, & + ! use nbytes to compute RecSize. + ! this prevents integer overflow with large global grids using nbits + ! nx*ny*nbits > 2^31 -1 (i.e., global grid 9000x7054x64) + nbytes = nbits/bits_per_byte + RecSize = nx*ny*nbytes + open(nu,file=filename,recl=RecSize, & form='unformatted',access='direct') endif ! nbits = 0 diff --git a/cicecore/shared/ice_init_column.F90 b/cicecore/shared/ice_init_column.F90 index 1c3dcfacd..19deb0159 100644 --- a/cicecore/shared/ice_init_column.F90 +++ b/cicecore/shared/ice_init_column.F90 @@ -254,8 +254,6 @@ subroutine init_shortwave allocate(ztrcr_sw(nbtrcr_sw, ncat)) - !!$OMP PARALLEL DO PRIVATE(iblk,i,j,n,ilo,ihi,jlo,jhi,this_block, & - !!$OMP l_print_point,debug,ipoint) do iblk=1,nblocks ! Initialize @@ -387,12 +385,18 @@ subroutine init_shortwave enddo endif + enddo ! i + enddo ! j + !----------------------------------------------------------------- ! Aggregate albedos + ! Match loop order in coupling_prep for same order of operations !----------------------------------------------------------------- - do n = 1, ncat - + do n = 1, ncat + do j = jlo, jhi + do i = ilo, ihi + if (aicen(i,j,n,iblk) > puny) then alvdf(i,j,iblk) = alvdf(i,j,iblk) & @@ -422,7 +426,12 @@ subroutine init_shortwave endif ! aicen > puny - enddo ! ncat + enddo ! i + enddo ! j + enddo ! ncat + + do j = 1, ny_block + do i = 1, nx_block !---------------------------------------------------------------- ! Store grid box mean albedos and fluxes before scaling by aice @@ -432,14 +441,14 @@ subroutine init_shortwave alidf_ai (i,j,iblk) = alidf (i,j,iblk) alvdr_ai (i,j,iblk) = alvdr (i,j,iblk) alidr_ai (i,j,iblk) = alidr (i,j,iblk) - + ! for history averaging !echmod? cszn = c0 !echmod if (coszen(i,j,iblk) > puny) cszn = c1 !echmod do n = 1, nstreams !echmod albcnt(i,j,iblk,n) = albcnt(i,j,iblk,n) + cszn !echmod enddo - + !---------------------------------------------------------------- ! Save net shortwave for scaling factor in scale_factor !---------------------------------------------------------------- diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index f321ff38b..406e8ec91 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -118,6 +118,9 @@ e_ratio = 2. basalstress = .false. k1 = 8. + k2 = 15. + alphab = 20. + threshold_hw = 30. coriolis = 'latitude' kridge = 1 ktransport = 1 diff --git a/configuration/scripts/options/set_nml.jra55 b/configuration/scripts/options/set_nml.jra55 index d0112a857..b35e99a6d 100755 --- a/configuration/scripts/options/set_nml.jra55 +++ b/configuration/scripts/options/set_nml.jra55 @@ -11,7 +11,7 @@ maskhalo_bound = .true. fyear_init = 2005 atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' atm_data_format = 'nc' -atm_data_type = 'JRA55' +atm_data_type = 'JRA55_gx1' precip_units = 'mks' ocn_data_dir = 'default' bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx1 b/configuration/scripts/options/set_nml.jra55_gx1 new file mode 100755 index 000000000..b35e99a6d --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55_gx1 @@ -0,0 +1,17 @@ +year_init = 2005 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2005 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx1' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx1_2008 b/configuration/scripts/options/set_nml.jra55_gx1_2008 new file mode 100755 index 000000000..bf92e93d1 --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55_gx1_2008 @@ -0,0 +1,17 @@ +year_init = 2008 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2008 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx1' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx3 b/configuration/scripts/options/set_nml.jra55_gx3 new file mode 100755 index 000000000..c2d1eefdc --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55_gx3 @@ -0,0 +1,17 @@ +year_init = 2005 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2005 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx3' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_gx3_2008 b/configuration/scripts/options/set_nml.jra55_gx3_2008 new file mode 100755 index 000000000..89c7dc55d --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55_gx3_2008 @@ -0,0 +1,17 @@ +year_init = 2008 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx3/iced_gx3_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2008 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx3/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55_gx3' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index b467b2fd3..4b0d35e5a 100755 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -37,8 +37,10 @@ smoke gx3 8x1 bgcskl,debug #smoke gx3 4x1 bgcz,thread smoke_gx3_8x2_bgcz restart gx1 4x2 bgcsklclim,medium restart gx1 8x1 bgczclim,medium -smoke gx1 24x1 jra55_2008,medium,run90day -restart gx1 24x1 jra55,short +smoke gx1 24x1 jra55_gx1_2008,medium,run90day +smoke gx3 8x1 jra55_gx3_2008,medium,run90day +restart gx1 24x1 jra55_gx1,short +restart gx3 8x1 jra55_gx3,short smoke gx3 4x2 fsd1,diag24,run5day,debug smoke gx3 8x2 fsd12,diag24,run5day,short restart gx3 4x2 fsd12,debug,short diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index 8117fd58d..4c9b6d502 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -163,7 +163,7 @@ The parameterization for the seabed stress is described in :cite:`Lemieux16`. Th :label: Cb where :math:`k_2` determines the maximum seabed stress that can be sustained by the grounded parameterized ridge(s), :math:`u_0` -is a small residual velocity and :math:`\alpha_b=20` is a parameter to ensure that the seabed stress quickly drops when +is a small residual velocity and :math:`\alpha_b` is a parameter to ensure that the seabed stress quickly drops when the ice concentration is smaller than 1. In the code, :math:`k_2 \max [0,(h_u - h_{cu})] e^{-\alpha_b * (1 - a_u)}` is defined as :math:`T_b`. The quantities :math:`h_u`, :math:`a_{u}` and :math:`h_{cu}` are calculated at the 'u' point based on local ice conditions (surrounding tracer points). They are respectively given by diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index f57cfd569..411b22fb8 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -275,6 +275,9 @@ Table of namelist options "","", "``-1``", "Transport Disabled", "" "","``basalstress``", "true/false", "use basal stress parameterization for landfast ice", "" "","``k1``", "real", "1st free parameter for landfast parameterization", "8." + "","``k2``", "real", "2nd free parameter (N/m^3) for landfast parameterization", "15." + "","``alphab``", "real", ":math:`\alpha_{b}` factor in :cite:`Lemieux16`", "20." + "","``threshold_hw``", "real", "Max water depth for grounding (see :cite:`Amundrud04`)", "30." "","``e_ratio``", "real", "EVP ellipse aspect ratio", "2.0" "","``Ktens``", "real", "Tensile strength factor (see :cite:`Konig10`)", "0.0" "","", "", "", "" @@ -347,7 +350,8 @@ Table of namelist options "","``nfreq``", "25", "number of frequencies in ocean surface wave spectral forcing", "" "\*","``atm_data_type``", "``default``", "constant values defined in the code", "" "","", "``LYq``", "COREII Large-Yeager (AOMIP) forcing data", ":cite:`Large09`" - "","", "``JRA55``", "JRA55 forcing data :cite:`Tsujino18`", "" + "","", "``JRA55_gx1``", "JRA55 forcing data for gx1 grid :cite:`Tsujino18`", "" + "","", "``JRA55_gx3``", "JRA55 forcing data for gx3 grid :cite:`Tsujino18`", "" "","", "``monthly``", "monthly forcing data", "" "","", "``ncar``", "NCAR bulk forcing data", "" "","", "``box2001``", "forcing data for :cite:`Hunke01` box problem", "" diff --git a/icepack b/icepack index 725955623..6a6bb3fe6 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 7259556232cebf37f4c9b42b4ba0b866a07d196b +Subproject commit 6a6bb3fe62b80727e9b31145dcaf1bf224a359ef