Skip to content

Commit

Permalink
Merge pull request NOAA-EMC#8 from ESCOMP/mvertens/nuopc_updates
Browse files Browse the repository at this point in the history 
introduction of mask file for generating mask when appropriate
  • Loading branch information
@dabail10
dabail10 authored Mar 22, 2021
2 parents d78c3a3 + 0f90ce1 commit 1139e7f
Show file tree
Hide file tree
Showing 6 changed files with 1,655 additions and 1,521 deletions.
2 changes: 1 addition & 1 deletion cicecore/cicedynB/general/ice_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1604,7 +1604,7 @@ subroutine input_data
grid_type /= 'rectangular' .and. &
grid_type /= 'cpom_grid' .and. &
grid_type /= 'regional' .and. &
grid_type /= 'latlon' ) then
grid_type /= 'setmask' ) then
if (my_task == master_task) write(nu_diag,*) subname//' ERROR: unknown grid_type=',trim(grid_type)
abort_list = trim(abort_list)//":20"
endif
Expand Down
287 changes: 1 addition & 286 deletions cicecore/cicedynB/infrastructure/ice_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ module ice_grid
private
public :: init_grid1, init_grid2, &
t2ugrid_vector, u2tgrid_vector, &
to_ugrid, to_tgrid, alloc_grid
to_ugrid, to_tgrid, alloc_grid, makemask

character (len=char_len_long), public :: &
grid_format , & ! file format ('bin'=binary or 'nc'=netcdf)
Expand Down Expand Up @@ -368,11 +368,6 @@ subroutine init_grid2
else
call popgrid ! read POP grid lengths directly
endif
#ifdef CESMCOUPLED
elseif (trim(grid_type) == 'latlon') then
call latlongrid ! lat lon grid for sequential CESM (CAM mode)
return
#endif
elseif (trim(grid_type) == 'cpom_grid') then
call cpomgrid ! cpom model orca1 type grid
else
Expand Down Expand Up @@ -882,286 +877,6 @@ subroutine popgrid_nc

end subroutine popgrid_nc

#ifdef CESMCOUPLED
!=======================================================================

! Read in kmt file that matches CAM lat-lon grid and has single column
! functionality
! author: Mariana Vertenstein
! 2007: Elizabeth Hunke upgraded to netcdf90 and cice ncdf calls

subroutine latlongrid

! use ice_boundary
use ice_domain_size
use ice_scam, only : scmlat, scmlon, single_column
use ice_constants, only: c0, c1, p5, p25, &
field_loc_center, field_type_scalar, radius
#ifdef USE_NETCDF
use netcdf
#endif

integer (kind=int_kind) :: &
i, j, iblk

integer (kind=int_kind) :: &
ni, nj, ncid, dimid, varid, ier

type (block) :: &
this_block ! block information for current block

integer (kind=int_kind) :: &
ilo,ihi,jlo,jhi ! beginning and end of physical domain

real (kind=dbl_kind) :: &
closelat, & ! Single-column latitude value
closelon, & ! Single-column longitude value
closelatidx, & ! Single-column latitude index to retrieve
closelonidx ! Single-column longitude index to retrieve

integer (kind=int_kind) :: &
start(2), & ! Start index to read in
count(2) ! Number of points to read in

integer (kind=int_kind) :: &
start3(3), & ! Start index to read in
count3(3) ! Number of points to read in

integer (kind=int_kind) :: &
status ! status flag

real (kind=dbl_kind), allocatable :: &
lats(:),lons(:),pos_lons(:), glob_grid(:,:) ! temporaries

real (kind=dbl_kind) :: &
pos_scmlon,& ! temporary
pi, &
puny, &
scamdata ! temporary

character(len=*), parameter :: subname = '(lonlatgrid)'

#ifdef USE_NETCDF
!-----------------------------------------------------------------
! - kmt file is actually clm fractional land file
! - Determine consistency of dimensions
! - Read in lon/lat centers in degrees from kmt file
! - Read in ocean from "kmt" file (1 for ocean, 0 for land)
!-----------------------------------------------------------------

call icepack_query_parameters(pi_out=pi, puny_out=puny)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

! Determine dimension of domain file and check for consistency

if (my_task == master_task) then
call ice_open_nc(kmt_file, ncid)

status = nf90_inq_dimid (ncid, 'ni', dimid)
status = nf90_inquire_dimension(ncid, dimid, len=ni)
status = nf90_inq_dimid (ncid, 'nj', dimid)
status = nf90_inquire_dimension(ncid, dimid, len=nj)
end if

! Determine start/count to read in for either single column or global lat-lon grid
! If single_column, then assume that only master_task is used since there is only one task

if (single_column) then
! Check for consistency
if (my_task == master_task) then
if ((nx_global /= 1).or. (ny_global /= 1)) then
write(nu_diag,*) 'Because you have selected the column model flag'
write(nu_diag,*) 'Please set nx_global=ny_global=1 in file'
write(nu_diag,*) 'ice_domain_size.F and recompile'
call abort_ice (subname//'ERROR: check nx_global, ny_global')
endif
end if

! Read in domain file for single column
allocate(lats(nj))
allocate(lons(ni))
allocate(pos_lons(ni))
allocate(glob_grid(ni,nj))

start3=(/1,1,1/)
count3=(/ni,nj,1/)
status = nf90_inq_varid(ncid, 'xc' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc')
status = nf90_get_var(ncid, varid, glob_grid, start3, count3)
if (status /= nf90_noerr) call abort_ice (subname//' get_var xc')
do i = 1,ni
lons(i) = glob_grid(i,1)
end do

status = nf90_inq_varid(ncid, 'yc' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc')
status = nf90_get_var(ncid, varid, glob_grid, start3, count3)
if (status /= nf90_noerr) call abort_ice (subname//' get_var yc')
do j = 1,nj
lats(j) = glob_grid(1,j)
end do

! convert lons array and scmlon to 0,360 and find index of value closest to 0
! and obtain single-column longitude/latitude indices to retrieve

pos_lons(:)= mod(lons(:) + 360._dbl_kind,360._dbl_kind)
pos_scmlon = mod(scmlon + 360._dbl_kind,360._dbl_kind)
start(1) = (MINLOC(abs(pos_lons-pos_scmlon),dim=1))
start(2) = (MINLOC(abs(lats -scmlat ),dim=1))

deallocate(lats)
deallocate(lons)
deallocate(pos_lons)
deallocate(glob_grid)

status = nf90_inq_varid(ncid, 'xc' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc')
status = nf90_get_var(ncid, varid, scamdata, start)
if (status /= nf90_noerr) call abort_ice (subname//' get_var xc')
TLON = scamdata
status = nf90_inq_varid(ncid, 'yc' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc')
status = nf90_get_var(ncid, varid, scamdata, start)
if (status /= nf90_noerr) call abort_ice (subname//' get_var yc')
TLAT = scamdata
status = nf90_inq_varid(ncid, 'area' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid area')
status = nf90_get_var(ncid, varid, scamdata, start)
if (status /= nf90_noerr) call abort_ice (subname//' get_var are')
tarea = scamdata
status = nf90_inq_varid(ncid, 'mask' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid mask')
status = nf90_get_var(ncid, varid, scamdata, start)
if (status /= nf90_noerr) call abort_ice (subname//' get_var mask')
hm = scamdata
status = nf90_inq_varid(ncid, 'frac' , varid)
if (status /= nf90_noerr) call abort_ice (subname//' inq_varid frac')
status = nf90_get_var(ncid, varid, scamdata, start)
if (status /= nf90_noerr) call abort_ice (subname//' get_var frac')
ocn_gridcell_frac = scamdata
else
! Check for consistency
if (my_task == master_task) then
if (nx_global /= ni .and. ny_global /= nj) then
write(nu_diag,*) 'latlongrid: ni,nj = ',ni,nj
write(nu_diag,*) 'latlongrid: nx_g,ny_g = ',nx_global, ny_global
call abort_ice (subname//'ERROR: ni,nj not equal to nx_global,ny_global')
end if
end if

! Read in domain file for global lat-lon grid
call ice_read_nc(ncid, 1, 'xc' , TLON , diag=.true.)
call ice_read_nc(ncid, 1, 'yc' , TLAT , diag=.true.)
call ice_read_nc(ncid, 1, 'area', tarea , diag=.true., &
field_loc=field_loc_center,field_type=field_type_scalar)
call ice_read_nc(ncid, 1, 'mask', hm , diag=.true.)
call ice_read_nc(ncid, 1, 'frac', ocn_gridcell_frac, diag=.true.)
end if

if (my_task == master_task) then
call ice_close_nc(ncid)
end if

!$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
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
! Convert from degrees to radians
TLON(i,j,iblk) = pi*TLON(i,j,iblk)/180._dbl_kind

! Convert from degrees to radians
TLAT(i,j,iblk) = pi*TLAT(i,j,iblk)/180._dbl_kind

! Convert from radians^2 to m^2
! (area in domain file is in radians^2 and tarea is in m^2)
tarea(i,j,iblk) = tarea(i,j,iblk) * (radius*radius)
end do
end do
end do
!$OMP END PARALLEL DO

!-----------------------------------------------------------------
! Calculate various geometric 2d arrays
! The U grid (velocity) is not used when run with sequential CAM
! because we only use thermodynamic sea ice. However, ULAT is used
! in the default initialization of CICE so we calculate it here as
! a "dummy" so that CICE will initialize with ice. If a no ice
! initialization is OK (or desired) this can be commented out and
! ULAT will remain 0 as specified above. ULAT is located at the
! NE corner of the grid cell, TLAT at the center, so here ULAT is
! hacked by adding half the latitudinal spacing (in radians) to
! TLAT.
!-----------------------------------------------------------------

!$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
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

if (ny_global == 1) then
uarea(i,j,iblk) = tarea(i,j, iblk)
else
uarea(i,j,iblk) = p25* &
(tarea(i,j, iblk) + tarea(i+1,j, iblk) &
+ tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk))
endif
tarear(i,j,iblk) = c1/tarea(i,j,iblk)
uarear(i,j,iblk) = c1/uarea(i,j,iblk)
tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)

if (single_column) then
ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj)
else
if (ny_global == 1) then
ULAT (i,j,iblk) = TLAT(i,j,iblk)
else
ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global)
endif
endif
ULON (i,j,iblk) = c0
ANGLE (i,j,iblk) = c0

ANGLET(i,j,iblk) = c0
HTN (i,j,iblk) = 1.e36_dbl_kind
HTE (i,j,iblk) = 1.e36_dbl_kind
dxt (i,j,iblk) = 1.e36_dbl_kind
dyt (i,j,iblk) = 1.e36_dbl_kind
dxu (i,j,iblk) = 1.e36_dbl_kind
dyu (i,j,iblk) = 1.e36_dbl_kind
dxhy (i,j,iblk) = 1.e36_dbl_kind
dyhx (i,j,iblk) = 1.e36_dbl_kind
cyp (i,j,iblk) = 1.e36_dbl_kind
cxp (i,j,iblk) = 1.e36_dbl_kind
cym (i,j,iblk) = 1.e36_dbl_kind
cxm (i,j,iblk) = 1.e36_dbl_kind
enddo
enddo
enddo
!$OMP END PARALLEL DO

call makemask
#else
call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', &
file=__FILE__, line=__LINE__)
#endif

end subroutine latlongrid
#endif

!=======================================================================

! Regular rectangular grid and mask
Expand Down
Loading

0 comments on commit 1139e7f

Please sign in to comment.