! SVN:$Id: CICE_RunMod.F90 1304 2015-11-12 19:28:14Z afrobert@nps.edu $ !======================================================================= ! ! Main driver for time stepping of CICE. ! ! authors Elizabeth C. Hunke, LANL ! Philip W. Jones, LANL ! William H. Lipscomb, LANL ! ! 2006 ECH: moved exit timeLoop to prevent execution of unnecessary timestep ! 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 module CICE_RunMod use ice_kinds_mod implicit none private public :: CICE_Run, ice_step save !======================================================================= contains !======================================================================= ! ! This is the main driver routine for advancing CICE forward in time. ! ! author Elizabeth C. Hunke, LANL ! Philip W. Jones, LANL ! William H. Lipscomb, LANL subroutine CICE_Run use ice_aerosol, only: faero_default use ice_algae, only: get_forcing_bgc use ice_calendar, only: istep, istep1, time, dt, stop_now, calendar use ice_forcing, only: get_forcing_atmo, get_forcing_ocn use ice_flux, only: init_flux_atm, init_flux_ocn use ice_state, only: tr_aero use ice_timers, only: ice_timer_start, ice_timer_stop, & timer_couple, timer_step use ice_zbgc_shared, only: skl_bgc !-------------------------------------------------------------------- ! initialize error code and step timer !-------------------------------------------------------------------- call ice_timer_start(timer_step) ! start timing entire run #ifndef CICE_IN_NEMO !-------------------------------------------------------------------- ! timestep loop !-------------------------------------------------------------------- timeLoop: do #endif call ice_step istep = istep + 1 ! update time step counters istep1 = istep1 + 1 time = time + dt ! determine the time and date call calendar(time) ! at the end of the timestep #ifndef CICE_IN_NEMO if (stop_now >= 1) exit timeLoop #endif #ifndef coupled call ice_timer_start(timer_couple) ! atm/ocn coupling call get_forcing_atmo ! atmospheric forcing from data call get_forcing_ocn(dt) ! ocean forcing from data ! if (tr_aero) call faero_data ! aerosols if (tr_aero) call faero_default ! aerosols if (skl_bgc) call get_forcing_bgc ! biogeochemistry call ice_timer_stop(timer_couple) ! atm/ocn coupling #endif call init_flux_atm ! initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler #ifndef CICE_IN_NEMO enddo timeLoop #endif !-------------------------------------------------------------------- ! end of timestep loop !-------------------------------------------------------------------- call ice_timer_stop(timer_step) ! end timestepping loop timer end subroutine CICE_Run !======================================================================= ! ! Calls drivers for physics components, some initialization, and output ! ! author Elizabeth C. Hunke, LANL ! William H. Lipscomb, LANL subroutine ice_step use ice_age, only: write_restart_age use ice_aerosol, only: write_restart_aero use ice_boundary, only: ice_HaloUpdate use ice_brine, only: hbrine_diags, write_restart_hbrine use ice_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep use ice_constants, only: field_loc_center, field_type_scalar use ice_diagnostics, only: init_mass_diags, runtime_diags use ice_domain, only: halo_info, nblocks use ice_domain_size, only: nslyr use ice_dyn_eap, only: write_restart_eap use ice_dyn_shared, only: kdyn use ice_firstyear, only: write_restart_FY use ice_flux, only: scale_factor, init_history_therm use ice_history, only: accum_hist use ice_lvl, only: write_restart_lvl use ice_restart, only: final_restart use ice_restart_driver, only: dumpfile use ice_meltpond_cesm, only: write_restart_pond_cesm use ice_meltpond_lvl, only: write_restart_pond_lvl use ice_meltpond_topo, only: write_restart_pond_topo use ice_restoring, only: restore_ice, ice_HaloRestore use ice_state, only: nt_qsno, trcrn, tr_iage, tr_FY, tr_lvl, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo, tr_brine, tr_aero use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, & post_thermo, step_dynamics, step_radiation use ice_therm_shared, only: calc_Tsfc use ice_timers, only: ice_timer_start, ice_timer_stop, & timer_diags, timer_column, timer_thermo, timer_bound, & timer_hist, timer_readwrite use ice_algae, only: bgc_diags, write_restart_bgc use ice_zbgc, only: init_history_bgc, biogeochemistry use ice_zbgc_shared, only: skl_bgc integer (kind=int_kind) :: & iblk , & ! block index k ! dynamics supercycling index #ifdef ICE_DA !--------------------------------------------------------------- ! Update CICE state variables using data assimilation increments !--------------------------------------------------------------- call da_state_update #endif !----------------------------------------------------------------- ! restoring on grid boundaries !----------------------------------------------------------------- if (restore_ice) call ice_HaloRestore !----------------------------------------------------------------- ! initialize diagnostics !----------------------------------------------------------------- call ice_timer_start(timer_diags) ! diagnostics/history call init_mass_diags ! diagnostics per timestep call init_history_therm call init_history_bgc call ice_timer_stop(timer_diags) ! diagnostics/history call ice_timer_start(timer_column) ! column physics call ice_timer_start(timer_thermo) ! thermodynamics !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks !----------------------------------------------------------------- ! Scale radiation fields !----------------------------------------------------------------- if (calc_Tsfc) call prep_radiation (dt, iblk) !----------------------------------------------------------------- ! thermodynamics !----------------------------------------------------------------- call step_therm1 (dt, iblk) ! vertical thermodynamics call biogeochemistry (dt, iblk) ! biogeochemistry call step_therm2 (dt, iblk) ! ice thickness distribution thermo enddo ! iblk !$OMP END PARALLEL DO call post_thermo (dt) ! finalize thermo update call ice_timer_stop(timer_thermo) ! thermodynamics call ice_timer_stop(timer_column) ! column physics !----------------------------------------------------------------- ! dynamics, transport, ridging !----------------------------------------------------------------- do k = 1, ndtd call step_dynamics (dt_dyn, ndtd) enddo !----------------------------------------------------------------- ! albedo, shortwave radiation !----------------------------------------------------------------- call ice_timer_start(timer_column) ! column physics call ice_timer_start(timer_thermo) ! thermodynamics !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call step_radiation (dt, iblk) !----------------------------------------------------------------- ! get ready for coupling and the next time step !----------------------------------------------------------------- call coupling_prep (iblk) enddo ! iblk !$OMP END PARALLEL DO call ice_timer_start(timer_bound) call ice_HaloUpdate (scale_factor, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) call ice_timer_stop(timer_thermo) ! thermodynamics call ice_timer_stop(timer_column) ! column physics !----------------------------------------------------------------- ! write data !----------------------------------------------------------------- call ice_timer_start(timer_diags) ! diagnostics if (mod(istep,diagfreq) == 0) then call runtime_diags(dt) ! log file if (skl_bgc) call bgc_diags (dt) if (tr_brine) call hbrine_diags (dt) endif call ice_timer_stop(timer_diags) ! diagnostics call ice_timer_start(timer_hist) ! history call accum_hist (dt) ! history file call ice_timer_stop(timer_hist) ! history call ice_timer_start(timer_readwrite) ! reading/writing if (write_restart == 1) then call dumpfile ! core variables for restarting if (tr_iage) call write_restart_age if (tr_FY) call write_restart_FY if (tr_lvl) call write_restart_lvl if (tr_pond_cesm) call write_restart_pond_cesm if (tr_pond_lvl) call write_restart_pond_lvl if (tr_pond_topo) call write_restart_pond_topo if (tr_aero) call write_restart_aero if (skl_bgc) call write_restart_bgc if (tr_brine) call write_restart_hbrine if (kdyn == 2) call write_restart_eap call final_restart endif call ice_timer_stop(timer_readwrite) ! reading/writing end subroutine ice_step !======================================================================= ! ! Prepare for coupling ! ! authors: Elizabeth C. Hunke, LANL subroutine coupling_prep (iblk) use ice_blocks, only: block, nx_block, ny_block use ice_calendar, only: dt, nstreams use ice_constants, only: c0, c1, puny, rhofresh use ice_domain_size, only: ncat use ice_flux, only: alvdf, alidf, alvdr, alidr, albice, albsno, & albpnd, albcnt, apeff_ai, coszen, fpond, fresh, l_mpond_fresh, & alvdf_ai, alidf_ai, alvdr_ai, alidr_ai, fhocn_ai, & fresh_ai, fsalt_ai, fsalt, & fswthru_ai, fhocn, fswthru, scale_factor, & swvdr, swidr, swvdf, swidf, Tf, Tair, Qa, strairxT, strairyt, & fsens, flat, fswabs, flwout, evap, Tref, Qref, faero_ocn, & fsurfn_f, flatn_f, scale_fluxes, frzmlt_init, frzmlt use ice_grid, only: tmask use ice_ocean, only: oceanmixed_ice, ocean_mixed_layer use ice_shortwave, only: alvdfn, alidfn, alvdrn, alidrn, & albicen, albsnon, albpndn, apeffn use ice_state, only: aicen, aice, aice_init, nbtrcr use ice_therm_shared, only: calc_Tsfc use ice_timers, only: timer_couple, ice_timer_start, ice_timer_stop use ice_zbgc_shared, only: flux_bio, flux_bio_ai integer (kind=int_kind), intent(in) :: & iblk ! block index ! local variables type (block) :: & this_block ! block information for current block integer (kind=int_kind) :: & n , & ! thickness category index i,j , & ! horizontal indices k , & ! tracer index ilo,ihi,jlo,jhi ! beginning and end of physical domain real (kind=dbl_kind) :: cszn ! counter for history averaging !----------------------------------------------------------------- ! Save current value of frzmlt for diagnostics. ! Update mixed layer with heat and radiation from ice. !----------------------------------------------------------------- do j = 1, ny_block do i = 1, nx_block frzmlt_init (i,j,iblk) = frzmlt(i,j,iblk) enddo enddo call ice_timer_start(timer_couple) ! atm/ocn coupling if (oceanmixed_ice) & call ocean_mixed_layer (dt,iblk) ! ocean surface fluxes and sst !----------------------------------------------------------------- ! Aggregate albedos !----------------------------------------------------------------- do j = 1, ny_block do i = 1, nx_block alvdf(i,j,iblk) = c0 alidf(i,j,iblk) = c0 alvdr(i,j,iblk) = c0 alidr(i,j,iblk) = c0 albice(i,j,iblk) = c0 albsno(i,j,iblk) = c0 albpnd(i,j,iblk) = c0 apeff_ai(i,j,iblk) = c0 ! for history averaging cszn = c0 if (coszen(i,j,iblk) > puny) cszn = c1 do n = 1, nstreams albcnt(i,j,iblk,n) = albcnt(i,j,iblk,n) + cszn enddo enddo enddo do n = 1, ncat do j = 1, ny_block do i = 1, nx_block alvdf(i,j,iblk) = alvdf(i,j,iblk) & + alvdfn(i,j,n,iblk)*aicen(i,j,n,iblk) alidf(i,j,iblk) = alidf(i,j,iblk) & + alidfn(i,j,n,iblk)*aicen(i,j,n,iblk) alvdr(i,j,iblk) = alvdr(i,j,iblk) & + alvdrn(i,j,n,iblk)*aicen(i,j,n,iblk) alidr(i,j,iblk) = alidr(i,j,iblk) & + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk) if (coszen(i,j,iblk) > puny) then ! sun above horizon albice(i,j,iblk) = albice(i,j,iblk) & + albicen(i,j,n,iblk)*aicen(i,j,n,iblk) albsno(i,j,iblk) = albsno(i,j,iblk) & + albsnon(i,j,n,iblk)*aicen(i,j,n,iblk) albpnd(i,j,iblk) = albpnd(i,j,iblk) & + albpndn(i,j,n,iblk)*aicen(i,j,n,iblk) endif apeff_ai(i,j,iblk) = apeff_ai(i,j,iblk) & ! for history + apeffn(i,j,n,iblk)*aicen(i,j,n,iblk) enddo enddo enddo do j = 1, ny_block do i = 1, nx_block !----------------------------------------------------------------- ! reduce fresh by fpond for coupling !----------------------------------------------------------------- if (l_mpond_fresh) then fpond(i,j,iblk) = fpond(i,j,iblk) * rhofresh/dt fresh(i,j,iblk) = fresh(i,j,iblk) - fpond(i,j,iblk) endif !---------------------------------------------------------------- ! Store grid box mean albedos and fluxes before scaling by aice !---------------------------------------------------------------- alvdf_ai (i,j,iblk) = alvdf (i,j,iblk) alidf_ai (i,j,iblk) = alidf (i,j,iblk) alvdr_ai (i,j,iblk) = alvdr (i,j,iblk) alidr_ai (i,j,iblk) = alidr (i,j,iblk) fresh_ai (i,j,iblk) = fresh (i,j,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) if (nbtrcr > 0) then do k = 1, nbtrcr flux_bio_ai (i,j,k,iblk) = flux_bio (i,j,k,iblk) enddo endif !----------------------------------------------------------------- ! Save net shortwave for scaling factor in scale_factor !----------------------------------------------------------------- scale_factor(i,j,iblk) = & swvdr(i,j,iblk)*(c1 - alvdr_ai(i,j,iblk)) & + swvdf(i,j,iblk)*(c1 - alvdf_ai(i,j,iblk)) & + swidr(i,j,iblk)*(c1 - alidr_ai(i,j,iblk)) & + swidf(i,j,iblk)*(c1 - alidf_ai(i,j,iblk)) enddo enddo !----------------------------------------------------------------- ! Divide fluxes by ice area ! - the CCSM coupler assumes fluxes are per unit ice area ! - also needed for global budget in diagnostics !----------------------------------------------------------------- call scale_fluxes (nx_block, ny_block, & tmask (:,:,iblk), nbtrcr, & aice (:,:,iblk), Tf (:,:,iblk), & Tair (:,:,iblk), Qa (:,:,iblk), & strairxT (:,:,iblk), strairyT(:,:,iblk), & fsens (:,:,iblk), flat (:,:,iblk), & fswabs (:,:,iblk), flwout (:,:,iblk), & evap (:,:,iblk), & Tref (:,:,iblk), Qref (:,:,iblk), & fresh (:,:,iblk), fsalt (:,:,iblk), & fhocn (:,:,iblk), fswthru (:,:,iblk), & faero_ocn(:,:,:,iblk), & alvdr (:,:,iblk), alidr (:,:,iblk), & alvdf (:,:,iblk), alidf (:,:,iblk), & flux_bio(:,:,:,iblk)) !echmod - comment this out for efficiency, if .not. calc_Tsfc if (.not. calc_Tsfc) then !--------------------------------------------------------------- ! If surface fluxes were provided, conserve these fluxes at ice ! free points by passing to ocean. !--------------------------------------------------------------- call sfcflux_to_ocn & (nx_block, ny_block, & tmask (:,:,iblk), aice_init(:,:,iblk), & fsurfn_f (:,:,:,iblk), flatn_f(:,:,:,iblk), & fresh (:,:,iblk), fhocn (:,:,iblk)) endif !echmod call ice_timer_stop(timer_couple) ! atm/ocn coupling 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 ! 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, ! conserve energy and water by passing these fluxes to the ocean. ! ! author: A. McLaren, Met Office subroutine sfcflux_to_ocn(nx_block, ny_block, & tmask, aice, & fsurfn_f, flatn_f, & fresh, fhocn) use ice_constants, only: c1, puny, lsub use ice_domain_size, only: ncat integer (kind=int_kind), intent(in) :: & nx_block, ny_block ! block dimensions logical (kind=log_kind), dimension (nx_block,ny_block), & intent(in) :: & tmask ! land/boundary mask, thickness (T-cell) real (kind=dbl_kind), dimension(nx_block,ny_block), & intent(in):: & aice ! initial ice concentration real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), & intent(in) :: & fsurfn_f, & ! net surface heat flux (provided as forcing) flatn_f ! latent heat flux (provided as forcing) real (kind=dbl_kind), dimension(nx_block,ny_block), & intent(inout):: & fresh , & ! fresh water flux to ocean (kg/m2/s) fhocn ! actual ocn/ice heat flx (W/m**2) #ifdef CICE_IN_NEMO ! local variables integer (kind=int_kind) :: & i, j, n ! horizontal indices real (kind=dbl_kind) :: & rLsub ! 1/Lsub rLsub = c1 / Lsub do n = 1, ncat do j = 1, ny_block do i = 1, nx_block if (tmask(i,j) .and. aice(i,j) <= puny) then fhocn(i,j) = fhocn(i,j) & + fsurfn_f(i,j,n) + flatn_f(i,j,n) fresh(i,j) = fresh(i,j) & + flatn_f(i,j,n) * rLsub endif enddo ! i enddo ! j enddo ! n #endif end subroutine sfcflux_to_ocn !======================================================================= ! ! Update the ice state variables using the ice concentration increment rate ! calculated in the NEMO data assimilation (DA) scheme. ! Ice area is added by adding ADDITIONAL ice with thickness hi_da_new. ! This implies the ADDITIONAL volume added is hi_da_new*daice, where ! daice is the change in ice area due to DA. ! Ice area is subtracted by removing ice area with the current category ! thickness. Ice area is first removed from the lowest category, and then ! removed from higher categories as needed. ! ! authors: D. Peterson, Met Office ! A. McLaren, Met Office subroutine da_state_update use ice_constants, only: c1, puny #ifdef ICE_DA integer (kind=int_kind) :: & i, j, ij , & ! horizontal indices iblk , & ! block index ilo,ihi,jlo,jhi, & ! beginning and end of physical domain n ! thickness category index integer (kind=int_kind) :: & nelevate ! number of elevations of increments to higher ! category (diagnostic) integer (kind=int_kind) :: & icells ! number of ocean cells integer (kind=int_kind), dimension(nx_block*ny_block) :: & indxi, indxj ! indirect indices for cells with aicen > puny type (block) :: & this_block ! block information for current block real (kind=dbl_kind) :: & hi_da_new , & ! specified ice thickness for new ice created by DA hicen , & ! ice thickness hsnon , & ! snow thickness daice , & ! change in ice concentration (for first category) dvice , & ! change in ice volume dvsno ! change in snow volume real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: & vsno_init , & ! initial snow volume vice_init ! initial ice volume !---------------------------------------------------------------- ! This routine will only work under certain circumstances!! ! Note, if any optional tracers are used in the run, they will not ! be conserved here. !---------------------------------------------------------------- if (nilyr /= 1 .or. nslyr /= 1 .or. ntrcr /= 1) & call abort_ice & ('da_state_update: only works for 1 cat, 1 layer, 1 tracer runs') !------------------------------------------------------------------ ! Set thickness for new ice ! (Currently using value of 0.5m, which was value thin ice was ! incremented toward in LIM ice model). !----------------------------------------------------------------- hi_da_new = 0.50_dbl_kind ! if ncat>1, this has to be less than ! the 1st category thickness limit ! Initialise various fields vsno_init(:,:,:) = c0 vice_init(:,:,:) = c0 fresh_da(:,:,:) = c0 fsalt_da(:,:,:) = c0 !---------------------------------------------------------------- ! Update category state variables !---------------------------------------------------------------- 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 !---------------------------------------------------------------- ! Find ocean points where data assimilation abs(increment) > puny ! (Note, daice_da is the RATE of change of ice concentration due ! to DA) !---------------------------------------------------------------- icells = 0 do j = jlo, jhi do i = ilo, ihi if (tmask(i,j,iblk) .and. abs(daice_da(i,j,iblk)*dt) > puny) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif vsno_init(i,j,iblk) = vsno(i,j,iblk) ! used for salinity changes vice_init(i,j,iblk) = vice(i,j,iblk) ! used for salinity changes enddo ! i enddo ! j if (icells > 0) then n = 1 ! only ever add increment to 1st category do ij = 1, icells i = indxi(ij) j = indxj(ij) !--------------------------------------------------- ! Apply concentration increment and associated ! volume change !--------------------------------------------------- if (aicen(i,j,n,iblk) > puny) then ! if decreasing concentration, subtract ice volume at ! current thickness hicen = vicen(i,j,n,iblk) / aicen(i,j,n,iblk) ! if increasing concentration, add ice volume at hi_da_new ! thickness if ( daice_da(i,j,iblk)*dt > puny) hicen = hi_da_new ! whether in/decreasing concentration, add/subtract snow ! volume at current thickness hsnon = vsnon(i,j,n,iblk) / aicen(i,j,n,iblk) daice = & min( ( c1 - aice(i,j,iblk) ), ( daice_da(i,j,iblk)*dt ) ) aicen(i,j,n,iblk) = aicen(i,j,n,iblk) + daice vicen(i,j,n,iblk) = vicen(i,j,n,iblk) + hicen*daice vsnon(i,j,n,iblk) = aicen(i,j,n,iblk) * hsnon !--------------------------------------------------- ! Create new ice points with specified thickness !--------------------------------------------------- else aicen(i,j,n,iblk) = & min( ( c1 - aice(i,j,iblk) ), ( daice_da(i,j,iblk)*dt ) ) ! note aicen/vicen < c0 taken care below vicen(i,j,n,iblk) = aicen(i,j,n,iblk) * hi_da_new vsnon(i,j,n,iblk) = c0 endif enddo ! ij do n = 1,ncat nelevate=0 do ij = 1,icells i = indxi(ij) j = indxj(ij) !---------------------------------------------------- ! Check is aicen < puny ! - remove from next category if necessary ! - otherwise just remove it ! Ignoring conservation issues here !---------------------------------------------------- if (aicen(i,j,n,iblk) < puny) then if ( n < ncat ) then if (aicen(i,j,n,iblk) < -1.0*puny ) then nelevate=nelevate+1 endif ! take concentration from next category -- constant thickness if ( aicen(i,j,n+1,iblk) > puny ) then hicen = vicen(i,j,n+1,iblk)/aicen(i,j,n+1,iblk) hsnon = vsnon(i,j,n+1,iblk)/aicen(i,j,n+1,iblk) else hicen = c0 hsnon = c0 endif ! aicen(n+1) > puny aicen(i,j,n+1,iblk) = aicen(i,j,n+1,iblk) + aicen(i,j,n,iblk) vicen(i,j,n+1,iblk) = aicen(i,j,n+1,iblk) * hicen vsnon(i,j,n+1,iblk) = aicen(i,j,n+1,iblk) * hsnon endif ! n < ncat aicen(i,j,n,iblk) = c0 vicen(i,j,n,iblk) = c0 vsnon(i,j,n,iblk) = c0 eicen(i,j,n,iblk) = c0 esnon(i,j,n,iblk) = c0 endif ! aicen(n) < puny !--------------------------------------------------- ! Update energies !--------------------------------------------------- ! Would need vertical layers here in the future if (aicen(i,j,n,iblk) > puny) then esnon(i,j,n,iblk) = -rhos*Lfresh*vsnon(i,j,n,iblk) eicen(i,j,n,iblk) = -rhoi*Lfresh*vicen(i,j,n,iblk) endif enddo ! ij !write(nu_diag,*) 'Elevated ', nelevate, ' incs to category ', n+1 enddo ! n endif ! icells enddo ! nblocks !------------------------------------------------------------------- ! Ghost cell updates for state variables. (Can't be called within ! block do loop). !------------------------------------------------------------------- call bound_state (aicen, trcrn, vicen, vsnon, eicen, esnon) 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 !----------------------------------------------------------- ! Find data assimilation points again !----------------------------------------------------------- icells = 0 do j = jlo, jhi do i = ilo, ihi if (tmask(i,j,iblk) .and. abs(daice_da(i,j,iblk)*dt) > puny) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif enddo ! i enddo ! j !------------------------------------------------------------- ! Update aggregate values !------------------------------------------------------------- call aggregate (nx_block, ny_block, & aicen(:,:,:,iblk), & trcrn(:,:,1:ntrcr,:,iblk), & vicen(:,:,:,iblk), vsnon(:,:, :,iblk), & eicen(:,:,:,iblk), esnon(:,:, :,iblk), & aice (:,:, iblk), & trcr (:,:,1:ntrcr, iblk), & vice (:,:, iblk), vsno (:,:, iblk), & eice (:,:, iblk), esno (:,:, iblk), & aice0(:,:, iblk), tmask(:,:, iblk), & ntrcr , trcr_depend(1:ntrcr) ) !------------------------------------------------------------- ! Calculate implied freshwater and salt fluxes !------------------------------------------------------------- if (icells > 0) then do ij = 1, icells i = indxi(ij) j = indxj(ij) dvice = vice(i,j,iblk) - vice_init(i,j,iblk) dvsno = vsno(i,j,iblk) - vsno_init(i,j,iblk) fresh_da(i,j,iblk) = - (rhoi * dvice + rhos * dvsno)/dt fsalt_da(i,j,iblk) = - rhoi*ice_ref_salinity*p001*dvice/dt enddo ! ij endif ! icells enddo ! iblk #endif end subroutine da_state_update !======================================================================= end module CICE_RunMod !=======================================================================