diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index e1d0cee61..634bb759a 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -1263,7 +1263,8 @@ subroutine input_data if (trim(ice_data_type) == 'default') ice_data_type = 'latsst' ! For backward compatibility - if (grid_format == 'nc') grid_format = 'pop_nc' + if (grid_format == 'nc' ) grid_format = 'pop_nc' + if (grid_format == 'nc_ext') grid_format = 'pop_nc_ext' !----------------------------------------------------------------- ! verify inputs @@ -2487,7 +2488,7 @@ subroutine input_data elseif (trim(wave_height_type) == 'none') then tmpstr2 = ' : no wave height data available, default==0' endif - write(nu_diag,1010) ' wave_height_type = ', trim(wave_height_type),trim(tmpstr2) + write(nu_diag,1030) ' wave_height_type = ', trim(wave_height_type),trim(tmpstr2) endif write(nu_diag,*) ' ' @@ -2780,6 +2781,7 @@ subroutine input_data endif ! my_task = master_task if (grid_format /= 'pop_nc' .and. & + grid_format /= 'pop_nc_ext' .and. & grid_format /= 'mom_nc' .and. & grid_format /= 'geosnc' .and. & grid_format /= 'meshnc' .and. & diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index f185da3c5..205f2150b 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -783,7 +783,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !----------------------------------------------------------------------- if (.not. ltripoleOnly) then - ! tripoleOnly skip fill, do not overwrite any values in interior as they may + ! tripoleOnly skip fill, do not overwrite any values in interior as they may ! already be set and filling tripole is not necessary ! fill outer boundary as needed diff --git a/cicecore/cicedyn/infrastructure/ice_blocks.F90 b/cicecore/cicedyn/infrastructure/ice_blocks.F90 index 513f3c06f..245f77bc7 100644 --- a/cicecore/cicedyn/infrastructure/ice_blocks.F90 +++ b/cicecore/cicedyn/infrastructure/ice_blocks.F90 @@ -161,7 +161,7 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & ! the global index will go from -nghost+1:0 on the lower boundary ! and n*_global+1:n*_global+nghost on the upper boundary and the ! haloUpdate and scatter, for instance, will not fill those values -! in those cases. Other boundary condition methods will fill the +! in those cases. Other boundary condition methods will fill the ! outer halo values in cases where ice exists on those boundaries. ! !---------------------------------------------------------------------- diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 406cff6d2..13a1aa098 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -36,7 +36,7 @@ module ice_grid ew_boundary_type, ns_boundary_type, init_domain_distribution use ice_fileunits, only: nu_diag, nu_grid, nu_kmt, & get_fileunit, release_fileunit, flush_fileunit - use ice_gather_scatter, only: gather_global, scatter_global, gather_global_ext + use ice_gather_scatter, only: gather_global, scatter_global, gather_global_ext, scatter_global_ext use ice_read_write, only: ice_read, ice_read_nc, ice_read_global, & ice_read_global_nc, ice_open, ice_open_nc, ice_close_nc, ice_check_nc use ice_timers, only: timer_bound, ice_timer_start, ice_timer_stop @@ -366,7 +366,7 @@ subroutine init_grid1 fieldname ! field name in netCDF file real (kind=dbl_kind), dimension(:,:), allocatable :: & - work_g1, work_g2, work_mom + work_g1, work_g2, work_g1x, work_mom integer (kind=int_kind) :: & max_blocks_min, & ! min value of max_blocks across procs @@ -453,6 +453,17 @@ subroutine init_grid1 call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.) call ice_close_nc(fid_grid) + case('pop_nc_ext') + + fieldname='ulat' + call ice_open_nc(grid_file,fid_grid) + allocate(work_g1x(nx_global+2*nghost, ny_global+2*nghost), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,.true.) + work_g1(1:nx_global,1:ny_global) = work_g1x(1+nghost:nx_global+nghost,1+nghost:ny_global+nghost) + deallocate(work_g1x) + call ice_close_nc(fid_grid) + case default call ice_open(nu_grid,grid_file,64) @@ -468,7 +479,7 @@ subroutine init_grid1 ! Fill kmt if (trim(kmt_type) =='file') then select case(trim(grid_format)) - case ('mom_nc', 'pop_nc', 'geosnc') + case ('mom_nc', 'pop_nc', 'pop_nc_ext', 'geosnc') ! mask variable name might be kmt or mask, check both call ice_open_nc(kmt_file,fid_kmt) @@ -487,7 +498,15 @@ subroutine init_grid1 #endif call broadcast_scalar(mask_fieldname, master_task) - call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g2,.true.) + if (trim(grid_format) == 'pop_nc_ext') then + allocate(work_g1x(nx_global+2*nghost, ny_global+2*nghost), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g1x,.true.) + work_g2(1:nx_global,1:ny_global) = work_g1x(1+nghost:nx_global+nghost,1+nghost:ny_global+nghost) + deallocate(work_g1x) + else + call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g2,.true.) + endif call ice_close_nc(fid_kmt) case default @@ -593,6 +612,8 @@ subroutine init_grid2 call mom_grid ! derive cice grid from MOM supergrid nc file case ('pop_nc') call popgrid_nc ! read POP grid lengths from nc file + case ('pop_nc_ext') + call popgrid_nc_ext ! read POP extended grid lengths from nc file case ('geosnc') call geosgrid_nc ! read GEOS MOM grid used from nc file case default @@ -614,6 +635,8 @@ subroutine init_grid2 select case (trim(grid_format)) case('mom_nc', 'pop_nc' ,'geosnc') call kmtmask('nc') + case('pop_nc_ext') + call kmtmask('nc_ext') case default call kmtmask('bin') end select @@ -922,8 +945,14 @@ subroutine kmtmask(filetype) elseif (filetype == 'nc') then call ice_open_nc(kmt_file,fid_kmt) call ice_read_nc(fid_kmt,1,mask_fieldname,kmt,diag, & - field_loc=field_loc_center, & - field_type=field_type_scalar) + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_close_nc(fid_kmt) + elseif (filetype == 'nc_ext') then + call ice_open_nc(kmt_file,fid_kmt) + call ice_read_nc(fid_kmt,1,mask_fieldname,kmt,diag,restart_ext=.true., & + field_loc=field_loc_center, & + field_type=field_type_scalar) call ice_close_nc(fid_kmt) else call abort_ice(subname//' ERROR: invalid filetype='//trim(filetype), file=__FILE__, line=__LINE__) @@ -936,7 +965,6 @@ subroutine kmtmask(filetype) ihi = this_block%ihi jlo = this_block%jlo jhi = this_block%jhi - do j = jlo, jhi do i = ilo, ihi ! force grid cells to land if ocn_gridcell_frac is defined @@ -944,6 +972,10 @@ subroutine kmtmask(filetype) ocn_gridcell_frac(i,j,iblk) < puny) then kmt(i,j,iblk) = c0 endif + enddo + enddo + do j = 1,ny_block + do i = 1,nx_block if (kmt(i,j,iblk) >= p5) hm(i,j,iblk) = c1 enddo enddo @@ -988,18 +1020,18 @@ subroutine popgrid if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) call ice_read_global(nu_grid,1,work_g1,'rda8',.true.) ! ULAT - call gridbox_verts(work_g1,latt_bounds) call scatter_global(ULAT, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULAT, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULAT,latt_bounds) call ice_read_global(nu_grid,2,work_g1,'rda8',.true.) ! ULON - call gridbox_verts(work_g1,lont_bounds) call scatter_global(ULON, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULON, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULON,lont_bounds) call ice_read_global(nu_grid,7,work_g1,'rda8',.true.) ! ANGLE call scatter_global(ANGLE, work_g1, master_task, distrb_info, & @@ -1090,19 +1122,19 @@ subroutine popgrid_nc fieldname='ulat' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULAT - call gridbox_verts(work_g1,latt_bounds) call scatter_global(ULAT, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULAT, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULAT,latt_bounds) fieldname='ulon' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULON - call gridbox_verts(work_g1,lont_bounds) call scatter_global(ULON, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULON, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULON,lont_bounds) fieldname='angle' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE @@ -1172,6 +1204,210 @@ subroutine popgrid_nc end subroutine popgrid_nc +!======================================================================= +! POP extended displaced pole grid and land mask. +! Grid record number, field and units are: \\ +! (1) ULAT (radians) \\ +! (2) ULON (radians) \\ +! (3) HTN (cm) \\ +! (4) HTE (cm) \\ +! (5) HUS (cm) \\ +! (6) HUW (cm) \\ +! (7) ANGLE (radians) +! +! author: T. Craig +! Revised for netcdf input: Ann Keen, Met Office, May 2007 + + subroutine popgrid_nc_ext + +#ifdef USE_NETCDF + use netcdf, only : nf90_inq_varid , nf90_inq_dimid, & + nf90_inquire_dimension, nf90_get_var, nf90_noerr +#endif + + integer (kind=int_kind) :: & + i, j, iblk, & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + fid_grid , & ! file id for netCDF grid file + ierr + + logical (kind=log_kind) :: diag + + character (char_len) :: & + fieldname ! field name in netCDF file + + real (kind=dbl_kind) :: & + pi + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + work_g1x ! temporary global extended array + + integer(kind=int_kind) :: & + varid, status + + type (block) :: & + this_block ! block information for current block + + character(len=*), parameter :: subname = '(popgrid_nc_ext)' + +#ifdef USE_NETCDF + call icepack_query_parameters(pi_out=pi) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call ice_open_nc(grid_file,fid_grid) + + diag = .true. ! write diagnostic info + + !----------------------------------------------------------------- + ! lat, lon, angle + !----------------------------------------------------------------- + + allocate(work_g1x(nx_global+2*nghost,ny_global+2*nghost), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + work_g1x = c0 + + fieldname='ulat' + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) ! ULAT + call scatter_global_ext(ULAT, work_g1x, master_task, distrb_info) + call gridbox_verts(ULAT,latt_bounds) + + fieldname='ulon' + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) ! ULON + call scatter_global_ext(ULON, work_g1x, master_task, distrb_info) + call gridbox_verts(ULON,lont_bounds) + + fieldname='angle' + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) ! ANGLE + call scatter_global_ext(ANGLE, work_g1x, master_task, distrb_info) + ! fix ANGLE: roundoff error due to single precision + where (ANGLE > pi) ANGLE = pi + where (ANGLE < -pi) ANGLE = -pi + + ! if grid file includes anglet then read instead + fieldname='anglet' + if (my_task == master_task) then + status = nf90_inq_varid(fid_grid, fieldname , varid) + if (status /= nf90_noerr) then + write(nu_diag,*) subname//' CICE will calculate angleT, TLON and TLAT' + else + write(nu_diag,*) subname//' angleT, TLON and TLAT is read from grid file' + l_readCenter = .true. + endif + endif + call broadcast_scalar(l_readCenter,master_task) + if (l_readCenter) then + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) + call scatter_global_ext(ANGLET, work_g1x, master_task, distrb_info) + where (ANGLET > pi) ANGLET = pi + where (ANGLET < -pi) ANGLET = -pi + + fieldname="tlon" + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) + call scatter_global_ext(TLON, work_g1x, master_task, distrb_info) + + fieldname="tlat" + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) + call scatter_global_ext(TLAT, work_g1x, master_task, distrb_info) + endif + !----------------------------------------------------------------- + ! cell dimensions + ! calculate derived quantities from global arrays to preserve + ! information on boundaries + !----------------------------------------------------------------- + + fieldname='htn' + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) ! HTN + if (my_task == master_task) then + work_g1x(:,:) = work_g1x(:,:) * cm_to_m ! HTN + endif +! call primary_grid_lengths_HTN(work_g1) ! dxU, dxT, dxN, dxE + if (save_ghte_ghtn) then + if (my_task == master_task) then + G_HTN = work_g1x + endif + endif + call scatter_global_ext(HTN, work_g1x, master_task, distrb_info) + + dxN(:,:,:) = HTN(:,:,:) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + dxU(i,j,iblk) = p5*(HTN(i,j,iblk)+HTN(i+1,j,iblk)) + dxT(i,j,iblk) = p5*(HTN(i,j,iblk)+HTN(i,j-1,iblk)) + dxE(i,j,iblk) = p25*(HTN(i,j,iblk)+HTN(i+1,j,iblk)+HTN(i,j-1,iblk)+HTN(i+1,j-1,iblk)) + enddo + enddo + enddo + call ice_HaloUpdate (dxU, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxU, distrb_info, & + ew_boundary_type, ns_boundary_type) + call ice_HaloUpdate (dxT, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxT, distrb_info, & + ew_boundary_type, ns_boundary_type) + call ice_HaloUpdate (dxE, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxE, distrb_info, & + ew_boundary_type, ns_boundary_type) + + fieldname='hte' + call ice_read_global_nc(fid_grid,1,fieldname,work_g1x,diag) ! HTE + if (my_task == master_task) then + work_g1x(:,:) = work_g1x(:,:) * cm_to_m ! HTN + endif +! call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE + if (save_ghte_ghtn) then + G_HTE = work_g1x + endif + call scatter_global_ext(HTE, work_g1x, master_task, distrb_info) + dyE(:,:,:) = HTE(:,:,:) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + dyU(i,j,iblk) = p5*(HTE(i,j,iblk)+HTE(i,j+1,iblk)) + dyT(i,j,iblk) = p5*(HTE(i,j,iblk)+HTE(i-1,j,iblk)) + dyN(i,j,iblk) = p25*(HTE(i,j,iblk)+HTE(i-1,j,iblk)+HTE(i,j+1,iblk)+HTE(i-1,j+1,iblk)) + enddo + enddo + enddo + call ice_HaloUpdate (dyU, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyU, distrb_info, & + ew_boundary_type, ns_boundary_type) + call ice_HaloUpdate (dyT, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyT, distrb_info, & + ew_boundary_type, ns_boundary_type) + call ice_HaloUpdate (dyN, halo_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyN, distrb_info, & + ew_boundary_type, ns_boundary_type) + + deallocate(work_g1x, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + call ice_close_nc(fid_grid) + +#else + call abort_ice(subname//' ERROR: USE_NETCDF cpp not defined', & + file=__FILE__, line=__LINE__) +#endif + + end subroutine popgrid_nc_ext + #ifdef CESMCOUPLED !======================================================================= ! Read in kmt file that matches CAM lat-lon grid and has single column @@ -2234,19 +2470,19 @@ subroutine geosgrid_nc fieldname='ulat' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULAT - call gridbox_verts(work_g1,latt_bounds) call scatter_global(ULAT, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULAT, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULAT,latt_bounds) fieldname='ulon' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULON - call gridbox_verts(work_g1,lont_bounds) call scatter_global(ULON, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) call ice_HaloExtrapolate(ULON, distrb_info, & ew_boundary_type, ns_boundary_type) + call gridbox_verts(ULON,lont_bounds) fieldname='angle' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE @@ -2845,7 +3081,7 @@ subroutine primary_grid_lengths_HTN(work_g) enddo do j = 1, ny_global do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate ip1 = i+1 if (i == nx_global) ip1 = 1 work_g2(i,j) = p5*(work_g(i,j) + work_g(ip1,j)) ! dxU @@ -2891,7 +3127,7 @@ subroutine primary_grid_lengths_HTN(work_g) if (my_task == master_task) then do j = 2, ny_global do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate ip1 = i+1 if (i == nx_global) ip1 = 1 work_g2(i,j) = p25*(work_g(i,j)+work_g(ip1,j)+work_g(i,j-1)+work_g(ip1,j-1)) ! dxE @@ -2899,7 +3135,7 @@ subroutine primary_grid_lengths_HTN(work_g) enddo ! extrapolate to obtain dxT along j=1 do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate ip1 = i+1 if (i == nx_global) ip1 = 1 work_g2(i,1) = p5*(c2*work_g(i ,2) - work_g(i ,3) + & @@ -2983,7 +3219,7 @@ subroutine primary_grid_lengths_HTE(work_g) if (my_task == master_task) then do j = 1, ny_global do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate im1 = i-1 if (i == 1) im1 = nx_global work_g2(i,j) = p5*(work_g(i,j) + work_g(im1,j)) ! dyT @@ -3000,7 +3236,7 @@ subroutine primary_grid_lengths_HTE(work_g) if (my_task == master_task) then do j = 1, ny_global-1 do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate im1 = i-1 if (i == 1) im1 = nx_global work_g2(i,j) = p25*(work_g(i,j) + work_g(im1,j) + work_g(i,j+1) + work_g(im1,j+1)) ! dyN @@ -3009,7 +3245,7 @@ subroutine primary_grid_lengths_HTE(work_g) ! extrapolate to obtain dyN along j=ny_global if (ny_global > 1) then do i = 1, nx_global - ! assume cyclic; noncyclic will be handled during scatter + ! assume cyclic; noncyclic will be handled during extrapolate im1 = i-1 if (i == 1) im1 = nx_global work_g2(i,ny_global) = p5*(c2*work_g(i ,ny_global-1) - work_g(i ,ny_global-2) + & @@ -5267,27 +5503,25 @@ end subroutine gridbox_edges ! ! authors: A. McLaren, Met Office ! E. Hunke, LANL +! T. Craig - subroutine gridbox_verts(work_g,vbounds) + subroutine gridbox_verts(work,vbounds) - real (kind=dbl_kind), dimension(:,:), intent(in) :: & - work_g + real (kind=dbl_kind), dimension(:,:,:), intent(in) :: & + work real (kind=dbl_kind), dimension(4,nx_block,ny_block,max_blocks), intent(out) :: & - vbounds + vbounds integer (kind=int_kind) :: & - i,j , & ! index counters + iblk,i,j,ilo,ihi,jlo,jhi, & ! index counters ierr - real (kind=dbl_kind) :: & - rad_to_deg - - real (kind=dbl_kind), dimension(:,:), allocatable :: & - work_g2 + type (block) :: & + this_block ! block information for current block - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - work1 + real (kind=dbl_kind) :: & + rad_to_deg character(len=*), parameter :: subname = '(gridbox_verts)' @@ -5296,87 +5530,22 @@ subroutine gridbox_verts(work_g,vbounds) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global), stat=ierr) - else - allocate(work_g2(1,1), stat=ierr) - endif - if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) - - !------------------------------------------------------------- - ! Get coordinates of grid boxes for each block as follows: - ! (1) SW corner, (2) SE corner, (3) NE corner, (4) NW corner - !------------------------------------------------------------- - - work_g2(:,:) = c0 - if (my_task == master_task) then - do j = 2, ny_global - do i = 2, nx_global - work_g2(i,j) = work_g(i-1,j-1) * rad_to_deg - enddo - enddo - ! extrapolate - do j = 1, ny_global - work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j) - enddo - do i = 1, nx_global - work_g2(i,1) = c2*work_g2(i,2) - work_g2(i,3) - enddo - endif - call scatter_global(work1, work_g2, & - master_task, distrb_info, & - field_loc_NEcorner, field_type_scalar) - vbounds(1,:,:,:) = work1(:,:,:) - - work_g2(:,:) = c0 - if (my_task == master_task) then - do j = 2, ny_global - do i = 1, nx_global - work_g2(i,j) = work_g(i,j-1) * rad_to_deg - enddo - enddo - ! extrapolate - do i = 1, nx_global - work_g2(i,1) = (c2*work_g2(i,2) - work_g2(i,3)) - enddo - endif - call scatter_global(work1, work_g2, & - master_task, distrb_info, & - field_loc_NEcorner, field_type_scalar) - vbounds(2,:,:,:) = work1(:,:,:) - - work_g2(:,:) = c0 - if (my_task == master_task) then - do j = 1, ny_global - do i = 1, nx_global - work_g2(i,j) = work_g(i,j) * rad_to_deg - enddo - enddo - endif - call scatter_global(work1, work_g2, & - master_task, distrb_info, & - field_loc_NEcorner, field_type_scalar) - vbounds(3,:,:,:) = work1(:,:,:) - - work_g2(:,:) = c0 - if (my_task == master_task) then - do j = 1, ny_global - do i = 2, nx_global - work_g2(i,j) = work_g(i-1,j ) * rad_to_deg - enddo + vbounds = c0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + vbounds(1,i,j,iblk) = work(i-1,j-1,iblk)*rad_to_deg + vbounds(2,i,j,iblk) = work(i ,j-1,iblk)*rad_to_deg + vbounds(3,i,j,iblk) = work(i ,j ,iblk)*rad_to_deg + vbounds(4,i,j,iblk) = work(i-1,j ,iblk)*rad_to_deg enddo - ! extrapolate - do j = 1, ny_global - work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j) enddo - endif - call scatter_global(work1, work_g2, & - master_task, distrb_info, & - field_loc_NEcorner, field_type_scalar) - vbounds(4,:,:,:) = work1(:,:,:) - - deallocate (work_g2, stat=ierr) - if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + enddo end subroutine gridbox_verts diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 index 37cf4d985..0186d73f2 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 @@ -224,7 +224,7 @@ subroutine init_restart_write(filename_spec) ! write pointer (path/file) if (my_task == master_task) then -#ifdef CESMCOUPLED +#ifdef CESMCOUPLED lpointer_file = 'rpointer.ice'//trim(inst_suffix) #else lpointer_file = pointer_file diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 index 2fecaf31e..6e7ff4173 100644 --- a/cicecore/drivers/unittest/halochk/halochk.F90 +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -450,15 +450,15 @@ program halochk k1m = 1 k2m = 1 halofld = '2DL1' - where (darrayi1 == fillval) + where (darrayi1 == fillval) larrayi1 = .false. elsewhere - larrayi1 = (mod(nint(darrayi1),2) == 1) + larrayi1 = (mod(nint(darrayi1),2) == 1) endwhere - where (darrayj1 == fillval) + where (darrayj1 == fillval) larrayj1 = .true. elsewhere - larrayj1 = (mod(nint(darrayj1),2) == 1) + larrayj1 = (mod(nint(darrayj1),2) == 1) endwhere if (halofill) then call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=0) diff --git a/cicecore/drivers/unittest/opticep/ice_step_mod.F90 b/cicecore/drivers/unittest/opticep/ice_step_mod.F90 index 6d8fff7f9..6e775a3c0 100644 --- a/cicecore/drivers/unittest/opticep/ice_step_mod.F90 +++ b/cicecore/drivers/unittest/opticep/ice_step_mod.F90 @@ -566,7 +566,7 @@ subroutine step_therm1 (dt, iblk) !opt dpnd_initialn=dpnd_initialn(i,j,:,iblk), & !opt dpnd_dlidn = dpnd_dlidn (i,j,:,iblk), & yday=yday) -!opt prescribed_ice=prescribed_ice) +!opt prescribed_ice=prescribed_ice) !----------------------------------------------------------------- ! handle per-category i2x fields, no merging @@ -674,12 +674,16 @@ subroutine step_therm2 (dt, iblk) tr_fsd, & ! floe size distribution tracers z_tracers ! vertical biogeochemistry + character (len=char_len) :: & + wave_height_type ! type of significant wave height forcing + type (block) :: & this_block ! block information for current block character(len=*), parameter :: subname = '(step_therm2)' - call icepack_query_parameters(z_tracers_out=z_tracers) + call icepack_query_parameters(z_tracers_out=z_tracers, & + wave_height_type_out=wave_height_type) call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) @@ -704,57 +708,59 @@ subroutine step_therm2 (dt, iblk) if (tmask(i,j,iblk) .or. opmask(i,j,iblk)) then - ! significant wave height for FSD - if (tr_fsd) & - wave_sig_ht(i,j,iblk) = c4*SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:))) - - call icepack_step_therm2(dt=dt, & - hin_max = hin_max (:), & - aicen = aicen (i,j,:,iblk), & - vicen = vicen (i,j,:,iblk), & - vsnon = vsnon (i,j,:,iblk), & - aicen_init = aicen_init(i,j,:,iblk), & - vicen_init = vicen_init(i,j,:,iblk), & - trcrn = trcrn (i,j,:,:,iblk), & - aice0 = aice0 (i,j, iblk), & - aice = aice (i,j, iblk), & - trcr_depend= trcr_depend(:), & - trcr_base = trcr_base(:,:), & - n_trcr_strata = n_trcr_strata(:), & - nt_strata = nt_strata(:,:), & - Tf = Tf (i,j, iblk), & - sss = sss (i,j, iblk), & - salinz = salinz (i,j,:,iblk), & - rsiden = rsiden (i,j,:,iblk), & - meltl = meltl (i,j, iblk), & -!opt wlat = wlat (i,j, iblk), & - frzmlt = frzmlt (i,j, iblk), & - frazil = frazil (i,j, iblk), & - frain = frain (i,j, iblk), & - fpond = fpond (i,j, iblk), & - fresh = fresh (i,j, iblk), & - fsalt = fsalt (i,j, iblk), & - fhocn = fhocn (i,j, iblk), & - faero_ocn = faero_ocn (i,j,:,iblk), & - first_ice = first_ice (i,j,:,iblk), & - flux_bio = flux_bio (i,j,1:nbtrcr,iblk), & - ocean_bio = ocean_bio (i,j,1:nbtrcr,iblk), & - frazil_diag= frazil_diag(i,j,iblk) & -!opt frz_onset = frz_onset (i,j, iblk), & -!opt yday = yday, & -!opt fiso_ocn = fiso_ocn (i,j,:,iblk), & -!opt HDO_ocn = HDO_ocn (i,j, iblk), & -!opt H2_16O_ocn = H2_16O_ocn(i,j, iblk), & -!opt H2_18O_ocn = H2_18O_ocn(i,j, iblk), & -!opt wave_sig_ht= wave_sig_ht(i,j, iblk), & -!opt wave_spectrum=wave_spectrum(i,j,:,iblk), & -!opt wavefreq = wavefreq (:), & -!opt dwavefreq = dwavefreq (:), & -!opt d_afsd_latg=d_afsd_latg(i,j,:,iblk), & -!opt d_afsd_newi=d_afsd_newi(i,j,:,iblk), & -!opt d_afsd_latm=d_afsd_latm(i,j,:,iblk), & -!opt d_afsd_weld=d_afsd_weld(i,j,:,iblk), & -!opt dpnd_melt = dpnd_melt( i,j, iblk)) + ! significant wave height + if (tr_fsd .and. trim(wave_height_type) == 'internal') then + wave_sig_ht(i,j,iblk) = c4*SQRT(SUM(wave_spectrum(i,j,:,iblk)*dwavefreq(:))) + ! else wave_sig_ht = 0 unless provided by coupler or other external data + endif + + call icepack_step_therm2(dt=dt, & + hin_max = hin_max (:), & + aicen = aicen (i,j,:,iblk), & + vicen = vicen (i,j,:,iblk), & + vsnon = vsnon (i,j,:,iblk), & + aicen_init = aicen_init (i,j,:,iblk), & + vicen_init = vicen_init (i,j,:,iblk), & + trcrn = trcrn (i,j,:,:,iblk), & + aice0 = aice0 (i,j, iblk), & + aice = aice (i,j, iblk), & + trcr_depend = trcr_depend(:), & + trcr_base = trcr_base (:,:), & + n_trcr_strata = n_trcr_strata(:), & + nt_strata = nt_strata (:,:), & + Tf = Tf (i,j, iblk), & + sss = sss (i,j, iblk), & + salinz = salinz (i,j,:,iblk), & + rsiden = rsiden (i,j,:,iblk), & + meltl = meltl (i,j, iblk), & +!opt wlat = wlat (i,j, iblk), & + frzmlt = frzmlt (i,j, iblk), & + frazil = frazil (i,j, iblk), & + frain = frain (i,j, iblk), & + fpond = fpond (i,j, iblk), & + fresh = fresh (i,j, iblk), & + fsalt = fsalt (i,j, iblk), & + fhocn = fhocn (i,j, iblk), & + faero_ocn = faero_ocn (i,j,:,iblk), & + first_ice = first_ice (i,j,:,iblk), & + flux_bio = flux_bio (i,j,1:nbtrcr,iblk), & + ocean_bio = ocean_bio (i,j,1:nbtrcr,iblk), & + frazil_diag = frazil_diag(i,j, iblk) & +!opt frz_onset = frz_onset (i,j, iblk), & +!opt yday = yday, & +!opt fiso_ocn = fiso_ocn (i,j,:,iblk), & +!opt HDO_ocn = HDO_ocn (i,j, iblk), & +!opt H2_16O_ocn = H2_16O_ocn (i,j, iblk), & +!opt H2_18O_ocn = H2_18O_ocn (i,j, iblk), & +!opt wave_sig_ht = wave_sig_ht(i,j, iblk), & +!opt wave_spectrum = & +!opt wave_spectrum(i,j,:,iblk), & +!opt wavefreq = wavefreq (:), & +!opt d_afsd_latg = d_afsd_latg(i,j,:,iblk), & +!opt d_afsd_newi = d_afsd_newi(i,j,:,iblk), & +!opt d_afsd_latm = d_afsd_latm(i,j,:,iblk), & +!opt d_afsd_weld = d_afsd_weld(i,j,:,iblk), & +!opt dpnd_melt = dpnd_melt (i,j, iblk)) ) endif ! tmask @@ -891,7 +897,7 @@ end subroutine update_state subroutine step_dyn_wave (dt) - use ice_arrays_column, only: wave_spectrum, & + use ice_arrays_column, only: wave_spectrum, wave_sig_ht, & d_afsd_wave, wavefreq, dwavefreq use ice_domain_size, only: ncat, nfsd, nfreq use ice_state, only: trcrn, aicen, aice, vice @@ -911,14 +917,11 @@ subroutine step_dyn_wave (dt) iblk, & ! block index i, j ! horizontal indices - character (len=char_len) :: wave_spec_type - character(len=*), parameter :: subname = '(step_dyn_wave)' call ice_timer_start(timer_column) call ice_timer_start(timer_fsd) - call icepack_query_parameters(wave_spec_type_out=wave_spec_type) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -935,8 +938,7 @@ subroutine step_dyn_wave (dt) do j = jlo, jhi do i = ilo, ihi d_afsd_wave(i,j,:,iblk) = c0 - call icepack_step_wavefracture(wave_spec_type = wave_spec_type, & - dt = dt, nfreq = nfreq, & + call icepack_step_wavefracture(dt = dt, nfreq = nfreq, & aice = aice (i,j, iblk), & vice = vice (i,j, iblk), & aicen = aicen (i,j,:, iblk), & @@ -944,7 +946,8 @@ subroutine step_dyn_wave (dt) wavefreq = wavefreq (:), & dwavefreq = dwavefreq (:), & trcrn = trcrn (i,j,:,:,iblk), & - d_afsd_wave = d_afsd_wave (i,j,:, iblk)) + d_afsd_wave = d_afsd_wave (i,j,:, iblk), & + wave_height = wave_sig_ht (i,j, iblk)) end do ! i end do ! j end do ! iblk @@ -1156,8 +1159,8 @@ subroutine step_dyn_ridge (dt, ndtd, iblk) fsalt = fsalt (i,j, iblk), & first_ice = first_ice(i,j,:,iblk), & flux_bio = flux_bio (i,j,1:nbtrcr,iblk), & - Tf = Tf (i,j, iblk) & -!opt dpnd_ridge=dpnd_ridge(i,j,iblk)) + Tf = Tf (i,j, iblk) & +!opt dpnd_ridge=dpnd_ridge(i,j, iblk)) ) endif ! tmask @@ -1446,7 +1449,7 @@ subroutine step_radiation (dt, iblk) albpndn =albpndn (i,j,: ,iblk), apeffn =apeffn (i,j,: ,iblk), & snowfracn=snowfracn(i,j,: ,iblk), & dhsn =dhsn (i,j,: ,iblk), ffracn =ffracn(i,j,:,iblk), & -!opt rsnow =rsnow (:,:), +!opt rsnow =rsnow (:,:), & l_print_point=l_print_point) endif diff --git a/configuration/scripts/machines/Macros.carpenter_gnu b/configuration/scripts/machines/Macros.carpenter_gnu index 61efa80c2..a1fe056fb 100644 --- a/configuration/scripts/machines/Macros.carpenter_gnu +++ b/configuration/scripts/machines/Macros.carpenter_gnu @@ -12,7 +12,8 @@ FFLAGS := -fconvert=big-endian -fbacktrace -ffree-line-length-none -fallow-a FFLAGS_NOOPT:= -O0 ifeq ($(ICE_BLDDEBUG), true) - FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow +# FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow + FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=zero,overflow CFLAGS += -O0 endif diff --git a/configuration/scripts/machines/Macros.carpenter_gnuimpi b/configuration/scripts/machines/Macros.carpenter_gnuimpi index ef0c5e96a..d57dfda2b 100644 --- a/configuration/scripts/machines/Macros.carpenter_gnuimpi +++ b/configuration/scripts/machines/Macros.carpenter_gnuimpi @@ -12,7 +12,8 @@ FFLAGS := -fconvert=big-endian -fbacktrace -ffree-line-length-none -fallow-a FFLAGS_NOOPT:= -O0 ifeq ($(ICE_BLDDEBUG), true) - FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow +# FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow + FFLAGS += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=zero,overflow CFLAGS += -O0 endif diff --git a/configuration/scripts/machines/env.carpenter_cray b/configuration/scripts/machines/env.carpenter_cray index d2c832d8f..0d70e7f54 100644 --- a/configuration/scripts/machines/env.carpenter_cray +++ b/configuration/scripts/machines/env.carpenter_cray @@ -13,21 +13,21 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-cray/8.4.0 +module load PrgEnv-cray/8.6.0 module unload cce -module load cce/16.0.0 +module load cce/19.0.0 module unload cray-mpich -module load cray-mpich/8.1.26 +module load cray-mpich/8.1.32 module unload cray-hdf5 module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.3 -module load cray-hdf5/1.12.2.3 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -41,7 +41,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME cray -setenv ICE_MACHINE_ENVINFO "Cray Fortran/Clang 16.0.0, cray-mpich/8.1.26, netcdf/4.9.0.3" +setenv ICE_MACHINE_ENVINFO "Cray Fortran/Clang 19.0.0, cray-mpich/8.1.32, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/configuration/scripts/machines/env.carpenter_gnu b/configuration/scripts/machines/env.carpenter_gnu index 96a04072f..9c8e4df4e 100644 --- a/configuration/scripts/machines/env.carpenter_gnu +++ b/configuration/scripts/machines/env.carpenter_gnu @@ -13,15 +13,15 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-gnu/8.4.0 +module load PrgEnv-gnu/8.6.0 module unload gcc -module load gcc/12.2.0 +module load gcc/14.2.0 module unload cray-mpich module unload mpi module unload openmpi -module load cray-mpich/8.1.26 +module load cray-mpich/8.1.32 #module load openmpi/4.1.6 #module load mpi/2021.11 @@ -30,8 +30,8 @@ module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.3 -module load cray-hdf5/1.12.2.3 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -45,7 +45,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME gnu -setenv ICE_MACHINE_ENVINFO "gnu gcc 12.2.0 20220819, mpich/8.1.26, netcdf/4.9.0.3" +setenv ICE_MACHINE_ENVINFO "gnu gcc 14.3.0, mpich/8.1.32, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/configuration/scripts/machines/env.carpenter_gnuimpi b/configuration/scripts/machines/env.carpenter_gnuimpi index f21bf97a5..900494b3b 100644 --- a/configuration/scripts/machines/env.carpenter_gnuimpi +++ b/configuration/scripts/machines/env.carpenter_gnuimpi @@ -13,25 +13,25 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-gnu/8.4.0 +module load PrgEnv-gnu/8.6.0 module unload gcc -module load gcc/11.2.0 +module load gcc/14.2.0 module unload cray-mpich module unload mpi module unload openmpi #module load cray-mpich/8.1.26 #module load openmpi/4.1.6 -module load mpi/2021.11 +module load mpi/2021.16 module unload cray-hdf5 module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.3 -module load cray-hdf5/1.12.2.3 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -45,7 +45,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME gnuimpi -setenv ICE_MACHINE_ENVINFO "gnu gcc 11.2.0 20210728, intel mpi 2021.11, netcdf/4.9.0.3" +setenv ICE_MACHINE_ENVINFO "gnu gcc 14.2.0, intel mpi 2021.16, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/configuration/scripts/machines/env.carpenter_intel b/configuration/scripts/machines/env.carpenter_intel index 141a5a82a..577d9642f 100644 --- a/configuration/scripts/machines/env.carpenter_intel +++ b/configuration/scripts/machines/env.carpenter_intel @@ -12,15 +12,15 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-intel/8.5.0 +module load PrgEnv-intel/8.6.0 module unload intel -module load intel/2023.0.0 +module load intel/2024.2 module unload cray-mpich module unload mpi module unload openmpi -module load cray-mpich/8.1.30 +module load cray-mpich/8.1.32 #module load mpi/2021.11 #module load openmpi/4.1.6 @@ -29,8 +29,8 @@ module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.3 -module load cray-hdf5/1.14.3.1 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -44,7 +44,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME intel -setenv ICE_MACHINE_ENVINFO "ifort 2021.8.0 20221119, cray-mpich/8.1.30, netcdf/4.9.0.3" +setenv ICE_MACHINE_ENVINFO "ifort 2021.13.0 20240602, cray-mpich/8.1.32, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/configuration/scripts/machines/env.carpenter_intelimpi b/configuration/scripts/machines/env.carpenter_intelimpi index ef43363cd..80f696772 100644 --- a/configuration/scripts/machines/env.carpenter_intelimpi +++ b/configuration/scripts/machines/env.carpenter_intelimpi @@ -12,16 +12,16 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-intel/8.5.0 +module load PrgEnv-intel/8.6.0 module unload intel -module load intel/2023.0.0 +module load intel/2023.2.0 module unload cray-mpich module unload mpi module unload openmpi #module load cray-mpich/8.1.26 -module load mpi/2021.13 +module load mpi/2021.16 #module load openmpi/4.1.6 module unload cray-hdf5 @@ -29,8 +29,8 @@ module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.3 -module load cray-hdf5/1.14.3.1 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -44,7 +44,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME intelimpi -setenv ICE_MACHINE_ENVINFO "ifort 2021.8.0 20221119, intel mpi 2021.13, netcdf/4.9.0.3" +setenv ICE_MACHINE_ENVINFO "ifort 2021.13.0 20240602, intel mpi 2021.16, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/configuration/scripts/machines/env.carpenter_inteloneapi b/configuration/scripts/machines/env.carpenter_inteloneapi index 2480918c1..62ba8b3c8 100644 --- a/configuration/scripts/machines/env.carpenter_inteloneapi +++ b/configuration/scripts/machines/env.carpenter_inteloneapi @@ -12,15 +12,15 @@ module unload PrgEnv-cray module unload PrgEnv-gnu module unload PrgEnv-intel module unload PrgEnv-pgi -module load PrgEnv-intel/8.5.0 +module load PrgEnv-intel/8.6.0 module unload intel -module load intel-oneapi/2024.2 +module load intel-oneapi/2025.2 module unload cray-mpich module unload mpi module unload openmpi -module load cray-mpich/8.1.30 +module load cray-mpich/8.1.32 #module load mpi/2021.11 #module load openmpi/4.1.6 @@ -29,8 +29,8 @@ module unload cray-hdf5-parallel module unload cray-netcdf-hdf5parallel module unload cray-parallel-netcdf module unload netcdf -module load cray-netcdf/4.9.0.13 -module load cray-hdf5/1.14.3.1 +module load cray-netcdf/4.9.0.17 +module load cray-hdf5/1.14.3.5 setenv NETCDF_PATH ${NETCDF_DIR} limit coredumpsize unlimited @@ -44,7 +44,7 @@ endif setenv ICE_MACHINE_MACHNAME carpenter setenv ICE_MACHINE_MACHINFO "Cray EX4000 AMD 9654 Genoa 2.1GHz, Slingshot Interconnect" setenv ICE_MACHINE_ENVNAME inteloneapi -setenv ICE_MACHINE_ENVINFO "Intel oneAPI DPC++/C++/icx/ifx 2024.2.0 20240602, cray-mpich/8.1.30, netcdf/4.9.0.13" +setenv ICE_MACHINE_ENVINFO "Intel oneAPI DPC++/C++/icx/ifx 2025.2.1 20250806, cray-mpich/8.1.32, netcdf/4.9.0.17" setenv ICE_MACHINE_MAKE gmake setenv ICE_MACHINE_WKDIR $WORKDIR/CICE_RUNS setenv ICE_MACHINE_INPUTDATA /p/app/unsupported/RASM/cice_consortium diff --git a/doc/source/developer_guide/dg_other.rst b/doc/source/developer_guide/dg_other.rst index a8f6e8b15..2819fa0eb 100644 --- a/doc/source/developer_guide/dg_other.rst +++ b/doc/source/developer_guide/dg_other.rst @@ -81,6 +81,247 @@ This is very likely to be bfb, but is not as fast or accurate as the reprosum implementation. See :cite:`He01` +.. _averages: + +Averages +----------------- + +Coupling and history output quantities may be averaged in different forms, depending on +whether the quantity represents a value averaged over the entire grid cell, the sea ice fraction, +or a subset of the sea ice fraction such as a thickness category or the ponded area. These +distinctions must also be considered for time averaging. + +The SIMIP Project :cite:`Notz16` +categorizes output variables as 'intensive' and 'extensive' based on their characteristics +relative to ice area. Extensive variables are proportional to area fraction, and their time +averages include zeroes when and where there is no ice. Intensive variables are not +proportional to area fraction, and their time averages should not include zeroes when and +where there is no ice. This is accomplished by summing area-weighted intensive values across categories +then dividing by the sum of the category areas. Tracers such as ice thickness, surface temperature, +and biogeochemical tracers are examples of intensive variables. + +The following formulas ignore subtleties such as some fluxes being computed on the initial ice area, which then +changes due to frazil ice formation, lateral melting and transport. The ice area used for both averaging and coupling should be carefully +considered in light of the model timestepping. Edge cases such as the complete disappearance or new appearance of ice +cause averaging errors. To address these cases, we could consider interpolating all quantities to the middle of the +timestep, but that is not currently done. + +Ice area +~~~~~~~~~~~~~~~~~ + +If :math:`\mathbf{X}=(x,y)`, :math:`A` is the cell area (:math:`m^2`) and :math:`g` represents +the ice thickness distribution discretized as :math:`a_n` for :math:`n=1,\, ncat`, then the +ice area (:math:`m^2`) is the sum of the thickness category areas :math:`a_n A`: + +.. math:: + A_{i}(t) = \int_{ice} g(\mathbf{X},t) \, d\mathbf{X} \sim \sum_{n=1}^{ncat} a_n(t) \, A + +and the (unitless) ice area fraction is + +.. math:: + a_{ice}(t) = {\int_{ice} g(\mathbf{X},t) \, d\mathbf{X} \over \int_{cell} d\mathbf{X} \, dt} \sim \sum_{n=1}^{ncat} a_n(t). + + +The time-averaged ice area over an interval of length :math:`N\Delta t` is + +.. math:: + \bar{A}_{i} = {\int_t A_{i}(t) \, dt \over \int_t \, dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} a_n \, A \, \Delta t \over N \, \Delta t} + = {A \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} a_n + +and the time-averaged ice area fraction is extensive (by definition): + +.. math:: + \bar{a}_{ice} = {\int_t \int_{ice} g(\mathbf{X},t) \, d\mathbf{X} \, dt \over \int_t \int_{cell} d\mathbf{X} \, dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} a_n \, A \Delta t \over A \, N \, \Delta t} + = {1 \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} a_n. + +Ice volume +~~~~~~~~~~~~~~~~~ + +Likewise for time averages of ice volume :math:`V_i` (:math:`m^3`), + +.. math:: + \bar{V}_{i} = {\int_t \int_{cell} \int_{0}^{h} g(\mathbf{X},t) \, dz \, d\mathbf{X} \, dt \over \int_{t} dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n \, A \, \Delta t \over N \, \Delta t} + = {A \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n + +for ice thickness :math:`h` assumed to be 0 in open water. Then the time-average ice volume per square meter of grid cell (:math:`m`) is + +.. math:: + \bar{v}_{ice} = {\int_t \int_{cell} \int_{0}^{h} g(\mathbf{X},t) \, dz \, d\mathbf{X} \, dt \over \int_{t} \int_{cell} d\mathbf{X} \, dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n \, A \, \Delta t \over A \, N \, \Delta t} + = {1 \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n = {1 \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} v_n. + +where :math:`v_n = h_n a_n`. :math:`v_{ice}` is the quantity labeled `hi` in history, which can be thought of as the mean ice thickness averaged over the entire +grid cell. The time-averaged ice volume per square meter of ice (mean 'actual' ice thickness, :math:`m`) is + +.. math:: + \bar{h}_{i} = {\int_t \int_{ice} \int_{0}^{h} g(\mathbf{X},t) \, dz \, d\mathbf{X} \, dt \over \int_{t} \int_{ice} d\mathbf{X} \, dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n \, A \, \Delta t \over \sum_{\Delta t} \sum_{n=1}^{ncat} a_n \, A \, \Delta t} + = {\sum_{\Delta t} \sum_{n=1}^{ncat} v_n \over \sum_{\Delta t} \sum_{n=1}^{ncat} a_n}. + +Snow volume is treated similarly. Ice and snow volumes are extensive, while thicknesses are +intensive. + +The form used here for time-averaging the average 'actual' thickness produces the average over all ice present +during the averaging interval. For intensive variables in particular, this form is slightly different from +the time-average of the category-averaged quantity per time step. The latter, two-step averaging process +requires additional divisions and re-multiplications by ice area, introducing errors where ice areas +are very small or cells change from ice-free to having ice or vice versa. The same is true for other tracers +and intensive variables. While both approaches are valid, averages as written here are preferred when +conservation is important. + +Volume content +~~~~~~~~~~~~~~~~~ + +Total content of tracers such as salt and enthalpy are necessary for conservative coupling. The time-average content +of a volume tracer :math:`b` (with units per :math:`m^3`) is + +.. math:: + \bar{B}_{i} = {\int_t \int_{cell} \int_{0}^{h} b(\mathbf{X},z,t) g(\mathbf{X},t) \, dz \, d\mathbf{X} \, dt \over \int_{t} dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} b_n \, h_n \, a_n \, A \, \Delta t \over N \, \Delta t} + = {A \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} b_n \, v_n + +and the time-averaged content per square meter of grid cell is + +.. math:: + \bar{b}_{ice} \sim {1 \over N} \sum_{\Delta t} \sum_{n=1}^{ncat} b_n \, v_n. + +The mean tracer value in sea ice is + +.. math:: + \bar{b}_{i} = {\int_t \int_{cell} \int_{0}^{h} b(\mathbf{X},z,t) g(\mathbf{X},t) \, dz \, d\mathbf{X} \, dt \over \int_{t} \int_{cell} \int_{0}^{h} dz \, d\mathbf{X} \, dt} + \sim {\sum_{\Delta t} \sum_{n=1}^{ncat} b_n \, h_n \, a_n \, A \, \Delta t \over \sum_{\Delta t} \sum_{n=1}^{ncat} h_n \, a_n \, A \, \Delta t} + = {\sum_{\Delta t} \sum_{n=1}^{ncat} b_n \, v_n \over \sum_{\Delta t} \sum_{n=1}^{ncat} v_n} + +Thus, volume content variables are extensive, while the tracers themselves are intensive. + +Surface quantities +~~~~~~~~~~~~~~~~~ + +Surface quantities such as temperature are intensive and treated similarly to volume tracers, with integrals taken over +the desired surface area rather than the volume. For example, + +.. math:: + T_{ice}(t) = {\int_{ice} T(\mathbf{X},t) g(\mathbf{X},t) \, d\mathbf{X} \over \int_{ice} g(\mathbf{X},t) \, d\mathbf{X}} + +and the time average is simply + +.. math:: + \bar{T}_{ice} = {\sum_{\Delta t} \sum_{n=1}^{ncat} T_n \, a_n \over \sum_{\Delta t} \sum_{n=1}^{ncat} \, a_n}. + +Note that since :math:`\sum_{n=0}^{ncat} \, a_n \,=\, 1`, a category-merged quantity can be considered the average over the cell area, assuming +the quantity is zero over open water: + +.. math:: + T_{cell} = {\sum_{n=0}^{ncat} T_n \, a_n \over \sum_{n=0}^{ncat} \, a_n} = \sum_{n=1}^{ncat} \, T_n \, a_n, + +and the average value over the ice is then + +.. math:: + T_{ice} = {\sum_{n=1}^{ncat} T_n \, a_n \over \sum_{n=1}^{ncat} \, a_n} = {T_{cell} \over a_{ice}}. + +This simplification is applicable for tracers carried on the ice area (or volume, similarly), which are zero over open water by definition. +When time-averaging CICE's history fields, the category-merged value in the numerator is saved (usually in Icepack), then accumulated in time and +later divided by the accumulated ice area fraction (or volume) in CICE. + + + +Tracer hierarchies +~~~~~~~~~~~~~~~~~ + +For tracers that are carried on other tracers, such as melt ponds, averages over different areas of a given cell differ in the denominator. +For melt ponds not carried on the level-ice area, for example, the average pond depths over the grid cell area, the ice area, and the ponded +area are, respectively, + +.. math:: + h_{p\,cell} = \frac{ \int_{cell} h_p \, a_p \, g \, d\mathbf{X} } + { \int_{cell} d\mathbf{X} } + \sim \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_n + +.. math:: + h_{p\,ice} = \frac{ \int_{ice} h_p \, a_p \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + = \frac{ \int_{cell} h_p \, a_p \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_n }{ \sum_{n=1}^{ncat} a_n } + +.. math:: + h_{p\,pond} = \frac{ \int_{pond} h_p \, a_p \, g \, d\mathbf{X} } + { \int_{pond} a_p \, g \, d\mathbf{X} } + = \frac{ \int_{cell} h_p \, a_p \, g \, d\mathbf{X} } + { \int_{ice} a_p \, g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_n }{ \sum_{n=1}^{ncat} a_{pn} \, a_n }. + +For level-ice ponds, there is an extra factor of :math:`a_{lvl}`. The level-ice pond depth averaged over the grid cell area, total ice area, level ice area and pond area are + +.. math:: + h_{p\,cell} = \frac{ \int_{cell} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{cell} d\mathbf{X} } + \sim \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_{lvln} \, a_n + +.. math:: + h_{p\,ice} = \frac{ \int_{ice} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + = \frac{ \int_{cell} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_{lvln} \, a_n }{ \sum_{n=1}^{ncat} a_n } + +.. math:: + h_{p\,lvl} = \frac{ \int_{lvl} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{lvl} a_{lvl} \, g \, d\mathbf{X} } + = \frac{ \int_{cell} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{ice} a_{lvl} \, a_{pn} \, g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_{lvln} \, a_n }{ \sum_{n=1}^{ncat} a_{lvln} \, a_n } + +.. math:: + h_{p\,pond} = \frac{ \int_{pond} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{pond} a_p \, a_{lvl} \, g \, d\mathbf{X} } + = \frac{ \int_{cell} h_p \, a_p \, a_{lvl} \, g \, d\mathbf{X} } + { \int_{ice} a_p \, a_{lvl} \, g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} h_{pn} \, a_{pn} \, a_{lvln} \, a_n }{ \sum_{n=1}^{ncat} a_{pn} \, a_{lvln} \, a_n }. + +Time averages follow analogously as above. + +Ridged-ice area and volume are handled slightly differently, since they are diagnostic based on +level-ice area and volume. Level-ice area is a tracer on ice area, and level-ice volume is a +tracer on ice volume. The tracer values are fractions of the total ice, and ridged (deformed) ice is +diagnosed as the remainder of the ice fraction or volume: +:math:`T_{ardg} = 1 - T_{alvl}` and :math:`T_{vrdg} = 1 - T_{vlvl}` for the area and volume tracers. +Thus the mean level and ridged ice area fractions of the ice area are + +.. math:: + a_{lvl\,ice} = \frac{ \int_{ice} T_{alvl} \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} a_{lvln} \, a_n }{ \sum_{n=1}^{ncat} a_n } + +.. math:: + a_{rdg\,ice} = \frac{ \int_{ice} (1 - T_{alvl}) \, g \, d\mathbf{X} } + { \int_{ice} g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} (1 - a_{lvln}) \, a_n }{ \sum_{n=1}^{ncat} a_n }. + +The mean thickness of level ice, averaging over just the level-ice areas from all categories, is + +.. math:: + h_{lvl} = \frac{ \int_{ice} \int_{0}^{h} T_{vlvl} \, g \, dz \, d\mathbf{X} } + { \int_{ice} T_{alvl} \, g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} T_{vlvln} \, a_n \, h_n } + { \sum_{n=1}^{ncat} T_{alvln} \, a_n } + \sim \frac{ \sum_{n=1}^{ncat} T_{vlvln} \, v_n } + { \sum_{n=1}^{ncat} T_{alvln} \, a_n } + +and the mean thickness of deformed ice (averaging over just the ridged-ice +areas from all categories) is + +.. math:: + h_{rdg} = \frac{ \int_{ice} \int_{0}^{h} (1 - T_{vlvl}) \, g \, dz \, d\mathbf{X} } + { \int_{ice} (1 - T_{alvl}) \, g \, d\mathbf{X} } + \sim \frac{ \sum_{n=1}^{ncat} (1 - T_{vlvln}) \, a_n \, h_n } + { \sum_{n=1}^{ncat} (1 - T_{alvln}) \, a_n } + \sim \frac{ \sum_{n=1}^{ncat} (1 - T_{vlvln}) \, v_n } + { \sum_{n=1}^{ncat} (1 - T_{alvln}) \, a_n }. + .. _addtimer: Adding Timers diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 3882073de..486dfdfc1 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -311,6 +311,7 @@ grid_nml "``grid_format``", "``bin``", "read direct access grid and kmt files", "``bin``" "", "``geosnc``", "read grid and kmt file in GEOS netcdf format", "" "", "``pop_nc``", "read grid and kmt files in POP netcdf format", "" + "", "``pop_nc_ext``", "read extended grid and kmt files in POP netcdf format", "" "", "``meshnc``", "coupled model grid option, no CICE code support", "" "", "``mom_nc``", "read grid in MOM (supergrid) format and kmt files", "" "``grid_ice``", "``B``", "use B grid structure with T at center and U at NE corner", "``B``" diff --git a/icepack b/icepack index eedb51924..0bcde2556 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit eedb519247f48bce9fa1d1b275b5cc3dc07643e4 +Subproject commit 0bcde255637a5947b1b7a4e4fc8dccd77803cb65