!======================================================================= !BOP ! ! !MODULE: ice_itd - initialize and redistribute ice in the ITD ! ! !DESCRIPTION: ! ! 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. ! ! !REVISION HISTORY: ! SVN:$Id: ice_itd.F90 138 2008-07-08 20:39:37Z eclare $ ! ! 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) ! ! !INTERFACE: ! module ice_itd ! ! !USES: ! use ice_kinds_mod use ice_communicate, only: my_task, master_task use ice_domain_size use ice_constants use ice_fileunits use ice_exit ! !EOP ! implicit none save integer (kind=int_kind) :: & kitd , & ! type of itd conversions ! 0 = delta function ! 1 = linear remap kcatbound , & ! 0 = old category boundary formula ! 1 = new formula giving round numbers ! 2 = WMO standard ilyr1 (ncat), & ! array position of top ice layer in each cat ilyrn (ncat), & ! array position of bottom ice layer in each cat slyr1 (ncat), & ! array position of top snow layer in each cat slyrn (ncat) ! array position of bottom snow layer in each cat real (kind=dbl_kind), parameter :: & hi_min = p01 ! minimum ice thickness allowed (m) real (kind=dbl_kind) :: & hin_max(0:ncat) ! category limits (m) character (len=35) :: c_hi_range(ncat) !------------------------------------------------------------------- ! a note regarding hi_min and hin_max(0): ! both represent a minimum ice thickness. hin_max(0) is ! intended to be used for particular numerical implementations ! of category conversions in the ice thickness distribution. ! hi_min is a more general purpose parameter, but is specifically ! for maintaining stability in the thermodynamics. ! hin_max(0) = 0.1 m for the delta function itd ! hin_max(0) = 0.0 m for linear remapping ! ! Also note that the upper limit on the thickest category ! is only used for the linear remapping scheme ! and it is not a true upper limit on the thickness !------------------------------------------------------------------- !======================================================================= contains !======================================================================= !BOP ! ! !IROUTINE: init_itd - initalize area fraction and thickness boundaries for ITD ! ! !INTERFACE: ! subroutine init_itd ! ! !DESCRIPTION: ! ! Initialize area fraction and thickness boundaries for the itd model ! ! !REVISION HISTORY: ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL ! C. M. Bitz, UW ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! !EOP ! 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 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=8) :: c_hinmax1,c_hinmax2 character(len=2) :: c_nc rncat = real(ncat, kind=dbl_kind) d1 = 3.0_dbl_kind / rncat d2 = 0.5_dbl_kind / rncat !----------------------------------------------------------------- ! Choose category boundaries based on one of three 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. !----------------------------------------------------------------- if (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 cc1 = max(1.1_dbl_kind/rncat,c1*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 ! 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 / hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo5(n) enddo elseif (ncat == 6) then ! 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 / hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo6(n) enddo elseif (ncat == 7) then ! 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 / hin_max(0) = c0 do n = 1, ncat hin_max(n) = wmo7(n) enddo else write (nu_diag,*) 'kcatbound=3 (WMO) must have ncat=5, 6 or 7' stop endif endif ! kcatbound if (my_task == master_task) then write (nu_diag,*) ' ' write (nu_diag,*) 'hin_max(n-1) < Cat n < hin_max(n)' do n = 1, ncat write (nu_diag,*) hin_max(n-1),' < Cat ',n, ' < ',hin_max(n) ! Write integer n to character string write (c_nc, '(i2)') n ! Write hin_max to character string write (c_hinmax1, '(f5.3)') hin_max(n-1) write (c_hinmax2, '(f5.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 (nu_diag,*) ' ' endif !----------------------------------------------------------------- ! vectors identifying first and last layer in each category !----------------------------------------------------------------- ilyr1(1) = 1 ! if nilyr = 4 ilyrn(1) = nilyr ! ilyr1 = { 1,5,9 } do n = 2, ncat ! ilyrn = { 4,8,12} etc ilyr1(n) = ilyrn(n-1) + 1 ilyrn(n) = ilyrn(n-1) + nilyr enddo slyr1(1) = 1 slyrn(1) = nslyr do n = 2, ncat slyr1(n) = slyrn(n-1) + 1 slyrn(n) = slyrn(n-1) + nslyr enddo end subroutine init_itd !======================================================================= !BOP ! ! !IROUTINE: aggregate - aggregate ice state variables ! ! !INTERFACE: ! subroutine aggregate (nx_block, ny_block, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & aice, trcr, & vice, vsno, & eice, esno, & aice0, tmask, & ntrcr, trcr_depend) ! ! !DESCRIPTION: ! ! Aggregate ice state variables over thickness categories. ! ! !REVISION HISTORY: ! ! authors: C. M. Bitz, UW ! W. H. Lipscomb, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ntrcr ! number of tracers in use real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & 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 (nx_block,ny_block,max_ntrcr,ncat), & intent(in) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), & intent(in) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(in) :: & esnon ! energy of melting for each snow layer (J/m^2) logical (kind=log_kind), dimension (nx_block,ny_block), & intent(in) :: & tmask ! land/boundary mask, thickness (T-cell) integer (kind=int_kind), dimension (max_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(out) :: & aice , & ! concentration of ice vice , & ! volume per unit area of ice (m) vsno , & ! volume per unit area of snow (m) eice , & ! energy of melt. of ice (J/m^2) esno , & ! energy of melt. of snow layer (J/m^2) aice0 ! concentration of open water real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr), & intent(out) :: & trcr ! ice tracers ! !EOP ! integer (kind=int_kind) :: & icells ! number of ocean/ice cells integer (kind=int_kind), dimension (nx_block*ny_block) :: & indxi, & ! compressed indices in i/j directions indxj integer (kind=int_kind) :: & i, j, k, n, it, & ij ! combined i/j horizontal index real (kind=dbl_kind), dimension (:,:), allocatable :: & atrcr ! sum of aicen*trcrn or vicen*trcrn or vsnon*trcrn !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- icells = 0 do j = 1, ny_block do i = 1, nx_block if (tmask(i,j)) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif ! tmask aice0(i,j) = c1 aice (i,j) = c0 vice (i,j) = c0 vsno (i,j) = c0 eice (i,j) = c0 esno (i,j) = c0 enddo enddo allocate (atrcr(icells,ntrcr)) !----------------------------------------------------------------- ! Aggregate !----------------------------------------------------------------- atrcr(:,:) = c0 do n = 1, ncat !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) aice(i,j) = aice(i,j) + aicen(i,j,n) vice(i,j) = vice(i,j) + vicen(i,j,n) vsno(i,j) = vsno(i,j) + vsnon(i,j,n) enddo ! ij 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, icells i = indxi(ij) j = indxj(ij) atrcr(ij,it) = atrcr(ij,it) & + trcrn(i,j,it,n)*aicen(i,j,n) enddo ! ij elseif (trcr_depend(it) == 1) then ! ice volume tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcr(ij,it) = atrcr(ij,it) & + trcrn(i,j,it,n)*vicen(i,j,n) enddo ! ij elseif (trcr_depend(it) ==2) then ! snow volume tracer !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) atrcr(ij,it) = atrcr(ij,it) & + trcrn(i,j,it,n)*vsnon(i,j,n) enddo ! ij endif ! trcr_depend enddo ! ntrcr do k = 1, nilyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) eice(i,j) = eice(i,j) + eicen(i,j,ilyr1(n)+k-1) enddo enddo ! nilyr do k = 1, nslyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) esno(i,j) = esno(i,j) + esnon(i,j,slyr1(n)+k-1) enddo enddo ! nslyr enddo ! ncat ! Open water fraction do ij = 1, icells i = indxi(ij) j = indxj(ij) aice0(i,j) = max (c1 - aice(i,j), c0) enddo ! ij ! Tracers call compute_tracers (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & atrcr(:,:), aice(:,:), & vice (:,:), vsno(:,:), & trcr(:,:,:)) deallocate (atrcr) end subroutine aggregate !======================================================================= !BOP ! ! !IROUTINE: aggregate_area - aggregate ice area ! ! !INTERFACE: ! subroutine aggregate_area (nx_block, ny_block, & aicen, aice, aice0) ! ! !DESCRIPTION: ! ! Aggregate ice area (but not other state variables) over thickness ! categories. ! ! !REVISION HISTORY: ! ! authors: William H. Lipscomb, LANL ! modified Jan 2004 by Clifford Chen, Fujitsu ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block ! block dimensions real (kind=dbl_kind), dimension (:,:,:), intent(in) :: & aicen ! concentration of ice real (kind=dbl_kind), dimension (:,:), intent(inout) :: & aice, & ! concentration of ice aice0 ! concentration of open water ! !EOP ! integer (kind=int_kind) :: i, j, n !----------------------------------------------------------------- ! Aggregate !----------------------------------------------------------------- aice(:,:) = c0 do n = 1, ncat do j = 1, ny_block do i = 1, nx_block aice(i,j) = aice(i,j) + aicen(i,j,n) enddo ! i enddo ! j enddo ! n do j = 1, ny_block do i = 1, nx_block ! open water fraction aice0(i,j) = max (c1 - aice(i,j), c0) enddo ! i enddo ! j end subroutine aggregate_area !======================================================================= !BOP ! ! !IROUTINE: rebin - rebins thicknesses into defined categories ! ! !INTERFACE: ! subroutine rebin (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & l_stop, & istop, jstop) ! ! !DESCRIPTION: ! ! Rebins thicknesses into defined categories ! ! !REVISION HISTORY: ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells , & ! number of grid cells with ice ntrcr ! number of tracers in use integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i/j indices integer (kind=int_kind), dimension (max_ntrcr), intent(in) :: & trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon 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,max_ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), & intent(inout) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(inout) :: & esnon ! energy of melting for each snow layer (J/m^2) 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 ! !EOP ! integer (kind=int_kind) :: & i,j , & ! horizontal indices n , & ! category index ij ! combined horizontal index logical (kind=log_kind) :: & shiftflag ! = .true. if ice must be shifted integer (kind=int_kind), dimension (icells,ncat) :: & donor ! donor category index real (kind=dbl_kind), dimension (icells,ncat) :: & daice , & ! ice area transferred dvice , & ! ice volume transferred hicen ! ice thickness for each cat (m) !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- l_stop = .false. istop = 0 jstop = 0 do n = 1, ncat do ij = 1, icells ! aice(i,j) > puny i = indxi(ij) j = indxj(ij) donor(ij,n) = 0 daice(ij,n) = c0 dvice(ij,n) = c0 !----------------------------------------------------------------- ! Compute ice thickness. !----------------------------------------------------------------- if (aicen(i,j,n) > puny) then hicen(ij,n) = vicen(i,j,n) / aicen(i,j,n) else hicen(ij,n) = c0 endif enddo ! ij enddo ! n !----------------------------------------------------------------- ! make sure thickness of cat 1 is at least hin_max(0) !----------------------------------------------------------------- do ij = 1, icells ! aice(i,j) > puny i = indxi(ij) j = indxj(ij) if (aicen(i,j,1) > puny) then if (hicen(ij,1) <= hin_max(0) .and. hin_max(0) > c0 ) then aicen(i,j,1) = vicen(i,j,1) / hin_max(0) hicen(ij,1) = hin_max(0) endif endif enddo ! ij !----------------------------------------------------------------- ! 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. do ij = 1, icells ! aice(i,j) > puny i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny .and. & hicen(ij,n) > hin_max(n)) then shiftflag = .true. donor(ij,n) = n daice(ij,n) = aicen(i,j,n) dvice(ij,n) = vicen(i,j,n) endif enddo ! ij if (shiftflag) then !----------------------------------------------------------------- ! shift ice between categories !----------------------------------------------------------------- call shift_ice (nx_block, ny_block, & indxi, indxj, & icells, & ntrcr, trcr_depend, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & hicen, donor, & daice, dvice, & l_stop, & istop, jstop) if (l_stop) return !----------------------------------------------------------------- ! reset shift parameters !----------------------------------------------------------------- do ij = 1, icells ! aice(i,j) > puny donor(ij,n) = 0 daice(ij,n) = c0 dvice(ij,n) = c0 enddo 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. do ij = 1, icells ! aice(i,j) > puny i = indxi(ij) j = indxj(ij) if (aicen(i,j,n+1) > puny .and. & hicen(ij,n+1) <= hin_max(n)) then shiftflag = .true. donor(ij,n) = n+1 daice(ij,n) = aicen(i,j,n+1) dvice(ij,n) = vicen(i,j,n+1) endif enddo ! ij if (shiftflag) then !----------------------------------------------------------------- ! shift ice between categories !----------------------------------------------------------------- call shift_ice (nx_block, ny_block, & indxi, indxj, & icells, & ntrcr, trcr_depend, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & hicen, donor, & daice, dvice, & l_stop, & istop, jstop) if (l_stop) return !----------------------------------------------------------------- ! reset shift parameters !----------------------------------------------------------------- do ij = 1, icells ! aice(i,j) > puny donor(ij,n) = 0 daice(ij,n) = c0 dvice(ij,n) = c0 enddo endif ! shiftflag enddo ! n end subroutine rebin !======================================================================= !BOP ! ! !IROUTINE: reduce_area - reduce area when ice melts for special case ncat=1 ! ! !INTERFACE: ! subroutine reduce_area (nx_block, ny_block, & ilo, ihi, jlo, jhi, & tmask, & aicen, vicen, & aicen_init,vicen_init) ! ! !DESCRIPTION: ! ! 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 ! ! !REVISION HISTORY: ! ! authors: C. M. Bitz, UW ! modified by: Elizabeth C. Hunke, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ilo,ihi,jlo,jhi ! beginning and end of physical domain logical (kind=log_kind), dimension (nx_block,ny_block), & intent(in) :: & tmask ! land/boundary mask, thickness (T-cell) real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & aicen , & ! concentration of ice vicen ! volume per unit area of ice (m) real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & aicen_init, & ! old ice area for category 1 (m) vicen_init ! old ice volume for category 1 (m) ! !EOP ! integer (kind=int_kind) :: & i, j ! horizontal indices real (kind=dbl_kind) :: & hi0 , & ! initial hi hi1 , & ! current hi dhi ! hi1 - hi0 do j = jlo, jhi do i = ilo, ihi if (tmask(i,j)) then hi0 = c0 if (aicen_init(i,j) > c0) & hi0 = vicen_init(i,j) / aicen_init(i,j) hi1 = c0 if (aicen(i,j) > c0) & hi1 = vicen(i,j) / aicen(i,j) ! make sure thickness of cat 1 is at least hin_max(0) if (hi1 <= hin_max(0) .and. hin_max(0) > c0 ) then aicen(i,j) = vicen(i,j) / hin_max(0) hi1 = hin_max(0) endif if (aicen(i,j) > c0) then dhi = hi1 - hi0 if (dhi < c0) then hi1 = vicen(i,j) / aicen(i,j) aicen(i,j) = c2 * vicen(i,j) / (hi1 + hi0) endif endif endif ! tmask enddo ! i enddo ! j end subroutine reduce_area !======================================================================= !BOP ! ! !IROUTINE: shift_ice - shift ice across category boundaries ! ! !INTERFACE: ! subroutine shift_ice (nx_block, ny_block, & indxi, indxj, & icells, & ntrcr, trcr_depend, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & hicen, donor, & daice, dvice, & l_stop, & istop, jstop) ! ! !DESCRIPTION: ! ! Shift ice across category boundaries, conserving area, volume, and ! energy. ! ! !REVISION HISTORY: ! ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL ! ! !USES: ! ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells , & ! number of ocean/ice cells ntrcr ! number of tracers in use integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi , & ! compressed indices in i/j directions indxj integer (kind=int_kind), dimension (max_ntrcr), intent(in) :: & trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon 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,max_ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), & intent(inout) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(inout) :: & esnon ! energy of melting for each snow layer (J/m^2) ! 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(icells,ncat), & intent(in) :: & donor ! donor category index real (kind=dbl_kind), dimension(icells,ncat), & intent(inout) :: & daice , & ! ice area transferred across boundary dvice , & ! ice volume transferred across boundary hicen ! ice thickness for each cat (m) 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 ! !EOP ! integer (kind=int_kind) :: & i, j, m , & ! horizontal indices n , & ! thickness category index nr , & ! receiver category nd , & ! donor category k , & ! ice layer index it , & ! tracer index ilo,ihi,jlo,jhi ! beginning and end of physical domain real (kind=dbl_kind), dimension(icells,max_ntrcr,ncat) :: & atrcrn ! aicen*trcrn real (kind=dbl_kind) :: & dvsnow , & ! snow volume transferred desnow , & ! snow energy transferred deice , & ! ice energy transferred datrcr ! aicen*train transferred integer (kind=int_kind), dimension (icells) :: & indxii , & ! compressed indices for i/j directions indxjj , & indxij integer (kind=int_kind) :: & ishift , & ! number of cells with ice to transfer ij ! combined i/j horizontal index 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), dimension (nx_block,ny_block) :: & worka, & workb !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- l_stop = .false. istop = 0 jstop = 0 worka(:,:) = c0 workb(:,:) = c0 !----------------------------------------------------------------- ! 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 endif enddo enddo !----------------------------------------------------------------- ! 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. do ij = 1, icells i = indxi(ij) j = indxj(ij) if (donor(ij,n) > 0) then nd = donor(ij,n) if (daice(ij,n) < c0) then if (daice(ij,n) > -puny*aicen(i,j,nd)) then daice(ij,n) = c0 ! shift no ice dvice(ij,n) = c0 else daice_negative = .true. endif endif if (dvice(ij,n) < c0) then if (dvice(ij,n) > -puny*vicen(i,j,nd)) then daice(ij,n) = c0 ! shift no ice dvice(ij,n) = c0 else dvice_negative = .true. endif endif if (daice(ij,n) > aicen(i,j,nd)*(c1-puny)) then if (daice(ij,n) < aicen(i,j,nd)*(c1+puny)) then daice(ij,n) = aicen(i,j,nd) dvice(ij,n) = vicen(i,j,nd) else daice_greater_aicen = .true. endif endif if (dvice(ij,n) > vicen(i,j,nd)*(c1-puny)) then if (dvice(ij,n) < vicen(i,j,nd)*(c1+puny)) then daice(ij,n) = aicen(i,j,nd) dvice(ij,n) = vicen(i,j,nd) else dvice_greater_vicen = .true. endif endif endif ! donor > 0 enddo ! ij !----------------------------------------------------------------- ! error messages !----------------------------------------------------------------- if (daice_negative) then do ij = 1, icells i = indxi(ij) j = indxj(ij) if (donor(ij,n) > 0 .and. & daice(ij,n) <= -puny*aicen(i,j,nd)) then write(nu_diag,*) ' ' write(nu_diag,*) 'shift_ice: negative daice' write(nu_diag,*) 'i, j:', i, j write(nu_diag,*) 'boundary, donor cat:', n, nd write(nu_diag,*) 'daice =', daice(ij,n) write(nu_diag,*) 'dvice =', dvice(ij,n) l_stop = .true. istop = i jstop = j endif enddo endif if (l_stop) return if (dvice_negative) then do ij = 1, icells i = indxi(ij) j = indxj(ij) if (donor(ij,n) > 0 .and. & dvice(ij,n) <= -puny*vicen(i,j,nd)) then write(nu_diag,*) ' ' write(nu_diag,*) 'shift_ice: negative dvice' write(nu_diag,*) 'i, j:', i, j write(nu_diag,*) 'boundary, donor cat:', n, nd write(nu_diag,*) 'daice =', daice(ij,n) write(nu_diag,*) 'dvice =', dvice(ij,n) l_stop = .true. istop = i jstop = j endif enddo endif if (l_stop) return if (daice_greater_aicen) then do ij = 1, icells i = indxi(ij) j = indxj(ij) if (donor(ij,n) > 0) then nd = donor(ij,n) if (daice(ij,n) >= aicen(i,j,nd)*(c1+puny)) then write(nu_diag,*) ' ' write(nu_diag,*) 'shift_ice: daice > aicen' write(nu_diag,*) 'i, j:', i, j write(nu_diag,*) 'boundary, donor cat:', n, nd write(nu_diag,*) 'daice =', daice(ij,n) write(nu_diag,*) 'aicen =', aicen(i,j,nd) l_stop = .true. istop = i jstop = j endif endif enddo endif if (l_stop) return if (dvice_greater_vicen) then do ij = 1, icells i = indxi(ij) j = indxj(ij) if (donor(ij,n) > 0) then nd = donor(ij,n) if (dvice(ij,n) >= vicen(i,j,nd)*(c1+puny)) then write(nu_diag,*) ' ' write(nu_diag,*) 'shift_ice: dvice > vicen' write(nu_diag,*) 'i, j:', i, j write(nu_diag,*) 'boundary, donor cat:', n, nd write(nu_diag,*) 'dvice =', dvice(ij,n) write(nu_diag,*) 'vicen =', vicen(i,j,nd) l_stop = .true. istop = i jstop = j endif endif enddo endif if (l_stop) return !----------------------------------------------------------------- ! transfer volume and energy between categories !----------------------------------------------------------------- ishift = 0 do ij = 1, icells i = indxi(ij) j = indxj(ij) if (daice(ij,n) > c0) then ! daice(n) can be < puny ishift = ishift + 1 indxii(ishift) = i indxjj(ishift) = j indxij(ishift) = ij endif ! tmask enddo !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, ishift i = indxii(ij) j = indxjj(ij) m = indxij(ij) nd = donor(m,n) worka(i,j) = dvice(m,n) / vicen(i,j,nd) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif aicen(i,j,nd) = aicen(i,j,nd) - daice(m,n) aicen(i,j,nr) = aicen(i,j,nr) + daice(m,n) vicen(i,j,nd) = vicen(i,j,nd) - dvice(m,n) vicen(i,j,nr) = vicen(i,j,nr) + dvice(m,n) dvsnow = vsnon(i,j,nd) * worka(i,j) vsnon(i,j,nd) = vsnon(i,j,nd) - dvsnow vsnon(i,j,nr) = vsnon(i,j,nr) + dvsnow workb(i,j) = dvsnow enddo ! ij do it = 1, ntrcr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, ishift i = indxii(ij) j = indxjj(ij) m = indxij(ij) nd = donor(m,n) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif if (trcr_depend(it) == 0) then datrcr = daice(m,n)*trcrn(i,j,it,nd) elseif (trcr_depend(it) == 1) then datrcr = dvice(m,n)*trcrn(i,j,it,nd) elseif (trcr_depend(it) == 2) then datrcr = workb(i,j) *trcrn(i,j,it,nd) endif atrcrn(m,it,nd) = atrcrn(m,it,nd) - datrcr atrcrn(m,it,nr) = atrcrn(m,it,nr) + datrcr enddo ! ij enddo ! ntrcr do k = 1, nilyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, ishift i = indxii(ij) j = indxjj(ij) m = indxij(ij) nd = donor(m,n) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif deice = eicen(i,j,ilyr1(nd)+k-1) * worka(i,j) eicen(i,j,ilyr1(nd)+k-1) = & eicen(i,j,ilyr1(nd)+k-1) - deice eicen(i,j,ilyr1(nr)+k-1) = & eicen(i,j,ilyr1(nr)+k-1) + deice enddo ! ij enddo ! nilyr do k = 1, nslyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, ishift i = indxii(ij) j = indxjj(ij) m = indxij(ij) nd = donor(m,n) if (nd == n) then nr = nd+1 else ! nd = n+1 nr = n endif desnow = esnon(i,j,slyr1(nd)+k-1) * worka(i,j) esnon(i,j,slyr1(nd)+k-1) = & esnon(i,j,slyr1(nd)+k-1) - desnow esnon(i,j,slyr1(nr)+k-1) = & esnon(i,j,slyr1(nr)+k-1) + desnow enddo ! ij enddo ! nslyr enddo ! boundaries, 1 to ncat-1 !----------------------------------------------------------------- ! Update ice thickness and tracers !----------------------------------------------------------------- do n = 1, ncat do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j,n) > puny) then hicen(ij,n) = vicen (i,j,n) / aicen(i,j,n) else hicen(ij,n) = c0 endif enddo call compute_tracers (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & atrcrn(:,:,n), aicen(:,:, n), & vicen (:,:, n), vsnon(:,:, n), & trcrn(:,:,:,n)) enddo ! ncat end subroutine shift_ice !======================================================================= !BOP ! ! !IROUTINE: column_sum - sum field over all ice categories ! ! !INTERFACE: ! subroutine column_sum (nx_block, ny_block, & icells, indxi, indxj, & nsum, & xin, xout) ! ! !DESCRIPTION: ! ! For each grid cell, sum field over all ice categories. ! ! !REVISION HISTORY: ! ! author: William H. Lipscomb, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions nsum , & ! number of categories/layers icells ! number of ice/ocean grid cells integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i/j indices real (kind=dbl_kind), dimension (nx_block,ny_block,nsum), & intent(in) :: & xin ! input field real (kind=dbl_kind), dimension (icells), intent(out) :: & xout ! output field ! !EOP ! integer (kind=int_kind) :: & i, j, ij , & ! horizontal indices n ! category/layer index do ij = 1, icells xout(ij) = c0 enddo do n = 1, nsum do ij = 1, icells i = indxi(ij) j = indxj(ij) xout(ij) = xout(ij) + xin(i,j,n) enddo ! ij enddo ! n end subroutine column_sum !======================================================================= !BOP ! ! !IROUTINE: column_conservation_check ! ! !INTERFACE: ! subroutine column_conservation_check (nx_block, ny_block, & icells, indxi, indxj, & fieldid, & x1, x2, & max_err, l_stop, & istop, jstop) ! ! !DESCRIPTION: ! ! For each physical grid cell, check that initial and final values ! of a conserved field are equal to within a small value. ! ! !REVISION HISTORY: ! ! author: William H. Lipscomb, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of ice/ocean grid cells integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i/j indices real (kind=dbl_kind), dimension(icells), 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 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 ! !EOP ! integer (kind=int_kind) :: & ij ! horizontal indices do ij = 1, icells if (abs (x2(ij)-x1(ij)) > max_err) then l_stop = .true. istop = indxi(ij) jstop = indxj(ij) write (nu_diag,*) ' ' write (nu_diag,*) 'Conservation error: ', trim(fieldid) write (nu_diag,*) 'i, j =', istop, jstop write (nu_diag,*) 'Initial value =', x1(ij) write (nu_diag,*) 'Final value =', x2(ij) write (nu_diag,*) 'Difference =', x2(ij) - x1(ij) endif enddo end subroutine column_conservation_check !======================================================================= !BOP ! ! !IROUTINE: compute_tracers - compute tracer fields ! ! !INTERFACE: ! subroutine compute_tracers (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & atrcrn, aicen, & vicen, vsnon, & trcrn) ! ! !DESCRIPTION: ! ! Compute tracer fields. ! Given atrcrn = aicen*trcrn (or vicen*trcrn, vsnon*trcrn), compute trcrn. ! ! !REVISION HISTORY: ! ! author: William H. Lipscomb, LANL ! ! !USES: ! use ice_state, only: nt_Tsfc ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells , & ! number of ice/ocean grid cells ntrcr ! number of tracers in use integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i/j indices integer (kind=int_kind), dimension (max_ntrcr), intent(in) :: & trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon real (kind=dbl_kind), dimension (icells,max_ntrcr), & intent(in) :: & atrcrn ! aicen*trcrn or vicen*trcrn or vsnon*trcrn real (kind=dbl_kind), dimension (nx_block,ny_block), & 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 (nx_block,ny_block,max_ntrcr), & intent(out) :: & trcrn ! ice tracers ! !EOP ! integer (kind=int_kind) :: & i, j, it, ij ! counting indices trcrn(:,:,:) = c0 !----------------------------------------------------------------- ! Compute new tracers !----------------------------------------------------------------- do it = 1, ntrcr if (it == nt_Tsfc) then ! surface temperature do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j) > puny) then trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j) else trcrn(i,j,it) = Tocnfrz endif enddo elseif (trcr_depend(it) == 0) then ! ice area tracers do ij = 1, icells i = indxi(ij) j = indxj(ij) if (aicen(i,j) > puny) then trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j) else trcrn(i,j,it) = c0 endif enddo elseif (trcr_depend(it) == 1) then ! ice volume tracers do ij = 1, icells i = indxi(ij) j = indxj(ij) if (vicen(i,j) > puny) then trcrn(i,j,it) = atrcrn(ij,it) / vicen(i,j) else trcrn(i,j,it) = c0 endif enddo elseif (trcr_depend(it) == 2) then ! snow volume tracers do ij = 1, icells i = indxi(ij) j = indxj(ij) if (vsnon(i,j) > puny) then trcrn(i,j,it) = atrcrn(ij,it) / vsnon(i,j) else trcrn(i,j,it) = c0 endif enddo endif ! trcr_depend enddo ! ntrcr end subroutine compute_tracers !======================================================================= !BOP ! ! !IROUTINE: cleanup_itd - rebin if needed, eliminate small ice areas, ! and aggregate over categories ! ! !INTERFACE: ! subroutine cleanup_itd (nx_block, ny_block, & ilo, ihi, jlo, jhi, & dt, ntrcr, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & aice0, aice, & trcr_depend, fresh, & fsalt, fhocn, & fsoot, tr_aero, & heat_capacity, l_stop, & istop, jstop, & limit_aice_in) ! ! !DESCRIPTION: ! ! 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). ! ! !REVISION HISTORY: ! ! author: William H. Lipscomb, LANL ! ! !USES: ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ntrcr ! number of tracers in use real (kind=dbl_kind), intent(in) :: & dt ! time step 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,max_ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), & intent(inout) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(inout) :: & esnon ! energy of melting for each snow layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & aice , & ! total ice concentration aice0 ! concentration of open water integer (kind=int_kind), dimension(max_ntrcr), intent(in) :: & trcr_depend ! tracer dependency information logical (kind=log_kind), intent(in) :: & tr_aero, & heat_capacity ! if false, ice and snow have zero heat capacity 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 ! ice-ocean fluxes (required for strict conservation) real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout), optional :: & 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) real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx), & intent(inout), optional :: & fsoot ! soot 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 ! !EOP ! integer (kind=int_kind) :: & i, j , & ! horizontal indices n , & ! category index icells ! number of grid cells with ice integer (kind=int_kind), dimension (nx_block*ny_block) :: & indxi, indxj ! compressed i/j indices real (kind=dbl_kind), dimension (nx_block,ny_block) :: & dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx) :: & dfsoot ! zapped soot flux (kg/m^2/s) logical (kind=log_kind) :: & limit_aice ! if true, check for aice out of bounds !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- if (present(limit_aice_in)) then limit_aice = limit_aice_in else limit_aice = .true. endif l_stop = .false. istop = 0 jstop = 0 !----------------------------------------------------------------- ! Compute total ice area. !----------------------------------------------------------------- call aggregate_area (nx_block, ny_block, & aicen(:,:,:), & aice, aice0) if (limit_aice) then ! check for aice out of bounds do j = jlo,jhi do i = ilo,ihi if (aice(i,j) > c1+puny .or. aice(i,j) < -puny) then l_stop = .true. istop = i jstop = j endif enddo enddo if (l_stop) then ! area out of bounds i = istop j = jstop write(nu_diag,*) ' ' write(nu_diag,*) 'aggregate ice area out of bounds' write(nu_diag,*) 'my_task, i, j, aice:', & my_task, i, j, aice(i,j) do n = 1, ncat write(nu_diag,*) 'n, aicen:', n, aicen(i,j,n) enddo return endif ! l_stop endif ! limit_aice !----------------------------------------------------------------- ! Identify grid cells with ice. !----------------------------------------------------------------- icells = 0 do j = jlo,jhi do i = ilo,ihi if (aice(i,j) > puny) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif enddo enddo !----------------------------------------------------------------- ! 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 (nx_block, ny_block, & icells, indxi, indxj, & ntrcr, trcr_depend, & aicen(:,:,:), trcrn(:,:,:,:), & vicen(:,:,:), vsnon(:,:,:), & eicen(:,:,:), esnon(:,:,:), & l_stop, & istop, jstop) if (l_stop) return !----------------------------------------------------------------- ! Zero out ice categories with very small areas. !----------------------------------------------------------------- if (limit_aice) then call zap_small_areas (nx_block, ny_block, & ilo, ihi, jlo, jhi, & dt, ntrcr, & aice, aice0, & aicen(:,:,:), trcrn(:,:,:,:), & vicen(:,:,:), vsnon(:,:,:), & eicen(:,:,:), esnon(:,:,:), & dfresh, dfsalt, & dfhocn, dfsoot, & tr_aero, & l_stop, & istop, jstop) if (l_stop) return endif ! l_limit_aice !------------------------------------------------------------------- ! Update ice-ocean fluxes for strict conservation !------------------------------------------------------------------- if (present(fresh)) & fresh (:,:) = fresh(:,:) + dfresh(:,:) if (present(fsalt)) & fsalt (:,:) = fsalt(:,:) + dfsalt(:,:) if (present(fhocn)) & fhocn (:,:) = fhocn(:,:) + dfhocn(:,:) if (present(fsoot)) & fsoot (:,:,:) = fsoot(:,:,:) + dfsoot(:,:,:) !---------------------------------------------------------------- ! If using zero-layer model (no heat capacity), check that the ! energy of snow and ice is correct. !---------------------------------------------------------------- if (.not. heat_capacity) then call zerolayer_check(nx_block, ny_block, & icells, indxi, indxj, & aicen, & vicen, vsnon, & eicen, esnon, & l_stop, & istop, jstop) endif end subroutine cleanup_itd !======================================================================= !BOP ! ! !IROUTINE: zap_small_areas - eliminate very small ice areas ! ! !INTERFACE: ! subroutine zap_small_areas (nx_block, ny_block, & ilo, ihi, jlo, jhi, & dt, ntrcr, & aice, aice0, & aicen, trcrn, & vicen, vsnon, & eicen, esnon, & dfresh, dfsalt, & dfhocn, dfsoot, & tr_aero, & l_stop, & istop, jstop) ! ! !DESCRIPTION: ! ! For each ice category in each grid cell, remove ice if the fractional ! area is less than puny. ! ! !REVISION HISTORY: ! ! author: William H. Lipscomb, LANL ! ! !USES: ! use ice_state, only: nt_Tsfc, nt_aero ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ilo,ihi,jlo,jhi , & ! beginning and end of physical domain ntrcr ! number of tracers in use real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & aice , & ! total ice concentration 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,ntilyr), & intent(inout) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(inout) :: & esnon ! energy of melting for each snow layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), & intent(inout) :: & trcrn ! ice tracers real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(out) :: & dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx), & intent(out) :: & dfsoot ! zapped soot flux (kg/m^2/s) logical (kind=log_kind), intent(in) :: & tr_aero 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 ! !EOP ! integer (kind=int_kind) :: & i,j, n, k, it , & ! counting indices icells , & ! number of cells with ice to zap ij ! combined i/j horizontal index integer (kind=int_kind), dimension (nx_block*ny_block) :: & indxi , & ! compressed indices for i/j directions indxj real (kind=dbl_kind) :: xtmp ! temporary variable !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- l_stop = .false. istop = 0 jstop = 0 dfresh(:,:) = c0 dfsalt(:,:) = c0 dfhocn(:,:) = c0 dfsoot(:,:,:) = c0 !----------------------------------------------------------------- ! Zap categories with very small areas. !----------------------------------------------------------------- do n = 1, ncat !----------------------------------------------------------------- ! Count categories to be zapped. !----------------------------------------------------------------- icells = 0 do j = jlo, jhi do i = ilo, ihi if (aicen(i,j,n) < -puny) then write (nu_diag,*) 'Zap ice: negative ice area' write (nu_diag,*) 'i, j, n, aicen =', & i, j, n, aicen(i,j,n) l_stop = .true. istop = i jstop = j return elseif ((aicen(i,j,n) >= -puny .and. aicen(i,j,n) < c0) .or. & (aicen(i,j,n) > c0 .and. aicen(i,j,n) <= puny)) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif enddo enddo !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- do k = 1, nilyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = eicen(i,j,ilyr1(n)+k-1) / dt ! < 0 dfhocn(i,j) = dfhocn(i,j) + xtmp eicen(i,j,ilyr1(n)+k-1) = c0 enddo ! ij enddo ! k !----------------------------------------------------------------- ! Zap snow energy and use ocean heat to melt snow !----------------------------------------------------------------- do k = 1, nslyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = esnon(i,j,slyr1(n)+k-1) / dt ! < 0 dfhocn(i,j) = dfhocn(i,j) + xtmp esnon(i,j,slyr1(n)+k-1) = c0 enddo ! ij enddo ! k !----------------------------------------------------------------- ! Zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt dfresh(i,j) = dfresh(i,j) + xtmp xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 / dt dfsalt(i,j) = dfsalt(i,j) + xtmp aice0(i,j) = aice0(i,j) + aicen(i,j,n) aicen(i,j,n) = c0 vicen(i,j,n) = c0 vsnon(i,j,n) = c0 trcrn(i,j,nt_Tsfc,n) = Tocnfrz enddo ! ij if (tr_aero) then !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) do it=1,n_aero xtmp & = (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) & +trcrn(i,j,nt_aero+1+4*(it-1),n)) & + vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) & +trcrn(i,j,nt_aero+3+4*(it-1),n))) & / dt dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp enddo ! n enddo ! ij endif !----------------------------------------------------------------- ! Zap tracers !----------------------------------------------------------------- if (ntrcr >= 2) then do it = 1, ntrcr ! this assumes nt_Tsfc = 1 do ij = 1, icells i = indxi(ij) j = indxj(ij) trcrn(i,j,it,n) = c0 enddo enddo endif enddo ! n !----------------------------------------------------------------- ! Count cells with excess ice (aice > c1) due to roundoff errors. ! Zap a little ice in each category so that aice = c1. !----------------------------------------------------------------- icells = 0 do j = jlo, jhi do i = ilo, ihi if (aice(i,j) > (c1+puny)) then write (nu_diag,*) 'Zap ice: excess ice area' write (nu_diag,*) 'i, j, aice =', & i, j, aice(i,j) l_stop = .true. istop = i jstop = j return elseif (aice(i,j) > c1 .and. aice(i,j) < (c1+puny)) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif enddo enddo do n = 1, ncat !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- do k = 1, nilyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = eicen(i,j,ilyr1(n)+k-1) & * (aice(i,j)-c1)/aice(i,j) / dt ! < 0 dfhocn(i,j) = dfhocn(i,j) + xtmp eicen(i,j,ilyr1(n)+k-1) = eicen(i,j,ilyr1(n)+k-1) & * (c1/aice(i,j)) enddo ! ij enddo ! k !----------------------------------------------------------------- ! Zap snow energy and use ocean heat to melt snow !----------------------------------------------------------------- do k = 1, nslyr !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = esnon(i,j,slyr1(n)+k-1) & * (aice(i,j)-c1)/aice(i,j) / dt ! < 0 dfhocn(i,j) = dfhocn(i,j) + xtmp esnon(i,j,slyr1(n)+k-1) = esnon(i,j,slyr1(n)+k-1) & *(c1/aice(i,j)) enddo ! ij enddo ! k !----------------------------------------------------------------- ! Zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) & * (aice(i,j)-c1)/aice(i,j) / dt dfresh(i,j) = dfresh(i,j) + xtmp xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 & * (aice(i,j)-c1)/aice(i,j) / dt dfsalt(i,j) = dfsalt(i,j) + xtmp aicen(i,j,n) = aicen(i,j,n) * (c1/aice(i,j)) vicen(i,j,n) = vicen(i,j,n) * (c1/aice(i,j)) vsnon(i,j,n) = vsnon(i,j,n) * (c1/aice(i,j)) enddo ! ij ! Note: Tracers are unchanged. !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu if (tr_aero) then do ij = 1, icells i = indxi(ij) j = indxj(ij) do it=1,n_aero xtmp & = (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) & +trcrn(i,j,nt_aero+1+4*(it-1),n)) & + vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) & +trcrn(i,j,nt_aero+3+4*(it-1),n))) & * (aice(i,j)-c1)/aice(i,j) / dt dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp enddo ! n enddo ! ij endif enddo ! n !----------------------------------------------------------------- ! Correct aice !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) aice(i,j) = c1 aice0(i,j) = c0 enddo end subroutine zap_small_areas !======================================================================= !BOP ! ! !IROUTINE: zerolayer_check - check that snow and ice energy is ! correct when using zero layer thermodynamics ! ! !INTERFACE: ! subroutine zerolayer_check (nx_block, ny_block, & icells, indxi, indxj, & aicen, & vicen, vsnon, & eicen, esnon, & l_stop, & istop, jstop) ! ! !DESCRIPTION: ! ! 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. ! ! !REVISION HISTORY: ! ! author: Alison McLaren, Met Office ! ! !USES: ! ! ! !INPUT/OUTPUT PARAMETERS: ! integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of grid cells with ice integer (kind=int_kind), dimension (nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i/j indices 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,ntilyr), & intent(in) :: & eicen ! energy of melting for each ice layer (J/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), & intent(in) :: & esnon ! energy of melting for each snow layer (J/m^2) 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 ! !EOP ! integer (kind=int_kind) :: & i, j , & ! horizontal indices n , & ! category index ij ! combined horizontal index real (kind=dbl_kind), parameter :: & max_error = puny*Lfresh*rhos ! max error in zero layer energy check ! (so max volume error = puny) logical (kind=log_kind) :: & ice_energy_correct , & ! zero layer ice energy check snow_energy_correct ! zero layer snow energy check real (kind=dbl_kind), dimension (nx_block,ny_block) :: & worka, & workb !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- l_stop = .false. istop = 0 jstop = 0 worka(:,:) = c0 workb(:,:) = c0 !---------------------------------------------------------------- ! 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. do n=1,ncat do ij=1,icells i=indxi(ij) j=indxj(ij) worka(i,j) = eicen(i,j,n) + rhoi * Lfresh * vicen(i,j,n) workb(i,j) = esnon(i,j,n) + rhos * Lfresh * vsnon(i,j,n) if(abs(worka(i,j)) > max_error) then ice_energy_correct = .false. endif if(abs(workb(i,j)) > max_error) then snow_energy_correct = .false. endif enddo !---------------------------------------------------------------- ! If there is a problem, abort with error message !---------------------------------------------------------------- if (.not. ice_energy_correct) then do ij=1,icells i=indxi(ij) j=indxj(ij) if(abs(worka(i,j)) > max_error) then write(nu_diag,*) ' ' write(nu_diag,*) & 'zerolayer check - wrong ice energy' write(nu_diag,*) 'i, j, n:', i,j,n write(nu_diag,*) 'eicen =', eicen(i,j,n) write(nu_diag,*) 'error=', worka(i,j) write(nu_diag,*) 'vicen =', vicen(i,j,n) write(nu_diag,*) 'aicen =', aicen(i,j,n) l_stop = .true. istop = i jstop = j endif enddo endif if (l_stop) return if (.not. snow_energy_correct) then do ij=1,icells i=indxi(ij) j=indxj(ij) if(abs(workb(i,j)) > max_error) then write(nu_diag,*) ' ' write(nu_diag,*) & 'zerolayer_check - wrong snow energy' write(nu_diag,*) 'i, j, n:', i,j,n write(nu_diag,*) 'esnon =', esnon(i,j,n) write(nu_diag,*) 'error=', workb(i,j) write(nu_diag,*) 'vsnon =', vsnon(i,j,n) write(nu_diag,*) 'aicen =', aicen(i,j,n) l_stop = .true. istop = i jstop = j return endif enddo endif enddo ! ncat end subroutine zerolayer_check !======================================================================= end module ice_itd !=======================================================================