diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index f1c0b8c19..26f282ea8 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -297,12 +297,10 @@ subroutine input_data dumpfreq='y' ! restart frequency option dumpfreq_n = 1 ! restart frequency dump_last = .false. ! write restart on last time step - 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_ext = .false. ! if true, read/write ghost cells restart_coszen = .false. ! if true, read/write coszen - use_restart_time = .true. ! if true, use time info written in file pointer_file = 'ice.restart_file' restart_format = 'default' ! restart file format lcdf64 = .false. ! 64 bit offset for netCDF @@ -452,6 +450,8 @@ subroutine input_data #ifndef CESMCOUPLED runid = 'unknown' ! run ID used in CESM and for machine 'bering' runtype = 'initial' ! run type: 'initial', 'continue' + restart = .false. ! if true, read restart files for initialization + use_restart_time = .true. ! if true, use time info written in file #endif ! extra tracers @@ -1775,7 +1775,8 @@ subroutine input_data grid_type /= 'rectangular' .and. & grid_type /= 'cpom_grid' .and. & grid_type /= 'regional' .and. & - grid_type /= 'latlon' ) then + grid_type /= 'latlon' .and. & + grid_type /= 'setmask' ) then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: unknown grid_type=',trim(grid_type) abort_list = trim(abort_list)//":20" endif diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 2d660af81..1c7937b5d 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -40,7 +40,7 @@ module ice_grid private public :: init_grid1, init_grid2, & t2ugrid_vector, u2tgrid_vector, & - to_ugrid, to_tgrid, alloc_grid + to_ugrid, to_tgrid, alloc_grid, makemask character (len=char_len_long), public :: & grid_format , & ! file format ('bin'=binary or 'nc'=netcdf) diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index cb70c9b4a..cfca994c3 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -1,445 +1,451 @@ -!======================================================================= -! -! This module contains the CICE initialization routine that sets model -! parameters and initializes the grid and CICE state variables. -! -! authors Elizabeth C. Hunke, LANL -! William H. Lipscomb, LANL -! Philip W. Jones, LANL -! -! 2006: Converted to free form source (F90) by Elizabeth Hunke -! 2008: E. Hunke moved ESMF code to its own driver - - module CICE_InitMod - - use ice_kinds_mod - use ice_exit, only: abort_ice - use ice_fileunits, only: init_fileunits, nu_diag - use icepack_intfc, only: icepack_aggregate - use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist - use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave - use icepack_intfc, only: icepack_configure - use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted - use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, & - icepack_query_tracer_indices, icepack_query_tracer_sizes - - implicit none - private - public :: cice_init +module CICE_InitMod -!======================================================================= + ! Initialize CICE model. + + use ice_kinds_mod + use ice_exit , only: abort_ice + use ice_fileunits, only: init_fileunits, nu_diag + use icepack_intfc, only: icepack_aggregate + use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist + use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_configure + use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted + use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags + use icepack_intfc, only: icepack_query_tracer_indices, icepack_query_tracer_sizes + + implicit none + private + public :: cice_init1 + public :: cice_init2 - contains + private :: init_restart !======================================================================= -! -! Initialize CICE model. - - subroutine cice_init - - ! Initialize the basic state, grid and all necessary parameters for - ! running the CICE model. - - use ice_arrays_column, only: hin_max, c_hi_range, alloc_arrays_column - use ice_arrays_column, only: floe_rad_l, floe_rad_c, & - floe_binwidth, c_fsd_range - use ice_state, only: alloc_state - use ice_flux_bgc, only: alloc_flux_bgc - use ice_calendar, only: dt, dt_dyn, time, istep, istep1, write_ic, & - init_calendar, calendar - use ice_communicate, only: my_task, master_task - use ice_diagnostics, only: init_diags - use ice_domain, only: init_domain_blocks - use ice_domain_size, only: ncat, nfsd - use ice_dyn_eap, only: init_eap, alloc_dyn_eap - use ice_dyn_shared, only: kdyn, init_dyn, alloc_dyn_shared - use ice_flux, only: init_coupler_flux, init_history_therm, & - init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux - use ice_forcing, only: init_forcing_ocn - use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & - faero_default, faero_optics, alloc_forcing_bgc, fiso_default - use ice_grid, only: init_grid1, init_grid2, alloc_grid - use ice_history, only: init_hist, accum_hist - use ice_restart_shared, only: restart, runtype - use ice_init, only: input_data, init_state - use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers - use ice_kinds_mod - use ice_restoring, only: ice_HaloRestore_init - use ice_timers, only: timer_total, init_ice_timers, ice_timer_start - use ice_transport_driver, only: init_transport - - logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, & - tr_iso, tr_fsd, wave_spec - character(len=*), parameter :: subname = '(cice_init)' - - call init_fileunits ! unit numbers - - call icepack_configure() ! initialize icepack - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - call input_data ! namelist variables - call input_zbgc ! vertical biogeochemistry namelist - call count_tracers ! count tracers - - call init_domain_blocks ! set up block decomposition - call init_grid1 ! domain distribution - call alloc_grid ! allocate grid arrays - call alloc_arrays_column ! allocate column arrays - call alloc_state ! allocate state arrays - call alloc_dyn_shared ! allocate dyn shared arrays - call alloc_flux_bgc ! allocate flux_bgc arrays - call alloc_flux ! allocate flux arrays - call init_ice_timers ! initialize all timers - call ice_timer_start(timer_total) ! start timing entire run - call init_grid2 ! grid variables - call init_zbgc ! vertical biogeochemistry initialization - call init_calendar ! initialize some calendar stuff - call init_hist (dt) ! initialize output history file - - call init_dyn (dt_dyn) ! define dynamics parameters, variables - if (kdyn == 2) then - call alloc_dyn_eap ! allocate dyn_eap arrays - call init_eap ! define eap dynamics parameters, variables - else if (kdyn == 3) then - call init_vp ! define vp dynamics parameters, variables - endif - - call init_coupler_flux ! initialize fluxes exchanged with coupler - call init_thermo_vertical ! initialize vertical thermodynamics - - call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution - if (my_task == master_task) then - call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output - endif - - call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution +contains +!======================================================================= + + subroutine cice_init1() + + ! Initialize the basic state, grid and all necessary parameters for + ! running the CICE model. + + use ice_init , only: input_data + use ice_init_column , only: input_zbgc, count_tracers + use ice_grid , only: init_grid1, alloc_grid + use ice_domain , only: init_domain_blocks + use ice_arrays_column , only: alloc_arrays_column + use ice_state , only: alloc_state + use ice_dyn_shared , only: alloc_dyn_shared + use ice_flux_bgc , only: alloc_flux_bgc + use ice_flux , only: alloc_flux + use ice_timers , only: timer_total, init_ice_timers, ice_timer_start + + character(len=*), parameter :: subname = '(cice_init1)' + !---------------------------------------------------- + + call init_fileunits ! unit numbers + call icepack_configure() ! initialize icepack + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + call input_data ! namelist variables + call input_zbgc ! vertical biogeochemistry namelist + call count_tracers ! count tracers + + call init_domain_blocks ! set up block decomposition + call init_grid1 ! domain distribution + call alloc_grid ! allocate grid arrays + call alloc_arrays_column ! allocate column arrays + call alloc_state ! allocate state arrays + call alloc_dyn_shared ! allocate dyn shared arrays + call alloc_flux_bgc ! allocate flux_bgc arrays + call alloc_flux ! allocate flux arrays + call init_ice_timers ! initialize all timers + call ice_timer_start(timer_total) ! start timing entire run + + end subroutine cice_init1 + + !======================================================================= + subroutine cice_init2() + + ! Initialize the basic state, and all necessary parameters for + ! running the CICE model. + + use ice_arrays_column , only: hin_max, c_hi_range + use ice_arrays_column , only: floe_rad_l, floe_rad_c, floe_binwidth, c_fsd_range + use ice_calendar , only: dt, dt_dyn, istep, istep1, write_ic, init_calendar, calendar + use ice_communicate , only: my_task, master_task + use ice_diagnostics , only: init_diags + use ice_domain_size , only: ncat, nfsd + use ice_dyn_eap , only: init_eap, alloc_dyn_eap + use ice_dyn_shared , only: kdyn, init_dyn + use ice_dyn_vp , only: init_vp + use ice_flux , only: init_coupler_flux, init_history_therm + use ice_flux , only: init_history_dyn, init_flux_atm, init_flux_ocn + use ice_forcing , only: init_forcing_ocn + use ice_forcing_bgc , only: get_forcing_bgc, get_atm_bgc + use ice_forcing_bgc , only: faero_default, faero_optics, alloc_forcing_bgc, fiso_default + use ice_history , only: init_hist, accum_hist + use ice_restart_shared , only: restart, runtype + use ice_init , only: input_data, init_state + use ice_init_column , only: init_thermo_vertical, init_shortwave, init_zbgc + use ice_restoring , only: ice_HaloRestore_init + use ice_timers , only: timer_total, init_ice_timers, ice_timer_start + use ice_transport_driver , only: init_transport + + logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers + logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec + character(len=*), parameter :: subname = '(cice_init2)' + !---------------------------------------------------- + + call init_zbgc ! vertical biogeochemistry initialization + call init_calendar ! initialize some calendar stuff + call init_hist (dt) ! initialize output history file + + call init_dyn (dt_dyn) ! define dynamics parameters, variables + if (kdyn == 2) then + call alloc_dyn_eap ! allocate dyn_eap arrays + call init_eap ! define eap dynamics parameters, variables + else if (kdyn == 3) then + call init_vp ! define vp dynamics parameters, variables + endif + + call init_coupler_flux ! initialize fluxes exchanged with coupler + call init_thermo_vertical ! initialize vertical thermodynamics + + call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution + if (my_task == master_task) then + call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output + endif + + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution floe_rad_l, & ! fsd size lower bound in m (radius) floe_rad_c, & ! fsd size bin centre in m (radius) floe_binwidth, & ! fsd size bin width in m (radius) c_fsd_range, & ! string for history output write_diags=(my_task == master_task)) ! write diag on master only - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - call calendar(time) ! determine the initial date - - ! TODO: - why is this being called when you are using CMEPS? - call init_forcing_ocn(dt) ! initialize sss and sst from data - - call init_state ! initialize the ice state - call init_transport ! initialize horizontal transport - call ice_HaloRestore_init ! restored boundary conditions - - call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & - wave_spec_out=wave_spec) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays - - call init_restart ! initialize restart variables - call init_diags ! initialize diagnostic output points - call init_history_therm ! initialize thermo history variables - call init_history_dyn ! initialize dynamic history variables - - call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) - call icepack_query_tracer_flags(tr_iso_out=tr_iso) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (tr_aero .or. tr_zaero) then - call faero_optics !initialize aerosol optical property tables - end if - - ! Initialize shortwave components using swdn from previous timestep - ! if restarting. These components will be scaled to current forcing - ! in prep_radiation. - - if (trim(runtype) == 'continue' .or. restart) then - call init_shortwave ! initialize radiative transfer - end if - - !-------------------------------------------------------------------- - ! coupler communication or forcing data initialization - !-------------------------------------------------------------------- - - if (z_tracers) call get_atm_bgc ! biogeochemistry - - if (runtype == 'initial' .and. .not. restart) then - call init_shortwave ! initialize radiative transfer using current swdn - end if - - call init_flux_atm ! initialize atmosphere fluxes sent to coupler - call init_flux_ocn ! initialize ocean fluxes sent to coupler - - end subroutine cice_init - -!======================================================================= - - subroutine init_restart - - use ice_arrays_column, only: dhsn - use ice_blocks, only: nx_block, ny_block - use ice_calendar, only: time, calendar - use ice_constants, only: c0 - use ice_domain, only: nblocks - use ice_domain_size, only: ncat, n_iso, n_aero, nfsd - use ice_dyn_eap, only: read_restart_eap - use ice_dyn_shared, only: kdyn - use ice_grid, only: tmask - use ice_init, only: ice_ic - use ice_init_column, only: init_age, init_FY, init_lvl, & - init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & - init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd - use ice_restart_column, only: restart_age, read_restart_age, & - restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & - restart_pond_cesm, read_restart_pond_cesm, & - restart_pond_lvl, read_restart_pond_lvl, & - restart_pond_topo, read_restart_pond_topo, & - restart_fsd, read_restart_fsd, & - restart_iso, read_restart_iso, & - restart_aero, read_restart_aero, & - restart_hbrine, read_restart_hbrine, & - restart_zsal, restart_bgc - use ice_restart_driver, only: restartfile - use ice_restart_shared, only: runtype, restart - use ice_state ! almost everything - - integer(kind=int_kind) :: & + call calendar() ! determine the initial date + + !TODO: - why is this being called when you are using CMEPS? + call init_forcing_ocn(dt) ! initialize sss and sst from data + + call init_state ! initialize the ice state + call init_transport ! initialize horizontal transport + call ice_HaloRestore_init ! restored boundary conditions + + call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & + wave_spec_out=wave_spec) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays + + call init_restart ! initialize restart variables + call init_diags ! initialize diagnostic output points + call init_history_therm ! initialize thermo history variables + call init_history_dyn ! initialize dynamic history variables + + call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) + call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_aero .or. tr_zaero) then + call faero_optics !initialize aerosol optical property tables + end if + + ! Initialize shortwave components using swdn from previous timestep + ! if restarting. These components will be scaled to current forcing + ! in prep_radiation. + + if (trim(runtype) == 'continue' .or. restart) then + call init_shortwave ! initialize radiative transfer + end if + + !-------------------------------------------------------------------- + ! coupler communication or forcing data initialization + !-------------------------------------------------------------------- + + if (z_tracers) call get_atm_bgc ! biogeochemistry + + if (runtype == 'initial' .and. .not. restart) then + call init_shortwave ! initialize radiative transfer using current swdn + end if + + call init_flux_atm ! initialize atmosphere fluxes sent to coupler + call init_flux_ocn ! initialize ocean fluxes sent to coupler + + end subroutine cice_init2 + + !======================================================================= + + subroutine init_restart() + + use ice_arrays_column, only: dhsn + use ice_blocks, only: nx_block, ny_block + use ice_calendar, only: calendar + use ice_constants, only: c0 + use ice_domain, only: nblocks + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_dyn_eap, only: read_restart_eap + use ice_dyn_shared, only: kdyn + use ice_grid, only: tmask + use ice_init, only: ice_ic + use ice_init_column, only: init_age, init_FY, init_lvl, & + init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & + init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd + use ice_restart_column, only: restart_age, read_restart_age, & + restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & + restart_pond_cesm, read_restart_pond_cesm, & + restart_pond_lvl, read_restart_pond_lvl, & + restart_pond_topo, read_restart_pond_topo, & + restart_fsd, read_restart_fsd, & + restart_iso, read_restart_iso, & + restart_aero, read_restart_aero, & + restart_hbrine, read_restart_hbrine, & + restart_zsal, restart_bgc + use ice_restart_driver, only: restartfile + use ice_restart_shared, only: runtype, restart + use ice_state ! almost everything + + integer(kind=int_kind) :: & i, j , & ! horizontal indices iblk ! block index - logical(kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & - tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & - skl_bgc, z_tracers, solve_zsal - integer(kind=int_kind) :: & - ntrcr - integer(kind=int_kind) :: & - nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & - nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice - - character(len=*), parameter :: subname = '(init_restart)' - - call icepack_query_tracer_sizes(ntrcr_out=ntrcr) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - call icepack_query_parameters(skl_bgc_out=skl_bgc, & - z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) - call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & - tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & - tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) - call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & - nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & - nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & - nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + logical(kind=log_kind) :: & + tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & + tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + skl_bgc, z_tracers, solve_zsal + integer(kind=int_kind) :: & + ntrcr + integer(kind=int_kind) :: & + nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice + + character(len=*), parameter :: subname = '(init_restart)' + !---------------------------------------------------- + + call icepack_query_tracer_sizes(ntrcr_out=ntrcr) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (trim(runtype) == 'continue') then - ! start from core restart file - call restartfile() ! given by pointer in ice_in - call calendar(time) ! update time parameters - if (kdyn == 2) call read_restart_eap ! EAP - else if (restart) then ! ice_ic = core restart file - call restartfile (ice_ic) ! or 'default' or 'none' - !!! uncomment to create netcdf - ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file - !!! uncomment if EAP restart data exists - ! if (kdyn == 2) call read_restart_eap - endif - - ! tracers - ! ice age tracer - if (tr_iage) then - if (trim(runtype) == 'continue') & - restart_age = .true. - if (restart_age) then - call read_restart_age - else - do iblk = 1, nblocks - call init_age(trcrn(:,:,nt_iage,:,iblk)) - enddo ! iblk - endif - endif - ! first-year area tracer - if (tr_FY) then - if (trim(runtype) == 'continue') restart_FY = .true. - if (restart_FY) then - call read_restart_FY - else - do iblk = 1, nblocks - call init_FY(trcrn(:,:,nt_FY,:,iblk)) - enddo ! iblk - endif - endif - ! level ice tracer - if (tr_lvl) then - if (trim(runtype) == 'continue') restart_lvl = .true. - if (restart_lvl) then - call read_restart_lvl - else - do iblk = 1, nblocks - call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & - trcrn(:,:,nt_vlvl,:,iblk)) - enddo ! iblk - endif - endif - ! CESM melt ponds - if (tr_pond_cesm) then - if (trim(runtype) == 'continue') & - restart_pond_cesm = .true. - if (restart_pond_cesm) then - call read_restart_pond_cesm - else - do iblk = 1, nblocks - call init_meltponds_cesm(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk)) - enddo ! iblk - endif - endif - ! level-ice melt ponds - if (tr_pond_lvl) then - if (trim(runtype) == 'continue') & - restart_pond_lvl = .true. - if (restart_pond_lvl) then - call read_restart_pond_lvl - else - do iblk = 1, nblocks - call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk), & - trcrn(:,:,nt_ipnd,:,iblk), & - dhsn(:,:,:,iblk)) - enddo ! iblk - endif - endif - ! topographic melt ponds - if (tr_pond_topo) then - if (trim(runtype) == 'continue') & - restart_pond_topo = .true. - if (restart_pond_topo) then - call read_restart_pond_topo - else - do iblk = 1, nblocks - call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk), & - trcrn(:,:,nt_ipnd,:,iblk)) - enddo ! iblk - endif ! .not. restart_pond - endif - ! floe size distribution - if (tr_fsd) then - if (trim(runtype) == 'continue') restart_fsd = .true. - if (restart_fsd) then - call read_restart_fsd - else - call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) - endif - endif - - ! isotopes - if (tr_iso) then - if (trim(runtype) == 'continue') restart_iso = .true. - if (restart_iso) then - call read_restart_iso - else - do iblk = 1, nblocks - call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & - trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) - enddo ! iblk - endif - endif - - if (tr_aero) then ! ice aerosol - if (trim(runtype) == 'continue') restart_aero = .true. - if (restart_aero) then - call read_restart_aero - else - do iblk = 1, nblocks - call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) - enddo ! iblk - endif ! .not. restart_aero - endif - - if (trim(runtype) == 'continue') then - if (tr_brine) & - restart_hbrine = .true. - if (solve_zsal) & - restart_zsal = .true. - if (skl_bgc .or. z_tracers) & - restart_bgc = .true. - endif - - if (tr_brine .or. skl_bgc) then ! brine height tracer - call init_hbrine - if (tr_brine .and. restart_hbrine) call read_restart_hbrine - endif - - if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry - if (tr_fsd) then - write (nu_diag,*) 'FSD implementation incomplete for use with BGC' - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + call icepack_query_parameters(skl_bgc_out=skl_bgc, & + z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & + tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (trim(runtype) == 'continue') then + ! start from core restart file + call restartfile() ! given by pointer in ice_in + call calendar() ! update time parameters + if (kdyn == 2) call read_restart_eap ! EAP + else if (restart) then ! ice_ic = core restart file + call restartfile (ice_ic) ! or 'default' or 'none' +!!! uncomment to create netcdf + ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file +!!! uncomment if EAP restart data exists + ! if (kdyn == 2) call read_restart_eap + endif + + ! tracers + ! ice age tracer + if (tr_iage) then + if (trim(runtype) == 'continue') & + restart_age = .true. + if (restart_age) then + call read_restart_age + else + do iblk = 1, nblocks + call init_age(trcrn(:,:,nt_iage,:,iblk)) + enddo ! iblk + endif + endif + ! first-year area tracer + if (tr_FY) then + if (trim(runtype) == 'continue') restart_FY = .true. + if (restart_FY) then + call read_restart_FY + else + do iblk = 1, nblocks + call init_FY(trcrn(:,:,nt_FY,:,iblk)) + enddo ! iblk + endif + endif + ! level ice tracer + if (tr_lvl) then + if (trim(runtype) == 'continue') restart_lvl = .true. + if (restart_lvl) then + call read_restart_lvl + else + do iblk = 1, nblocks + call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & + trcrn(:,:,nt_vlvl,:,iblk)) + enddo ! iblk + endif + endif + ! CESM melt ponds + if (tr_pond_cesm) then + if (trim(runtype) == 'continue') & + restart_pond_cesm = .true. + if (restart_pond_cesm) then + call read_restart_pond_cesm + else + do iblk = 1, nblocks + call init_meltponds_cesm(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk)) + enddo ! iblk + endif + endif + ! level-ice melt ponds + if (tr_pond_lvl) then + if (trim(runtype) == 'continue') & + restart_pond_lvl = .true. + if (restart_pond_lvl) then + call read_restart_pond_lvl + else + do iblk = 1, nblocks + call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk), & + dhsn(:,:,:,iblk)) + enddo ! iblk + endif + endif + ! topographic melt ponds + if (tr_pond_topo) then + if (trim(runtype) == 'continue') & + restart_pond_topo = .true. + if (restart_pond_topo) then + call read_restart_pond_topo + else + do iblk = 1, nblocks + call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk)) + enddo ! iblk + endif ! .not. restart_pond + endif + ! floe size distribution + if (tr_fsd) then + if (trim(runtype) == 'continue') restart_fsd = .true. + if (restart_fsd) then + call read_restart_fsd + else + call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) + endif + endif + + ! isotopes + if (tr_iso) then + if (trim(runtype) == 'continue') restart_iso = .true. + if (restart_iso) then + call read_restart_iso + else + do iblk = 1, nblocks + call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & + trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) + enddo ! iblk + endif + endif + + if (tr_aero) then ! ice aerosol + if (trim(runtype) == 'continue') restart_aero = .true. + if (restart_aero) then + call read_restart_aero + else + do iblk = 1, nblocks + call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) + enddo ! iblk + endif ! .not. restart_aero + endif + + if (trim(runtype) == 'continue') then + if (tr_brine) & + restart_hbrine = .true. + if (solve_zsal) & + restart_zsal = .true. + if (skl_bgc .or. z_tracers) & + restart_bgc = .true. + endif + + if (tr_brine .or. skl_bgc) then ! brine height tracer + call init_hbrine + if (tr_brine .and. restart_hbrine) call read_restart_hbrine + endif + + if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry + if (tr_fsd) then + write (nu_diag,*) 'FSD implementation incomplete for use with BGC' + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - endif - call init_bgc - endif - - !----------------------------------------------------------------- - ! aggregate tracers - !----------------------------------------------------------------- - - !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - if (tmask(i,j,iblk)) then - call icepack_aggregate(ncat = ncat, & - aicen = aicen(i,j,:,iblk), & - trcrn = trcrn(i,j,:,:,iblk), & - vicen = vicen(i,j,:,iblk), & - vsnon = vsnon(i,j,:,iblk), & - aice = aice (i,j, iblk), & - trcr = trcr (i,j,:,iblk), & - vice = vice (i,j, iblk), & - vsno = vsno (i,j, iblk), & - aice0 = aice0(i,j, iblk), & - ntrcr = ntrcr, & - trcr_depend = trcr_depend, & - trcr_base = trcr_base, & - n_trcr_strata = n_trcr_strata, & - nt_strata = nt_strata) - else - ! tcraig, reset all tracer values on land to zero - trcrn(i,j,:,:,iblk) = c0 - endif - enddo - enddo - enddo - !$OMP END PARALLEL DO - - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + endif + call init_bgc + endif + + !----------------------------------------------------------------- + ! aggregate tracers + !----------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + if (tmask(i,j,iblk)) then + call icepack_aggregate(ncat = ncat, & + aicen = aicen(i,j,:,iblk), & + trcrn = trcrn(i,j,:,:,iblk), & + vicen = vicen(i,j,:,iblk), & + vsnon = vsnon(i,j,:,iblk), & + aice = aice (i,j, iblk), & + trcr = trcr (i,j,:,iblk), & + vice = vice (i,j, iblk), & + vsno = vsno (i,j, iblk), & + aice0 = aice0(i,j, iblk), & + ntrcr = ntrcr, & + trcr_depend = trcr_depend, & + trcr_base = trcr_base, & + n_trcr_strata = n_trcr_strata, & + nt_strata = nt_strata) + else + ! tcraig, reset all tracer values on land to zero + trcrn(i,j,:,:,iblk) = c0 + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - end subroutine init_restart + end subroutine init_restart -!======================================================================= + !======================================================================= - end module CICE_InitMod +end module CICE_InitMod !======================================================================= diff --git a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 index 9cb81e45f..81fa367c1 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_RunMod.F90 @@ -7,7 +7,7 @@ ! William H. Lipscomb, LANL ! ! 2006 ECH: moved exit timeLoop to prevent execution of unnecessary timestep -! 2006 ECH: Streamlined for efficiency +! 2006 ECH: Streamlined for efficiency ! 2006 ECH: Converted to free source form (F90) ! 2007 BPB: Modified Delta-Eddington shortwave interface ! 2008 ECH: moved ESMF code to its own driver @@ -44,7 +44,7 @@ module CICE_RunMod subroutine CICE_Run - use ice_calendar, only: istep, istep1, time, dt, calendar + use ice_calendar, only: istep, istep1, dt, calendar, advance_timestep use ice_forcing, only: get_forcing_atmo, get_forcing_ocn, & get_wave_spec use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & @@ -77,19 +77,15 @@ subroutine CICE_Run ! timestep loop !-------------------------------------------------------------------- - istep = istep + 1 ! update time step counters - istep1 = istep1 + 1 - time = time + dt ! determine the time and date - call ice_timer_start(timer_couple) ! atm/ocn coupling + call advance_timestep() ! advance timestep and update calendar data + if (z_tracers) call get_atm_bgc ! biogeochemistry call init_flux_atm ! Initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler - call calendar(time) ! at the end of the timestep - call ice_timer_stop(timer_couple) ! atm/ocn coupling call ice_step @@ -98,7 +94,7 @@ subroutine CICE_Run ! end of timestep loop !-------------------------------------------------------------------- - call ice_timer_stop(timer_step) ! end timestepping loop timer + call ice_timer_stop(timer_step) ! end timestepping loop timer end subroutine CICE_Run @@ -140,7 +136,7 @@ subroutine ice_step use ice_prescribed_mod integer (kind=int_kind) :: & - iblk , & ! block index + iblk , & ! block index k , & ! dynamics supercycling index ktherm ! thermodynamics is off when ktherm = -1 @@ -317,7 +313,7 @@ subroutine ice_step if (tr_iso) call write_restart_iso if (tr_aero) call write_restart_aero if (solve_zsal .or. skl_bgc .or. z_tracers) & - call write_restart_bgc + call write_restart_bgc if (tr_brine) call write_restart_hbrine if (kdyn == 2) call write_restart_eap call final_restart @@ -361,12 +357,12 @@ subroutine coupling_prep (iblk) use ice_step_mod, only: ocean_mixed_layer use ice_timers, only: timer_couple, ice_timer_start, ice_timer_stop - integer (kind=int_kind), intent(in) :: & - iblk ! block index + integer (kind=int_kind), intent(in) :: & + iblk ! block index ! local variables - integer (kind=int_kind) :: & + integer (kind=int_kind) :: & ilo,ihi,jlo,jhi, & ! beginning and end of physical domain n , & ! thickness category index i,j , & ! horizontal indices @@ -506,8 +502,8 @@ subroutine coupling_prep (iblk) fsalt_ai (i,j,iblk) = fsalt (i,j,iblk) fhocn_ai (i,j,iblk) = fhocn (i,j,iblk) fswthru_ai(i,j,iblk) = fswthru(i,j,iblk) - fzsal_ai (i,j,iblk) = fzsal (i,j,iblk) - fzsal_g_ai(i,j,iblk) = fzsal_g(i,j,iblk) + fzsal_ai (i,j,iblk) = fzsal (i,j,iblk) + fzsal_g_ai(i,j,iblk) = fzsal_g(i,j,iblk) if (nbtrcr > 0) then do k = 1, nbtrcr @@ -528,7 +524,7 @@ subroutine coupling_prep (iblk) enddo !----------------------------------------------------------------- - ! Divide fluxes by ice area + ! Divide fluxes by ice area ! - the CESM coupler assumes fluxes are per unit ice area ! - also needed for global budget in diagnostics !----------------------------------------------------------------- @@ -580,16 +576,16 @@ subroutine coupling_prep (iblk) if (.not. calc_Tsfc) then !--------------------------------------------------------------- - ! If surface fluxes were provided, conserve these fluxes at ice - ! free points by passing to ocean. + ! If surface fluxes were provided, conserve these fluxes at ice + ! free points by passing to ocean. !--------------------------------------------------------------- - call sfcflux_to_ocn & + call sfcflux_to_ocn & (nx_block, ny_block, & tmask (:,:,iblk), aice_init(:,:,iblk), & fsurfn_f (:,:,:,iblk), flatn_f(:,:,:,iblk), & fresh (:,:,iblk), fhocn (:,:,iblk)) - endif + endif !echmod call ice_timer_stop(timer_couple,iblk) ! atm/ocn coupling @@ -598,10 +594,10 @@ end subroutine coupling_prep !======================================================================= ! ! If surface heat fluxes are provided to CICE instead of CICE calculating -! them internally (i.e. .not. calc_Tsfc), then these heat fluxes can +! them internally (i.e. .not. calc_Tsfc), then these heat fluxes can ! be provided at points which do not have ice. (This is could be due to ! the heat fluxes being calculated on a lower resolution grid or the -! heat fluxes not recalculated at every CICE timestep.) At ice free points, +! heat fluxes not recalculated at every CICE timestep.) At ice free points, ! conserve energy and water by passing these fluxes to the ocean. ! ! author: A. McLaren, Met Office diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 8d2e23740..ec409495b 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -15,30 +15,21 @@ module ice_comp_nuopc use NUOPC_Model , only : model_label_SetRunClock => label_SetRunClock use NUOPC_Model , only : model_label_Finalize => label_Finalize use NUOPC_Model , only : NUOPC_ModelGet, SetVM - use ice_constants , only : ice_init_constants + use ice_constants , only : ice_init_constants, c0 use ice_shr_methods , only : chkerr, state_setscalar, state_getscalar, state_diagnose, alarmInit - use ice_shr_methods , only : set_component_logging, get_component_instance - use ice_shr_methods , only : state_flddebug - use ice_import_export , only : ice_import, ice_export - use ice_import_export , only : ice_advertise_fields, ice_realize_fields + use ice_shr_methods , only : set_component_logging, get_component_instance, state_flddebug + use ice_import_export , only : ice_import, ice_export, ice_advertise_fields, ice_realize_fields use ice_domain_size , only : nx_global, ny_global - use ice_domain , only : nblocks, blocks_ice, distrb_info - use ice_blocks , only : block, get_block, nx_block, ny_block, nblocks_x, nblocks_y - use ice_blocks , only : nblocks_tot, get_block_parameter - use ice_distribution , only : ice_distributiongetblockloc - use ice_grid , only : tlon, tlat, hm, tarea, ULON, ULAT + use ice_grid , only : grid_type, init_grid2 use ice_communicate , only : init_communicate, my_task, master_task, mpi_comm_ice use ice_calendar , only : force_restart_now, write_ic - use ice_calendar , only : idate, mday, time, mmonth, time2sec, year_init + use ice_calendar , only : idate, mday, mmonth, year_init, timesecs use ice_calendar , only : msec, dt, calendar, calendar_type, nextsw_cday, istep use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long - use ice_scam , only : scmlat, scmlon, single_column use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit - use ice_restart_shared , only : runid, runtype, restart_dir, restart_file + use ice_restart_shared , only : runid, runtype, restart, use_restart_time, restart_dir, restart_file use ice_history , only : accum_hist - use CICE_InitMod , only : cice_init - use CICE_RunMod , only : cice_run use ice_exit , only : abort_ice use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc , only : icepack_init_orbit, icepack_init_parameters, icepack_query_orbit @@ -48,9 +39,15 @@ module ice_comp_nuopc #ifdef CESMCOUPLED use shr_const_mod use shr_orb_mod , only : shr_orb_decl, shr_orb_params, SHR_ORB_UNDEF_REAL, SHR_ORB_UNDEF_INT + use ice_scam , only : scmlat, scmlon, scol_mask, scol_frac, scol_ni, scol_nj #endif use ice_timers + use CICE_InitMod , only : cice_init1, cice_init2 + use CICE_RunMod , only : cice_run + use ice_mesh_mod , only : ice_mesh_set_distgrid, ice_mesh_setmask_from_maskfile, ice_mesh_check + use ice_mesh_mod , only : ice_mesh_init_tlon_tlat_area_hm, ice_mesh_create_scolumn use ice_prescribed_mod , only : ice_prescribed_init + use ice_scam , only : scol_valid, single_column implicit none private @@ -86,6 +83,10 @@ module ice_comp_nuopc character(len=*),parameter :: shr_cal_noleap = 'NO_LEAP' character(len=*),parameter :: shr_cal_gregorian = 'GREGORIAN' + type(ESMF_Mesh) :: ice_mesh + + integer :: nthrds ! Number of threads to use in this component + integer :: dbug = 0 integer , parameter :: debug_import = 0 ! internal debug level integer , parameter :: debug_export = 0 ! internal debug level @@ -179,8 +180,51 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Local variables character(len=char_len_long) :: cvalue - character(len=char_len_long) :: logmsg + character(len=char_len_long) :: ice_meshfile + character(len=char_len_long) :: ice_maskfile + character(len=char_len_long) :: errmsg logical :: isPresent, isSet + real(dbl_kind) :: eccen, obliqr, lambm0, mvelpp + type(ESMF_DistGrid) :: ice_distGrid + real(kind=dbl_kind) :: atmiter_conv + real(kind=dbl_kind) :: atmiter_conv_driver + integer (kind=int_kind) :: natmiter + integer (kind=int_kind) :: natmiter_driver + character(len=char_len) :: tfrz_option_driver ! tfrz_option from driver attributes + character(len=char_len) :: tfrz_option ! tfrz_option from cice namelist + integer(int_kind) :: ktherm + integer :: localPet + integer :: npes + logical :: mastertask + type(ESMF_VM) :: vm + integer :: lmpicom ! local communicator + type(ESMF_Time) :: currTime ! Current time + type(ESMF_Time) :: startTime ! Start time + type(ESMF_Time) :: stopTime ! Stop time + type(ESMF_Time) :: refTime ! Ref time + type(ESMF_TimeInterval) :: timeStep ! Model timestep + type(ESMF_Calendar) :: esmf_calendar ! esmf calendar + type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type + integer :: start_ymd ! Start date (YYYYMMDD) + integer :: start_tod ! start time of day (s) + integer :: curr_ymd ! Current date (YYYYMMDD) + integer :: curr_tod ! Current time of day (s) + integer :: stop_ymd ! stop date (YYYYMMDD) + integer :: stop_tod ! stop time of day (sec) + integer :: ref_ymd ! Reference date (YYYYMMDD) + integer :: ref_tod ! reference time of day (s) + integer :: yy,mm,dd ! Temporaries for time query + integer :: iyear ! yyyy + integer :: dtime ! time step + integer :: shrlogunit ! original log unit + character(len=char_len) :: starttype ! infodata start type + integer :: lsize ! local size of coupling array + integer :: n,c,g,i,j,m ! indices + integer :: iblk, jblk ! indices + integer :: ig, jg ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + character(len=char_len_long) :: diag_filename = 'unset' + character(len=char_len_long) :: logmsg character(len=*), parameter :: subname=trim(modName)//':(InitializeAdvertise) ' !-------------------------------- @@ -244,102 +288,26 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) write(logmsg,'(i6)') dbug call ESMF_LogWrite('CICE_cap: dbug = '//trim(logmsg), ESMF_LOGMSG_INFO) - call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - end subroutine InitializeAdvertise - - !=============================================================================== - - subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) - - ! Arguments - type(ESMF_GridComp) :: gcomp - type(ESMF_State) :: importState - type(ESMF_State) :: exportState - type(ESMF_Clock) :: clock - integer, intent(out) :: rc - - ! Local variables - real(dbl_kind) :: eccen, obliqr, lambm0, mvelpp - type(ESMF_DistGrid) :: distGrid - type(ESMF_Mesh) :: Emesh, EmeshTemp - integer :: spatialDim - integer :: numOwnedElements - real(dbl_kind), pointer :: ownedElemCoords(:) - real(dbl_kind), pointer :: lat(:), latMesh(:) - real(dbl_kind), pointer :: lon(:), lonMesh(:) - integer , allocatable :: gindex_ice(:) - integer , allocatable :: gindex_elim(:) - integer , allocatable :: gindex(:) - integer :: globalID - character(ESMF_MAXSTR) :: cvalue - character(len=char_len) :: tfrz_option - character(ESMF_MAXSTR) :: convCIM, purpComp - type(ESMF_VM) :: vm - type(ESMF_Time) :: currTime ! Current time - type(ESMF_Time) :: startTime ! Start time - type(ESMF_Time) :: stopTime ! Stop time - type(ESMF_Time) :: refTime ! Ref time - type(ESMF_TimeInterval) :: timeStep ! Model timestep - type(ESMF_Calendar) :: esmf_calendar ! esmf calendar - type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type - integer :: start_ymd ! Start date (YYYYMMDD) - integer :: start_tod ! start time of day (s) - integer :: curr_ymd ! Current date (YYYYMMDD) - integer :: curr_tod ! Current time of day (s) - integer :: stop_ymd ! stop date (YYYYMMDD) - integer :: stop_tod ! stop time of day (sec) - integer :: ref_ymd ! Reference date (YYYYMMDD) - integer :: ref_tod ! reference time of day (s) - integer :: yy,mm,dd ! Temporaries for time query - integer :: iyear ! yyyy - integer :: dtime ! time step - integer :: lmpicom - integer :: shrlogunit ! original log unit - character(len=char_len) :: starttype ! infodata start type - integer :: lsize ! local size of coupling array - logical :: isPresent - logical :: isSet - integer :: localPet - integer :: n,c,g,i,j,m ! indices - integer :: iblk, jblk ! indices - integer :: ig, jg ! indices - integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain - type(block) :: this_block ! block information for current block - integer :: compid ! component id - character(len=char_len_long) :: tempc1,tempc2 - real(dbl_kind) :: diff_lon - integer :: npes - integer :: num_elim_global - integer :: num_elim_local - integer :: num_elim - integer :: num_ice - integer :: num_elim_gcells ! local number of eliminated gridcells - integer :: num_elim_blocks ! local number of eliminated blocks - integer :: num_total_blocks - integer :: my_elim_start, my_elim_end - real(dbl_kind) :: rad_to_deg - integer(int_kind) :: ktherm - logical :: mastertask - character(len=char_len_long) :: diag_filename = 'unset' - character(len=*), parameter :: F00 = "('(ice_comp_nuopc) ',2a,1x,d21.14)" - character(len=*), parameter :: subname=trim(modName)//':(InitializeRealize) ' - !-------------------------------- - - rc = ESMF_SUCCESS - if (dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - !---------------------------------------------------------------------------- ! generate local mpi comm !---------------------------------------------------------------------------- call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMGet(vm, mpiCommunicator=lmpicom, localPet=localPet, PetCount=npes, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return +#ifdef CESMCOUPLED + call ESMF_VMGet(vm, pet=localPet, peCount=nthrds, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (nthrds==1) then + call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=u_FILE_u)) return + read(cvalue,*) nthrds + endif +!$ call omp_set_num_threads(nthrds) +#endif + !---------------------------------------------------------------------------- ! Initialize cice communicators !---------------------------------------------------------------------------- @@ -371,6 +339,8 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call ice_init_constants(omega_in=SHR_CONST_OMEGA, radius_in=SHR_CONST_REARTH, & spval_dbl_in=SHR_CONST_SPVAL) + ! TODO: get tfrz_option from driver + call icepack_init_parameters( & secday_in = SHR_CONST_CDAY, & rhoi_in = SHR_CONST_RHOICE, & @@ -395,7 +365,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) Tocnfrz_in = -34.0_dbl_kind*0.054_dbl_kind, & pi_in = SHR_CONST_PI, & snowpatch_in = 0.005_dbl_kind, & - dragio_in = 0.00962_dbl_kind) + dragio_in = 0.00536_dbl_kind) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & @@ -422,8 +392,12 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) runtype = "initial" else if (trim(starttype) == trim('continue') ) then runtype = "continue" + restart = .true. + use_restart_time = .true. else if (trim(starttype) == trim('branch')) then runtype = "continue" + restart = .true. + use_restart_time = .true. else call abort_ice( subname//' ERROR: unknown starttype' ) end if @@ -443,23 +417,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) runtype = 'initial' ! determined from the namelist in ice_init if CESMCOUPLED is not defined end if - ! Determine if single column - call NUOPC_CompAttributeGet(gcomp, name='single_column', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - read(cvalue,*) single_column - if (single_column) then - call NUOPC_CompAttributeGet(gcomp, name='scmlon', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) scmlon - call NUOPC_CompAttributeGet(gcomp, name='scmlat', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) scmlat - end if - else - single_column = .false. - end if - ! Determine runid call NUOPC_CompAttributeGet(gcomp, name='case_name', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) if (isPresent .and. isSet) then @@ -514,12 +471,14 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call shr_file_setLogUnit (shrlogunit) - call NUOPC_CompAttributeGet(gcomp, name="diro", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + call NUOPC_CompAttributeGet(gcomp, name="diro", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then diag_filename = trim(cvalue) end if - call NUOPC_CompAttributeGet(gcomp, name="logfile", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + call NUOPC_CompAttributeGet(gcomp, name="logfile", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then diag_filename = trim(diag_filename) // '/' // trim(cvalue) @@ -532,15 +491,125 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end if !---------------------------------------------------------------------------- - ! Initialize cice + ! First cice initialization phase - before initializing grid info !---------------------------------------------------------------------------- - ! Note that cice_init also sets time manager info as well as mpi communicator info, - ! including master_task and my_task + ! Read the cice namelist as part of the call to cice_init1 + call t_startf ('cice_init1') + call cice_init1 + call t_stopf ('cice_init1') + +#ifdef CESMCOUPLED + ! Form of ocean freezing temperature + ! 'minus1p8' = -1.8 C + ! 'linear_salt' = -depressT * sss + ! 'mushy' conforms with ktherm=2 + call NUOPC_CompAttributeGet(gcomp, name="tfreeze_option", value=tfrz_option_driver, & + isPresent=isPresent, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (.not. isPresent) then + tfrz_option_driver = 'linear_salt' + end if + call icepack_query_parameters( tfrz_option_out=tfrz_option) + if (tfrz_option_driver /= tfrz_option) then + write(errmsg,'(a)') trim(subname)//'error: tfrz_option from driver '//trim(tfrz_option_driver)//& + ' must be the same as tfrz_option from cice namelist '//trim(tfrz_option) + call abort_ice(trim(errmsg)) + endif + + ! Flux convergence tolerance - always use the driver attribute value + call NUOPC_CompAttributeGet(gcomp, name="flux_convergence", value=cvalue, & + isPresent=isPresent, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent) then + read(cvalue,*) atmiter_conv_driver + call icepack_query_parameters( atmiter_conv_out=atmiter_conv) + if (atmiter_conv_driver /= atmiter_conv) then + write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'warning: atmiter_ from driver ',& + atmiter_conv_driver,' is overwritting atmiter_conv from cice namelist ',atmiter_conv + write(nu_diag,*) trim(errmsg) + call icepack_warnings_flush(nu_diag) + call icepack_init_parameters(atmiter_conv_in=atmiter_conv_driver) + end if + end if + + ! Number of iterations for boundary layer calculations + call NUOPC_CompAttributeGet(gcomp, name="flux_max_iteration", value=cvalue, isPresent=isPresent, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent) then + read(cvalue,*) natmiter_driver + else + natmiter_driver = 5 + end if + call icepack_query_parameters( natmiter_out=natmiter) + if (natmiter_driver /= natmiter) then + write(errmsg,'(a,i8,a,i8)') trim(subname)//'error: natmiter_driver ',natmiter_driver, & + ' must be the same as natmiter from cice namelist ',natmiter + call abort_ice(trim(errmsg)) + endif +#endif + !---------------------------------------------------------------------------- + ! Initialize grid info + !---------------------------------------------------------------------------- + + ! Initialize cice mesh and mask if appropriate + + if (single_column .and. scol_valid) then + call ice_mesh_init_tlon_tlat_area_hm() + else + ! Determine mesh input file + call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=ice_meshfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Determine mask input file + call NUOPC_CompAttributeGet(gcomp, name='mesh_mask', value=cvalue, isPresent=isPresent, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + ice_maskfile = trim(cvalue) + else + ice_maskfile = ice_meshfile + end if + if (my_task == master_task) then + write(nu_diag,*)'mesh file for cice domain is ',trim(ice_meshfile) + write(nu_diag,*)'mask file for cice domain is ',trim(ice_maskfile) + end if + + ! Determine the model distgrid using the decomposition obtained in + ! call to init_grid1 called from cice_init1 + call ice_mesh_set_distgrid(localpet, npes, ice_distgrid, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Read in the ice mesh on the cice distribution + ice_mesh = ESMF_MeshCreate(filename=trim(ice_meshfile), fileformat=ESMF_FILEFORMAT_ESMFMESH, & + elementDistGrid=ice_distgrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - call t_startf ('cice_init') - call cice_init - call t_stopf ('cice_init') + ! Initialize the cice mesh and the cice mask + if (trim(grid_type) == 'setmask') then + ! In this case cap code determines the mask file + call ice_mesh_setmask_from_maskfile(ice_maskfile, ice_mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ice_mesh_init_tlon_tlat_area_hm() + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + ! In this case init_grid2 will initialize tlon, tlat, area and hm + call init_grid2() + call ice_mesh_check(gcomp,ice_mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end if + + !---------------------------------------------------------------------------- + ! Second cice initialization phase -after initializing grid info + !---------------------------------------------------------------------------- + + ! Note that cice_init2 also sets time manager info as well as mpi communicator info, + ! including master_task and my_task + ! Note that cice_init2 calls ice_init() which in turn calls icepack_init_parameters + ! which sets the tfrz_option + call t_startf ('cice_init2') + call cice_init2() + call t_stopf ('cice_init2') !---------------------------------------------------------------------------- ! reset shr logging to my log file @@ -554,14 +623,14 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! Now write output to nu_diag - this must happen AFTER call to cice_init if (mastertask) then - write(nu_diag,F00) trim(subname),' cice init nextsw_cday = ',nextsw_cday - write(nu_diag,*) trim(subname),' tfrz_option = ',trim(tfrz_option) + write(nu_diag,'(a,d21.14)') trim(subname)//' cice init nextsw_cday = ',nextsw_cday + write(nu_diag,'(a)') trim(subname)//' tfrz_option = '//trim(tfrz_option) if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then write(nu_diag,*) trim(subname),' Warning: Using ktherm = 2 and tfrz_option = ', trim(tfrz_option) endif - write(nu_diag,*) trim(subname),' inst_name = ',trim(inst_name) - write(nu_diag,*) trim(subname),' inst_index = ',inst_index - write(nu_diag,*) trim(subname),' inst_suffix = ',trim(inst_suffix) + write(nu_diag,'(a )') trim(subname)//' inst_name = '//trim(inst_name) + write(nu_diag,'(a,i8 )') trim(subname)//' inst_index = ',inst_index + write(nu_diag,'(a )') trim(subname)//' inst_suffix = ',trim(inst_suffix) endif !--------------------------------------------------------------------------- @@ -620,247 +689,153 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call time2sec(iyear-year_init,mmonth,mday,time) endif #endif - time = time+start_tod + timesecs = timesecs+start_tod end if - call calendar(time) ! update calendar info + call calendar() ! update calendar info if (write_ic) then call accum_hist(dt) ! write initial conditions end if - !--------------------------------------------------------------------------- - ! Determine the global index space needed for the distgrid - !--------------------------------------------------------------------------- - - ! number the local grid to get allocation size for gindex_ice - lsize = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - lsize = lsize + 1 - enddo - enddo - enddo - - ! set global index array - allocate(gindex_ice(lsize)) - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - gindex_ice(n) = (jg-1)*nx_global + ig - enddo - enddo - enddo - - ! Determine total number of eliminated blocks globally - globalID = 0 - num_elim_global = 0 ! number of eliminated blocks - num_total_blocks = 0 - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - num_total_blocks = num_total_blocks + 1 - if (distrb_info%blockLocation(globalID) == 0) then - num_elim_global = num_elim_global + 1 - end if - end do - end do + !----------------------------------------------------------------- + ! Prescribed ice initialization + !----------------------------------------------------------------- - if (num_elim_global > 0) then + call ice_prescribed_init(clock, ice_mesh, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Distribute the eliminated blocks in a round robin fashion amoung processors - num_elim_local = num_elim_global / npes - my_elim_start = num_elim_local*localPet + min(localPet, mod(num_elim_global, npes)) + 1 - if (localPet < mod(num_elim_global, npes)) then - num_elim_local = num_elim_local + 1 - end if - my_elim_end = my_elim_start + num_elim_local - 1 - - ! Determine the number of eliminated gridcells locally - globalID = 0 - num_elim_blocks = 0 ! local number of eliminated blocks - num_elim_gcells = 0 - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - if (distrb_info%blockLocation(globalID) == 0) then - num_elim_blocks = num_elim_blocks + 1 - if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then - this_block = get_block(globalID, globalID) - num_elim_gcells = num_elim_gcells + & - (this_block%jhi-this_block%jlo+1) * (this_block%ihi-this_block%ilo+1) - end if - end if - end do - end do - - ! Determine the global index space of the eliminated gridcells - allocate(gindex_elim(num_elim_gcells)) - globalID = 0 - num_elim_gcells = 0 ! local number of eliminated gridcells - num_elim_blocks = 0 ! local number of eliminated blocks - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - if (distrb_info%blockLocation(globalID) == 0) then - this_block = get_block(globalID, globalID) - num_elim_blocks = num_elim_blocks + 1 - if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - num_elim_gcells = num_elim_gcells + 1 - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - gindex_elim(num_elim_gcells) = (jg-1)*nx_global + ig - end do - end do - end if - end if - end do - end do - - ! create a global index that includes both active and eliminated gridcells - num_ice = size(gindex_ice) - num_elim = size(gindex_elim) - allocate(gindex(num_elim + num_ice)) - do n = 1,num_ice - gindex(n) = gindex_ice(n) - end do - do n = num_ice+1,num_ice+num_elim - gindex(n) = gindex_elim(n-num_ice) - end do - - deallocate(gindex_elim) + !----------------------------------------------------------------- + ! Advertise fields + !----------------------------------------------------------------- - else + ! NOTE: the advertise phase needs to be called after the ice + ! initialization since the number of ice categories is needed for + ! ice_fraction_n and mean_sw_pen_to_ocn_ifrac_n + call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! No eliminated land blocks - num_ice = size(gindex_ice) - allocate(gindex(num_ice)) - do n = 1,num_ice - gindex(n) = gindex_ice(n) - end do + call t_stopf ('cice_init_total') - end if + end subroutine InitializeAdvertise - !--------------------------------------------------------------------------- - ! Create distGrid from global index array - !--------------------------------------------------------------------------- + !=============================================================================== - DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) - !--------------------------------------------------------------------------- - ! Create the CICE mesh - !--------------------------------------------------------------------------- + ! Arguments + type(ESMF_GridComp) :: gcomp + type(ESMF_State) :: importState + type(ESMF_State) :: exportState + type(ESMF_Clock) :: clock + integer, intent(out) :: rc - ! read in the mesh - call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Local variables + integer :: n + integer :: fieldcount + type(ESMF_Field) :: lfield + character(len=char_len_long) :: cvalue + real(dbl_kind) :: scol_lon + real(dbl_kind) :: scol_lat + real(dbl_kind) :: scol_spval + real(dbl_kind), pointer :: fldptr1d(:) + real(dbl_kind), pointer :: fldptr2d(:,:) + integer :: rank + character(len=char_len_long) :: single_column_lnd_domainfile + character(len=char_len_long) , pointer :: lfieldnamelist(:) => null() + character(len=*), parameter :: subname=trim(modName)//':(InitializeRealize) ' + !-------------------------------- - EMeshTemp = ESMF_MeshCreate(filename=trim(cvalue), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (my_task == master_task) then - write(nu_diag,*)'mesh file for cice domain is ',trim(cvalue) - end if + rc = ESMF_SUCCESS + if (dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - ! recreate the mesh using the above distGrid - EMesh = ESMF_MeshCreate(EMeshTemp, elementDistgrid=Distgrid, rc=rc) +#ifdef CESMCOUPLED + call NUOPC_CompAttributeGet(gcomp, name='scol_lon', value=cvalue, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! obtain mesh lats and lons - call ESMF_MeshGet(Emesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) + read(cvalue,*) scmlon + call NUOPC_CompAttributeGet(gcomp, name='scol_lat', value=cvalue, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - allocate(ownedElemCoords(spatialDim*numOwnedElements)) - allocate(lonMesh(numOwnedElements), latMesh(numOwnedElements)) - call ESMF_MeshGet(Emesh, ownedElemCoords=ownedElemCoords) + read(cvalue,*) scmlat + call NUOPC_CompAttributeGet(gcomp, name='scol_spval', value=cvalue, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_spval - do n = 1,numOwnedElements - lonMesh(n) = ownedElemCoords(2*n-1) - latMesh(n) = ownedElemCoords(2*n) - end do + if (scmlon > scol_spval .and. scmlat > scol_spval) then + call NUOPC_CompAttributeGet(gcomp, name='single_column_lnd_domainfile', & + value=single_column_lnd_domainfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (trim(single_column_lnd_domainfile) /= 'UNSET') then + single_column = .true. + else + call abort_ice('single_column_domainfile cannot be null for single column mode') + end if + call NUOPC_CompAttributeGet(gcomp, name='scol_ocnmask', value=cvalue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_mask + call NUOPC_CompAttributeGet(gcomp, name='scol_ocnfrac', value=cvalue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_frac + call NUOPC_CompAttributeGet(gcomp, name='scol_ni', value=cvalue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_ni + call NUOPC_CompAttributeGet(gcomp, name='scol_nj', value=cvalue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_nj - call icepack_query_parameters(rad_to_deg_out=rad_to_deg) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) + call ice_mesh_create_scolumn(scmlon, scmlat, ice_mesh, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! obtain internally generated cice lats and lons for error checks - allocate(lon(lsize)) - allocate(lat(lsize)) - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - lon(n) = tlon(i,j,iblk)*rad_to_deg - lat(n) = tlat(i,j,iblk)*rad_to_deg + scol_valid = (scol_mask == 1) + if (.not. scol_valid) then + ! if single column is not valid - set all export state fields to zero and return + write(nu_diag,'(a)')' (ice_comp_nuopc) single column mode point does not contain any ocn/ice '& + //' - setting all export data to 0' + call ice_realize_fields(gcomp, mesh=ice_mesh, & + flds_scalar_name=flds_scalar_name, flds_scalar_num=flds_scalar_num, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_StateGet(exportState, itemCount=fieldCount, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(lfieldnamelist(fieldCount)) + call ESMF_StateGet(exportState, itemNameList=lfieldnamelist, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + do n = 1, fieldCount + if (trim(lfieldnamelist(n)) /= flds_scalar_name) then + call ESMF_StateGet(exportState, itemName=trim(lfieldnamelist(n)), field=lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, rank=rank, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (rank == 2) then + call ESMF_FieldGet(lfield, farrayPtr=fldptr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + fldptr2d(:,:) = 0._dbl_kind + else + call ESMF_FieldGet(lfield, farrayPtr=fldptr1d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + fldptr1d(:) = 0._dbl_kind + end if + end if enddo - enddo - enddo - - ! error check differences between internally generated lons and those read in - do n = 1,lsize - diff_lon = abs(lonMesh(n) - lon(n)) - if ( (diff_lon > 1.e2 .and. abs(diff_lon - 360_dbl_kind) > 1.e-1) .or.& - (diff_lon > 1.e-3 .and. diff_lon < 1._dbl_kind) ) then - !write(6,100)n,lonMesh(n),lon(n), diff_lon -100 format('ERROR: CICE n, lonmesh(n), lon(n), diff_lon = ',i6,2(f21.13,3x),d21.5) - !call abort_ice() - end if - if (abs(latMesh(n) - lat(n)) > 1.e-1) then - !write(6,101)n,latMesh(n),lat(n), abs(latMesh(n)-lat(n)) -101 format('ERROR: CICE n, latmesh(n), lat(n), diff_lat = ',i6,2(f21.13,3x),d21.5) - !call abort_ice() + deallocate(lfieldnamelist) + ! ******************* + ! *** RETURN HERE *** + ! ******************* + RETURN + else + write(nu_diag,'(a,3(f10.5,2x))')' (ice_comp_nuopc) single column mode lon/lat/frac is ',& + scmlon,scmlat,scol_frac end if - end do - - ! deallocate memory - deallocate(ownedElemCoords) - deallocate(lon, lonMesh) - deallocate(lat, latMesh) + else + single_column = .false. + end if +#endif !----------------------------------------------------------------- ! Realize the actively coupled fields !----------------------------------------------------------------- - call ice_realize_fields(gcomp, mesh=Emesh, & + call ice_realize_fields(gcomp, mesh=ice_mesh, & flds_scalar_name=flds_scalar_name, flds_scalar_num=flds_scalar_num, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - !----------------------------------------------------------------- - ! Prescribed ice initialization - first get compid - !----------------------------------------------------------------- - - call NUOPC_CompAttributeGet(gcomp, name='MCTID', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - read(cvalue,*) compid ! convert from string to integer - else - compid = 0 - end if - call ice_prescribed_init(lmpicom, compid, gindex_ice) - !----------------------------------------------------------------- ! Create cice export state !----------------------------------------------------------------- @@ -875,16 +850,16 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) flds_scalar_name, flds_scalar_num, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !-------------------------------- + ! diagnostics + !-------------------------------- + ! TODO (mvertens, 2018-12-21): fill in iceberg_prognostic as .false. if (debug_export > 0 .and. my_task==master_task) then call State_fldDebug(exportState, flds_scalar_name, 'cice_export:', & idate, msec, nu_diag, rc=rc) end if - !-------------------------------- - ! diagnostics - !-------------------------------- - if (dbug > 0) then call state_diagnose(exportState,subname//':ES',rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -892,11 +867,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if (dbug > 5) call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) - call t_stopf ('cice_init_total') - - deallocate(gindex_ice) - deallocate(gindex) - call flush_fileunit(nu_diag) end subroutine InitializeRealize @@ -939,7 +909,6 @@ subroutine ModelAdvance(gcomp, rc) character(char_len_long) :: restart_date character(char_len_long) :: restart_filename logical :: isPresent, isSet - character(*) , parameter :: F00 = "('(ice_comp_nuopc) ',2a,i8,d21.14)" character(len=*),parameter :: subname=trim(modName)//':(ModelAdvance) ' character(char_len_long) :: msgString !-------------------------------- @@ -1005,7 +974,7 @@ subroutine ModelAdvance(gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if if (my_task == master_task) then - write(nu_diag,F00) trim(subname),' cice istep, nextsw_cday = ',istep, nextsw_cday + write(nu_diag,'(a,2x,i8,2x,d24.14)') trim(subname)//' cice istep, nextsw_cday = ',istep, nextsw_cday end if !-------------------------------- @@ -1283,28 +1252,26 @@ end subroutine ModelSetRunClock !=============================================================================== subroutine ModelFinalize(gcomp, rc) + + !-------------------------------- + ! Finalize routine + !-------------------------------- + type(ESMF_GridComp) :: gcomp integer, intent(out) :: rc ! local variables - character(*), parameter :: F00 = "('(ice_comp_nuopc) ',8a)" - character(*), parameter :: F91 = "('(ice_comp_nuopc) ',73('-'))" + character(*), parameter :: F91 = "('(ice_comp_nuopc) ',73('-'))" character(len=*),parameter :: subname=trim(modName)//':(ModelFinalize) ' !-------------------------------- - !-------------------------------- - ! Finalize routine - !-------------------------------- - rc = ESMF_SUCCESS if (dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - if (my_task == master_task) then write(nu_diag,F91) - write(nu_diag,F00) 'CICE: end of main integration loop' + write(nu_diag,'(a)') 'CICE: end of main integration loop' write(nu_diag,F91) end if - if (dbug > 5) call ESMF_LogWrite(subname//' done', ESMF_LOGMSG_INFO) end subroutine ModelFinalize @@ -1469,7 +1436,4 @@ subroutine ice_cal_ymd2date(year, month, day, date) end subroutine ice_cal_ymd2date - !=============================================================================== - - end module ice_comp_nuopc diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index b32085143..7f394dd61 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -4,7 +4,7 @@ module ice_import_export use NUOPC use NUOPC_Model use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind - use ice_constants , only : c0, c1, spval_dbl + use ice_constants , only : c0, c1, spval_dbl, radius use ice_constants , only : field_loc_center, field_type_scalar, field_type_vector use ice_blocks , only : block, get_block, nx_block, ny_block use ice_domain , only : nblocks, blocks_ice, halo_info, distrb_info @@ -21,10 +21,12 @@ module ice_import_export use ice_flux , only : fresh, fsalt, zlvl, uatm, vatm, potT, Tair, Qa use ice_flux , only : rhoa, swvdr, swvdf, swidr, swidf, flw, frain use ice_flux , only : fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt + use ice_flux , only : send_i2x_per_cat use ice_flux , only : sss, Tf, wind, fsw use ice_state , only : vice, vsno, aice, aicen_init, trcr - use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm, ocn_gridcell_frac + use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm use ice_grid , only : grid_type, t2ugrid_vector + use ice_mesh_mod , only : ocn_gridcell_frac use ice_boundary , only : ice_HaloUpdate use ice_fileunits , only : nu_diag, flush_fileunit use ice_communicate , only : my_task, master_task, MPI_COMM_ICE @@ -34,9 +36,10 @@ module ice_import_export use icepack_intfc , only : icepack_query_parameters, icepack_query_tracer_flags use icepack_intfc , only : icepack_liquidus_temperature use icepack_intfc , only : icepack_sea_freezing_temperature - use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf + use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED use shr_frz_mod , only : shr_frz_freezetemp + use shr_mpi_mod , only : shr_mpi_min, shr_mpi_max #endif implicit none @@ -54,20 +57,18 @@ module ice_import_export interface state_getfldptr module procedure state_getfldptr_1d module procedure state_getfldptr_2d - module procedure state_getfldptr_3d - module procedure state_getfldptr_4d end interface state_getfldptr private :: state_getfldptr interface state_getimport - module procedure state_getimport_4d_output - module procedure state_getimport_3d_output + module procedure state_getimport_4d + module procedure state_getimport_3d end interface state_getimport private :: state_getimport interface state_setexport - module procedure state_setexport_4d_input - module procedure state_setexport_3d_input + module procedure state_setexport_4d + module procedure state_setexport_3d end interface state_setexport private :: state_setexport @@ -79,12 +80,15 @@ module ice_import_export integer :: ungridded_ubound = 0 end type fld_list_type + ! area correction factors for fluxes send and received from mediator + real(dbl_kind), allocatable :: mod2med_areacor(:) ! ratios of model areas to input mesh areas + real(dbl_kind), allocatable :: med2mod_areacor(:) ! ratios of input mesh areas to model areas + integer, parameter :: fldsMax = 100 integer :: fldsToIce_num = 0 integer :: fldsFrIce_num = 0 type (fld_list_type) :: fldsToIce(fldsMax) type (fld_list_type) :: fldsFrIce(fldsMax) - type(ESMF_GeomType_Flag) :: geomtype integer , parameter :: io_dbug = 10 ! i/o debug messages character(*), parameter :: u_FILE_u = & @@ -108,7 +112,6 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam character(char_len) :: stdname character(char_len) :: cvalue logical :: flds_wiso ! use case - logical :: flds_i2o_per_cat ! .true. => select per ice thickness category logical :: isPresent, isSet character(len=*), parameter :: subname='(ice_import_export:ice_advertise_fields)' !------------------------------------------------------------------------------- @@ -116,21 +119,28 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam rc = ESMF_SUCCESS if (io_dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - ! Determine if the following attributes are sent by the driver and if so read them in - flds_wiso = .false. - call NUOPC_CompAttributeGet(gcomp, name='flds_wiso', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + ! Determine if ice sends multiple ice category info back to mediator + send_i2x_per_cat = .false. + call NUOPC_CompAttributeGet(gcomp, name='flds_i2o_per_cat', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then - read(cvalue,*) flds_wiso - call ESMF_LogWrite('flds_wiso = '// trim(cvalue), ESMF_LOGMSG_INFO) + read(cvalue,*) send_i2x_per_cat + end if + if (my_task == master_task) then + write(nu_diag,*)'send_i2x_per_cat = ',send_i2x_per_cat end if - flds_i2o_per_cat = .false. - call NUOPC_CompAttributeGet(gcomp, name='flds_i2o_per_cat', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + ! Determine if the following attributes are sent by the driver and if so read them in + flds_wiso = .false. + call NUOPC_CompAttributeGet(gcomp, name='flds_wiso', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then - read(cvalue,*) send_i2x_per_cat - call ESMF_LogWrite('flds_i2o_per_cat = '// trim(cvalue), ESMF_LOGMSG_INFO) + read(cvalue,*) flds_wiso + end if + if (my_task == master_task) then + write(nu_diag,*)'flds_wiso = ',flds_wiso end if !----------------- @@ -262,21 +272,35 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam end subroutine ice_advertise_fields -!============================================================================== - - subroutine ice_realize_fields(gcomp, mesh, grid, flds_scalar_name, flds_scalar_num, rc) + !============================================================================== + subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc) ! input/output variables - type(ESMF_GridComp) :: gcomp - type(ESMF_Mesh) , optional , intent(in) :: mesh - type(ESMF_Grid) , optional , intent(in) :: grid - character(len=*) , intent(in) :: flds_scalar_name - integer , intent(in) :: flds_scalar_num - integer , intent(out) :: rc + type(ESMF_GridComp) :: gcomp + type(ESMF_Mesh) , intent(in) :: mesh + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(in) :: flds_scalar_num + integer , intent(out) :: rc ! local variables - type(ESMF_State) :: importState - type(ESMF_State) :: exportState + type(ESMF_State) :: importState + type(ESMF_State) :: exportState + type(ESMF_Field) :: lfield + integer :: numOwnedElements + integer :: i, j, iblk, n + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type(block) :: this_block ! block information for current block + real(dbl_kind), allocatable :: mesh_areas(:) + real(dbl_kind), allocatable :: model_areas(:) + real(dbl_kind), pointer :: dataptr(:) + real(dbl_kind) :: max_mod2med_areacor + real(dbl_kind) :: max_med2mod_areacor + real(dbl_kind) :: min_mod2med_areacor + real(dbl_kind) :: min_med2mod_areacor + real(dbl_kind) :: max_mod2med_areacor_glob + real(dbl_kind) :: max_med2mod_areacor_glob + real(dbl_kind) :: min_mod2med_areacor_glob + real(dbl_kind) :: min_med2mod_areacor_glob character(len=*), parameter :: subname='(ice_import_export:realize_fields)' !--------------------------------------------------------------------------- @@ -285,60 +309,86 @@ subroutine ice_realize_fields(gcomp, mesh, grid, flds_scalar_name, flds_scalar_n call NUOPC_ModelGet(gcomp, importState=importState, exportState=exportState, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (present(mesh)) then - - geomtype = ESMF_GEOMTYPE_MESH - - call fldlist_realize( & - state=ExportState, & - fldList=fldsFrIce, & - numflds=fldsFrIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Export',& - mesh=mesh, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call fldlist_realize( & - state=importState, & - fldList=fldsToIce, & - numflds=fldsToIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Import',& - mesh=mesh, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - else if (present(grid)) then - - geomtype = ESMF_GEOMTYPE_GRID + call fldlist_realize( & + state=ExportState, & + fldList=fldsFrIce, & + numflds=fldsFrIce_num, & + flds_scalar_name=flds_scalar_name, & + flds_scalar_num=flds_scalar_num, & + tag=subname//':CICE_Export',& + mesh=mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - call fldlist_realize( & - state=ExportState, & - fldList=fldsFrIce, & - numflds=fldsFrIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Export',& - grid=grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + call fldlist_realize( & + state=importState, & + fldList=fldsToIce, & + numflds=fldsToIce_num, & + flds_scalar_name=flds_scalar_name, & + flds_scalar_num=flds_scalar_num, & + tag=subname//':CICE_Import',& + mesh=mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - call fldlist_realize( & - state=importState, & - fldList=fldsToIce, & - numflds=fldsToIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Import',& - grid=grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return +#ifdef CESMCOUPLED + ! Get mesh areas from second field - using second field since the + ! first field is the scalar field + call ESMF_MeshGet(mesh, numOwnedElements=numOwnedElements, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_StateGet(exportState, itemName=trim(fldsFrIce(2)%stdname), field=lfield, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldRegridGetArea(lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, farrayPtr=dataptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(mesh_areas(numOwnedElements)) + mesh_areas(:) = dataptr(:) + + ! Determine flux correction factors (module variables) + allocate(model_areas(numOwnedElements)) + allocate(mod2med_areacor(numOwnedElements)) + allocate(med2mod_areacor(numOwnedElements)) + mod2med_areacor(:) = 1._dbl_kind + med2mod_areacor(:) = 1._dbl_kind + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + model_areas(n) = tarea(i,j,iblk)/(radius*radius) + mod2med_areacor(n) = model_areas(n) / mesh_areas(n) + med2mod_areacor(n) = mesh_areas(n) / model_areas(n) + enddo + enddo + enddo + deallocate(model_areas) + deallocate(mesh_areas) + + min_mod2med_areacor = minval(mod2med_areacor) + max_mod2med_areacor = maxval(mod2med_areacor) + min_med2mod_areacor = minval(med2mod_areacor) + max_med2mod_areacor = maxval(med2mod_areacor) + call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpi_comm_ice) + call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpi_comm_ice) + call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpi_comm_ice) + call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpi_comm_ice) + + if (my_task == master_task) then + write(nu_diag,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& + min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'CICE6' + write(nu_diag,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& + min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'CICE6' end if +#endif end subroutine ice_realize_fields !============================================================================== - subroutine ice_import( importState, rc ) ! input/output variables @@ -355,7 +405,11 @@ subroutine ice_import( importState, rc ) real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh - real (kind=dbl_kind) :: inst_pres_height_lowest + real (kind=dbl_kind) :: inst_pres_height_lowest + real (kind=dbl_kind), pointer :: dataptr2d(:,:) + real (kind=dbl_kind), pointer :: dataptr1d(:) + real (kind=dbl_kind), pointer :: dataptr2d_dstwet(:,:) + real (kind=dbl_kind), pointer :: dataptr2d_dstdry(:,:) character(len=char_len) :: tfrz_option integer(int_kind) :: ktherm character(len=*), parameter :: subname = 'ice_import' @@ -365,17 +419,20 @@ subroutine ice_import( importState, rc ) call icepack_query_parameters(Tffresh_out=Tffresh) call icepack_query_parameters(tfrz_option_out=tfrz_option) call icepack_query_parameters(ktherm_out=ktherm) - if (io_dbug > 5) then - write(msgString,'(A,i8)')trim(subname)//' tfrz_option = ' & - // trim(tfrz_option)//', ktherm = ',ktherm - call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) + + if (io_dbug > 5) then + write(msgString,'(A,i8)')trim(subname)//' tfrz_option = ' & + // trim(tfrz_option)//', ktherm = ',ktherm + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) end if -! call icepack_query_parameters(tfrz_option_out=tfrz_option, & -! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & -! Tffresh_out=Tffresh) -! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & -! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & -! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + + ! call icepack_query_parameters(tfrz_option_out=tfrz_option, & + ! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & + ! Tffresh_out=Tffresh) + ! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & + ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & + ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -429,30 +486,38 @@ subroutine ice_import( importState, rc ) ! import ocn/ice fluxes - call state_getimport(importState, 'freezing_melting_potential', output=aflds, index=9, rc=rc) + call state_getimport(importState, 'freezing_melting_potential', output=aflds, index=9, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! import atm fluxes - call state_getimport(importState, 'mean_down_sw_vis_dir_flx', output=aflds, index=10, rc=rc) + call state_getimport(importState, 'mean_down_sw_vis_dir_flx', output=aflds, index=10, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_ir_dir_flx', output=aflds, index=11, rc=rc) + call state_getimport(importState, 'mean_down_sw_ir_dir_flx', output=aflds, index=11, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_vis_dif_flx', output=aflds, index=12, rc=rc) + call state_getimport(importState, 'mean_down_sw_vis_dif_flx', output=aflds, index=12, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_ir_dif_flx', output=aflds, index=13, rc=rc) + call state_getimport(importState, 'mean_down_sw_ir_dif_flx', output=aflds, index=13, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_lw_flx', output=aflds, index=14, rc=rc) + call state_getimport(importState, 'mean_down_lw_flx', output=aflds, index=14, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_prec_rate', output=aflds, index=15, rc=rc) + call state_getimport(importState, 'mean_prec_rate', output=aflds, index=15, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_fprec_rate', output=aflds, index=16, rc=rc) + call state_getimport(importState, 'mean_fprec_rate', output=aflds, index=16, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! perform a halo update @@ -488,7 +553,7 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO - if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then + if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1,ny_block @@ -518,7 +583,7 @@ subroutine ice_import( importState, rc ) endif end do !i end do !j - end do !iblk + end do !iblk !$OMP END PARALLEL DO end if @@ -577,34 +642,45 @@ subroutine ice_import( importState, rc ) ! bcphodry ungridded_index=2 ! bcphiwet ungridded_index=3 - ! bcphodry - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=1, ungridded_index=2, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! bcphidry + bcphiwet - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=2, ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=2, do_sum=.true., ungridded_index=3, rc=rc) + call state_getfldptr(importState, 'Faxa_bcph', fldptr=dataPtr2d, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + faero_atm(i,j,1,iblk) = dataPtr2d(2,n) * med2mod_areacor(n) ! bcphodry + faero_atm(i,j,2,iblk) = (dataptr2d(1,n) + dataPtr2d(3,n)) * med2mod_areacor(n) ! bcphidry + bcphiwet + end do + end do + end do end if ! Sum over all dry and wet dust fluxes from ath atmosphere if (State_FldChk(importState, 'Faxa_dstwet') .and. State_FldChk(importState, 'Faxa_dstdry')) then - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=2, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=2, rc=rc) + call state_getfldptr(importState, 'Faxa_dstwet', fldptr=dataPtr2d_dstwet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=3, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=3, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=4, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=4, rc=rc) + call state_getfldptr(importState, 'Faxa_dstdry', fldptr=dataPtr2d_dstdry, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + faero_atm(i,j,3,iblk) = dataPtr2d_dstwet(1,n) + dataptr2d_dstdry(1,n) + & + dataPtr2d_dstwet(2,n) + dataptr2d_dstdry(2,n) + & + dataPtr2d_dstwet(3,n) + dataptr2d_dstdry(3,n) + & + dataPtr2d_dstwet(4,n) + dataptr2d_dstdry(4,n) + faero_atm(i,j,3,iblk) = faero_atm(i,j,3,iblk) * med2mod_areacor(n) + end do + end do + end do end if !------------------------------------------------------- @@ -623,12 +699,15 @@ subroutine ice_import( importState, rc ) call state_getimport(importState, 'inst_spec_humid_height_lowest_wiso', output=Qa_iso, index=3, ungridded_index=2, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=1, ungridded_index=3, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=2, ungridded_index=1, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=3, ungridded_index=2, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=1, ungridded_index=3, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=2, ungridded_index=1, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=3, ungridded_index=2, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getimport(importState, 'mean_fprec_rate_wiso', output=fiso_atm, index=1, ungridded_index=3, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -637,11 +716,14 @@ subroutine ice_import( importState, rc ) call state_getimport(importState, 'mean_fprec_rate_wiso', output=fiso_atm, index=3, ungridded_index=2, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=HDO_ocn , ungridded_index=3, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=HDO_ocn , ungridded_index=3, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=H2_16O_ocn, ungridded_index=1, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=H2_16O_ocn, ungridded_index=1, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=H2_18O_ocn, ungridded_index=2, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=H2_18O_ocn, ungridded_index=2, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -690,9 +772,11 @@ subroutine ice_import( importState, rc ) #ifdef CESMCOUPLED ! Use shr_frz_mod for this - Tf(:,:,iblk) = shr_frz_freezetemp(sss(:,:,iblk)) -#else - !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky) + do iblk = 1, nblocks + Tf(:,:,iblk) = shr_frz_freezetemp(sss(:,:,iblk)) + end do +#else + !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1,ny_block do i = 1,nx_block @@ -747,7 +831,6 @@ subroutine ice_import( importState, rc ) end subroutine ice_import !=============================================================================== - subroutine ice_export( exportState, rc ) ! input/output variables @@ -770,8 +853,10 @@ subroutine ice_export( exportState, rc ) real (kind=dbl_kind) :: tauxo (nx_block,ny_block,max_blocks) ! ice/ocean stress real (kind=dbl_kind) :: tauyo (nx_block,ny_block,max_blocks) ! ice/ocean stress real (kind=dbl_kind) :: ailohi(nx_block,ny_block,max_blocks) ! fractional ice area - real (kind=dbl_kind), allocatable :: tempfld(:,:,:) real (kind=dbl_kind) :: Tffresh + real (kind=dbl_kind), allocatable :: tempfld(:,:,:) + real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:) + real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:) character(len=*),parameter :: subname = 'ice_export' !----------------------------------------------------- @@ -779,12 +864,13 @@ subroutine ice_export( exportState, rc ) if (io_dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) call icepack_query_parameters(Tffresh_out=Tffresh) -! call icepack_query_parameters(tfrz_option_out=tfrz_option, & -! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & -! Tffresh_out=Tffresh) -! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & -! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & -! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + ! call icepack_query_parameters(tfrz_option_out=tfrz_option, & + ! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & + ! Tffresh_out=Tffresh) + ! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & + ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & + ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -880,7 +966,7 @@ subroutine ice_export( exportState, rc ) call state_setexport(exportState, 'ice_fraction', input=ailohi, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (trim(grid_type) == 'latlon') then + if (trim(grid_type) == 'setmask') then call state_setexport(exportState, 'ice_mask', input=ocn_gridcell_frac, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return else @@ -967,31 +1053,38 @@ subroutine ice_export( exportState, rc ) ! ------ ! Zonal air/ice stress - call state_setexport(exportState, 'stress_on_air_ice_zonal' , input=tauxa, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_air_ice_zonal' , input=tauxa, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Meridional air/ice stress - call state_setexport(exportState, 'stress_on_air_ice_merid' , input=tauya, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_air_ice_merid' , input=tauya, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Latent heat flux (atm into ice) - call state_setexport(exportState, 'mean_laten_heat_flx_atm_into_ice' , input=flat, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_laten_heat_flx_atm_into_ice' , input=flat, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Sensible heat flux (atm into ice) - call state_setexport(exportState, 'mean_sensi_heat_flx_atm_into_ice' , input=fsens, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sensi_heat_flx_atm_into_ice' , input=fsens, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! longwave outgoing (upward), average over ice fraction only - call state_setexport(exportState, 'mean_up_lw_flx_ice' , input=flwout, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_up_lw_flx_ice' , input=flwout, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Evaporative water flux (kg/m^2/s) - call state_setexport(exportState, 'mean_evap_rate_atm_into_ice' , input=evap, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_evap_rate_atm_into_ice' , input=evap, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Shortwave flux absorbed in ice and ocean (W/m^2) - call state_setexport(exportState, 'Faii_swnet' , input=fswabs, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Faii_swnet' , input=fswabs, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------ @@ -999,43 +1092,53 @@ subroutine ice_export( exportState, rc ) ! ------ ! flux of shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn' , input=fswthru, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn' , input=fswthru, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of vis dir shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dir_flx' , input=fswthru_vdr, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dir_flx' , input=fswthru_vdr, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of vis dif shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dif_flx' , input=fswthru_vdf, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dif_flx' , input=fswthru_vdf, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of ir dir shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dir_flx' , input=fswthru_idr, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dir_flx' , input=fswthru_idr, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of ir dif shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dif_flx' , input=fswthru_idf, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dif_flx' , input=fswthru_idf, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! heat exchange with ocean - call state_setexport(exportState, 'net_heat_flx_to_ocn' , input=fhocn, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux of heat exchange with ocean + call state_setexport(exportState, 'net_heat_flx_to_ocn' , input=fhocn, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! fresh water to ocean (h2o flux from melting) - call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate' , input=fresh, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux fresh water to ocean (h2o flux from melting) + call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate' , input=fresh, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! salt to ocean (salt flux from melting) - call state_setexport(exportState, 'mean_salt_rate' , input=fsalt, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux of salt to ocean (salt flux from melting) + call state_setexport(exportState, 'mean_salt_rate' , input=fsalt, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! stress n i/o zonal - call state_setexport(exportState, 'stress_on_ocn_ice_zonal' , input=tauxo, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_ocn_ice_zonal' , input=tauxo, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! stress n i/o meridional - call state_setexport(exportState, 'stress_on_ocn_ice_merid' , input=tauyo, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_ocn_ice_merid' , input=tauyo, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------ @@ -1044,19 +1147,22 @@ subroutine ice_export( exportState, rc ) ! hydrophobic bc if (State_FldChk(exportState, 'Fioi_bcpho')) then - call state_setexport(exportState, 'Fioi_bcpho' , input=faero_ocn, index=1, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_bcpho' , input=faero_ocn, index=1, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if ! hydrophilic bc if (State_FldChk(exportState, 'Fioi_bcphi')) then - call state_setexport(exportState, 'Fioi_bcphi' , input=faero_ocn, index=2, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_bcphi' , input=faero_ocn, index=2, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if ! dust if (State_FldChk(exportState, 'Fioi_flxdst')) then - call state_setexport(exportState, 'Fioi_flxdst' , input=faero_ocn, index=3, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_flxdst' , input=faero_ocn, index=3, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -1070,13 +1176,13 @@ subroutine ice_export( exportState, rc ) ! HDO => ungridded_index=3 call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=1, & - lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=3, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=2, & - lmask=tmask, ifrac=ailohi, ungridded_index=1, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=1, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=3, & - lmask=tmask, ifrac=ailohi, ungridded_index=2, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=2, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -1087,16 +1193,16 @@ subroutine ice_export( exportState, rc ) if (State_FldChk(exportState, 'mean_evap_rate_atm_into_ice_wiso')) then ! Isotope evap to atm call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=1, & - lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=3, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=2, & - lmask=tmask, ifrac=ailohi, ungridded_index=1, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=1, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=3, & - lmask=tmask, ifrac=ailohi, ungridded_index=2, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=2, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Isotope evap to atm + ! qref to atm call state_setexport(exportState, 'Si_qref_wiso' , input=Qref_iso, index=1, & lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1124,7 +1230,7 @@ subroutine ice_export( exportState, rc ) ! Note: no need zero out pass-through fields over land for benefit of x2oacc fields in cpl hist files since ! the export state has been zeroed out at the beginning call state_setexport(exportState, 'mean_sw_pen_to_ocn_ifrac_n', input=fswthrun_ai, index=n, & - lmask=tmask, ifrac=ailohi, ungridded_index=n, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=n, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end do end if @@ -1132,7 +1238,6 @@ subroutine ice_export( exportState, rc ) end subroutine ice_export !=============================================================================== - subroutine fldlist_add(num, fldlist, stdname, ungridded_lbound, ungridded_ubound) ! input/output variables @@ -1162,7 +1267,6 @@ subroutine fldlist_add(num, fldlist, stdname, ungridded_lbound, ungridded_ubound end subroutine fldlist_add !=============================================================================== - subroutine fldlist_realize(state, fldList, numflds, flds_scalar_name, flds_scalar_num, mesh, grid, tag, rc) use NUOPC, only : NUOPC_IsConnected, NUOPC_Realize @@ -1187,6 +1291,7 @@ subroutine fldlist_realize(state, fldList, numflds, flds_scalar_name, flds_scala integer :: n type(ESMF_Field) :: field character(len=80) :: stdname + character(ESMF_MAXSTR) :: msg character(len=*),parameter :: subname='(ice_import_export:fld_list_realize)' ! ---------------------------------------------- @@ -1203,8 +1308,6 @@ subroutine fldlist_realize(state, fldList, numflds, flds_scalar_name, flds_scala if (ChkErr(rc,__LINE__,u_FILE_u)) return else if (present(mesh)) then - call ESMF_LogWrite(trim(subname)//trim(tag)//" Field = "//trim(stdname)//" is connected using mesh", & - ESMF_LOGMSG_INFO) ! Create the field if (fldlist(n)%ungridded_lbound > 0 .and. fldlist(n)%ungridded_ubound > 0) then field = ESMF_FieldCreate(mesh, ESMF_TYPEKIND_R8, name=stdname, meshloc=ESMF_MESHLOC_ELEMENT, & @@ -1212,9 +1315,16 @@ subroutine fldlist_realize(state, fldList, numflds, flds_scalar_name, flds_scala ungriddedUbound=(/fldlist(n)%ungridded_ubound/), & gridToFieldMap=(/2/), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(msg, '(a,i4,2x,i4)') trim(subname)//trim(tag)//" Field = "//trim(stdname)//& + " is connected using mesh with lbound, ubound = ",& + fldlist(n)%ungridded_lbound,fldlist(n)%ungridded_ubound + call ESMF_LogWrite(msg, ESMF_LOGMSG_INFO) else field = ESMF_FieldCreate(mesh, ESMF_TYPEKIND_R8, name=stdname, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(msg, '(a,i4,a,i4)') trim(subname)//trim(tag)//" Field = "//trim(stdname)//& + " is connected using mesh without ungridded dimension" + call ESMF_LogWrite(msg, ESMF_LOGMSG_INFO) end if else if (present(grid)) then call ESMF_LogWrite(trim(subname)//trim(tag)//" Field = "//trim(stdname)//" is connected using grid", & @@ -1287,7 +1397,6 @@ end subroutine SetScalarField end subroutine fldlist_realize !=============================================================================== - logical function State_FldChk(State, fldname) ! ---------------------------------------------- ! Determine if field is in state @@ -1302,27 +1411,25 @@ logical function State_FldChk(State, fldname) ! ---------------------------------------------- call ESMF_StateGet(State, trim(fldname), itemType) - State_FldChk = (itemType /= ESMF_STATEITEM_NOTFOUND) end function State_FldChk !=============================================================================== - - subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungridded_index, rc) + subroutine state_getimport_4d(state, fldname, output, index, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array ! ---------------------------------------------- ! input/output variables - type(ESMF_State) , intent(in) :: state - character(len=*) , intent(in) :: fldname - real (kind=dbl_kind) , intent(inout) :: output(:,:,:,:) - integer , intent(in) :: index - logical, optional , intent(in) :: do_sum - integer, optional , intent(in) :: ungridded_index - integer , intent(out) :: rc + type(ESMF_State) , intent(in) :: state + character(len=*) , intent(in) :: fldname + real (kind=dbl_kind) , intent(inout) :: output(:,:,:,:) + integer , intent(in) :: index + integer, optional , intent(in) :: ungridded_index + real(kind=dbl_kind), optional , intent(in) :: areacor(:) + integer , intent(out) :: rc ! local variables type(block) :: this_block ! block information for current block @@ -1330,9 +1437,7 @@ subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungr integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid - character(len=*), parameter :: subname='(ice_import_export:state_getimport)' + character(len=*), parameter :: subname='(ice_import_export:state_getimport_4d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1340,103 +1445,65 @@ subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungr ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! set values of output array - n=0 + ! set values of output array + n=0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(ungridded_index)) then + output(i,j,index,iblk) = dataPtr2d(ungridded_index,n) + else + output(i,j,index,iblk) = dataPtr1d(n) + end if + end do + end do + end do + if (present(areacor)) then + n = 0 do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi do j = jlo, jhi do i = ilo, ihi - n = n+1 - if (present(do_sum)) then ! do sum - if (present(ungridded_index)) then - output(i,j,index,iblk) = output(i,j,index, iblk) + dataPtr2d(ungridded_index,n) - else - output(i,j,index,iblk) = output(i,j,index, iblk) + dataPtr1d(n) - end if - else ! do not do sum - if (present(ungridded_index)) then - output(i,j,index,iblk) = dataPtr2d(ungridded_index,n) - else - output(i,j,index,iblk) = dataPtr1d(n) - end if - end if + n = n + 1 + output(i,j,index,iblk) = output(i,j,index,iblk) * areacor(n) end do end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! set values of output array - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(do_sum)) then - if (present(ungridded_index)) then - output(i,j,index,iblk) = output(i,j,index,iblk) + dataPtr4d(i1,j1,iblk,ungridded_index) - else - output(i,j,index,iblk) = output(i,j,index,iblk) + dataPtr3d(i1,j1,iblk) - end if - else - if (present(ungridded_index)) then - output(i,j,index,iblk) = dataPtr4d(i1,j1,iblk,ungridded_index) - else - output(i,j,index,iblk) = dataPtr3d(i1,j1,iblk) - end if - end if - end do - end do - end do - end if - end subroutine state_getimport_4d_output + end subroutine state_getimport_4d !=============================================================================== - - subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_index, rc) + subroutine state_getimport_3d(state, fldname, output, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array ! ---------------------------------------------- ! input/output variables - type(ESMF_State) , intent(in) :: state - character(len=*) , intent(in) :: fldname - real (kind=dbl_kind) , intent(inout) :: output(:,:,:) - logical, optional , intent(in) :: do_sum - integer, optional , intent(in) :: ungridded_index - integer , intent(out) :: rc + type(ESMF_State) , intent(in) :: state + character(len=*) , intent(in) :: fldname + real (kind=dbl_kind) , intent(inout) :: output(:,:,:) + integer, optional , intent(in) :: ungridded_index + real(kind=dbl_kind), optional , intent(in) :: areacor(:) + integer , intent(out) :: rc ! local variables type(block) :: this_block ! block information for current block @@ -1444,9 +1511,7 @@ subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_i integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid - character(len=*) , parameter :: subname='(ice_import_export:state_getimport)' + character(len=*) , parameter :: subname='(ice_import_export:state_getimport_3d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1454,83 +1519,53 @@ subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_i ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! determine output array - n=0 + ! determine output array + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(ungridded_index)) then + output(i,j,iblk) = dataPtr2d(ungridded_index,n) + else + output(i,j,iblk) = dataPtr1d(n) + end if + end do + end do + end do + if (present(areacor)) then + n = 0 do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi do j = jlo, jhi do i = ilo, ihi - n = n+1 - if (present(do_sum) .and. present(ungridded_index)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr2d(ungridded_index,n) - else if (present(do_sum)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr1d(n) - else if (present(ungridded_index)) then - output(i,j,iblk) = dataPtr2d(ungridded_index,n) - else - output(i,j,iblk) = dataPtr1d(n) - end if - end do - end do - end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! set values of output array - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(do_sum) .and. present(ungridded_index)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr4d(i1,j1,iblk,ungridded_index) - else if (present(do_sum)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr3d(i1,j1,iblk) - else if (present(ungridded_index)) then - output(i,j,iblk) = dataPtr4d(i1,j1,iblk, ungridded_index) - else - output(i,j,iblk) = dataPtr3d(i1,j1,iblk) - end if + n = n + 1 + output(i,j,iblk) = output(i,j,iblk) * areacor(n) end do end do end do - end if - end subroutine state_getimport_3d_output + end subroutine state_getimport_3d !=============================================================================== - - subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ungridded_index, rc) + subroutine state_setexport_4d(state, fldname, input, index, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 4d input array to export state field @@ -1544,6 +1579,7 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, logical , optional, intent(in) :: lmask(:,:,:) real(kind=dbl_kind) , optional, intent(in) :: ifrac(:,:,:) integer , optional, intent(in) :: ungridded_index + real(kind=dbl_kind) , optional, intent(in) :: areacor(:) integer , intent(out) :: rc ! local variables @@ -1552,9 +1588,8 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, integer :: i, j, iblk, n, i1, j1 ! indices real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid - character(len=*), parameter :: subname='(ice_import_export:state_setexport)' + integer :: ice_num + character(len=*), parameter :: subname='(ice_import_export:state_setexport_4d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1562,93 +1597,81 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (ungridded_index == 1) then + dataptr2d(:,:) = c0 end if - - ! set values of field pointer n = 0 do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - if (present(lmask) .and. present(ifrac)) then + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + if (present(lmask) .and. present(ifrac)) then + do j = jlo, jhi + do i = ilo, ihi + n = n+1 if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) - else - dataPtr1d(n) = input(i,j,index,iblk) - end if - end if - else - if (present(ungridded_index)) then dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) else - dataPtr1d(n) = input(i,j,index,iblk) + dataPtr2d(ungridded_index,n) = c0 end if - end if + end do end do - end do + else + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) + end do + end do + end if end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ice_num = n + if (present(areacor)) then + do n = 1,ice_num + dataPtr2d(ungridded_index,n) = dataPtr2d(ungridded_index,n) * areacor(n) + end do end if - - do iblk = 1,nblocks + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + dataptr1d(:) = c0 + n = 0 + do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(lmask) .and. present(ifrac)) then + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + if (present(lmask) .and. present(ifrac)) then + do j = jlo, jhi + do i = ilo, ihi + n = n+1 if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,index,iblk) - end if - else - dataPtr3d(i1,j1,iblk) = input(i,j,index,iblk) - end if - else - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,index,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,index,iblk) + dataPtr1d(n) = input(i,j,index,iblk) end if - end if + end do end do - end do + else + do i = ilo, ihi + n = n+1 + dataPtr1d(n) = input(i,j,index,iblk) + end do + end if end do - + ice_num = n + if (present(areacor)) then + do n = 1,ice_num + dataPtr1d(n) = dataPtr1d(n) * areacor(n) + end do + end if end if - end subroutine state_setexport_4d_input + end subroutine state_setexport_4d !=============================================================================== - - subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridded_index, rc) + subroutine state_setexport_3d(state, fldname, input, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 3d input array to export state field @@ -1661,6 +1684,7 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd logical , optional , intent(in) :: lmask(:,:,:) real(kind=dbl_kind) , optional , intent(in) :: ifrac(:,:,:) integer , optional , intent(in) :: ungridded_index + real(kind=dbl_kind) , optional , intent(in) :: areacor(:) integer , intent(out) :: rc ! local variables @@ -1669,9 +1693,8 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid - character(len=*), parameter :: subname='(ice_import_export:state_setexport)' + integer :: num_ice + character(len=*), parameter :: subname='(ice_import_export:state_setexport_3d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1679,92 +1702,59 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr2d(ungridded_index,n) = input(i,j,iblk) - else - dataPtr1d(n) = input(i,j,iblk) - end if - end if - else + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(lmask) .and. present(ifrac)) then + if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then if (present(ungridded_index)) then dataPtr2d(ungridded_index,n) = input(i,j,iblk) else dataPtr1d(n) = input(i,j,iblk) end if end if - end do + else + if (present(ungridded_index)) then + dataPtr2d(ungridded_index,n) = input(i,j,iblk) + else + dataPtr1d(n) = input(i,j,iblk) + end if + end if end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer + end do + num_ice = n + if (present(areacor)) then if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,num_ice + dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) + end do else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,iblk) - end if - end if - else - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,iblk) - end if - end if - end do + do n = 1,num_ice + dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do - end do - + end if end if - end subroutine state_setexport_3d_input + end subroutine state_setexport_3d !=============================================================================== - subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc) ! ---------------------------------------------- ! Get pointer to a state field @@ -1788,10 +1778,10 @@ subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc) call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + end subroutine State_GetFldPtr_1d !=============================================================================== - subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) ! ---------------------------------------------- ! Get pointer to a state field @@ -1815,60 +1805,7 @@ subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_2d - - !=============================================================================== - - subroutine State_GetFldPtr_3d(State, fldname, fldptr, rc) - ! ---------------------------------------------- - ! Get pointer to a state field - ! ---------------------------------------------- - - ! input/output variables - type(ESMF_State) , intent(in) :: State - character(len=*) , intent(in) :: fldname - real(kind=dbl_kind) , pointer , intent(inout) :: fldptr(:,:,:) - integer , optional , intent(out) :: rc - - ! local variables - type(ESMF_Field) :: lfield - character(len=*),parameter :: subname='(ice_import_export:State_GetFldPtr_3d)' - ! ---------------------------------------------- - - rc = ESMF_SUCCESS - - call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_3d - !=============================================================================== - - subroutine State_GetFldPtr_4d(State, fldname, fldptr, rc) - ! ---------------------------------------------- - ! Get pointer to a state field - ! ---------------------------------------------- - - ! input/output variables - type(ESMF_State) , intent(in) :: State - character(len=*) , intent(in) :: fldname - real(kind=dbl_kind) , pointer , intent(inout) :: fldptr(:,:,:,:) - integer , optional , intent(out) :: rc - - ! local variables - type(ESMF_Field) :: lfield - character(len=*),parameter :: subname='(ice_import_export:State_GetFldPtr_3d)' - ! ---------------------------------------------- - - rc = ESMF_SUCCESS - - call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_4d + end subroutine State_GetFldPtr_2d end module ice_import_export diff --git a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 new file mode 100644 index 000000000..17941435d --- /dev/null +++ b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 @@ -0,0 +1,666 @@ +module ice_mesh_mod + + use ESMF + use NUOPC , only : NUOPC_CompAttributeGet + use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long + use ice_domain_size , only : nx_global, ny_global, max_blocks + use ice_domain , only : nblocks, blocks_ice, distrb_info + use ice_blocks , only : block, get_block, nx_block, ny_block, nblocks_x, nblocks_y + use ice_shr_methods , only : chkerr + use ice_fileunits , only : nu_diag + use ice_communicate , only : my_task, master_task + use ice_exit , only : abort_ice + use icepack_intfc , only : icepack_query_parameters + use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted + + implicit none + private + + public :: ice_mesh_set_distgrid + public :: ice_mesh_setmask_from_maskfile + public :: ice_mesh_create_scolumn + public :: ice_mesh_init_tlon_tlat_area_hm + public :: ice_mesh_check + + ! Only relevant for lat-lon grids gridcell value of [1 - (land fraction)] (T-cell) + real (dbl_kind), allocatable, public :: ocn_gridcell_frac(:,:,:) + + character(*), parameter :: u_FILE_u = & + __FILE__ + +!======================================================================= +contains +!======================================================================= + + subroutine ice_mesh_set_distgrid(localpet, npes, distgrid, rc) + + ! Determine the global index space needed for the distgrid + + ! input/output variables + integer , intent(in) :: localpet + integer , intent(in) :: npes + type(ESMF_DistGrid) , intent(inout) :: distgrid + integer , intent(out) :: rc + + ! local variables + integer :: n,c,g,i,j,m ! indices + integer :: iblk, jblk ! indices + integer :: ig, jg ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + integer :: lsize ! local size of coupling array + type(block) :: this_block ! block information for current block + integer :: num_elim_global + integer :: num_elim_local + integer :: num_elim + integer :: num_ice + integer :: num_elim_gcells ! local number of eliminated gridcells + integer :: num_elim_blocks ! local number of eliminated blocks + integer :: num_total_blocks + integer :: my_elim_start, my_elim_end + integer , allocatable :: gindex(:) + integer , allocatable :: gindex_ice(:) + integer , allocatable :: gindex_elim(:) + integer :: globalID + character(len=*), parameter :: subname = ' ice_mesh_set_distgrid: ' + !---------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! number the local grid to get allocation size for gindex_ice + lsize = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + lsize = lsize + 1 + enddo + enddo + enddo + + ! set global index array + allocate(gindex_ice(lsize)) + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + ig = this_block%i_glob(i) + jg = this_block%j_glob(j) + gindex_ice(n) = (jg-1)*nx_global + ig + enddo + enddo + enddo + + ! Determine total number of eliminated blocks globally + globalID = 0 + num_elim_global = 0 ! number of eliminated blocks + num_total_blocks = 0 + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + num_total_blocks = num_total_blocks + 1 + if (distrb_info%blockLocation(globalID) == 0) then + num_elim_global = num_elim_global + 1 + end if + end do + end do + + if (num_elim_global > 0) then + + ! Distribute the eliminated blocks in a round robin fashion amoung processors + num_elim_local = num_elim_global / npes + my_elim_start = num_elim_local*localPet + min(localPet, mod(num_elim_global, npes)) + 1 + if (localPet < mod(num_elim_global, npes)) then + num_elim_local = num_elim_local + 1 + end if + my_elim_end = my_elim_start + num_elim_local - 1 + + ! Determine the number of eliminated gridcells locally + globalID = 0 + num_elim_blocks = 0 ! local number of eliminated blocks + num_elim_gcells = 0 + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + if (distrb_info%blockLocation(globalID) == 0) then + num_elim_blocks = num_elim_blocks + 1 + if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then + this_block = get_block(globalID, globalID) + num_elim_gcells = num_elim_gcells + & + (this_block%jhi-this_block%jlo+1) * (this_block%ihi-this_block%ilo+1) + end if + end if + end do + end do + + ! Determine the global index space of the eliminated gridcells + allocate(gindex_elim(num_elim_gcells)) + globalID = 0 + num_elim_gcells = 0 ! local number of eliminated gridcells + num_elim_blocks = 0 ! local number of eliminated blocks + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + if (distrb_info%blockLocation(globalID) == 0) then + this_block = get_block(globalID, globalID) + num_elim_blocks = num_elim_blocks + 1 + if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then + do j=this_block%jlo,this_block%jhi + do i=this_block%ilo,this_block%ihi + num_elim_gcells = num_elim_gcells + 1 + ig = this_block%i_glob(i) + jg = this_block%j_glob(j) + gindex_elim(num_elim_gcells) = (jg-1)*nx_global + ig + end do + end do + end if + end if + end do + end do + + ! create a global index that includes both active and eliminated gridcells + num_ice = size(gindex_ice) + num_elim = size(gindex_elim) + allocate(gindex(num_elim + num_ice)) + do n = 1,num_ice + gindex(n) = gindex_ice(n) + end do + do n = num_ice+1,num_ice+num_elim + gindex(n) = gindex_elim(n-num_ice) + end do + + deallocate(gindex_elim) + + else + + ! No eliminated land blocks + num_ice = size(gindex_ice) + allocate(gindex(num_ice)) + do n = 1,num_ice + gindex(n) = gindex_ice(n) + end do + + end if + + !--------------------------------------------------------------------------- + ! Create distGrid from global index array + !--------------------------------------------------------------------------- + + DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + deallocate(gindex_ice) + deallocate(gindex) + + end subroutine ice_mesh_set_distgrid + + !======================================================================= + subroutine ice_mesh_setmask_from_maskfile(ice_maskfile, ice_mesh, rc) + + use ice_grid , only : tlon, tlat, hm, tarea + use ice_constants , only : c0, c1, c2, p25, radius + + ! input/output variables + character(len=*) , intent(in) :: ice_maskfile + type(ESMF_Mesh) , intent(inout) :: ice_mesh + integer , intent(out) :: rc + + ! local variables + integer :: i, j, n + integer (int_kind) :: ni, nj + integer :: iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type (block) :: this_block ! block information for current block + real(dbl_kind) , pointer :: ice_frac(:) + type(ESMF_Field) :: areaField + type(ESMF_Mesh) :: mesh_mask + type(ESMF_Field) :: field_mask + type(ESMF_Field) :: field_dst + type(ESMF_RouteHandle) :: rhandle + integer :: srcMaskValue = 0 + integer :: dstMaskValue = -987987 ! spval for RH mask values + integer :: srcTermProcessing_Value = 0 + logical :: checkflag = .false. + integer, pointer :: ice_mask(:) + real(dbl_kind) , pointer :: mask_src(:) ! on mesh created from ice_maskfile + real(dbl_kind) , pointer :: dataptr1d(:) + type(ESMF_DistGrid) :: distgrid_mask + type(ESMF_Array) :: elemMaskArray + integer :: lsize_mask, lsize_dst + integer :: spatialDim + real(dbl_kind) :: fminval = 0.001_dbl_kind ! TODO: make this a share constant + real(dbl_kind) :: fmaxval = 1._dbl_kind + real(dbl_kind) :: lfrac + real(dbl_kind) , pointer :: mesh_areas(:) + integer :: numownedelements + real(dbl_kind) , pointer :: ownedElemCoords(:) + real(dbl_kind) :: pi + real(dbl_kind) :: c180 + real(dbl_kind) :: puny + real(dbl_kind) :: deg_to_rad + character(len=*), parameter :: subname = ' ice_mesh_setmask_from_maskfile' + !--------------------------------------------------- + + rc = ESMF_SUCCESS + + mesh_mask = ESMF_MeshCreate(trim(ice_maskfile), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=lsize_dst, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(ice_mask(lsize_dst)) + allocate(ice_frac(lsize_dst)) + + ! create fields on source and destination meshes + field_mask = ESMF_FieldCreate(mesh_mask, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + field_dst = ESMF_FieldCreate(ice_mesh, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create route handle to map source mask (assume ocean) to destination mesh (assume atm/lnd) + call ESMF_FieldRegridStore(field_mask, field_dst, routehandle=rhandle, & + srcMaskValues=(/srcMaskValue/), dstMaskValues=(/dstMaskValue/), & + regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=ESMF_NORMTYPE_DSTAREA, & + srcTermProcessing=srcTermProcessing_Value, & + ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! fill in values for field_mask with mask on source mesh + call ESMF_MeshGet(mesh_mask, elementdistGrid=distgrid_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distgrid_mask, localDe=0, elementCount=lsize_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(mask_src(lsize_mask)) + elemMaskArray = ESMF_ArrayCreate(distgrid_mask, mask_src, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! The following call fills in the values of mask_src + call ESMF_MeshGet(mesh_mask, elemMaskArray=elemMaskArray, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! The following call fills in the values of field_mask + call ESMF_FieldGet(field_mask, farrayptr=dataptr1d, rc=rc) + dataptr1d(:) = mask_src(:) + + ! map source mask to destination mesh - to obtain destination mask and frac + call ESMF_FieldRegrid(field_mask, field_dst, routehandle=rhandle, & + termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(field_dst, farrayptr=dataptr1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! now determine ice_mask and ice_frac + do n = 1,size(dataptr1d) + lfrac = c1 - dataptr1d(n) + if (lfrac > fmaxval) lfrac = c1 + if (lfrac < fminval) lfrac = c0 + ice_frac(n) = c1 - lfrac + if (ice_frac(n) == c0) then + ice_mask(n) = 0 + else + ice_mask(n) = 1 + end if + enddo + + ! reset the model mesh mask + call ESMF_MeshSet(ice_mesh, elementMask=ice_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! deallocate memory + call ESMF_RouteHandleDestroy(rhandle, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldDestroy(field_mask, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldDestroy(field_dst, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + deallocate(mask_src) + + ! Allocate module variable ocn_gridcell_frac + allocate(ocn_gridcell_frac(nx_block,ny_block,max_blocks)) + + ! Obtain mesh areas in radians^2 + areaField = ESMF_FieldCreate(ice_mesh, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldRegridGetArea(areaField, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(areaField, farrayPtr=mesh_areas, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Obtain mesh lons and lats in degrees + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(ownedElemCoords(spatialDim*numownedelements)) + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Get required constants + call icepack_query_parameters(pi_out=pi, c180_out=c180) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + deg_to_rad = pi/c180 + + ! Set tlon, tlat, tarea, hm + ! Convert mesh areas from radians^2 to m^2 (tarea is in m^2) + ! Convert lons and lats from degrees to radians + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n + 1 + tlon(i,j,iblk) = ownedElemCoords(2*n-1) * deg_to_rad + tlat(i,j,iblk) = ownedElemCoords(2*n) * deg_to_rad + tarea(i,j,iblk) = mesh_areas(n) * (radius*radius) + hm(i,j,iblk) = real(ice_mask(n),kind=dbl_kind) + ocn_gridcell_frac(i,j,iblk) = ice_frac(n) + enddo + enddo + enddo + + ! Dealocate memory + deallocate(ownedElemCoords) + call ESMF_FieldDestroy(areaField) + + end subroutine ice_mesh_setmask_from_maskfile + + !=============================================================================== + subroutine ice_mesh_create_scolumn(scol_lon, scol_lat, ice_mesh, rc) + + use ice_constants , only : c0, c1 + use ice_scam , only : scmlat, scmlon, scol_area, scol_mask, scol_frac, scol_nj + use netcdf + + ! Create the model mesh from the domain file - for either single column mode + ! or for a regional grid + + ! input/output variables + real(dbl_kind) , intent(in) :: scol_lon + real(dbl_kind) , intent(in) :: scol_lat + type(ESMF_Mesh) , intent(inout) :: ice_mesh + integer , intent(out) :: rc + + ! local variables + type(ESMF_Grid) :: lgrid + integer :: maxIndex(2) + real(dbl_kind) :: mincornerCoord(2) + real(dbl_kind) :: maxcornerCoord(2) + integer :: i, j,iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type (block) :: this_block ! block information for current block + character(len=*), parameter :: subname = ' ice_mesh_create_scolumn' + ! ---------------------------------------------- + + rc = ESMF_SUCCESS + + ! Use center and come up with arbitrary area delta lon and lat = .1 degree + maxIndex(1) = 1 ! number of lons + maxIndex(2) = 1 ! number of lats + mincornerCoord(1) = scol_lon - .1_dbl_kind ! min lon + mincornerCoord(2) = scol_lat - .1_dbl_kind ! min lat + maxcornerCoord(1) = scol_lon + .1_dbl_kind ! max lon + maxcornerCoord(2) = scol_lat + .1_dbl_kind ! max lat + + ! create the ESMF grid + lgrid = ESMF_GridCreateNoPeriDimUfrm (maxindex=maxindex, & + mincornercoord=mincornercoord, maxcornercoord= maxcornercoord, & + staggerloclist=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create the mesh from the lgrid + ice_mesh = ESMF_MeshCreate(lgrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Allocate module variable ocn_gridcell_frac + allocate(ocn_gridcell_frac(nx_block,ny_block,max_blocks)) + ocn_gridcell_frac(:,:,:) = scol_frac + + end subroutine ice_mesh_create_scolumn + + !=============================================================================== + subroutine ice_mesh_init_tlon_tlat_area_hm() + + use ice_grid , only : tlon, tlat, hm, tarea, ULON, ULAT, HTN, HTE, ANGLE, ANGLET + use ice_grid , only : uarea, uarear, tarear, tinyarea + use ice_grid , only : dxt, dyt, dxu, dyu, dyhx, dxhy, cyp, cxp, cym, cxm + use ice_grid , only : makemask + use ice_boundary , only : ice_HaloUpdate + use ice_domain , only : blocks_ice, nblocks, halo_info, distrb_info + use ice_constants , only : c0, c1, p25 + use ice_constants , only : field_loc_center, field_type_scalar + use ice_scam , only : scmlat, scmlon, scol_area, scol_mask, scol_frac, scol_nj, single_column + + ! local variables + integer :: i,j,n + integer :: iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type (block) :: this_block ! block information for current block + real(dbl_kind) :: puny + real(dbl_kind) :: pi + character(len=*), parameter :: subname = ' ice_mesh_init_tlon_tlat_area_hm' + ! ---------------------------------------------- + + ! Get required constants + call icepack_query_parameters(pi_out=pi, puny_out=puny) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + ! Check for consistency + if (single_column) then + if ((nx_global /= 1).or. (ny_global /= 1)) then + write(nu_diag,*) 'nx_global = ',nx_global + write(nu_diag,*) 'ny_global = ',ny_global + write(nu_diag,*) 'Because you have selected the column model flag' + write(nu_diag,*) 'then require nx_global=ny_global=1 in file ice_domain_size.F' + call abort_ice(' ice_mesh_init_tlon_tlat_area_hm: nx_global and ny_global need to be 1 for single column') + else + write(nu_diag,'(a,f10.5)')' single column mode lon/lat does contain ocn with ocn fraction ',scol_frac + end if + + TLON = scmlon + TLAT = scmlat + tarea = scol_area + hm = scol_mask + ULAT = TLAT + pi/scol_nj + end if + + call ice_HaloUpdate (TLON , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (TLAT , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (tarea , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (hm , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + + !----------------------------------------------------------------- + ! CALCULATE various geometric 2d arrays + ! The U grid (velocity) is not used when run with sequential CAM + ! because we only use thermodynamic sea ice. However, ULAT is used + ! in the default initialization of CICE so we calculate it here as + ! a "dummy" so that CICE will initialize with ice. If a no ice + ! initialization is OK (or desired) this can be commented out and + ! ULAT will remain 0 as specified above. ULAT is located at the + ! NE corner of the grid cell, TLAT at the center, so here ULAT is + ! hacked by adding half the latitudinal spacing (in radians) to TLAT. + !----------------------------------------------------------------- + + ANGLET(:,:,:) = c0 + + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + do j = jlo, jhi + do i = ilo, ihi + + if (ny_global == 1) then + uarea(i,j,iblk) = tarea(i,j, iblk) + else + uarea(i,j,iblk) = p25* & + (tarea(i,j, iblk) + tarea(i+1,j, iblk) & + + tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk)) + endif + tarear(i,j,iblk) = c1/tarea(i,j,iblk) + uarear(i,j,iblk) = c1/uarea(i,j,iblk) + tinyarea(i,j,iblk) = puny*tarea(i,j,iblk) + + if (.not. single_column) then + if (ny_global == 1) then + ULAT(i,j,iblk) = TLAT(i,j,iblk) + else + ULAT(i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global) + endif + endif + ULON (i,j,iblk) = c0 + ANGLE (i,j,iblk) = c0 + + HTN (i,j,iblk) = 1.e36_dbl_kind + HTE (i,j,iblk) = 1.e36_dbl_kind + dxt (i,j,iblk) = 1.e36_dbl_kind + dyt (i,j,iblk) = 1.e36_dbl_kind + dxu (i,j,iblk) = 1.e36_dbl_kind + dyu (i,j,iblk) = 1.e36_dbl_kind + dxhy (i,j,iblk) = 1.e36_dbl_kind + dyhx (i,j,iblk) = 1.e36_dbl_kind + cyp (i,j,iblk) = 1.e36_dbl_kind + cxp (i,j,iblk) = 1.e36_dbl_kind + cym (i,j,iblk) = 1.e36_dbl_kind + cxm (i,j,iblk) = 1.e36_dbl_kind + enddo + enddo + enddo + + call ice_HaloUpdate (ULAT, halo_info, field_loc_center, field_type_scalar, fillValue=c1) + + ! Set the boundary values for the T cell land mask (hm) and + ! make the logical land masks for T and U cells (tmask, umask). + ! Also create hemisphere masks (mask-n northern, mask-s southern) + call makemask() + + end subroutine ice_mesh_init_tlon_tlat_area_hm + + !=============================================================================== + subroutine ice_mesh_check(gcomp, ice_mesh, rc) + + ! Check CICE mesh + + use ice_constants, only : c1,c0,c360 + use ice_grid , only : tlon, tlat + + ! input/output parameters + type(ESMF_GridComp) , intent(inout) :: gcomp + type(ESMF_Mesh) , intent(inout) :: ice_mesh + integer , intent(out) :: rc + + ! local variables + type(ESMF_DistGrid) :: distGrid + integer :: n,c,g,i,j,m ! indices + integer :: iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type(block) :: this_block ! block information for current block + integer :: spatialDim + integer :: numOwnedElements + real(dbl_kind), pointer :: ownedElemCoords(:) + real(dbl_kind), pointer :: lat(:), latMesh(:) + real(dbl_kind), pointer :: lon(:), lonMesh(:) + real(dbl_kind) :: diff_lon + real(dbl_kind) :: diff_lat + real(dbl_kind) :: rad_to_deg + real(dbl_kind) :: tmplon, eps_imesh + logical :: isPresent, isSet + character(len=char_len_long) :: cvalue + character(len=char_len_long) :: logmsg + character(len=*), parameter :: subname = ' ice_mesh_check: ' + !--------------------------------------------------- + + ! Determine allowed mesh error + call NUOPC_CompAttributeGet(gcomp, name='eps_imesh', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) eps_imesh + else + eps_imesh = 1.0e-1_dbl_kind + end if + write(logmsg,*) eps_imesh + call ESMF_LogWrite(trim(subname)//' eps_imesh = '//trim(logmsg), ESMF_LOGMSG_INFO) + + ! error check differences between internally generated lons and those read in + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(ownedElemCoords(spatialDim*numownedelements)) + allocate(lonmesh(numOwnedElements)) + allocate(latmesh(numOwnedElements)) + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,numOwnedElements + lonMesh(n) = ownedElemCoords(2*n-1) + latMesh(n) = ownedElemCoords(2*n) + end do + + ! obtain internally generated cice lats and lons for error checks + call icepack_query_parameters(rad_to_deg_out=rad_to_deg) + allocate(lon(numOwnedElements)) + allocate(lat(numOwnedElements)) + lon(:) = 0. + lat(:) = 0. + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n + 1 + lon(n) = tlon(i,j,iblk)*rad_to_deg + lat(n) = tlat(i,j,iblk)*rad_to_deg + + tmplon = lon(n) + if(tmplon < c0)tmplon = tmplon + c360 + + ! error check differences between internally generated lons and those read in + diff_lon = abs(mod(lonMesh(n) - tmplon,360.0)) + if (diff_lon > eps_imesh ) then + write(6,100)n,lonMesh(n),tmplon, diff_lon + call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + end if + diff_lat = abs(latMesh(n) - lat(n)) + if (diff_lat > eps_imesh) then + write(6,101)n,latMesh(n),lat(n), diff_lat + call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + end if + + enddo + enddo + enddo + +100 format('ERROR: CICE n, lonmesh, lon, diff_lon = ',i6,2(f21.13,3x),d21.5) +101 format('ERROR: CICE n, latmesh, lat, diff_lat = ',i6,2(f21.13,3x),d21.5) + + ! deallocate memory + deallocate(ownedElemCoords) + deallocate(lon, lonMesh) + deallocate(lat, latMesh) + + end subroutine ice_mesh_check + +end module ice_mesh_mod diff --git a/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90 index 6eca4f2b4..dc40177d8 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90 @@ -7,39 +7,33 @@ module ice_prescribed_mod ! Ice/ocean fluxes are set to zero, and ice dynamics are not calculated. ! Regridding and data cycling capabilities are included. + use ESMF + #ifndef CESMCOUPLED use ice_kinds_mod - implicit none private ! except - public :: ice_prescribed_init ! initialize input data stream logical(kind=log_kind), parameter, public :: prescribed_ice = .false. ! true if prescribed ice - contains ! This is a stub routine for now - subroutine ice_prescribed_init(mpicom, compid, gindex) - integer(kind=int_kind), intent(in) :: mpicom - integer(kind=int_kind), intent(in) :: compid - integer(kind=int_kind), intent(in) :: gindex(:) + subroutine ice_prescribed_init(clock, mesh, rc) + type(ESMF_Clock) , intent(in) :: clock + type(ESMF_Mesh) , intent(in) :: mesh + integer , intent(out) :: rc ! do nothing end subroutine ice_prescribed_init -#else - - use shr_nl_mod , only : shr_nl_find_group_name - use shr_strdata_mod - use shr_dmodel_mod - use shr_string_mod - use shr_ncread_mod - use shr_sys_mod - use shr_mct_mod - use mct_mod - use pio +#else + + use ice_kinds_mod + use shr_nl_mod , only : shr_nl_find_group_name + use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_print + use dshr_strdata_mod , only : shr_strdata_init_from_inline, shr_strdata_advance + use dshr_methods_mod , only : dshr_fldbun_getfldptr use ice_broadcast use ice_communicate , only : my_task, master_task, MPI_COMM_ICE - use ice_kinds_mod use ice_fileunits use ice_exit , only : abort_ice use ice_domain_size , only : nx_global, ny_global, ncat, nilyr, nslyr, max_blocks @@ -54,306 +48,288 @@ end subroutine ice_prescribed_init use icepack_intfc , only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc , only: icepack_query_tracer_indices, icepack_query_tracer_sizes use icepack_intfc , only: icepack_query_parameters + use ice_shr_methods , only: chkerr implicit none private ! except - ! MEMBER FUNCTIONS: - public :: ice_prescribed_init ! initialize input data stream - public :: ice_prescribed_run ! get time slices and time interp - public :: ice_prescribed_phys ! set prescribed ice state and fluxes - - ! !PUBLIC DATA MEMBERS: - logical(kind=log_kind), public :: prescribed_ice ! true if prescribed ice - integer(kind=int_kind),parameter :: nFilesMaximum = 400 ! max number of files - integer(kind=int_kind) :: stream_year_first ! first year in stream to use - integer(kind=int_kind) :: stream_year_last ! last year in stream to use - integer(kind=int_kind) :: model_year_align ! align stream_year_first with this model year - character(len=char_len_long) :: stream_fldVarName - character(len=char_len_long) :: stream_fldFileName(nFilesMaximum) - character(len=char_len_long) :: stream_domTvarName - character(len=char_len_long) :: stream_domXvarName - character(len=char_len_long) :: stream_domYvarName - character(len=char_len_long) :: stream_domAreaName - character(len=char_len_long) :: stream_domMaskName - character(len=char_len_long) :: stream_domFileName - character(len=char_len_long) :: stream_mapread - logical(kind=log_kind) :: prescribed_ice_fill ! true if data fill required - type(shr_strdata_type) :: sdat ! prescribed data stream - character(len=char_len_long) :: fldList ! list of fields in data stream - real(kind=dbl_kind),allocatable :: ice_cov(:,:,:) ! ice cover + ! public member functions: + public :: ice_prescribed_init ! initialize input data stream + public :: ice_prescribed_run ! get time slices and time interp + public :: ice_prescribed_phys ! set prescribed ice state and fluxes -contains + ! public data members: + logical(kind=log_kind), public :: prescribed_ice ! true if prescribed ice - subroutine ice_prescribed_init(mpicom, compid, gindex) + ! private data members: + type(shr_strdata_type) :: sdat ! prescribed data stream + real(kind=dbl_kind),allocatable :: ice_cov(:,:,:) ! ice cover - ! Prescribed ice initialization - needed to - ! work with new shr_strdata module derived type + character(*), parameter :: u_FILE_u = & + __FILE__ - use shr_pio_mod, only : shr_pio_getiotype, shr_pio_getiosys, shr_pio_getioformat +!======================================================================= +contains +!=============================================================================== + + subroutine ice_prescribed_init(clock, mesh, rc) + + ! Prescribed ice initialization - implicit none include 'mpif.h' - ! !nput/output parameters: - integer(kind=int_kind), intent(in) :: mpicom - integer(kind=int_kind), intent(in) :: compid - integer(kind=int_kind), intent(in) :: gindex(:) + ! input/output parameters + type(ESMF_Clock) , intent(in) :: clock + type(ESMF_Mesh) , intent(in) :: mesh + integer , intent(out) :: rc + + ! local parameters + integer(kind=int_kind),parameter :: nFilesMaximum = 400 ! max number of files + integer(kind=int_kind) :: n, nFile, ierr + integer(kind=int_kind) :: nml_error ! namelist i/o error flag + character(len=char_len_long) :: stream_meshFile + character(len=char_len_long) :: stream_dataFiles(nFilesMaximum) + character(len=char_len_long) :: stream_varname + character(len=char_len_long) :: stream_mapalgo + integer(kind=int_kind) :: stream_yearfirst ! first year in stream to use + integer(kind=int_kind) :: stream_yearlast ! last year in stream to use + integer(kind=int_kind) :: stream_yearalign ! align stream_year_first + integer(kind=int_kind) :: nu_nml + logical :: prescribed_ice_mode + character(*),parameter :: subName = "('ice_prescribed_init')" + character(*),parameter :: F00 = "('(ice_prescribed_init) ',4a)" + character(*),parameter :: F01 = "('(ice_prescribed_init) ',a,i0)" + character(*),parameter :: F02 = "('(ice_prescribed_init) ',2a,i0,)" + !-------------------------------- - !----- Local ------ - type(mct_gsMap) :: gsmap_ice - type(mct_gGrid) :: dom_ice - integer(kind=int_kind) :: lsize - integer(kind=int_kind) :: gsize - integer(kind=int_kind) :: nml_error ! namelist i/o error flag - integer(kind=int_kind) :: n, nFile, ierr - character(len=8) :: fillalgo - character(*),parameter :: subName = '(ice_prescribed_init)' - - namelist /ice_prescribed_nml/ & - prescribed_ice, & - model_year_align, & - stream_year_first , & - stream_year_last , & - stream_fldVarName , & - stream_fldFileName, & - stream_domTvarName, & - stream_domXvarName, & - stream_domYvarName, & - stream_domAreaName, & - stream_domMaskName, & - stream_domFileName, & - stream_mapread, & - prescribed_ice_fill + namelist /ice_prescribed_nml/ & + prescribed_ice_mode, & + stream_meshfile, & + stream_varname , & + stream_datafiles, & + stream_mapalgo, & + stream_yearalign, & + stream_yearfirst , & + stream_yearlast + + rc = ESMF_SUCCESS ! default values for namelist - prescribed_ice = .false. ! if true, prescribe ice - stream_year_first = 1 ! first year in pice stream to use - stream_year_last = 1 ! last year in pice stream to use - model_year_align = 1 ! align stream_year_first with this model year - stream_fldVarName = 'ice_cov' - stream_fldFileName(:) = ' ' - stream_domTvarName = 'time' - stream_domXvarName = 'lon' - stream_domYvarName = 'lat' - stream_domAreaName = 'area' - stream_domMaskName = 'mask' - stream_domFileName = ' ' - stream_mapread = 'NOT_SET' - prescribed_ice_fill = .false. ! true if pice data fill required - - ! read from input file - call get_fileunit(nu_nml) + prescribed_ice_mode = .false. ! if true, prescribe ice + stream_yearfirst = 1 ! first year in pice stream to use + stream_yearlast = 1 ! last year in pice stream to use + stream_yearalign = 1 ! align stream_year_first with this model year + stream_varname = 'ice_cov' + stream_meshfile = ' ' + stream_datafiles(:) = ' ' + stream_mapalgo = 'bilinear' + + ! read namelist on master task if (my_task == master_task) then - open (nu_nml, file=nml_filename, status='old',iostat=nml_error) + open (newunit=nu_nml, file=nml_filename, status='old',iostat=nml_error) call shr_nl_find_group_name(nu_nml, 'ice_prescribed_nml', status=nml_error) - if (nml_error == 0) then - read(nu_nml, ice_prescribed_nml, iostat=nml_error) - if (nml_error > 0) then - call shr_sys_abort( 'problem on read of ice_prescribed namelist in ice_prescribed_mod' ) - endif + if (nml_error /= 0) then + write(nu_diag,F00) "ERROR: problem on read of ice_prescribed_nml namelist" + call abort_ice(subName) endif + read(nu_nml, ice_prescribed_nml, iostat=nml_error) + close(nu_nml) end if - call release_fileunit(nu_nml) - call broadcast_scalar(prescribed_ice, master_task) - - ! *** If not prescribed ice then return *** - if (.not. prescribed_ice) RETURN - - call broadcast_scalar(model_year_align,master_task) - call broadcast_scalar(stream_year_first,master_task) - call broadcast_scalar(stream_year_last,master_task) - call broadcast_scalar(stream_fldVarName,master_task) - call broadcast_scalar(stream_domTvarName,master_task) - call broadcast_scalar(stream_domXvarName,master_task) - call broadcast_scalar(stream_domYvarName,master_task) - call broadcast_scalar(stream_domAreaName,master_task) - call broadcast_scalar(stream_domMaskName,master_task) - call broadcast_scalar(stream_domFileName,master_task) - call broadcast_scalar(stream_mapread,master_task) - call broadcast_scalar(prescribed_ice_fill,master_task) - call mpi_bcast(stream_fldFileName, len(stream_fldFileName(1))*NFilesMaximum, & - MPI_CHARACTER, 0, MPI_COMM_ICE, ierr) - - nFile = 0 - do n=1,nFilesMaximum - if (stream_fldFileName(n) /= ' ') nFile = nFile + 1 - end do - ! Read shr_strdata_nml namelist - if (prescribed_ice_fill) then - fillalgo='nn' - else - fillalgo='none' - endif + ! broadcast namelist input + call broadcast_scalar(prescribed_ice_mode, master_task) - if (my_task == master_task) then - write(nu_diag,*) ' ' - write(nu_diag,*) 'This is the prescribed ice coverage option.' - write(nu_diag,*) ' stream_year_first = ',stream_year_first - write(nu_diag,*) ' stream_year_last = ',stream_year_last - write(nu_diag,*) ' model_year_align = ',model_year_align - write(nu_diag,*) ' stream_fldVarName = ',trim(stream_fldVarName) - do n = 1,nFile - write(nu_diag,*) ' stream_fldFileName = ',trim(stream_fldFileName(n)),n + ! set module variable 'prescribed_ice' + prescribed_ice = prescribed_ice_mode + + ! -------------------------------------------------- + ! only do the following if prescribed ice mode is on + ! -------------------------------------------------- + + if (prescribed_ice_mode) then + + call broadcast_scalar(stream_yearalign , master_task) + call broadcast_scalar(stream_yearfirst , master_task) + call broadcast_scalar(stream_yearlast , master_task) + call broadcast_scalar(stream_meshfile , master_task) + call broadcast_scalar(stream_mapalgo , master_task) + call broadcast_scalar(stream_varname , master_task) + call mpi_bcast(stream_dataFiles, len(stream_datafiles(1))*NFilesMaximum, MPI_CHARACTER, 0, MPI_COMM_ICE, ierr) + + nFile = 0 + do n = 1,nFilesMaximum + if (stream_datafiles(n) /= ' ') nFile = nFile + 1 end do - write(nu_diag,*) ' stream_domTvarName = ',trim(stream_domTvarName) - write(nu_diag,*) ' stream_domXvarName = ',trim(stream_domXvarName) - write(nu_diag,*) ' stream_domYvarName = ',trim(stream_domYvarName) - write(nu_diag,*) ' stream_domFileName = ',trim(stream_domFileName) - write(nu_diag,*) ' stream_mapread = ',trim(stream_mapread) - write(nu_diag,*) ' stream_fillalgo = ',trim(fillalgo) - write(nu_diag,*) ' ' - endif - - gsize = nx_global*ny_global - lsize = size(gindex) - call mct_gsMap_init( gsmap_ice, gindex, MPI_COMM_ICE, compid, lsize, gsize) - call ice_prescribed_set_domain( lsize, MPI_COMM_ICE, gsmap_ice, dom_ice ) - - call shr_strdata_create(sdat,name="prescribed_ice", & - mpicom=MPI_COMM_ICE, compid=compid, & - gsmap=gsmap_ice, ggrid=dom_ice, & - nxg=nx_global,nyg=ny_global, & - yearFirst=stream_year_first, & - yearLast=stream_year_last, & - yearAlign=model_year_align, & - offset=0, & - domFilePath='', & - domFileName=trim(stream_domFileName), & - domTvarName=stream_domTvarName, & - domXvarName=stream_domXvarName, & - domYvarName=stream_domYvarName, & - domAreaName=stream_domAreaName, & - domMaskName=stream_domMaskName, & - filePath='', & - filename=stream_fldFileName(1:nFile), & - fldListFile=stream_fldVarName, & - fldListModel=stream_fldVarName, & - fillalgo=trim(fillalgo), & - calendar=trim(calendar_type), & - mapread=trim(stream_mapread)) - if (my_task == master_task) then - call shr_strdata_print(sdat,'SPRESICE data') - endif + if (my_task == master_task) then + write(nu_diag,*) ' ' + write(nu_diag,F00) 'This is the prescribed ice coverage option.' + write(nu_diag,F01) ' stream_yearfirst = ',stream_yearfirst + write(nu_diag,F01) ' stream_yearlast = ',stream_yearlast + write(nu_diag,F01) ' stream_yearalign = ',stream_yearalign + write(nu_diag,F00) ' stream_meshfile = ',trim(stream_meshfile) + write(nu_diag,F00) ' stream_varname = ',trim(stream_varname) + write(nu_diag,F00) ' stream_mapalgo = ',trim(stream_mapalgo) + do n = 1,nFile + write(nu_diag,F00) ' stream_datafiles = ',trim(stream_dataFiles(n)) + end do + write(nu_diag,*) ' ' + endif + + ! initialize sdat + call shr_strdata_init_from_inline(sdat, & + my_task = my_task, & + logunit = nu_diag, & + compname = 'ICE', & + model_clock = clock, & + model_mesh = mesh, & + stream_meshfile = stream_meshfile, & + stream_lev_dimname = 'null', & + stream_mapalgo = trim(stream_mapalgo), & + stream_filenames = stream_datafiles(1:nfile), & + stream_fldlistFile = (/'ice_cov'/), & + stream_fldListModel = (/'ice_cov'/), & + stream_yearFirst = stream_yearFirst, & + stream_yearLast = stream_yearLast, & + stream_yearAlign = stream_yearAlign , & + stream_offset = 0, & + stream_taxmode = 'cycle', & + stream_dtlimit = 1.5_dbl_kind, & + stream_tintalgo = 'linear', & + rc = rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! print out sdat info + if (my_task == master_task) then + call shr_strdata_print(sdat,'ice coverage prescribed data') + endif + + ! For one ice category, set hin_max(1) to something big + if (ncat == 1) then + hin_max(1) = 999._dbl_kind + end if + + end if ! end of if prescribed ice mode - !----------------------------------------------------------------- - ! For one ice category, set hin_max(1) to something big - !----------------------------------------------------------------- - if (ncat == 1) then - hin_max(1) = 999._dbl_kind - end if end subroutine ice_prescribed_init !======================================================================= subroutine ice_prescribed_run(mDateIn, secIn) - ! !DESCRIPTION: - ! Finds two time slices bounding current model time, remaps if necessary - - implicit none + ! Finds two time slices bounding current model time, remaps if necessary + ! Interpolate to new ice coverage - ! !INPUT/OUTPUT PARAMETERS: - integer(kind=int_kind), intent(in) :: mDateIn ! Current model date (yyyymmdd) - integer(kind=int_kind), intent(in) :: secIn ! Elapsed seconds on model date + ! input/output parameters: + integer(kind=int_kind), intent(in) :: mDateIn ! Current model date (yyyymmdd) + integer(kind=int_kind), intent(in) :: secIn ! Elapsed seconds on model date + + ! local variables + integer(kind=int_kind) :: i,j,n,iblk ! loop indices and counter + integer(kind=int_kind) :: ilo,ihi,jlo,jhi ! beginning and end of physical domain + type (block) :: this_block + real(kind=dbl_kind) :: aice_max ! maximun ice concentration + real(kind=dbl_kind), pointer :: dataptr(:) + integer :: rc ! ESMF return code + character(*),parameter :: subName = "('ice_prescribed_run')" + character(*),parameter :: F00 = "('(ice_prescribed_run) ',a,2g20.13)" + logical :: first_time = .true. + !------------------------------------------------------------------------ - ! local varaibles - integer(kind=int_kind) :: i,j,n,iblk ! loop indices and counter - integer(kind=int_kind) :: ilo,ihi,jlo,jhi ! beginning and end of physical domain - type (block) :: this_block - real(kind=dbl_kind) :: aice_max ! maximun ice concentration - logical, save :: first_time = .true. - character(*),parameter :: subName = '(ice_prescribed_run)' - character(*),parameter :: F00 = "(a,2g20.13)" + rc = ESMF_SUCCESS - !------------------------------------------------------------------------ - ! Interpolate to new ice coverage - !------------------------------------------------------------------------ + ! Advance sdat stream + call shr_strdata_advance(sdat, ymd=mDateIn, tod=SecIn, logunit=nu_diag, istr='cice_pice', rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if - call shr_strdata_advance(sdat,mDateIn,SecIn,MPI_COMM_ICE,'cice_pice') + ! Get pointer for stream data that is time and spatially interpolate to model time and grid + call dshr_fldbun_getFldPtr(sdat%pstrm(1)%fldbun_model, 'ice_cov', dataptr, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if - if (first_time) then + ! Fill in module ice_cov array + if (.not. allocated(ice_cov)) then allocate(ice_cov(nx_block,ny_block,max_blocks)) - endif - + end if ice_cov(:,:,:) = c0 ! This initializes ghost cells as well - - n=0 + n = 0 do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo ihi = this_block%ihi jlo = this_block%jlo jhi = this_block%jhi - do j = jlo, jhi do i = ilo, ihi n = n+1 - ice_cov(i,j,iblk) = sdat%avs(1)%rAttr(1,n) + ice_cov(i,j,iblk) = dataptr(n) end do end do end do - !-------------------------------------------------------------------- ! Check to see that ice concentration is in fraction, not percent - !-------------------------------------------------------------------- if (first_time) then aice_max = maxval(ice_cov) - if (aice_max > c10) then - write(nu_diag,F00) subname//" ERROR: Ice conc data must be in fraction, aice_max= ",& - aice_max - call abort_ice(subName) + write(nu_diag,F00) "ERROR: Ice conc data must be in fraction, aice_max= ", aice_max + rc = ESMF_FAILURE + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if end if first_time = .false. end if - !----------------------------------------------------------------- ! Set prescribed ice state and fluxes - !----------------------------------------------------------------- - call ice_prescribed_phys() end subroutine ice_prescribed_run - !=============================================================================== - subroutine ice_prescribed_phys + !======================================================================= + subroutine ice_prescribed_phys() ! Set prescribed ice state using input ice concentration; ! set surface ice temperature to atmospheric value; use ! linear temperature gradient in ice to ocean temperature. - ! !USES: use ice_flux use ice_state use icepack_intfc, only : icepack_aggregate use ice_dyn_evp - implicit none !----- Local ------ integer(kind=int_kind) :: layer ! level index integer(kind=int_kind) :: nc ! ice category index integer(kind=int_kind) :: i,j,k ! longitude, latitude and level indices integer(kind=int_kind) :: iblk - integer(kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, ntrcr - - real(kind=dbl_kind) :: slope ! diff in underlying ocean tmp and ice surface tmp - real(kind=dbl_kind) :: Ti ! ice level temperature - real(kind=dbl_kind) :: Tmlt ! ice level melt temperature - real(kind=dbl_kind) :: qin_save(nilyr) - real(kind=dbl_kind) :: qsn_save(nslyr) - real(kind=dbl_kind) :: hi ! ice prescribed (hemispheric) ice thickness - real(kind=dbl_kind) :: hs ! snow thickness - real(kind=dbl_kind) :: zn ! normalized ice thickness - real(kind=dbl_kind) :: salin(nilyr) ! salinity (ppt) - real(kind=dbl_kind) :: rad_to_deg, pi, puny - real(kind=dbl_kind) :: rhoi, rhos, cp_ice, cp_ocn, lfresh, depressT - + integer(kind=int_kind) :: nt_Tsfc + integer(kind=int_kind) :: nt_sice + integer(kind=int_kind) :: nt_qice + integer(kind=int_kind) :: nt_qsno + integer(kind=int_kind) :: ntrcr + real(kind=dbl_kind) :: slope ! diff in underlying ocean tmp and ice surface tmp + real(kind=dbl_kind) :: Ti ! ice level temperature + real(kind=dbl_kind) :: Tmlt ! ice level melt temperature + real(kind=dbl_kind) :: qin_save(nilyr) + real(kind=dbl_kind) :: qsn_save(nslyr) + real(kind=dbl_kind) :: hi ! ice prescribed (hemispheric) ice thickness + real(kind=dbl_kind) :: hs ! snow thickness + real(kind=dbl_kind) :: zn ! normalized ice thickness + real(kind=dbl_kind) :: salin(nilyr) ! salinity (ppt) + real(kind=dbl_kind) :: rad_to_deg, pi, puny + real(kind=dbl_kind) :: rhoi + real(kind=dbl_kind) :: rhos + real(kind=dbl_kind) :: cp_ice + real(kind=dbl_kind) :: cp_ocn + real(kind=dbl_kind) :: lfresh + real(kind=dbl_kind) :: depressT real(kind=dbl_kind), parameter :: nsal = 0.407_dbl_kind real(kind=dbl_kind), parameter :: msal = 0.573_dbl_kind real(kind=dbl_kind), parameter :: saltmax = 3.2_dbl_kind ! max salinity at ice base (ppm) character(*),parameter :: subName = '(ice_prescribed_phys)' + !----------------------------------------------------------------- call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, & nt_qice_out=nt_qice, nt_qsno_out=nt_qsno) @@ -458,7 +434,7 @@ subroutine ice_prescribed_phys trcrn(i,j,nt_sice:nt_sice+nilyr-1,:,iblk) = c0 trcrn(i,j,nt_qice:nt_qice+nilyr-1,:,iblk) = c0 trcrn(i,j,nt_qsno:nt_qsno+nslyr-1,:,iblk) = c0 - end if ! ice_cov >= eps04 + end if ! ice_cov >= eps04 !-------------------------------------------------------------------- ! compute aggregate ice state and open water area @@ -478,10 +454,11 @@ subroutine ice_prescribed_phys trcr_base = trcr_base(1:ntrcr,:), & n_trcr_strata = n_trcr_strata(1:ntrcr), & nt_strata = nt_strata(1:ntrcr,:)) - end if ! tmask - enddo ! i - enddo ! j - enddo ! iblk + + end if ! tmask + enddo ! i + enddo ! j + enddo ! iblk do iblk = 1, nblocks do j = 1, ny_block @@ -509,105 +486,6 @@ subroutine ice_prescribed_phys end subroutine ice_prescribed_phys - !=============================================================================== - subroutine ice_prescribed_set_domain( lsize, mpicom, gsmap_i, dom_i ) - - ! Arguments - integer , intent(in) :: lsize - integer , intent(in) :: mpicom - type(mct_gsMap), intent(in) :: gsMap_i - type(mct_ggrid), intent(inout) :: dom_i - - ! Local Variables - integer :: i, j, iblk, n ! indices - integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain - real(dbl_kind), pointer :: data1(:) ! temporary - real(dbl_kind), pointer :: data2(:) ! temporary - real(dbl_kind), pointer :: data3(:) ! temporary - real(dbl_kind), pointer :: data4(:) ! temporary - real(dbl_kind), pointer :: data5(:) ! temporary - real(dbl_kind), pointer :: data6(:) ! temporary - integer , pointer :: idata(:) ! temporary - real(kind=dbl_kind) :: rad_to_deg - type(block) :: this_block ! block information for current block - character(*),parameter :: subName = '(ice_prescribed_set_domain)' - !-------------------------------- - - call icepack_query_parameters(rad_to_deg_out=rad_to_deg) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - ! Initialize mct domain type - call mct_gGrid_init(GGrid=dom_i, & - CoordChars='lat:lon:hgt', OtherChars='area:aream:mask:frac', lsize=lsize ) - call mct_aVect_zero(dom_i%data) - - ! Determine global gridpoint number attribute, GlobGridNum, which is set automatically by MCT - call mct_gsMap_orderedPoints(gsMap_i, my_task, idata) - call mct_gGrid_importIAttr(dom_i,'GlobGridNum',idata,lsize) - deallocate(idata) - - ! Determine domain (numbering scheme is: West to East and South to North to South pole) - ! Initialize attribute vector with special value - - allocate(data1(lsize)) - allocate(data2(lsize)) - allocate(data3(lsize)) - allocate(data4(lsize)) - allocate(data5(lsize)) - allocate(data6(lsize)) - - data1(:) = -9999.0_dbl_kind - data2(:) = -9999.0_dbl_kind - data3(:) = -9999.0_dbl_kind - data4(:) = -9999.0_dbl_kind - call mct_gGrid_importRAttr(dom_i,"lat" ,data1,lsize) - call mct_gGrid_importRAttr(dom_i,"lon" ,data2,lsize) - call mct_gGrid_importRAttr(dom_i,"area" ,data3,lsize) - call mct_gGrid_importRAttr(dom_i,"aream",data4,lsize) - data5(:) = 0.0_dbl_kind - data6(:) = 0.0_dbl_kind - call mct_gGrid_importRAttr(dom_i,"mask" ,data5,lsize) - call mct_gGrid_importRAttr(dom_i,"frac" ,data6,lsize) - - ! Fill in correct values for domain components - ! lat/lon in degrees, area in radians^2, mask is 1 (ocean), 0 (non-ocean) - n=0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - - data1(n) = TLON(i,j,iblk)*rad_to_deg - data2(n) = TLAT(i,j,iblk)*rad_to_deg - data3(n) = tarea(i,j,iblk)/(radius*radius) - - data5(n) = real(nint(hm(i,j,iblk)),kind=dbl_kind) - if (trim(grid_type) == 'latlon') then - data6(n) = ocn_gridcell_frac(i,j,iblk) - else - data6(n) = real(nint(hm(i,j,iblk)),kind=dbl_kind) - end if - - enddo !i - enddo !j - enddo !iblk - call mct_gGrid_importRattr(dom_i,"lon" ,data1,lsize) - call mct_gGrid_importRattr(dom_i,"lat" ,data2,lsize) - call mct_gGrid_importRattr(dom_i,"area",data3,lsize) - call mct_gGrid_importRattr(dom_i,"mask",data5,lsize) - call mct_gGrid_importRattr(dom_i,"frac",data6,lsize) - - deallocate(data1, data2, data3, data4, data5, data6) - - end subroutine ice_prescribed_set_domain - #endif end module ice_prescribed_mod diff --git a/cicecore/drivers/nuopc/cmeps/ice_scam.F90 b/cicecore/drivers/nuopc/cmeps/ice_scam.F90 index f5280b259..b92900e4f 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_scam.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_scam.F90 @@ -6,9 +6,15 @@ module ice_scam ! single column control variables (only used for latlon grid) - logical :: single_column ! true => single column mode - real (kind=dbl_kind) scmlat ! single column latitude (degrees) - real (kind=dbl_kind) scmlon ! single column longitude (degrees) + logical :: single_column = .false. ! true => single column mode + real (kind=dbl_kind) :: scmlat ! single column latitude (degrees) + real (kind=dbl_kind) :: scmlon ! single column longitude (degrees) + real (kind=dbl_kind) :: scol_frac ! single column ocn fraction + real (kind=dbl_kind) :: scol_mask ! single column ocn mask + real (kind=dbl_kind) :: scol_area ! single column ocn area + integer :: scol_ni ! ni size of single column domain file + integer :: scol_nj ! nj size of single column domain file + logical :: scol_valid = .false. ! true => single column mask is 1 end module ice_scam diff --git a/cicecore/drivers/nuopc/cmeps/ice_shr_methods.F90 b/cicecore/drivers/nuopc/cmeps/ice_shr_methods.F90 index 323cba9a4..1144568b4 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_shr_methods.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_shr_methods.F90 @@ -11,8 +11,8 @@ module ice_shr_methods use ESMF , only : ESMF_GeomType_Flag, ESMF_FieldStatus_Flag use ESMF , only : ESMF_Mesh, ESMF_MeshGet use ESMF , only : ESMF_GEOMTYPE_MESH, ESMF_GEOMTYPE_GRID, ESMF_FIELDSTATUS_COMPLETE - use ESMF , only : ESMF_Clock, ESMF_ClockCreate, ESMF_ClockGet, ESMF_ClockSet - use ESMF , only : ESMF_ClockPrint, ESMF_ClockAdvance + use ESMF , only : ESMF_Clock, ESMF_ClockCreate, ESMF_ClockGet, ESMF_ClockSet + use ESMF , only : ESMF_ClockPrint, ESMF_ClockAdvance use ESMF , only : ESMF_Alarm, ESMF_AlarmCreate, ESMF_AlarmGet, ESMF_AlarmSet use ESMF , only : ESMF_Calendar, ESMF_CALKIND_NOLEAP, ESMF_CALKIND_GREGORIAN use ESMF , only : ESMF_Time, ESMF_TimeGet, ESMF_TimeSet @@ -38,7 +38,7 @@ module ice_shr_methods public :: state_reset public :: state_flddebug public :: state_diagnose - public :: alarmInit + public :: alarmInit public :: chkerr private :: timeInit @@ -65,7 +65,7 @@ module ice_shr_methods optMonthly = "monthly" , & optYearly = "yearly" , & optDate = "date" , & - optIfdays0 = "ifdays0" + optIfdays0 = "ifdays0" ! Module data integer, parameter :: SecPerDay = 86400 ! Seconds per day @@ -588,7 +588,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) call ESMF_MeshGet(lmesh, numOwnedNodes=nnodes, numOwnedElements=nelements, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (nnodes == 0 .and. nelements == 0) lrank = 0 - else + else call ESMF_LogWrite(trim(subname)//": ERROR geomtype not supported ", & ESMF_LOGMSG_INFO, rc=rc) rc = ESMF_FAILURE @@ -949,7 +949,7 @@ end subroutine alarmInit subroutine timeInit( Time, ymd, cal, tod, rc) - ! Create the ESMF_Time object corresponding to the given input time, + ! Create the ESMF_Time object corresponding to the given input time, ! given in YMD (Year Month Day) and TOD (Time-of-day) format. ! Set the time by an integer as YYYYMMDD and integer seconds in the day diff --git a/configuration/scripts/forapps/ufs/comp_ice.backend.clean b/configuration/scripts/forapps/ufs/comp_ice.backend.clean deleted file mode 100755 index 823f1f586..000000000 --- a/configuration/scripts/forapps/ufs/comp_ice.backend.clean +++ /dev/null @@ -1,42 +0,0 @@ -#! /bin/csh -f - -### Expect to find the following environment variables set on entry: -# SITE -# SYSTEM_USERDIR -# SRCDIR -# EXEDIR - -setenv OBJDIR $EXEDIR/compile ; if !(-d $OBJDIR) mkdir -p $OBJDIR - -if (${SITE} =~ cheyenne*) then - setenv ARCH cheyenne_intel -else if (${SITE} =~ orion*) then - setenv ARCH orion_intel -else if (${SITE} =~ hera*) then - setenv ARCH hera_intel -else - echo "CICE6 ${0}: ERROR in ARCH setup, ${hname}" - exit -2 -endif - -echo "CICE6 ${0}: ARCH = $ARCH" - -cd $OBJDIR - -setenv MAKENAME gmake -setenv MAKETHRDS 1 -setenv MAKEFILE ${SRCDIR}/configuration/scripts/Makefile -setenv MACROSFILE ${SRCDIR}/configuration/scripts/machines/Macros.$ARCH - -echo "CICE6 ${0}: EXEDIR = ${EXEDIR}" -echo "CICE6 ${0}: OBJDIR = ${OBJDIR}" -echo "CICE6 ${0}: MAKEFILE = ${MAKEFILE}" -echo "CICE6 ${0}: MACROSFILE = ${MACROSFILE}" -echo "CICE6 ${0}: ESMFMKFILE = ${ESMFMKFILE}" - -#clean -${MAKENAME} EXEC=${OBJDIR}/libcice6.a \ - -f ${MAKEFILE} MACFILE=${MACROSFILE} clean - -#clean install -rm -r -f ${BINDIR} diff --git a/configuration/scripts/forapps/ufs/comp_ice.backend.libcice b/configuration/scripts/forapps/ufs/comp_ice.backend.libcice deleted file mode 100755 index ea38e048b..000000000 --- a/configuration/scripts/forapps/ufs/comp_ice.backend.libcice +++ /dev/null @@ -1,145 +0,0 @@ -#! /bin/csh -f - -### Expect to find the following environment variables set on entry: -# SITE -# SYSTEM_USERDIR -# SRCDIR -# EXEDIR - -### local variable that begin with ICE_ are needed in the Macros file -# ICE_COMMDIR -# ICE_BLDDEBUG -# ICE_THREADED -# ICE_CPPDEFS - -setenv OBJDIR $EXEDIR/compile ; if !(-d $OBJDIR) mkdir -p $OBJDIR - -setenv THRD no # set to yes for OpenMP threading - -if (${SITE} =~ cheyenne*) then - setenv ARCH cheyenne_intel -else if (${SITE} =~ orion*) then - setenv ARCH orion_intel -else if (${SITE} =~ hera*) then - setenv ARCH hera_intel -else - echo "CICE6 ${0}: ERROR in ARCH setup, ${hname}" - exit -2 -endif - -echo "CICE6 ${0}: ARCH = $ARCH" - -cd $OBJDIR - -setenv SHRDIR csm_share # location of CCSM shared code -setenv DRVDIR nuopc/cmeps - -#if ($NTASK == 1) then -# setenv ICE_COMMDIR serial -#else - setenv ICE_COMMDIR mpi -#endif - -if ($THRD == 'yes') then - setenv ICE_THREADED true -else - setenv ICE_THREADED false -endif - -if ($?ICE_CPPDEFS) then - setenv ICE_CPPDEFS "${ICE_CPPDEFS} -Dcoupled" -else - setenv ICE_CPPDEFS "-Dcoupled" -endif - -if !($?IO_TYPE) then - setenv IO_TYPE netcdf4 # set to none if netcdf library is unavailable -endif -if ($IO_TYPE == 'netcdf3' || $IO_TYPE == 'netcdf4') then - setenv IODIR io_netcdf - setenv ICE_CPPDEFS "${ICE_CPPDEFS} -DUSE_NETCDF" -else if ($IO_TYPE == 'pio') then - setenv IODIR io_pio - setenv ICE_CPPDEFS "${ICE_CPPDEFS} -DUSE_NETCDF" -else - setenv IODIR io_binary -endif - -# Build in debug mode. If DEBUG=Y, enable DEBUG compilation. This -# flag is set in ${ROOTDIR}/coupledFV3_MOM6_CICE_debug.appBuilder file. -if (! $?DEBUG) then - setenv ICE_BLDDEBUG false -else - if ($DEBUG == "Y") then - setenv ICE_BLDDEBUG true - else - setenv ICE_BLDDEBUG false - endif -endif -echo "CICE6 ${0}: DEBUG = ${ICE_BLDDEBUG}" - -### List of source code directories (in order of importance). -cat >! Filepath << EOF -${SRCDIR}/cicecore/drivers/${DRVDIR} -${SRCDIR}/cicecore/cicedynB/dynamics -${SRCDIR}/cicecore/cicedynB/general -${SRCDIR}/cicecore/cicedynB/analysis -${SRCDIR}/cicecore/cicedynB/infrastructure -${SRCDIR}/cicecore/cicedynB/infrastructure/io/${IODIR} -${SRCDIR}/cicecore/cicedynB/infrastructure/comm/${ICE_COMMDIR} -${SRCDIR}/cicecore/shared -${SRCDIR}/icepack/columnphysics -${SRCDIR}/$SHRDIR -EOF - -setenv MAKENAME gmake -setenv MAKETHRDS 1 -setenv MAKEFILE ${SRCDIR}/configuration/scripts/Makefile -setenv MACROSFILE ${SRCDIR}/configuration/scripts/machines/Macros.$ARCH -setenv DEPFILE ${SRCDIR}/configuration/scripts/makdep.c - -echo "CICE6 ${0}: EXEDIR = ${EXEDIR}" -echo "CICE6 ${0}: OBJDIR = ${OBJDIR}" -echo "CICE6 ${0}: MAKEFILE = ${MAKEFILE}" -echo "CICE6 ${0}: MACROSFILE = ${MACROSFILE}" -echo "CICE6 ${0}: DEPFILE = ${DEPFILE}" -echo "CICE6 ${0}: ESMFMKFILE = ${ESMFMKFILE}" - -#diagnostics -#${MAKENAME} -j ${MAKETHRDS} VPFILE=Filepath EXEC=${OBJDIR}/cice \ -# -f ${MAKEFILE} MACFILE=${MACROSFILE} DEPFILE=${DEPFILE} db_files -#${MAKENAME} -j ${MAKETHRDS} VPFILE=Filepath EXEC=${OBJDIR}/cice \ -# -f ${MAKEFILE} MACFILE=${MACROSFILE} DEPFILE=${DEPFILE} db_flags - -#clean -#${MAKENAME} VPFILE=Filepath EXEC=${OBJDIR}/cice \ -# -f ${MAKEFILE} MACFILE=${MACROSFILE} DEPFILE=${DEPFILE} clean - -#needed to trigger a failed build to rest of system -rm ${BINDIR}/cice6.mk - -#build lib (includes dependencies) -${MAKENAME} -j ${MAKETHRDS} VPFILE=Filepath EXEC=${OBJDIR}/libcice6.a \ - -f ${MAKEFILE} MACFILE=${MACROSFILE} DEPFILE=${DEPFILE} libcice - -if ($status != 0) then - echo "CICE6 ${0}: gmake failed, exiting" - exit -2 -endif - -#install -mkdir -p ${BINDIR} -cp -f ${OBJDIR}/libcice6.a ${BINDIR}/ -cp -f ${OBJDIR}/ice_comp_nuopc.mod ${BINDIR}/ -cp -f ${OBJDIR}/ice_timers.mod ${BINDIR}/ - -cat >! ${BINDIR}/cice6.mk << EOF -# ESMF self-describing build dependency makefile fragment - -ESMF_DEP_FRONT = ice_comp_nuopc -ESMF_DEP_INCPATH = ${BINDIR} -ESMF_DEP_CMPL_OBJS = -ESMF_DEP_LINK_OBJS = ${BINDIR}/libcice6.a - -EOF - diff --git a/configuration/scripts/machines/Macros.hera_intel b/configuration/scripts/machines/Macros.hera_intel index 230f43e70..caad25ead 100644 --- a/configuration/scripts/machines/Macros.hera_intel +++ b/configuration/scripts/machines/Macros.hera_intel @@ -4,11 +4,11 @@ CPP := fpp CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} -CFLAGS := -c -O2 -fp-model precise -xHost +CFLAGS := -c -O2 -fp-model precise -xHost FIXEDFLAGS := -132 FREEFLAGS := -FR -FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost +FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost FFLAGS_NOOPT:= -O0 ifeq ($(ICE_BLDDEBUG), true) @@ -48,9 +48,9 @@ INCLDIR := $(INCLDIR) -I$(INC_NETCDF) #SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -L$(LIB_PNETCDF) -lpnetcdf -lgptl SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -ifeq ($(ICE_THREADED), true) - LDFLAGS += -qopenmp - CFLAGS += -qopenmp - FFLAGS += -qopenmp +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp endif diff --git a/configuration/scripts/machines/Macros.orion_intel b/configuration/scripts/machines/Macros.orion_intel index 6dffdd0a2..fa6745e03 100644 --- a/configuration/scripts/machines/Macros.orion_intel +++ b/configuration/scripts/machines/Macros.orion_intel @@ -4,11 +4,11 @@ CPP := fpp CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} -CFLAGS := -c -O2 -fp-model precise -xHost +CFLAGS := -c -O2 -fp-model precise -xHost FIXEDFLAGS := -132 FREEFLAGS := -FR -FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost +FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost FFLAGS_NOOPT:= -O0 ifeq ($(ICE_BLDDEBUG), true) @@ -48,9 +48,9 @@ INCLDIR := $(INCLDIR) -I$(INC_NETCDF) #SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -L$(LIB_PNETCDF) -lpnetcdf -lgptl SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -ifeq ($(ICE_THREADED), true) - LDFLAGS += -qopenmp - CFLAGS += -qopenmp - FFLAGS += -qopenmp +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp endif diff --git a/configuration/scripts/machines/Macros.stampede_intel b/configuration/scripts/machines/Macros.stampede_intel new file mode 100644 index 000000000..14bbc7a4a --- /dev/null +++ b/configuration/scripts/machines/Macros.stampede_intel @@ -0,0 +1,56 @@ +#============================================================================== +# Makefile macros for TACC stampede, intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise -xHost + +FIXEDFLAGS := -132 +FREEFLAGS := -FR +FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost +FFLAGS_NOOPT:= -O0 + +ifeq ($(ICE_BLDDEBUG), true) + FFLAGS += -O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created -link_mpi=dbg +else + FFLAGS += -O2 +endif + +SCC := icc +SFC := ifort +MPICC := mpiicc +MPIFC := mpiifort + +ifeq ($(ICE_COMMDIR), mpi) + FC := $(MPIFC) + CC := $(MPICC) +else + FC := $(SFC) + CC := $(SCC) +endif +LD:= $(FC) + +NETCDF_PATH := $(NETCDF_ROOT) + +PIO_CONFIG_OPTS:= --enable-filesystem-hints=gpfs + +#PNETCDF_PATH := $(PNETCDF) +#PNETCDF_PATH := /glade/u/apps/ch/opt/pio/2.2/mpt/2.15f/intel/17.0.1/lib + +INC_NETCDF := $(NETCDF_PATH)/include +LIB_NETCDF := $(NETCDF_PATH)/lib + +#LIB_PNETCDF := $(PNETCDF_PATH)/lib +#LIB_MPI := $(IMPILIBDIR) + +INCLDIR := $(INCLDIR) -I$(INC_NETCDF) +#SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -L$(LIB_PNETCDF) -lpnetcdf -lgptl +SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp +endif + diff --git a/configuration/scripts/machines/Macros.wcoss_dell_p3_intel b/configuration/scripts/machines/Macros.wcoss_dell_p3_intel new file mode 100644 index 000000000..a835be424 --- /dev/null +++ b/configuration/scripts/machines/Macros.wcoss_dell_p3_intel @@ -0,0 +1,49 @@ +#============================================================================== +# Makefile macros for wcoss phase3 machine, intel compiler +#============================================================================== + +CPP := fpp +CPPDEFS := -DFORTRANUNDERSCORE ${ICE_CPPDEFS} +CFLAGS := -c -O2 -fp-model precise -xHost + +FIXEDFLAGS := -132 +FREEFLAGS := -FR +FFLAGS := -fp-model precise -convert big_endian -assume byterecl -ftz -traceback -align array64byte -xHost +FFLAGS_NOOPT:= -O0 + +ifeq ($(ICE_BLDDEBUG), true) + FFLAGS += -O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created -link_mpi=dbg +else + FFLAGS += -O2 +endif + +SCC := icc +SFC := ifort +MPICC := mpiicc +MPIFC := mpiifort + +ifeq ($(ICE_COMMDIR), mpi) + FC := $(MPIFC) + CC := $(MPICC) +else + FC := $(SFC) + CC := $(SCC) +endif +LD:= $(FC) + +NETCDF_PATH := $(NETCDF) + +PIO_CONFIG_OPTS:= --enable-filesystem-hints=gpfs + +INC_NETCDF := $(NETCDF_PATH)/include +LIB_NETCDF := $(NETCDF_PATH)/lib + +INCLDIR := $(INCLDIR) -I$(INC_NETCDF) +#SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -L$(LIB_PNETCDF) -lpnetcdf -lgptl +SLIBS := -L$(LIB_NETCDF) -lnetcdf -lnetcdff + +ifeq ($(ICE_THREADED), true) + LDFLAGS += -qopenmp + CFLAGS += -qopenmp + FFLAGS += -qopenmp +endif