From f96ebd4231d732e981b780cec6c1c11b2da41067 Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Fri, 6 Dec 2024 09:19:27 +1100 Subject: [PATCH] Add implementation for CICE to read the mom supergrid file directly, and use this as the grid file. --- cicecore/cicedyn/general/ice_init.F90 | 7 +- cicecore/cicedyn/infrastructure/ice_grid.F90 | 1057 ++++++++++++++++-- doc/source/user_guide/ug_case_settings.rst | 1 + doc/source/user_guide/ug_implementation.rst | 15 +- icepack | 2 +- 5 files changed, 960 insertions(+), 122 deletions(-) diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index 5cbaedcb6..5552b0f2d 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -1948,6 +1948,7 @@ subroutine input_data write(nu_diag,*) ' ' write(nu_diag,*) ' Grid, Discretization' write(nu_diag,*) '--------------------------------' + write(nu_diag,1030) ' grid_format = ',trim(grid_format) tmpstr2 = ' ' if (trim(grid_type) == 'rectangular') tmpstr2 = ' : internally defined, rectangular grid' if (trim(grid_type) == 'regional') tmpstr2 = ' : user-defined, regional grid' @@ -2692,14 +2693,16 @@ subroutine input_data kmt_type /= 'channel_onenorth' .and. & kmt_type /= 'wall' .and. & kmt_type /= 'default' .and. & - kmt_type /= 'boxislands') then + kmt_type /= 'boxislands'.and. & + kmt_type /= 'none' ) then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: unknown kmt_type=',trim(kmt_type) abort_list = trim(abort_list)//":27" endif if (grid_type /= 'column' .and. & grid_type /= 'rectangular' .and. & - kmt_type /= 'file') then + kmt_type /= 'file' .and. & + kmt_type /= 'none') then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: need kmt file, kmt_type=',trim(kmt_type) abort_list = trim(abort_list)//":28" endif diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 53082aeaf..047f27a2a 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -25,7 +25,7 @@ module ice_grid use ice_kinds_mod use ice_broadcast, only: broadcast_scalar, broadcast_array use ice_boundary, only: ice_HaloUpdate, ice_HaloExtrapolate - use ice_constants, only: c0, c1, c1p5, c2, c4, c20, c360, & + use ice_constants, only: c0, c1, c1p5, c2, c4, c20, c180, c360, & p5, p25, radius, cm_to_m, m_to_cm, & field_loc_center, field_loc_NEcorner, field_loc_Nface, field_loc_Eface, & field_type_scalar, field_type_vector, field_type_angle @@ -53,7 +53,7 @@ module ice_grid grid_neighbor_min, grid_neighbor_max character (len=char_len_long), public :: & - grid_format , & ! file format ('bin'=binary or 'nc'=netcdf) + grid_format , & ! file format ('bin'=binary or 'nc'=netcdf or 'mom_nc'=mom (supergrid) netcdf) gridcpl_file , & ! input file for POP coupling grid info grid_file , & ! input file for POP grid info kmt_file , & ! input file for POP grid info @@ -172,7 +172,6 @@ module ice_grid character (len=char_len), private :: & mask_fieldname !field/var name for the mask variable (in nc files) - interface grid_average_X2Y module procedure grid_average_X2Y_base , & grid_average_X2Y_userwghts, & @@ -302,17 +301,19 @@ subroutine init_grid1 integer (kind=int_kind) :: & fid_grid, & ! file id for netCDF grid file - fid_kmt ! file id for netCDF kmt file + fid_kmt ! file id for netCDF kmt file character (char_len) :: & fieldname ! field name in netCDF file real (kind=dbl_kind), dimension(:,:), allocatable :: & - work_g1, work_g2 + work_g1, work_g2, work_mom integer (kind=int_kind) :: & max_blocks_min, & ! min value of max_blocks across procs - max_blocks_max ! max value of max_blocks across procs + max_blocks_max, & ! max value of max_blocks across procs + i, j, im, jm , & + ierr real (kind=dbl_kind) :: & rad_to_deg @@ -328,75 +329,120 @@ subroutine init_grid1 if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - allocate(work_g1(nx_global,ny_global)) - allocate(work_g2(nx_global,ny_global)) + allocate( & + work_g1(nx_global,ny_global), & + work_g2(nx_global,ny_global), & + stat=ierr & + ) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) ! check tripole flags here ! can't check in init_data because ns_boundary_type is not yet read ! can't check in init_domain_blocks because grid_type is not accessible due to circular logic - if (grid_type == 'tripole' .and. ns_boundary_type /= 'tripole' .and. & + if (trim(grid_type) == 'tripole' .and. ns_boundary_type /= 'tripole' .and. & ns_boundary_type /= 'tripoleT') then call abort_ice(subname//' ERROR: grid_type tripole needs tripole ns_boundary_type', & file=__FILE__, line=__LINE__) endif - if (grid_type == 'tripole' .and. (mod(nx_global,2)/=0)) then + if (trim(grid_type) == 'tripole' .and. (mod(nx_global,2)/=0)) then call abort_ice(subname//' ERROR: grid_type tripole requires even nx_global number', & file=__FILE__, line=__LINE__) endif + if (trim(grid_format) == 'mom_nc' .and. ns_boundary_type == 'tripoleT') then + call abort_ice(subname//" ERROR: ns_boundary_type='tripoleT' not implemented "& + "for grid_format='mom_nc' tripole needs tripole ns_boundary_type", & + file=__FILE__, line=__LINE__) + endif + if (trim(grid_type) == 'displaced_pole' .or. & trim(grid_type) == 'tripole' .or. & trim(grid_type) == 'regional' ) then - if (trim(grid_format) == 'nc') then - - fieldname='ulat' - call ice_open_nc(grid_file,fid_grid) - call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.) - call ice_close_nc(fid_grid) + ! Fill ULAT + select case(trim(grid_format)) + case ('mom_nc') + + if (my_task == master_task) then + allocate(work_mom(nx_global*2+1, ny_global*2+1), stat=ierr) + else + allocate(work_mom(1, 1), stat=ierr) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + fieldname='y' ! use mom y field to fill cice ULAT + call ice_open_nc(grid_file,fid_grid) + call ice_read_global_nc(fid_grid,1,fieldname,work_mom,.true.) + call ice_close_nc(fid_grid) + im = 3 + do i = 1, nx_global + jm = 3 + do j = 1, ny_global + work_g1(i,j) = work_mom(im, jm) + jm = jm + 2 + enddo + im = im + 2 + enddo + + deallocate(work_mom, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + case('nc') + + fieldname='ulat' + call ice_open_nc(grid_file,fid_grid) + call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.) + call ice_close_nc(fid_grid) + + case default + + call ice_open(nu_grid,grid_file,64) + call ice_read_global(nu_grid,1,work_g1,'rda8',.true.) + if (my_task == master_task) close (nu_grid) + + end select + + else ! rectangular grid + work_g1(:,:) = 75._dbl_kind/rad_to_deg ! arbitrary polar latitude + endif + + ! Fill kmt + if (trim(kmt_type) =='file') then + select case(trim(grid_format)) + case ('mom_nc', 'nc') ! mask variable name might be kmt or mask, check both - call ice_open_nc(kmt_file,fid_kmt) + call ice_open_nc(kmt_file,fid_kmt) #ifdef USE_NETCDF - if ( my_task==master_task ) then - status = nf90_inq_varid(fid_kmt, 'kmt', varid) - if (status == nf90_noerr) then - mask_fieldname = 'kmt' - else - status = nf90_inq_varid(fid_kmt, 'mask', varid) - call ice_check_nc(status, subname//' ERROR: does '//trim(kmt_file)//& - ' contain "kmt" or "mask" variable?', file=__FILE__, line=__LINE__) - mask_fieldname = 'mask' + if ( my_task==master_task ) then + status = nf90_inq_varid(fid_kmt, 'kmt', varid) + if (status == nf90_noerr) then + mask_fieldname = 'kmt' + else + status = nf90_inq_varid(fid_kmt, 'mask', varid) + call ice_check_nc(status, subname//' ERROR: does '//trim(kmt_file)//& + ' contain "kmt" or "mask" variable?', file=__FILE__, line=__LINE__) + mask_fieldname = 'mask' + endif endif - endif #endif - call broadcast_scalar(mask_fieldname, master_task) - - call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g2,.true.) - call ice_close_nc(fid_kmt) + call broadcast_scalar(mask_fieldname, master_task) - else - - call ice_open(nu_grid,grid_file,64) ! ULAT - call ice_open(nu_kmt, kmt_file, 32) ! KMT + call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g2,.true.) + call ice_close_nc(fid_kmt) - call ice_read_global(nu_grid,1,work_g1,'rda8',.true.) ! ULAT - call ice_read_global(nu_kmt, 1,work_g2,'ida4',.true.) ! KMT + case default - if (my_task == master_task) then - close (nu_grid) - close (nu_kmt) - endif + call ice_open(nu_kmt, kmt_file, 32) ! KMT + call ice_read_global(nu_kmt, 1,work_g2,'ida4',.true.) ! KMT + if (my_task == master_task) close (nu_kmt) - endif - - else ! rectangular grid + end select - work_g1(:,:) = 75._dbl_kind/rad_to_deg ! arbitrary polar latitude + else work_g2(:,:) = c1 - endif call broadcast_array(work_g1, master_task) ! ULAT @@ -408,8 +454,8 @@ subroutine init_grid1 call init_domain_distribution(work_g2, work_g1, grid_ice) ! KMT, ULAT - deallocate(work_g1) - deallocate(work_g2) + deallocate(work_g1, work_g2, stat = ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) !----------------------------------------------------------------- ! write additional domain information @@ -449,8 +495,9 @@ subroutine init_grid2 #endif integer (kind=int_kind) :: & - i, j, iblk, & - ilo,ihi,jlo,jhi ! beginning and end of physical domain + i, j, iblk, & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + ierr real (kind=dbl_kind) :: & angle_0, angle_w, angle_s, angle_sw, & @@ -486,13 +533,14 @@ subroutine init_grid2 if (trim(grid_type) == 'displaced_pole' .or. & trim(grid_type) == 'tripole' .or. & trim(grid_type) == 'regional' ) then - if (trim(grid_format) == 'nc') then - call kmtmask_nc ! read mask from nc file - call popgrid_nc ! read POP grid lengths from nc file - else - call kmtmask ! read kmt directly - call popgrid ! read POP grid lengths directly - endif + select case (trim(grid_format)) + case('mom_nc') + call mom_grid ! derive cice grid from mom supergrid nc file + case ('nc') + call popgrid_nc ! read POP grid lengths from nc file + case default + call popgrid ! read POP grid lengths directly + end select #ifdef CESMCOUPLED elseif (trim(grid_type) == 'latlon') then call latlongrid ! lat lon grid for sequential CESM (CAM mode) @@ -502,6 +550,18 @@ subroutine init_grid2 call rectgrid ! regular rectangular grid endif + if (trim(kmt_type) =='none') then + kmt(:,:,:) = c1 + hm(:,:,:) = c1 + else if (trim(kmt_type) =='file') then + select case (trim(grid_format)) + case('mom_nc', 'nc') + call kmtmask_nc + case default + call kmtmask + end select + endif !the other types are handled by rectgrid + !----------------------------------------------------------------- ! Diagnose OpenMP thread schedule, force order in output !----------------------------------------------------------------- @@ -533,6 +593,26 @@ subroutine init_grid2 ! at halos. !----------------------------------------------------------------- + if (trim(grid_format) /= 'mom_nc') then + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + 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 = 1,ny_block + do i = 1,nx_block + tarea(i,j,iblk) = dxT(i,j,iblk)*dyT(i,j,iblk) + uarea(i,j,iblk) = dxU(i,j,iblk)*dyU(i,j,iblk) + narea(i,j,iblk) = dxN(i,j,iblk)*dyN(i,j,iblk) + earea(i,j,iblk) = dxE(i,j,iblk)*dyE(i,j,iblk) + enddo + enddo + enddo + endif + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) @@ -543,11 +623,6 @@ subroutine init_grid2 do j = 1,ny_block do i = 1,nx_block - tarea(i,j,iblk) = dxT(i,j,iblk)*dyT(i,j,iblk) - uarea(i,j,iblk) = dxU(i,j,iblk)*dyU(i,j,iblk) - narea(i,j,iblk) = dxN(i,j,iblk)*dyN(i,j,iblk) - earea(i,j,iblk) = dxE(i,j,iblk)*dyE(i,j,iblk) - if (tarea(i,j,iblk) > c0) then tarear(i,j,iblk) = c1/tarea(i,j,iblk) else @@ -616,8 +691,7 @@ subroutine init_grid2 call ice_timer_stop(timer_bound) !----------------------------------------------------------------- - ! Calculate ANGLET to be compatible with POP ocean model - ! First, ensure that -pi <= ANGLE <= pi + ! Ensure that -pi <= ANGLE <= pi !----------------------------------------------------------------- out_of_range = .false. @@ -625,7 +699,17 @@ subroutine init_grid2 if (count(out_of_range) > 0) then write(nu_diag,*) subname,' angle = ',minval(ANGLE),maxval(ANGLE),count(out_of_range) call abort_ice (subname//' ANGLE out of expected range', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__) + endif + + if (l_readCenter) then + out_of_range = .false. + where (ANGLET < -pi .or. ANGLET > pi) out_of_range = .true. + if (count(out_of_range) > 0) then + write(nu_diag,*) subname,' angle = ',minval(ANGLET),maxval(ANGLET),count(out_of_range) + call abort_ice (subname//' ANGLET out of expected range', & + file=__FILE__, line=__LINE__) + endif endif if (l_readCenter) then @@ -701,10 +785,21 @@ subroutine init_grid2 call ice_timer_stop(timer_bound) call makemask ! velocity mask, hemisphere masks - if (.not. (l_readCenter)) then - call Tlatlon ! get lat, lon on the T grid + + !---------------------------------------------------------------- + ! Coordinates for all T/N/E cells + !---------------------------------------------------------------- + + if (trim(grid_format) /= 'mom_nc') then + if (.not. (l_readCenter)) then + call Tlatlon ! get lat, lon on the T grid + endif + call NElatlon ! get lat, lon on the N, E grid + + ! corners for CF-compliant output + call gridbox_corners + call gridbox_edges endif - call NElatlon ! get lat, lon on the N, E grid !----------------------------------------------------------------- ! bathymetry @@ -719,33 +814,29 @@ subroutine init_grid2 file=__FILE__, line=__LINE__) endif - !---------------------------------------------------------------- - ! Corner coordinates for CF compliant history files - !---------------------------------------------------------------- - - call gridbox_corners - call gridbox_edges - !----------------------------------------------------------------- ! Compute global index (used for unpacking messages from coupler) !----------------------------------------------------------------- if (my_task==master_task) then - allocate(work_g1(nx_global,ny_global)) + allocate(work_g1(nx_global,ny_global), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) do j=1,ny_global do i=1,nx_global work_g1(i,j) = real((j-1)*nx_global + i,kind=dbl_kind) enddo enddo else - allocate(work_g1(1,1)) ! to save memory + allocate(work_g1(1,1), stat=ierr) ! to save memory + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) endif call scatter_global(rndex_global, work_g1, & master_task, distrb_info, & field_loc_center, field_type_scalar) - deallocate(work_g1) + deallocate(work_g1, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) end subroutine init_grid2 @@ -827,6 +918,8 @@ subroutine popgrid real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g1 + integer (int_kind) :: ierr + character(len=*), parameter :: subname = '(popgrid)' call ice_open(nu_grid,grid_file,64) @@ -837,7 +930,8 @@ subroutine popgrid ! lat, lon, angle !----------------------------------------------------------------- - allocate(work_g1(nx_global,ny_global)) + allocate(work_g1(nx_global,ny_global), stat=ierr) + 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) @@ -869,7 +963,8 @@ subroutine popgrid call ice_read_global(nu_grid,4,work_g1,'rda8',.true.) ! HTE call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE - deallocate(work_g1) + deallocate(work_g1, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) if (my_task == master_task) then close (nu_grid) @@ -948,13 +1043,15 @@ end subroutine kmtmask_nc subroutine popgrid_nc #ifdef USE_NETCDF - 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 + fid_grid , & ! file id for netCDF grid file + ierr logical (kind=log_kind) :: diag @@ -986,7 +1083,8 @@ subroutine popgrid_nc ! lat, lon, angle !----------------------------------------------------------------- - allocate(work_g1(nx_global,ny_global)) + allocate(work_g1(nx_global,ny_global), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) fieldname='ulat' call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULAT @@ -1015,7 +1113,7 @@ subroutine popgrid_nc ! if grid file includes anglet then read instead fieldname='anglet' if (my_task == master_task) then - status = nf90_inq_varid(fid_grid, trim(fieldname) , varid) + 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 @@ -1052,7 +1150,8 @@ subroutine popgrid_nc call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTE call primary_grid_lengths_HTE(work_g1) ! dyU, dyT, dyN, dyE - deallocate(work_g1) + deallocate(work_g1, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) call ice_close_nc(fid_grid) @@ -1075,7 +1174,8 @@ subroutine latlongrid use ice_scam, only : scmlat, scmlon, single_column #ifdef USE_NETCDF - use netcdf + use netcdf, only : nf90_inq_varid , nf90_inq_dimid, & + nf90_inquire_dimension, nf90_get_var #endif integer (kind=int_kind) :: & @@ -1344,6 +1444,706 @@ subroutine latlongrid end subroutine latlongrid #endif + +!======================================================================= +! Create the CICE grid from the MOM supergrid netcdf file. +! CICE Fields and units are: +! ULAT, ULON, TLAT, TLON, ELAT, ELON, NLAT, NLON (radians) +! HTN, HTE (m) +! dxT, dyT, dxU, dyU, dxN, dyN, dxE, dyE, (m) +! ANGLE, ANGLET (radians) +! tarea, uarea, narea, earea (m^2) +! lont_bounds, latt_bounds, etc (degrees) + + subroutine mom_grid + + integer (kind=int_kind) :: & + fid_grid, & ! file id for netCDF grid file + varid, & ! netcdf varid + ierr + + logical (kind=log_kind) :: diag + + character (char_len) :: & + fieldname ! field name in netCDF file + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + G_TLON, work_gE, G_ULON, work_gN, work_mom, G_ULAT, G_TLAT, work_g1 + + character(len=*), parameter :: subname = '(mom_grid)' + + call ice_open_nc(grid_file,fid_grid) + + !----------------------------------------------------------------- + ! lat, lon, angle + !----------------------------------------------------------------- + + if (my_task == master_task) then + allocate( & + work_mom(nx_global*2+1, ny_global*2+1), & + work_gE(nx_global+1,ny_global+1) , & + work_gN(nx_global+1,ny_global+1) , & + G_ULAT(nx_global+1,ny_global+1) , & !1 bigger to include left and bottom + G_TLAT(nx_global+1,ny_global+1) , & !1 bigger to include top and right + G_TLON(nx_global+1,ny_global+1) , & !1 bigger to include left and bottom + G_ULON(nx_global+1,ny_global+1) , & !1 bigger to include top and right + stat = ierr & + ) + else + allocate(work_mom(1,1), work_gE(1,1), work_gN(1,1), & + G_ULAT(1,1), G_TLAT(1,1), G_TLON(1,1), G_ULON(1,1), & + stat=ierr) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + ! populate all LAT fields + fieldname='y' + call ice_read_global_nc(fid_grid,1,fieldname,work_mom,diag) + call mom_corners_global(work_mom, G_ULAT, G_TLAT, work_gE, work_gN) + ! create bounds fields for cf-compliant output + call mom_bounds(G_ULAT, latt_bounds) ! u points define corners for t-cells + call mom_bounds(G_TLAT, latu_bounds) + call mom_bounds(work_gN, late_bounds) + call mom_bounds(work_gE, latn_bounds) + !distribute global array to local + call mom_corners_scatter(G_ULAT, G_TLAT, work_gE, work_gN, & + ULAT, TLAT, ELAT, NLAT) + + ! populate all LON fields + fieldname='x' + call ice_read_global_nc(fid_grid,1,fieldname,work_mom,diag) + call mom_corners_global(work_mom, G_ULON, G_TLON, work_gE, work_gN) + call mom_bounds(G_ULON, lont_bounds) + call mom_bounds(G_TLON, lonu_bounds) + call mom_bounds(work_gN, lone_bounds) + call mom_bounds(work_gE, lonn_bounds) + call mom_corners_scatter(G_ULON, G_TLON, work_gE, work_gN, & + ULON, TLON, ELON, NLON) + + deallocate(work_gE, work_gN, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + if (my_task == master_task) then + allocate(work_g1(nx_global, ny_global), stat=ierr) !array for angle field + else + allocate(work_g1(1, 1), stat=ierr) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + ! populate angle fields, angle is u-points, angleT is t-points + ! even though mom supergrid files contain angle_dx, mom6 calculates internally + call grid_rotation_angle(G_ULON, G_ULAT, G_TLON(1:nx_global,1:ny_global), work_g1) ! anglet + call scatter_global(ANGLET, work_g1, master_task, distrb_info, & + field_loc_center, field_type_angle) + call grid_rotation_angle(G_TLON, G_TLAT, G_ULON(2:nx_global+1,2:ny_global+1), work_g1) ! angle + call scatter_global(ANGLE, work_g1, master_task, distrb_info, & + field_loc_NEcorner, field_type_angle) + + deallocate(work_g1, G_ULAT, G_TLAT, G_TLON, G_ULON, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + !----------------------------------------------------------------- + ! cell dimensions + !----------------------------------------------------------------- + fieldname='dx' + ! dx uses the cells in x, edges in y, reallocate work_mom to this size + deallocate(work_mom, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + if (my_task == master_task) then + allocate(work_mom(nx_global*2, ny_global*2+1), stat=ierr) + else + allocate(work_mom(1, 1), stat=ierr) + endif + 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_mom,diag) + call mom_dx(work_mom) + deallocate(work_mom, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + fieldname='dy' + ! dy uses the edges in x, cells in y, reallocate work_mom to this size + if (my_task == master_task) then + allocate(work_mom(nx_global*2+1, ny_global*2), stat=ierr) + else + allocate(work_mom(1, 1), stat=ierr) + endif + 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_mom,diag) + call mom_dy(work_mom) + deallocate(work_mom, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + + !----------------------------------------------------------------- + ! cell areas + !----------------------------------------------------------------- + fieldname = 'area' + if (my_task == master_task) then + allocate(work_mom(nx_global*2, ny_global*2), stat=ierr) + else + allocate(work_mom(1, 1), stat=ierr) + endif + 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_mom,diag) + call mom_area(work_mom) + deallocate(work_mom, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc', file=__FILE__, line=__LINE__) + + !----------------------------------------------------------------- + ! fin + !----------------------------------------------------------------- + call ice_close_nc(fid_grid) + l_readCenter = .true. ! we have read t quantities + + end subroutine mom_grid + + subroutine mom_corners_global(work_mom, G_U, G_T, G_E, G_N) + + ! mom supergrid has four cells for every model cell + ! we need to select the correct edges to get lat & lon for a model cell + ! we include left/bottom edges for U-points, and top/right edges for T-points + ! and close per ew_boundary_type & ns_boundary_type + + real (kind=dbl_kind), dimension(:,:), intent(in) :: work_mom + ! supergrid array of x or y + + real (kind=dbl_kind), dimension(:,:), intent(out) :: G_U, G_T, G_E, G_N + ! global grids + + integer (kind=int_kind) :: & + i, j, & + im1, im2, jm1, jm2 ! i & j for mom supergrid + + character(len=*), parameter :: subname = '(mom_corners_global)' + + if (my_task == master_task) then + + im1 = 1 ; im2 = 2 ! lh , middle hand edge of first col + do i = 1, nx_global + jm1 = 1; jm2 = 2 ! bottom, middle of first row + do j = 1, ny_global + G_U(i,j) = work_mom(im1, jm1) ! ULAT/LON + G_N(i,j) = work_mom(im2, jm1) ! NLAT/LON + G_E(i,j) = work_mom(im1, jm2) ! ELAT/LON + G_T(i,j) = work_mom(im2, jm2) ! TLAT/LON + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + + ! fill last col + jm1 = 1; jm2 = 2 ! bottom, middle of first row + do j = 1, ny_global + G_U(nx_global+1,j) = work_mom(2*nx_global+1, jm1) + G_E(nx_global+1,j) = work_mom(2*nx_global+1, jm2) + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + select case (trim(ew_boundary_type)) + case('cyclic') + G_T(nx_global+1,:) = G_T(1,:) + G_N(nx_global+1,:) = G_N(1,:) + case('open') + do j=1, ny_global+1 + G_T(nx_global+1,j) = 2 * G_T(nx_global, j) - G_T(nx_global-1, j) + G_N(nx_global+1,j) = 2 * G_N(nx_global, j) - G_N(nx_global-1, j) + enddo + end select + + ! fill last row + im1 = 1 ; im2 = 2 + do i = 1, nx_global+1 + G_U(i,ny_global + 1) = work_mom(im1, 2*ny_global+1) + G_N(i,ny_global + 1) = work_mom(im2, 2*ny_global+1) + im1 = im1 + 2 + enddo + select case (trim(ns_boundary_type)) + case ('tripole') + do i = 1, nx_global+1 + G_T(i,ny_global+1) = G_T(nx_global+1-i, ny_global) + G_E(i,ny_global+1) = G_E(nx_global+1-i, ny_global) + enddo + case ('cyclic') + G_T(:,ny_global+1) = G_T(:,1) + G_E(:,ny_global+1) = G_E(:,1) + case ('open') + do i = 1, nx_global+1 + G_T(i,ny_global+1) = 2 * G_T(i, ny_global) - G_T(i, ny_global-1) + G_E(i,ny_global+1) = 2 * G_E(i, ny_global) - G_E(i, ny_global-1) + enddo + end select + + endif + + end subroutine mom_corners_global + + + subroutine mom_bounds(G_corners, bounds) + + ! with an global array of corner point, subset and distribute + ! into the cice bounds variables + + real (kind=dbl_kind), dimension(:,:), intent(in) :: G_corners + real (kind=dbl_kind), dimension(:,:,:,:), intent(out) :: bounds + + ! local vars + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + work_bounds + + character(len=*), parameter :: subname = '(mom_bounds)' + + ! Get bounds of grid boxes for each block as follows: + ! (1) SW corner, (2) SE corner, (3) NE corner, (4) NW corner + call scatter_global(work_bounds, G_corners(1:nx_global, 1:ny_global), & + master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + bounds(1,:,:,:) = work_bounds(:,:,:) + call scatter_global(work_bounds, G_corners(2:nx_global+1, 1:ny_global), & + master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + bounds(2,:,:,:) = work_bounds(:,:,:) + call scatter_global(work_bounds, G_corners(2:nx_global+1, 2:ny_global+1), & + master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + bounds(3,:,:,:) = work_bounds(:,:,:) + call scatter_global(work_bounds, G_corners(1:nx_global, 2:ny_global+1), & + master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + bounds(4,:,:,:) = work_bounds(:,:,:) + + end subroutine mom_bounds + + subroutine mom_corners_scatter(G_U, G_T, G_E, G_N, U, T, E, N ) + + ! with a global array of corner points in degrees, convert to rad and scatter to workers + + real (kind=dbl_kind), dimension(:,:), intent(inout) :: G_U, G_T, G_E, G_N + ! global grids + + real (kind=dbl_kind), dimension(:,:,:), intent(out) :: U, T, E, N ! local grids + + real (kind=dbl_kind) :: deg_to_rad , pi + + character(len=*), parameter :: subname = '(mom_corners_scatter)' + + 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__) + + deg_to_rad = pi/c180 + + ! convert to rad + G_T = G_T * deg_to_rad + G_U = G_U * deg_to_rad + G_N = G_N * deg_to_rad + G_E = G_E * deg_to_rad + + ! distribute to processors + ! subset G_T to active cells by dropping right/top halo + call scatter_global(T, G_T(1:nx_global, 1:ny_global), & + master_task, distrb_info, & + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(T, distrb_info, & + ew_boundary_type, ns_boundary_type) + ! subset G_U/G_E/G_N to active cells by dropping left/bottom edge + call scatter_global(U, G_U(2:nx_global+1, 2:ny_global+1), & + master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + call ice_HaloExtrapolate(U, distrb_info, & + ew_boundary_type, ns_boundary_type) + call scatter_global(N, G_N(1:nx_global, 2:ny_global+1), master_task, distrb_info, & + field_loc_Nface, field_type_scalar) + call ice_HaloExtrapolate(N, distrb_info, & + ew_boundary_type, ns_boundary_type) + call scatter_global(E, G_E(2:nx_global+1, 1:ny_global), master_task, distrb_info, & + field_loc_Eface, field_type_scalar) + call ice_HaloExtrapolate(E, distrb_info, & + ew_boundary_type, ns_boundary_type) + + end subroutine mom_corners_scatter + + + subroutine mom_dx(work_mom) + + ! mom supergrid has four cells for every model cell, + ! we just need to sum the correct sidelengths to get dx + + real (kind=dbl_kind), dimension(:,:) :: work_mom + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + G_dxT, G_dxN, G_dxE, G_dxU + + integer (kind=int_kind) :: & + i, j , & + im1, im2, jm1, jm2, im3, jm3 , & ! i & j for mom supergrid + ierr + + character(len=*), parameter :: subname = '(mom_dx)' + + if (my_task == master_task) then + allocate( & + G_dxT(nx_global,ny_global), & + G_dxN(nx_global,ny_global), & + G_dxE(nx_global,ny_global), & + G_dxU(nx_global,ny_global), & + stat=ierr & + ) + else + allocate(G_dxT(1,1), G_dxE(1,1), G_dxU(1,1), G_dxN(1,1), stat=ierr) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + if (my_task == master_task) then + work_mom(:,:) = work_mom(:,:) * m_to_cm ! convert to cm + + im1 = 1 ; im2 = 2 ! left ; right coloum of first t-cell + im3 = 3 ! left column of second T-cell, (ie right column of U-cell) + do i = 1, nx_global - 1 + jm1 = 2 ; jm2 = 3 ! middle , top of first row + do j = 1, ny_global + G_dxT(i,j) = work_mom(im1, jm1) + work_mom(im2, jm1) !dxT + G_dxN(i,j) = work_mom(im1, jm2) + work_mom(im2, jm2) !dxN + G_dxE(i,j) = work_mom(im2, jm1) + work_mom(im3, jm1) !dxE + G_dxU(i,j) = work_mom(im2, jm2) + work_mom(im3, jm2) !dxU + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + im1 = im1 + 2 ; im2 = im2 + 2 ; im3 = im3 + 2 + enddo + + ! fill the last col + jm1 = 2 ; jm2 = 3 ! middle , top of first row + do j = 1, ny_global + G_dxT(nx_global,j) = work_mom(2*nx_global - 1, jm1) + work_mom(2*nx_global, jm1) !dxT + G_dxN(nx_global,j) = work_mom(2*nx_global - 1, jm2) + work_mom(2*nx_global, jm2) !dxN + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + jm1 = 2 ; jm2 = 3 ! middle , top of first row + if (trim(ew_boundary_type) == 'cyclic') then + jm1 = 2 ; jm2 = 3 ! middle , top of first row + do j = 1, ny_global + G_dxE(nx_global,j) = work_mom(2*nx_global, jm1) + work_mom(1, jm1) !dxE + G_dxU(nx_global,j) = work_mom(2*nx_global, jm2) + work_mom(1, jm2) !dxU + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + else if (trim(ew_boundary_type) == 'open') then + do j = 1, ny_global + G_dxE(nx_global,j) = 4*work_mom(2*nx_global, jm1) - 2*work_mom(2*nx_global-1, jm1) !dxE + G_dxU(nx_global,j) = 4*work_mom(2*nx_global, jm2) - 2*work_mom(2*nx_global-1, jm2) !dxU + jm1 = jm1 + 2 ; jm2 = jm2 + 2 + enddo + endif + + if (save_ghte_ghtn) then + do j = 1, ny_global + do i = 1, nx_global + G_HTN(i+nghost,j+nghost) = G_dxN(i,j) + enddo + enddo + call global_ext_halo(G_HTN) + endif + endif + + call scatter_global(dxT, G_dxT, master_task, distrb_info, & + field_loc_center, field_type_scalar) + call scatter_global(HTN, G_dxN, master_task, distrb_info, & + field_loc_Nface, field_type_scalar) + dxN(:,:,:) = HTN(:,:,:) + call scatter_global(dxE, G_dxE, master_task, distrb_info, & + field_loc_center, field_type_scalar) + call scatter_global(dxU, G_dxU, master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + + deallocate(G_dxT, G_dxE, G_dxU, G_dxN, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + end subroutine mom_dx + + + subroutine mom_dy(work_mom) + + ! mom sueprgrid has four cells for every model cell, + ! we just need to sum the correct sidelengths to get dy + + real (kind=dbl_kind), dimension(:,:) :: work_mom + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + G_dyT, G_dyN, G_dyE, G_dyU + + integer (kind=int_kind) :: & + i, j, & + im1, im2, jm1, jm2, im3, jm3 , & ! i & j for mom supergrid + ierr + + character(len=*), parameter :: subname = '(mom_dy)' + + if (my_task == master_task) then + allocate( & + G_dyT(nx_global,ny_global), & + G_dyN(nx_global,ny_global), & + G_dyE(nx_global,ny_global), & + G_dyU(nx_global,ny_global), & + stat=ierr & + ) + else + allocate(G_dyT(1,1), G_dyE(1,1), G_dyU(1,1), G_dyN(1,1), stat=ierr ) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + if (my_task == master_task) then + work_mom(:,:) = work_mom(:,:) * m_to_cm ! convert to cm + + im1 = 2 ; im2 = 3 ! middle , right edge of first T-cell + do i = 1, nx_global + jm1 = 1 ; jm2 = 2 ; jm3 = 3 + do j = 1, ny_global - 1 + G_dyT(i,j) = work_mom(im1, jm1) + work_mom(im1, jm2) !dyT + G_dyN(i,j) = work_mom(im1, jm2) + work_mom(im1, jm3) !dyN + G_dyE(i,j) = work_mom(im2, jm1) + work_mom(im2, jm2) !dyE + G_dyU(i,j) = work_mom(im2, jm2) + work_mom(im2, jm3) !dyU + jm1 = jm1 + 2 ; jm2 = jm2 + 2 ; jm3 = jm3 + 2 + enddo + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + + ! fill the top row + im1 = 2 ; im2 = 3 ! middle , right edge of first column + do i = 1, nx_global + G_dyT(i,ny_global) = work_mom(im1, 2*ny_global - 1) + work_mom(im1, 2*ny_global) !dyT + G_dyE(i,ny_global) = work_mom(im2, 2*ny_global - 1) + work_mom(im2, 2*ny_global) !dyE + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + im1 = 2 ; im2 = 3 + if (trim(ns_boundary_type) == 'tripole') then + do i = 1, nx_global + G_dyN(i,ny_global) = work_mom(im1, 2*ny_global) + work_mom(2*nx_global+2-im1, 2*ny_global) !dyN + G_dyU(i,ny_global) = work_mom(im2, 2*ny_global) + work_mom(2*nx_global+2-im2, 2*ny_global) !dyU + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + else if (trim(ns_boundary_type) == 'cyclic') then + do i = 1, nx_global + G_dyN(i,ny_global) = work_mom(im1, 2*ny_global) + work_mom(im1, 1) !dyN + G_dyU(i,ny_global) = work_mom(im2, 2*ny_global) + work_mom(im2, 1) !dyU + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + else if (trim(ns_boundary_type) == 'open') then + do i = 1, nx_global + G_dyN(i,ny_global) = 4*work_mom(im1, 2*ny_global) - 2*work_mom(im1, 2*ny_global-1) !dyN + G_dyU(i,ny_global) = 4*work_mom(im2, 2*ny_global) - 2*work_mom(im2, 2*ny_global-1) !dyU + im1 = im1 + 2 ; im2 = im2 + 2 + enddo + endif + + if (save_ghte_ghtn) then + do j = 1, ny_global + do i = 1, nx_global + G_HTE(i+nghost,j+nghost) = G_dyE(i,j) + enddo + enddo + call global_ext_halo(G_HTE) + endif + endif + + call scatter_global(dyT, G_dyT, master_task, distrb_info, & + field_loc_center, field_type_scalar) + call scatter_global(dyN, G_dyN, master_task, distrb_info, & + field_loc_Nface, field_type_scalar) + call scatter_global(HTE, G_dyE, master_task, distrb_info, & + field_loc_center, field_type_scalar) + dyE(:,:,:) = HTE(:,:,:) + call scatter_global(dyU, G_dyU, master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + + deallocate(G_dyT, G_dyN, G_dyE, G_dyU) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + end subroutine mom_dy + + + subroutine mom_area(work_mom) + + ! mom sueprgrid has four cells for every model cell, we need to sum these + ! to get model cell areas + ! however, earea and narea are calculated from dx & dy - see https://github.com/NOAA-GFDL/MOM6/issues/740 + + real (kind=dbl_kind), dimension(:,:), intent(in) :: work_mom + + integer (kind=int_kind) :: & + i, j, iblk, & + im1, im2, jm1, jm2, im3, jm3 , & ! i & j for mom supergrid + ilo,ihi,jlo,jhi , & ! beginning and end of physical domain + ierr + + type (block) :: & + this_block ! block information for current block + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + G_tarea, G_uarea + + character(len=*), parameter :: subname = '(mom_area)' + + ! calculate narea and earea + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + 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 = 1,ny_block + do i = 1,nx_block + narea(i,j,iblk) = dxN(i,j,iblk)*dyN(i,j,iblk)*cm_to_m*cm_to_m + earea(i,j,iblk) = dxE(i,j,iblk)*dyE(i,j,iblk)*cm_to_m*cm_to_m + enddo + enddo + enddo + + if (my_task == master_task) then + allocate( & + G_tarea(nx_global,ny_global), & + G_uarea(nx_global,ny_global), & + stat=ierr & + ) + else + allocate(G_tarea(1,1), G_uarea(1,1), stat=ierr ) + endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) + + ! load tarea and uarea + if (my_task == master_task) then + im1 = 1 ; im2 = 2 ! left/right -half of first column + im3 = 3 ! right of first U - cell + do i = 1, nx_global - 1 + jm1 = 1 ; jm2 = 2 ! bottom/top -half of first row + jm3 = 3 ! top of first U - cell + do j = 1, ny_global - 1 + G_tarea(i,j) = work_mom(im1, jm1) + work_mom(im1, jm2) & + + work_mom(im2, jm1) + work_mom(im2, jm2) + G_uarea(i,j) = work_mom(im2, jm2) + work_mom(im2, jm3) & + + work_mom(im3, jm2) + work_mom(im3, jm3) + jm1 = jm1 + 2 ; jm2 = jm2 + 2 ; jm3 = jm3 + 2 + enddo + im1 = im1 + 2 ; im2 = im2 + 2 ; im3 = im3 + 2 + enddo + + ! fill last column + jm1 = 1 ; jm2 = 2 ; jm3 = 3 + im1 = 2*nx_global - 1 ; im2 = 2*nx_global ; im3 = 1 + do j = 1, ny_global - 1 + G_tarea(nx_global,j) = work_mom(im1, jm1) + work_mom(im1, jm2) & + + work_mom(im2, jm1) + work_mom(im2, jm2) + if (trim(ew_boundary_type) == 'cyclic') then + G_uarea(nx_global,j) = work_mom(im2, jm2) + work_mom(im2, jm3) & + + work_mom(im3, jm2) + work_mom(im3, jm3) + else if (trim(ew_boundary_type) == 'open') then + G_uarea(nx_global,j) = 4*work_mom(im2, jm2) + 4*work_mom(im2, jm3) & + - 2*work_mom(im1, jm2) - 2*work_mom(im1, jm3) + endif + jm1 = jm1 + 2 ; jm2 = jm2 + 2 ; jm3 = jm3 + 2 + enddo + + ! fill last row + jm1 = ny_global*2 - 1 ; jm2 = ny_global*2 + im1 = 1 ; im2 = 2 ; im3 = 3 + do i = 1, nx_global -1 + G_tarea(i,ny_global) = work_mom(im1, jm1) + work_mom(im1, jm2) & + + work_mom(im2, jm1) + work_mom(im2, jm2) + if (trim(ns_boundary_type) == 'tripole') then + G_uarea(i,ny_global) = work_mom(im2, jm2) + work_mom(2*nx_global+1-im2, jm2) & + + work_mom(im3, jm2) + work_mom(2*nx_global+1-im3, jm2) + else if (trim(ns_boundary_type) == 'cyclic') then + G_uarea(i,ny_global) = work_mom(im2, jm2) + work_mom(im2, jm3) & + + work_mom(im3, jm2) + work_mom(im3, jm3) + else if (trim(ns_boundary_type) == 'open') then + G_uarea(i,ny_global) = 4*work_mom(im2, jm2) + 4*work_mom(im3, jm2) & + - 2*work_mom(im2, jm1) - 2*work_mom(im3, jm1) + endif + im1 = im1 + 2 ; im2 = im2 + 2 ; im3 = im3 + 2 + enddo + + ! the top right corner + im1 = nx_global*2-1 ; im2 = nx_global*2 + jm1 = ny_global*2-1 ; jm2 = ny_global*2 + G_tarea(nx_global,ny_global) = work_mom(im1, jm1) + work_mom(im1, jm2) & + + work_mom(im2, jm1) + work_mom(im2, jm2) + if (trim(ns_boundary_type) == 'tripole') then + G_uarea(nx_global,ny_global) = 2*(work_mom(im2, jm2) + work_mom(1, jm2)) + else if (trim(ns_boundary_type) == 'cyclic' & + .and. trim(ew_boundary_type) == 'cyclic') then + G_uarea(nx_global,ny_global) = work_mom(im2, jm2) + work_mom(1, jm2) & + + work_mom(im2, 1) + work_mom(1, 1) + else if (trim(ns_boundary_type) == 'cyclic' & + .and. trim(ew_boundary_type) == 'open') then + G_uarea(nx_global,ny_global) = 4*work_mom(im2, jm2) + 4*work_mom(im2, 1) & + - 2*work_mom(im1, jm2) - 2*work_mom(im1, 1) + else if (trim(ns_boundary_type) == 'open' & + .and. trim(ew_boundary_type) == 'cyclic') then + G_uarea(nx_global,ny_global) = 4*work_mom(im2, jm2) + 4*work_mom(1, jm2) & + - 2*work_mom(im2, jm1) - 2*work_mom(1, jm1) + else if (trim(ns_boundary_type) == 'open' & + .and. trim(ew_boundary_type) == 'open') then + G_uarea(nx_global,ny_global) = 8*work_mom(im2, jm2) & + - 2*work_mom(im2, jm1) - 2*work_mom(im1, jm2) + endif + + endif + + call scatter_global(tarea, G_tarea, master_task, distrb_info, & + field_loc_center, field_type_scalar) + call scatter_global(uarea, G_uarea, master_task, distrb_info, & + field_loc_NEcorner, field_type_scalar) + deallocate(G_tarea, G_uarea, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) + + end subroutine mom_area + + + subroutine grid_rotation_angle(lon_cnr, lat_cnr, lon_cen, angle) + ! create angles in the same way mom6 creates the angle + ! based on https://github.com/mom-ocean/MOM6/blob/129e1bda02d454fb280819d1d87ae16347fd044c/src/initialization/MOM_shared_initialization.F90#L535 + ! the angle is between logical north on the grid and true north. + + ! global lat/lons/angles + real (kind=dbl_kind), dimension(:,:), intent(in) :: & + lon_cnr, & ! array of lon corner points + lat_cnr, & ! array of lat corner points + lon_cen ! array of lon centre points (i.e. the location the angle is calculated for) + real (kind=dbl_kind), dimension(:,:), intent(out) :: angle + + ! local vars + real (kind=dbl_kind) :: & + lon_scale, & ! The trigonometric scaling factor converting changes in longitude to equivalent distances in latitudes [nondim] + len_lon, & + lon_adj, & + lonB(2,2) + integer (kind=int_kind) :: i, j, m, n + + character(len=*), parameter :: subname = '(grid_rotation_angle)' + + if (my_task == master_task) then + len_lon = maxval(lon_cnr)-minval(lon_cnr) ! The periodic range of longitudes, usually 2pi. + + do j=1,ny_global + do i=1,nx_global + lon_adj = lon_cen(i,j)-p5*len_lon + do n=1,2 ; do m=1,2 + ! shift 4 lon corner points to be similar range to centre point + ! e.g. upper limit of 0 might be shifted to 2*pi + lonB(m,n) = modulo(lon_cnr(i+m-1,j+n-1)-lon_adj, len_lon) & + + lon_adj + enddo ; enddo + lon_scale = cos(p25*(lat_cnr(I,J) + lat_cnr(I+1,J+1) + lat_cnr(I+1,J) + lat_cnr(I,J+1))) + angle(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1) + lonB(2,2) - lonB(1,1))), & + (lat_cnr(I,J+1) - lat_cnr(I+1,J) + lat_cnr(I+1,J+1) - lat_cnr(I,J)) ) + enddo + enddo + endif + + end subroutine grid_rotation_angle + !======================================================================= ! Regular rectangular grid and mask @@ -1846,7 +2646,8 @@ subroutine primary_grid_lengths_HTN(work_g) integer (kind=int_kind) :: & i, j, & - ip1 ! i+1 + ip1 , & ! i+1 + ierr real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g2 @@ -1854,10 +2655,11 @@ subroutine primary_grid_lengths_HTN(work_g) character(len=*), parameter :: subname = '(primary_grid_lengths_HTN)' if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global)) + allocate(work_g2(nx_global,ny_global), stat=ierr) else - allocate(work_g2(1,1)) + allocate(work_g2(1,1), stat=ierr) endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) ! HTN, dxU = average of 2 neighbor HTNs in i @@ -1932,7 +2734,8 @@ subroutine primary_grid_lengths_HTN(work_g) call scatter_global(dxE, work_g2, master_task, distrb_info, & field_loc_center, field_type_scalar) - deallocate(work_g2) + deallocate(work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) end subroutine primary_grid_lengths_HTN @@ -1951,7 +2754,8 @@ subroutine primary_grid_lengths_HTE(work_g) integer (kind=int_kind) :: & i, j, & - im1 ! i-1 + im1, & ! i-1 + ierr real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g2 @@ -1959,10 +2763,11 @@ subroutine primary_grid_lengths_HTE(work_g) character(len=*), parameter :: subname = '(primary_grid_lengths_HTE)' if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global)) + allocate(work_g2(nx_global,ny_global), stat=ierr) else - allocate(work_g2(1,1)) + allocate(work_g2(1,1), stat=ierr) endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) ! HTE, dyU = average of 2 neighbor HTE in j @@ -2041,7 +2846,8 @@ subroutine primary_grid_lengths_HTE(work_g) dyE(:,:,:) = HTE(:,:,:) - deallocate(work_g2) + deallocate(work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc errro', file=__FILE__, line=__LINE__) end subroutine primary_grid_lengths_HTE @@ -2100,7 +2906,8 @@ subroutine makemask integer (kind=int_kind) :: & i, j, iblk, & - ilo,ihi,jlo,jhi ! beginning and end of physical domain + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + ierr real (kind=dbl_kind) :: & puny @@ -2130,7 +2937,8 @@ subroutine makemask !----------------------------------------------------------------- bm = c0 - allocate(uvmCD(nx_block,ny_block,max_blocks)) + allocate(uvmCD(nx_block,ny_block,max_blocks), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) do iblk = 1, nblocks @@ -2235,7 +3043,8 @@ subroutine makemask enddo ! iblk !$OMP END PARALLEL DO - deallocate(uvmCD) + deallocate(uvmCD, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc errro', file=__FILE__, line=__LINE__) end subroutine makemask @@ -3790,8 +4599,9 @@ end function grid_neighbor_max subroutine gridbox_corners integer (kind=int_kind) :: & - i,j,iblk,icorner,& ! index counters - ilo,ihi,jlo,jhi ! beginning and end of physical domain + i,j,iblk,icorner,& ! index counters + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + ierr real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g2 @@ -3851,10 +4661,11 @@ subroutine gridbox_corners !---------------------------------------------------------------- if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global)) + allocate(work_g2(nx_global,ny_global), stat=ierr) else - allocate(work_g2(1,1)) + allocate(work_g2(1,1), stat=ierr) endif + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) work1(:,:,:) = latu_bounds(2,:,:,:) ! work_g2 = c0 @@ -3944,13 +4755,15 @@ subroutine gridbox_corners field_loc_NEcorner, field_type_scalar) lonu_bounds(4,:,:,:) = work1(:,:,:) - deallocate(work_g2) + deallocate(work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc errro', file=__FILE__, line=__LINE__) !---------------------------------------------------------------- ! Convert longitude to Degrees East >0 for history output !---------------------------------------------------------------- - allocate(work_g2(nx_block,ny_block)) ! not used as global here + allocate(work_g2(nx_block,ny_block), stat=ierr) ! not used as global here + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) !OMP fails in this loop do iblk = 1, nblocks do icorner = 1, 4 @@ -3964,7 +4777,8 @@ subroutine gridbox_corners lonu_bounds(icorner,:,:,iblk) = work_g2(:,:) enddo enddo - deallocate(work_g2) + deallocate(work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) end subroutine gridbox_corners @@ -3981,8 +4795,9 @@ end subroutine gridbox_corners subroutine gridbox_edges integer (kind=int_kind) :: & - i,j,iblk,icorner,& ! index counters - ilo,ihi,jlo,jhi ! beginning and end of physical domain + i,j,iblk,icorner,& ! index counters + ilo,ihi,jlo,jhi , & ! beginning and end of physical domain + ierr real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g2 @@ -4055,10 +4870,11 @@ subroutine gridbox_edges !---------------------------------------------------------------- if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global)) + allocate(work_g2(nx_global,ny_global), stat=ierr) else - allocate(work_g2(1,1)) + allocate(work_g2(1,1), stat=ierr) endif + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) ! latn_bounds @@ -4240,13 +5056,15 @@ subroutine gridbox_edges field_loc_NEcorner, field_type_scalar) lone_bounds(3,:,:,:) = work1(:,:,:) - deallocate(work_g2) + deallocate(work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) !---------------------------------------------------------------- ! Convert longitude to Degrees East >0 for history output !---------------------------------------------------------------- - allocate(work_g2(nx_block,ny_block)) ! not used as global here + allocate(work_g2(nx_block,ny_block), stat=ierr) ! not used as global here + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) !OMP fails in this loop do iblk = 1, nblocks do icorner = 1, 4 @@ -4260,7 +5078,7 @@ subroutine gridbox_edges lone_bounds(icorner,:,:,iblk) = work_g2(:,:) enddo enddo - deallocate(work_g2) + deallocate(work_g2, stat=ierr) end subroutine gridbox_edges @@ -4282,7 +5100,8 @@ subroutine gridbox_verts(work_g,vbounds) vbounds integer (kind=int_kind) :: & - i,j ! index counters + i,j , & ! index counters + ierr real (kind=dbl_kind) :: & rad_to_deg @@ -4301,10 +5120,11 @@ subroutine gridbox_verts(work_g,vbounds) file=__FILE__, line=__LINE__) if (my_task == master_task) then - allocate(work_g2(nx_global,ny_global)) + allocate(work_g2(nx_global,ny_global), stat=ierr) else - allocate(work_g2(1,1)) + 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: @@ -4378,7 +5198,8 @@ subroutine gridbox_verts(work_g,vbounds) field_loc_NEcorner, field_type_scalar) vbounds(4,:,:,:) = work1(:,:,:) - deallocate (work_g2) + deallocate (work_g2, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) end subroutine gridbox_verts @@ -4495,7 +5316,8 @@ subroutine get_bathymetry_popfile write(nu_diag,*) subname,' KMT max = ',nlevel endif - allocate(depth(nlevel),thick(nlevel)) + allocate(depth(nlevel),thick(nlevel), stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__) thick = -999999. depth = -999999. @@ -4564,7 +5386,8 @@ subroutine get_bathymetry_popfile if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - deallocate(depth,thick) + deallocate(depth,thick, stat=ierr) + if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) end subroutine get_bathymetry_popfile diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 9bfd3a65b..76c334642 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -307,6 +307,7 @@ grid_nml "``grid_file``", "string", "name of grid file to be read", "'unknown_grid_file'" "``grid_format``", "``bin``", "read direct access grid and kmt files", "``bin``" "", "``nc``", "read grid and kmt files", "" + "", "``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``" "", "``C``", "use C grid structure with T at center, U at E edge, V at N edge", "" "``grid_ocn``", "``A``", "ocn forcing/coupling grid, all fields on T grid", "``A``" diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 0e03c2f13..0dc3156ef 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -133,8 +133,10 @@ The user has several ways to initialize the grid, which can be read from files or created internally. The *rectgrid* code creates a regular rectangular grid (use the namelist option ``grid_type='rectangular'``). The *popgrid* and *popgrid_nc* code reads grid lengths and other parameters for a nonuniform grid (including tripole -and regional grids). -The input files **grid_gx3.bin** and **kmt_gx3.bin** contain the +and regional grids). The *mom_mosaic* code reads grids for a nonuniform grid defined +in the mom mosaic (supergrid) format. + +For the *popgrid* formats, the input files **grid_gx3.bin** and **kmt_gx3.bin** contain the :math:`\left<3^\circ\right>` POP grid and land mask; **grid_gx1.bin** and **kmt_gx1.bin** contain the :math:`\left<1^\circ\right>` grid and land mask, and **grid_tx1.bin** @@ -149,6 +151,15 @@ and ANGLE variables. If the variables exist, ANGLET, TLON and TLAT will also be read from a netcdf grid file. The kmt (mask) netcdf file needs a variable named kmt or mask, set to 0 for land and 1 for ocean. +For the *mom_nc* formats (``grid_format='mom_nc'``), a grid file following +the MOM supergrid convention as generated by `FRE-NCTools `_ +and described in `this short summary `_ +is required. CICE derives the required grid variables to match the grid implementation +in MOM6. All lat/lon and grid lengths are read from the file. Cell areas for A & B grid +cells are loaded from the file, and cell areas for C grid cells are calculated internally. +Rotation angles for cells are also calculated interally. All calculations match the methods +used in MOM6 for consistency. + The input grid file for the B-grid and CD-grid is identical. That file contains each cells' HTN, HTE, ULON, ULAT, and kmt value. From those variables, the longitude, latitude, grid lengths (dx and dy), areas, diff --git a/icepack b/icepack index 379252056..05ac0ec3e 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 3792520561cf9419082ef41f9f0dffd03edf2e43 +Subproject commit 05ac0ec3ea666080eed36e67f6cf8ce1255b243f