diff --git a/cicecore/cicedynB/general/ice_state.F90 b/cicecore/cicedynB/general/ice_state.F90 index d5c115a0c..4a3f83ce2 100644 --- a/cicecore/cicedynB/general/ice_state.F90 +++ b/cicecore/cicedynB/general/ice_state.F90 @@ -186,6 +186,8 @@ subroutine alloc_state n_trcr_strata = 0 nt_strata = 0 trcr_base = c0 + aicen = c0 + aicen_init = c0 end subroutine alloc_state diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index f9b5116d0..3d5e5cc2a 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -78,7 +78,7 @@ subroutine cice_init2() 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_domain_size , only: ncat, nfsd, nfreq 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 @@ -94,10 +94,12 @@ subroutine cice_init2() 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 + use ice_arrays_column , only: wavefreq, dwavefreq logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow character(len=char_len) :: snw_aging_table + real(kind=dbl_kind), dimension(25) :: wave_spectrum_profile ! hardwire for now character(len=*), parameter :: subname = '(cice_init2)' !---------------------------------------------------- @@ -177,6 +179,11 @@ subroutine cice_init2() endif endif + if (wave_spec) then + call icepack_init_wave(nfreq=nfreq, & + wave_spectrum_profile=wave_spectrum_profile, wavefreq=wavefreq, dwavefreq=dwavefreq) + end if + ! Initialize shortwave components using swdn from previous timestep ! if restarting. These components will be scaled to current forcing ! in prep_radiation. diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 8920ea386..182308973 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -25,6 +25,7 @@ module ice_comp_nuopc use ice_calendar , only : force_restart_now, write_ic use ice_calendar , only : idate, mday, mmonth, myear, year_init use ice_calendar , only : msec, dt, calendar, calendar_type, nextsw_cday, istep + use ice_calendar , only : ice_calendar_noleap, ice_calendar_gregorian use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit @@ -80,9 +81,6 @@ module ice_comp_nuopc character(len=*) , parameter :: orb_variable_year = 'variable_year' character(len=*) , parameter :: orb_fixed_parameters = 'fixed_parameters' - 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 @@ -216,7 +214,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) 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) @@ -339,7 +336,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call get_component_instance(gcomp, inst_suffix, inst_index, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - inst_name = "ICE"//trim(inst_suffix) +! inst_name = "ICE"//trim(inst_suffix) + inst_name = "ICE" !---------------------------------------------------------------------------- ! start cice timers @@ -470,9 +468,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ESMF_TimeGet( currTime, calkindflag=esmf_caltype, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (esmf_caltype == ESMF_CALKIND_NOLEAP) then - calendar_type = shr_cal_noleap + calendar_type = ice_calendar_noleap else if (esmf_caltype == ESMF_CALKIND_GREGORIAN) then - calendar_type = shr_cal_gregorian + calendar_type = ice_calendar_gregorian else call abort_ice( subname//'ERROR:: bad calendar for ESMF' ) end if @@ -581,9 +579,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) 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)) + write(errmsg,'(a)') trim(subname)//'WARNING: tfrz_option from driver '//trim(tfrz_option_driver)//& + ' is overwriting tfrz_option from cice namelist '//trim(tfrz_option) + write(nu_diag,*) trim(errmsg) + call icepack_warnings_flush(nu_diag) + call icepack_init_parameters(tfrz_option_in=tfrz_option_driver) endif ! Flux convergence tolerance - always use the driver attribute value @@ -594,7 +594,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) 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 ',& + 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) diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 7bfc53f45..f5a00cdf8 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -3,12 +3,13 @@ module ice_import_export use ESMF use NUOPC use NUOPC_Model - use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind + use ice_kinds_mod , only : int_kind, dbl_kind, char_len, char_len_long, log_kind 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 use ice_domain_size , only : nx_global, ny_global, block_size_x, block_size_y, max_blocks, ncat + use ice_domain_size , only : nfreq, nfsd use ice_exit , only : abort_ice use ice_flux , only : strairxT, strairyT, strocnxT, strocnyT use ice_flux , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref @@ -23,9 +24,10 @@ module ice_import_export 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_arrays_column , only : floe_rad_c, wave_spectrum + use ice_state , only : vice, vsno, aice, aicen_init, trcr, trcrn use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm - use ice_grid , only : grid_type, grid_average_X2Y + use ice_grid , only : grid_type use ice_mesh_mod , only : ocn_gridcell_frac use ice_boundary , only : ice_HaloUpdate use ice_fileunits , only : nu_diag, flush_fileunit @@ -34,8 +36,10 @@ module ice_import_export use ice_shr_methods , only : chkerr, state_reset 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 use icepack_intfc , only : icepack_liquidus_temperature use icepack_intfc , only : icepack_sea_freezing_temperature + use icepack_parameters , only : puny, c2 use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED use shr_frz_mod , only : shr_frz_freezetemp @@ -112,6 +116,7 @@ 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_wave ! use case logical :: isPresent, isSet character(len=*), parameter :: subname='(ice_import_export:ice_advertise_fields)' !------------------------------------------------------------------------------- @@ -148,6 +153,17 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam write(nu_diag,*)'flds_wiso = ',flds_wiso end if + flds_wave = .false. + call NUOPC_CompAttributeGet(gcomp, name='wav_coupling_to_cice', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) flds_wave + end if + if (my_task == master_task) then + write(nu_diag,*)'flds_wave = ',flds_wave + end if + !----------------- ! advertise import fields !----------------- @@ -192,6 +208,14 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam ! from atm - dry dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstdry', ungridded_lbound=1, ungridded_ubound=4) + ! the following are advertised but might not be connected if they are not advertised in the + ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific + ! from wave + if (flds_wave) then + call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & + ungridded_ubound=25) + end if + do n = 1,fldsToIce_num call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, & TransferOfferGeomObject='will provide', rc=rc) @@ -225,6 +249,10 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam call fldlist_add(fldsFrIce_num, fldsFrIce, 'ice_fraction_n', & ungridded_lbound=1, ungridded_ubound=ncat) end if + if (flds_wave) then + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_thick' ) + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_floediam' ) + end if ! ice/atm fluxes computed by ice call fldlist_add(fldsFrIce_num, fldsFrIce, 'stress_on_air_ice_zonal' ) @@ -292,7 +320,7 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc type(ESMF_State) :: exportState type(ESMF_Field) :: lfield integer :: numOwnedElements - integer :: i, j, iblk, n + integer :: i, j, iblk, n, k 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(:) @@ -403,11 +431,10 @@ subroutine ice_import( importState, rc ) ! local variables integer,parameter :: nflds=16 integer,parameter :: nfldv=6 - integer :: i, j, iblk, n + integer :: i, j, iblk, n, k integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain type(block) :: this_block ! block information for current block real (kind=dbl_kind),allocatable :: aflds(:,:,:,:) - real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: work real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh @@ -559,6 +586,29 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO + ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) + if (State_FldChk(importState, 'Sw_elevation_spectrum')) then + if (nfreq /= 25) then + call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") + end if + call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do k = 1,nfreq + 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 + wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) + end do + end do + end do + end do + end if + 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 @@ -845,7 +895,7 @@ subroutine ice_export( exportState, rc ) ! local variables type(block) :: this_block ! block information for current block - integer :: i, j, iblk, n ! incides + integer :: i, j, iblk, n, k ! indices integer :: n2 ! thickness category index integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain real (kind=dbl_kind) :: workx, worky ! tmps for converting grid @@ -859,7 +909,11 @@ 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) :: floediam(nx_block,ny_block,max_blocks) + real (kind=dbl_kind) :: floethick(nx_block,ny_block,max_blocks) ! ice thickness real (kind=dbl_kind) :: Tffresh + logical (kind=log_kind) :: tr_fsd + integer (kind=int_kind) :: nt_fsd real (kind=dbl_kind), allocatable :: tempfld(:,:,:) real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:) real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:) @@ -877,6 +931,9 @@ subroutine ice_export( exportState, rc ) ! 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_tracer_indices(nt_fsd_out=nt_fsd) + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -890,8 +947,10 @@ subroutine ice_export( exportState, rc ) tauya(:,:,:) = c0 tauxo(:,:,:) = c0 tauyo(:,:,:) = c0 + floediam(:,:,:) = c0 + floethick(:,:,:) = c0 - !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky, this_block, ilo, ihi, jlo, jhi) + !$OMP PARALLEL DO PRIVATE(iblk,i,j,k,workx,worky, this_block, ilo, ihi, jlo, jhi) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -904,6 +963,27 @@ subroutine ice_export( exportState, rc ) ! ice fraction ailohi(i,j,iblk) = min(aice(i,j,iblk), c1) + if (tr_fsd) then + ! floe thickness (m) + if (aice(i,j,iblk) > puny) then + floethick(i,j,iblk) = vice(i,j,iblk) / aice(i,j,iblk) + else + floethick(i,j,iblk) = c0 + end if + + ! floe diameter (m) + workx = c0 + worky = c0 + do n = 1, ncat + do k = 1, nfsd + workx = workx + floe_rad_c(k) * aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + worky = worky + aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + end do + end do + if (worky > c0) workx = c2*workx / worky + floediam(i,j,iblk) = MAX(c2*floe_rad_c(1),workx) + endif + ! surface temperature Tsrf(i,j,iblk) = Tffresh + trcr(i,j,1,iblk) !Kelvin (original ???) @@ -1054,6 +1134,22 @@ subroutine ice_export( exportState, rc ) call state_setexport(exportState, 'Si_snowh' , input=tempfld , lmask=tmask, ifrac=ailohi, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------ + ! optional floe diameter and ice thickness to wave + ! ------ + + ! Sea ice thickness (m) + if (State_FldChk(exportState, 'Si_thick')) then + call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Sea ice floe diameter (m) + if (State_FldChk(exportState, 'Si_floediam')) then + call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + ! ------ ! ice/atm fluxes computed by ice ! ------