! SVN:$Id: ice_mechred.F90 1304 2015-11-12 19:28:14Z afrobert@nps.edu $ !======================================================================= ! Driver for ice mechanical redistribution (ridging) ! ! See these references: ! ! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength ! in modeling the thickness distribution of Arctic sea ice, ! J. Geophys. Res., 100, 18,611-18,626. ! ! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice ! cover, Mon. Wea. Rev., 108, 1943-1973, 1980. ! ! Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki, 2007: ! Improving ridging schemes for high-resolution sea ice models. ! J. Geophys. Res. 112, C03S91, doi:10.1029/2005JC003355. ! ! Rothrock, D. A., 1975: The energetics of the plastic deformation of ! pack ice by ridging, J. Geophys. Res., 80, 4514-4519. ! ! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, ! 1975: The thickness distribution of sea ice, J. Geophys. Res., ! 80, 4501-4513. ! ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL ! ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb ! 2004: Block structure added by William Lipscomb ! 2006: New options for participation and redistribution (WHL) ! 2006: Streamlined for efficiency by Elizabeth Hunke ! Converted to free source form (F90) module ice_mechred use ice_kinds_mod use ice_constants use ice_domain_size, only: ncat, max_aero, n_aero, nilyr, nslyr, nblyr use ice_fileunits, only: nu_diag use ice_itd, only: hin_max, column_sum, & column_conservation_check, compute_tracers implicit none save private public :: ice_strength, ridge_ice !----------------------------------------------------------------------- ! Ridging parameters !----------------------------------------------------------------------- integer (kind=int_kind), public :: & ! defined in namelist kstrength , & ! 0 for simple Hibler (1979) formulation ! 1 for Rothrock (1975) pressure formulation krdg_partic , & ! 0 for Thorndike et al. (1975) formulation ! 1 for exponential participation function krdg_redist ! 0 for Hibler (1980) formulation ! 1 for exponential redistribution function real (kind=dbl_kind), public :: & mu_rdg, & ! gives e-folding scale of ridged ice (m^.5) (krdg_redist = 1) Cf ! ratio of ridging work to PE change in ridging (kstrength = 1) real (kind=dbl_kind), parameter :: & Cs = p25 , & ! fraction of shear energy contrbtng to ridging Cp = p5*gravit*(rhow-rhoi)*rhoi/rhow, & ! proport const for PE fsnowrdg = p5 , & ! snow fraction that survives in ridging Gstar = p15 , & ! max value of G(h) that participates ! (krdg_partic = 0) astar = p05 , & ! e-folding scale for G(h) participation !echmod astar = p1 , & ! e-folding scale for G(h) participation ! (krdg_partic = 1) maxraft= c1 , & ! max value of hrmin - hi = max thickness ! of ice that rafts (m) Hstar = c25 , & ! determines mean thickness of ridged ice (m) ! (krdg_redist = 0) ! Flato & Hibler (1995) have Hstar = 100 Pstar = 2.75e4_dbl_kind, & ! constant in Hibler strength formula ! (kstrength = 0) Cstar = c20 ! constant in Hibler strength formula ! (kstrength = 0) logical (kind=log_kind), parameter :: & ! l_conservation_check = .true. ! if true, check conservation l_conservation_check = .false. ! if true, check conservation ! (useful for debugging) !======================================================================= contains !======================================================================= ! Compute changes in the ice thickness distribution due to divergence ! and shear. ! ! author: William H. Lipscomb, LANL subroutine ridge_ice (nx_block, ny_block, & dt, ndtd, & ntrcr, icells, & indxi, indxj, & rdg_conv, rdg_shear, & aicen, trcrn, & vicen, vsnon, & aice0, & trcr_depend, l_stop, & istop, jstop, & dardg1dt, dardg2dt, & dvirdgdt, opening, & fpond, & fresh, fhocn, & faero_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & dvirdgndt, & araftn, vraftn) use ice_state, only: nt_qice, nt_qsno, tr_brine, nt_fbri, nt_sice integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells , & ! number of cells with ice present ndtd , & ! number of dynamics subcycles ntrcr ! number of tracers in use integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & rdg_conv, & ! normalized energy dissipation due to convergence (1/s) rdg_shear ! normalized energy dissipation due to shear (1/s) real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(inout) :: & aicen , & ! concentration of ice vicen , & ! volume per unit area of ice (m) vsnon ! volume per unit area of snow (m) real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & aice0 ! concentration of open water integer (kind=int_kind), dimension(ntrcr), intent(in) :: & trcr_depend logical (kind=log_kind), intent(out) :: & l_stop ! if true, abort on return integer (kind=int_kind), intent(out) :: & istop, jstop ! indices of grid cell where model aborts ! optional history fields real (kind=dbl_kind), dimension(nx_block,ny_block), & intent(inout), optional :: & dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) dardg2dt , & ! rate of fractional area gain by new ridges (1/s) dvirdgdt , & ! rate of ice volume ridged (m/s) opening , & ! rate of opening due to divergence/shear (1/s) fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fhocn ! net heat flux to ocean (W/m^2) real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), & intent(out), optional :: & dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) dvirdgndt , & ! rate of ice volume ridged (m/s) aparticn , & ! participation function krdgn , & ! mean ridge thickness/thickness of ridging ice araftn , & ! rafting ice area vraftn , & ! rafting ice volume aredistn , & ! redistribution function: fraction of new ridge area vredistn ! redistribution function: fraction of new ridge volume real (kind=dbl_kind), dimension(nx_block,ny_block,max_aero), & intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) ! local variables real (kind=dbl_kind), dimension (nx_block,ny_block,ncat) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ncat) :: & esnon , & ! energy of melting for each snow layer (J/m^2) vbrin, & ! ice volume with defined by brine height (m) sicen ! Bulk salt in h ice (ppt*m) real (kind=dbl_kind), dimension (icells) :: & asum , & ! sum of ice and open water area aksum , & ! ratio of area removed to area ridged msnow_mlt , & ! mass of snow added to ocean (kg m-2) esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) mpond , & ! mass of pond added to ocean (kg m-2) closing_net, & ! net rate at which area is removed (1/s) ! (ridging ice area - area of new ridges) / dt divu_adv , & ! divu as implied by transport scheme (1/s) opning , & ! rate of opening due to divergence/shear ! opning is a local variable; ! opening is the history diagnostic variable ardg1 , & ! fractional area loss by ridging ice ardg2 , & ! fractional area gain by new ridges virdg , & ! ice volume ridged aopen ! area opening due to divergence/shear real (kind=dbl_kind), dimension (icells,max_aero) :: & maero ! aerosol mass added to ocean (kg m-2) real (kind=dbl_kind), dimension (icells,0:ncat) :: & apartic ! participation function; fraction of ridging ! and closing associated w/ category n real (kind=dbl_kind), dimension (icells,ncat) :: & hrmin , & ! minimum ridge thickness hrmax , & ! maximum ridge thickness (krdg_redist = 0) hrexp , & ! ridge e-folding thickness (krdg_redist = 1) krdg , & ! mean ridge thickness/thickness of ridging ice ardg1n , & ! area of ice ridged ardg2n , & ! area of new ridges virdgn , & ! ridging ice volume mraftn ! rafting ice mask real (kind=dbl_kind), dimension (icells) :: & vice_init, vice_final, & ! ice volume summed over categories vsno_init, vsno_final, & ! snow volume summed over categories eice_init, eice_final, & ! ice energy summed over layers vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories Shi_init , Shi_final , & ! ice bulk salinity summed over categories esno_init, esno_final ! snow energy summed over layers integer (kind=int_kind), parameter :: & nitermax = 20 ! max number of ridging iterations integer (kind=int_kind) :: & i,j , & ! horizontal indices n , & ! thickness category index niter , & ! iteration counter ij , & ! horizontal index, combines i and j loops k ! vertical index real (kind=dbl_kind) :: & dti ! 1 / dt logical (kind=log_kind) :: & iterate_ridging ! if true, repeat the ridging character (len=char_len) :: & fieldid ! field identifier !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- l_stop = .false. istop = 0 jstop = 0 do ij = 1, icells msnow_mlt(ij) = c0 esnow_mlt(ij) = c0 maero (ij,:) = c0 mpond (ij) = c0 ardg1 (ij) = c0 ardg2 (ij) = c0 virdg (ij) = c0 ardg1n (ij,:) = c0 ardg2n (ij,:) = c0 virdgn (ij,:) = c0 mraftn (ij,:) = c0 ! aopen (ij) = c0 enddo !----------------------------------------------------------------- ! Compute area of ice plus open water before ridging. !----------------------------------------------------------------- call asum_ridging (nx_block, ny_block, & icells, indxi, indxj, & aicen, aice0, & asum) !----------------------------------------------------------------- ! Compute the area opening and closing. !----------------------------------------------------------------- call ridge_prep (nx_block, ny_block, & icells, indxi, indxj, & dt, & rdg_conv, rdg_shear, & asum, closing_net, & divu_adv, opning) !----------------------------------------------------------------- ! Compute initial values of conserved quantities. !----------------------------------------------------------------- if (l_conservation_check) then eicen(:,:,:) = c0 esnon(:,:,:) = c0 vbrin(:,:,:) = c0 sicen(:,:,:) = c0 do n = 1, ncat do j = 1, ny_block do i = 1, nx_block vbrin(i,j,n) = vicen(i,j,n) if (tr_brine) vbrin(i,j,n) = trcrn(i,j,nt_fbri,n) * vicen(i,j,n) enddo enddo do k = 1, nilyr do j = 1, ny_block do i = 1, nx_block eicen(i,j,n) = eicen(i,j,n) + trcrn(i,j,nt_qice+k-1,n) & * vicen(i,j,n)/real(nilyr,kind=dbl_kind) enddo enddo enddo do k = 1, nslyr do j = 1, ny_block do i = 1, nx_block esnon(i,j,n) = esnon(i,j,n) + trcrn(i,j,nt_qsno+k-1,n) & * vsnon(i,j,n)/real(nslyr,kind=dbl_kind) enddo enddo enddo do k = 1, nilyr do j = 1, ny_block do i = 1, nx_block sicen(i,j,n) = sicen(i,j,n) + trcrn(i,j,nt_sice+k-1,n) & * vicen(i,j,n)/real(nilyr,kind=dbl_kind) enddo enddo enddo enddo call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vicen, vice_init) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vsnon, vsno_init) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & eicen, eice_init) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & esnon, esno_init) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vbrin, vbri_init) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & sicen, Shi_init) endif do niter = 1, nitermax !----------------------------------------------------------------- ! Compute the thickness distribution of ridging ice ! and various quantities associated with the new ridged ice. !----------------------------------------------------------------- call ridge_itd (nx_block, ny_block, & icells, indxi, indxj, & aicen, vicen, & aice0, & aksum, apartic, & hrmin, hrmax, & hrexp, krdg, & aparticn, krdgn, mraftn) !----------------------------------------------------------------- ! Redistribute area, volume, and energy. !----------------------------------------------------------------- call ridge_shift (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, dt, & aicen, trcrn, & vicen, vsnon, & aice0, trcr_depend, & aksum, apartic, & hrmin, hrmax, & hrexp, krdg, & closing_net, opning, & ardg1, ardg2, & virdg, aopen, & ardg1n, ardg2n, & virdgn, & msnow_mlt, esnow_mlt, & maero, mpond, & l_stop, & istop, jstop, & aredistn, vredistn) if (l_stop) return !----------------------------------------------------------------- ! Make sure the new area = 1. If not (because the closing ! and opening rates were reduced above), prepare to ridge again ! with new rates. !----------------------------------------------------------------- call asum_ridging (nx_block, ny_block, & icells, indxi, indxj, & aicen, aice0, & asum) call ridge_check (nx_block, ny_block, & icells, indxi, indxj, & dt, & asum, closing_net, & divu_adv, opning, & iterate_ridging) !----------------------------------------------------------------- ! If done, exit. If not, prepare to ridge again. !----------------------------------------------------------------- if (iterate_ridging) then write(nu_diag,*) 'Repeat ridging, niter =', niter else exit endif if (niter == nitermax) then write(nu_diag,*) ' ' write(nu_diag,*) 'Exceeded max number of ridging iterations' write(nu_diag,*) 'max =',nitermax l_stop = .true. return endif enddo ! niter !----------------------------------------------------------------- ! Compute final values of conserved quantities. ! Check for conservation (allowing for snow thrown into ocean). !----------------------------------------------------------------- if (l_conservation_check) then eicen(:,:,:) = c0 esnon(:,:,:) = c0 sicen(:,:,:) = c0 vbrin(:,:,:) = c0 do n = 1, ncat do k = 1, nilyr do j = 1, ny_block do i = 1, nx_block eicen(i,j,n) = eicen(i,j,n) + trcrn(i,j,nt_qice+k-1,n) & * vicen(i,j,n)/real(nilyr,kind=dbl_kind) enddo enddo enddo do k = 1, nslyr do j = 1, ny_block do i = 1, nx_block esnon(i,j,n) = esnon(i,j,n) + trcrn(i,j,nt_qsno+k-1,n) & * vsnon(i,j,n)/real(nslyr,kind=dbl_kind) enddo enddo enddo do j = 1, ny_block do i = 1, nx_block vbrin(i,j,n) = vicen(i,j,n) if (tr_brine) vbrin(i,j,n) = trcrn(i,j,nt_fbri,n) * vbrin(i,j,n) enddo enddo do k = 1, nilyr do j = 1, ny_block do i = 1, nx_block sicen(i,j,n) = sicen(i,j,n) + trcrn(i,j,nt_sice+k-1,n) & * vicen(i,j,n)/real(nilyr,kind=dbl_kind) enddo enddo enddo enddo call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vicen, vice_final) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vsnon, vsno_final) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & eicen, eice_final) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & esnon, esno_final) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & sicen, Shi_final) call column_sum (nx_block, ny_block, & icells, indxi, indxj, & ncat, & vbrin, vbri_final) do ij = 1, icells vsno_final(ij) = vsno_final(ij) + msnow_mlt(ij)/rhos esno_final(ij) = esno_final(ij) + esnow_mlt(ij) enddo fieldid = 'vice, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & vice_init, vice_final, & puny, l_stop, & istop, jstop) if (l_stop) return fieldid = 'vsno, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & vsno_init, vsno_final, & puny, l_stop, & istop, jstop) if (l_stop) return fieldid = 'eice, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & eice_init, eice_final, & puny*Lfresh*rhoi, & l_stop, & istop, jstop) if (l_stop) return fieldid = 'esno, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & esno_init, esno_final, & puny*Lfresh*rhos, & l_stop, & istop, jstop) if (l_stop) return fieldid = 'sice, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & Shi_init, Shi_final, & puny, l_stop, & istop, jstop) if (l_stop) return fieldid = 'vbrin, ridging' call column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & vbri_init, vbri_final, & puny*c10, l_stop, & istop, jstop) if (l_stop) return endif ! l_conservation_check !----------------------------------------------------------------- ! Compute ridging diagnostics. !----------------------------------------------------------------- dti = c1/dt if (present(dardg1dt)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) dardg1dt(i,j) = ardg1(ij)*dti enddo endif if (present(dardg2dt)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) dardg2dt(i,j) = ardg2(ij)*dti enddo endif if (present(dvirdgdt)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) dvirdgdt(i,j) = virdg(ij)*dti enddo endif if (present(opening)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) opening(i,j) = aopen(ij)*dti enddo endif if (present(dardg1ndt)) then do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) dardg1ndt(i,j,n) = ardg1n(ij,n)*dti enddo enddo endif if (present(dardg2ndt)) then do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) dardg2ndt(i,j,n) = ardg2n(ij,n)*dti enddo enddo endif if (present(dvirdgndt)) then do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) dvirdgndt(i,j,n) = virdgn(ij,n)*dti enddo enddo endif if (present(araftn)) then do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) araftn(i,j,n) = mraftn(ij,n)*ardg2n(ij,n) ! araftn(i,j,n) = mraftn(ij,n)*ardg1n(ij,n)*p5 enddo enddo endif if (present(vraftn)) then do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) vraftn(i,j,n) = mraftn(ij,n)*virdgn(ij,n) enddo enddo endif !----------------------------------------------------------------- ! Update fresh water and heat fluxes due to snow melt. !----------------------------------------------------------------- ! use thermodynamic time step (ndtd*dt here) to average properly dti = c1/(ndtd*dt) if (present(fresh)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) fresh(i,j) = fresh(i,j) + msnow_mlt(ij)*dti enddo endif if (present(fhocn)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) fhocn(i,j) = fhocn(i,j) + esnow_mlt(ij)*dti enddo endif if (present(faero_ocn)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) faero_ocn(i,j,:) = faero_ocn(i,j,:) + maero(ij,:)*dti enddo endif if (present(fpond)) then do ij = 1, icells i = indxi(ij) j = indxj(ij) fpond(i,j) = fpond(i,j) - mpond(ij) ! units change later enddo endif !----------------------------------------------------------------- ! Check for fractional ice area > 1. !----------------------------------------------------------------- do ij = 1, icells i = indxi(ij) j = indxj(ij) if (abs(asum(ij) - c1) > puny) then l_stop = .true. istop = i jstop = j write(nu_diag,*) ' ' write(nu_diag,*) 'Ridging error: total area > 1' write(nu_diag,*) 'i, j, area:', i, j, asum(ij) write(nu_diag,*) 'n, aicen:' write(nu_diag,*) 0, aice0(i,j) do n = 1, ncat write(nu_diag,*) n, aicen(i,j,n) enddo return endif enddo end subroutine ridge_ice !======================================================================= ! Find the total area of ice plus open water in each grid cell. ! ! This is similar to the aggregate_area subroutine except that the ! total area can be greater than 1, so the open water area is ! included in the sum instead of being computed as a residual. ! ! author: William H. Lipscomb, LANL subroutine asum_ridging (nx_block, ny_block, & icells, indxi, indxj, & aicen, aice0, & asum) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells with ice present integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(in) :: & aicen ! concentration of ice in each category real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & aice0 ! concentration of open water real (kind=dbl_kind), dimension (icells), intent(out):: & asum ! sum of ice and open water area ! local variables integer (kind=int_kind) :: & i, j, n, & ij ! horizontal index, combines i and j loops !----------------------------------------------------------------- ! open water !----------------------------------------------------------------- do ij = 1, icells i = indxi(ij) j = indxj(ij) asum(ij) = aice0(i,j) enddo !----------------------------------------------------------------- ! ice categories !----------------------------------------------------------------- do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) asum(ij) = asum(ij) + aicen(i,j,n) enddo enddo end subroutine asum_ridging !======================================================================= ! Initialize arrays, compute area of closing and opening ! ! author: William H. Lipscomb, LANL subroutine ridge_prep (nx_block, ny_block, & icells, indxi, indxj, & dt, & rdg_conv, rdg_shear, & asum, closing_net, & divu_adv, opning) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells with ice present integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), intent(in) :: & dt ! time step (s) real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & rdg_conv, & ! normalized energy dissipation due to convergence (1/s) rdg_shear ! normalized energy dissipation due to shear (1/s) real (kind=dbl_kind), dimension(icells), & intent(inout):: & asum ! sum of ice and open water area real (kind=dbl_kind), dimension(icells), & intent(out):: & closing_net, & ! net rate at which area is removed (1/s) divu_adv , & ! divu as implied by transport scheme (1/s) opning ! rate of opening due to divergence/shear ! local variables real (kind=dbl_kind), parameter :: & big = 1.0e+8_dbl_kind integer (kind=int_kind) :: & i,j, & ! horizontal indices ij ! horizontal index, combines i and j loops ! Set hin_max(ncat) to a big value to ensure that all ridged ice ! is thinner than hin_max(ncat). hin_max(ncat) = big do ij = 1, icells i = indxi(ij) j = indxj(ij) !----------------------------------------------------------------- ! Compute the net rate of closing due to convergence ! and shear, based on Flato and Hibler (1995). ! ! For the elliptical yield curve: ! rdg_conv = -min (divu, 0) ! rdg_shear = (1/2) * (Delta - abs(divu)) ! Note that the shear term also accounts for divergence. ! ! The energy dissipation rate is equal to the net closing rate ! times the ice strength. ! ! NOTE: The NET closing rate is equal to the rate that open water ! area is removed, plus the rate at which ice area is removed by ! ridging, minus the rate at which area is added in new ridges. ! The GROSS closing rate is equal to the first two terms (open ! water closing and thin ice ridging) without the third term ! (thick, newly ridged ice). ! ! rdg_conv is calculated differently in EAP (update_ice_rdg) and ! represents closing_net directly. In that case, rdg_shear=0. !----------------------------------------------------------------- closing_net(ij) = Cs*rdg_shear(i,j) + rdg_conv(i,j) !----------------------------------------------------------------- ! Compute divu_adv, the divergence rate given by the transport/ ! advection scheme, which may not be equal to divu as computed ! from the velocity field. ! ! If divu_adv < 0, make sure the closing rate is large enough ! to give asum = 1.0 after ridging. !----------------------------------------------------------------- divu_adv(ij) = (c1-asum(ij)) / dt if (divu_adv(ij) < c0) & closing_net(ij) = max(closing_net(ij), -divu_adv(ij)) !----------------------------------------------------------------- ! Compute the (non-negative) opening rate that will give ! asum = 1.0 after ridging. !----------------------------------------------------------------- opning(ij) = closing_net(ij) + divu_adv(ij) enddo end subroutine ridge_prep !======================================================================= ! Compute the thickness distribution of the ice and open water ! participating in ridging and of the resulting ridges. ! ! This version includes new options for ridging participation and ! redistribution. ! The new participation scheme (krdg_partic = 1) improves stability ! by increasing the time scale for large changes in ice strength. ! The new exponential redistribution function (krdg_redist = 1) improves ! agreement between ITDs of modeled and observed ridges. ! ! author: William H. Lipscomb, LANL ! ! 2006: Changed subroutine name to ridge_itd ! Added new options for ridging participation and redistribution. subroutine ridge_itd (nx_block, ny_block, & icells, indxi, indxj, & aicen, vicen, & aice0, & aksum, apartic, & hrmin, hrmax, & hrexp, krdg, & aparticn, krdgn, mraft) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells with ice present integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(in) :: & aicen , & ! concentration of ice vicen ! volume per unit area of ice (m) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & aice0 ! concentration of open water real (kind=dbl_kind), dimension (icells), intent(out):: & aksum ! ratio of area removed to area ridged real (kind=dbl_kind), dimension (icells,0:ncat), & intent(out) :: & apartic ! participation function; fraction of ridging ! and closing associated w/ category n real (kind=dbl_kind), dimension (icells,ncat), & intent(out) :: & hrmin , & ! minimum ridge thickness hrmax , & ! maximum ridge thickness (krdg_redist = 0) hrexp , & ! ridge e-folding thickness (krdg_redist = 1) krdg ! mean ridge thickness/thickness of ridging ice ! diagnostic, category values real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), & intent(out), optional :: & aparticn , & ! participation function krdgn ! mean ridge thickness/thickness of ridging ice real (kind=dbl_kind), dimension (icells,ncat), & intent(out), optional :: & mraft ! rafting ice mask ! local variables integer (kind=int_kind) :: & i,j , & ! horizontal indices n , & ! thickness category index ij ! horizontal index, combines i and j loops real (kind=dbl_kind), parameter :: & Gstari = c1/Gstar, & astari = c1/astar real (kind=dbl_kind), dimension(icells,-1:ncat) :: & Gsum ! Gsum(n) = sum of areas in categories 0 to n real (kind=dbl_kind), dimension(icells) :: & work ! temporary work array real (kind=dbl_kind) :: & hi , & ! ice thickness for each cat (m) hrmean , & ! mean ridge thickness (m) xtmp ! temporary variable !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- do ij = 1, icells Gsum (ij,-1) = c0 ! by definition Gsum (ij,0) = c1 ! to avoid divzero below apartic(ij,0) = c0 enddo do n = 1, ncat do ij = 1, icells Gsum (ij,n) = c1 ! to avoid divzero below apartic(ij,n) = c0 hrmin (ij,n) = c0 hrmax (ij,n) = c0 hrexp (ij,n) = c0 krdg (ij,n) = c1 enddo enddo !----------------------------------------------------------------- ! Compute the thickness distribution of ice participating in ridging. !----------------------------------------------------------------- !----------------------------------------------------------------- ! First compute the cumulative thickness distribution function Gsum, ! where Gsum(n) is the fractional area in categories 0 to n. ! Ignore categories with very small areas. !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aice0(i,j) > puny) then Gsum(ij,0) = aice0(i,j) else Gsum(ij,0) = Gsum(ij,-1) endif enddo do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny) then Gsum(ij,n) = Gsum(ij,n-1) + aicen(i,j,n) else Gsum(ij,n) = Gsum(ij,n-1) endif enddo enddo ! normalize do ij = 1, icells work(ij) = c1 / Gsum(ij,ncat) enddo do n = 0, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells Gsum(ij,n) = Gsum(ij,n) * work(ij) enddo enddo !----------------------------------------------------------------- ! Compute the participation function apartic; this is analogous to ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). ! ! area lost from category n due to ridging/closing ! apartic(n) = -------------------------------------------------- ! total area lost due to ridging/closing ! !----------------------------------------------------------------- if (krdg_partic == 0) then ! Thornike et al. 1975 formulation !----------------------------------------------------------------- ! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). ! The expressions for apartic are found by integrating b(h)g(h) between ! the category boundaries. !----------------------------------------------------------------- do n = 0, ncat do ij = 1, icells if (Gsum(ij,n) < Gstar) then apartic(ij,n) = Gstari*(Gsum(ij,n)-Gsum(ij,n-1)) * & (c2 - (Gsum(ij,n-1)+Gsum(ij,n))*Gstari) elseif (Gsum(ij,n-1) < Gstar) then apartic(ij,n) = Gstari * (Gstar-Gsum(ij,n-1)) * & (c2 - (Gsum(ij,n-1)+Gstar)*Gstari) endif enddo ! ij enddo ! n elseif (krdg_partic==1) then ! exponential dependence on G(h) !----------------------------------------------------------------- ! b(h) = exp(-G(h)/astar) ! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)]. ! The expression for apartic is found by integrating b(h)g(h) ! between the category boundaries. !----------------------------------------------------------------- ! precompute exponential terms using Gsum as work array xtmp = c1 / (c1 - exp(-astari)) do n = -1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells Gsum(ij,n) = exp(-Gsum(ij,n)*astari) * xtmp enddo ! ij enddo ! n do n = 0, ncat do ij = 1, icells apartic(ij,n) = Gsum(ij,n-1) - Gsum(ij,n) enddo ! ij enddo ! n endif ! krdg_partic !----------------------------------------------------------------- ! Compute variables related to ITD of ridged ice: ! ! krdg = mean ridge thickness / thickness of ridging ice ! hrmin = min ridge thickness ! hrmax = max ridge thickness (krdg_redist = 0) ! hrexp = ridge e-folding scale (krdg_redist = 1) !---------------------------------------------------------------- if (krdg_redist == 0) then ! Hibler 1980 formulation !----------------------------------------------------------------- ! Assume ridged ice is uniformly distributed between hrmin and hrmax. ! ! This parameterization is a modified version of Hibler (1980). ! In the original paper the min ridging thickness is hrmin = 2*hi, ! and the max thickness is hrmax = 2*sqrt(hi*Hstar). ! ! Here the min thickness is hrmin = min(2*hi, hi+maxraft), ! so thick ridging ice is not required to raft. ! !----------------------------------------------------------------- do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny) then hi = vicen(i,j,n) / aicen(i,j,n) hrmin(ij,n) = min(c2*hi, hi + maxraft) hrmax(ij,n) = c2*sqrt(Hstar*hi) hrmax(ij,n) = max(hrmax(ij,n), hrmin(ij,n)+puny) hrmean = p5 * (hrmin(ij,n) + hrmax(ij,n)) krdg(ij,n) = hrmean / hi ! diagnostic rafting mask not implemented endif enddo ! ij enddo ! n else ! krdg_redist = 1; exponential redistribution !----------------------------------------------------------------- ! The ridge ITD is a negative exponential: ! ! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin ! ! where hrmin is the minimum thickness of ridging ice and ! hrexp is the e-folding thickness. ! ! Here, assume as above that hrmin = min(2*hi, hi+maxraft). ! That is, the minimum ridge thickness results from rafting, ! unless the ice is thicker than maxraft. ! ! Also, assume that hrexp = mu_rdg*sqrt(hi). ! The parameter mu_rdg is tuned to give e-folding scales mostly ! in the range 2-4 m as observed by upward-looking sonar. ! ! Values of mu_rdg in the right column give ice strengths ! roughly equal to values of Hstar in the left column ! (within ~10 kN/m for typical ITDs): ! ! Hstar mu_rdg ! ! 25 3.0 ! 50 4.0 ! 75 5.0 ! 100 6.0 !----------------------------------------------------------------- do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny) then hi = vicen(i,j,n) / aicen(i,j,n) hi = max(hi,puny) hrmin(ij,n) = min(c2*hi, hi + maxraft) hrexp(ij,n) = mu_rdg * sqrt(hi) krdg(ij,n) = (hrmin(ij,n) + hrexp(ij,n)) / hi !echmod: check computational efficiency ! diagnostic rafting mask if (present(mraft)) then mraft(ij,n) = max(c0, sign(c1, hi+maxraft-hrmin(ij,n))) xtmp = mraft(ij,n)*((c2*hi+hrexp(ij,n))/hi - krdg(ij,n)) mraft(ij,n) = max(c0, sign(c1, puny-abs(xtmp))) endif endif enddo enddo endif ! krdg_redist !---------------------------------------------------------------- ! Compute aksum = net ice area removed / total area participating. ! For instance, if a unit area of ice with h = 1 participates in ! ridging to form a ridge with a = 1/3 and h = 3, then ! aksum = 1 - 1/3 = 2/3. !---------------------------------------------------------------- do ij = 1, icells aksum(ij) = apartic(ij,0) ! area participating = area removed enddo do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells ! area participating > area removed aksum(ij) = aksum(ij) & + apartic(ij,n) * (c1 - c1/krdg(ij,n)) enddo enddo ! diagnostics if (present(aparticn)) then do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) aparticn(i,j,n) = apartic(ij,n) enddo enddo endif if (present(krdgn)) then do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) krdgn(i,j,n) = krdg(ij,n) enddo enddo endif end subroutine ridge_itd !======================================================================= ! Remove area, volume, and energy from each ridging category ! and add to thicker ice categories. ! ! Tracers: Ridging conserves ice volume and therefore conserves volume ! tracers. It does not conserve ice area, and therefore a portion of area ! tracers are lost (corresponding to the net closing). Area tracers on ! ice that participates in ridging are carried onto the resulting ridged ! ice (except the portion that are lost due to closing). Therefore, ! tracers must be decremented if they are lost to the ocean during ridging ! (e.g. snow, ponds) or if they are being carried only on the level ice ! area. ! ! author: William H. Lipscomb, LANL subroutine ridge_shift (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, dt, & aicen, trcrn, & vicen, vsnon, & aice0, trcr_depend, & aksum, apartic, & hrmin, hrmax, & hrexp, krdg, & closing_net, opning, & ardg1, ardg2, & virdg, aopen, & ardg1nn, ardg2nn, & virdgnn, & msnow_mlt, esnow_mlt, & maero, mpond, & l_stop, & istop, jstop, & aredistn, vredistn) use ice_state, only: nt_qsno, & nt_alvl, nt_vlvl, nt_aero, tr_lvl, tr_aero, & nt_apnd, nt_hpnd, nt_ipnd, tr_pond, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo, & nt_fbri integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells , & ! number of cells with ice present ntrcr ! number of tracers in use integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), intent(in) :: & dt ! time step (s) integer (kind=int_kind), dimension (ntrcr), intent(in) :: & trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & aice0 ! concentration of open water real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(inout) :: & aicen , & ! concentration of ice vicen , & ! volume per unit area of ice (m) vsnon ! volume per unit area of snow (m) real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (icells), intent(in) :: & aksum ! ratio of area removed to area ridged real (kind=dbl_kind), dimension (icells,0:ncat), intent(in) :: & apartic ! participation function; fraction of ridging ! and closing associated w/ category n real (kind=dbl_kind), dimension (icells,ncat), intent(in) :: & hrmin , & ! minimum ridge thickness hrmax , & ! maximum ridge thickness (krdg_redist = 0) hrexp , & ! ridge e-folding thickness (krdg_redist = 1) krdg ! mean ridge thickness/thickness of ridging ice real (kind=dbl_kind), dimension(icells), intent(inout) :: & closing_net, & ! net rate at which area is removed (1/s) opning , & ! rate of opening due to divergence/shear (1/s) ardg1 , & ! fractional area loss by ridging ice ardg2 , & ! fractional area gain by new ridges virdg , & ! ice volume ridged (m) aopen ! area opened due to divergence/shear real (kind=dbl_kind), dimension(icells,ncat), intent(inout) :: & ardg1nn , & ! area of ice ridged ardg2nn , & ! area of new ridges virdgnn ! ridging ice volume real (kind=dbl_kind), dimension(icells), intent(inout) :: & msnow_mlt, & ! mass of snow added to ocean (kg m-2) esnow_mlt, & ! energy needed to melt snow in ocean (J m-2) mpond ! mass of pond added to ocean (kg m-2) real (kind=dbl_kind), dimension(icells,max_aero), intent(inout) :: & maero ! aerosol mass added to ocean (kg m-2) logical (kind=log_kind), intent(inout) :: & l_stop ! if true, abort on return integer (kind=int_kind), intent(inout) :: & istop, jstop ! indices of grid cell where model aborts real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(inout), optional :: & aredistn , & ! redistribution function: fraction of new ridge area vredistn ! redistribution function: fraction of new ridge volume ! local variables integer (kind=int_kind) :: & i,j , & ! horizontal indices n, nr , & ! thickness category indices k , & ! ice layer index it , & ! tracer index ij, m , & ! horizontal indices, combine i and j loops iridge ! number of cells with nonzero ridging integer (kind=int_kind), dimension (icells) :: & indxii, indxjj, & ! compressed indices indxij ! compressed indices real (kind=dbl_kind), dimension (icells,ncat) :: & aicen_init , & ! ice area before ridging vicen_init , & ! ice volume before ridging vsnon_init ! snow volume before ridging real (kind=dbl_kind), dimension(icells,ntrcr,ncat) :: & atrcrn ! aicen*trcrn real (kind=dbl_kind), dimension (icells) :: & closing_gross ! rate at which area removed, not counting ! area of new ridges ! ECH note: the following arrays only need be defined on iridge cells real (kind=dbl_kind), dimension (icells) :: & afrac , & ! fraction of category area ridged ardg1n , & ! area of ice ridged ardg2n , & ! area of new ridges virdgn , & ! ridging ice volume vsrdgn , & ! ridging snow volume dhr , & ! hrmax - hrmin dhr2 , & ! hrmax^2 - hrmin^2 farea , & ! fraction of new ridge area going to nr fvol ! fraction of new ridge volume going to nr real (kind=dbl_kind) :: & esrdgn ! ridging snow energy real (kind=dbl_kind) :: & hi1 , & ! thickness of ridging ice hexp , & ! ridge e-folding thickness hL, hR , & ! left and right limits of integration expL, expR , & ! exponentials involving hL, hR tmpfac , & ! factor by which opening/closing rates are cut wk1 ! work variable !----------------------------------------------------------------- ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn !----------------------------------------------------------------- do n = 1, ncat do it = 1, ntrcr if (trcr_depend(it) == 0) then ! ice area tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = aicen(i,j,n)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 1) then ! ice volume tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = vicen(i,j,n)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2) then ! snow volume tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = vsnon(i,j,n)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = vicen(i,j,n) & * trcrn(i,j,nt_fbri,n) & * trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_alvl) then ! level ice tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = aicen(i,j,n) & * trcrn(i,j,nt_alvl,n) & * trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_apnd .and. & (tr_pond_cesm .or. tr_pond_topo)) then ! CESM or topo pond area tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = aicen(i,j,n) & * trcrn(i,j,nt_apnd,n) & * trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_apnd .and. & tr_pond_lvl) then ! level-ice pond area tracer do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcrn(ij,it,n) = aicen(i,j,n) & * trcrn(i,j,nt_alvl,n) & * trcrn(i,j,nt_apnd,n) & * trcrn(i,j,it,n) enddo endif enddo enddo !----------------------------------------------------------------- ! Based on the ITD of ridging and ridged ice, convert the net ! closing rate to a gross closing rate. ! NOTE: 0 < aksum <= 1 !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) closing_gross(ij) = closing_net(ij) / aksum(ij) !----------------------------------------------------------------- ! Reduce the closing rate if more than 100% of the open water ! would be removed. Reduce the opening rate proportionately. !----------------------------------------------------------------- if (apartic(ij,0) > c0) then wk1 = apartic(ij,0) * closing_gross(ij) * dt if (wk1 > aice0(i,j)) then tmpfac = aice0(i,j) / wk1 closing_gross(ij) = closing_gross(ij) * tmpfac opning(ij) = opning(ij) * tmpfac endif endif enddo ! ij !----------------------------------------------------------------- ! Reduce the closing rate if more than 100% of any ice category ! would be removed. Reduce the opening rate proportionately. !----------------------------------------------------------------- do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0) then wk1 = apartic(ij,n) * closing_gross(ij) * dt if (wk1 > aicen(i,j,n)) then tmpfac = aicen(i,j,n) / wk1 closing_gross(ij) = closing_gross(ij) * tmpfac opning(ij) = opning(ij) * tmpfac endif endif enddo ! ij enddo ! n !----------------------------------------------------------------- ! Compute change in open water area due to closing and opening. !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) aice0(i,j) = aice0(i,j) & - apartic(ij,0)*closing_gross(ij)*dt & + opning(ij)*dt if (aice0(i,j) < -puny) then l_stop = .true. istop = i jstop = j write (nu_diag,*) ' ' write (nu_diag,*) 'Ridging error: aice0 < 0' write (nu_diag,*) 'i, j, aice0:', i, j, aice0(i,j) return elseif (aice0(i,j) < c0) then ! roundoff error aice0(i,j) = c0 endif aopen(ij) = opning(ij)*dt ! optional diagnostic enddo !----------------------------------------------------------------- ! Save initial state variables !----------------------------------------------------------------- do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) aicen_init(ij,n) = aicen(i,j,n) vicen_init(ij,n) = vicen(i,j,n) vsnon_init(ij,n) = vsnon(i,j,n) enddo enddo !----------------------------------------------------------------- ! Compute the area, volume, and energy of ice ridging in each ! category, along with the area of the resulting ridge. !----------------------------------------------------------------- do n = 1, ncat !----------------------------------------------------------------- ! Identify grid cells with nonzero ridging !----------------------------------------------------------------- iridge = 0 do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen_init(ij,n) > puny .and. apartic(ij,n) > c0 & .and. closing_gross(ij) > c0) then iridge = iridge + 1 indxii(iridge) = i indxjj(iridge) = j indxij(iridge) = ij endif enddo ! ij !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) !----------------------------------------------------------------- ! Compute area of ridging ice (ardg1n) and of new ridge (ardg2n). ! Make sure ridging fraction <=1. (Roundoff errors can give ! ardg1 slightly greater than aicen.) !----------------------------------------------------------------- ardg1n(ij) = apartic(m,n)*closing_gross(m)*dt if (ardg1n(ij) > aicen_init(m,n) + puny) then l_stop = .true. istop = i jstop = j write (nu_diag,*) ' ' write (nu_diag,*) 'Ridging error: ardg > aicen' write (nu_diag,*) 'i, j, n:', i, j, n write (nu_diag,*) 'ardg, aicen:', & ardg1n(ij), aicen_init(m,n) return else ardg1n(ij) = min(aicen_init(m,n), ardg1n(ij)) endif ardg2n(ij) = ardg1n(ij) / krdg(m,n) afrac(ij) = ardg1n(ij) / aicen_init(m,n) !----------------------------------------------------------------- ! Subtract area, volume, and energy from ridging category n. ! Note: Tracer values are unchanged. !----------------------------------------------------------------- virdgn(ij) = vicen_init(m,n) * afrac(ij) vsrdgn(ij) = vsnon_init(m,n) * afrac(ij) aicen(i,j,n) = aicen(i,j,n) - ardg1n(ij) vicen(i,j,n) = vicen(i,j,n) - virdgn(ij) vsnon(i,j,n) = vsnon(i,j,n) - vsrdgn(ij) !----------------------------------------------------------------- ! Increment ridging diagnostics !----------------------------------------------------------------- ardg1(m) = ardg1(m) + ardg1n(ij) ardg2(m) = ardg2(m) + ardg2n(ij) virdg(m) = virdg(m) + virdgn(ij) ardg1nn(m,n) = ardg1n(ij) ardg2nn(m,n) = ardg2n(ij) virdgnn(m,n) = virdgn(ij) !----------------------------------------------------------------- ! Place part of the snow and tracer lost by ridging into the ocean. !----------------------------------------------------------------- msnow_mlt(m) = msnow_mlt(m) + rhos*vsrdgn(ij)*(c1-fsnowrdg) if (tr_aero) then do it = 1, n_aero maero(m,it) = maero(m,it) & + vsrdgn(ij)*(c1-fsnowrdg) & *(trcrn(i,j,nt_aero +4*(it-1),n) & + trcrn(i,j,nt_aero+1+4*(it-1),n)) enddo endif if (tr_pond_topo) then mpond(m) = mpond(m) + ardg1n(ij) & * trcrn(i,j,nt_apnd,n) & * trcrn(i,j,nt_hpnd,n) endif !----------------------------------------------------------------- ! Compute quantities used to apportion ice among categories ! in the nr loop below !----------------------------------------------------------------- dhr(ij) = hrmax(ij,n) - hrmin(m,n) dhr2(ij) = hrmax(ij,n) * hrmax(ij,n) - hrmin(m,n) * hrmin(m,n) enddo ! ij !----------------------------------------------------------------- ! Increment energy needed to melt snow in ocean. ! Note that esnow_mlt < 0; the ocean must cool to melt snow. !----------------------------------------------------------------- do k = 1, nslyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) esrdgn = vsrdgn(ij) * trcrn(i,j,nt_qsno+k-1,n) & / real(nslyr,kind=dbl_kind) esnow_mlt(m) = esnow_mlt(m) + esrdgn*(c1-fsnowrdg) enddo enddo !----------------------------------------------------------------- ! Subtract area- and volume-weighted tracers from category n. !----------------------------------------------------------------- do it = 1, ntrcr if (trcr_depend(it) == 0) then ! ice area tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - ardg1n(ij)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 1) then ! ice volume tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - virdgn(ij)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2) then ! snow volume tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - vsrdgn(ij)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_alvl) then ! level ice tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - ardg1n(ij)*trcrn(i,j,nt_alvl,n)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_apnd .and. & (tr_pond_cesm .or. tr_pond_topo)) then ! CESM or topo pond area tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - ardg1n(ij)*trcrn(i,j,nt_apnd,n)*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_apnd .and. & tr_pond_lvl) then ! level-ice pond area tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - ardg1n(ij) & * trcrn(i,j,nt_alvl,n) & * trcrn(i,j,nt_apnd,n) & * trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,n) = atrcrn(m,it,n) & - virdgn(ij)*trcrn(i,j,it,n)*trcrn(i,j,nt_fbri,n) enddo endif ! trcr_depend enddo ! ntrcr !----------------------------------------------------------------- ! Add area, volume, and energy of new ridge to each category nr. !----------------------------------------------------------------- do nr = 1, ncat if (krdg_redist == 0) then ! Hibler 1980 formulation do ij = 1, iridge m = indxij(ij) !----------------------------------------------------------------- ! Compute the fraction of ridged ice area and volume going to ! thickness category nr. !----------------------------------------------------------------- if (hrmin(m,n) >= hin_max(nr) .or. & hrmax(ij,n) <= hin_max(nr-1)) then hL = c0 hR = c0 else hL = max (hrmin(m,n), hin_max(nr-1)) hR = min (hrmax(ij,n), hin_max(nr)) endif farea(ij) = (hR-hL) / dhr(ij) fvol (ij) = (hR*hR - hL*hL) / dhr2(ij) enddo ! ij else ! krdg_redist = 1; 2005 exponential formulation !----------------------------------------------------------------- ! Compute the fraction of ridged ice area and volume going to ! thickness category nr. !----------------------------------------------------------------- if (nr < ncat) then do ij = 1, iridge m = indxij(ij) hi1 = hrmin(m,n) hexp = hrexp(m,n) if (hi1 >= hin_max(nr)) then farea(ij) = c0 fvol (ij) = c0 else hL = max (hi1, hin_max(nr-1)) hR = hin_max(nr) expL = exp(-(hL-hi1)/hexp) expR = exp(-(hR-hi1)/hexp) farea(ij) = expL - expR fvol (ij) = ((hL + hexp)*expL & - (hR + hexp)*expR) / (hi1 + hexp) endif enddo ! ij else ! nr = ncat do ij = 1, iridge m = indxij(ij) hi1 = hrmin(m,n) hexp = hrexp(m,n) hL = max (hi1, hin_max(nr-1)) expL = exp(-(hL-hi1)/hexp) farea(ij) = expL fvol (ij) = (hL + hexp)*expL / (hi1 + hexp) enddo endif ! nr < ncat ! diagnostics if (n ==1) then ! only for thinnest ridging ice if (present(aredistn)) then do ij = 1, iridge i = indxii(ij) j = indxjj(ij) aredistn(i,j,nr) = farea(ij)*ardg2n(ij) enddo endif if (present(vredistn)) then do ij = 1, iridge i = indxii(ij) j = indxjj(ij) vredistn(i,j,nr) = fvol(ij)*virdgn(ij) enddo endif endif endif ! krdg_redist !----------------------------------------------------------------- ! Transfer ice area, ice volume, and snow volume to category nr. !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) aicen(i,j,nr) = aicen(i,j,nr) + farea(ij)*ardg2n(ij) vicen(i,j,nr) = vicen(i,j,nr) + fvol(ij) *virdgn(ij) vsnon(i,j,nr) = vsnon(i,j,nr) & + fvol(ij)*vsrdgn(ij)*fsnowrdg enddo !----------------------------------------------------------------- ! Transfer area-weighted and volume-weighted tracers to category nr. ! Note: The global sum aicen*trcrn of ice area tracers ! (trcr_depend = 0) is not conserved by ridging. ! However, ridging conserves the global sum of volume ! tracers (trcr_depend = 1 or 2). ! Tracers associated with level ice, or that are otherwise lost ! from ridging ice, are not transferred. ! We assume that all pond water is lost from ridging ice. !----------------------------------------------------------------- do it = 1, ntrcr if (trcr_depend(it) == 0) then ! ice area tracer if (it /= nt_alvl) then !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,nr) = atrcrn(m,it,nr) & + farea(ij)*ardg2n(ij)*trcrn(i,j,it,n) enddo endif elseif (trcr_depend(it) == 1) then ! ice volume tracer if (it /= nt_vlvl) then !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,nr) = atrcrn(m,it,nr) & + fvol(ij)*virdgn(ij)*trcrn(i,j,it,n) enddo endif elseif (trcr_depend(it) == 2) then ! snow volume tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,nr) = atrcrn(m,it,nr) & + fvol(ij)*vsrdgn(ij)*fsnowrdg*trcrn(i,j,it,n) enddo elseif (trcr_depend(it) == 2+nt_fbri) then ! brine tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, iridge i = indxii(ij) j = indxjj(ij) m = indxij(ij) atrcrn(m,it,nr) = atrcrn(m,it,nr) & ! + fvol(ij)*virdgn(ij)*trcrn(i,j,nt_fbri,n)*trcrn(i,j,it,n) enddo endif ! trcr_depend enddo ! ntrcr enddo ! nr (new ridges) enddo ! n (ridging categories) !----------------------------------------------------------------- ! Compute new tracers !----------------------------------------------------------------- do n = 1, ncat call compute_tracers (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & atrcrn(:,:,n), aicen(:,:, n), & vicen (:,:, n), vsnon(:,:, n), & trcrn(:,:,:,n)) enddo end subroutine ridge_shift !======================================================================= ! Make sure ice area <=1. If not, prepare to repeat ridging. ! ! authors William H. Lipscomb, LANL subroutine ridge_check (nx_block, ny_block, & icells, indxi, indxj, & dt, & asum, closing_net, & divu_adv, opning, & iterate_ridging) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells with ice present integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed indices for cells with ice real (kind=dbl_kind), intent(in) :: & dt ! time step (s) real (kind=dbl_kind), dimension (icells), intent(in) :: & asum ! sum of ice and open water area real (kind=dbl_kind), dimension (icells), & intent(inout) :: & closing_net, & ! net rate at which area is removed (1/s) divu_adv , & ! divu as implied by transport scheme (1/s) opning ! rate of opening due to divergence/shear logical (kind=log_kind), intent(out) :: & iterate_ridging ! if true, repeat the ridging ! local variables integer (kind=int_kind) :: & ij ! horizontal index, combines i and j loops iterate_ridging = .false. do ij = 1, icells if (abs(asum(ij) - c1) < puny) then closing_net(ij) = c0 opning (ij) = c0 else iterate_ridging = .true. divu_adv(ij) = (c1 - asum(ij)) / dt closing_net(ij) = max(c0, -divu_adv(ij)) opning(ij) = max(c0, divu_adv(ij)) endif enddo end subroutine ridge_check !======================================================================= ! Compute the strength of the ice pack, defined as the energy (J m-2) ! dissipated per unit area removed from the ice pack under compression, ! and assumed proportional to the change in potential energy caused ! by ridging. ! ! See Rothrock (1975) and Hibler (1980). ! ! For simpler strength parameterization, see this reference: ! Hibler, W. D. III, 1979: A dynamic-thermodynamic sea ice model, ! J. Phys. Oceanog., 9, 817-846. ! ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL subroutine ice_strength (nx_block, ny_block, & ilo, ihi, jlo, jhi, & icells, & indxi, indxj, & aice, vice, & aice0, aicen, & vicen, strength) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ilo,ihi,jlo,jhi ! beg and end of physical domain integer (kind=int_kind), intent(in) :: & icells ! no. of cells where icetmask = 1 integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi , & ! compressed index in i-direction indxj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & aice , & ! concentration of ice vice , & ! volume per unit area of ice (m) aice0 ! concentration of open water real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(in) :: & aicen , & ! concentration of ice vicen ! volume per unit area of ice (m) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & strength ! ice strength (N/m) ! local variables real (kind=dbl_kind), dimension (icells) :: & asum , & ! sum of ice and open water area aksum ! ratio of area removed to area ridged real (kind=dbl_kind), dimension (icells,0:ncat) :: & apartic ! participation function; fraction of ridging ! and closing associated w/ category n real (kind=dbl_kind), dimension (icells,ncat) :: & hrmin , & ! minimum ridge thickness hrmax , & ! maximum ridge thickness (krdg_redist = 0) hrexp , & ! ridge e-folding thickness (krdg_redist = 1) krdg ! mean ridge thickness/thickness of ridging ice integer (kind=int_kind) :: & i,j , & ! horizontal indices n , & ! thickness category index ij ! horizontal index, combines i and j loops real (kind=dbl_kind) :: & hi , & ! ice thickness (m) h2rdg , & ! mean value of h^2 for new ridge dh2rdg ! change in mean value of h^2 per unit area ! consumed by ridging !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- strength(:,:) = c0 if (kstrength == 1) then ! Rothrock '75 formulation !----------------------------------------------------------------- ! Compute thickness distribution of ridging and ridged ice. !----------------------------------------------------------------- call asum_ridging (nx_block, ny_block, & icells, indxi, indxj, & aicen, aice0, & asum) call ridge_itd (nx_block, ny_block, & icells, indxi, indxj, & aicen, vicen, & aice0, & aksum, apartic, & hrmin, hrmax, & hrexp, krdg) !----------------------------------------------------------------- ! Compute ice strength based on change in potential energy, ! as in Rothrock (1975) !----------------------------------------------------------------- if (krdg_redist==0) then ! Hibler 1980 formulation do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0)then hi = vicen(i,j,n) / aicen(i,j,n) h2rdg = p333 * (hrmax(ij,n)**3 - hrmin(ij,n)**3) & / (hrmax(ij,n) - hrmin(ij,n)) dh2rdg = -hi*hi + h2rdg/krdg(ij,n) strength(i,j) = strength(i,j) & + apartic(ij,n) * dh2rdg endif ! aicen > puny enddo ! ij enddo ! n elseif (krdg_redist==1) then ! exponential formulation do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny .and. apartic(ij,n) > c0)then hi = vicen(i,j,n) / aicen(i,j,n) h2rdg = hrmin(ij,n)*hrmin(ij,n) & + c2*hrmin(ij,n)*hrexp(ij,n) & + c2*hrexp(ij,n)*hrexp(ij,n) dh2rdg = -hi*hi + h2rdg/krdg(ij,n) strength(i,j) = strength(i,j) & + apartic(ij,n) * dh2rdg endif enddo ! ij enddo ! n endif ! krdg_redist !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) strength(i,j) = Cf * Cp * strength(i,j) / aksum(ij) ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) ! Cf accounts for frictional dissipation enddo ! ij else ! kstrength /= 1: Hibler (1979) form !----------------------------------------------------------------- ! Compute ice strength as in Hibler (1979) !----------------------------------------------------------------- do j = jlo, jhi do i = ilo, ihi strength(i,j) = Pstar*vice(i,j)*exp(-Cstar*(c1-aice(i,j))) enddo ! j enddo ! i endif ! kstrength end subroutine ice_strength !======================================================================= end module ice_mechred !=======================================================================