From d4888eee48ec117555bea439132ba80314148bb0 Mon Sep 17 00:00:00 2001 From: mhrib Date: Tue, 6 Nov 2018 10:49:28 +0000 Subject: [PATCH 1/5] Record number for popgrid_nc (NetCDF) is always one --- cicecore/cicedynB/infrastructure/ice_grid.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 6baddd63f..65a1ebccd 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -339,7 +339,7 @@ subroutine init_grid2 type (block) :: & this_block ! block information for current block - + character(len=*), parameter :: subname = '(init_grid2)' !----------------------------------------------------------------- @@ -719,13 +719,13 @@ subroutine popgrid_nc 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_kmt ! file id for netCDF kmt file + fid_grid, & ! file id for netCDF grid file + fid_kmt ! file id for netCDF kmt file logical (kind=log_kind) :: diag character (char_len) :: & - fieldname ! field name in netCDF file + fieldname ! field name in netCDF file real (kind=dbl_kind) :: & pi @@ -794,7 +794,7 @@ subroutine popgrid_nc ew_boundary_type, ns_boundary_type) fieldname='ulon' - call ice_read_global_nc(fid_grid,2,fieldname,work_g1,diag) ! ULON + call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ULON call gridbox_verts(work_g1,lont_bounds) call scatter_global(ULON, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) @@ -802,7 +802,7 @@ subroutine popgrid_nc ew_boundary_type, ns_boundary_type) fieldname='angle' - call ice_read_global_nc(fid_grid,7,fieldname,work_g1,diag) ! ANGLE + 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) @@ -817,11 +817,11 @@ subroutine popgrid_nc !----------------------------------------------------------------- fieldname='htn' - call ice_read_global_nc(fid_grid,3,fieldname,work_g1,diag) ! 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,4,fieldname,work_g1,diag) ! HTE + call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! HTE call primary_grid_lengths_HTE(work_g1) ! dyu, dyt deallocate(work_g1) From 36b5de516f2d65a83835a61aa17200a4b4464131 Mon Sep 17 00:00:00 2001 From: mhrib Date: Tue, 6 Nov 2018 10:51:57 +0000 Subject: [PATCH 2/5] HYCOM (NetCDF) I/O routines used by DMI --- cicecore/cicedynB/general/ice_forcing.F90 | 369 +++++++++++++++++- .../infrastructure/ice_read_write.F90 | 137 ++++++- cicecore/shared/ice_calendar.F90 | 48 ++- 3 files changed, 545 insertions(+), 9 deletions(-) diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index 6eaf4151b..20c367a12 100644 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -21,10 +21,11 @@ module ice_forcing use ice_communicate, only: my_task, master_task use ice_calendar, only: istep, istep1, time, time_forc, & sec, mday, month, nyr, yday, daycal, dayyr, & - daymo, days_per_year + daymo, days_per_year, hc_jday use ice_fileunits, only: nu_diag, nu_forcing use ice_exit, only: abort_ice use ice_read_write, only: ice_open, ice_read, & + ice_get_ncvarsize, ice_read_vec_nc, & ice_open_nc, ice_read_nc, ice_close_nc use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, & timer_bound @@ -118,7 +119,7 @@ module ice_forcing bgc_data_type, & ! 'default', 'clim', 'ncar', 'oned' ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned', ! 'hadgem_sst' or 'hadgem_sst_uvocn' - precip_units ! 'mm_per_month', 'mm_per_sec', 'mks' + precip_units ! 'mm_per_month', 'mm_per_sec', 'mks','m_per_sec' character(char_len_long), public :: & atm_data_dir , & ! top directory for atmospheric data @@ -151,6 +152,12 @@ module ice_forcing logical (kind=log_kind), public :: & dbug ! prints debugging output if true + real (dbl_kind), dimension(:), allocatable, public :: & + jday_atm ! jday time vector from atm forcing files + + integer (kind=int_kind), public :: & + Njday_atm ! Number of atm forcing timesteps + !======================================================================= contains @@ -242,6 +249,8 @@ subroutine init_forcing_atmo call ISPOL_files elseif (trim(atm_data_type) == 'box') then call box_data + elseif (trim(atm_data_type) == 'hycom_atm') then + call hycom_atm_files endif end subroutine init_forcing_atmo @@ -441,6 +450,11 @@ subroutine init_forcing_ocn(dt) call ocn_data_ncar_init_3D endif + if (trim(ocn_data_type) == 'hycom' .or. & + trim(bgc_data_type) == 'hycom') then + call ocn_data_hycom_init + endif + end subroutine init_forcing_ocn !======================================================================= @@ -534,6 +548,8 @@ subroutine get_forcing_atmo call oned_data elseif (trim(atm_data_type) == 'box') then call box_data + elseif (trim(atm_data_type) == 'hycom_atm') then + call hycom_atm_data else ! default values set in init_flux return endif @@ -621,6 +637,10 @@ subroutine get_forcing_ocn (dt) elseif (trim(ocn_data_type) == 'oned' .or. & trim(bgc_data_type) == 'oned') then call ocn_data_oned + elseif (trim(ocn_data_type) == 'hycom' .or. & + trim(bgc_data_type) == 'hycom') then +! call ocn_data_hycom(dt) +!MHRI: NOT IMPLEMENTED YET endif end subroutine get_forcing_ocn @@ -825,7 +845,7 @@ subroutine read_data_nc (flag, recd, yr, ixm, ixx, ixp, & character(len=*), parameter :: subname = '(read_data_nc)' -#ifdef ncdf +#ifdef ncdf integer (kind=int_kind) :: & nrec , & ! record number to read n2, n4 , & ! like ixm and ixp, but @@ -931,6 +951,81 @@ subroutine read_data_nc (flag, recd, yr, ixm, ixx, ixp, & #endif end subroutine read_data_nc +!======================================================================= + + subroutine read_data_nc_hycom (flag, recd, & + data_file, fieldname, field_data, & + field_loc, field_type) + +! Data is assumed to cover the entire time period of simulation. +! It is not bounded by start of year nor end of year +! Data must be accesible both before and after (or on) the point in time +! Assume increasing timeaxis within the forcing files, but they do not +! have to be equal spaced. Read time vector from "MT" in "init_hycom" +! +! Adapted by Mads Hvid Ribergaard, DMI from read_data_nc + + use ice_diagnostics, only: check_step + use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite + + logical (kind=log_kind), intent(in) :: flag + + integer (kind=int_kind), intent(in) :: & + recd ! baseline record number + + character (char_len_long) :: & + data_file ! data file to be read + character (char_len), intent(in) :: & + fieldname ! field name in netCDF file + + integer (kind=int_kind), intent(in) :: & + field_loc, & ! location of field on staggered grid + field_type ! type of field (scalar, vector, angle) + + real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), & + intent(out) :: & + field_data ! 2 values needed for interpolation + +#ifdef ncdf + ! local variables + integer (kind=int_kind) :: & + fid ! file id for netCDF routines + + call ice_timer_start(timer_readwrite) ! reading/writing + + if (istep1 > check_step) dbug = .true. !! debugging + + if (my_task==master_task .and. (dbug)) then + write(nu_diag,*) ' ', trim(data_file) + endif + + if (flag) then + + call ice_open_nc (data_file, fid) + !----------------------------------------------------------------- + ! read data + !----------------------------------------------------------------- + call ice_read_nc & + (fid, recd , fieldname, field_data(:,:,1,:), dbug, & + field_loc, field_type) + + call ice_read_nc & + (fid, recd+1, fieldname, field_data(:,:,2,:), dbug, & + field_loc, field_type) + + call ice_close_nc(fid) + + endif ! flag + + call ice_timer_stop(timer_readwrite) ! reading/writing + +#else + field_data = c0 ! to satisfy intent(out) attribute + write(*,*)'ERROR: CICE not compiled with NetCDF' + stop +#endif + end subroutine read_data_nc_hycom + !======================================================================= subroutine read_clim_data (readflag, recd, ixm, ixx, ixp, & @@ -1216,6 +1311,24 @@ subroutine interp_coeff (recnum, recslot, secint, dataloc) end subroutine interp_coeff +!======================================================================= + + subroutine interp_coeff2 (tt, t1, t2) + +! Compute coefficients for interpolating data to current time step. +! Works for any data interval using decimal daynumbers + + ! local variables + real (kind=dbl_kind), intent(in) :: & + tt , & ! current decimal daynumber + t1, t2 ! first+last decimal daynumber + + ! Compute coefficients + c1intp = abs((t2 - tt) / (t2 - t1)) + c2intp = c1 - c1intp + + end subroutine interp_coeff2 + !======================================================================= subroutine interpolate_data (field_data, field) @@ -1450,6 +1563,8 @@ subroutine prepare_forcing (nx_block, ny_block, & elseif (trim(precip_units) == 'mm_per_sec' .or. & trim(precip_units) == 'mks') then precip_factor = c1 ! mm/sec = kg/m^2 s + elseif (trim(precip_units) == 'm_per_sec') then + precip_factor = c1000 endif do j = jlo, jhi @@ -2913,7 +3028,7 @@ subroutine oned_data fieldname ! field name in netcdf file integer (kind=int_kind) :: & - fid ! file id for netCDF file + fid ! file id for netCDF file real (kind=dbl_kind):: & work ! temporary variable @@ -3881,6 +3996,252 @@ subroutine ocn_data_hadgem(dt) end subroutine ocn_data_hadgem +!======================================================================= + + subroutine ocn_data_hycom_init + ! Read SSS+SST from a HYCOM file converted to NetCDF format. + ! HYCOM binary2NetCDF: hcdata2ncdf2d (or hcdata2ncdf3z) + ! + rename/link file + use ice_blocks, only: nx_block, ny_block + use ice_domain, only: nblocks + use ice_domain_size, only: max_blocks + use ice_flux, only: sss, sst, Tf +#ifdef ncdf + use netcdf +#endif + + integer (kind=int_kind) :: & + i, j, iblk , & ! horizontal indices + fid ! file id for netCDF file + + character (char_len) :: & + fieldname ! field name in netcdf file + + if (trim(bgc_data_type) == 'hycom') then + sss_file = trim(ocn_data_dir)//'ice.restart.surf.nc' + + if (my_task == master_task) then + write (nu_diag,*)' ' + write (nu_diag,*)'Initial SSS file: ',trim(sss_file) + endif + + fieldname = 'sss' + call ice_open_nc (sss_file, fid) + call ice_read_nc (fid, 1 , fieldname, sss, dbug, & + field_loc_center, field_type_scalar) + call ice_close_nc(fid) + + call ocn_freezing_temperature + endif + + if (trim(ocn_data_type) == 'hycom') then + sst_file = trim(ocn_data_dir)//'ice.restart.surf.nc' + + if (my_task == master_task) then + write (nu_diag,*)' ' + write (nu_diag,*)'Initial SST file: ',trim(sst_file) + endif + + fieldname = 'sst' + call ice_open_nc (sst_file, fid) + call ice_read_nc (fid, 1 , fieldname, sst, dbug, & + field_loc_center, field_type_scalar) + call ice_close_nc(fid) + + ! Make sure sst is not less than freezing temperature Tf + !$OMP PARALLEL DO PRIVATE(iblk,i,j) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk)) + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif + + end subroutine ocn_data_hycom_init + +!======================================================================= + subroutine hycom_atm_files + + use ice_broadcast, only: broadcast_array, broadcast_scalar + + integer (kind = int_kind) :: & + fid ! File id + character (char_len) :: & + varname ! variable name in netcdf file + + fsw_file = trim(atm_data_dir)//'/forcing.shwflx.nc' + flw_file = trim(atm_data_dir)//'/forcing.radflx.nc' + rain_file = trim(atm_data_dir)//'/forcing.precip.nc' + uwind_file = trim(atm_data_dir)//'/forcing.ewndsp.nc' !actually Xward, not Eward + vwind_file = trim(atm_data_dir)//'/forcing.nwndsp.nc' !actually Yward, not Nward + tair_file = trim(atm_data_dir)//'/forcing.airtmp.nc' + humid_file = trim(atm_data_dir)//'/forcing.vapmix.nc' + + ! Read time vector from "tair_file" + call ice_open_nc(tair_file, fid) + varname='MT' + call ice_get_ncvarsize(fid,varname,Njday_atm) + call broadcast_scalar(Njday_atm,master_task) + allocate(jday_atm(Njday_atm)) + call ice_read_vec_nc(fid,Njday_atm, varname,jday_atm, .true.) + call ice_close_nc(fid) + call broadcast_array(jday_atm ,master_task) + + ! Write diag info + if (my_task == master_task) then + write (nu_diag,*) ' ' + write (nu_diag,*) 'CICE: Atm. (hycomdate) Start = ',jday_atm(1) + write (nu_diag,*) 'CICE: Atm. (hycomdate) End = ',jday_atm(Njday_atm) + write (nu_diag,*) 'CICE: Total Atm timesteps = ',Njday_atm + write (nu_diag,*) 'CICE: Atmospheric forcing files:' + write (nu_diag,*) trim(fsw_file) + write (nu_diag,*) trim(flw_file) + write (nu_diag,*) trim(rain_file) + write (nu_diag,*) trim(uwind_file) + write (nu_diag,*) trim(vwind_file) + write (nu_diag,*) trim(tair_file) + write (nu_diag,*) trim(humid_file) + endif ! master_task + + end subroutine hycom_atm_files + +!======================================================================= + + subroutine hycom_atm_data + + use ice_flux, only: fsw, fsnow, Tair, uatm, vatm, Qa, flw + use ice_domain, only: nblocks + use ice_calendar, only: year_init + + integer (kind=int_kind) :: & + recnum ! record number + + real (kind=dbl_kind) :: & + hcdate ! current time in HYCOM jday units + + logical (kind=log_kind) :: read6 + + character (char_len) :: & + fieldname ! field name in netcdf file + + integer (kind=int_kind) :: & + i, j, iblk ! horizontal indices + + real (kind=dbl_kind) :: Tffresh, secday + + character(len=*), parameter :: subname = '(hycom_atm_data)' + + call icepack_query_parameters(Tffresh_out=Tffresh) + call icepack_query_parameters(secday_out=secday) + + ! current time in HYCOM jday units + hcdate = hc_jday(nyr+year_init-1,0,0)+ yday+sec/secday + + ! Init recnum try + recnum=min(max(oldrecnum,1),Njday_atm-1) + + ! Find correct time in ATM data ... assume cont. incr. time-axis + do while ( recnum=jday_atm(recnum) .and. & + hcdate<=jday_atm(recnum+1) ) exit + if ( abs(hcdate-jday_atm(recnum))jday_atm(recnum+1)+p001 ) then + 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 + call abort_ice ('ERROR: CICE stopped') + endif + + ! Compute interpolation coefficients + call interp_coeff2 (hcdate, jday_atm(recnum), jday_atm(recnum+1) ) + + + ! Read + read6 = .false. + if (istep==1 .or. oldrecnum /= recnum) read6 = .true. + + if (trim(atm_data_format) == 'nc') then + + if (read6 .and. my_task == master_task) write(nu_diag,*) & + 'CICE: Atm. read: = ',jday_atm(recnum), jday_atm(recnum+1) + + fieldname = 'airtmp' + call read_data_nc_hycom (read6, recnum, & + tair_file, fieldname, Tair_data, & + field_loc_center, field_type_scalar) + fieldname = 'ewndsp' + call read_data_nc_hycom (read6, recnum, & + uwind_file, fieldname, uatm_data, & + field_loc_center, field_type_vector) + fieldname = 'nwndsp' + call read_data_nc_hycom (read6, recnum, & + vwind_file, fieldname, vatm_data, & + field_loc_center, field_type_vector) + fieldname = 'vapmix' + call read_data_nc_hycom (read6, recnum, & + humid_file, fieldname, Qa_data, & + field_loc_center, field_type_scalar) + fieldname = 'shwflx' + call read_data_nc_hycom (read6, recnum, & + fsw_file, fieldname, fsw_data, & + field_loc_center, field_type_scalar) + fieldname = 'radflx' + call read_data_nc_hycom (read6, recnum, & + flw_file, fieldname, flw_data, & + field_loc_center, field_type_scalar) + fieldname = 'precip' + call read_data_nc_hycom (read6, recnum, & + rain_file, fieldname, fsnow_data,& + field_loc_center, field_type_scalar) + + else + call abort_ice(subname//'ERROR: atm_data_format unavailable for hycom') + endif + + ! Interpolate + if (dbug) then + if (my_task == master_task) then + write(nu_diag,*)'CICE: Atm. interpolate: = ',& + hcdate,c1intp,c2intp + endif + endif + call interpolate_data (Tair_data, Tair) + call interpolate_data (uatm_data, uatm) + call interpolate_data (vatm_data, vatm) + call interpolate_data ( fsw_data, fsw) + call interpolate_data ( flw_data, flw) + call interpolate_data ( Qa_data, Qa) + call interpolate_data (fsnow_data, fsnow) + + ! Adjust data forcing to CICE units + !$OMP PARALLEL DO PRIVATE(iblk,i,j) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + ! Air temperature: Degrees --> Kelvin + Tair(i,j,iblk) = Tair(i,j,iblk) + Tffresh + enddo ! i + enddo ! j + enddo ! nblocks + !$OMP END PARALLEL DO + + ! Save record number for next time step + oldrecnum = recnum + + end subroutine hycom_atm_data + !======================================================================= ! subroutine read_data_nc_point (flag, recd, yr, ixm, ixx, ixp, & diff --git a/cicecore/cicedynB/infrastructure/ice_read_write.F90 b/cicecore/cicedynB/infrastructure/ice_read_write.F90 index 20eff9765..bf6e850fe 100644 --- a/cicecore/cicedynB/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedynB/infrastructure/ice_read_write.F90 @@ -41,6 +41,8 @@ module ice_read_write ice_write, & ice_write_nc, & ice_write_ext, & + ice_read_vec_nc, & + ice_get_ncvarsize, & ice_close_nc interface ice_write @@ -76,12 +78,15 @@ module ice_read_write ! ! author: Tony Craig, NCAR - subroutine ice_open(nu, filename, nbits) + subroutine ice_open(nu, filename, nbits, algn) integer (kind=int_kind), intent(in) :: & nu , & ! unit number nbits ! no. of bits per variable (0 for sequential access) + integer (kind=int_kind), intent(in), optional :: algn + integer (kind=int_kind) :: RecSize, Remnant + character (*) :: filename character(len=*), parameter :: subname = '(ice_open)' @@ -93,7 +98,18 @@ subroutine ice_open(nu, filename, nbits) open(nu,file=filename,form='unformatted') else ! direct access - open(nu,file=filename,recl=nx_global*ny_global*nbits/8, & + RecSize = nx_global*ny_global*nbits/8 + if (present(algn)) then + ! If data is keept in blocks using given sizes (=algn) + ! Used in eg. HYCOM binary files, which are stored as "blocks" dividable by 16384 bit (=algn) + if (algn /= 0) then + Remnant = modulo(RecSize,algn) + if (Remnant /= 0) then + RecSize = RecSize + (algn - Remnant) + endif + endif + endif + open(nu,file=filename,recl=RecSize, & form='unformatted',access='direct') endif ! nbits = 0 @@ -1842,7 +1858,7 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) nrec ! record number character (char_len), intent(in) :: & - varname ! field name in netcdf file + varname ! field name in netcdf file real (kind=dbl_kind), dimension(nx_global,ny_global), & intent(out) :: & @@ -2102,6 +2118,121 @@ subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, & #endif end subroutine ice_read_nc_uv +!======================================================================= +! Read a vector in a netcdf file. +! Just like ice_read_global_nc except that it returns a vector. +! work_g is a real vector +! +! Adapted by William Lipscomb, LANL, from ice_read +! Adapted by Ann Keen, Met Office, to read from a netcdf file + + subroutine ice_read_vec_nc (fid, nrec, varname, work_g, diag) + + integer (kind=int_kind), intent(in) :: & + fid , & ! file id + nrec ! record number + + character (char_len), intent(in) :: & + varname ! field name in netcdf file + + real (kind=dbl_kind), dimension(nrec), & + intent(out) :: & + work_g ! output array (real, 8-byte) + + logical (kind=log_kind) :: & + diag ! if true, write diagnostic output + + ! local variables + + character(len=*), parameter :: subname = '(ice_read_vec_nc)' + +#ifdef ncdf +! netCDF file diagnostics: + integer (kind=int_kind) :: & + varid, & ! netcdf id for field + status, & ! status output from netcdf routines + nvar ! sizes of netcdf vector + + real (kind=dbl_kind) :: & + amin, amax ! min, max values of input vector + + character (char_len) :: & + dimname ! dimension name +! + work_g(:) = c0 + + if (my_task == master_task) then + + !------------------------------------------------------------- + ! Find out ID of required variable + !------------------------------------------------------------- + + status = nf90_inq_varid(fid, trim(varname), varid) + + if (status /= nf90_noerr) then + call abort_ice (subname//'ERROR: Cannot find variable '//trim(varname) ) + endif + + !-------------------------------------------------------------- + ! Read global array + !-------------------------------------------------------------- + + status = nf90_get_var( fid, varid, work_g, & + start=(/1/), & + count=(/nrec/) ) + endif ! my_task = master_task + + !------------------------------------------------------------------- + ! optional diagnostics + !------------------------------------------------------------------- + + if (my_task == master_task .and. diag) then + amin = minval(work_g) + amax = maxval(work_g) + write(nu_diag,*) 'min, max, nrec = ', amin, amax, nrec + endif + +#else + write(*,*) 'ERROR: ncdf not defined during compilation' + work_g = c0 ! to satisfy intent(out) attribute +#endif + end subroutine ice_read_vec_nc + +!======================================================================= +! Get number of variables of a given variable + subroutine ice_get_ncvarsize(fid,varname,recsize) + + integer (kind=int_kind), intent(in) :: & + fid ! file id + character (char_len), intent(in) :: & + varname ! field name in netcdf file + integer (kind=int_kind), intent(out) :: & + recsize ! Number of records in file + integer (kind=int_kind) :: & + ndims, i, status + character (char_len) :: & + cvar + character(len=*), parameter :: subname = '(ice_get_ncvarsize)' + + if (my_task == master_task) then + status=nf90_inquire(fid, nDimensions = nDims) + if (status /= nf90_noerr) then + call abort_ice (subname//'ERROR: inquire nDimensions' ) + endif + do i=1,nDims + status = nf90_inquire_dimension(fid,i,name=cvar,len=recsize) + if (status /= nf90_noerr) then + call abort_ice (subname//'ERROR: inquire len for variable '//trim(cvar) ) + endif + call flush(nu_diag) + if (trim(cvar) == trim(varname)) exit + enddo + if (trim(cvar) .ne. trim(varname)) then + call abort_ice (subname//'ERROR: Did not find variable '//trim(varname) ) + endif + endif + end subroutine ice_get_ncvarsize + !======================================================================= end module ice_read_write diff --git a/cicecore/shared/ice_calendar.F90 b/cicecore/shared/ice_calendar.F90 index c721281ba..85b827038 100644 --- a/cicecore/shared/ice_calendar.F90 +++ b/cicecore/shared/ice_calendar.F90 @@ -25,7 +25,7 @@ module ice_calendar implicit none private - public :: init_calendar, calendar, time2sec, sec2time + public :: init_calendar, calendar, time2sec, sec2time, hc_jday integer (kind=int_kind), public :: & days_per_year , & ! number of days in one year @@ -64,7 +64,7 @@ module ice_calendar istep0 , & ! counter, number of steps taken in previous run istep1 , & ! counter, number of steps at current timestep mday , & ! day of the month - hour , & ! hour of the year + hour , & ! hour of the day month , & ! month number, 1 to 12 monthp , & ! last month year_init, & ! initial year @@ -311,6 +311,9 @@ subroutine calendar(ttime) case ("d", "D") if (new_day .and. mod(elapsed_days, dumpfreq_n)==0) & write_restart = 1 + case ("h", "H") + if (new_hour .and. mod(elapsed_hours, dumpfreq_n)==0) & + write_restart = 1 end select if (force_restart_now) write_restart = 1 @@ -530,6 +533,47 @@ subroutine set_calendar(year) end subroutine set_calendar +!======================================================================= + + real(kind=dbl_kind) function hc_jday(iyear,imm,idd,ihour) +!-------------------------------------------------------------------- +! converts "calendar" date to HYCOM julian day: +! 1) year,month,day,hour (4 arguments) +! 2) year,doy,hour (3 arguments) +! +! HYCOM model day is calendar days since 31/12/1900 +!-------------------------------------------------------------------- + real(kind=dbl_kind) :: dtime + integer(kind=int_kind) :: iyear,iyr,imm,idd,idoy,ihr + integer(kind=int_kind), optional :: ihour + + if (present(ihour)) then + !----------------- + ! yyyy mm dd HH + !----------------- + iyr=iyear-1901 + if (mod(iyr,4)==3) then + dtime = floor(365.25_dbl_kind*iyr)*c1 + daycal366(imm)*c1 + idd*c1 + ihour/24._dbl_kind + else + dtime = floor(365.25_dbl_kind*iyr)*c1 + daycal365(imm)*c1 + idd*c1 + ihour/24._dbl_kind + endif + + else + !----------------- + ! yyyy DOY HH + !----------------- + ihr = idd ! redefine input + idoy = imm ! redefine input + iyr = iyear - 1901 + dtime = floor(365.25_dbl_kind*iyr)*c1 + idoy*c1 + ihr/24._dbl_kind + + endif + + hc_jday=dtime + + return + end function hc_jday + !======================================================================= end module ice_calendar From f25632015ac7da85e1360493e6b41b4265965c43 Mon Sep 17 00:00:00 2001 From: mhrib Date: Tue, 6 Nov 2018 11:15:49 +0000 Subject: [PATCH 3/5] Updated documentation for HYCOM forcing, precip units and restart dump in hours interval --- doc/source/user_guide/ug_case_settings.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 42f73f9e9..669dbb450 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -110,6 +110,7 @@ Table of namelist options "``dumpfreq``", "``y``", "write restart every ``dumpfreq_n`` years", "y" "", "``m``", "write restart every ``dumpfreq_n`` months", "" "", "``d``", "write restart every ``dumpfreq_n`` days", "" + "", "``h``", "write restart every ``dumpfreq_n`` hours", "" "``dumpfreq_n``", "integer", "frequency restart data is written", "1" "``dump_last``", "true/false", "if true, write restart on last time step of simulation", "" "", "", "*Model Output*", "" @@ -342,6 +343,7 @@ Table of namelist options "``precip_units``", "``mks``", "liquid precipitation data units", "" "", "``mm_per_month``", "", "" "", "``mm_per_sec``", "(same as MKS units)", "" + "", "``m_per_sec``", "", "" "``tfrz_option``", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8^{\circ} C`)", "" "", "``linear_salt``", "linear function of salinity (ktherm=1)", "" "", "``mushy_layer``", "matches mushy-layer thermo (ktherm=2)", "" @@ -362,12 +364,15 @@ Table of namelist options "", "``monthly``", "monthly forcing data", "" "", "``ncar``", "NCAR bulk forcing data", "" "", "``oned``", "column forcing data", "" + "", "``hycom_atm``", "HYCOM atm forcing in NetCDF format", "" "``ocn_data_type``", "``default``", "constant values defined in the code", "" "", "``clim``", "climatological data", "" "", "``ncar``", "POP ocean forcing data", "" + "", "``hycom``", "HYCOM ocean forcing data.", "Constant initial forcing" "``bgc_data_type``", "``default``", "constant values defined in the code", "" "", "``clim``", "climatological data", "" "", "``near``", "POP ocean forcing data", "" + "", "``hycom``", "HYCOM ocean forcing data.", "Constant initial forcing" "``sil_data_type``", "``default``", "default forcing value for silicate", "" "", "``clim``", "silicate forcing from ocean climatology :cite:`Garcia06`", "" "``nit_data_type``", "``default``", "default forcing value for nitrate", "" From 788c6a84d8d36aabf18aff695818905d16c7f365 Mon Sep 17 00:00:00 2001 From: mhrib Date: Tue, 6 Nov 2018 19:44:32 +0000 Subject: [PATCH 4/5] "hycom_atm" option changed to "hycom" --- cicecore/cicedynB/general/ice_forcing.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index 20c367a12..d4ce5a8db 100644 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -249,7 +249,7 @@ subroutine init_forcing_atmo call ISPOL_files elseif (trim(atm_data_type) == 'box') then call box_data - elseif (trim(atm_data_type) == 'hycom_atm') then + elseif (trim(atm_data_type) == 'hycom') then call hycom_atm_files endif @@ -548,7 +548,7 @@ subroutine get_forcing_atmo call oned_data elseif (trim(atm_data_type) == 'box') then call box_data - elseif (trim(atm_data_type) == 'hycom_atm') then + elseif (trim(atm_data_type) == 'hycom') then call hycom_atm_data else ! default values set in init_flux return From 7b4dbbad0527c2946c7b281341b173f6f22041bc Mon Sep 17 00:00:00 2001 From: mhrib Date: Tue, 6 Nov 2018 19:45:05 +0000 Subject: [PATCH 5/5] option "HYCOM_atm" -> "hycom". Correct "near" -> "ncar" --- doc/source/user_guide/ug_case_settings.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 669dbb450..f39301866 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -364,15 +364,15 @@ Table of namelist options "", "``monthly``", "monthly forcing data", "" "", "``ncar``", "NCAR bulk forcing data", "" "", "``oned``", "column forcing data", "" - "", "``hycom_atm``", "HYCOM atm forcing in NetCDF format", "" + "", "``hycom``", "HYCOM atm forcing data in NetCDF format", "" "``ocn_data_type``", "``default``", "constant values defined in the code", "" "", "``clim``", "climatological data", "" "", "``ncar``", "POP ocean forcing data", "" - "", "``hycom``", "HYCOM ocean forcing data.", "Constant initial forcing" + "", "``hycom``", "HYCOM ocean forcing data in NetCDF format", "Constant initial forcing" "``bgc_data_type``", "``default``", "constant values defined in the code", "" "", "``clim``", "climatological data", "" - "", "``near``", "POP ocean forcing data", "" - "", "``hycom``", "HYCOM ocean forcing data.", "Constant initial forcing" + "", "``ncar``", "POP ocean forcing data", "" + "", "``hycom``", "HYCOM ocean forcing data in NetCDF format", "Constant initial forcing" "``sil_data_type``", "``default``", "default forcing value for silicate", "" "", "``clim``", "silicate forcing from ocean climatology :cite:`Garcia06`", "" "``nit_data_type``", "``default``", "default forcing value for nitrate", ""