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.

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

Update gx1 initial condition, Icepack, tests, and version. Fix hmix default value #586

Merged
merged 9 commits into from
Apr 6, 2021
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
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
245 changes: 2 additions & 243 deletions cicecore/cicedynB/general/ice_forcing.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2668,247 +2668,6 @@ subroutine JRA55_data

end subroutine JRA55_data

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

subroutine Jra55_data_old (yr)

use ice_blocks, only: block, get_block
use ice_global_reductions, only: global_minval, global_maxval
use ice_domain, only: nblocks, distrb_info, blocks_ice
use ice_flux, only: fsnow, Tair, uatm, vatm, Qa, fsw, flw
use ice_grid, only: hm, tlon, tlat, tmask, umask
use ice_state, only: aice
use ice_calendar, only: days_per_year, use_leap_years

integer (kind=int_kind), intent(in) :: &
yr ! current forcing year

integer (kind=int_kind) :: &
ncid , & ! netcdf file id
i, j, n1, iblk, &
yrp , & ! year after yr in forcing cycle
recnum , & ! record number
maxrec , & ! maximum record number
recslot , & ! spline slot for current record
dataloc ! = 1 for data located in middle of time interval
! = 2 for date located at end of time interval

real (kind=dbl_kind) :: &
sec3hr , & ! number of seconds in 3 hours
secday , & ! number of seconds in day
eps, tt , & ! interpolation coeff calc
Tffresh , &
vmin, vmax

logical (kind=log_kind) :: debug_n_d = .false.

character (char_len_long) :: uwind_file_old
character(len=64) :: fieldname !netcdf field name
character(len=*), parameter :: subname = '(Jra55_data_old)'

if (forcing_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'

call icepack_query_parameters(Tffresh_out=Tffresh)
call icepack_query_parameters(secday_out=secday)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

sec3hr = secday/c8 ! seconds in 3 hours
maxrec = days_per_year*8

if (debug_n_d .and. my_task == master_task) then
write (nu_diag,*) subname,'recnum',recnum
write (nu_diag,*) subname,'maxrec',maxrec
write (nu_diag,*) subname,'days_per_year', days_per_year
endif

!-------------------------------------------------------------------
! 3-hourly data
! states are instantaneous, 1st record is 00z Jan 1
! fluxes are 3 hour averages, 1st record is 00z-03z Jan 1
! Both states and fluxes have 1st record defined as 00z Jan 1
! interpolate states, do not interpolate fluxes
! fluxes are held constant from [init period, end period)
!-------------------------------------------------------------------
! File is NETCDF with winds in NORTH and EAST direction
! file variable names are:
! glbrad (shortwave W/m^2)
! dlwsfc (longwave W/m^2)
! wndewd (eastward wind m/s)
! wndnwd (northward wind m/s)
! airtmp (air temperature K)
! spchmd (specific humidity kg/kg)
! ttlpcp (precipitation kg/m s-1)
!-------------------------------------------------------------------

uwind_file_old = uwind_file
call file_year(uwind_file,yr)
if (uwind_file /= uwind_file_old .and. my_task == master_task) then
write(nu_diag,*) subname,' reading forcing file = ',trim(uwind_file)
endif

call ice_open_nc(uwind_file,ncid)

do n1 = 1,2

if (n1 == 1) then
recnum = 8*int(yday) - 7 + int(real(msec,kind=dbl_kind)/sec3hr)
if (my_task == master_task .and. (recnum <= 2 .or. recnum >= maxrec-1)) then
write(nu_diag,*) subname,' reading forcing file 1st ts = ',trim(uwind_file)
endif
elseif (n1 == 2) then
recnum = 8*int(yday) - 7 + int(real(msec,kind=dbl_kind)/sec3hr) + 1
if (recnum > maxrec) then
yrp = fyear_init + mod(myear,ycycle) ! next year
recnum = 1
call file_year(uwind_file,yrp)
if (my_task == master_task) then
write(nu_diag,*) subname,' reading forcing file 2nd ts = ',trim(uwind_file)
endif
call ice_close_nc(ncid)
call ice_open_nc(uwind_file,ncid)
endif
endif

if (debug_n_d .and. my_task == master_task) then
write(nu_diag,*) subname,' read recnum = ',recnum,n1
endif

fieldname = 'airtmp'
call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

fieldname = 'wndewd'
call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

fieldname = 'wndnwd'
call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

fieldname = 'spchmd'
call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

! only read one timestep for fluxes, 3 hr average, no interpolation
if (n1 == 1) then
fieldname = 'glbrad'
call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

fieldname = 'dlwsfc'
call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

fieldname = 'ttlpcp'
call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,n1,:),debug_n_d, &
field_loc=field_loc_center, &
field_type=field_type_scalar)
endif

enddo

call ice_close_nc(ncid)

! reset uwind_file to original year
call file_year(uwind_file,yr)

! Compute interpolation coefficients
eps = 1.0e-6
tt = real(mod(msec,nint(sec3hr)),kind=dbl_kind)
c2intp = tt / sec3hr
if (c2intp < c0 .and. c2intp > c0-eps) c2intp = c0
if (c2intp > c1 .and. c2intp < c1+eps) c2intp = c1
c1intp = 1.0_dbl_kind - c2intp
if (c2intp < c0 .or. c2intp > c1) then
write(nu_diag,*) subname,' ERROR: c2intp = ',c2intp
call abort_ice (error_message=subname//' ERROR: c2intp out of range', &
file=__FILE__, line=__LINE__)
endif
if (debug_n_d .and. my_task == master_task) then
write(nu_diag,*) subname,' c12intp = ',c1intp,c2intp
endif

! Interpolate
call interpolate_data (Tair_data, Tair)
call interpolate_data (uatm_data, uatm)
call interpolate_data (vatm_data, vatm)
call interpolate_data (Qa_data, Qa)
! use 3 hr average for heat flux and precip fields
! call interpolate_data (fsw_data, fsw)
! call interpolate_data (flw_data, flw)
! call interpolate_data (fsnow_data, fsnow)
fsw(:,:,:) = fsw_data(:,:,1,:)
flw(:,:,:) = flw_data(:,:,1,:)
fsnow(:,:,:) = fsnow_data(:,:,1,:)

!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
! limit summer Tair values where ice is present
do j = 1, ny_block
do i = 1, nx_block
if (aice(i,j,iblk) > p1) Tair(i,j,iblk) = min(Tair(i,j,iblk), Tffresh+p1)
enddo
enddo

do j = 1, ny_block
do i = 1, nx_block
Qa (i,j,iblk) = Qa (i,j,iblk) * hm(i,j,iblk)
Tair(i,j,iblk) = Tair(i,j,iblk) * hm(i,j,iblk)
uatm(i,j,iblk) = uatm(i,j,iblk) * hm(i,j,iblk)
vatm(i,j,iblk) = vatm(i,j,iblk) * hm(i,j,iblk)
fsw (i,j,iblk) = fsw (i,j,iblk) * hm(i,j,iblk)
flw (i,j,iblk) = flw (i,j,iblk) * hm(i,j,iblk)
fsnow(i,j,iblk) = fsnow (i,j,iblk) * hm(i,j,iblk)
enddo
enddo

enddo ! iblk
!$OMP END PARALLEL DO

if (debug_n_d .or. dbug) then
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'JRA55_bulk_data'
vmin = global_minval(fsw,distrb_info,tmask)
vmax = global_maxval(fsw,distrb_info,tmask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'fsw',vmin,vmax
vmin = global_minval(flw,distrb_info,tmask)
vmax = global_maxval(flw,distrb_info,tmask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'flw',vmin,vmax
vmin =global_minval(fsnow,distrb_info,tmask)
vmax =global_maxval(fsnow,distrb_info,tmask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'fsnow',vmin,vmax
vmin = global_minval(Tair,distrb_info,tmask)
vmax = global_maxval(Tair,distrb_info,tmask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'Tair',vmin,vmax
vmin = global_minval(uatm,distrb_info,umask)
vmax = global_maxval(uatm,distrb_info,umask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'uatm',vmin,vmax
vmin = global_minval(vatm,distrb_info,umask)
vmax = global_maxval(vatm,distrb_info,umask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'vatm',vmin,vmax
vmin = global_minval(Qa,distrb_info,tmask)
vmax = global_maxval(Qa,distrb_info,tmask)
if (my_task.eq.master_task) &
write (nu_diag,*) subname,'Qa',vmin,vmax

endif ! dbug

end subroutine Jra55_data_old

!=======================================================================
!
! AOMIP shortwave forcing
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a note, mainly as a reminder: We are making JRA the default, but the older forcing options are still available. Those will need to be deprecated via our process (not in this release).

Expand Down Expand Up @@ -4394,7 +4153,7 @@ subroutine ocn_data_ncar(dt)
do j = 1, ny_block
do i = 1, nx_block
if (n == 2) sss (i,j,:) = c0
if (n == 3) hmix (i,j,:) = c0
if (n == 3) hmix (i,j,:) = c20
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this just temporary because of NaNs in the file? If not, then it would be better to set this as a parameter at the top of the module, to make it more explicit. 'default mixed layer depth'

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this avoids a divide by zero. I assume there are also some science reasons for setting it to 20. I can create a parameter in the code as suggested and will do that later today. I think that's a good idea.

Copy link
Contributor

@TillRasmussen TillRasmussen Apr 6, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a way of avoiding divison by zero. @dabail10 suggested max(c20,work1(i,j,iblk)). I dont have a preference. C20 is set to match the only other fixed setting in subroutine init_coupler_flux.

if (n == 4) uocn (i,j,:) = c0
if (n == 5) vocn (i,j,:) = c0
if (n == 6) ss_tltx(i,j,:) = c0
Expand All @@ -4403,7 +4162,7 @@ subroutine ocn_data_ncar(dt)
do iblk = 1, nblocks
if (hm(i,j,iblk) == c1) then
if (n == 2) sss (i,j,iblk) = work1(i,j,iblk)
if (n == 3) hmix (i,j,iblk) = work1(i,j,iblk)
if (n == 3) hmix (i,j,iblk) = max(c20,work1(i,j,iblk))
if (n == 4) uocn (i,j,iblk) = work1(i,j,iblk)
if (n == 5) vocn (i,j,iblk) = work1(i,j,iblk)
if (n == 6) ss_tltx(i,j,iblk) = work1(i,j,iblk)
Expand Down
2 changes: 1 addition & 1 deletion cicecore/version.txt
Original file line number Diff line number Diff line change
@@ -1 +1 @@
CICE 6.1.4
CICE 6.2.0
3 changes: 3 additions & 0 deletions configuration/scripts/ice_in
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,9 @@
days_per_year = 365
use_leap_years = .false.
year_init = 1997
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't year_init match the new forcing data? maybe it doesn't matter, as long as it's properly set in the -s options. I see below that fyear_init=2005, but I thought that one wasn't used now. It's a little confusing, would be better to be consistent.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I saw this and wondered the same things. I was reluctant to make this change as I didn't understand the importance or history. The other thing is that it will change all the gx3 results as well in our regression test. But I agree 100% that it would be nice if all of this were consistent. Are there other defaults that need to change. Should use_leap_years be always true, for instance? I can make this change as part of this PR, but it will change answers.

month_init = 1
day_init = 1
sec_init = 0
istep0 = 0
dt = 3600.0
npt_unit = '1'
Expand Down
4 changes: 2 additions & 2 deletions configuration/scripts/options/set_nml.gx1
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ runtype = 'initial'
year_init = 2005
use_leap_years = .true.
use_restart_time = .false.
ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc'
ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced.2005-01-01-00000.nc'
grid_format = 'bin'
grid_type = 'displaced_pole'
grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin'
Expand All @@ -17,5 +17,5 @@ atm_data_format = 'nc'
atm_data_type = 'JRA55_gx1'
atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55'
precip_units = 'mks'
ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/COREII'
ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/CESM/MONTHLY'
bgc_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/WOA/MONTHLY'
5 changes: 5 additions & 0 deletions configuration/scripts/options/set_nml.gx1apr
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
year_init = 2005
month_init = 4
day_init = 1
sec_init = 0
ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced.2005-04-01-00000.nc'
1 change: 1 addition & 0 deletions configuration/scripts/options/set_nml.gx1coreii
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
year_init = 1997
use_leap_years = .false.
use_restart_time = .true.
ice_ic = '/glade/p/cesm/pcwg_dev/CICE_data/ic/gx1/iced_gx1_v5.nc'
fyear_init = 2005
ycycle = 1
atm_data_format = 'bin'
Expand Down
21 changes: 15 additions & 6 deletions configuration/scripts/options/set_nml.gx1prod
Original file line number Diff line number Diff line change
@@ -1,7 +1,16 @@
year_init = 1958
dt = 3600
npt = 87600
year_init = 2005
npt_unit = 'y'
npt = 1
dumpfreq = 'm'
fyear_init = 1958
ycycle = 52
ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1'
fyear_init = 2005
ycycle = 5
ocn_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/CESM/MONTHLY/'
use_bathymetry = .true.
seabed_stress = .true.
seabed_stress_method = 'LKD'
ocn_data_type = 'ncar'
ocn_data_format = 'nc'
oceanmixed_file = 'ocean_forcing_clim_2D_gx1.nc'
tr_brine = .true.
f_taubx = 'm'
f_tauby = 'm'
6 changes: 6 additions & 0 deletions configuration/scripts/options/set_nml.seabedLKD
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
use_bathymetry = .true.
seabed_stress = .true.
seabed_stress_method = 'LKD'
histfreq = 'm','d','x','x','x'
f_taubx = 'md'
f_tauby = 'md'
6 changes: 6 additions & 0 deletions configuration/scripts/options/set_nml.seabedprob
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
use_bathymetry = .true.
seabed_stress = .true.
seabed_stress_method = 'probabilistic'
histfreq = 'm','d','x','x','x'
f_taubx = 'md'
f_tauby = 'md'
6 changes: 4 additions & 2 deletions configuration/scripts/tests/base_suite.ts
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,10 @@ restart gx1 8x1 bgczclim,medium
smoke gx1 24x1 medium,run90day,yi2008
smoke gx3 8x1 medium,run90day,yi2008
restart gx1 24x1 short
restart gx1 24x1 short,ml
restart gx3 8x1 short
restart gx1 16x2 seabedLKD,gx1apr,medium,debug
restart gx1 16x2 seabedLKD,gx1apr,medium
restart gx1 15x2 seabedprob,medium
restart gx1 32x1 gx1prod,medium
smoke gx3 4x2 fsd1,diag24,run5day,debug
smoke gx3 8x2 fsd12,diag24,run5day,short
restart gx3 4x2 fsd12,debug,short
Expand Down
4 changes: 2 additions & 2 deletions doc/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,9 @@
# built documents.
#
# The short X.Y version.
version = u'6.1.4'
version = u'6.2.0'
# The full version, including alpha/beta/rc tags.
version = u'6.1.4'
version = u'6.2.0'

# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
Expand Down