!=======================================================================
!BOP
!
! !MODULE: ice_init - parameter and variable initializations
!
! !DESCRIPTION:
!
! parameter and variable initializations
!
! !REVISION HISTORY:
! SVN:$Id: ice_init.F90 143 2008-08-08 23:08:22Z eclare $
!
! authors Elizabeth C. Hunke and William H. Lipscomb, LANL
! C. M. Bitz, UW
!
! 2004 WHL: Block structure added
! 2006 ECH: Added namelist variables, warnings.
! Replaced old default initial ice conditions with 3.14 version.
! Converted to free source form (F90).
!
! !INTERFACE:
!
module ice_init 1,4
!
! !USES:
!
use ice_kinds_mod
use ice_communicate, only: my_task, master_task, lprint_stats
use ice_domain_size
use ice_constants
!
!EOP
!
implicit none
save
!=======================================================================
contains
!=======================================================================
!BOP
!
! !IROUTINE: input_data - namelist variables
!
! !INTERFACE:
!
subroutine input_data 1,133
!
! !DESCRIPTION:
!
! Namelist variables, set to default values; may be altered
! at run time
!
! !REVISION HISTORY:
!
! author Elizabeth C. Hunke, LANL
!
! !USES:
!
use ice_broadcast
use ice_diagnostics
use ice_fileunits
use ice_calendar, only: year_init, istep0, histfreq, histfreq_n, &
dumpfreq, dumpfreq_n, diagfreq, &
npt, dt, xndt_dyn, days_per_year, write_ic
use ice_restart, only: &
restart, restart_dir, restart_file, pointer_file, &
runid, runtype, ice_ic, resttype, restart_format, lcdf64
use ice_history_fields, only: hist_avg, &
history_format, history_dir, history_file, &
incond_dir, incond_file
use ice_exit
use ice_itd, only: kitd, kcatbound
use ice_ocean, only: oceanmixed_ice
use ice_flux, only: Tfrzpt, update_ocn_f
use ice_forcing, only: &
ycycle, fyear_init, dbug, &
atm_data_type, atm_data_dir, precip_units, &
atm_data_format, ocn_data_format, &
sss_data_type, sst_data_type, ocn_data_dir, &
oceanmixed_file, restore_sst, trestore
use ice_grid, only: grid_file, kmt_file, grid_type, grid_format
use ice_mechred, only: kstrength, krdg_partic, krdg_redist
use ice_dyn_evp, only: ndte, kdyn, evp_damping, yield_curve
use ice_shortwave, only: albicev, albicei, albsnowv, albsnowi, &
shortwave, albedo_type, R_ice, R_pnd, &
R_snw, dT_mlt_in, rsnw_melt_in
use ice_atmo, only: atmbndy, calc_strair
use ice_transport_driver, only: advection
use ice_state, only: nt_Tsfc, nt_iage, nt_FY, nt_volpn, nt_aero, &
tr_aero, tr_iage, tr_FY, tr_pond, &
nt_alvl, nt_vlvl, tr_lvl
use ice_aerosol, only: restart_aero
use ice_age, only: restart_age
use ice_FY, only: restart_FY
use ice_lvl, only: restart_lvl
use ice_meltpond, only: restart_pond
use ice_therm_vertical, only: calc_Tsfc, heat_capacity
use ice_restoring
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
integer (kind=int_kind) :: &
nml_error, & ! namelist i/o error flag
ntr , & ! counter for number of tracers turned on
ns
character (len=6) :: chartmp
#ifdef CCSMCOUPLED
logical :: exists
#endif
!-----------------------------------------------------------------
! Namelist variables.
!-----------------------------------------------------------------
namelist /setup_nml/ &
days_per_year, year_init, istep0, dt, &
npt, xndt_dyn, &
runtype, runid, &
ice_ic, restart, &
restart_dir, restart_file, restart_format, lcdf64, &
pointer_file, dumpfreq, dumpfreq_n, &
diagfreq, diag_type, diag_file, &
print_global, print_points, latpnt, lonpnt, &
dbug, histfreq, histfreq_n, hist_avg, &
history_dir, history_file, history_format, &
write_ic, incond_dir, incond_file, lprint_stats
namelist /grid_nml/ &
grid_format, grid_type, grid_file, kmt_file, &
kcatbound
namelist /ice_nml/ &
kitd, kdyn, ndte, &
evp_damping, yield_curve, &
kstrength, krdg_partic, krdg_redist, advection, &
heat_capacity, shortwave, albedo_type, &
albicev, albicei, albsnowv, albsnowi, &
R_ice, R_pnd, R_snw, &
dT_mlt_in, rsnw_melt_in, &
atmbndy, fyear_init, ycycle, atm_data_format,&
atm_data_type, atm_data_dir, calc_strair, calc_Tsfc, &
precip_units, Tfrzpt, update_ocn_f, &
oceanmixed_ice, ocn_data_format, sss_data_type, sst_data_type, &
ocn_data_dir, oceanmixed_file, restore_sst, trestore, &
restore_ice
namelist /tracer_nml/ &
tr_iage, restart_age, &
tr_FY, restart_FY, &
tr_lvl, restart_lvl, &
tr_pond, restart_pond, &
tr_aero, restart_aero !MH for soot
!-----------------------------------------------------------------
! default values
!-----------------------------------------------------------------
days_per_year = 365 ! number of days in a year
year_init = 0 ! initial year
istep0 = 0 ! no. of steps taken in previous integrations,
! real (dumped) or imagined (to set calendar)
#ifndef CCSMCOUPLED
dt = 3600.0_dbl_kind ! time step, s
#endif
npt = 99999 ! total number of time steps (dt)
diagfreq = 24 ! how often diag output is written
print_points = .false. ! if true, print point data
print_global = .true. ! if true, print global diagnostic data
lprint_stats = .false. ! if true, prints decomposition statistics
diag_type = 'stdout'
diag_file = 'ice_diag.d'
histfreq(:) = 'd' ! output frequency option for different streams **** note this is only for testing - needs to get changed back'
histfreq_n(:) = 1 ! output frequency of histfreq
hist_avg = .true. ! if true, write time-averages (not snapshots)
history_dir = ' ' ! write to executable dir for default
history_file = 'iceh' ! history file name prefix
history_format = 'nc' ! file format ('bin'=binary or 'nc'=netcdf)
write_ic = .false. ! write out initial condition
incond_dir = history_dir ! write to history dir for default
incond_file = 'iceh_ic'! file prefix
dumpfreq='y' ! restart frequency option
dumpfreq_n = 1 ! restart frequency
restart = .false. ! if true, read restart files for initialization
restart_dir = ' ' ! write to executable dir for default
restart_file = 'iced' ! restart file name prefix
restart_format = 'nc' ! file format ('bin'=binary or 'nc'=netcdf)
lcdf64 = .false. ! 64 bit offset for netCDF
pointer_file = 'ice.restart_file'
ice_ic = 'default' ! latitude and sst-dependent
grid_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
grid_type = 'rectangular' ! define rectangular grid internally
grid_file = 'unknown_grid_file'
kmt_file = 'unknown_kmt_file'
kitd = 1 ! type of itd conversions (0 = delta, 1 = linear)
kcatbound = 1 ! category boundary formula (0 = old, 1 = new)
kdyn = 1 ! type of dynamics (1 = evp)
xndt_dyn = c1 ! dynamic time steps per thermodynamic time step
ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte
evp_damping = .false. ! if true, use damping procedure in evp dynamics
yield_curve = 'ellipse'
kstrength = 1 ! 1 = Rothrock 75 strength, 0 = Hibler 79
krdg_partic = 1 ! 1 = new participation, 0 = Thorndike et al 75
krdg_redist = 1 ! 1 = new redistribution, 0 = Hibler 80
advection = 'remap' ! incremental remapping transport scheme
shortwave = 'default' ! or 'dEdd' (delta-Eddington)
albedo_type = 'default'! or 'constant'
heat_capacity = .true. ! nonzero heat capacity (F => 0-layer thermo)
calc_Tsfc = .true. ! calculate surface temperature
Tfrzpt = 'linear_S' ! ocean freezing temperature, 'constant'=-1.8C
update_ocn_f = .false. ! include fresh water and salt fluxes for frazil
R_ice = 0.00_dbl_kind ! tuning parameter for sea ice
R_pnd = 0.00_dbl_kind ! tuning parameter for ponded sea ice
R_snw = 0.00_dbl_kind ! tuning parameter for snow over sea ice
dT_mlt_in = 1.50_dbl_kind ! melt transition (tuning)
rsnw_melt_in = 1500._dbl_kind ! max snow grain radius (tuning)
albicev = 0.78_dbl_kind ! visible ice albedo for h > ahmax
albicei = 0.36_dbl_kind ! near-ir ice albedo for h > ahmax
albsnowv = 0.98_dbl_kind ! cold snow albedo, visible
albsnowi = 0.70_dbl_kind ! cold snow albedo, near IR
atmbndy = 'default' ! or 'constant'
fyear_init = 1900 ! first year of forcing cycle
ycycle = 1 ! number of years in forcing cycle
atm_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
atm_data_type = 'default'
atm_data_dir = ' '
calc_strair = .true. ! calculate wind stress
precip_units = 'mks' ! 'mm_per_month' or
! 'mm_per_sec' = 'mks' = kg/m^2 s
oceanmixed_ice = .false. ! if true, use internal ocean mixed layer
ocn_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
sss_data_type = 'default'
sst_data_type = 'default'
ocn_data_dir = ' '
oceanmixed_file = 'unknown_oceanmixed_file' ! ocean forcing data
restore_sst = .false. ! restore sst if true
trestore = 90 ! restoring timescale, days (0 instantaneous)
restore_ice = .false. ! restore ice state on grid edges if true
dbug = .false. ! true writes diagnostics for input forcing
latpnt(1) = 90._dbl_kind ! latitude of diagnostic point 1 (deg)
lonpnt(1) = 0._dbl_kind ! longitude of point 1 (deg)
latpnt(2) = -65._dbl_kind ! latitude of diagnostic point 2 (deg)
lonpnt(2) = -45._dbl_kind ! longitude of point 2 (deg)
#ifndef CCSMCOUPLED
runid = 'unknown' ! run ID, only used in CCSM
runtype = 'initial' ! run type: 'initial', 'continue'
#endif
!-----------------------------------------------------------------
! extra tracers (no longer namelist variables set in ice_domain_size)
!-----------------------------------------------------------------
tr_iage = .false. ! ice age
restart_age = .false. ! ice age restart
filename_iage = 'none'
tr_FY = .false. ! FY ice area
restart_FY = .false. ! FY ice area restart
filename_FY = 'none'
tr_lvl = .false. ! FY ice area
restart_lvl = .false. ! FY ice area restart
filename_lvl = 'none'
tr_pond = .false. ! explicit melt ponds
restart_pond = .false. ! melt ponds restart
filename_volpn = 'none'
tr_aero = .false. ! aerosols MH
restart_aero = .false. ! aerosol restart MH
filename_aero = 'none'
resttype = 'old'
!-----------------------------------------------------------------
! read from input file
!-----------------------------------------------------------------
call get_fileunit(nu_nml)
if (my_task == master_task) then
open (nu_nml, file=nml_filename, status='old',iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
do while (nml_error > 0)
print*,'Reading setup_nml'
read(nu_nml, nml=setup_nml,iostat=nml_error)
print*,'Reading grid_nml'
read(nu_nml, nml=grid_nml,iostat=nml_error)
print*,'Reading ice_nml'
read(nu_nml, nml=ice_nml,iostat=nml_error)
print*,'Reading tracer_nml'
read(nu_nml, nml=tracer_nml,iostat=nml_error)
if (nml_error > 0) read(nu_nml,*) ! for Nagware compiler
end do
if (nml_error == 0) close(nu_nml)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call abort_ice('ice: error reading namelist')
endif
call release_fileunit(nu_nml)
! Need to broadcast ice_ic at this point.
call broadcast_scalar(ice_ic, master_task)
!-----------------------------------------------------------------
! set up diagnostics output and resolve conflicts
!-----------------------------------------------------------------
#ifdef CCSMCOUPLED
! Note that diag_file is not utilized in CCSMCOUPLED mode
if (my_task == master_task) then
inquire(file='ice_modelio.nml',exist=exists)
if (exists) then
nu_diag = shr_file_getUnit()
call shr_file_setIO('ice_modelio.nml',nu_diag)
end if
! Note in CCSMCOUPLED mode the runid and runtype flag are obtained from
! the sequential driver - not from the cice namelist
history_file = trim(runid) // ".cice.h"
restart_file = trim(runid) // ".cice.r"
incond_file = trim(runid) // ".cice.i."
end if
if (trim(ice_ic) /= 'default' .and. trim(ice_ic) /= 'none') &
restart = .true.
#else
if (trim(diag_type) == 'file') call get_fileunit(nu_diag)
#endif
if (my_task == master_task) then
if (trim(diag_type) == 'file') then
write(ice_stdout,*) 'Diagnostic output will be in file ',diag_file
open (nu_diag, file=diag_file, status='unknown')
endif
write(nu_diag,*) '--------------------------------'
write(nu_diag,*) ' CICE model diagnostic output '
write(nu_diag,*) '--------------------------------'
write(nu_diag,*) ' '
endif
if (runtype == 'continue') restart = .true.
if (runtype /= 'continue' .and. (restart)) then
if (ice_ic == 'none' .or. ice_ic == 'default') then
if (my_task == master_task) then
write(nu_diag,*) &
'WARNING: runtype, restart, ice_ic are inconsistent:'
write(nu_diag,*) runtype, restart, ice_ic
write(nu_diag,*) &
'WARNING: Need ice_ic = <filename>.'
write(nu_diag,*) &
'WARNING: Initializing using ice_ic conditions'
endif
restart = .false.
endif
endif
if (runtype == 'initial' .and. .not.(restart)) then
if (ice_ic /= 'none' .and. ice_ic /= 'default') then
if (my_task == master_task) then
write(nu_diag,*) &
'WARNING: runtype, restart, ice_ic are inconsistent:'
write(nu_diag,*) runtype, restart, ice_ic
write(nu_diag,*) &
'WARNING: Initializing with NO ICE: '
endif
ice_ic = 'none'
endif
endif
#ifndef ncdf
! netcdf is unavailable
history_format = 'bin'
restart_format = 'bin'
grid_format = 'bin'
atm_data_format = 'bin'
ocn_data_format = 'bin'
#endif
if (days_per_year /= 365) shortwave = 'default' ! definite conflict
chartmp = advection(1:6)
if (chartmp /= 'upwind' .and. chartmp /= 'remap ') advection = 'remap'
if (ncat == 1 .and. kitd == 1) then
write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1'
write (nu_diag,*) 'Using the delta function ITD option instead'
kitd = 0
endif
if (trim(atm_data_type) == 'monthly' .and. calc_strair) &
calc_strair = .false.
if (trim(atm_data_type) == 'hadgem' .and. &
trim(precip_units) /= 'mks') then
if (my_task == master_task) &
write (nu_diag,*) &
'WARNING: HadGEM atmospheric data chosen with wrong precip_units'
write (nu_diag,*) &
'WARNING: Changing precip_units to mks (i.e. kg/m2 s).'
precip_units='mks'
endif
call broadcast_scalar(days_per_year, master_task)
call broadcast_scalar(year_init, master_task)
call broadcast_scalar(istep0, master_task)
call broadcast_scalar(dt, master_task)
call broadcast_scalar(npt, master_task)
call broadcast_scalar(diagfreq, master_task)
call broadcast_scalar(print_points, master_task)
call broadcast_scalar(print_global, master_task)
call broadcast_scalar(diag_type, master_task)
call broadcast_scalar(diag_file, master_task)
call broadcast_scalar(history_format, master_task)
! call broadcast_scalar(histfreq, master_task)
! call broadcast_scalar(histfreq_n, master_task)
do ns=1,max_nstrm
call broadcast_scalar(histfreq(ns), master_task)
call broadcast_scalar(histfreq_n(ns), master_task)
enddo
! call broadcast_array(histfreq_n(:), master_task)
call broadcast_scalar(hist_avg, master_task)
call broadcast_scalar(history_dir, master_task)
call broadcast_scalar(history_file, master_task)
call broadcast_scalar(write_ic, master_task)
call broadcast_scalar(incond_dir, master_task)
call broadcast_scalar(incond_file, master_task)
call broadcast_scalar(dumpfreq, master_task)
call broadcast_scalar(dumpfreq_n, master_task)
call broadcast_scalar(restart_file, master_task)
call broadcast_scalar(restart, master_task)
call broadcast_scalar(restart_dir, master_task)
call broadcast_scalar(restart_format, master_task)
call broadcast_scalar(lcdf64, master_task)
call broadcast_scalar(pointer_file, master_task)
call broadcast_scalar(grid_format, master_task)
call broadcast_scalar(grid_type, master_task)
call broadcast_scalar(grid_file, master_task)
call broadcast_scalar(kmt_file, master_task)
call broadcast_scalar(kitd, master_task)
call broadcast_scalar(kcatbound, master_task)
call broadcast_scalar(kdyn, master_task)
call broadcast_scalar(xndt_dyn, master_task)
call broadcast_scalar(ndte, master_task)
call broadcast_scalar(evp_damping, master_task)
call broadcast_scalar(yield_curve, master_task)
call broadcast_scalar(kstrength, master_task)
call broadcast_scalar(krdg_partic, master_task)
call broadcast_scalar(krdg_redist, master_task)
call broadcast_scalar(advection, master_task)
call broadcast_scalar(shortwave, master_task)
call broadcast_scalar(albedo_type, master_task)
call broadcast_scalar(heat_capacity, master_task)
call broadcast_scalar(R_ice, master_task)
call broadcast_scalar(R_pnd, master_task)
call broadcast_scalar(R_snw, master_task)
call broadcast_scalar(dT_mlt_in, master_task)
call broadcast_scalar(rsnw_melt_in, master_task)
call broadcast_scalar(albicev, master_task)
call broadcast_scalar(albicei, master_task)
call broadcast_scalar(albsnowv, master_task)
call broadcast_scalar(albsnowi, master_task)
call broadcast_scalar(atmbndy, master_task)
call broadcast_scalar(fyear_init, master_task)
call broadcast_scalar(ycycle, master_task)
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(calc_strair, master_task)
call broadcast_scalar(calc_Tsfc, master_task)
call broadcast_scalar(Tfrzpt, master_task)
call broadcast_scalar(update_ocn_f, master_task)
call broadcast_scalar(precip_units, master_task)
call broadcast_scalar(oceanmixed_ice, master_task)
call broadcast_scalar(ocn_data_format, master_task)
call broadcast_scalar(sss_data_type, master_task)
call broadcast_scalar(sst_data_type, master_task)
call broadcast_scalar(ocn_data_dir, master_task)
call broadcast_scalar(oceanmixed_file, master_task)
call broadcast_scalar(restore_sst, master_task)
call broadcast_scalar(trestore, master_task)
call broadcast_scalar(restore_ice, master_task)
call broadcast_scalar(dbug, master_task)
call broadcast_array (latpnt(1:2), master_task)
call broadcast_array (lonpnt(1:2), master_task)
call broadcast_scalar(runid, master_task)
call broadcast_scalar(runtype, master_task)
if (dbug) & ! else only master_task writes to file
call broadcast_scalar(nu_diag, master_task)
! tracers
call broadcast_scalar(tr_iage, master_task)
call broadcast_scalar(restart_age, master_task)
call broadcast_scalar(tr_FY, master_task)
call broadcast_scalar(restart_FY, master_task)
call broadcast_scalar(tr_lvl, master_task)
call broadcast_scalar(restart_lvl, master_task)
call broadcast_scalar(tr_pond, master_task)
call broadcast_scalar(restart_pond, master_task)
call broadcast_scalar(tr_aero, master_task)
call broadcast_scalar(restart_aero, master_task)
!-----------------------------------------------------------------
! spew
!-----------------------------------------------------------------
if (my_task == master_task) then
write(nu_diag,*) ' Document ice_in namelist parameters:'
write(nu_diag,*) ' ==================================== '
write(nu_diag,*) ' '
if (trim(runid) /= 'unknown') &
write(nu_diag,*) ' runid = ', &
trim(runid)
write(nu_diag,1030) ' runtype = ', &
trim(runtype)
write(nu_diag,1020) ' days_per_year = ', days_per_year
write(nu_diag,1020) ' year_init = ', year_init
write(nu_diag,1020) ' istep0 = ', istep0
write(nu_diag,1000) ' dt = ', dt
write(nu_diag,1020) ' npt = ', npt
write(nu_diag,1020) ' diagfreq = ', diagfreq
write(nu_diag,1010) ' print_global = ', &
print_global
write(nu_diag,1010) ' print_points = ', &
print_points
write(nu_diag,1050) ' histfreq = ', histfreq(:)
write(nu_diag,1040) ' histfreq_n = ', histfreq_n(:)
write(nu_diag,1010) ' hist_avg = ', hist_avg
if (hist_avg) then
write (nu_diag,*) 'History data will be averaged over ', &
histfreq_n,' ',histfreq
else
write (nu_diag,*) 'History data will be snapshots'
endif
write(nu_diag,*) ' history_dir = ', &
trim(history_dir)
write(nu_diag,*) ' history_file = ', &
trim(history_file)
write(nu_diag,*) ' history_format = ', &
trim(history_format)
if (write_ic) then
write (nu_diag,*) 'Initial condition will be written in ', &
trim(incond_dir)
endif
write(nu_diag,1030) ' dumpfreq = ', &
trim(dumpfreq)
write(nu_diag,1020) ' dumpfreq_n = ', dumpfreq_n
write(nu_diag,1010) ' restart = ', restart
write(nu_diag,*) ' restart_dir = ', &
trim(restart_dir)
write(nu_diag,*) ' restart_file = ', &
trim(restart_file)
write(nu_diag,*) ' restart_format = ', &
trim(restart_format)
write(nu_diag,*) ' lcdf64 = ', &
lcdf64
write(nu_diag,*) ' pointer_file = ', &
trim(pointer_file)
write(nu_diag,* ) ' ice_ic = ', &
trim(ice_ic)
write(nu_diag,*) ' grid_type = ', &
trim(grid_type)
if (trim(grid_type) /= 'rectangular' .or. &
trim(grid_type) /= 'column') then
write(nu_diag,*) ' grid_file = ', &
trim(grid_file)
write(nu_diag,*) ' kmt_file = ', &
trim(kmt_file)
endif
write(nu_diag,1020) ' kitd = ', kitd
write(nu_diag,1020) ' kcatbound = ', &
kcatbound
write(nu_diag,1020) ' kdyn = ', kdyn
write(nu_diag,1000) ' xndt_dyn = ', xndt_dyn
write(nu_diag,1020) ' ndte = ', ndte
write(nu_diag,1010) ' evp_damping = ', &
evp_damping
write(nu_diag,*) ' yield_curve = ', &
trim(yield_curve)
write(nu_diag,1020) ' kstrength = ', kstrength
write(nu_diag,1020) ' krdg_partic = ', &
krdg_partic
write(nu_diag,1020) ' krdg_redist = ', &
krdg_redist
write(nu_diag,1030) ' advection = ', &
trim(advection)
write(nu_diag,1030) ' shortwave = ', &
trim(shortwave)
write(nu_diag,1030) ' albedo_type = ', &
trim(albedo_type)
write(nu_diag,1000) ' R_ice = ', R_ice
write(nu_diag,1000) ' R_pnd = ', R_pnd
write(nu_diag,1000) ' R_snw = ', R_snw
write(nu_diag,1000) ' dT_mlt_in = ', dT_mlt_in
write(nu_diag,1000) ' rsnw_melt_in = ', rsnw_melt_in
write(nu_diag,1000) ' albicev = ', albicev
write(nu_diag,1000) ' albicei = ', albicei
write(nu_diag,1000) ' albsnowv = ', albsnowv
write(nu_diag,1000) ' albsnowi = ', albsnowi
write(nu_diag,1010) ' heat_capacity = ', &
heat_capacity
write(nu_diag,1030) ' atmbndy = ', &
trim(atmbndy)
write(nu_diag,1020) ' fyear_init = ', &
fyear_init
write(nu_diag,1020) ' ycycle = ', ycycle
write(nu_diag,*) ' atm_data_type = ', &
trim(atm_data_type)
write(nu_diag,1010) ' calc_strair = ', calc_strair
write(nu_diag,1010) ' calc_Tsfc = ', calc_Tsfc
write(nu_diag,*) ' Tfrzpt = ', trim(Tfrzpt)
write(nu_diag,1010) ' update_ocn_f = ', update_ocn_f
if (trim(atm_data_type) /= 'default') then
write(nu_diag,*) ' atm_data_dir = ', &
trim(atm_data_dir)
write(nu_diag,*) ' precip_units = ', &
trim(precip_units)
endif
write(nu_diag,1010) ' oceanmixed_ice = ', &
oceanmixed_ice
write(nu_diag,*) ' sss_data_type = ', &
trim(sss_data_type)
write(nu_diag,*) ' sst_data_type = ', &
trim(sst_data_type)
if (trim(sss_data_type) /= 'default' .or. &
trim(sst_data_type) /= 'default') then
write(nu_diag,*) ' ocn_data_dir = ', &
trim(ocn_data_dir)
endif
if (trim(sss_data_type) == 'ncar' .or. &
trim(sst_data_type) == 'ncar') then
write(nu_diag,*) ' oceanmixed_file = ', &
trim(oceanmixed_file)
endif
#ifdef coupled
if( oceanmixed_ice ) then
write (nu_diag,*) 'WARNING WARNING WARNING WARNING '
write (nu_diag,*) '*Coupled and oceanmixed flags are *'
write (nu_diag,*) '*BOTH ON. Ocean data received from*'
write (nu_diag,*) '*coupler will be altered by mixed *'
write (nu_diag,*) '*layer routine! *'
write (nu_diag,*) ' '
endif
#endif
write (nu_diag,*) ' '
write (nu_diag,'(a30,2f8.2)') 'Diagnostic point 1: lat, lon =', &
latpnt(1), lonpnt(1)
write (nu_diag,'(a30,2f8.2)') 'Diagnostic point 2: lat, lon =', &
latpnt(2), lonpnt(2)
! tracers
write(nu_diag,1010) ' tr_iage = ', tr_iage
write(nu_diag,1010) ' restart_age = ', restart_age
write(nu_diag,1010) ' tr_FY = ', tr_FY
write(nu_diag,1010) ' restart_FY = ', restart_FY
write(nu_diag,1010) ' tr_lvl = ', tr_lvl
write(nu_diag,1010) ' restart_lvl = ', restart_lvl
write(nu_diag,1010) ' tr_pond = ', tr_pond
write(nu_diag,1010) ' restart_pond = ', restart_pond
write(nu_diag,1010) ' tr_aero = ', tr_aero !MH
write(nu_diag,1010) ' restart_aero = ', restart_aero !MH
nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1
ntrcr = 1 ! count tracers, starting with Tsfc = 1
if (tr_iage) then
nt_iage = ntrcr + 1
ntrcr = ntrcr + 1
else
nt_iage = max_ntrcr
endif
if (tr_FY) then
nt_FY = ntrcr + 1
ntrcr = ntrcr + 1
else
nt_FY = max_ntrcr
endif
if (tr_lvl) then
nt_alvl = ntrcr + 1
ntrcr = ntrcr + 1
nt_vlvl = ntrcr + 1
ntrcr = ntrcr + 1
else
nt_alvl = max_ntrcr
nt_vlvl = max_ntrcr
endif
if (tr_pond) then
nt_volpn = ntrcr + 1
ntrcr = ntrcr + 1
else
nt_volpn = max_ntrcr
endif
if (tr_aero) then
nt_aero = ntrcr + 1
ntrcr = ntrcr + n_aero*4 !MH 2 for snow soot and 2 for ice
!MH for multiple (n_aero) aerosols
else
nt_aero = max_ntrcr
endif
if (ntrcr > max_ntrcr) then
write(nu_diag,*) 'max_ntrcr < number of namelist tracers'
call abort_ice('max_ntrcr < number of namelist tracers')
endif
1000 format (a30,2x,f9.2) ! a30 to align formatted, unformatted statements
1010 format (a30,2x,l6) ! logical
1020 format (a30,2x,i6) ! integer
1030 format (a30, a8) ! character
1040 format (a30,2x,6i6) ! integer
1050 format (a30, 6a8) ! character
write (nu_diag,*) ' '
if (grid_type /= 'displaced_pole' .and. &
grid_type /= 'tripole' .and. &
grid_type /= 'column' .and. &
grid_type /= 'rectangular' .and. &
grid_type /= 'panarctic' .and. &
grid_type /= 'latlon' ) then
call abort_ice('ice_init: unknown grid_type')
endif
endif ! my_task = master_task
call broadcast_scalar(ntrcr, master_task)
call broadcast_scalar(nt_Tsfc, master_task)
call broadcast_scalar(nt_iage, master_task)
call broadcast_scalar(nt_FY, master_task)
call broadcast_scalar(nt_alvl, master_task)
call broadcast_scalar(nt_vlvl, master_task)
call broadcast_scalar(nt_volpn, master_task)
call broadcast_scalar(nt_aero, master_task)
end subroutine input_data
!=======================================================================
!BOP
!
! !IROUTINE: init_state - initialize ice state variables
!
! !INTERFACE:
!
subroutine init_state 2,35
!
! !DESCRIPTION:
!
! Initialize state for the itd model
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz, UW
! William H. Lipscomb, LANL
!
! !USES:
!
use ice_blocks
use ice_domain
use ice_flux, only: sst, Tf, Tair
use ice_grid
use ice_state
use ice_itd
use ice_exit
use ice_therm_vertical, only: heat_capacity
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
integer (kind=int_kind) :: &
i, j , & ! horizontal indices
it, n , & ! tracer index
iblk ! block index
!-----------------------------------------------------------------
! Check number of layers in ice and snow.
!-----------------------------------------------------------------
if (my_task == master_task) then
if (nilyr < 1) then
write (nu_diag,*) 'nilyr =', nilyr
write (nu_diag,*) 'Must have at least one ice layer'
call abort_ice ('ice_init: Not enough ice layers')
endif
if (nslyr < 1) then
write (nu_diag,*) 'nslyr =', nslyr
write (nu_diag,*) 'Must have at least one snow layer'
call abort_ice('ice_init: Not enough snow layers')
endif
if (.not.heat_capacity) then
write (nu_diag,*) 'WARNING - Zero-layer thermodynamics'
if (nilyr > 1) then
write (nu_diag,*) 'nilyr =', nilyr
write (nu_diag,*) &
'Must have nilyr = 1 if heat_capacity = F'
call abort_ice('ice_init: Too many ice layers')
endif
if (nslyr > 1) then
write (nu_diag,*) 'nslyr =', nslyr
write (nu_diag,*) &
'Must have nslyr = 1 if heat_capacity = F'
call abort_ice('ice_init: Too many snow layers')
endif
endif ! heat_capacity = F
endif ! my_task
!-----------------------------------------------------------------
! Set tracer types
!-----------------------------------------------------------------
trcr_depend(nt_Tsfc) = 0 ! ice/snow surface temperature
if (tr_iage) trcr_depend(nt_iage) = 1 ! volume-weighted ice age
if (tr_FY) trcr_depend(nt_FY) = 0 ! area-weighted FY conc
if (tr_lvl) trcr_depend(nt_alvl) = 0 ! level ice area
if (tr_lvl) trcr_depend(nt_vlvl) = 1 ! level ice volume
if (tr_pond) trcr_depend(nt_volpn) = 0 ! melt pond volume
if (tr_aero) then
do n=1,n_aero
trcr_depend(nt_aero+(n-1)*4 :nt_aero+(n-1)*4+1) = 2 ! snow volume-weighted MH
trcr_depend(nt_aero+(n-1)*4+2:nt_aero+(n-1)*4+3) = 1 ! volume-weighted MH
enddo
endif
!$OMP PARALLEL DO PRIVATE(iblk,it)
do iblk = 1, nblocks
!-----------------------------------------------------------------
! Set state variables
!-----------------------------------------------------------------
call set_state_var (nx_block, ny_block, &
tmask(:,:, iblk), &
ULON (:,:, iblk), ULAT (:,:, iblk), &
Tair (:,:, iblk), sst (:,:, iblk), &
Tf (:,:, iblk), trcr_depend, &
aicen(:,:, :,iblk), trcrn(:,:,:,:,iblk), &
vicen(:,:, :,iblk), vsnon(:,:, :,iblk), &
eicen(:,:, :,iblk), esnon(:,:, :,iblk))
!-----------------------------------------------------------------
! compute aggregate ice state and open water area
!-----------------------------------------------------------------
aice(:,:,iblk) = c0
vice(:,:,iblk) = c0
vsno(:,:,iblk) = c0
eice(:,:,iblk) = c0
esno(:,:,iblk) = c0
do it = 1, max_ntrcr
trcr(:,:,it,iblk) = c0
enddo
call aggregate (nx_block, ny_block, &
aicen(:,:,:,iblk), &
trcrn(:,:,:,:,iblk), &
vicen(:,:,:,iblk), &
vsnon(:,:,:,iblk), &
eicen(:,:,:,iblk), &
esnon(:,:,:,iblk), &
aice (:,:, iblk), &
trcr (:,:,:,iblk), &
vice (:,:, iblk), &
vsno (:,:, iblk), &
eice (:,:, iblk), &
esno (:,:, iblk), &
aice0(:,:, iblk), &
tmask(:,:, iblk), &
ntrcr, &
trcr_depend)
aice_init(:,:,iblk) = aice(:,:,iblk)
enddo ! iblk
!$OMP END PARALLEL DO
!-----------------------------------------------------------------
! ghost cell updates
!-----------------------------------------------------------------
call bound_state (aicen, trcrn, &
vicen, vsnon, &
eicen, esnon)
end subroutine init_state
!=======================================================================
!BOP
!
! !IROUTINE: set_state_var - initialize single-category state variables
!
! !INTERFACE:
!
subroutine set_state_var (nx_block, ny_block, & 1,5
tmask, ULON, ULAT, &
Tair, sst, &
Tf, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon)
!
! !DESCRIPTION:
!
! Initialize state in each ice thickness category
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz
! William H. Lipscomb, LANL
!
! !USES:
!
use ice_state, only: nt_Tsfc
use ice_therm_vertical, only: heat_capacity, calc_Tsfc, Tmlt
use ice_itd, only: ilyr1, slyr1, hin_max
use ice_grid, only: grid_type
use ice_restart, only: ice_ic
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ! block dimensions
logical (kind=log_kind), dimension (nx_block,ny_block), &
intent(in) :: &
tmask ! true for ice/ocean cells
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
ULON , & ! latitude of velocity pts (radians)
ULAT ! latitude of velocity pts (radians)
real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
Tair , & ! air temperature (K)
Tf , & ! freezing temperature (C)
sst ! sea surface temperature (C)
integer (kind=int_kind), dimension (max_ntrcr), intent(inout) :: &
trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(out) :: &
aicen , & ! concentration of ice
vicen , & ! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(out) :: &
trcrn ! ice tracers
! 1: surface temperature of ice/snow (C)
real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), &
intent(out) :: &
eicen ! energy of melting for each ice layer (J/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), &
intent(out) :: &
esnon ! energy of melting for each ice layer (J/m^2)
!
!EOP
!
integer (kind=int_kind) :: &
i, j , & ! horizontal indices
ij , & ! horizontal index, combines i and j loops
k , & ! ice layer index
n , & ! thickness category index
it , & ! tracer index
icells ! number of cells initialized with ice
integer (kind=int_kind), dimension(nx_block*ny_block) :: &
indxi, indxj ! compressed indices for cells with aicen > puny
real (kind=dbl_kind) :: &
slope, Ti, sum, hbar, &
ainit(ncat), &
hinit(ncat)
real (kind=dbl_kind), parameter :: &
hsno_init = 0.20_dbl_kind , & ! initial snow thickness (m)
edge_init_nh = 70._dbl_kind, & ! initial ice edge, N.Hem. (deg)
edge_init_sh = -60._dbl_kind ! initial ice edge, S.Hem. (deg)
indxi(:) = 0
indxj(:) = 0
! Initialize state variables.
! If restarting, these values are overwritten.
do n = 1, ncat
do j = 1, ny_block
do i = 1, nx_block
aicen(i,j,n) = c0
vicen(i,j,n) = c0
vsnon(i,j,n) = c0
trcrn(i,j,nt_Tsfc,n) = Tf(i,j) ! surface temperature
if (max_ntrcr >= 2) then
do it = 2, max_ntrcr
trcrn(i,j,it,n) = c0
enddo
endif
enddo
enddo
enddo
eicen(:,:,:) = c0
esnon(:,:,:) = c0
if (trim(ice_ic) == 'default') then
!-----------------------------------------------------------------
! Place ice where ocean surface is cold.
! Note: If SST is not read from a file, then the ocean is assumed
! to be at its freezing point everywhere, and ice will
! extend to the prescribed edges.
!-----------------------------------------------------------------
! initial category areas in cells with ice
hbar = c3 ! initial ice thickness with greatest area
! Note: the resulting average ice thickness
! tends to be less than hbar due to the
! nonlinear distribution of ice thicknesses
sum = c0
do n = 1, ncat
if (n < ncat) then
hinit(n) = p5*(hin_max(n-1) + hin_max(n)) ! m
else ! n=ncat
hinit(n) = (hin_max(n-1) + c1) ! m
endif
! parabola, max at h=hbar, zero at h=0, 2*hbar
ainit(n) = max(c0, (c2*hbar*hinit(n) - hinit(n)**2))
sum = sum + ainit(n)
enddo
do n = 1, ncat
ainit(n) = ainit(n) / (sum + puny/ncat) ! normalize
enddo
if (trim(grid_type) == 'rectangular') then
! place ice on left side of domain
icells = 0
do j = 1, ny_block
do i = 1, nx_block
if (tmask(i,j)) then
if (ULON(i,j) < -50./rad_to_deg) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif ! ULON
endif ! tmask
enddo ! i
enddo ! j
else
! place ice at high latitudes where ocean sfc is cold
icells = 0
do j = 1, ny_block
do i = 1, nx_block
if (tmask(i,j)) then
! place ice in high latitudes where ocean sfc is cold
if ( (sst (i,j) <= Tf(i,j)+p2) .and. &
(ULAT(i,j) < edge_init_sh/rad_to_deg .or. &
ULAT(i,j) > edge_init_nh/rad_to_deg) ) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif ! cold surface
endif ! tmask
enddo ! i
enddo ! j
endif ! rectgrid
do n = 1, ncat
! ice volume, snow volume
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
aicen(i,j,n) = ainit(n)
vicen(i,j,n) = hinit(n) * ainit(n) ! m
vsnon(i,j,n) =min(aicen(i,j,n)*hsno_init,p2*vicen(i,j,n))
enddo ! ij
! surface temperature
if (calc_Tsfc) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
trcrn(i,j,nt_Tsfc,n) = min(Tsmelt, Tair(i,j) - Tffresh) !deg C
enddo
else ! Tsfc is not calculated by the ice model
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
trcrn(i,j,nt_Tsfc,n) = Tf(i,j) ! not used
enddo
endif ! calc_Tsfc
! other tracers (none at present)
if (heat_capacity) then
! ice energy
do k = 1, nilyr
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
! assume linear temp profile and compute enthalpy
slope = Tf(i,j) - trcrn(i,j,nt_Tsfc,n)
Ti = trcrn(i,j,nt_Tsfc,n) &
+ slope*(real(k,kind=dbl_kind)-p5) &
/real(nilyr,kind=dbl_kind)
eicen(i,j,ilyr1(n)+k-1) = &
-(rhoi * (cp_ice*(Tmlt(k)-Ti) &
+ Lfresh*(c1-Tmlt(k)/Ti) - cp_ocn*Tmlt(k))) &
* vicen(i,j,n)/real(nilyr,kind=dbl_kind)
enddo ! ij
enddo ! nilyr
! snow energy
do k = 1, nslyr
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
Ti = min(c0, trcrn(i,j,nt_Tsfc,n))
esnon(i,j,slyr1(n)+k-1) = -rhos*(Lfresh - cp_ice*Ti) &
*vsnon(i,j,n) &
/real(nslyr,kind=dbl_kind)
enddo ! ij
enddo ! nslyr
else ! one layer with zero heat capacity
! ice energy
k = 1
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
eicen(i,j,ilyr1(n)+k-1) = &
- rhoi * Lfresh * vicen(i,j,n)
enddo ! ij
! snow energy
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
esnon(i,j,slyr1(n)+k-1) = &
- rhos * Lfresh * vsnon(i,j,n)
enddo ! ij
endif ! heat_capacity
enddo ! ncat
endif ! ice_ic
end subroutine set_state_var
!=======================================================================
end module ice_init
!=======================================================================