Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

retrieve updates #9

Merged
merged 2 commits into from
Jun 24, 2020
Merged

retrieve updates #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f2dc42c
Optionally rotate wind/stress from to computational grid (from true n…
mhrib Jun 23, 2020
3fedc78
Allow for read of tlat, tlon, anglet with popgrid (#463)
TillRasmussen Jun 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions cicecore/cicedynB/dynamics/ice_transport_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -609,8 +609,8 @@ subroutine transport_remap (dt)
asum_init(0), asum_final(0))

if (l_stop) then
write (nu_diag,*) 'istep1, my_task, iblk =', &
istep1, my_task, iblk
write (nu_diag,*) 'istep1, my_task =', &
istep1, my_task
write (nu_diag,*) 'transport: conservation error, cat 0'
call abort_ice(subname//'ERROR: conservation error1')
endif
Expand All @@ -623,8 +623,8 @@ subroutine transport_remap (dt)
atsum_init(:,n), atsum_final(:,n))

if (l_stop) then
write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
istep1, my_task, iblk, n
write (nu_diag,*) 'istep1, my_task, cat =', &
istep1, my_task, n
write (nu_diag,*) 'transport: conservation error, cat ',n
call abort_ice(subname//'ERROR: conservation error2')
endif
Expand Down
82 changes: 46 additions & 36 deletions cicecore/cicedynB/general/ice_forcing.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,9 @@ module ice_forcing
! 'hadgem_sst' or 'hadgem_sst_uvocn'
ice_data_type, & ! 'default', 'box2001', 'boxslotcyl'
precip_units ! 'mm_per_month', 'mm_per_sec', 'mks','m_per_sec'

logical (kind=log_kind), public :: &
rotate_wind ! rotate wind/stress to computational grid from true north directed

character(char_len_long), public :: &
atm_data_dir , & ! top directory for atmospheric data
Expand Down Expand Up @@ -1628,11 +1631,10 @@ subroutine prepare_forcing (nx_block, ny_block, &

if (calc_strair) then

do j = jlo, jhi
do i = ilo, ihi

wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)

if (rotate_wind) then
do j = jlo, jhi
do i = ilo, ihi
wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
!-----------------------------------------------------------------
! Rotate zonal/meridional vectors to local coordinates.
! Velocity comes in on T grid, but is oriented geographically ---
Expand All @@ -1644,30 +1646,38 @@ subroutine prepare_forcing (nx_block, ny_block, &
! atmo_boundary_layer, and are interpolated to the U grid later as
! necessary.
!-----------------------------------------------------------------
workx = uatm(i,j) ! wind velocity, m/s
worky = vatm(i,j)
uatm (i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+ worky*sin(ANGLET(i,j)) ! note uatm, vatm, wind
vatm (i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
- workx*sin(ANGLET(i,j))

enddo ! i
enddo ! j

else ! strax, stray, wind are read from files

do j = jlo, jhi
do i = ilo, ihi

workx = strax(i,j) ! wind stress
worky = stray(i,j)
strax(i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+ worky*sin(ANGLET(i,j)) ! note strax, stray, wind
stray(i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
- workx*sin(ANGLET(i,j))

enddo ! i
enddo ! j
workx = uatm(i,j) ! wind velocity, m/s
worky = vatm(i,j)
uatm (i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+ worky*sin(ANGLET(i,j)) ! note uatm, vatm, wind
vatm (i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
- workx*sin(ANGLET(i,j))
enddo ! i
enddo ! j
else ! not rotated
do j = jlo, jhi
do i = ilo, ihi
wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
enddo ! i
enddo ! j
endif ! rotated

else ! strax, stray, wind are read from files

if (rotate_wind) then
do j = jlo, jhi
do i = ilo, ihi
workx = strax(i,j) ! wind stress
worky = stray(i,j)
strax(i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+ worky*sin(ANGLET(i,j)) ! note strax, stray, wind
stray(i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
- workx*sin(ANGLET(i,j))
enddo ! i
enddo ! j
else ! not rotated
! wind (speed) is already read from file, so all is in place
endif ! rotated

endif ! calc_strair

Expand Down Expand Up @@ -2050,11 +2060,11 @@ subroutine JRA55_gx1_files(yr)
uwind_file = &
trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc'
call file_year(uwind_file,yr)
if (my_task == master_task) then
if (my_task == master_task) then
write (nu_diag,*) ' '
write (nu_diag,*) 'Atmospheric data files:'
write (nu_diag,*) trim(uwind_file)
endif
endif
end subroutine JRA55_gx1_files
subroutine JRA55_tx1_files(yr)
!
Expand All @@ -2066,11 +2076,11 @@ subroutine JRA55_tx1_files(yr)
uwind_file = &
trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_tx1_2005.nc'
call file_year(uwind_file,yr)
if (my_task == master_task) then
if (my_task == master_task) then
write (nu_diag,*) ' '
write (nu_diag,*) 'Atmospheric data files:'
write (nu_diag,*) trim(uwind_file)
endif
endif
end subroutine JRA55_tx1_files
subroutine JRA55_gx3_files(yr)
!
Expand All @@ -2082,11 +2092,11 @@ subroutine JRA55_gx3_files(yr)
uwind_file = &
trim(atm_data_dir)//'/8XDAILY/JRA55_gx3_03hr_forcing_2005.nc'
call file_year(uwind_file,yr)
if (my_task == master_task) then
if (my_task == master_task) then
write (nu_diag,*) ' '
write (nu_diag,*) 'Atmospheric data files:'
write (nu_diag,*) trim(uwind_file)
endif
endif
end subroutine JRA55_gx3_files
!=======================================================================
!
Expand Down Expand Up @@ -4471,7 +4481,7 @@ subroutine hycom_atm_data
write (nu_diag,*) &
'ERROR: CICE: Atm forcing not available at hcdate =',hcdate
write (nu_diag,*) &
'ERROR: CICE: nyr, year_init, yday = ',nyr, year_init, yday
'ERROR: CICE: nyr, year_init, yday ,sec = ',nyr, year_init, yday, sec
call abort_ice ('ERROR: CICE stopped')
endif

Expand Down
9 changes: 6 additions & 3 deletions cicecore/cicedynB/general/ice_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ subroutine input_data
use ice_flux_bgc, only: cpl_bgc
use ice_forcing, only: &
ycycle, fyear_init, dbug, &
atm_data_type, atm_data_dir, precip_units, &
atm_data_type, atm_data_dir, precip_units, rotate_wind, &
atm_data_format, ocn_data_format, &
bgc_data_type, &
ocn_data_type, ocn_data_dir, wave_spec_file, &
Expand Down Expand Up @@ -206,15 +206,15 @@ subroutine input_data
namelist /forcing_nml/ &
formdrag, atmbndy, calc_strair, calc_Tsfc, &
highfreq, natmiter, atmiter_conv, &
ustar_min, emissivity, &
ustar_min, emissivity, &
fbot_xfer_type, update_ocn_f, l_mpond_fresh, tfrz_option, &
oceanmixed_ice, restore_ice, restore_ocn, trestore, &
precip_units, default_season, wave_spec_type,nfreq, &
atm_data_type, ocn_data_type, bgc_data_type, fe_data_type, &
ice_data_type, wave_spec_file, &
fyear_init, ycycle, &
atm_data_dir, ocn_data_dir, bgc_data_dir, &
atm_data_format, ocn_data_format, &
atm_data_format, ocn_data_format, rotate_wind, &
oceanmixed_file

!-----------------------------------------------------------------
Expand Down Expand Up @@ -357,6 +357,7 @@ subroutine input_data
atm_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
atm_data_type = 'default'
atm_data_dir = ' '
rotate_wind = .true. ! rotate wind/stress composants to computational grid orientation
calc_strair = .true. ! calculate wind stress
formdrag = .false. ! calculate form drag
highfreq = .false. ! calculate high frequency RASM coupling
Expand Down Expand Up @@ -631,6 +632,7 @@ subroutine input_data
call broadcast_scalar(atm_data_format, master_task)
call broadcast_scalar(atm_data_type, master_task)
call broadcast_scalar(atm_data_dir, master_task)
call broadcast_scalar(rotate_wind, master_task)
call broadcast_scalar(calc_strair, master_task)
call broadcast_scalar(calc_Tsfc, master_task)
call broadcast_scalar(formdrag, master_task)
Expand Down Expand Up @@ -1279,6 +1281,7 @@ subroutine input_data
write(nu_diag,*) '--------------------------------'
write(nu_diag,1012) ' calc_Tsfc = ', calc_Tsfc,' calculate surface temperature as part of thermo'
write(nu_diag,1012) ' calc_strair = ', calc_strair,' calculate wind stress and speed'
write(nu_diag,1012) ' rotate_wind = ', rotate_wind,' rotate wind/stress to computational grid'
write(nu_diag,1012) ' formdrag = ', formdrag,' use form drag parameterization'
if (trim(atmbndy) == 'constant') then
tmpstr2 = ': stability-based boundary layer'
Expand Down
6 changes: 4 additions & 2 deletions cicecore/cicedynB/general/ice_step_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,15 +182,17 @@ subroutine step_therm1 (dt, iblk)
logical (kind=log_kind) :: &
prescribed_ice ! if .true., use prescribed ice instead of computed
#endif

real (kind=dbl_kind), intent(in) :: &
dt ! time step

integer (kind=int_kind), intent(in) :: &
iblk ! block index

! local variables

#ifdef CICE_IN_NEMO
real (kind=dbl_kind) :: &
raice ! temporary reverse ice concentration
#endif
integer (kind=int_kind) :: &
ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
i, j , & ! horizontal indices
Expand Down
61 changes: 46 additions & 15 deletions cicecore/cicedynB/infrastructure/ice_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,10 @@ module ice_grid
real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
rndex_global ! global index for local subdomain (dbl)

logical (kind=log_kind), private :: &
l_readCenter ! If anglet exist in grid file read it otherwise calculate it


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

contains
Expand Down Expand Up @@ -332,7 +336,6 @@ subroutine init_grid2
real (kind=dbl_kind) :: &
angle_0, angle_w, angle_s, angle_sw, &
pi, pi2, puny
! real (kind=dbl_kind) :: ANGLET_dum
logical (kind=log_kind), dimension(nx_block,ny_block,max_blocks):: &
out_of_range

Expand Down Expand Up @@ -470,11 +473,10 @@ subroutine init_grid2
!-----------------------------------------------------------------
! Compute ANGLE on T-grid
!-----------------------------------------------------------------
ANGLET = c0

if (trim(grid_type) == 'cpom_grid') then
ANGLET(:,:,:) = ANGLE(:,:,:)
else
else if (.not. (l_readCenter)) then
ANGLET = c0

!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP angle_0,angle_w,angle_s,angle_sw)
Expand Down Expand Up @@ -504,7 +506,8 @@ subroutine init_grid2
enddo
!$OMP END PARALLEL DO
endif ! cpom_grid
if (trim(grid_type) == 'regional') then
if (trim(grid_type) == 'regional' .and. &
(.not. (l_readCenter))) then
! for W boundary extrapolate from interior
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
do iblk = 1, nblocks
Expand All @@ -531,9 +534,9 @@ subroutine init_grid2
call ice_timer_stop(timer_bound)

call makemask ! velocity mask, hemisphere masks

call Tlatlon ! get lat, lon on the T grid

if (.not. (l_readCenter)) then
call Tlatlon ! get lat, lon on the T grid
endif
!-----------------------------------------------------------------
! bathymetry
!-----------------------------------------------------------------
Expand Down Expand Up @@ -716,6 +719,7 @@ subroutine popgrid_nc
field_loc_center, field_loc_NEcorner, &
field_type_scalar, field_type_angle
use ice_domain_size, only: max_blocks
use netcdf

integer (kind=int_kind) :: &
i, j, iblk, &
Expand All @@ -739,6 +743,12 @@ subroutine popgrid_nc

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

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


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

Expand All @@ -751,7 +761,7 @@ subroutine popgrid_nc
call ice_open_nc(kmt_file,fid_kmt)

diag = .true. ! write diagnostic info

l_readCenter = .false.
!-----------------------------------------------------------------
! topography
!-----------------------------------------------------------------
Expand Down Expand Up @@ -806,11 +816,37 @@ subroutine popgrid_nc
call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE
call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
field_loc_NEcorner, field_type_angle)

! 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, trim(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_g1,diag)
call scatter_global(ANGLET, work_g1, master_task, distrb_info, &
field_loc_center, field_type_angle)
where (ANGLET > pi) ANGLET = pi
where (ANGLET < -pi) ANGLET = -pi
fieldname="tlon"
call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
call scatter_global(TLON, work_g1, master_task, distrb_info, &
field_loc_center, field_type_scalar)
fieldname="tlat"
call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
call scatter_global(TLAT, work_g1, master_task, distrb_info, &
field_loc_center, field_type_scalar)
endif
!-----------------------------------------------------------------
! cell dimensions
! calculate derived quantities from global arrays to preserve
Expand All @@ -820,7 +856,6 @@ subroutine popgrid_nc
fieldname='htn'
call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTN
call primary_grid_lengths_HTN(work_g1) ! dxu, dxt

fieldname='hte'
call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTE
call primary_grid_lengths_HTE(work_g1) ! dyu, dyt
Expand All @@ -831,7 +866,6 @@ subroutine popgrid_nc
call ice_close_nc(fid_grid)
call ice_close_nc(fid_kmt)
endif

#endif
end subroutine popgrid_nc

Expand Down Expand Up @@ -1737,7 +1771,6 @@ subroutine Tlatlon
enddo ! j
enddo ! iblk
!$OMP END PARALLEL DO

if (trim(grid_type) == 'regional') then
! for W boundary extrapolate from interior
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
Expand Down Expand Up @@ -2463,7 +2496,6 @@ subroutine read_basalstress_bathy

! use module
use ice_read_write
use ice_communicate, only: my_task, master_task
use ice_constants, only: field_loc_center, field_type_scalar

! local variables
Expand Down Expand Up @@ -2491,7 +2523,6 @@ subroutine read_basalstress_bathy

if (my_task == master_task) then
write(nu_diag,*) 'reading ',TRIM(fieldname)
write(*,*) 'reading ',TRIM(fieldname)
call icepack_warnings_flush(nu_diag)
endif
call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &
Expand Down
Loading