diff --git a/cicecore/cicedyn/analysis/ice_diagnostics.F90 b/cicecore/cicedyn/analysis/ice_diagnostics.F90 index 8879d6632..53631b2d4 100644 --- a/cicecore/cicedyn/analysis/ice_diagnostics.F90 +++ b/cicecore/cicedyn/analysis/ice_diagnostics.F90 @@ -87,6 +87,8 @@ module ice_diagnostics totms , & ! total ice/snow water mass (sh) totmin , & ! total ice water mass (nh) totmis , & ! total ice water mass (sh) + totsn , & ! total salt mass (nh) + totss , & ! total salt mass (sh) toten , & ! total ice/snow energy (J) totes ! total ice/snow energy (J) @@ -154,7 +156,7 @@ subroutine runtime_diags (dt) rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh character (len=char_len) :: & - snwredist + snwredist, saltflux_option ! hemispheric state quantities real (kind=dbl_kind) :: & @@ -162,6 +164,8 @@ subroutine runtime_diags (dt) umaxs, hmaxs, shmaxs, areas, snwmxs, extents, shmaxst, & etotn, mtotn, micen, msnwn, pmaxn, ketotn, & etots, mtots, mices, msnws, pmaxs, ketots, & + stotn, & + stots, & urmsn, albtotn, arean_alb, mpndn, ptotn, spondn, & urmss, albtots, areas_alb, mpnds, ptots, sponds @@ -226,7 +230,7 @@ subroutine runtime_diags (dt) awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, & rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, & ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, & - snwgrain_out=snwgrain) + snwgrain_out=snwgrain, saltflux_option_out=saltflux_option) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) @@ -512,6 +516,15 @@ subroutine runtime_diags (dt) etots = global_sum(work1, distrb_info, & field_loc_center, tareas) + ! total salt volume + call total_salt (work2) + + stotn = global_sum(work2, distrb_info, & + field_loc_center, tarean) + stots = global_sum(work2, distrb_info, & + field_loc_center, tareas) + + !----------------------------------------------------------------- ! various fluxes !----------------------------------------------------------------- @@ -785,12 +798,22 @@ subroutine runtime_diags (dt) swerrs = (fswnets - fswdns) / (fswnets - c1) ! salt mass - msltn = micen*ice_ref_salinity*p001 - mslts = mices*ice_ref_salinity*p001 + if (saltflux_option == 'prognostic') then + ! compute the total salt mass + msltn = stotn*rhoi*p001 + mslts = stots*rhoi*p001 + + ! change in salt mass + delmsltn = rhoi*(stotn-totsn)*p001 + delmslts = rhoi*(stots-totss)*p001 + else + msltn = micen*ice_ref_salinity*p001 + mslts = mices*ice_ref_salinity*p001 - ! change in salt mass - delmsltn = delmxn*ice_ref_salinity*p001 - delmslts = delmxs*ice_ref_salinity*p001 + ! change in salt mass + delmsltn = delmxn*ice_ref_salinity*p001 + delmslts = delmxs*ice_ref_salinity*p001 + endif ! salt error serrn = (sfsaltn + delmsltn) / (msltn + c1) @@ -1275,7 +1298,7 @@ subroutine init_mass_diags rhoi, rhos, rhofresh real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - work1 + work1, work2 character(len=*), parameter :: subname = '(init_mass_diags)' @@ -1310,6 +1333,12 @@ subroutine init_mass_diags toten = global_sum(work1, distrb_info, field_loc_center, tarean) totes = global_sum(work1, distrb_info, field_loc_center, tareas) + ! north/south salt + call total_salt (work2) + totsn = global_sum(work2, distrb_info, field_loc_center, tarean) + totss = global_sum(work2, distrb_info, field_loc_center, tareas) + + if (print_points) then do n = 1, npnt if (my_task == pmloc(n)) then diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90 index 2142310b9..9ba5cf4d4 100644 --- a/cicecore/cicedyn/analysis/ice_history.F90 +++ b/cicecore/cicedyn/analysis/ice_history.F90 @@ -420,10 +420,6 @@ subroutine init_hist (dt) f_taubyE = f_tauby endif -#ifndef ncdf - f_bounds = .false. -#endif - ! write dimensions for 3D or 4D history variables ! note: list of variables checked here is incomplete if (f_aicen(1:1) /= 'x' .or. f_vicen(1:1) /= 'x' .or. & @@ -2168,12 +2164,13 @@ subroutine accum_hist (dt) real (kind=dbl_kind) :: awtvdr, awtidr, awtvdf, awtidf, puny, secday, rad_to_deg real (kind=dbl_kind) :: Tffresh, rhoi, rhos, rhow, ice_ref_salinity - real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt + real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt, sicen logical (kind=log_kind) :: formdrag, skl_bgc logical (kind=log_kind) :: tr_pond, tr_aero, tr_brine, tr_snow integer (kind=int_kind) :: ktherm integer (kind=int_kind) :: nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY, nt_Tsfc, & nt_alvl, nt_vlvl + character (len=char_len) :: saltflux_option type (block) :: & this_block ! block information for current block @@ -2185,6 +2182,7 @@ subroutine accum_hist (dt) call icepack_query_parameters(Tffresh_out=Tffresh, rhoi_out=rhoi, rhos_out=rhos, & rhow_out=rhow, ice_ref_salinity_out=ice_ref_salinity) call icepack_query_parameters(formdrag_out=formdrag, skl_bgc_out=skl_bgc, ktherm_out=ktherm) + call icepack_query_parameters(saltflux_option_out=saltflux_option) call icepack_query_tracer_flags(tr_pond_out=tr_pond, tr_aero_out=tr_aero, & tr_brine_out=tr_brine, tr_snow_out=tr_snow) call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_qice_out=nt_qice, & @@ -2269,7 +2267,7 @@ subroutine accum_hist (dt) !--------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, & - !$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt, & + !$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt,sicen, & !$OMP worka,workb,worka3,Tinz4d,Sinz4d,Tsnz4d) do iblk = 1, nblocks @@ -3228,7 +3226,16 @@ subroutine accum_hist (dt) dfresh = -rhoi*frazil(i,j,iblk)/dt endif endif - dfsalt = ice_ref_salinity*p001*dfresh + if (saltflux_option == 'prognostic') then + sicen = c0 + do k = 1, nzilyr + sicen = sicen + trcr(i,j,nt_sice+k-1,iblk)*vice(i,j,iblk) & + / real(nzilyr,kind=dbl_kind) + enddo + dfsalt = sicen*p001*dfresh + else + dfsalt = ice_ref_salinity*p001*dfresh + endif worka(i,j) = aice(i,j,iblk)*(fsalt(i,j,iblk)+dfsalt) endif enddo diff --git a/cicecore/cicedyn/analysis/ice_history_shared.F90 b/cicecore/cicedyn/analysis/ice_history_shared.F90 index d9c62edde..70aa5e14c 100644 --- a/cicecore/cicedyn/analysis/ice_history_shared.F90 +++ b/cicecore/cicedyn/analysis/ice_history_shared.F90 @@ -40,7 +40,7 @@ module ice_history_shared logical (kind=log_kind), public :: & hist_avg ! if true, write averaged data instead of snapshots - character (len=char_len), public :: & + character (len=char_len_long), public :: & history_file , & ! output file for history incond_file ! output file for snapshot initial conditions diff --git a/cicecore/cicedyn/dynamics/ice_dyn_evp.F90 b/cicecore/cicedyn/dynamics/ice_dyn_evp.F90 index 69305e131..cf111cccf 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_evp.F90 @@ -945,7 +945,8 @@ subroutine evp (dt) uarea (:,:,iblk), DminTarea (:,:,iblk), & strength (:,:,iblk), shearU (:,:,iblk), & zetax2T (:,:,iblk), etax2T (:,:,iblk), & - stresspT (:,:,iblk), stressmT (:,:,iblk)) + stresspT (:,:,iblk), stressmT (:,:,iblk), & + stress12T (:,:,iblk)) enddo !$OMP END PARALLEL DO @@ -1730,7 +1731,8 @@ subroutine stressC_T (nx_block, ny_block , & uarea , DminTarea , & strength , shearU , & zetax2T , etax2T , & - stressp , stressm ) + stresspT , stressmT , & + stress12T) use ice_dyn_shared, only: strain_rates_T, capping, & visc_replpress, e_factor @@ -1744,24 +1746,25 @@ subroutine stressC_T (nx_block, ny_block , & indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & - uvelE , & ! x-component of velocity (m/s) at the E point - vvelE , & ! y-component of velocity (m/s) at the E point - uvelN , & ! x-component of velocity (m/s) at the N point - vvelN , & ! y-component of velocity (m/s) at the N point - dxN , & ! width of N-cell through the middle (m) - dyE , & ! height of E-cell through the middle (m) - dxT , & ! width of T-cell through the middle (m) - dyT , & ! height of T-cell through the middle (m) - strength , & ! ice strength (N/m) - shearU , & ! shearU local for this routine - uarea , & ! area of u cell - DminTarea ! deltaminEVP*tarea + uvelE , & ! x-component of velocity (m/s) at the E point + vvelE , & ! y-component of velocity (m/s) at the E point + uvelN , & ! x-component of velocity (m/s) at the N point + vvelN , & ! y-component of velocity (m/s) at the N point + dxN , & ! width of N-cell through the middle (m) + dyE , & ! height of E-cell through the middle (m) + dxT , & ! width of T-cell through the middle (m) + dyT , & ! height of T-cell through the middle (m) + strength , & ! ice strength (N/m) + shearU , & ! shearU local for this routine + uarea , & ! area of u cell + DminTarea ! deltaminEVP*tarea real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & - zetax2T , & ! zetax2 = 2*zeta (bulk viscosity) - etax2T , & ! etax2 = 2*eta (shear viscosity) - stressp , & ! sigma11+sigma22 - stressm ! sigma11-sigma22 + zetax2T , & ! zetax2 = 2*zeta (bulk viscosity) + etax2T , & ! etax2 = 2*eta (shear viscosity) + stresspT , & ! sigma11+sigma22 + stressmT , & ! sigma11-sigma22 + stress12T ! sigma12 ! local variables @@ -1769,12 +1772,14 @@ subroutine stressC_T (nx_block, ny_block , & i, j, ij real (kind=dbl_kind), dimension (nx_block,ny_block) :: & - divT , & ! divergence at T point - tensionT ! tension at T point + divT , & ! divergence at T point + tensionT ! tension at T point real (kind=dbl_kind) :: & shearTsqr , & ! strain rates squared at T point + shearT , & ! strain rate at T point DeltaT , & ! delt at T point + uareaavgr , & ! 1 / uarea avg rep_prsT ! replacement pressure at T point character(len=*), parameter :: subname = '(stressC_T)' @@ -1801,11 +1806,19 @@ subroutine stressC_T (nx_block, ny_block , & ! U point values (Bouillon et al., 2013, Kimmritz et al., 2016 !----------------------------------------------------------------- + uareaavgr = c1/(uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j)) + shearTsqr = (shearU(i ,j )**2 * uarea(i ,j ) & + shearU(i ,j-1)**2 * uarea(i ,j-1) & + shearU(i-1,j-1)**2 * uarea(i-1,j-1) & + shearU(i-1,j )**2 * uarea(i-1,j )) & - / (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j)) + * uareaavgr + + shearT = (shearU(i ,j ) * uarea(i ,j ) & + + shearU(i ,j-1) * uarea(i ,j-1) & + + shearU(i-1,j-1) * uarea(i-1,j-1) & + + shearU(i-1,j ) * uarea(i-1,j )) & + * uareaavgr DeltaT = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr)) @@ -1822,11 +1835,14 @@ subroutine stressC_T (nx_block, ny_block , & ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - stressp(i,j) = (stressp(i,j)*(c1-arlx1i*revp) & - + arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1 + stresspT(i,j) = (stresspT (i,j)*(c1-arlx1i*revp) & + + arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1 - stressm(i,j) = (stressm(i,j)*(c1-arlx1i*revp) & - + arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1 + stressmT(i,j) = (stressmT (i,j)*(c1-arlx1i*revp) & + + arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1 + + stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) & + + arlx1i*p5*etax2T(i,j)*shearT ) * denom1 enddo ! ij @@ -1851,7 +1867,7 @@ subroutine stressC_U (nx_block , ny_block ,& uarea , & etax2U , deltaU ,& strengthU, shearU ,& - stress12) + stress12U) use ice_dyn_shared, only: visc_replpress, & visc_method, deltaminEVP, capping @@ -1872,7 +1888,7 @@ subroutine stressC_U (nx_block , ny_block ,& strengthU ! ice strength at the U point real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & - stress12 ! sigma12 + stress12U ! sigma12 ! local variables @@ -1891,15 +1907,15 @@ subroutine stressC_U (nx_block , ny_block ,& ! viscosities and replacement pressure at U point ! avg_zeta: Bouillon et al. 2013, C1 method of Kimmritz et al. 2016 ! avg_strength: C2 method of Kimmritz et al. 2016 - ! if outside do and stress12 equation repeated in each loop for performance + ! if outside do and stress12U equation repeated in each loop for performance !----------------------------------------------------------------- if (visc_method == 'avg_zeta') then do ij = 1, icellU i = indxUi(ij) j = indxUj(ij) - stress12(i,j) = (stress12(i,j)*(c1-arlx1i*revp) & - + arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1 + stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) & + + arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1 enddo elseif (visc_method == 'avg_strength') then @@ -1911,8 +1927,8 @@ subroutine stressC_U (nx_block , ny_block ,& ! minimal extra calculations here even though it seems like there is call visc_replpress (strengthU(i,j), DminUarea, deltaU(i,j), & lzetax2U , letax2U , lrep_prsU , capping) - stress12(i,j) = (stress12(i,j)*(c1-arlx1i*revp) & - + arlx1i*p5*letax2U*shearU(i,j)) * denom1 + stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) & + + arlx1i*p5*letax2U*shearU(i,j)) * denom1 enddo endif @@ -1976,7 +1992,7 @@ subroutine stressCD_T (nx_block, ny_block , & real (kind=dbl_kind), dimension (nx_block,ny_block) :: & divT , & ! divergence at T point tensionT , & ! tension at T point - shearT , & ! sheat at T point + shearT , & ! shear at T point DeltaT ! delt at T point real (kind=dbl_kind) :: & @@ -1991,7 +2007,7 @@ subroutine stressCD_T (nx_block, ny_block , & call strain_rates_T (nx_block , ny_block , & icellT , & - indxTi (:), indxTj (:) , & + indxTi(:) , indxTj (:) , & uvelE (:,:), vvelE (:,:), & uvelN (:,:), vvelN (:,:), & dxN (:,:), dyE (:,:), & @@ -2046,8 +2062,7 @@ subroutine stressCD_U (nx_block, ny_block, & stresspU , stressmU, & stress12U) - use ice_dyn_shared, only: strain_rates_U, & - visc_replpress, & + use ice_dyn_shared, only: visc_replpress, & visc_method, deltaminEVP, capping integer (kind=int_kind), intent(in) :: & @@ -2059,7 +2074,7 @@ subroutine stressCD_U (nx_block, ny_block, & indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & - uarea , & ! area of U-cell (m^2) + uarea , & ! area of U-cell (m^2) zetax2U , & ! 2*zeta at U point etax2U , & ! 2*eta at U point strengthU, & ! ice strength at U point diff --git a/cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90 b/cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90 index e874611bd..b7daab0a0 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90 @@ -889,7 +889,7 @@ subroutine evp1d_halo_update(NAVEL_len, lb, ub, uvel, vvel, & #ifdef _OPENACC !$acc parallel & - !$acc present(uvel, vvel) & + !$acc present(uvel, vvel) !$acc loop do iw = 1, NAVEL_len if (halo_parent(iw) == 0) cycle diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index a1f5aafc3..f51425780 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -17,11 +17,11 @@ module ice_init use ice_constants, only: c0, c1, c2, c3, c5, c12, p01, p2, p3, p5, p75, p166, & cm_to_m use ice_exit, only: abort_ice - use ice_fileunits, only: nu_nml, nu_diag, nu_diag_set, nml_filename, diag_type, & + use ice_fileunits, only: nu_nml, nu_diag, nml_filename, diag_type, & ice_stdout, get_fileunit, release_fileunit, bfbflag, flush_fileunit, & ice_IOUnitsMinUnit, ice_IOUnitsMaxUnit #ifdef CESMCOUPLED - use ice_fileunits, only: inst_suffix + use ice_fileunits, only: inst_suffix, nu_diag_set #endif use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_aggregate @@ -151,7 +151,7 @@ subroutine input_data kitd, kcatbound, ktransport character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & - tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table, & + tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table, & capping_method, snw_ssp_table logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, & @@ -163,7 +163,7 @@ subroutine input_data integer (kind=int_kind) :: numin, numax ! unit number limits integer (kind=int_kind) :: rplvl, rptopo - real (kind=dbl_kind) :: Cf, ksno, puny, Tocnfrz + real (kind=dbl_kind) :: Cf, ksno, puny, ice_ref_salinity, Tocnfrz character (len=char_len) :: abort_list character (len=128) :: tmpstr2 @@ -260,6 +260,7 @@ subroutine input_data highfreq, natmiter, atmiter_conv, calc_dragio, & ustar_min, emissivity, iceruf, iceruf_ocn, & fbot_xfer_type, update_ocn_f, l_mpond_fresh, tfrz_option, & + saltflux_option,ice_ref_salinity, & oceanmixed_ice, restore_ice, restore_ocn, trestore, & precip_units, default_season, wave_spec_type,nfreq, & atm_data_type, ocn_data_type, bgc_data_type, fe_data_type, & @@ -499,6 +500,8 @@ subroutine input_data precip_units = 'mks' ! 'mm_per_month' or ! 'mm_per_sec' = 'mks' = kg/m^2 s tfrz_option = 'mushy' ! freezing temp formulation + saltflux_option = 'constant' ! saltflux calculation + ice_ref_salinity = 4.0_dbl_kind ! Ice reference salinity for coupling oceanmixed_ice = .false. ! if true, use internal ocean mixed layer wave_spec_type = 'none' ! type of wave spectrum forcing nfreq = 25 ! number of wave frequencies @@ -761,8 +764,8 @@ subroutine input_data ! each task gets unique ice log filename when if test is true, for debugging if (1 == 0) then call get_fileUnit(nu_diag) - write(tmpstr,'(a,i4.4)') "ice.log.task_",my_task - open(nu_diag,file=tmpstr) + write(tmpstr2,'(a,i4.4)') "ice.log.task_",my_task + open(nu_diag,file=tmpstr2) endif end if if (trim(ice_ic) /= 'default' .and. & @@ -984,6 +987,8 @@ subroutine input_data call broadcast_scalar(wave_spec_file, master_task) call broadcast_scalar(nfreq, master_task) call broadcast_scalar(tfrz_option, master_task) + call broadcast_scalar(saltflux_option, master_task) + call broadcast_scalar(ice_ref_salinity, master_task) call broadcast_scalar(ocn_data_format, master_task) call broadcast_scalar(bgc_data_type, master_task) call broadcast_scalar(fe_data_type, master_task) @@ -1421,6 +1426,12 @@ subroutine input_data write(nu_diag,*) subname//' WARNING: For consistency, set tfrz_option = mushy' endif endif + if (ktherm == 1 .and. trim(saltflux_option) /= 'constant') then + if (my_task == master_task) then + write(nu_diag,*) subname//' WARNING: ktherm = 1 and saltflux_option = ',trim(saltflux_option) + write(nu_diag,*) subname//' WARNING: For consistency, set saltflux_option = constant' + endif + endif !tcraig if (ktherm == 1 .and. .not.sw_redist) then if (my_task == master_task) then @@ -1987,6 +1998,10 @@ subroutine input_data write(nu_diag,*) ' WARNING: will impact ocean forcing interaction' write(nu_diag,*) ' WARNING: coupled forcing will be modified by mixed layer routine' endif + write(nu_diag,1030) ' saltflux_option = ', trim(saltflux_option) + if (trim(saltflux_option) == 'constant') then + write(nu_diag,1002) ' ice_ref_salinity = ',ice_ref_salinity + endif if (trim(tfrz_option(1:8)) == 'constant') then tmpstr2 = ' : constant ocean freezing temperature (Tocnfrz)' elseif (trim(tfrz_option(1:8)) == 'minus1p8') then @@ -2397,6 +2412,7 @@ subroutine input_data wave_spec_type_in = wave_spec_type, & wave_spec_in=wave_spec, nfreq_in=nfreq, & tfrz_option_in=tfrz_option, kalg_in=kalg, fbot_xfer_type_in=fbot_xfer_type, & + saltflux_option_in=saltflux_option, ice_ref_salinity_in=ice_ref_salinity, & Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, iceruf_ocn_in=iceruf_ocn, calc_dragio_in=calc_dragio, & windmin_in=windmin, drhosdwind_in=drhosdwind, & rsnw_fall_in=rsnw_fall, rsnw_tmax_in=rsnw_tmax, rhosnew_in=rhosnew, & @@ -2813,7 +2829,7 @@ subroutine set_state_var (nx_block, ny_block, & indxi, indxj ! compressed indices for cells with aicen > puny real (kind=dbl_kind) :: & - Tsfc, sum, hbar, abar, puny, rhos, Lfresh, rad_to_deg, rsnw_fall, dist_ratio + Tsfc, sum, hbar, abar, puny, rhos, Lfresh, rad_to_deg, rsnw_fall, dist_ratio, Tffresh real (kind=dbl_kind), dimension(ncat) :: & ainit, hinit ! initial area, thickness @@ -2855,7 +2871,7 @@ subroutine set_state_var (nx_block, ny_block, & nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & nt_rhos_out=nt_rhos, nt_rsnw_out=nt_rsnw) call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny, & - rad_to_deg_out=rad_to_deg, rsnw_fall_out=rsnw_fall) + rad_to_deg_out=rad_to_deg, rsnw_fall_out=rsnw_fall, Tffresh_out=Tffresh) call icepack_query_parameters(secday_out=secday, pi_out=pi) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & @@ -3093,7 +3109,12 @@ subroutine set_state_var (nx_block, ny_block, & do i = ilo, ihi if (tmask(i,j)) then ! place ice in high latitudes where ocean sfc is cold +#ifdef CESMCOUPLED + ! Option to use Tair instead. + if ( (Tair (i,j) <= Tffresh) .and. & +#else if ( (sst (i,j) <= Tf(i,j)+p2) .and. & +#endif (TLAT(i,j) < edge_init_sh/rad_to_deg .or. & TLAT(i,j) > edge_init_nh/rad_to_deg) ) then icells = icells + 1 diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index dfccdd413..d193eca02 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -507,23 +507,26 @@ subroutine init_grid2 ! Diagnose OpenMP thread schedule, force order in output !----------------------------------------------------------------- +! This code does not work in CESM. Needs to be investigated further. +#ifndef CESMCOUPLED #if defined (_OPENMP) - !$OMP PARALLEL DO ORDERED PRIVATE(iblk) SCHEDULE(runtime) - do iblk = 1, nblocks - if (my_task == master_task) then - !$OMP ORDERED - if (iblk == 1) then - call omp_get_schedule(ompsk,ompcs) - write(nu_diag,*) '' - write(nu_diag,*) subname,' OpenMP runtime thread schedule:' - write(nu_diag,*) subname,' omp schedule = ',ompsk,ompcs - endif - write(nu_diag,*) subname,' block, thread = ',iblk,OMP_GET_THREAD_NUM() - call flush_fileunit(nu_diag) - !$OMP END ORDERED - endif - enddo - !$OMP END PARALLEL DO + !$OMP PARALLEL DO ORDERED PRIVATE(iblk) SCHEDULE(runtime) + do iblk = 1, nblocks + if (my_task == master_task) then + !$OMP ORDERED + if (iblk == 1) then + call omp_get_schedule(ompsk,ompcs) +! write(nu_diag,*) '' + write(nu_diag,*) subname,' OpenMP runtime thread schedule:' + write(nu_diag,*) subname,' omp schedule = ',ompsk,ompcs + endif + write(nu_diag,*) subname,' block, thread = ',iblk,OMP_GET_THREAD_NUM() + call flush_fileunit(nu_diag) + !$OMP END ORDERED + endif + enddo + !$OMP END PARALLEL DO +#endif #endif !----------------------------------------------------------------- diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 index 10fcf8b81..7019f7128 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 @@ -749,10 +749,6 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, & ! if (ndim3 == ncat .and. ncat>1) then if (ndim3 == ncat .and. ndims == 3) then call pio_read_darray(File, vardesc, iodesc3d_ncat, work, status) -!#ifndef CESM1_PIO -!! This only works for PIO2 -! where (work == PIO_FILL_DOUBLE) work = c0 -!#endif if (present(field_loc)) then do n=1,ndim3 call ice_HaloUpdate (work(:,:,n,:), halo_info, & @@ -762,10 +758,6 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, & ! elseif (ndim3 == 1) then elseif (ndim3 == 1 .and. ndims == 2) then call pio_read_darray(File, vardesc, iodesc2d, work, status) -!#ifndef CESM1_PIO -!! This only works for PIO2 -! where (work == PIO_FILL_DOUBLE) work = c0 -!#endif if (present(field_loc)) then call ice_HaloUpdate (work(:,:,1,:), halo_info, & field_loc, field_type) diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index 59d13e91a..3934212b5 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -107,8 +107,8 @@ subroutine cice_init2() call init_hist (dt) ! initialize output history file if (kdyn == 1) then - call init_evp ! allocate dyn_evp arrays - if (kdyn == 2) then + call init_evp ! define evp dynamics parameters, variables + elseif (kdyn == 2) then call init_eap ! define eap dynamics parameters, variables else if (kdyn == 3) then call init_vp ! define vp dynamics parameters, variables diff --git a/cicecore/shared/ice_fileunits.F90 b/cicecore/shared/ice_fileunits.F90 index 1854dda64..c8ca3a937 100644 --- a/cicecore/shared/ice_fileunits.F90 +++ b/cicecore/shared/ice_fileunits.F90 @@ -80,8 +80,10 @@ module ice_fileunits integer (kind=int_kind), public :: & nu_diag = ice_stdout ! diagnostics output file, unit number may be overwritten +#ifdef CESMCOUPLED logical (kind=log_kind), public :: & nu_diag_set = .false. ! flag to indicate whether nu_diag is already set +#endif integer (kind=int_kind), public :: & ice_IOUnitsMinUnit = 11, & ! do not use unit numbers below @@ -116,7 +118,11 @@ subroutine init_fileunits ice_IOUnitsInUse(ice_stdout) = .true. ! reserve unit 6 ice_IOUnitsInUse(ice_stderr) = .true. if (nu_diag >= 1 .and. nu_diag <= ice_IOUnitsMaxUnit) & - ice_IOUnitsInUse(nu_diag) = .true. ! reserve unit nu_diag + ice_IOUnitsInUse(nu_diag) = .true. ! reserve unit nu_diag +#ifdef CESMCOUPLED + ! CESM can have negative unit numbers. + if (nu_diag < 0) nu_diag_set = .true. +#endif call get_fileunit(nu_grid) call get_fileunit(nu_kmt) @@ -239,7 +245,12 @@ subroutine release_all_fileunits call release_fileunit(nu_rst_pointer) call release_fileunit(nu_history) call release_fileunit(nu_hdr) +#ifdef CESMCOUPLED + ! CESM can have negative unit numbers + if (nu_diag > 0 .and. nu_diag /= ice_stdout) call release_fileunit(nu_diag) +#else if (nu_diag /= ice_stdout) call release_fileunit(nu_diag) +#endif end subroutine release_all_fileunits diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 1a6123fff..d51f114e0 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -257,6 +257,8 @@ update_ocn_f = .false. l_mpond_fresh = .false. tfrz_option = 'mushy' + saltflux_option = 'constant' + ice_ref_salinity = 4.0 oceanmixed_ice = .true. wave_spec_type = 'none' wave_spec_file = 'unknown_wave_spec_file' diff --git a/configuration/scripts/options/set_nml.saltflux b/configuration/scripts/options/set_nml.saltflux new file mode 100644 index 000000000..d50ddc4e3 --- /dev/null +++ b/configuration/scripts/options/set_nml.saltflux @@ -0,0 +1,2 @@ + ktherm = 2 + saltflux_option = 'prognostic' diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index b659fcb19..8685ab9a8 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -22,6 +22,7 @@ restart gx3 4x4 alt04 restart gx3 4x4 alt05 restart gx3 8x2 alt06 restart gx3 8x3 alt07 +restart gx3 8x3 saltflux restart gx3 18x2 debug,maskhalo restart gx3 6x2 alt01,debug,short restart gx3 8x2 alt02,debug,short diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 1829b0f1e..664437cb0 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -350,7 +350,7 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "ice_ic", "choice of initial conditions (see :ref:`tab-ic`)", "" "ice_stdout", "unit number for standard output", "" "ice_stderr", "unit number for standard error output", "" - "ice_ref_salinity", "reference salinity for ice–ocean exchanges", "4. ppt" + "ice_ref_salinity", "reference salinity for ice–ocean exchanges", "" "icells", "number of grid cells with specified property (for vectorization)", "" "iceruf", "ice surface roughness at atmosphere interface", "5.\ :math:`\times`\ 10\ :math:`^{-4}` m" "iceruf_ocn", "under-ice roughness (at ocean interface)", "0.03 m" @@ -675,6 +675,7 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "Tf", "freezing temperature", "C" "Tffresh", "freezing temp of fresh ice", "273.15 K" "tfrz_option", "form of ocean freezing temperature", "" + "saltflux_option", "form of coupled salt flux ", "" "thinS", "minimum ice thickness for brine tracer", "" "timer_stats", "logical to turn on extra timer statistics", ".false." "timesecs", "total elapsed time in seconds", "s" diff --git a/doc/source/science_guide/sg_horiztrans.rst b/doc/source/science_guide/sg_horiztrans.rst index f85f13ee5..d66046465 100644 --- a/doc/source/science_guide/sg_horiztrans.rst +++ b/doc/source/science_guide/sg_horiztrans.rst @@ -35,8 +35,11 @@ versions but have not yet been implemented. Two transport schemes are available: upwind and the incremental remapping scheme of :cite:`Dukowicz00` as modified for sea ice by -:cite:`Lipscomb04`. The upwind scheme is naturally suited for a C grid discretization. As such, the C grid velocity components (i.e. :math:`uvelE=u` at the E point and :math:`vvelN=v` at the N point) are directly passed to the upwind transport scheme. On the other hand, if the B grid is used, :math:`uvel` and :math:`vvel` (respectively :math:`u` and :math:`v` at the U point) are interpolated to the E and N points such that the upwind advection can be performed. Conversely, as the remapping scheme was originally developed for B grid applications, :math:`uvel` and :math:`vvel` are directly used for the advection. If the remapping scheme is used for the C grid, :math:`uvelE` and :math:`vvelN` are first interpolated to the U points before performing the advection. +:cite:`Lipscomb04`. +- The upwind scheme uses velocity points at the East and North face (i.e. :math:`uvelE=u` at the E point and :math:`vvelN=v` at the N point) of a T gridcell. As such, the prognostic C grid velocity components (:math:`uvelE` and :math:`vvelN`) can be passed directly to the upwind transport scheme. If the upwind scheme is used with the B grid, the B grid velocities, :math:`uvelU` and :math:`vvelU` (respectively :math:`u` and :math:`v` at the U point) are interpolated to the E and N points first. (Note however that the upwind scheme does not transport all potentially available tracers.) + +- The remapping scheme uses :math:`uvelU` and :math:`vvelU` if l_fixed_area is false and :math:`uvelE` and :math:`vvelN` if l_fixed_area is true. l_fixed_area is hardcoded to false by default and further described below. As such, the B grid velocities (:math:`uvelU` and :math:`vvelU`) are used directly in the remapping scheme, while the C grid velocities (:math:`uvelE` and :math:`vvelN`) are interpolated to U points first. If l_fixed_area is changed to true, then the reverse is true. The C grid velocities are used directly and the B grid velocities are interpolated. The remapping scheme has several desirable features: @@ -464,14 +467,14 @@ In general, the fluxes in this expression are not equal to those implied by the above scheme for locating departure regions. For some applications it may be desirable to prescribe the divergence by prescribing the area of the departure region for each edge. This can be -done in CICE 4.0 by setting `l\_fixed\_area` = true in +done by setting `l\_fixed\_area` = true in **ice\_transport\_driver.F90** and passing the prescribed departure areas (`edgearea\_e` and `edgearea\_n`) into the remapping routine. An extra triangle is then constructed for each departure region to ensure that the total area is equal to the prescribed value. This idea was suggested and first implemented by Mats Bentsen of the Nansen Environmental and Remote Sensing Center (Norway), who applied an earlier version of the -CICE remapping scheme to an ocean model. The implementation in CICE v4.0 +CICE remapping scheme to an ocean model. The implementation in CICE is somewhat more general, allowing for departure regions lying on both sides of a cell edge. The extra triangle is constrained to lie in one but not both of the grid cells that share the edge. Since this option diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index dd6ec4bcd..d8fa4c7c7 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -651,6 +651,7 @@ forcing_nml "", "``eastblock``", "ice block covering about 25 percent of domain at the east edge of the domain", "" "", "``latsst``", "ice dependent on latitude and ocean temperature", "" "", "``uniform``", "ice defined at all grid points", "" + "``ice_ref_salinity``", "real", "sea ice salinity for coupling fluxes (ppt)", "4.0" "``iceruf``", "real", "ice surface roughness at atmosphere interface in meters", "0.0005" "``l_mpond_fresh``", "``.false.``", "release pond water immediately to ocean", "``.false.``" "", "``true``", "retain (topo) pond water until ponds drain", "" @@ -673,6 +674,8 @@ forcing_nml "``restore_ocn``", "logical", "restore sst to data", "``.false.``" "``restore_ice``", "logical", "restore ice state along lateral boundaries", "``.false.``" "``rotate_wind``", "logical", "rotate wind from east/north to computation grid", "``.true.``" + "``saltflux_option``", "``constant``", "computed using ice_ref_salinity", "``constant``" + "", "``prognostic``", "computed using prognostic salinity", "" "``tfrz_option``","``constant``", "constant ocean freezing temperature (Tocnfrz)","``mushy``" "", "``linear_salt``", "linear function of salinity (ktherm=1)", "", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8^{\circ} C`)", "" diff --git a/icepack b/icepack index 4d131da97..a9fc4fafa 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 4d131da975b17e8df50c227d5272840b39a11483 +Subproject commit a9fc4fafa83460a8b8582e5a016d97259f0bada8