! SVN:$Id: ice_atmo.F90 1623 2017-05-15 00:18:33Z tcraig@ucar.edu $ !======================================================================= ! Atmospheric boundary interface (stability based flux calculations) ! author: Elizabeth C. Hunke, LANL ! ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb ! 2004: Block structure added by William Lipscomb ! 2006: Converted to free source form (F90) by Elizabeth Hunke ! 2013: Form drag routine added (neutral_drag_coeffs) by David Schroeder ! 2014: Adjusted form drag and added high frequency coupling by Andrew Roberts module ice_atmo use ice_kinds_mod use ice_blocks, only: nx_block, ny_block use ice_constants use ice_domain_size, only: max_blocks implicit none save private public :: atmo_boundary_layer, atmo_boundary_const, neutral_drag_coeffs character (len=char_len), public :: & atmbndy ! atmo boundary method, 'default' ('ccsm3') or 'constant' logical (kind=log_kind), public :: & calc_strair, & ! if true, calculate wind stress components formdrag, & ! if true, calculate form drag highfreq ! if true, use high frequency coupling integer (kind=int_kind), public :: & natmiter ! number of iterations for boundary layer calculations real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: & Cdn_atm, & ! neutral atm drag coefficient at zrefa Cdn_ocn, & ! neutral ocn drag coefficient at zrefo Cd_atm, & ! atm drag coefficient at zlvla Cd_ocn, & ! ocn drag coefficient at zlvlo hfreebd, & ! formdrag freeboard (m) hdraft, & ! formdrag draft of ice + snow column hridge, & ! formdrag ridge height distrdg, & ! formdrag distance between ridges hkeel, & ! formdrag keel depth dkeel, & ! formdrag distance between keels lfloe, & ! formdrag floe length dfloe, & ! formdrag distance between floes Cdn_atm_skin, & ! formdrag neutral skin drag coefficient at zrefa Cdn_atm_floe, & ! formdrag neutral floe edge drag coefficient at zrefa Cdn_atm_pond, & ! formdrag neutral pond edge drag coefficient at zrefa Cdn_atm_rdg, & ! formdrag neutral ridge drag coefficient at zrefa Cdn_ocn_skin, & ! formdrag neutral skin drag coefficient at zrefo Cdn_ocn_floe, & ! formdrag neutral floe edge drag coefficient at zrefo Cdn_ocn_keel, & ! formdrag neutral keel drag coefficient at zrefo Cdn_atm_ocn, & ! ratio neutral atm drag at zrefa / ocn drag at zrefo Cd_atm_ocn ! ratio atm drag at zlvla / ocn drag at zlvlo !======================================================================= contains !======================================================================= ! Compute coefficients for atm/ice fluxes, stress, and reference ! temperature and humidity. NOTE: \\ ! (1) All fluxes are positive downward, \\ ! (2) Here, tstar = (WT)/U*, and qstar = (WQ)/U*, \\ ! (3) wind speeds should all be above a minimum speed (eg. 1.0 m/s). \\ ! (4) highfreq is from Roberts et al. (2015), Annals, doi:10.3189/2015AoG69A760 ! ! ASSUME: ! The saturation humidity of air at T(K): qsat(T) (kg/m**3) ! ! Code originally based on CSM1 subroutine atmo_boundary_layer (nx_block, ny_block, & sfctype, icells, & indxi, indxj, & Tsf, potT, & uatm, vatm, & wind, zlvla, & Qa, rhoa, & strx, stry, & Tref, Qref, & delt, delq, & lhcoef, shcoef, & Cdn_atm, Cd_atm_n, & uice, vice, & Uref ) use ice_fileunits, only: nu_diag use ice_communicate, only: my_task, master_task use ice_exit, only: abort_ice integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells that require atmo fluxes integer (kind=int_kind), dimension(nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i and j indices character (len=3), intent(in) :: & sfctype ! ice or ocean real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & Tsf , & ! surface temperature of ice or ocean potT , & ! air potential temperature (K) uatm , & ! x-direction wind speed (m/s) vatm , & ! y-direction wind speed (m/s) wind , & ! wind speed (m/s) zlvla , & ! atm level height of lowest model layer (m) Qa , & ! specific humidity (kg/kg) rhoa ! air density (kg/m^3) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & Cdn_atm ! neutral atm drag coefficient real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & Cd_atm_n ! atm drag coefficient real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(inout) :: & strx , & ! x surface stress (N) stry ! y surface stress (N) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & Tref , & ! reference height temperature (K) Qref , & ! reference height specific humidity (kg/kg) delt , & ! potential T difference (K) delq , & ! humidity difference (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat real (kind=dbl_kind), dimension (nx_block,ny_block), optional, intent(out) :: & Uref ! reference height wind speed (m/s) real (kind=dbl_kind), dimension (nx_block,ny_block), optional, intent(in) :: & uice , & ! x-direction ice speed (m/s) vice ! y-direction ice speed (m/s) ! local variables logical (kind=log_kind), save :: & firstpass=.true. ! high frequency checks on first pass integer (kind=int_kind) :: & k , & ! iteration index i, j , & ! horizontal indices ij ! combined ij index real (kind=dbl_kind) :: & TsfK , & ! surface temperature in Kelvin (K) xqq , & ! temporary variable psimh , & ! stability function at zlvla (momentum) tau , & ! stress at zlvla fac , & ! interpolation factor al2 , & ! ln(z10 /zTrf) psix2 , & ! stability function at zTrf (heat and water) psimhs, & ! stable profile ssq , & ! sat surface humidity (kg/kg) qqq , & ! for qsat, dqsfcdt TTT , & ! for qsat, dqsfcdt qsat , & ! the saturation humidity of air (kg/m^3) Lheat , & ! Lvap or Lsub, depending on surface type umin ! minimum wind speed (m/s) real (kind=dbl_kind), dimension (icells) :: & ustar , & ! ustar (m/s) tstar , & ! tstar qstar , & ! qstar rdn , & ! sqrt of neutral exchange coefficient (momentum) rhn , & ! sqrt of neutral exchange coefficient (heat) ren , & ! sqrt of neutral exchange coefficient (water) rd , & ! sqrt of exchange coefficient (momentum) re , & ! sqrt of exchange coefficient (water) rh , & ! sqrt of exchange coefficient (heat) vmag , & ! surface wind magnitude (m/s) alz , & ! ln(zlvla/z10) thva , & ! virtual temperature (K) cp , & ! specific heat of moist air hol , & ! H (at zlvla) over L stable, & ! stability factor psixh ! stability function at zlvla (heat and water) real (kind=dbl_kind), parameter :: & cpvir = cp_wv/cp_air-c1, & ! defined as cp_wv/cp_air - 1. zTrf = c2 ! reference height for air temp (m) ! local functions real (kind=dbl_kind) :: & xd , & ! dummy argument psimhu, & ! unstable part of psimh psixhu ! unstable part of psimx !------------------------------------------------------------ ! Define functions !------------------------------------------------------------ psimhu(xd) = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) & - c2*atan(xd) + pih !ech - c2*atan(xd) + 1.571_dbl_kind psixhu(xd) = c2 * log((c1 + xd*xd)/c2) al2 = log(zrefa/zTrf) !------------------------------------------------------------ ! Initialize !------------------------------------------------------------ if (highfreq) then ! high frequency coupling follows Roberts et al. (2015) if (my_task==master_task.and.firstpass.and.sfctype(1:3)=='ice') then if (present(uice) .and. present(vice)) then write(nu_diag,*)'Using high frequency RASM atmospheric coupling' else call abort_ice('High frequency RASM coupling missing uice and vice') endif firstpass = .false. endif umin = p5 ! minumum allowable wind-ice speed difference of 0.5 m/s else umin = c1 ! minumum allowable wind speed of 1m/s endif do j = 1, ny_block do i = 1, nx_block if (present(Uref)) then Uref(i,j) = c0 endif Tref(i,j) = c0 Qref(i,j) = c0 delt(i,j) = c0 delq(i,j) = c0 shcoef(i,j) = c0 lhcoef(i,j) = c0 enddo enddo !------------------------------------------------------------ ! Compute turbulent flux coefficients, wind stress, and ! reference temperature and humidity. !------------------------------------------------------------ !------------------------------------------------------------ ! define variables that depend on surface type !------------------------------------------------------------ if (sfctype(1:3)=='ice') then qqq = qqqice ! for qsat TTT = TTTice ! for qsat Lheat = Lsub ! ice to vapor do ij = 1, icells i = indxi(ij) j = indxj(ij) if (highfreq) then vmag(ij) = max(umin, sqrt( (uatm(i,j)-uice(i,j))**2 + & (vatm(i,j)-vice(i,j))**2) ) else vmag(ij) = max(umin, wind(i,j)) endif if (formdrag .and. Cdn_atm(i,j) > puny) then rdn(ij) = sqrt(Cdn_atm(i,j)) else rdn(ij) = vonkar/log(zrefa/icerufa) ! neutral coefficient Cdn_atm(i,j) = rdn(ij) * rdn(ij) endif enddo ! ij elseif (sfctype(1:3)=='ocn') then qqq = qqqocn TTT = TTTocn Lheat = Lvap ! liquid to vapor do ij = 1, icells i = indxi(ij) j = indxj(ij) vmag(ij) = max(umin, wind(i,j)) rdn(ij) = sqrt(0.0027_dbl_kind/vmag(ij) & + .000142_dbl_kind + .0000764_dbl_kind*vmag(ij)) enddo ! ij endif ! sfctype do ij = 1, icells i = indxi(ij) j = indxj(ij) !------------------------------------------------------------ ! define some more needed variables !------------------------------------------------------------ TsfK = Tsf(i,j) + Tffresh ! surface temp (K) qsat = qqq * exp(-TTT/TsfK) ! saturation humidity (kg/m^3) ssq = qsat / rhoa(i,j) ! sat surf hum (kg/kg) thva(ij) = potT(i,j) * (c1 + zvir * Qa(i,j)) ! virtual pot temp (K) delt(i,j) = potT(i,j) - TsfK ! pot temp diff (K) delq(i,j) = Qa(i,j) - ssq ! spec hum dif (kg/kg) alz(ij) = log(zlvla(i,j)/zrefa) cp(ij) = cp_air*(c1 + cpvir*ssq) !------------------------------------------------------------ ! first estimate of Z/L and ustar, tstar and qstar !------------------------------------------------------------ ! neutral coefficients, z/L = 0.0 rhn(ij) = rdn(ij) ren(ij) = rdn(ij) ! ustar,tstar,qstar ustar(ij) = rdn(ij) * vmag(ij) tstar(ij) = rhn(ij) * delt(i,j) qstar(ij) = ren(ij) * delq(i,j) enddo ! ij !------------------------------------------------------------ ! iterate to converge on Z/L, ustar, tstar and qstar !------------------------------------------------------------ do k = 1, natmiter do ij = 1, icells i = indxi(ij) j = indxj(ij) ! compute stability & evaluate all stability functions hol(ij) = vonkar * gravit * zlvla(i,j) & * (tstar(ij)/thva(ij) & + qstar(ij)/(c1/zvir+Qa(i,j))) & / ustar(ij)**2 hol(ij) = sign( min(abs(hol(ij)),c10), hol(ij) ) stable(ij) = p5 + sign(p5 , hol(ij)) xqq = max(sqrt(abs(c1 - c16*hol(ij))) , c1) xqq = sqrt(xqq) ! Jordan et al 1999 psimhs = -(0.7_dbl_kind*hol(ij) & + 0.75_dbl_kind*(hol(ij)-14.3_dbl_kind) & * exp(-0.35_dbl_kind*hol(ij)) + 10.7_dbl_kind) psimh = psimhs*stable(ij) & + (c1 - stable(ij))*psimhu(xqq) psixh(ij) = psimhs*stable(ij) & + (c1 - stable(ij))*psixhu(xqq) ! shift all coeffs to measurement height and stability rd(ij) = rdn(ij) / (c1+rdn(ij)/vonkar*(alz(ij)-psimh)) rh(ij) = rhn(ij) / (c1+rhn(ij)/vonkar*(alz(ij)-psixh(ij))) re(ij) = ren(ij) / (c1+ren(ij)/vonkar*(alz(ij)-psixh(ij))) ! update ustar, tstar, qstar using updated, shifted coeffs ustar(ij) = rd(ij) * vmag(ij) tstar(ij) = rh(ij) * delt(i,j) qstar(ij) = re(ij) * delq(i,j) enddo ! ij enddo ! end iteration if (calc_strair) then ! initialize do j = 1, ny_block do i = 1, nx_block strx(i,j) = c0 stry(i,j) = c0 enddo enddo !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) if (highfreq .and. sfctype(1:3)=='ice') then !------------------------------------------------------------ ! momentum flux for RASM !------------------------------------------------------------ ! tau = rhoa(i,j) * rd * rd ! strx = tau * |Uatm-U| * (uatm-u) ! stry = tau * |Uatm-U| * (vatm-v) !------------------------------------------------------------ tau = rhoa(i,j) * rd(ij) * rd(ij) ! not the stress at zlvla(i,j) ! high frequency momentum coupling following Roberts et al. (2015) strx(i,j) = tau * sqrt((uatm(i,j)-uice(i,j))**2 + & (vatm(i,j)-vice(i,j))**2) * & (uatm(i,j)-uice(i,j)) stry(i,j) = tau * sqrt((uatm(i,j)-uice(i,j))**2 + & (vatm(i,j)-vice(i,j))**2) * & (vatm(i,j)-vice(i,j)) else !------------------------------------------------------------ ! momentum flux !------------------------------------------------------------ ! tau = rhoa(i,j) * ustar * ustar ! strx = tau * uatm(i,j) / vmag ! stry = tau * vatm(i,j) / vmag !------------------------------------------------------------ tau = rhoa(i,j) * ustar(ij) * rd(ij) ! not the stress at zlvla(i,j) strx(i,j) = tau * uatm(i,j) stry(i,j) = tau * vatm(i,j) endif !------------------------------------------------------------ ! Stability and height corrected atm drag coefficient !------------------------------------------------------------ Cd_atm_n(i,j) = rd(ij)*rd(ij) enddo ! ij endif ! calc_strair !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) !------------------------------------------------------------ ! coefficients for turbulent flux calculation !------------------------------------------------------------ ! add windless coefficient for sensible heat flux ! as in Jordan et al (JGR, 1999) !------------------------------------------------------------ shcoef(i,j) = rhoa(i,j) * ustar(ij) * cp(ij) * rh(ij) + c1 lhcoef(i,j) = rhoa(i,j) * ustar(ij) * Lheat * re(ij) !------------------------------------------------------------ ! Compute diagnostics: 2m ref T, Q, U !------------------------------------------------------------ hol(ij) = hol(ij)*zTrf/zlvla(i,j) xqq = max( c1, sqrt(abs(c1-c16*hol(ij))) ) xqq = sqrt(xqq) psix2 = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psixhu(xqq) fac = (rh(ij)/vonkar) & * (alz(ij) + al2 - psixh(ij) + psix2) Tref(i,j)= potT(i,j) - delt(i,j)*fac Tref(i,j)= Tref(i,j) - p01*zTrf ! pot temp to temp correction fac = (re(ij)/vonkar) & * (alz(ij) + al2 - psixh(ij) + psix2) Qref(i,j)= Qa(i,j) - delq(i,j)*fac if (present(Uref)) then if (highfreq .and. sfctype(1:3)=='ice') then Uref(i,j) = sqrt((uatm(i,j)-uice(i,j))**2 + & (vatm(i,j)-vice(i,j))**2) * & rd(ij) / rdn(ij) else Uref(i,j) = vmag(ij) * rd(ij) / rdn(ij) endif endif ! (present(Uref)) enddo ! ij end subroutine atmo_boundary_layer !======================================================================= ! Compute coefficients for atm/ice fluxes, stress ! NOTE: \\ ! (1) all fluxes are positive downward, \\ ! (2) reference temperature and humidity are NOT computed subroutine atmo_boundary_const (nx_block, ny_block, & sfctype, icells, & indxi, indxj, & uatm, vatm, & wind, rhoa, & strx, stry, & Tsf, potT, & Qa, & delt, delq, & lhcoef, shcoef, & Cdn_atm) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icells ! number of cells that require atmo fluxes integer (kind=int_kind), dimension(nx_block*ny_block), & intent(in) :: & indxi, indxj ! compressed i and j indices character (len=3), intent(in) :: & sfctype ! ice or ocean real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & Tsf , & ! surface temperature of ice or ocean potT , & ! air potential temperature (K) Qa , & ! specific humidity (kg/kg) uatm , & ! x-direction wind speed (m/s) vatm , & ! y-direction wind speed (m/s) wind , & ! wind speed (m/s) rhoa ! air density (kg/m^3) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & Cdn_atm ! neutral drag coefficient real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout):: & strx , & ! x surface stress (N) stry ! y surface stress (N) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: & delt , & ! potential T difference (K) delq , & ! humidity difference (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat ! local variables integer (kind=int_kind) :: & i, j, & ! horizontal indices ij ! combined ij index real (kind=dbl_kind) :: & TsfK, & ! surface temperature in Kelvin (K) qsat, & ! the saturation humidity of air (kg/m^3) ssq , & ! sat surface humidity (kg/kg) tau, & ! stress at zlvla Lheat ! Lvap or Lsub, depending on surface type !------------------------------------------------------------ ! Initialize !------------------------------------------------------------ do j = 1, ny_block do i = 1, nx_block delt(i,j) = c0 delq(i,j) = c0 shcoef(i,j) = c0 lhcoef(i,j) = c0 enddo enddo if (calc_strair) then do j = 1, ny_block do i = 1, nx_block strx(i,j) = c0 stry(i,j) = c0 enddo enddo !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) !------------------------------------------------------------ ! momentum flux !------------------------------------------------------------ tau = rhoa(i,j) * 0.0012_dbl_kind * wind(i,j) !AOMIP tau = rhoa(i,j) * (1.10_dbl_kind + c4*p01*wind(i,j)) & !AOMIP * wind(i,j) * p001 strx(i,j) = tau * uatm(i,j) stry(i,j) = tau * vatm(i,j) enddo ! ij endif ! calc_strair !------------------------------------------------------------ ! define variables that depend on surface type !------------------------------------------------------------ if (sfctype(1:3)=='ice') then Lheat = Lsub ! ice to vapor elseif (sfctype(1:3)=='ocn') then Lheat = Lvap ! liquid to vapor endif ! sfctype !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) !------------------------------------------------------------ ! potential temperature and specific humidity differences !------------------------------------------------------------ TsfK = Tsf(i,j) + Tffresh ! surface temp (K) qsat = qqqocn * exp(-TTTocn/TsfK) ! sat humidity (kg/m^3) ssq = qsat / rhoa(i,j) ! sat surf hum (kg/kg) delt(i,j)= potT(i,j) - TsfK ! pot temp diff (K) delq(i,j)= Qa(i,j) - ssq ! spec hum dif (kg/kg) !------------------------------------------------------------ ! coefficients for turbulent flux calculation !------------------------------------------------------------ shcoef(i,j) = (1.20e-3_dbl_kind)*cp_air*rhoa(i,j)*wind(i,j) lhcoef(i,j) = (1.50e-3_dbl_kind)*Lheat *rhoa(i,j)*wind(i,j) enddo ! ij end subroutine atmo_boundary_const !======================================================================= ! Neutral drag coefficients for the atmosphere and ocean taking into account ! skin drag, and form drag from ridges and keels, floe edges and melt ponds. ! Equations and parameters cited in this subroutine are from: ! 1) Tsamados et al. (2014), JPO, DOI: 10.1175/JPO-D-13-0215.1 ! 2) Lupkes et al. (2012), JGR, DOI:10.1029/2012JD017630 ! 3) Shaw et al. (2008), JGR, DOI: 10.1029/2007JC004550 ! 4) Stoessel & Claussen (1993), Clim. Dyn. DOI:10.1007/BF00210010 ! 5) Steele 1989 ! 6) Hanssen1988 ! 7) Lu 2011 ! ! authors: Michel Tsamados, CPOM ! David Schroeder, CPOM ! Andrew Roberts, NPS ! Elizabeth Hunke, LANL subroutine neutral_drag_coeffs (nx_block, ny_block, & ilo, ihi, jlo, jhi, & apnd, hpnd, & ipnd, & alvl, vlvl, & aice, vice, & vsno, z0ocn, & aicen, & vicen, vsnon, & Cdn_ocn, Cdn_ocn_skin, & Cdn_ocn_floe, Cdn_ocn_keel,& Cdn_atm, Cdn_atm_skin, & Cdn_atm_floe, Cdn_atm_pond,& Cdn_atm_rdg, hfreebd, & hdraft, hridge, & distrdg, hkeel, & dkeel, lfloe, & dfloe, ncat) !----------------------------------------------------------------- ! START TEMPORARY MODULE USAGE use ice_communicate, only: my_task, master_task use ice_fileunits, only: nu_diag ! END TEMPORARY MODULE USAGE !----------------------------------------------------------------- integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ncat integer (kind=int_kind), dimension(nx_block*ny_block) :: & indxi, indxj ! compressed i and j indices real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & intent(in) :: & apnd ,& ! melt pond fraction of sea ice hpnd ,& ! mean melt pond depth over sea ice ipnd ,& ! mean ice pond depth over sea ice in cat n alvl ,& ! level ice area fraction (of grid cell ?) vlvl ! level ice mean thickness real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & aice , & ! concentration of ice vice , & ! volume per unit area of ice vsno , & ! volume per unit area of snow z0ocn ! ocean roughness length 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), & intent(out) :: & hfreebd , & ! freeboard (m) hdraft , & ! draft of ice + snow column (Stössel & Clausse 1993) hridge , & ! ridge height distrdg , & ! distance between ridges hkeel , & ! keel depth dkeel , & ! distance between keels lfloe , & ! floe length (m) dfloe , & ! distance between floes Cdn_ocn , & ! ocean-ice neutral drag coefficient Cdn_ocn_skin , & ! drag coefficient due to skin drag Cdn_ocn_floe , & ! drag coefficient due to floe edges Cdn_ocn_keel , & ! drag coefficient dur to keels Cdn_atm , & ! ice-atmosphere drag coefficient Cdn_atm_skin , & ! drag coefficient due to skin drag Cdn_atm_floe , & ! drag coefficient due to floe edges Cdn_atm_pond , & ! drag coefficient dur to ponds Cdn_atm_rdg ! drag coefficient dur to ridges !----------------------------------------------------------------- ! START TEMPORARY COMPONENT SWITCHES logical, parameter :: & ridgkeeldrag = .true., & ! Form drag due to ridges and keels floeedgedrag = .true., & ! Form drag due to floe edges meltponddrag = .true. ! Form drag on upper surface due to melt ponds ! END TEMPORARY COMPONENT SWITCHES !----------------------------------------------------------------- real (kind=dbl_kind), parameter :: & ! [,] = range of values that can be tested mrdga = c20 ,&! atmos screening effect see Lu2011 [5,50] mrdgo = c10 ,&! ocn screening effect see Lu2011 [5,50] beta = p5 ,&! power exponent appearing in astar and ! L=Lmin(A*/(A*-A))**beta [0,1] Lmin = c8 ,&! min length of floe (m) [5,100] Lmax = 300._dbl_kind ,&! max length of floe (m) [30,3000] cfa = p2 ,&! Eq. 12 ratio of local form drag over ! geometrical parameter for floes [0,1] cfw = p2 ,&! Eq. 15 ratio of local form drag over ! geometrical parameter for floes [0,1] cpa = p2 ,&! Eq. 16 ratio of local form drag over ! geometrical parameter [0,1] cra = p2 ,&! Eq. 10 local form drag coefficient [0,1] crw = p2 ,&! Eq. 11 local form drag coefficient [0,1] sl = 22._dbl_kind ,&! Sheltering parameter Lupkes et al 2012 [10,30] lpmin = 2.26_dbl_kind ,&! min pond length (m) Tsamados et al 2014 [1,10] lpmax = 24.63_dbl_kind ,&! max pond length (m) Tsamados et al 2014 [10,100] tanar = 0.404_dbl_kind ,&! tan ridge slope [0.4,1] tanak = 0.404_dbl_kind ,&! tan keel slope [0.4,1] phir = 0.8_dbl_kind ,&! porosity of ridges [0.4,1] phik = 0.8_dbl_kind ,&! porosity of keels [0.4,1] sHGB = 0.18_dbl_kind ! Lupkes et al 2012 Eq. 28, Hanssen1988, ! Steele1989 suggest instead 0.18 !----------------------------------------------------------------- ! START TEMPORARY SETTING logical, save :: & first=.true. ! END TEMPORARY SETTING !----------------------------------------------------------------- integer (kind=int_kind) :: & icells, & ! number of cells that require atmo fluxes n , & ! category index i, j , & ! horizontal indices ij ! combined ij index real (kind=dbl_kind) :: & astar, & ! new constant for form drag ctecaf, & ! constant ctecwf, & ! constant sca, & ! wind attenuation function scw, & ! ocean attenuation function lp, & ! pond length (m) ctecar, & ! WHAT ARE THESE? ctecwk, & ! WHAT ARE THESE? ai, aii ! ice area and its inverse real (kind=dbl_kind) :: & apond , & ! melt pond fraction of grid cell vpond , & ! mean melt pond depth over grid cell ipond , & ! mean melt pond ice depth over grid cell ardg , & ! ridged ice area fraction of grid cell vrdg , & ! ridged ice mean thickness adi , & ! mean deformed ice fraction of grid cell vds , & ! mean deformed snow thickness over grid cell vdi , & ! mean deformed ice thickness over grid cell hdf , & ! mean deformed ice freeboard hdd ! mean deformed ice draft !----------------------------------------------------------------- ! START TEMPORARY OUTPUT if (my_task==master_task.and.first) then write(nu_diag,*)'-------------------------------------------' write(nu_diag,*)'ridgkeeldrag=',ridgkeeldrag write(nu_diag,*)'floeedgedrag=',floeedgedrag write(nu_diag,*)'meltponddrag=',meltponddrag write(nu_diag,*)' phir=',phir write(nu_diag,*)' phik=',phik write(nu_diag,*)' tanar=',tanar write(nu_diag,*)' tanak=',tanak write(nu_diag,*)' icerufas=',icerufas write(nu_diag,*)' icerufos=',icerufos write(nu_diag,*)' sl=',sl write(nu_diag,*)' beta=',beta write(nu_diag,*)' Lmin=',Lmin write(nu_diag,*)' Lmax=',Lmax write(nu_diag,*)' lpmin=',lpmin write(nu_diag,*)' lpmax=',lpmax write(nu_diag,*)' mrdga=',mrdga write(nu_diag,*)' mrdgo=',mrdgo write(nu_diag,*)'-------------------------------------------' first=.false. endif ! END TEMPORARY OUTPUT !----------------------------------------------------------------- !----------------------------------------------------------------- ! Initialize parametric constants !----------------------------------------------------------------- astar = c1/(c1-(Lmin/Lmax)**(c1/beta)) !----------------------------------------------------------------- ! Initialize fundamental arrays, setting to background values !----------------------------------------------------------------- hfreebd(:,:) = c0 hdraft (:,:) = c0 hridge (:,:) = c0 distrdg(:,:) = c0 hkeel (:,:) = c0 dkeel (:,:) = c0 lfloe (:,:) = c0 dfloe (:,:) = c0 Cdn_ocn_skin(:,:) = c0 Cdn_ocn_floe(:,:) = c0 Cdn_ocn_keel(:,:) = c0 Cdn_atm_skin(:,:) = c0 Cdn_atm_floe(:,:) = c0 Cdn_atm_pond(:,:) = c0 Cdn_atm_rdg (:,:) = c0 Cdn_ocn(:,:) = (vonkar/log(zrefo/icerufo)) * (vonkar/log(zrefo/icerufo)) Cdn_atm(:,:) = (vonkar/log(zrefa/icerufa)) * (vonkar/log(zrefa/icerufa)) !----------------------------------------------------------------- ! Identify cells with ice area exceeding the ice dynamics minimum !----------------------------------------------------------------- icells = 0 do j = jlo, jhi do i = ilo, ihi if (aice(i,j) > p001) then icells = icells + 1 indxi(icells) = i indxj(icells) = j endif enddo ! i enddo ! j !----------------------------------------------------------------- ! Primary skin and form drag loop !----------------------------------------------------------------- !DIR$ CONCURRENT !Cray !cdir nodep !NEC !ocl novrec !Fujitsu do ij = 1, icells i = indxi(ij) j = indxj(ij) !------------------------------------------------------------ ! Set temporary concentration variable and its inverse !------------------------------------------------------------ ! concentration and its inverse ai = aice(i,j) aii = c1/ai !------------------------------------------------------------ ! Determine mean melt pond area, volume and depth !------------------------------------------------------------ apond = c0 vpond = c0 ipond = c0 do n = 1,ncat ! area of pond per unit area of grid cell apond = apond+apnd(i,j,n)*aicen(i,j,n) ! volume of pond per unit area of grid cell vpond = vpond+apnd(i,j,n)*hpnd(i,j,n)*aicen(i,j,n) ! volume of lid per unit area of grid cell ipond = ipond+apnd(i,j,n)*ipnd(i,j,n)*aicen(i,j,n) enddo !------------------------------------------------------------ ! Determine mean draft and freeboard !------------------------------------------------------------ ! mean draft and freeboard (Eq. 27 Tsamados et al. 2014) hdraft(i,j) = (rhoi*vice(i,j)+rhos*vsno(i,j))*aii/rhow ! without ponds hfreebd(i,j) = (vice(i,j)+vsno(i,j))*aii-hdraft(i,j) ! adjust freeboard for snow ice (Eq. 5 Stoessel & Claussen1993) if (hdraft(i,j) >= vice(i,j)*aii) then hfreebd(i,j) = (hdraft(i,j)*ai*(c1-rhoi/rhow) + & (vsno(i,j)-(vice(i,j)-hdraft(i,j)*ai)*rhoi/rhos) * & (c1-rhos/rhow))*aii endif !------------------------------------------------------------ ! Compute mean floe size metrics !------------------------------------------------------------ ! floe size parameterization (Eq. 13 Tsamados et al. 2014) lfloe(i,j) = Lmin * (astar / (astar - ai))**beta ! distance between floes parameterization (Eq. 14 Tsamados et al. 2014) dfloe(i,j) = lfloe(i,j) * (c1/sqrt(ai) - c1) !------------------------------------------------------------ ! Calculate deformed ice draft and freeboard !------------------------------------------------------------ adi = c0 vds = c0 vdi = c0 do n=1,ncat ! deformed ice area fraction over grid cell adi=adi+(c1-alvl(i,j,n))*aicen(i,j,n) ! total deformed snow volume per unit grid cell area vds=vds+(c1-alvl(i,j,n))*vsnon(i,j,n) ! total deformed ice volume per unit grid cell area vdi=vdi+(c1-vlvl(i,j,n))*vicen(i,j,n) enddo if (adi>puny) then hdd = (rhoi*vdi+rhos*vds)/(adi*rhow) hdf = ((vdi+vds)/adi)-hdd ! adjust freeboard for snow ice (Eq. 5 Stoessel & Claussen 1993) if (hdd >= vdi/adi) then hdf = (hdd*adi*(c1-rhoi/rhow) + & (vds-(vdi-hdf*adi)*rhoi/rhos) * & (c1-rhos/rhow))/adi endif else hdd = c0 hdf = c0 endif !------------------------------------------------------------ ! Isostatic ridge/keel structure !------------------------------------------------------------ if (hdd/phik>hdraft(i,j).and.hdf/phir>hfreebd(i,j).and.ridgkeeldrag) then hkeel(i,j) = (c2*hdd/phik) - hdraft(i,j) hridge(i,j) = hfreebd(i,j) + c2*sqrt(tanar*((hdd/phik)-hdraft(i,j)) & *((hdf/phir)-hfreebd(i,j))/tanak) dkeel(i,j) = c2*ai*(hkeel(i,j)-hdraft(i,j))/(tanak*adi) distrdg(i,j) = dkeel(i,j) !------------------------------------------------------------ ! Skin drag (atmosphere) !------------------------------------------------------------ ! Eq. X Lupkes et al. 2014 Cdn_atm_skin(i,j) = (c1 - mrdga*(hridge(i,j)-hfreebd(i,j))/distrdg(i,j)) * & (vonkar/log(zrefa/icerufas))**c2 !------------------------------------------------------------ ! Ridge effect (atmosphere) !------------------------------------------------------------ ! Eq. X Lupkes et al. 2014 sca = c1 - exp(-sHGB*distrdg(i,j)/(hridge(i,j)-hfreebd(i,j))) ctecar = cra*p5 Cdn_atm_rdg(i,j) = ctecar*(hridge(i,j)-hfreebd(i,j))/distrdg(i,j)*sca* & (log((hridge(i,j)-hfreebd(i,j))/icerufas)/log(zrefa/icerufas))**c2 !------------------------------------------------------------ ! Skin drag bottom ice (ocean) !------------------------------------------------------------ ! Eqs. X-Y Lupkes et al. 2014 Cdn_ocn_skin(i,j) = (c1 - mrdgo*(hkeel(i,j)-hdraft(i,j))/dkeel(i,j)) * & (vonkar/log(zrefo/icerufos))**c2 !------------------------------------------------------------ ! Keel effect (ocean) !------------------------------------------------------------ ! Eqs. X-Y Lupkes et al. 2014 scw = c1 - exp(-sHGB*dkeel(i,j)/(hkeel(i,j)-hdraft(i,j))) ctecwk = crw*p5 Cdn_ocn_keel(i,j) = ctecwk*(hkeel(i,j)-hdraft(i,j))/dkeel(i,j)*scw* & (log((hkeel(i,j)-hdraft(i,j))/icerufos)/log(zrefo/icerufos))**c2 else hkeel(i,j) = hdraft(i,j) hridge(i,j) = hfreebd(i,j) Cdn_atm_skin(i,j) = (vonkar/log(zrefa/icerufas))**2 Cdn_ocn_skin(i,j) = (vonkar/log(zrefo/icerufos))**2 endif ! (hdd/phir>hdraft(i,j) & hdf/phir>hfreebd(i,j)) !------------------------------------------------------------ ! Floe edge drag (atmosphere) !------------------------------------------------------------ if (hfreebd(i,j) > puny .and. floeedgedrag) then ! Eqs. X-Y Lupkes et al. 2014 sca = c1 - exp(-sl*beta*(c1-ai)) ctecaf = cfa*p5*((log(hfreebd(i,j)/z0ocn(i,j))/log(zrefa/z0ocn(i,j)))**c2)*sca Cdn_atm_floe(i,j) = ctecaf * hfreebd(i,j) / lfloe(i,j) endif !------------------------------------------------------------ ! Pond edge drag (atmosphere) !------------------------------------------------------------ if (hfreebd(i,j) > puny .and. meltponddrag) then ! Eqs. X-Y Lupkes et al. 2014 sca = (apond)**(c1/(zrefa*beta)) lp = lpmin*(1-apond)+lpmax*apond Cdn_atm_pond(i,j) = cpa*p5*sca*apond*hfreebd(i,j)/lp * & ((log(hfreebd(i,j)/z0ocn(i,j))/log(zrefa/z0ocn(i,j)))**c2) endif !------------------------------------------------------------ ! Floe edge drag (ocean) !------------------------------------------------------------ if (hdraft(i,j) > puny .and. floeedgedrag) then ! Eqs. X-Y Lupkes et al. 2014 scw = c1 - exp(-sl*beta*(c1-ai)) ctecwf = cfw*p5*((log(hdraft(i,j)/z0ocn(i,j))/log(zrefo/z0ocn(i,j)))**c2)*scw Cdn_ocn_floe(i,j) = ctecwf * hdraft(i,j) / lfloe(i,j) endif !------------------------------------------------------------ ! Total drag coefficient (atmosphere) !------------------------------------------------------------ ! Eqs. X Tsamados et al. 2014 Cdn_atm(i,j) = Cdn_atm_skin(i,j) + Cdn_atm_floe(i,j) + & Cdn_atm_pond(i,j) + Cdn_atm_rdg(i,j) !------------------------------------------------------------ ! Total drag coefficient (ocean) !------------------------------------------------------------ ! Eqs. X Tsamados et al. 2014 Cdn_ocn(i,j) = Cdn_ocn_skin(i,j) + Cdn_ocn_floe(i,j) + & Cdn_ocn_keel(i,j) enddo ! ij end subroutine neutral_drag_coeffs !======================================================================= end module ice_atmo !=======================================================================