diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 old mode 100644 new mode 100755 index f1a307ae0..53eb0ec49 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -30,7 +30,7 @@ module ice_forcing use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, & timer_bound use ice_arrays_column, only: oceanmixed_ice, restore_bgc - use ice_constants, only: c0, c1, c2, c3, c4, c5, c10, c12, c15, c20, & + use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c15, c20, & c180, c360, c365, c1000, c3600 use ice_constants, only: p001, p01, p1, p2, p25, p5, p6 use ice_constants, only: cm_to_m @@ -115,7 +115,8 @@ module ice_forcing atm_data_format, & ! 'bin'=binary or 'nc'=netcdf ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf atm_data_type, & ! 'default', 'monthly', 'ncar', - ! 'LYq' or 'hadgem' or 'oned' + ! 'LYq' or 'hadgem' or 'oned' or + ! 'JRA55' bgc_data_type, & ! 'default', 'clim' ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned', ! 'hadgem_sst' or 'hadgem_sst_uvocn' @@ -211,6 +212,8 @@ subroutine init_forcing_atmo ! Determine the current and final year of the forcing cycle based on ! namelist input; initialize the atmospheric forcing data filenames. + use ice_calendar, only: use_leap_years + character(len=*), parameter :: subname = '(init_forcing_atmo)' ! Allocate forcing arrays @@ -235,6 +238,14 @@ subroutine init_forcing_atmo file=__FILE__, line=__LINE__) endif + if (use_leap_years .and. (trim(atm_data_type) /= 'JRA55' .and. & + trim(atm_data_type) /= 'default' .and. & + trim(atm_data_type) /= 'box2001')) then + write(nu_diag,*) 'use_leap_years option is currently only supported for' + write(nu_diag,*) 'JRA55, default , and box2001 atmospheric data' + call abort_ice(error_message=subname, file=__FILE__, line=__LINE__) + endif + !------------------------------------------------------------------- ! Get filenames for input forcing data !------------------------------------------------------------------- @@ -244,6 +255,8 @@ subroutine init_forcing_atmo call NCAR_files(fyear) elseif (trim(atm_data_type) == 'LYq') then call LY_files(fyear) + elseif (trim(atm_data_type) == 'JRA55') then + call JRA55_files(fyear) elseif (trim(atm_data_type) == 'hadgem') then call hadgem_files(fyear) elseif (trim(atm_data_type) == 'monthly') then @@ -539,6 +552,8 @@ subroutine get_forcing_atmo call ncar_data elseif (trim(atm_data_type) == 'LYq') then call LY_data + elseif (trim(atm_data_type) == 'JRA55') then + call JRA55_data(fyear) elseif (trim(atm_data_type) == 'hadgem') then call hadgem_data elseif (trim(atm_data_type) == 'monthly') then @@ -1376,6 +1391,10 @@ subroutine file_year (data_file, yr) i = index(data_file,'.nc') - 5 tmpname = data_file write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' + elseif (trim(atm_data_type) == 'JRA55') then ! netcdf + i = index(data_file,'.nc') - 5 + tmpname = data_file + write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc' else ! LANL/NCAR naming convention i = index(data_file,'.dat') - 5 tmpname = data_file @@ -2002,6 +2021,22 @@ subroutine LY_files (yr) endif ! master_task end subroutine LY_files + subroutine JRA55_files(yr) +! + integer (kind=int_kind), intent(in) :: & + yr ! current forcing year + + character(len=*), parameter :: subname = '(JRA55_files)' + + uwind_file = & + trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc' + call file_year(uwind_file,yr) + if (my_task == master_task) then + write (nu_diag,*) ' ' + write (nu_diag,*) 'Atmospheric data files:' + write (nu_diag,*) trim(uwind_file) + endif + end subroutine JRA55_files !======================================================================= ! @@ -2225,6 +2260,239 @@ subroutine LY_data end subroutine LY_data + subroutine JRA55_data (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) :: & + ncid , & ! netcdf file id + i, j , & + ixm,ixx,ixp , & ! record numbers for neighboring months + recnum , & ! record number + maxrec , & ! maximum record number + recslot , & ! spline slot for current record + midmonth , & ! middle day of month + dataloc , & ! = 1 for data located in middle of time interval + ! = 2 for date located at end of time interval + iblk , & ! block index + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + yr ! current forcing year + + real (kind=dbl_kind) :: & + sec3hr , & ! number of seconds in 3 hours + secday , & ! number of seconds in day + Tffresh , & + vmin, vmax + + logical (kind=log_kind) :: readm, read6,debug_n_d + + type (block) :: & + this_block ! block information for current block + + character(len=64) :: fieldname !netcdf field name + character(len=*), parameter :: subname = '(JRA55_data)' + + debug_n_d = .false. !usually false + + 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__) + + !------------------------------------------------------------------- + ! 3-hourly data + ! + ! Assume that the 3-hourly value is located at the end of the + ! 3-hour period. This is the convention for NCEP reanalysis data. + ! E.g. record 1 gives conditions at 3 am GMT on 1 January. + !------------------------------------------------------------------- + + dataloc = 2 ! data located at end of interval + sec3hr = secday/c8 ! seconds in 3 hours + !maxrec = 2920 ! 365*8; for leap years = 366*8 + + if(use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar + maxrec = days_per_year*8 + + if(days_per_year == 365 .and. (mod(yr, 4) == 0)) then + call abort_ice('days_per_year should be set to 366 for leap years') + end if + + ! current record number + recnum = 8*int(yday) - 7 + int(real(sec,kind=dbl_kind)/sec3hr) + + ! Compute record numbers for surrounding data (2 on each side) + + ixm = mod(recnum+maxrec-2,maxrec) + 1 + ixx = mod(recnum-1, maxrec) + 1 + + ! Compute interpolation coefficients + ! If data is located at the end of the time interval, then the + ! data value for the current record goes in slot 2 + + recslot = 2 + ixp = -99 + call interp_coeff (recnum, recslot, sec3hr, dataloc) + + ! Read + read6 = .false. + if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true. + !------------------------------------------------------------------- + ! 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) + !------------------------------------------------------------------- + call ice_open_nc(uwind_file,ncid) + + fieldname = 'airtmp' + call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'wndewd' + call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'wndnwd' + call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'spchmd' + call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'glbrad' + call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'dlwsfc' + call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname = 'ttlpcp' + call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + call ice_close_nc(ncid) + + + ! 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) + call interpolate_data (fsw_data, fsw) + call interpolate_data (flw_data, flw) + call interpolate_data (fsnow_data, fsnow) + + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + 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 + + ! Save record number + oldrecnum = recnum + + if (dbug) then + if (my_task == master_task) write (nu_diag,*) '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,*) '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,*) '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,*) '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,*) '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,*) '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,*) '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,*) 'Qa',vmin,vmax + if (my_task.eq.master_task) & + write (nu_diag,*) 'maxrec',maxrec + write (nu_diag,*) 'days_per_year', days_per_year + + endif ! dbug + + end subroutine JRA55_data + + !======================================================================= ! ! AOMIP shortwave forcing diff --git a/configuration/scripts/options/set_nml.boxrestore b/configuration/scripts/options/set_nml.boxrestore index bc913a3dc..6789e1ff8 100644 --- a/configuration/scripts/options/set_nml.boxrestore +++ b/configuration/scripts/options/set_nml.boxrestore @@ -2,7 +2,7 @@ nilyr = 1 ice_ic = 'default' restart = .false. restart_ext = .true. -use_leap_years = .true. +use_leap_years = .false. ndtd = 2 kcatbound = 1 ew_boundary_type = 'cyclic' diff --git a/configuration/scripts/options/set_nml.jra55 b/configuration/scripts/options/set_nml.jra55 new file mode 100755 index 000000000..d0112a857 --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55 @@ -0,0 +1,17 @@ +year_init = 2005 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2005 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.jra55_2008 b/configuration/scripts/options/set_nml.jra55_2008 new file mode 100755 index 000000000..042431fc0 --- /dev/null +++ b/configuration/scripts/options/set_nml.jra55_2008 @@ -0,0 +1,17 @@ +year_init = 2008 +ice_ic = 'ICE_MACHINE_INPUTDATA/CICE_data/ic/gx1/iced_gx1_v5.nc' +grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin' +kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin' +bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc' +use_leap_years = .true. +use_restart_time = .false. +maskhalo_dyn = .true. +maskhalo_remap = .true. +maskhalo_bound = .true. +fyear_init = 2008 +atm_data_dir = 'ICE_MACHINE_INPUTDATA/CICE_data/forcing/gx1/JRA55' +atm_data_format = 'nc' +atm_data_type = 'JRA55' +precip_units = 'mks' +ocn_data_dir = 'default' +bgc_data_dir = 'default' diff --git a/configuration/scripts/options/set_nml.run90day b/configuration/scripts/options/set_nml.run90day new file mode 100644 index 000000000..06db1a3d8 --- /dev/null +++ b/configuration/scripts/options/set_nml.run90day @@ -0,0 +1,5 @@ +npt = 2160 +dumpfreq = 'd' +dumpfreq_n = 30 +histfreq = 'd','x','x','x','x' +f_aice = 'd' diff --git a/configuration/scripts/options/test_nml.restart1 b/configuration/scripts/options/test_nml.restart1 index 4874c25a5..82f934720 100644 --- a/configuration/scripts/options/test_nml.restart1 +++ b/configuration/scripts/options/test_nml.restart1 @@ -2,3 +2,4 @@ npt = 240 dumpfreq = 'd' dumpfreq_n = 5 runtype = 'initial' +use_restart_time = .false. diff --git a/configuration/scripts/options/test_nml.restart2 b/configuration/scripts/options/test_nml.restart2 index 5981918a8..4ae10c5a6 100644 --- a/configuration/scripts/options/test_nml.restart2 +++ b/configuration/scripts/options/test_nml.restart2 @@ -3,4 +3,4 @@ dumpfreq = 'd' dumpfreq_n = 5 runtype = 'continue' restart = .true. - +use_restart_time = .true. diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts old mode 100644 new mode 100755 index 413047e7a..39f68a5f2 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -37,3 +37,5 @@ smoke gx3 8x1 bgcskl,debug #smoke gx3 4x1 bgcz,thread smoke_gx3_8x2_bgcz restart gx1 4x2 bgcsklclim,medium restart gx1 8x1 bgczclim,medium +smoke gx1 24x1 jra55_2008,medium,run90day +restart gx1 24x1 jra55,short diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 9c4a4d905..11bb3e887 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -928,6 +928,15 @@ @Article{Roberts18 volume = {376}, url = {http://dx.doi.org/10.1098/rsta.2017.0344} } +@Article{Tsujino18, + author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.", + title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}", + journal = OM, + year = {2018}, + volume = {130}, + pages = {79-139}, + url = {http://dx.doi.org/10.1016/j.ocemod.2018.07.002} +} % ********************************************** % For new entries, see example entry in BIB_TEMPLATE.txt % ********************************************** diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index beb7c7f1d..c0ee9331d 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -338,6 +338,7 @@ Table of namelist options "","``restore_ice``", "true/false", "restore ice state along lateral boundaries", "" "\*","``atm_data_type``", "``default``", "constant values defined in the code", "" "","", "``LYq``", "COREII Large-Yeager (AOMIP) forcing data", ":cite:`Large09`" + "","", "``JRA55``", "JRA55 forcing data :cite:`Tsujino18`", "" "","", "``monthly``", "monthly forcing data", "" "","", "``ncar``", "NCAR bulk forcing data", "" "","", "``box2001``", "forcing data for :cite:`Hunke01` box problem", ""