!======================================================================= ! Routines to initialize the ice thickness distribution and ! utilities to redistribute ice among categories. These routines ! are not specific to a particular numerical implementation. ! ! See Bitz, C.M., and W.H. Lipscomb, 1999: ! An energy-conserving thermodynamic model of sea ice, ! J. Geophys. Res., 104, 15,669--15,677. ! ! See Bitz, C.M., M.M. Holland, A.J. Weaver, M. Eby, 2001: ! Simulating the ice-thickness distribution in a climate model, ! J. Geophys. Res., 106, 2441--2464. ! ! authors: C. M. Bitz, UW ! William H. Lipscomb and Elizabeth C. Hunke, LANL ! ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb ! ! 2004 WHL: Added multiple snow layers, block structure, cleanup_itd ! 2006 ECH: Added WMO standard ice thickness categories as kcatbound=2 ! Streamlined for efficiency ! Converted to free source form (F90) ! 2014 ECH: Converted to column package module icepack_itd use icepack_kinds use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers use icepack_parameters, only: kcatbound, kitd use icepack_therm_shared, only: Tmin, hi_min use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted use icepack_zbgc_shared, only: zap_small_bgc use icepack_tracers, only: icepack_compute_tracers implicit none private public :: aggregate_area, & shift_ice, & column_sum, & column_conservation_check, & cleanup_itd, & reduce_area, & icepack_init_itd, & icepack_init_itd_hist, & icepack_aggregate !======================================================================= contains !======================================================================= ! Aggregate ice area (but not other state variables) over thickness ! categories. ! ! authors: William H. Lipscomb, LANL subroutine aggregate_area (ncat, aicen, aice, aice0) integer (kind=int_kind), intent(in) :: & ncat ! number of thickness categories real (kind=dbl_kind), dimension(:), intent(in) :: & aicen ! concentration of ice real (kind=dbl_kind), intent(inout) :: & aice, & ! concentration of ice aice0 ! concentration of open water ! local variables integer (kind=int_kind) :: n character(len=*),parameter :: subname='(aggregate_area)' !----------------------------------------------------------------- ! Aggregate !----------------------------------------------------------------- aice = c0 do n = 1, ncat aice = aice + aicen(n) enddo ! n ! open water fraction aice0 = max (c1 - aice, c0) end subroutine aggregate_area !======================================================================= ! Rebins thicknesses into defined categories ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL subroutine rebin (ntrcr, trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata, & aicen, trcrn, & vicen, vsnon, & ncat, hin_max ) integer (kind=int_kind), intent(in) :: & ntrcr , & ! number of tracers in use ncat ! number of thickness categories integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon n_trcr_strata ! number of underlying tracer layers real (kind=dbl_kind), dimension (:,:), intent(in) :: & trcr_base ! = 0 or 1 depending on tracer dependency ! argument 2: (1) aice, (2) vice, (3) vsno integer (kind=int_kind), dimension (:,:), intent(in) :: & nt_strata ! indices of underlying tracer layers real (kind=dbl_kind), dimension (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 (:,:), intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension(0:ncat), intent(in) :: & hin_max ! category limits (m) ! local variables integer (kind=int_kind) :: & n ! category index logical (kind=log_kind) :: & shiftflag ! = .true. if ice must be shifted integer (kind=int_kind), dimension (ncat) :: & donor ! donor category index real (kind=dbl_kind), dimension (ncat) :: & daice , & ! ice area transferred dvice , & ! ice volume transferred hicen ! ice thickness for each cat (m) character(len=*),parameter :: subname='(rebin)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- do n = 1, ncat donor(n) = 0 daice(n) = c0 dvice(n) = c0 !----------------------------------------------------------------- ! Compute ice thickness. !----------------------------------------------------------------- if (aicen(n) > puny) then hicen(n) = vicen(n) / aicen(n) else hicen(n) = c0 endif enddo ! n !----------------------------------------------------------------- ! make sure thickness of cat 1 is at least hin_max(0) !----------------------------------------------------------------- if (aicen(1) > puny) then if (hicen(1) <= hin_max(0) .and. hin_max(0) > c0 ) then aicen(1) = vicen(1) / hin_max(0) hicen(1) = hin_max(0) endif endif !----------------------------------------------------------------- ! If a category thickness is not in bounds, shift the ! entire area, volume, and energy to the neighboring category !----------------------------------------------------------------- !----------------------------------------------------------------- ! Move thin categories up !----------------------------------------------------------------- do n = 1, ncat-1 ! loop over category boundaries !----------------------------------------------------------------- ! identify thicknesses that are too big !----------------------------------------------------------------- shiftflag = .false. if (aicen(n) > puny .and. & hicen(n) > hin_max(n)) then shiftflag = .true. donor(n) = n daice(n) = aicen(n) dvice(n) = vicen(n) endif if (shiftflag) then !----------------------------------------------------------------- ! shift ice between categories !----------------------------------------------------------------- call shift_ice (ntrcr, ncat, & trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata, & aicen, trcrn, & vicen, vsnon, & hicen, donor, & daice, dvice ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! reset shift parameters !----------------------------------------------------------------- donor(n) = 0 daice(n) = c0 dvice(n) = c0 endif ! shiftflag enddo ! n !----------------------------------------------------------------- ! Move thick categories down !----------------------------------------------------------------- do n = ncat-1, 1, -1 ! loop over category boundaries !----------------------------------------------------------------- ! identify thicknesses that are too small !----------------------------------------------------------------- shiftflag = .false. if (aicen(n+1) > puny .and. & hicen(n+1) <= hin_max(n)) then shiftflag = .true. donor(n) = n+1 daice(n) = aicen(n+1) dvice(n) = vicen(n+1) endif if (shiftflag) then !----------------------------------------------------------------- ! shift ice between categories !----------------------------------------------------------------- call shift_ice (ntrcr, ncat, & trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata, & aicen, trcrn, & vicen, vsnon, & hicen, donor, & daice, dvice ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! reset shift parameters !----------------------------------------------------------------- donor(n) = 0 daice(n) = c0 dvice(n) = c0 endif ! shiftflag enddo ! n end subroutine rebin !======================================================================= ! Reduce area when ice melts for special case of ncat=1 ! ! Use CSM 1.0-like method of reducing ice area ! when melting occurs: assume only half the ice volume ! change goes to thickness decrease, the other half ! to reduction in ice fraction ! ! authors: C. M. Bitz, UW ! modified by: Elizabeth C. Hunke, LANL subroutine reduce_area (hin_max, & aicen, vicen, & aicen_init,vicen_init) real (kind=dbl_kind), intent(in) :: & hin_max ! lowest category boundary real (kind=dbl_kind), intent(inout) :: & aicen , & ! concentration of ice vicen ! volume per unit area of ice (m) real (kind=dbl_kind), intent(in) :: & aicen_init, & ! old ice area for category 1 (m) vicen_init ! old ice volume for category 1 (m) ! local variables real (kind=dbl_kind) :: & hi0 , & ! initial hi hi1 , & ! current hi dhi ! hi1 - hi0 character(len=*),parameter :: subname='(reduce_area)' hi0 = c0 if (aicen_init > c0) & hi0 = vicen_init / aicen_init hi1 = c0 if (aicen > c0) & hi1 = vicen / aicen ! make sure thickness of cat 1 is at least hin_max(0) if (hi1 <= hin_max .and. hin_max > c0 ) then aicen = vicen / hin_max hi1 = hin_max endif if (aicen > c0) then dhi = hi1 - hi0 if (dhi < c0) then hi1 = vicen / aicen aicen = c2 * vicen / (hi1 + hi0) endif endif end subroutine reduce_area !======================================================================= ! Shift ice across category boundaries, conserving area, volume, and ! energy. ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL subroutine shift_ice (ntrcr, ncat, & trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata, & aicen, trcrn, & vicen, vsnon, & hicen, donor, & daice, dvice ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories ntrcr ! number of tracers in use integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon n_trcr_strata ! number of underlying tracer layers real (kind=dbl_kind), dimension (:,:), intent(in) :: & trcr_base ! = 0 or 1 depending on tracer dependency ! argument 2: (1) aice, (2) vice, (3) vsno integer (kind=int_kind), dimension (:,:), intent(in) :: & nt_strata ! indices of underlying tracer layers real (kind=dbl_kind), dimension (:), 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 (:,:), intent(inout) :: & trcrn ! ice tracers ! NOTE: Third index of donor, daice, dvice should be ncat-1, ! except that compilers would have trouble when ncat = 1 integer (kind=int_kind), dimension(:), intent(in) :: & donor ! donor category index real (kind=dbl_kind), dimension(:), intent(inout) :: & daice , & ! ice area transferred across boundary dvice , & ! ice volume transferred across boundary hicen ! ice thickness for each cat (m) ! local variables integer (kind=int_kind) :: & n , & ! thickness category index nr , & ! receiver category nd , & ! donor category it , & ! tracer index ntr , & ! tracer index itl ! loop index real (kind=dbl_kind), dimension(ntrcr,ncat) :: & atrcrn ! aicen*trcrn real (kind=dbl_kind) :: & dvsnow , & ! snow volume transferred datrcr ! aicen*train transferred logical (kind=log_kind) :: & daice_negative , & ! true if daice < -puny dvice_negative , & ! true if dvice < -puny daice_greater_aicen, & ! true if daice > aicen dvice_greater_vicen ! true if dvice > vicen real (kind=dbl_kind) :: & worka, workb real (kind=dbl_kind), dimension(ncat) :: aicen_init real (kind=dbl_kind), dimension(ncat) :: vicen_init real (kind=dbl_kind), dimension(ncat) :: vsnon_init character(len=*),parameter :: subname='(shift_ice)' !----------------------------------------------------------------- ! store initial snow and ice volume !----------------------------------------------------------------- aicen_init(:) = aicen(:) vicen_init(:) = vicen(:) vsnon_init(:) = vsnon(:) !----------------------------------------------------------------- ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn !----------------------------------------------------------------- do n = 1, ncat do it = 1, ntrcr atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) & + trcr_base(it,2) * vicen(n) & + trcr_base(it,3) * vsnon(n)) if (n_trcr_strata(it) > 0) then do itl = 1, n_trcr_strata(it) ntr = nt_strata(it,itl) atrcrn(it,n) = atrcrn(it,n) * trcrn(ntr,n) enddo endif enddo ! it enddo ! n !----------------------------------------------------------------- ! Check for daice or dvice out of range, allowing for roundoff error !----------------------------------------------------------------- do n = 1, ncat-1 daice_negative = .false. dvice_negative = .false. daice_greater_aicen = .false. dvice_greater_vicen = .false. if (donor(n) > 0) then nd = donor(n) if (daice(n) < c0) then if (daice(n) > -puny*aicen(nd)) then daice(n) = c0 ! shift no ice dvice(n) = c0 else daice_negative = .true. endif endif if (dvice(n) < c0) then if (dvice(n) > -puny*vicen(nd)) then daice(n) = c0 ! shift no ice dvice(n) = c0 else dvice_negative = .true. endif endif if (daice(n) > aicen(nd)*(c1-puny)) then if (daice(n) < aicen(nd)*(c1+puny)) then daice(n) = aicen(nd) dvice(n) = vicen(nd) else daice_greater_aicen = .true. endif endif if (dvice(n) > vicen(nd)*(c1-puny)) then if (dvice(n) < vicen(nd)*(c1+puny)) then daice(n) = aicen(nd) dvice(n) = vicen(nd) else dvice_greater_vicen = .true. endif endif endif ! donor > 0 !----------------------------------------------------------------- ! error messages !----------------------------------------------------------------- if (daice_negative) then if (donor(n) > 0 .and. & daice(n) <= -puny*aicen(nd)) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'shift_ice: negative daice' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'boundary, donor cat:', n, nd call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'daice =', daice(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dvice =', dvice(n) call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' shift_ice: negative daice') endif endif if (icepack_warnings_aborted(subname)) return if (dvice_negative) then if (donor(n) > 0 .and. & dvice(n) <= -puny*vicen(nd)) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'shift_ice: negative dvice' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'boundary, donor cat:', n, nd call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'daice =', daice(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dvice =', dvice(n) call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' shift_ice: negative dvice') endif endif if (icepack_warnings_aborted(subname)) return if (daice_greater_aicen) then if (donor(n) > 0) then nd = donor(n) if (daice(n) >= aicen(nd)*(c1+puny)) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'shift_ice: daice > aicen' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'boundary, donor cat:', n, nd call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'daice =', daice(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'aicen =', aicen(nd) call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' shift_ice: daice > aicen') endif endif endif if (icepack_warnings_aborted(subname)) return if (dvice_greater_vicen) then if (donor(n) > 0) then nd = donor(n) if (dvice(n) >= vicen(nd)*(c1+puny)) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'shift_ice: dvice > vicen' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'boundary, donor cat:', n, nd call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dvice =', dvice(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'vicen =', vicen(nd) call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' shift_ice: dvice > vicen') endif endif endif if (icepack_warnings_aborted(subname)) return enddo ! boundaries, 1 to ncat-1 !----------------------------------------------------------------- ! transfer volume and energy between categories !----------------------------------------------------------------- do n = 1, ncat-1 if (daice(n) > c0) then ! daice(n) can be < puny nd = donor(n) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif aicen(nd) = aicen(nd) - daice(n) aicen(nr) = aicen(nr) + daice(n) vicen(nd) = vicen(nd) - dvice(n) vicen(nr) = vicen(nr) + dvice(n) worka = daice(n) / aicen_init(nd) dvsnow = vsnon_init(nd) * worka vsnon(nd) = vsnon(nd) - dvsnow vsnon(nr) = vsnon(nr) + dvsnow workb = dvsnow do it = 1, ntrcr nd = donor(n) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif datrcr = trcrn(it,nd)*(trcr_base(it,1) * daice(n) & + trcr_base(it,2) * dvice(n) & + trcr_base(it,3) * workb) if (n_trcr_strata(it) > 0) then do itl = 1, n_trcr_strata(it) ntr = nt_strata(it,itl) datrcr = datrcr * trcrn(ntr,nd) enddo endif atrcrn(it,nd) = atrcrn(it,nd) - datrcr atrcrn(it,nr) = atrcrn(it,nr) + datrcr enddo ! ntrcr endif ! daice enddo ! boundaries, 1 to ncat-1 !----------------------------------------------------------------- ! Update ice thickness and tracers !----------------------------------------------------------------- do n = 1, ncat if (aicen(n) > puny) then hicen(n) = vicen (n) / aicen(n) else hicen(n) = c0 endif !----------------------------------------------------------------- ! Compute new tracers !----------------------------------------------------------------- call icepack_compute_tracers (ntrcr, trcr_depend, & atrcrn(:,n), aicen(n), & vicen(n), vsnon(n), & trcr_base, n_trcr_strata, & nt_strata, trcrn(:,n)) if (icepack_warnings_aborted(subname)) return enddo ! ncat end subroutine shift_ice !======================================================================= ! For each grid cell, sum field over all ice categories. ! ! author: William H. Lipscomb, LANL subroutine column_sum (nsum, xin, xout) integer (kind=int_kind), intent(in) :: & nsum ! number of categories/layers real (kind=dbl_kind), dimension (nsum), & intent(in) :: & xin ! input field real (kind=dbl_kind), intent(out) :: & xout ! output field ! local variables integer (kind=int_kind) :: & n ! category/layer index character(len=*),parameter :: subname='(column_sum)' xout = c0 do n = 1, nsum xout = xout + xin(n) enddo ! n end subroutine column_sum !======================================================================= ! For each physical grid cell, check that initial and final values ! of a conserved field are equal to within a small value. ! ! author: William H. Lipscomb, LANL subroutine column_conservation_check (fieldid, & x1, x2, & max_err ) real (kind=dbl_kind), intent(in) :: & x1 , & ! initial field x2 ! final field real (kind=dbl_kind), intent(in) :: & max_err ! max allowed error character (len=char_len), intent(in) :: & fieldid ! field identifier character(len=*),parameter :: subname='(column_conservation_check)' ! local variables if (abs (x2-x1) > max_err) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Conservation error: ', trim(fieldid) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Initial value =', x1 call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final value =', x2 call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Difference =', x2 - x1 call icepack_warnings_add(warnstr) endif end subroutine column_conservation_check !======================================================================= ! Cleanup subroutine that rebins thickness categories if necessary, ! eliminates very small ice areas while conserving mass and energy, ! aggregates state variables, and does a boundary call. ! It is a good idea to call this subroutine after the thermodynamics ! (thermo_vertical/thermo_itd) and again after the dynamics ! (evp/transport/ridging). ! ! author: William H. Lipscomb, LANL subroutine cleanup_itd (dt, ntrcr, & nilyr, nslyr, & ncat, hin_max, & aicen, trcrn, & vicen, vsnon, & aice0, aice, & n_aero, & nbtrcr, nblyr, & tr_aero, & tr_pond_topo, & heat_capacity, & first_ice, & trcr_depend, trcr_base, & n_trcr_strata,nt_strata, & fpond, fresh, & fsalt, fhocn, & faero_ocn, fzsal, & flux_bio, limit_aice_in) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nblyr , & ! number of bio layers nslyr , & ! number of snow layers ntrcr , & ! number of tracers in use nbtrcr, & ! number of bio tracers in use n_aero ! number of aerosol tracers real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), dimension(0:ncat), intent(in) :: & hin_max ! category boundaries (m) real (kind=dbl_kind), dimension (:), 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 (:,:), intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), intent(inout) :: & aice , & ! total ice concentration aice0 ! concentration of open water integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon n_trcr_strata ! number of underlying tracer layers real (kind=dbl_kind), dimension (:,:), intent(in) :: & trcr_base ! = 0 or 1 depending on tracer dependency ! argument 2: (1) aice, (2) vice, (3) vsno integer (kind=int_kind), dimension (:,:), intent(in) :: & nt_strata ! indices of underlying tracer layers logical (kind=log_kind), intent(in) :: & tr_aero, & ! aerosol flag tr_pond_topo, & ! topo pond flag heat_capacity ! if false, ice and snow have zero heat capacity logical (kind=log_kind), dimension(ncat),intent(inout) :: & first_ice ! For bgc and S tracers. set to true if zapping ice. ! ice-ocean fluxes (required for strict conservation) real (kind=dbl_kind), intent(inout), optional :: & fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) fzsal ! net salt flux to ocean from zsalinity (kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout), optional :: & flux_bio ! net tracer flux to ocean from biology (mmol/m^2/s) real (kind=dbl_kind), dimension (:), & intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) logical (kind=log_kind), intent(in), optional :: & limit_aice_in ! if false, allow aice to be out of bounds ! may want to allow this for unit tests ! local variables integer (kind=int_kind) :: & n , & ! category index it ! tracer index real (kind=dbl_kind) & dfpond , & ! zapped pond water flux (kg/m^2/s) dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) dfhocn , & ! zapped energy flux ( W/m^2) dfzsal ! zapped salt flux for zsalinity (kg/m^2/s) real (kind=dbl_kind), dimension (n_aero) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) real (kind=dbl_kind), dimension (ntrcr) :: & dflux_bio ! zapped biology flux (mmol/m^2/s) logical (kind=log_kind) :: & limit_aice ! if true, check for aice out of bounds character(len=*),parameter :: subname='(cleanup_itd)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- if (present(limit_aice_in)) then limit_aice = limit_aice_in else limit_aice = .true. endif dfpond = c0 dfresh = c0 dfsalt = c0 dfhocn = c0 dfaero_ocn(:) = c0 dflux_bio(:) = c0 dfzsal = c0 !----------------------------------------------------------------- ! Compute total ice area. !----------------------------------------------------------------- call aggregate_area (ncat, aicen, aice, aice0) if (icepack_warnings_aborted(subname)) return if (limit_aice) then ! check for aice out of bounds if (aice > c1+puny .or. aice < -puny) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' aggregate ice area out of bounds') write(warnstr,*) subname, 'aice:', aice call icepack_warnings_add(warnstr) do n = 1, ncat write(warnstr,*) subname, 'n, aicen:', n, aicen(n) call icepack_warnings_add(warnstr) enddo return endif endif ! limit_aice !----------------------------------------------------------------- ! Identify grid cells with ice. !----------------------------------------------------------------- if (aice > puny) then !----------------------------------------------------------------- ! Make sure ice in each category is within its thickness bounds. ! NOTE: The rebin subroutine is needed only in the rare cases ! when the linear_itd subroutine cannot transfer ice ! correctly (e.g., very fast ice growth). !----------------------------------------------------------------- call rebin (ntrcr, trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata, & aicen, trcrn, & vicen, vsnon, & ncat, hin_max ) if (icepack_warnings_aborted(subname)) return endif ! aice > puny !----------------------------------------------------------------- ! Zero out ice categories with very small areas. !----------------------------------------------------------------- if (limit_aice) then call zap_small_areas (dt, ntrcr, & ncat, n_aero, & nblyr, & nilyr, nslyr, & aice, aice0, & aicen, trcrn, & vicen, vsnon, & dfpond, & dfresh, dfsalt, & dfhocn, dfaero_ocn, & tr_aero, tr_pond_topo, & first_ice, nbtrcr, & dfzsal, dflux_bio ) if (icepack_warnings_aborted(subname)) then write(warnstr,*) subname, 'aice:', aice call icepack_warnings_add(warnstr) do n = 1, ncat write(warnstr,*) subname, 'n, aicen:', n, aicen(n) call icepack_warnings_add(warnstr) enddo return endif endif ! l_limit_aice !------------------------------------------------------------------- ! Zap snow that has out of bounds temperatures !------------------------------------------------------------------- call zap_snow_temperature(dt, ncat, & heat_capacity, nblyr, & nslyr, aicen, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & dflux_bio, nbtrcr, & n_aero) if (icepack_warnings_aborted(subname)) return !------------------------------------------------------------------- ! Update ice-ocean fluxes for strict conservation !------------------------------------------------------------------- if (present(fpond)) & fpond = fpond + dfpond if (present(fresh)) & fresh = fresh + dfresh if (present(fsalt)) & fsalt = fsalt + dfsalt if (present(fhocn)) & fhocn = fhocn + dfhocn if (present(faero_ocn)) then do it = 1, n_aero faero_ocn(it) = faero_ocn(it) + dfaero_ocn(it) enddo endif if (present(flux_bio)) then do it = 1, nbtrcr flux_bio (it) = flux_bio(it) + dflux_bio(it) enddo endif if (present(fzsal)) & fzsal = fzsal + dfzsal !---------------------------------------------------------------- ! If using zero-layer model (no heat capacity), check that the ! energy of snow and ice is correct. !---------------------------------------------------------------- if ((.not. heat_capacity) .and. aice > puny) then call zerolayer_check (ncat, nilyr, & nslyr, aicen, & vicen, vsnon, & trcrn) if (icepack_warnings_aborted(subname)) return endif end subroutine cleanup_itd !======================================================================= ! For each ice category in each grid cell, remove ice if the fractional ! area is less than puny. ! ! author: William H. Lipscomb, LANL subroutine zap_small_areas (dt, ntrcr, & ncat, n_aero, & nblyr, & nilyr, nslyr, & aice, aice0, & aicen, trcrn, & vicen, vsnon, & dfpond, & dfresh, dfsalt, & dfhocn, dfaero_ocn, & tr_aero, tr_pond_topo, & first_ice, nbtrcr, & dfzsal, dflux_bio ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nblyr , & ! number of bio layers nslyr , & ! number of snow layers ntrcr , & ! number of tracers in use n_aero , & ! number of aerosol tracers nbtrcr ! number of biology tracers real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), intent(inout) :: & aice , & ! total ice concentration aice0 ! concentration of open water real (kind=dbl_kind), dimension(:), 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 (:,:), intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), intent(inout) :: & dfpond , & ! zapped pond water flux (kg/m^2/s) dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) dfhocn , & ! zapped energy flux ( W/m^2) dfzsal ! zapped salt flux from zsalinity(kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout) :: & dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & tr_aero, & ! aerosol flag tr_pond_topo ! pond flag logical (kind=log_kind), dimension (:),intent(inout) :: & first_ice ! For bgc tracers. Set to true if zapping ice ! local variables integer (kind=int_kind) :: & n, k, it, & !counting indices blevels real (kind=dbl_kind) :: xtmp ! temporary variables real (kind=dbl_kind) , dimension (1):: trcr_skl real (kind=dbl_kind) , dimension (nblyr+1):: bvol character(len=*),parameter :: subname='(zap_small_areas)' !----------------------------------------------------------------- ! I. Zap categories with very small areas. !----------------------------------------------------------------- dfzsal = c0 do n = 1, ncat !----------------------------------------------------------------- ! Count categories to be zapped. !----------------------------------------------------------------- if (aicen(n) < -puny) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' Zap ice: negative ice area') return elseif (abs(aicen(n)) /= c0 .and. & abs(aicen(n)) <= puny) then !----------------------------------------------------------------- ! Account for tracers important for conservation !----------------------------------------------------------------- if (tr_pond_topo) then xtmp = aicen(n) & * trcrn(nt_apnd,n) * trcrn(nt_hpnd,n) dfpond = dfpond - xtmp endif if (tr_aero) then do it = 1, n_aero xtmp = (vicen(n)*(trcrn(nt_aero+2+4*(it-1),n) & + trcrn(nt_aero+3+4*(it-1),n)))/dt dfaero_ocn(it) = dfaero_ocn(it) + xtmp enddo endif if (solve_zsal) then do it = 1, nblyr xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& *trcrn(nt_bgc_S+it-1,n)/ & real(nblyr,kind=dbl_kind)/dt dfzsal = dfzsal + xtmp enddo ! n endif if (skl_bgc .and. nbtrcr > 0) then blevels = 1 bvol(1) = aicen(n)*sk_l it = 1 do it = 1, nbtrcr trcr_skl(1) = trcrn(bio_index(it),n) call zap_small_bgc(blevels, dflux_bio(it), & dt, bvol(1:blevels), trcr_skl(blevels)) if (icepack_warnings_aborted(subname)) return enddo elseif (z_tracers .and. nbtrcr > 0) then blevels = nblyr + 1 bvol(:) = vicen(n)/real(nblyr,kind=dbl_kind)*trcrn(nt_fbri,n) bvol(1) = p5*bvol(1) bvol(blevels) = p5*bvol(blevels) do it = 1, nbtrcr call zap_small_bgc(blevels, dflux_bio(it), & dt, bvol(1:blevels),trcrn(bio_index(it):bio_index(it)+blevels-1,n)) if (icepack_warnings_aborted(subname)) return enddo endif !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- do k = 1, nilyr xtmp = trcrn(nt_qice+k-1,n) / dt & * vicen(n)/real(nilyr,kind=dbl_kind) ! < 0 dfhocn = dfhocn + xtmp trcrn(nt_qice+k-1,n) = c0 enddo ! k !----------------------------------------------------------------- ! Zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- xtmp = (rhoi*vicen(n)) / dt dfresh = dfresh + xtmp xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt dfsalt = dfsalt + xtmp aice0 = aice0 + aicen(n) aicen(n) = c0 vicen(n) = c0 trcrn(nt_Tsfc,n) = Tocnfrz !----------------------------------------------------------------- ! Zap snow !----------------------------------------------------------------- call zap_snow(dt, nslyr, & trcrn(:,n), vsnon(n), & dfresh, dfhocn, & dfaero_ocn, tr_aero, & dflux_bio, nbtrcr, & n_aero, & aicen(n), nblyr) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Zap tracers !----------------------------------------------------------------- if (ntrcr >= 2) then do it = 2, ntrcr if (tr_brine .and. it == nt_fbri) then trcrn(it,n) = c1 else trcrn(it,n) = c0 endif enddo endif first_ice(n) = .true. endif ! aicen enddo ! n !----------------------------------------------------------------- ! II. Count cells with excess ice (aice > c1) due to roundoff errors. ! Zap a little ice in each category so that aice = c1. !----------------------------------------------------------------- if (aice > (c1+puny)) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' Zap ice: excess ice area') return elseif (aice > c1 .and. aice < (c1+puny)) then do n = 1, ncat !----------------------------------------------------------------- ! Account for tracers important for conservation !----------------------------------------------------------------- if (tr_pond_topo) then xtmp = aicen(n) & * trcrn(nt_apnd,n) * trcrn(nt_hpnd,n) & * (aice-c1)/aice dfpond = dfpond - xtmp endif if (tr_aero) then do it = 1, n_aero xtmp = (vsnon(n)*(trcrn(nt_aero +4*(it-1),n) & + trcrn(nt_aero+1+4*(it-1),n)) & + vicen(n)*(trcrn(nt_aero+2+4*(it-1),n) & + trcrn(nt_aero+3+4*(it-1),n))) & * (aice-c1)/aice / dt dfaero_ocn(it) = dfaero_ocn(it) + xtmp enddo ! it endif !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- do k = 1, nilyr xtmp = trcrn(nt_qice+k-1,n) & * vicen(n)/real(nilyr,kind=dbl_kind) & * (aice-c1)/aice / dt ! < 0 dfhocn = dfhocn + xtmp enddo ! k !----------------------------------------------------------------- ! Zap snow energy and use ocean heat to melt snow !----------------------------------------------------------------- do k = 1, nslyr xtmp = trcrn(nt_qsno+k-1,n) & * vsnon(n)/real(nslyr,kind=dbl_kind) & * (aice-c1)/aice / dt ! < 0 dfhocn = dfhocn + xtmp enddo ! k !----------------------------------------------------------------- ! Zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- xtmp = (rhoi*vicen(n) + rhos*vsnon(n)) & * (aice-c1)/aice / dt dfresh = dfresh + xtmp xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 & * (aice-c1)/aice / dt dfsalt = dfsalt + xtmp if (solve_zsal) then do k = 1,nblyr xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& /real(nblyr,kind=dbl_kind)*trcrn(nt_bgc_S+k-1,n) & * (aice-c1)/aice /dt dfzsal = dfzsal + xtmp enddo if (vicen(n) > vicen(n)*trcrn(nt_fbri,n)) then xtmp = (vicen(n)-vicen(n)*trcrn(nt_fbri,n))*(aice-c1)/& aice*p001*rhosi*min_salin/dt dfzsal = dfzsal + xtmp endif endif ! solve_zsal aicen(n) = aicen(n) * (c1/aice) vicen(n) = vicen(n) * (c1/aice) vsnon(n) = vsnon(n) * (c1/aice) ! Note: Tracers are unchanged. enddo ! n !----------------------------------------------------------------- ! Correct aice !----------------------------------------------------------------- aice = c1 aice0 = c0 endif ! aice end subroutine zap_small_areas !======================================================================= subroutine zap_snow(dt, nslyr, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & dflux_bio, nbtrcr, & n_aero, & aicen, nblyr) integer (kind=int_kind), intent(in) :: & nslyr , & ! number of snow layers n_aero , & ! number of aerosol tracers nblyr , & ! number of bio layers nbtrcr real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), dimension (:), intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), intent(in) :: & aicen ! ice area fraction real (kind=dbl_kind), intent(inout) :: & vsnon , & ! volume per unit area of snow (m) dfresh , & ! zapped fresh water flux (kg/m^2/s) dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout) :: & dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & tr_aero ! aerosol flag ! local variables integer (kind=int_kind) :: & k, it ! counting indices real (kind=dbl_kind) :: xtmp, dvssl, dvint character(len=*),parameter :: subname='(zap_snow)' ! aerosols if (tr_aero) then do it = 1, n_aero xtmp = (vsnon*(trcrn(nt_aero +4*(it-1)) & + trcrn(nt_aero+1+4*(it-1))))/dt dfaero_ocn(it) = dfaero_ocn(it) + xtmp enddo ! it endif ! tr_aero if (z_tracers) then dvssl = min(p5*vsnon, hs_ssl*aicen) !snow surface layer dvint = vsnon- dvssl !snow interior do it = 1, nbtrcr xtmp = (trcrn(bio_index(it)+nblyr+1)*dvssl + & trcrn(bio_index(it)+nblyr+2)*dvint)/dt dflux_bio(it) = dflux_bio(it) + xtmp enddo ! it endif ! z_tracers ! snow enthalpy tracer do k = 1, nslyr xtmp = trcrn(nt_qsno+k-1) / dt & * vsnon/real(nslyr,kind=dbl_kind) ! < 0 dfhocn = dfhocn + xtmp trcrn(nt_qsno+k-1) = c0 enddo ! k ! snow volume xtmp = (rhos*vsnon) / dt dfresh = dfresh + xtmp vsnon = c0 end subroutine zap_snow !======================================================================= subroutine zap_snow_temperature(dt, ncat, & heat_capacity, & nblyr, & nslyr, aicen, & trcrn, vsnon, & dfresh, dfhocn, & dfaero_ocn, tr_aero, & dflux_bio, nbtrcr, & n_aero) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nslyr , & ! number of snow layers n_aero, & ! number of aerosol tracers nbtrcr, & ! number of z-tracers in use nblyr ! number of bio layers in ice real (kind=dbl_kind), intent(in) :: & dt ! time step logical (kind=log_kind), intent(in) :: & heat_capacity ! if false, ice and snow have zero heat capacity real (kind=dbl_kind), dimension (:), intent(in) :: & aicen ! concentration of ice real (kind=dbl_kind), dimension(:), intent(inout) :: & vsnon ! volume per unit area of snow (m) real (kind=dbl_kind), dimension (:,:), intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), intent(inout) :: & dfresh , & ! zapped fresh water flux (kg/m^2/s) dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) real (kind=dbl_kind), dimension (:),intent(inout) :: & dflux_bio ! zapped biology flux (mmol/m^2/s) logical (kind=log_kind), intent(in) :: & tr_aero ! aerosol flag ! local variables integer (kind=int_kind) :: & n, k ! counting indices real (kind=dbl_kind) :: & rnslyr , & ! real(nslyr) hsn , & ! snow thickness (m) zqsn , & ! snow layer enthalpy (J m-2) zTsn , & ! snow layer temperature (C) Tmax ! maximum allowed snow temperature logical :: & l_zap ! logical whether zap snow character(len=*),parameter :: subname='(zap_snow_temperature)' rnslyr = real(nslyr,kind=dbl_kind) do n = 1, ncat !----------------------------------------------------------------- ! Determine cells to zap !----------------------------------------------------------------- l_zap = .false. if (aicen(n) > puny) then ! snow thickness hsn = vsnon(n) / aicen(n) ! check each snow layer - zap all if one is bad do k = 1, nslyr ! snow enthalpy and max temperature if (hsn > hs_min .and. heat_capacity) then ! zqsn < 0 zqsn = trcrn(nt_qsno+k-1,n) Tmax = -zqsn*puny*rnslyr / (rhos*cp_ice*vsnon(n)) else zqsn = -rhos * Lfresh Tmax = puny endif ! snow temperature zTsn = (Lfresh + zqsn/rhos)/cp_ice ! check for zapping if (zTsn < Tmin .or. zTsn > Tmax) then l_zap = .true. write(warnstr,*) subname, "zap_snow_temperature: temperature out of bounds!" call icepack_warnings_add(warnstr) write(warnstr,*) subname, "k:" , k call icepack_warnings_add(warnstr) write(warnstr,*) subname, "zTsn:", zTsn call icepack_warnings_add(warnstr) write(warnstr,*) subname, "Tmin:", Tmin call icepack_warnings_add(warnstr) write(warnstr,*) subname, "Tmax:", Tmax call icepack_warnings_add(warnstr) write(warnstr,*) subname, "zqsn:", zqsn call icepack_warnings_add(warnstr) endif enddo ! k endif ! aicen > puny !----------------------------------------------------------------- ! Zap the cells !----------------------------------------------------------------- if (l_zap) & call zap_snow(dt, nslyr, & trcrn(:,n), vsnon(n), & dfresh, dfhocn, & dfaero_ocn, tr_aero, & dflux_bio, nbtrcr, & n_aero, & aicen(n), nblyr) if (icepack_warnings_aborted(subname)) return enddo ! n end subroutine zap_snow_temperature !======================================================================= ! Checks that the snow and ice energy in the zero layer thermodynamics ! model still agrees with the snow and ice volume. ! If there is an error, the model will abort. ! This subroutine is only called if heat_capacity = .false. ! ! author: Alison McLaren, Met Office ! May 2010: ECH replaced eicen, esnon with trcrn but did not test ! the changes. The loop below runs over n=1,ncat and I added loops ! over k, making the test more stringent. subroutine zerolayer_check (ncat, nilyr, & nslyr, aicen, & vicen, vsnon, & trcrn) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), dimension (:), 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 (:,:), intent(inout) :: & trcrn ! ice tracers ! local variables integer (kind=int_kind) :: & k , & ! vertical index n ! category index real (kind=dbl_kind) :: & max_error ! max error in zero layer energy check ! (so max volume error = puny) real (kind=dbl_kind), dimension (ncat) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (ncat) :: & esnon ! energy of melting for each snow layer (J/m^2) logical (kind=log_kind) :: & ice_energy_correct , & ! zero layer ice energy check snow_energy_correct ! zero layer snow energy check real (kind=dbl_kind) :: & worka, workb character(len=*),parameter :: subname='(zerolayer_check)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- max_error = puny*Lfresh*rhos ! max error in zero layer energy check ! (so max volume error = puny) !---------------------------------------------------------------- ! Calculate difference between ice and snow energies and the ! energy values derived from the ice and snow volumes !---------------------------------------------------------------- ice_energy_correct = .true. snow_energy_correct = .true. worka = c0 workb = c0 do n = 1, ncat eicen(n) = c0 do k = 1, nilyr eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) & * vicen(n) / real(nilyr,kind=dbl_kind) enddo worka = eicen(n) + rhoi * Lfresh * vicen(n) esnon(n) = c0 do k = 1, nslyr esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) & * vsnon(n) / real(nslyr,kind=dbl_kind) enddo workb = esnon(n) + rhos * Lfresh * vsnon(n) if(abs(worka) > max_error) ice_energy_correct = .false. if(abs(workb) > max_error) snow_energy_correct = .false. !---------------------------------------------------------------- ! If there is a problem, abort with error message !---------------------------------------------------------------- if (.not. ice_energy_correct) then if (abs(worka) > max_error) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' zerolayer check - wrong ice energy') write(warnstr,*) subname, 'n:', n call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'eicen =', eicen(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'error=', worka call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'vicen =', vicen(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'aicen =', aicen(n) call icepack_warnings_add(warnstr) endif endif if (icepack_warnings_aborted(subname)) return if (.not. snow_energy_correct) then if (abs(workb) > max_error) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//' zerolayer check - wrong snow energy') write(warnstr,*) subname, 'n:', n call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'esnon =', esnon(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'error=', workb call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'vsnon =', vsnon(n) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'aicen =', aicen(n) call icepack_warnings_add(warnstr) return endif endif enddo ! ncat end subroutine zerolayer_check !======================================================================= ! Initialize area fraction and thickness boundaries for the itd model ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL ! C. M. Bitz, UW subroutine icepack_init_itd(ncat, hin_max) integer (kind=int_kind), intent(in) :: & ncat ! number of thickness categories real (kind=dbl_kind), intent(out) :: & hin_max(0:ncat) ! category limits (m) ! local variables integer (kind=int_kind) :: & n ! thickness category index real (kind=dbl_kind) :: & cc1, cc2, cc3, & ! parameters for kcatbound = 0 x1 , & rn , & ! real(n) rncat , & ! real(ncat) d1 , & ! parameters for kcatbound = 1 (m) d2 , & ! b1 , & ! parameters for kcatbound = 3 b2 , & ! b3 real (kind=dbl_kind), dimension(5) :: wmo5 ! data for wmo itd real (kind=dbl_kind), dimension(6) :: wmo6 ! data for wmo itd real (kind=dbl_kind), dimension(7) :: wmo7 ! data for wmo itd character(len=*),parameter :: subname='(icepack_init_itd)' ! thinnest 3 categories combined data wmo5 / 0.30_dbl_kind, 0.70_dbl_kind, & 1.20_dbl_kind, 2.00_dbl_kind, & 999._dbl_kind / ! thinnest 2 categories combined data wmo6 / 0.15_dbl_kind, & 0.30_dbl_kind, 0.70_dbl_kind, & 1.20_dbl_kind, 2.00_dbl_kind, & 999._dbl_kind / !echmod wmo6a ! data wmo6 /0.30_dbl_kind, 0.70_dbl_kind, & ! 1.20_dbl_kind, 2.00_dbl_kind, & ! 4.56729_dbl_kind, & ! 999._dbl_kind / ! all thickness categories data wmo7 / 0.10_dbl_kind, 0.15_dbl_kind, & 0.30_dbl_kind, 0.70_dbl_kind, & 1.20_dbl_kind, 2.00_dbl_kind, & 999._dbl_kind / rncat = real(ncat, kind=dbl_kind) d1 = 3.0_dbl_kind / rncat d2 = 0.5_dbl_kind / rncat b1 = p1 ! asymptotic category width (m) b2 = c3 ! thickness for which participation function is small (m) b3 = max(rncat*(rncat-1), c2*b2/b1) hi_min = p01 ! minimum ice thickness allowed (m) for thermo ! note hi_min is reset to 0.1 for kitd=0, below !----------------------------------------------------------------- ! Choose category boundaries based on one of four options. ! ! The first formula (kcatbound = 0) was used in Lipscomb (2001) ! and in CICE versions 3.0 and 3.1. ! ! The second formula is more user-friendly in the sense that it ! is easy to obtain round numbers for category boundaries: ! ! H(n) = n * [d1 + d2*(n-1)] ! ! Default values are d1 = 300/ncat, d2 = 50/ncat. ! For ncat = 5, boundaries in cm are 60, 140, 240, 360, which are ! close to the standard values given by the first formula. ! For ncat = 10, boundaries in cm are 30, 70, 120, 180, 250, 330, ! 420, 520, 630. ! ! The third option provides support for World Meteorological ! Organization classification based on thickness. The full ! WMO thickness distribution is used if ncat = 7; if ncat=5 ! or ncat = 6, some of the thinner categories are combined. ! For ncat = 5, boundaries are 30, 70, 120, 200, >200 cm. ! For ncat = 6, boundaries are 15, 30, 70, 120, 200, >200 cm. ! For ncat = 7, boundaries are 10, 15, 30, 70, 120, 200, >200 cm. ! ! The fourth formula asymptotes to a particular category width as ! the number of categories increases, given by the parameter b1. ! The parameter b3 is computed so that the category boundaries ! are even numbers. ! ! H(n) = b1 * [n + b3*n*(n+1)/(2*N*(N-1))] for N=ncat ! ! kcatbound=-1 is available only for 1-category runs, with ! boundaries 0 and 100 m. !----------------------------------------------------------------- if (kcatbound == -1) then ! single category hin_max(0) = c0 hin_max(1) = c100 elseif (kcatbound == 0) then ! original scheme if (kitd == 1) then ! linear remapping itd category limits cc1 = c3/rncat cc2 = c15*cc1 cc3 = c3 hin_max(0) = c0 ! minimum ice thickness, m else ! delta function itd category limits #ifndef CESMCOUPLED hi_min = p1 ! minimum ice thickness allowed (m) for thermo #endif cc1 = max(1.1_dbl_kind/rncat,hi_min) cc2 = c25*cc1 cc3 = 2.25_dbl_kind ! hin_max(0) should not be zero ! use some caution in making it less than 0.10 hin_max(0) = hi_min ! minimum ice thickness, m endif ! kitd do n = 1, ncat x1 = real(n-1,kind=dbl_kind) / rncat hin_max(n) = hin_max(n-1) & + cc1 + cc2*(c1 + tanh(cc3*(x1-c1))) enddo elseif (kcatbound == 1) then ! new scheme hin_max(0) = c0 do n = 1, ncat rn = real(n, kind=dbl_kind) hin_max(n) = rn * (d1 + (rn-c1)*d2) enddo elseif (kcatbound == 2) then ! WMO standard if (ncat == 5) then hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo5(n) enddo elseif (ncat == 6) then hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo6(n) enddo elseif (ncat == 7) then hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo7(n) enddo else call icepack_warnings_add(subname//' kcatbound=2 (WMO) must have ncat=5, 6 or 7') call icepack_warnings_setabort(.true.,__FILE__,__LINE__) return endif elseif (kcatbound == 3) then ! asymptotic scheme hin_max(0) = c0 do n = 1, ncat rn = real(n, kind=dbl_kind) hin_max(n) = b1 * (rn + b3*rn*(rn+c1)/(c2*rncat*(rncat-c1))) enddo endif ! kcatbound end subroutine icepack_init_itd !======================================================================= ! Initialize area fraction and thickness boundaries for the itd model ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL ! C. M. Bitz, UW subroutine icepack_init_itd_hist (ncat, hin_max, c_hi_range) integer (kind=int_kind), intent(in) :: & ncat ! number of thickness categories real (kind=dbl_kind), intent(in) :: & hin_max(0:ncat) ! category limits (m) character (len=35), intent(out) :: & c_hi_range(ncat) ! string for history output ! local variables integer (kind=int_kind) :: & n ! thickness category index character(len=8) :: c_hinmax1,c_hinmax2 character(len=2) :: c_nc character(len=*),parameter :: subname='(icepack_init_itd_hist)' write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname call icepack_warnings_add(warnstr) write(warnstr,*) 'hin_max(n-1) < Cat n < hin_max(n)' call icepack_warnings_add(warnstr) do n = 1, ncat write(warnstr,*) hin_max(n-1),' < Cat ',n, ' < ',hin_max(n) call icepack_warnings_add(warnstr) ! Write integer n to character string write (c_nc, '(i2)') n ! Write hin_max to character string write (c_hinmax1, '(f6.3)') hin_max(n-1) write (c_hinmax2, '(f6.3)') hin_max(n) ! Save character string to write to history file c_hi_range(n)=c_hinmax1//'m < hi Cat '//c_nc//' < '//c_hinmax2//'m' enddo write(warnstr,*) ' ' call icepack_warnings_add(warnstr) end subroutine icepack_init_itd_hist !======================================================================= ! Aggregate ice state variables over thickness categories. ! ! authors: C. M. Bitz, UW ! W. H. Lipscomb, LANL subroutine icepack_aggregate (ncat, & aicen, trcrn, & vicen, vsnon, & aice, trcr, & vice, vsno, & aice0, & ntrcr, & trcr_depend, & trcr_base, & n_trcr_strata, & nt_strata) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories ntrcr ! number of tracers in use real (kind=dbl_kind), dimension (:), intent(in) :: & 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 (:,:), & intent(inout) :: & trcrn ! ice tracers integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon n_trcr_strata ! number of underlying tracer layers real (kind=dbl_kind), dimension (:,:), intent(in) :: & trcr_base ! = 0 or 1 depending on tracer dependency ! argument 2: (1) aice, (2) vice, (3) vsno integer (kind=int_kind), dimension (:,:), intent(in) :: & nt_strata ! indices of underlying tracer layers real (kind=dbl_kind), intent(out) :: & aice , & ! concentration of ice vice , & ! volume per unit area of ice (m) vsno , & ! volume per unit area of snow (m) aice0 ! concentration of open water real (kind=dbl_kind), dimension (:), & intent(out) :: & trcr ! ice tracers ! local variables integer (kind=int_kind) :: & n, it, itl, & ! loop indices ntr ! tracer index real (kind=dbl_kind), dimension (:), allocatable :: & atrcr ! sum of aicen*trcrn or vicen*trcrn or vsnon*trcrn real (kind=dbl_kind) :: & atrcrn ! category value character(len=*),parameter :: subname='(icepack_aggregate)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- aice0 = c1 aice = c0 vice = c0 vsno = c0 allocate (atrcr(ntrcr)) !----------------------------------------------------------------- ! Aggregate !----------------------------------------------------------------- atrcr(:) = c0 do n = 1, ncat aice = aice + aicen(n) vice = vice + vicen(n) vsno = vsno + vsnon(n) do it = 1, ntrcr atrcrn = trcrn(it,n)*(trcr_base(it,1) * aicen(n) & + trcr_base(it,2) * vicen(n) & + trcr_base(it,3) * vsnon(n)) if (n_trcr_strata(it) > 0) then ! additional tracer layers do itl = 1, n_trcr_strata(it) ntr = nt_strata(it,itl) atrcrn = atrcrn * trcrn(ntr,n) enddo endif atrcr(it) = atrcr(it) + atrcrn enddo ! ntrcr enddo ! ncat ! Open water fraction aice0 = max (c1 - aice, c0) ! Tracers call icepack_compute_tracers (ntrcr, trcr_depend, & atrcr, aice, & vice , vsno, & trcr_base, n_trcr_strata, & nt_strata, trcr) if (icepack_warnings_aborted(subname)) return deallocate (atrcr) end subroutine icepack_aggregate !======================================================================= end module icepack_itd !=======================================================================