!========================================================================= ! ! Update ice and snow internal temperatures and compute ! thermodynamic growth rates and atmospheric fluxes. ! ! NOTE: The thermodynamic calculation is split in two for load balancing. ! First icepack_therm_vertical computes vertical growth rates and coupler ! fluxes. Then icepack_therm_itd does thermodynamic calculations not ! needed for coupling. ! ! authors: William H. Lipscomb, LANL ! C. M. Bitz, UW ! Elizabeth C. Hunke, LANL ! ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb ! 2004: Block structure added by William Lipscomb ! 2006: Streamlined for efficiency by Elizabeth Hunke ! Converted to free source form (F90) module icepack_therm_vertical use icepack_kinds use icepack_parameters, only: c0, c1, p001, p5, puny use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity use icepack_parameters, only: ktherm, heat_capacity, calc_Tsfc use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd use icepack_parameters, only: phi_i_mushy, floeshape, floediam use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo use icepack_tracers, only: n_aero, n_iso use icepack_therm_shared, only: ferrmax, l_brine use icepack_therm_shared, only: calculate_tin_from_qin, Tmin use icepack_therm_shared, only: hi_min use icepack_therm_bl99, only: temperature_changes use icepack_therm_0layer, only: zerolayer_temperature use icepack_therm_mushy, only: temperature_changes_salinity use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted use icepack_mushy_physics, only: icepack_mushy_temperature_mush use icepack_mushy_physics, only: liquidus_temperature_mush use icepack_mushy_physics, only: enthalpy_mush, enthalpy_of_melting use icepack_aerosol, only: update_aerosol use icepack_isotope, only: update_isotope use icepack_atmo, only: neutral_drag_coeffs, icepack_atm_boundary use icepack_age, only: increment_age use icepack_firstyear, only: update_FYarea use icepack_flux, only: set_sfcflux, merge_fluxes use icepack_meltpond_cesm, only: compute_ponds_cesm use icepack_meltpond_lvl, only: compute_ponds_lvl use icepack_meltpond_topo, only: compute_ponds_topo implicit none private public :: frzmlt_bottom_lateral, & thermo_vertical, & icepack_step_therm1 !======================================================================= contains !======================================================================= ! ! Driver for updating ice and snow internal temperatures and ! computing thermodynamic growth rates and atmospheric fluxes. ! ! authors: William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine thermo_vertical (nilyr, nslyr, & dt, aicen, & vicen, vsnon, & Tsf, zSin, & zqin, zqsn, & apond, hpond, & tr_pond_topo,& flw, potT, & Qa, rhoa, & fsnow, fpond, & fbot, Tbot, & Tsnice, sss, & lhcoef, shcoef, & fswsfc, fswint, & Sswabs, Iswabs, & fsurfn, fcondtopn, & fcondbotn, & fsensn, flatn, & flwoutn, evapn, & evapsn, evapin, & freshn, fsaltn, & fhocnn, meltt, & melts, meltb, & congel, snoice, & mlt_onset, frz_onset, & yday, dsnow, & prescribed_ice) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt ! time step ! ice state variables real (kind=dbl_kind), intent(inout) :: & aicen , & ! concentration of ice vicen , & ! volume per unit area of ice (m) vsnon ! volume per unit area of snow (m) ! tracers real (kind=dbl_kind), intent(inout) :: & Tsf , & ! ice/snow top surface temp, same as Tsfcn (deg C) apond , & ! melt pond area fraction hpond ! melt pond depth (m) ! iage ! ice age (s) logical (kind=log_kind), intent(in) :: & tr_pond_topo ! if .true., use melt pond tracer logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed real (kind=dbl_kind), dimension (:), intent(inout) :: & zqsn , & ! snow layer enthalpy, zqsn < 0 (J m-3) zqin , & ! ice layer enthalpy, zqin < 0 (J m-3) zSin ! internal ice layer salinities ! input from atmosphere real (kind=dbl_kind), & intent(in) :: & flw , & ! incoming longwave radiation (W/m^2) potT , & ! air potential temperature (K) Qa , & ! specific humidity (kg/kg) rhoa , & ! air density (kg/m^3) fsnow , & ! snowfall rate (kg m-2 s-1) shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat real (kind=dbl_kind), & intent(inout) :: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) fswint , & ! SW absorbed in ice interior, below surface (W m-2) fpond ! fresh water flux to ponds (kg/m^2/s) real (kind=dbl_kind), dimension (:), intent(inout) :: & Sswabs , & ! SW radiation absorbed in snow layers (W m-2) Iswabs ! SW radiation absorbed in ice layers (W m-2) ! input from ocean real (kind=dbl_kind), intent(in) :: & fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) Tbot , & ! ice bottom surface temperature (deg C) sss ! ocean salinity ! coupler fluxes to atmosphere real (kind=dbl_kind), intent(out):: & flwoutn , & ! outgoing longwave radiation (W/m^2) evapn , & ! evaporative water flux (kg/m^2/s) evapsn , & ! evaporative water flux over snow (kg/m^2/s) evapin ! evaporative water flux over ice (kg/m^2/s) ! Note: these are intent out if calc_Tsfc = T, otherwise intent in real (kind=dbl_kind), intent(inout):: & fsensn , & ! sensible heat flux (W/m^2) flatn , & ! latent heat flux (W/m^2) fsurfn , & ! net flux to top surface, excluding fcondtopn fcondtopn, & ! downward cond flux at top surface (W m-2) fcondbotn ! downward cond flux at bottom surface (W m-2) ! coupler fluxes to ocean real (kind=dbl_kind), intent(out):: & freshn , & ! fresh water flux to ocean (kg/m^2/s) fsaltn , & ! salt flux to ocean (kg/m^2/s) fhocnn ! net heat flux to ocean (W/m^2) ! diagnostic fields real (kind=dbl_kind), & intent(inout):: & Tsnice , & ! snow ice interface temperature (deg C) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) meltb , & ! basal ice melt (m/step-->cm/day) congel , & ! basal ice growth (m/step-->cm/day) snoice , & ! snow-ice formation (m/step-->cm/day) dsnow , & ! change in snow thickness (m/step-->cm/day) mlt_onset, & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) real (kind=dbl_kind), intent(in) :: & yday ! day of year ! local variables integer (kind=int_kind) :: & k ! ice layer index real (kind=dbl_kind) :: & dhi , & ! change in ice thickness dhs ! change in snow thickness ! 2D state variables (thickness, temperature) real (kind=dbl_kind) :: & hilyr , & ! ice layer thickness hslyr , & ! snow layer thickness hin , & ! ice thickness (m) hsn , & ! snow thickness (m) hsn_new , & ! thickness of new snow (m) worki , & ! local work array works ! local work array real (kind=dbl_kind), dimension (nilyr) :: & zTin ! internal ice layer temperatures real (kind=dbl_kind), dimension (nslyr) :: & zTsn ! internal snow layer temperatures ! other 2D flux and energy variables real (kind=dbl_kind) :: & einit , & ! initial energy of melting (J m-2) efinal , & ! final energy of melting (J m-2) einter ! intermediate energy real (kind=dbl_kind) :: & fadvocn ! advective heat flux to ocean character(len=*),parameter :: subname='(thermo_vertical)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- flwoutn = c0 evapn = c0 evapsn = c0 evapin = c0 freshn = c0 fsaltn = c0 fhocnn = c0 fadvocn = c0 meltt = c0 meltb = c0 melts = c0 congel = c0 snoice = c0 dsnow = c0 zTsn(:) = c0 zTin(:) = c0 if (calc_Tsfc) then fsensn = c0 flatn = c0 fsurfn = c0 fcondtopn = c0 endif !----------------------------------------------------------------- ! Compute variables needed for vertical thermo calculation !----------------------------------------------------------------- call init_vertical_profile (nilyr, nslyr, & aicen, & vicen, vsnon, & hin, hilyr, & hsn, hslyr, & zqin, zTin, & zqsn, zTsn, & zSin, & einit ) if (icepack_warnings_aborted(subname)) return ! Save initial ice and snow thickness (for fresh and fsalt) worki = hin works = hsn !----------------------------------------------------------------- ! Compute new surface temperature and internal ice and snow ! temperatures. !----------------------------------------------------------------- if (heat_capacity) then ! usual case if (ktherm == 2) then call temperature_changes_salinity(dt, & nilyr, nslyr, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & fswsfc, fswint, & Sswabs, Iswabs, & hilyr, hslyr, & apond, hpond, & zqin, zTin, & zqsn, zTsn, & zSin, & Tsf, Tbot, & sss, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn, fcondbotn, & fadvocn, snoice) if (icepack_warnings_aborted(subname)) return else ! ktherm call temperature_changes(dt, & nilyr, nslyr, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & fswsfc, fswint, & Sswabs, Iswabs, & hilyr, hslyr, & zqin, zTin, & zqsn, zTsn, & zSin, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn, fcondbotn, & einit ) if (icepack_warnings_aborted(subname)) return endif ! ktherm else if (calc_Tsfc) then call zerolayer_temperature(nilyr, nslyr, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & fswsfc, & hilyr, hslyr, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn, fcondbotn ) if (icepack_warnings_aborted(subname)) return else !------------------------------------------------------------ ! Set fcondbot = fcondtop for zero layer thermodynamics ! fcondtop is set in call to set_sfcflux in step_therm1 !------------------------------------------------------------ fcondbotn = fcondtopn ! zero layer endif ! calc_Tsfc endif ! heat_capacity ! intermediate energy for error check einter = c0 do k = 1, nslyr einter = einter + hslyr * zqsn(k) enddo ! k do k = 1, nilyr einter = einter + hilyr * zqin(k) enddo ! k Tsnice = c0 if ((hslyr+hilyr) > puny) then if (hslyr > puny) then Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr) else Tsnice = Tsf endif endif if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Compute growth and/or melting at the top and bottom surfaces. ! Add new snowfall. ! Repartition ice into equal-thickness layers, conserving energy. !----------------------------------------------------------------- call thickness_changes(nilyr, nslyr, & dt, yday, & efinal, & hin, hilyr, & hsn, hslyr, & zqin, zqsn, & fbot, Tbot, & flatn, fsurfn, & fcondtopn, fcondbotn, & fsnow, hsn_new, & fhocnn, evapn, & evapsn, evapin, & meltt, melts, & meltb, & congel, snoice, & mlt_onset, frz_onset, & zSin, sss, & dsnow) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Check for energy conservation by comparing the change in energy ! to the net energy input !----------------------------------------------------------------- call conservation_check_vthermo(dt, & fsurfn, flatn, & fhocnn, fswint, & fsnow, einit, & einter, efinal, & fcondtopn, fcondbotn, & fadvocn, fbot ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! If prescribed ice, set hi back to old values !----------------------------------------------------------------- #ifdef CESMCOUPLED if (present(prescribed_ice)) then if (prescribed_ice) then hin = worki fhocnn = c0 ! for diagnostics endif endif #endif !----------------------------------------------------------------- ! Compute fluxes of water and salt from ice to ocean. ! evapn < 0 => sublimation, evapn > 0 => condensation ! aerosol flux is accounted for in icepack_aerosol.F90 !----------------------------------------------------------------- dhi = hin - worki dhs = hsn - works - hsn_new freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt fhocnn = fhocnn + fadvocn ! for ktherm=2 if (hin == c0) then if (tr_pond_topo) fpond = fpond - aicen * apond * hpond endif !----------------------------------------------------------------- ! Given the vertical thermo state variables, compute the new ice ! state variables. !----------------------------------------------------------------- call update_state_vthermo(nilyr, nslyr, & Tbot, Tsf, & hin, hsn, & zqin, zSin, & zqsn, & aicen, & vicen, vsnon) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Reload passive tracer array !----------------------------------------------------------------- end subroutine thermo_vertical !======================================================================= ! ! Compute heat flux to bottom surface. ! Compute fraction of ice that melts laterally. ! ! authors C. M. Bitz, UW ! William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL subroutine frzmlt_bottom_lateral (dt, ncat, & nilyr, nslyr, & aice, frzmlt, & vicen, vsnon, & qicen, qsnon, & sst, Tf, & ustar_min, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & fside) 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), intent(in) :: & dt ! time step real (kind=dbl_kind), intent(in) :: & aice , & ! ice concentration frzmlt , & ! freezing/melting potential (W/m^2) sst , & ! sea surface temperature (C) Tf , & ! freezing temperature (C) ustar_min,& ! minimum friction velocity for ice-ocean heat flux Cdn_ocn , & ! ocean-ice neutral drag coefficient strocnxT, & ! ice-ocean stress, x-direction strocnyT ! ice-ocean stress, y-direction character (char_len), intent(in) :: & fbot_xfer_type ! transfer coefficient type for ice-ocean heat flux real (kind=dbl_kind), dimension(:), intent(in) :: & vicen , & ! ice volume (m) vsnon ! snow volume (m) real (kind=dbl_kind), dimension(:,:), intent(in) :: & qicen , & ! ice layer enthalpy (J m-3) qsnon ! snow layer enthalpy (J m-3) real (kind=dbl_kind), intent(out) :: & Tbot , & ! ice bottom surface temperature (deg C) fbot , & ! heat flux to ice bottom (W/m^2) rside , & ! fraction of ice that melts laterally fside ! lateral heat flux (W/m^2) ! local variables integer (kind=int_kind) :: & n , & ! thickness category index k ! layer index real (kind=dbl_kind) :: & etot , & ! total energy in column qavg ! average enthalpy in column (approximate) real (kind=dbl_kind) :: & deltaT , & ! SST - Tbot >= 0 ustar , & ! skin friction velocity for fbot (m/s) wlat , & ! lateral melt rate (m/s) xtmp ! temporary variable ! Parameters for bottom melting real (kind=dbl_kind) :: & cpchr ! -cp_ocn*rhow*exchange coefficient ! Parameters for lateral melting real (kind=dbl_kind), parameter :: & m1 = 1.6e-6_dbl_kind , & ! constant from Maykut & Perovich ! (m/s/deg^(-m2)) m2 = 1.36_dbl_kind ! constant from Maykut & Perovich ! (unitless) character(len=*),parameter :: subname='(frzmlt_bottom_lateral)' !----------------------------------------------------------------- ! Identify grid cells where ice can melt. !----------------------------------------------------------------- rside = c0 fside = c0 Tbot = Tf fbot = c0 wlat = c0 if (aice > puny .and. frzmlt < c0) then ! ice can melt !----------------------------------------------------------------- ! Use boundary layer theory for fbot. ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703. !----------------------------------------------------------------- deltaT = max((sst-Tbot),c0) ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2 ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow) ustar = max (ustar,ustar_min) if (trim(fbot_xfer_type) == 'Cdn_ocn') then ! Note: Cdn_ocn has already been used for calculating ustar ! (formdrag only) --- David Schroeder (CPOM) cpchr = -cp_ocn*rhow*Cdn_ocn else ! fbot_xfer_type == 'constant' ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut cpchr = -cp_ocn*rhow*0.006_dbl_kind endif fbot = cpchr * deltaT * ustar ! < 0 fbot = max (fbot, frzmlt) ! frzmlt < fbot < 0 !!! uncomment to use all frzmlt for standalone runs ! fbot = min (c0, frzmlt) !----------------------------------------------------------------- ! Compute rside. See these references: ! Maykut and Perovich (1987): JGR, 92, 7032-7044 ! Steele (1992): JGR, 97, 17,729-17,738 !----------------------------------------------------------------- wlat = m1 * deltaT**m2 ! Maykut & Perovich rside = wlat*dt*pi/(floeshape*floediam) ! Steele rside = max(c0,min(rside,c1)) !----------------------------------------------------------------- ! Compute heat flux associated with this value of rside. !----------------------------------------------------------------- do n = 1, ncat etot = c0 qavg = c0 ! melting energy/unit area in each column, etot < 0 do k = 1, nslyr etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind) qavg = qavg + qsnon(k,n) enddo do k = 1, nilyr etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind) qavg = qavg + qicen(k,n) enddo ! nilyr ! lateral heat flux, fside < 0 if (tr_fsd) then ! floe size distribution fside = fside + wlat*qavg else ! default floe size fside = fside + rside*etot/dt endif enddo ! n !----------------------------------------------------------------- ! Limit bottom and lateral heat fluxes if necessary. !----------------------------------------------------------------- xtmp = frzmlt/(fbot + fside + puny) xtmp = min(xtmp, c1) fbot = fbot * xtmp rside = rside * xtmp fside = fside * xtmp endif end subroutine frzmlt_bottom_lateral !======================================================================= ! ! Given the state variables (vicen, vsnon, zqin, etc.), ! compute variables needed for the vertical thermodynamics ! (hin, hsn, zTin, etc.) ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine init_vertical_profile(nilyr, nslyr, & aicen, vicen, & vsnon, & hin, hilyr, & hsn, hslyr, & zqin, zTin, & zqsn, zTsn, & zSin, & einit ) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), 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), intent(out):: & hilyr , & ! ice layer thickness hslyr , & ! snow layer thickness einit ! initial energy of melting (J m-2) real (kind=dbl_kind), intent(out):: & hin , & ! ice thickness (m) hsn ! snow thickness (m) real (kind=dbl_kind), dimension (:), intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) zTin ! internal ice layer temperatures real (kind=dbl_kind), dimension (:), intent(in) :: & zSin ! internal ice layer salinities real (kind=dbl_kind), dimension (:), & intent(inout) :: & zqsn , & ! snow enthalpy zTsn ! snow temperature ! local variables real (kind=dbl_kind), dimension(nilyr) :: & Tmlts ! melting temperature integer (kind=int_kind) :: & k ! ice layer index real (kind=dbl_kind) :: & rnslyr, & ! real(nslyr) Tmax ! maximum allowed snow/ice temperature (deg C) logical (kind=log_kind) :: & ! for vector-friendly error checks tsno_high , & ! flag for zTsn > Tmax tice_high , & ! flag for zTin > Tmlt tsno_low , & ! flag for zTsn < Tmin tice_low ! flag for zTin < Tmin character(len=*),parameter :: subname='(init_vertical_profile)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- rnslyr = real(nslyr,kind=dbl_kind) tsno_high = .false. tice_high = .false. tsno_low = .false. tice_low = .false. einit = c0 !----------------------------------------------------------------- ! Surface temperature, ice and snow thickness ! Initialize internal energy !----------------------------------------------------------------- hin = vicen / aicen hsn = vsnon / aicen hilyr = hin / real(nilyr,kind=dbl_kind) hslyr = hsn / rnslyr !----------------------------------------------------------------- ! Snow enthalpy and maximum allowed snow temperature ! If heat_capacity = F, zqsn and zTsn are used only for checking ! conservation. !----------------------------------------------------------------- do k = 1, nslyr !----------------------------------------------------------------- ! Tmax based on the idea that dT ~ dq / (rhos*cp_ice) ! dq ~ q dv / v ! dv ~ puny = eps11 ! where 'd' denotes an error due to roundoff. !----------------------------------------------------------------- if (hslyr > hs_min/rnslyr .and. heat_capacity) then ! zqsn < 0 Tmax = -zqsn(k)*puny*rnslyr / & (rhos*cp_ice*vsnon) else zqsn (k) = -rhos * Lfresh Tmax = puny endif !----------------------------------------------------------------- ! Compute snow temperatures from enthalpies. ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0. !----------------------------------------------------------------- zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice !----------------------------------------------------------------- ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin. !----------------------------------------------------------------- if (zTsn(k) > Tmax) then tsno_high = .true. elseif (zTsn(k) < Tmin) then tsno_low = .true. endif enddo ! nslyr !----------------------------------------------------------------- ! If zTsn is out of bounds, print diagnostics and exit. !----------------------------------------------------------------- if (tsno_high .and. heat_capacity) then do k = 1, nslyr if (hslyr > hs_min/rnslyr) then Tmax = -zqsn(k)*puny*rnslyr / & (rhos*cp_ice*vsnon) else Tmax = puny endif if (zTsn(k) > Tmax) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Starting thermo, zTsn > Tmax' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTsn=',zTsn(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Tmax=',Tmax call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zqsn',zqsn(k),-Lfresh*rhos,zqsn(k)+Lfresh*rhos call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn > Tmax" ) return endif enddo ! nslyr endif ! tsno_high if (tsno_low .and. heat_capacity) then do k = 1, nslyr if (zTsn(k) < Tmin) then ! allowing for roundoff error write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Starting thermo, zTsn < Tmin' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTsn=', zTsn(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Tmin=', Tmin call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zqsn', zqsn(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, hin call icepack_warnings_add(warnstr) write(warnstr,*) subname, hsn call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn < Tmin" ) return endif enddo ! nslyr endif ! tsno_low do k = 1, nslyr if (zTsn(k) > c0) then ! correct roundoff error zTsn(k) = c0 zqsn(k) = -rhos*Lfresh endif !----------------------------------------------------------------- ! initial energy per unit area of ice/snow, relative to 0 C !----------------------------------------------------------------- einit = einit + hslyr*zqsn(k) enddo ! nslyr do k = 1, nilyr !--------------------------------------------------------------------- ! Use initial salinity profile for thin ice !--------------------------------------------------------------------- if (ktherm == 1 .and. zSin(k) < min_salin-puny) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Starting zSin < min_salin, layer', k call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zSin =', zSin(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'min_salin =', min_salin call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting zSin < min_salin, layer" ) return endif if (ktherm == 2) then Tmlts(k) = liquidus_temperature_mush(zSin(k)) else Tmlts(k) = -zSin(k) * depressT endif !----------------------------------------------------------------- ! Compute ice enthalpy ! If heat_capacity = F, zqin and zTin are used only for checking ! conservation. !----------------------------------------------------------------- !----------------------------------------------------------------- ! Compute ice temperatures from enthalpies using quadratic formula !----------------------------------------------------------------- if (ktherm == 2) then zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k)) else zTin(k) = calculate_Tin_from_qin(zqin(k),Tmlts(k)) endif if (l_brine) then Tmax = Tmlts(k) else ! fresh ice Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen) endif !----------------------------------------------------------------- ! Check for zTin > Tmax and zTin < Tmin !----------------------------------------------------------------- if (zTin(k) > Tmax) then tice_high = .true. elseif (zTin(k) < Tmin) then tice_low = .true. endif !----------------------------------------------------------------- ! If zTin is out of bounds, print diagnostics and exit. !----------------------------------------------------------------- if (tice_high .and. heat_capacity) then if (l_brine) then Tmax = Tmlts(k) else ! fresh ice Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen) endif if (zTin(k) > Tmax) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Starting thermo, T > Tmax, layer', k call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'k:', k call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTin =',zTin(k),', Tmax=',Tmax call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zSin =',zSin(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'hin =',hin call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zqin =',zqin(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'qmlt=',enthalpy_of_melting(zSin(k)) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Tmlt=',Tmlts(k) call icepack_warnings_add(warnstr) if (ktherm == 2) then zqin(k) = enthalpy_of_melting(zSin(k)) - c1 zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k)) write(warnstr,*) subname, 'Corrected quantities' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zqin=',zqin(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTin=',zTin(k) call icepack_warnings_add(warnstr) else call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T > Tmax, layer" ) return endif endif endif ! tice_high if (tice_low .and. heat_capacity) then if (zTin(k) < Tmin) then write(warnstr,*) ' ' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Starting thermo T < Tmin, layer', k call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zTin =', zTin(k) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Tmin =', Tmin call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T < Tmin, layer" ) return endif endif ! tice_low !----------------------------------------------------------------- ! correct roundoff error !----------------------------------------------------------------- if (ktherm /= 2) then if (zTin(k) > c0) then zTin(k) = c0 zqin(k) = -rhoi*Lfresh endif endif ! echmod: is this necessary? ! if (ktherm == 1) then ! if (zTin(k)>= -zSin(k)*depressT) then ! zTin(k) = -zSin(k)*depressT - puny ! zqin(k) = -rhoi*cp_ocn*zSin(k)*depressT ! endif ! endif !----------------------------------------------------------------- ! initial energy per unit area of ice/snow, relative to 0 C !----------------------------------------------------------------- einit = einit + hilyr*zqin(k) enddo ! nilyr end subroutine init_vertical_profile !======================================================================= ! ! Compute growth and/or melting at the top and bottom surfaces. ! Convert snow to ice if necessary. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine thickness_changes (nilyr, nslyr, & dt, yday, & efinal, & hin, hilyr, & hsn, hslyr, & zqin, zqsn, & fbot, Tbot, & flatn, fsurfn, & fcondtopn, fcondbotn, & fsnow, hsn_new, & fhocnn, evapn, & evapsn, evapin, & meltt, melts, & meltb, & congel, snoice, & mlt_onset, frz_onset,& zSin, sss, & dsnow) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt , & ! time step yday ! day of the year real (kind=dbl_kind), intent(in) :: & fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) Tbot , & ! ice bottom surface temperature (deg C) fsnow , & ! snowfall rate (kg m-2 s-1) flatn , & ! surface downward latent heat (W m-2) fsurfn , & ! net flux to top surface, excluding fcondtopn fcondtopn ! downward cond flux at top surface (W m-2) real (kind=dbl_kind), intent(inout) :: & fcondbotn ! downward cond flux at bottom surface (W m-2) real (kind=dbl_kind), dimension (:), intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) zqsn ! snow layer enthalpy (J m-3) real (kind=dbl_kind), intent(inout) :: & hilyr , & ! ice layer thickness (m) hslyr ! snow layer thickness (m) real (kind=dbl_kind), intent(inout) :: & meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) meltb , & ! basal ice melt (m/step-->cm/day) congel , & ! basal ice growth (m/step-->cm/day) snoice , & ! snow-ice formation (m/step-->cm/day) dsnow , & ! snow formation (m/step-->cm/day) ! iage , & ! ice age (s) mlt_onset , & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) real (kind=dbl_kind), intent(inout) :: & hin , & ! total ice thickness (m) hsn ! total snow thickness (m) real (kind=dbl_kind), intent(out):: & efinal ! final energy of melting (J m-2) real (kind=dbl_kind), intent(out):: & fhocnn , & ! fbot, corrected for any surplus energy (W m-2) evapn , & ! ice/snow mass sublimated/condensed (kg m-2 s-1) evapsn , & ! ice/snow mass sublimated/condensed over snow (kg m-2 s-1) evapin ! ice/snow mass sublimated/condensed over ice (kg m-2 s-1) real (kind=dbl_kind), intent(out):: & hsn_new ! thickness of new snow (m) ! changes to zSin in this subroutine are not reloaded into the ! trcrn array for ktherm /= 2, so we could remove ktherm=2 conditionals real (kind=dbl_kind), dimension (:), intent(inout) :: & zSin ! ice layer salinity (ppt) real (kind=dbl_kind), intent(in) :: & sss ! ocean salinity (PSU) ! local variables integer (kind=int_kind) :: & k ! vertical index real (kind=dbl_kind) :: & esub , & ! energy for sublimation, > 0 (J m-2) econ , & ! energy for condensation, < 0 (J m-2) etop_mlt , & ! energy for top melting, > 0 (J m-2) ebot_mlt , & ! energy for bottom melting, > 0 (J m-2) ebot_gro , & ! energy for bottom growth, < 0 (J m-2) emlt_atm , & ! total energy of brine, from atmosphere (J m-2) emlt_ocn ! total energy of brine, to ocean (J m-2) real (kind=dbl_kind) :: & qbotmax , & ! max enthalpy of ice growing at bottom dhi , & ! change in ice thickness dhs , & ! change in snow thickness Ti , & ! ice temperature Ts , & ! snow temperature qbot , & ! enthalpy of ice growing at bottom surface (J m-3) qsub , & ! energy/unit volume to sublimate ice/snow (J m-3) hqtot , & ! sum of h*q for two layers wk1 , & ! temporary variable zqsnew , & ! enthalpy of new snow (J m-3) hstot , & ! snow thickness including new snow (m) Tmlts ! melting temperature real (kind=dbl_kind), dimension (nilyr+1) :: & zi1 , & ! depth of ice layer boundaries (m) zi2 ! adjusted depths, with equal hilyr (m) real (kind=dbl_kind), dimension (nslyr+1) :: & zs1 , & ! depth of snow layer boundaries (m) zs2 ! adjusted depths, with equal hslyr (m) real (kind=dbl_kind), dimension (nilyr) :: & dzi ! ice layer thickness after growth/melting real (kind=dbl_kind), dimension (nslyr) :: & dzs ! snow layer thickness after growth/melting real (kind=dbl_kind), dimension (nilyr) :: & qm , & ! energy of melting (J m-3) = zqin in BL99 formulation qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation real (kind=dbl_kind) :: & qbotm , & qbotp , & qbot0 character(len=*),parameter :: subname='(thickness_changes)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- dhi = c0 dhs = c0 hsn_new = c0 do k = 1, nilyr dzi(k) = hilyr enddo do k = 1, nslyr dzs(k) = hslyr enddo do k = 1, nilyr if (ktherm == 2) then qmlt(k) = enthalpy_of_melting(zSin(k)) else qmlt(k) = c0 endif qm(k) = zqin(k) - qmlt(k) emlt_atm = c0 emlt_ocn = c0 enddo !----------------------------------------------------------------- ! For l_brine = false (fresh ice), check for temperatures > 0. ! Melt ice or snow as needed to bring temperatures back to 0. ! For l_brine = true, this should not be necessary. !----------------------------------------------------------------- if (.not. l_brine) then do k = 1, nslyr Ts = (Lfresh + zqsn(k)/rhos) / cp_ice if (Ts > c0) then dhs = cp_ice*Ts*dzs(k) / Lfresh dzs(k) = dzs(k) - dhs zqsn(k) = -rhos*Lfresh endif enddo do k = 1, nilyr Ti = (Lfresh + zqin(k)/rhoi) / cp_ice if (Ti > c0) then dhi = cp_ice*Ti*dzi(k) / Lfresh dzi(k) = dzi(k) - dhi zqin(k) = -rhoi*Lfresh endif enddo ! k endif ! .not. l_brine !----------------------------------------------------------------- ! Compute energy available for sublimation/condensation, top melt, ! and bottom growth/melt. !----------------------------------------------------------------- wk1 = -flatn * dt esub = max(wk1, c0) ! energy for sublimation, > 0 econ = min(wk1, c0) ! energy for condensation, < 0 wk1 = (fsurfn - fcondtopn) * dt etop_mlt = max(wk1, c0) ! etop_mlt > 0 wk1 = (fcondbotn - fbot) * dt ebot_mlt = max(wk1, c0) ! ebot_mlt > 0 ebot_gro = min(wk1, c0) ! ebot_gro < 0 !-------------------------------------------------------------- ! Condensation (evapn > 0) ! Note: evapn here has unit of kg/m^2. Divide by dt later. ! This is the only case in which energy from the atmosphere ! is used for changes in the brine energy (emlt_atm). !-------------------------------------------------------------- evapn = c0 ! initialize evapsn = c0 ! initialize evapin = c0 ! initialize if (hsn > puny) then ! add snow with enthalpy zqsn(1) dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0 dzs(1) = dzs(1) + dhs evapn = evapn + dhs*rhos evapsn = evapsn + dhs*rhos else ! add ice with enthalpy zqin(1) dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0 dzi(1) = dzi(1) + dhi evapn = evapn + dhi*rhoi evapin = evapin + dhi*rhoi ! enthalpy of melt water emlt_atm = emlt_atm - qmlt(1) * dhi endif !-------------------------------------------------------------- ! Grow ice (bottom) !-------------------------------------------------------------- if (ktherm == 2) then qbotm = enthalpy_mush(Tbot, sss) qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy) qbot0 = qbotm - qbotp dhi = ebot_gro / qbotp ! dhi > 0 hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss emlt_ocn = emlt_ocn - qbot0 * dhi else Tmlts = -zSin(nilyr) * depressT ! enthalpy of new ice growing at bottom surface if (heat_capacity) then if (l_brine) then qbotmax = -p5*rhoi*Lfresh ! max enthalpy of ice growing at bottom qbot = -rhoi * (cp_ice * (Tmlts-Tbot) & + Lfresh * (c1-Tmlts/Tbot) & - cp_ocn * Tmlts) qbot = min (qbot, qbotmax) ! in case Tbot is close to Tmlt else qbot = -rhoi * (-cp_ice * Tbot + Lfresh) endif else ! zero layer qbot = -rhoi * Lfresh endif dhi = ebot_gro / qbot ! dhi > 0 hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbot hstot = c0 endif ! ktherm dzi(nilyr) = dzi(nilyr) + dhi if (dzi(nilyr) > puny) then zqin(nilyr) = hqtot / dzi(nilyr) if (ktherm == 2) then zSin(nilyr) = hstot / dzi(nilyr) qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr)) else qmlt(nilyr) = c0 endif qm(nilyr) = zqin(nilyr) - qmlt(nilyr) endif ! update ice age due to freezing (new ice age = dt) ! if (tr_iage) & ! iage = (iage*hin + dt*dhi) / (hin + dhi) ! history diagnostics congel = congel + dhi if (dhi > puny .and. frz_onset < puny) & frz_onset = yday do k = 1, nslyr !-------------------------------------------------------------- ! Remove internal snow melt !-------------------------------------------------------------- if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then dhs = max(-dzs(k), & -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k)) dzs(k) = dzs(k) + dhs zqsn(k) = -rhos * Lfresh melts = melts - dhs ! delta E = zqsn(k) + rhos * Lfresh endif !-------------------------------------------------------------- ! Sublimation of snow (evapn < 0) !-------------------------------------------------------------- qsub = zqsn(k) - rhos*Lvap ! qsub < 0 dhs = max (-dzs(k), esub/qsub) ! esub > 0, dhs < 0 dzs(k) = dzs(k) + dhs esub = esub - dhs*qsub esub = max(esub, c0) ! in case of roundoff error evapn = evapn + dhs*rhos evapsn = evapsn + dhs*rhos !-------------------------------------------------------------- ! Melt snow (top) !-------------------------------------------------------------- dhs = max(-dzs(k), etop_mlt/zqsn(k)) dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0 etop_mlt = etop_mlt - dhs*zqsn(k) etop_mlt = max(etop_mlt, c0) ! in case of roundoff error ! history diagnostics if (dhs < -puny .and. mlt_onset < puny) & mlt_onset = yday melts = melts - dhs enddo ! nslyr do k = 1, nilyr !-------------------------------------------------------------- ! Sublimation of ice (evapn < 0) !-------------------------------------------------------------- qsub = qm(k) - rhoi*Lvap ! qsub < 0 dhi = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0 dzi(k) = dzi(k) + dhi esub = esub - dhi*qsub esub = max(esub, c0) evapn = evapn + dhi*rhoi evapin = evapin + dhi*rhoi emlt_ocn = emlt_ocn - qmlt(k) * dhi !-------------------------------------------------------------- ! Melt ice (top) !-------------------------------------------------------------- if (qm(k) < c0) then dhi = max(-dzi(k), etop_mlt/qm(k)) else qm(k) = c0 dhi = -dzi(k) endif emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0 etop_mlt = max(etop_mlt - dhi*qm(k), c0) ! history diagnostics if (dhi < -puny .and. mlt_onset < puny) & mlt_onset = yday meltt = meltt - dhi enddo ! nilyr do k = nilyr, 1, -1 !-------------------------------------------------------------- ! Melt ice (bottom) !-------------------------------------------------------------- if (qm(k) < c0) then dhi = max(-dzi(k), ebot_mlt/qm(k)) else qm(k) = c0 dhi = -dzi(k) endif emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi dzi(k) = dzi(k) + dhi ! zqin < 0, dhi < 0 ebot_mlt = max(ebot_mlt - dhi*qm(k), c0) ! history diagnostics meltb = meltb -dhi enddo ! nilyr do k = nslyr, 1, -1 !-------------------------------------------------------------- ! Melt snow (only if all the ice has melted) !-------------------------------------------------------------- dhs = max(-dzs(k), ebot_mlt/zqsn(k)) dzs(k) = dzs(k) + dhs ! zqsn < 0, dhs < 0 ebot_mlt = ebot_mlt - dhs*zqsn(k) ebot_mlt = max(ebot_mlt, c0) ! Add this to the snow melt (J. Zhu) melts = melts - dhs enddo ! nslyr !----------------------------------------------------------------- ! Compute heat flux used by the ice (<=0). ! fhocn is the available ocean heat that is left after use by ice !----------------------------------------------------------------- fhocnn = fbot & + (esub + etop_mlt + ebot_mlt)/dt !---!----------------------------------------------------------------- !---! Add new snowfall at top surface. !---!----------------------------------------------------------------- !---------------------------------------------------------------- ! NOTE: If heat flux diagnostics are to work, new snow should ! have T = 0 (i.e. q = -rhos*Lfresh) and should not be ! converted to rain. !---------------------------------------------------------------- if (fsnow > c0) then hsn_new = fsnow/rhos * dt zqsnew = -rhos*Lfresh hstot = dzs(1) + hsn_new if (hstot > c0) then zqsn(1) = (dzs(1) * zqsn(1) & + hsn_new * zqsnew) / hstot ! avoid roundoff errors zqsn(1) = min(zqsn(1), -rhos*Lfresh) dzs(1) = hstot endif endif !----------------------------------------------------------------- ! Find the new ice and snow thicknesses. !----------------------------------------------------------------- hin = c0 hsn = c0 do k = 1, nilyr hin = hin + dzi(k) enddo ! k do k = 1, nslyr hsn = hsn + dzs(k) dsnow = dsnow + dzs(k) - hslyr enddo ! k !------------------------------------------------------------------- ! Convert snow to ice if snow lies below freeboard. !------------------------------------------------------------------- if (ktherm /= 2) & call freeboard (nslyr, & snoice, & hin, hsn, & zqin, zqsn, & dzi, dzs, & dsnow) if (icepack_warnings_aborted(subname)) return !---!------------------------------------------------------------------- !---! Repartition the ice and snow into equal-thickness layers, !---! conserving energy. !---!------------------------------------------------------------------- !----------------------------------------------------------------- ! Compute desired layer thicknesses. !----------------------------------------------------------------- if (hin > c0) then hilyr = hin / real(nilyr,kind=dbl_kind) else hin = c0 hilyr = c0 endif if (hsn > c0) then hslyr = hsn / real(nslyr,kind=dbl_kind) else hsn = c0 hslyr = c0 endif !----------------------------------------------------------------- ! Compute depths zi1 of old layers (unequal thickness). ! Compute depths zi2 of new layers (equal thickness). !----------------------------------------------------------------- zi1(1) = c0 zi1(1+nilyr) = hin zi2(1) = c0 zi2(1+nilyr) = hin if (heat_capacity) then do k = 1, nilyr-1 zi1(k+1) = zi1(k) + dzi(k) zi2(k+1) = zi2(k) + hilyr enddo !----------------------------------------------------------------- ! Conserving energy, compute the enthalpy of the new equal layers. !----------------------------------------------------------------- call adjust_enthalpy (nilyr, & zi1, zi2, & hilyr, hin, & zqin) if (icepack_warnings_aborted(subname)) return if (ktherm == 2) & call adjust_enthalpy (nilyr, & zi1, zi2, & hilyr, hin, & zSin) if (icepack_warnings_aborted(subname)) return else ! zero layer (nilyr=1) zqin(1) = -rhoi * Lfresh zqsn(1) = -rhos * Lfresh endif if (nslyr > 1) then !----------------------------------------------------------------- ! Compute depths zs1 of old layers (unequal thickness). ! Compute depths zs2 of new layers (equal thickness). !----------------------------------------------------------------- zs1(1) = c0 zs1(1+nslyr) = hsn zs2(1) = c0 zs2(1+nslyr) = hsn do k = 1, nslyr-1 zs1(k+1) = zs1(k) + dzs(k) zs2(k+1) = zs2(k) + hslyr enddo !----------------------------------------------------------------- ! Conserving energy, compute the enthalpy of the new equal layers. !----------------------------------------------------------------- call adjust_enthalpy (nslyr, & zs1, zs2, & hslyr, hsn, & zqsn) if (icepack_warnings_aborted(subname)) return endif ! nslyr > 1 !----------------------------------------------------------------- ! Remove very thin snow layers (ktherm = 2) !----------------------------------------------------------------- if (ktherm == 2) then do k = 1, nslyr if (hsn <= puny) then fhocnn = fhocnn & + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt) zqsn(k) = -rhos*Lfresh hslyr = c0 endif enddo endif !----------------------------------------------------------------- ! Compute final ice-snow energy, including the energy of ! sublimated/condensed ice. !----------------------------------------------------------------- efinal = -evapn*Lvap evapn = evapn/dt evapsn = evapsn/dt evapin = evapin/dt do k = 1, nslyr efinal = efinal + hslyr*zqsn(k) enddo do k = 1, nilyr efinal = efinal + hilyr*zqin(k) enddo ! k if (ktherm < 2) then emlt_atm = c0 emlt_ocn = c0 endif ! melt water is no longer zero enthalpy with ktherm=2 fhocnn = fhocnn + emlt_ocn/dt efinal = efinal + emlt_atm ! for conservation check end subroutine thickness_changes !======================================================================= ! ! If there is enough snow to lower the ice/snow interface below ! sea level, convert enough snow to ice to bring the interface back ! to sea level. ! ! authors William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL subroutine freeboard (nslyr, & snoice, & hin, hsn, & zqin, zqsn, & dzi, dzs, & dsnow) integer (kind=int_kind), intent(in) :: & nslyr ! number of snow layers ! real (kind=dbl_kind), intent(in) :: & ! dt ! time step real (kind=dbl_kind), & intent(inout) :: & snoice , & ! snow-ice formation (m/step-->cm/day) dsnow ! change in snow thickness after snow-ice formation (m) ! iage ! ice age (s) real (kind=dbl_kind), & intent(inout) :: & hin , & ! ice thickness (m) hsn ! snow thickness (m) real (kind=dbl_kind), dimension (:), intent(in) :: & zqsn ! snow layer enthalpy (J m-3) real (kind=dbl_kind), dimension (:), intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) dzi , & ! ice layer thicknesses (m) dzs ! snow layer thicknesses (m) ! local variables integer (kind=int_kind) :: & k ! vertical index real (kind=dbl_kind) :: & dhin , & ! change in ice thickness (m) dhsn , & ! change in snow thickness (m) hqs ! sum of h*q for snow (J m-2) real (kind=dbl_kind) :: & wk1 , & ! temporary variable dhs ! snow to remove from layer (m) character(len=*),parameter :: subname='(freeboard)' !----------------------------------------------------------------- ! Determine whether snow lies below freeboard. !----------------------------------------------------------------- dhin = c0 dhsn = c0 hqs = c0 wk1 = hsn - hin*(rhow-rhoi)/rhos if (wk1 > puny .and. hsn > puny) then ! snow below freeboard dhsn = min(wk1*rhoi/rhow, hsn) ! snow to remove dhin = dhsn * rhos/rhoi ! ice to add endif !----------------------------------------------------------------- ! Adjust snow layer thickness. ! Compute energy to transfer from snow to ice. !----------------------------------------------------------------- do k = nslyr, 1, -1 if (dhin > puny) then dhs = min(dhsn, dzs(k)) ! snow to remove from layer hsn = hsn - dhs dsnow = dsnow -dhs !new snow addition term dzs(k) = dzs(k) - dhs dhsn = dhsn - dhs dhsn = max(dhsn,c0) hqs = hqs + dhs * zqsn(k) endif ! dhin > puny enddo !----------------------------------------------------------------- ! Transfer volume and energy from snow to top ice layer. !----------------------------------------------------------------- if (dhin > puny) then ! update ice age due to freezing (new ice age = dt) ! if (tr_iage) & ! iage = (iage*hin+dt*dhin)/(hin+dhin) wk1 = dzi(1) + dhin hin = hin + dhin zqin(1) = (dzi(1)*zqin(1) + hqs) / wk1 dzi(1) = wk1 ! history diagnostic snoice = snoice + dhin endif ! dhin > puny end subroutine freeboard !======================================================================= ! ! Conserving energy, compute the new enthalpy of equal-thickness ice ! or snow layers. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine adjust_enthalpy (nlyr, & z1, z2, & hlyr, hn, & qn) integer (kind=int_kind), intent(in) :: & nlyr ! number of layers (nilyr or nslyr) real (kind=dbl_kind), dimension (:), intent(in) :: & z1 , & ! interface depth for old, unequal layers (m) z2 ! interface depth for new, equal layers (m) real (kind=dbl_kind), intent(in) :: & hlyr ! new layer thickness (m) real (kind=dbl_kind), intent(in) :: & hn ! total thickness (m) real (kind=dbl_kind), dimension (:), intent(inout) :: & qn ! layer quantity (enthalpy, salinity...) ! local variables integer (kind=int_kind) :: & k, k1, k2 ! vertical indices real (kind=dbl_kind) :: & hovlp ! overlap between old and new layers (m) real (kind=dbl_kind) :: & rhlyr ! 1./hlyr real (kind=dbl_kind), dimension (nlyr) :: & hq ! h * q for a layer character(len=*),parameter :: subname='(adjust_enthalpy)' !----------------------------------------------------------------- ! Compute reciprocal layer thickness. !----------------------------------------------------------------- rhlyr = c0 if (hn > puny) rhlyr = c1 / hlyr !----------------------------------------------------------------- ! Compute h*q for new layers (k2) given overlap with old layers (k1) !----------------------------------------------------------------- do k2 = 1, nlyr hq(k2) = c0 enddo ! k k1 = 1 k2 = 1 do while (k1 <= nlyr .and. k2 <= nlyr) hovlp = min (z1(k1+1), z2(k2+1)) & - max (z1(k1), z2(k2)) hovlp = max (hovlp, c0) hq(k2) = hq(k2) + hovlp*qn(k1) if (z1(k1+1) > z2(k2+1)) then k2 = k2 + 1 else k1 = k1 + 1 endif enddo ! while !----------------------------------------------------------------- ! Compute new enthalpies. !----------------------------------------------------------------- do k = 1, nlyr qn(k) = hq(k) * rhlyr enddo ! k end subroutine adjust_enthalpy !======================================================================= ! ! Check for energy conservation by comparing the change in energy ! to the net energy input. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW ! Adrian K. Turner, LANL subroutine conservation_check_vthermo(dt, & fsurfn, flatn, & fhocnn, fswint, & fsnow, & einit, einter, & efinal, & fcondtopn,fcondbotn, & fadvocn, fbot ) real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), intent(in) :: & fsurfn , & ! net flux to top surface, excluding fcondtopn flatn , & ! surface downward latent heat (W m-2) fhocnn , & ! fbot, corrected for any surplus energy fswint , & ! SW absorbed in ice interior, below surface (W m-2) fsnow , & ! snowfall rate (kg m-2 s-1) fcondtopn , & fadvocn , & fbot real (kind=dbl_kind), intent(in) :: & einit , & ! initial energy of melting (J m-2) einter , & ! intermediate energy of melting (J m-2) efinal , & ! final energy of melting (J m-2) fcondbotn ! local variables real (kind=dbl_kind) :: & einp , & ! energy input during timestep (J m-2) ferr ! energy conservation error (W m-2) character(len=*),parameter :: subname='(conservation_check_vthermo)' !---------------------------------------------------------------- ! If energy is not conserved, print diagnostics and exit. !---------------------------------------------------------------- !----------------------------------------------------------------- ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent ! heat is not included in the energy input, since (efinal - einit) ! is the energy change in the system ice + vapor, and the latent ! heat lost by the ice is equal to that gained by the vapor. !----------------------------------------------------------------- einp = (fsurfn - flatn + fswint - fhocnn & - fsnow*Lfresh - fadvocn) * dt ferr = abs(efinal-einit-einp) / dt if (ferr > ferrmax) then call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" ) write(warnstr,*) subname, 'Thermo energy conservation error' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Energy error (J) =', ferr*dt call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Initial energy =', einit call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final energy =', efinal call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'efinal - einit =', efinal-einit call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fsurfn,flatn,fswint,fhocn, fsnow*Lfresh:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn, fsnow*Lfresh call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Input energy =', einp call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fbot,fcondbot:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fbot,fcondbotn call icepack_warnings_add(warnstr) ! if (ktherm == 2) then write(warnstr,*) subname, 'Intermediate energy =', einter call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'efinal - einter =', & efinal-einter call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'einter - einit =', & einter-einit call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Conduction Error =', (einter-einit) & - (fcondtopn*dt - fcondbotn*dt + fswint*dt) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Melt/Growth Error =', (einter-einit) & + ferr*dt - (fcondtopn*dt - fcondbotn*dt + fswint*dt) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Advection Error =', fadvocn*dt call icepack_warnings_add(warnstr) ! endif ! write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn ! call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dt*(fsurfn, flatn, fswint, fhocn, fsnow*Lfresh, fadvocn):' call icepack_warnings_add(warnstr) write(warnstr,*) subname, fsurfn*dt, flatn*dt, & fswint*dt, fhocnn*dt, & fsnow*Lfresh*dt, fadvocn*dt call icepack_warnings_add(warnstr) return endif end subroutine conservation_check_vthermo !======================================================================= ! ! Given the vertical thermo state variables (hin, hsn), ! compute the new ice state variables (vicen, vsnon). ! Zero out state variables if ice has melted entirely. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW ! Elizabeth C. Hunke, LANL subroutine update_state_vthermo(nilyr, nslyr, & Tf, Tsf, & hin, hsn, & zqin, zSin, & zqsn, & aicen, vicen, & vsnon) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & Tf ! freezing temperature (C) real (kind=dbl_kind), intent(inout) :: & Tsf ! ice/snow surface temperature, Tsfcn real (kind=dbl_kind), intent(in) :: & hin , & ! ice thickness (m) hsn ! snow thickness (m) real (kind=dbl_kind), dimension (:), intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) zSin , & ! ice salinity (ppt) zqsn ! snow layer enthalpy (J m-3) real (kind=dbl_kind), intent(inout) :: & aicen , & ! concentration of ice vicen , & ! volume per unit area of ice (m) vsnon ! volume per unit area of snow (m) ! local variables integer (kind=int_kind) :: & k ! ice layer index character(len=*),parameter :: subname='(update_state_vthermo)' if (hin <= c0) then aicen = c0 vicen = c0 vsnon = c0 Tsf = Tf do k = 1, nilyr zqin(k) = c0 enddo if (ktherm == 2) then do k = 1, nilyr zSin(k) = c0 enddo endif do k = 1, nslyr zqsn(k) = c0 enddo else ! aicen is already up to date vicen = aicen * hin vsnon = aicen * hsn endif end subroutine update_state_vthermo !======================================================================= !autodocument_start icepack_step_therm1 ! Driver for thermodynamic changes not needed for coupling: ! transport in thickness space, lateral growth and melting. ! ! authors: William H. Lipscomb, LANL ! Elizabeth C. Hunke, LANL subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & aicen_init , & vicen_init , vsnon_init , & aice , aicen , & vice , vicen , & vsno , vsnon , & uvel , vvel , & Tsfc , zqsn , & zqin , zSin , & alvl , vlvl , & apnd , hpnd , & ipnd , & iage , FY , & aerosno , aeroice , & isosno , isoice , & uatm , vatm , & wind , zlvl , & Qa , rhoa , & Qa_iso , & Tair , Tref , & Qref , Uref , & Qref_iso , & Cdn_atm_ratio, & 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 , & strax , stray , & strairxT , strairyT , & potT , sst , & sss , Tf , & strocnxT , strocnyT , & fbot , & Tbot , Tsnice , & frzmlt , rside , & fside , & fsnow , frain , & fpond , & fsurf , fsurfn , & fcondtop , fcondtopn , & fcondbot , fcondbotn , & fswsfcn , fswintn , & fswthrun , & fswthrun_vdr, & fswthrun_vdf, & fswthrun_idr, & fswthrun_idf, & fswabs , & flwout , & Sswabsn , Iswabsn , & flw , & fsens , fsensn , & flat , flatn , & evap , & evaps , evapi , & fresh , fsalt , & fhocn , & fswthru , & fswthru_vdr , & fswthru_vdf , & fswthru_idr , & fswthru_idf , & flatn_f , fsensn_f , & fsurfn_f , fcondtopn_f , & faero_atm , faero_ocn , & fiso_atm , fiso_ocn , & fiso_evap , & HDO_ocn , H2_16O_ocn , & H2_18O_ocn , & dhsn , ffracn , & meltt , melttn , & meltb , meltbn , & melts , meltsn , & congel , congeln , & snoice , snoicen , & dsnown , & lmask_n , lmask_s , & mlt_onset , frz_onset , & yday , prescribed_ice,& zlvs) 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), intent(in) :: & dt , & ! time step uvel , & ! x-component of velocity (m/s) vvel , & ! y-component of velocity (m/s) strax , & ! wind stress components (N/m^2) stray , & ! yday ! day of year logical (kind=log_kind), intent(in) :: & lmask_n , & ! northern hemisphere mask lmask_s ! southern hemisphere mask logical (kind=log_kind), intent(in), optional :: & prescribed_ice ! if .true., use prescribed ice instead of computed real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration vice , & ! volume per unit area of ice (m) vsno , & ! volume per unit area of snow (m) zlvl , & ! atm level height (m) uatm , & ! wind velocity components (m/s) vatm , & wind , & ! wind speed (m/s) potT , & ! air potential temperature (K) Tair , & ! air temperature (K) Qa , & ! specific humidity (kg/kg) rhoa , & ! air density (kg/m^3) frain , & ! rainfall rate (kg/m^2 s) fsnow , & ! snowfall rate (kg/m^2 s) 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) fswthru , & ! shortwave penetrating to ocean (W/m^2) fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2) fcondtop , & ! top surface conductive flux (W/m^2) fcondbot , & ! bottom surface conductive flux (W/m^2) fsens , & ! sensible heat flux (W/m^2) flat , & ! latent heat flux (W/m^2) fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2) flw , & ! incoming longwave radiation (W/m^2) flwout , & ! outgoing longwave radiation (W/m^2) evap , & ! evaporative water flux (kg/m^2/s) evaps , & ! evaporative water flux over snow (kg/m^2/s) evapi , & ! evaporative water flux over ice (kg/m^2/s) congel , & ! basal ice growth (m/step-->cm/day) snoice , & ! snow-ice formation (m/step-->cm/day) Tref , & ! 2m atm reference temperature (K) Qref , & ! 2m atm reference spec humidity (kg/kg) Uref , & ! 10m atm reference wind speed (m/s) Cdn_atm , & ! atm drag coefficient Cdn_ocn , & ! ocn drag coefficient hfreebd , & ! freeboard (m) hdraft , & ! draft of ice + snow column (Stoessel1993) hridge , & ! ridge height distrdg , & ! distance between ridges hkeel , & ! keel depth dkeel , & ! distance between keels lfloe , & ! floe length dfloe , & ! distance between floes Cdn_atm_skin, & ! neutral skin drag coefficient Cdn_atm_floe, & ! neutral floe edge drag coefficient Cdn_atm_pond, & ! neutral pond edge drag coefficient Cdn_atm_rdg , & ! neutral ridge drag coefficient Cdn_ocn_skin, & ! skin drag coefficient Cdn_ocn_floe, & ! floe edge drag coefficient Cdn_ocn_keel, & ! keel drag coefficient Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm strairxT , & ! stress on ice by air, x-direction strairyT , & ! stress on ice by air, y-direction strocnxT , & ! ice-ocean stress, x-direction strocnyT , & ! ice-ocean stress, y-direction fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) frzmlt , & ! freezing/melting potential (W/m^2) rside , & ! fraction of ice that melts laterally fside , & ! lateral heat flux (W/m^2) sst , & ! sea surface temperature (C) Tf , & ! freezing temperature (C) Tbot , & ! ice bottom surface temperature (deg C) Tsnice , & ! snow ice interface temperature (deg C) sss , & ! sea surface salinity (ppt) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) meltb , & ! basal ice melt (m/step-->cm/day) mlt_onset , & ! day of year that sfc melting begins frz_onset ! day of year that freezing begins (congel or frazil) real (kind=dbl_kind), intent(inout), optional :: & fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2) fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) fswthru_idf ! nir dif shortwave penetrating to ocean (W/m^2) real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & Qa_iso , & ! isotope specific humidity (kg/kg) Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg) fiso_atm , & ! isotope deposition rate (kg/m^2 s) fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) fiso_evap ! isotope evaporation (kg/m^2/s) real (kind=dbl_kind), optional, intent(in) :: & HDO_ocn , & ! ocean concentration of HDO (kg/kg) H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice vicen_init , & ! volume per unit area of ice (m) vsnon_init , & ! volume per unit area of snow (m) aicen , & ! concentration of ice vicen , & ! volume per unit area of ice (m) vsnon , & ! volume per unit area of snow (m) Tsfc , & ! ice/snow surface temperature, Tsfcn alvl , & ! level ice area fraction vlvl , & ! level ice volume fraction apnd , & ! melt pond area fraction hpnd , & ! melt pond depth (m) ipnd , & ! melt pond refrozen lid thickness (m) iage , & ! volume-weighted ice age FY , & ! area-weighted first-year ice area fsurfn , & ! net flux to top surface, excluding fcondtop fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) flatn , & ! latent heat flux (W m-2) fsensn , & ! sensible heat flux (W m-2) fsurfn_f , & ! net flux to top surface, excluding fcondtop fcondtopn_f , & ! downward cond flux at top surface (W m-2) flatn_f , & ! latent heat flux (W m-2) fsensn_f , & ! sensible heat flux (W m-2) fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) fswthrun , & ! SW through ice to ocean (W/m^2) fswintn , & ! SW absorbed in ice interior, below surface (W m-2) faero_atm , & ! aerosol deposition rate (kg/m^2 s) faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond meltsn , & ! snow melt (m) melttn , & ! top ice melt (m) meltbn , & ! bottom ice melt (m) congeln , & ! congelation ice growth (m) snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) real (kind=dbl_kind), optional, dimension(:), intent(inout) :: & fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind), dimension(:,:), intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zqin , & ! ice layer enthalpy (J m-3) zSin , & ! internal ice layer salinities Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) Iswabsn ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & aerosno , & ! snow aerosol tracer (kg/m^2) aeroice ! ice aerosol tracer (kg/m^2) real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: & isosno , & ! snow isotope tracer (kg/m^2) isoice ! ice isotope tracer (kg/m^2) !autodocument_end ! local variables integer (kind=int_kind) :: & n ! category index real (kind=dbl_kind) :: & worka, workb ! temporary variables ! 2D coupler variables (computed for each category, then aggregated) real (kind=dbl_kind) :: & fswabsn , & ! shortwave absorbed by ice (W/m^2) flwoutn , & ! upward LW at surface (W/m^2) evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1) evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1) evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1) freshn , & ! flux of water, ice to ocean (kg/m^2/s) fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s) fhocnn , & ! fbot corrected for leftover energy (W/m^2) strairxn , & ! air/ice zonal stress, (N/m^2) strairyn , & ! air/ice meridional stress, (N/m^2) Cdn_atm_ratio_n, & ! drag coefficient ratio Trefn , & ! air tmp reference level (K) Urefn , & ! air speed reference level (m/s) Qrefn , & ! air sp hum reference level (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef , & ! transfer coefficient for latent heat rfrac ! water fraction retained for melt ponds real (kind=dbl_kind), dimension(n_iso) :: & Qrefn_iso , & ! isotope air sp hum reference level (kg/kg) fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s) fiso_evapn ! isotope evaporation (kg/m^2/s) real (kind=dbl_kind), allocatable, dimension(:,:) :: & l_isosno , & ! local snow isotope tracer (kg/m^2) l_isoice ! local ice isotope tracer (kg/m^2) real (kind=dbl_kind), allocatable, dimension(:) :: & l_Qa_iso , & ! local isotope specific humidity (kg/kg) l_Qref_iso , & ! local isotope 2m atm reference spec humidity (kg/kg) l_fiso_atm , & ! local isotope deposition rate (kg/m^2 s) l_fiso_ocn , & ! local isotope flux to ocean (kg/m^2/s) l_fiso_evap ! local isotope evaporation (kg/m^2/s) real (kind=dbl_kind) :: & l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg) l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg) l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg) real (kind=dbl_kind) :: & l_fswthru_vdr , & ! vis dir SW through ice to ocean (W/m^2) l_fswthru_vdf , & ! vis dif SW through ice to ocean (W/m^2) l_fswthru_idr , & ! nir dir SW through ice to ocean (W/m^2) l_fswthru_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind), dimension(:), allocatable :: & l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind) :: & pond ! water retained in ponds (m) character(len=*),parameter :: subname='(icepack_step_therm1)' !----------------------------------------------------------------- ! allocate local optional arguments !----------------------------------------------------------------- if (present(isosno) ) then allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2))) l_isosno = isosno else allocate(l_isosno(1,1)) l_isosno = c0 endif if (present(isoice) ) then allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2))) l_isoice = isoice else allocate(l_isoice(1,1)) l_isoice = c0 endif if (present(Qa_iso) ) then allocate(l_Qa_iso(size(Qa_iso))) l_Qa_iso = Qa_iso else allocate(l_Qa_iso(1)) l_Qa_iso = c0 endif if (present(Qref_iso) ) then allocate(l_Qref_iso(size(Qref_iso))) l_Qref_iso = Qref_iso else allocate(l_Qref_iso(1)) l_Qref_iso = c0 endif if (present(fiso_atm) ) then allocate(l_fiso_atm(size(fiso_atm))) l_fiso_atm = fiso_atm else allocate(l_fiso_atm(1)) l_fiso_atm = c0 endif if (present(fiso_ocn) ) then allocate(l_fiso_ocn(size(fiso_ocn))) l_fiso_ocn = fiso_ocn else allocate(l_fiso_ocn(1)) l_fiso_ocn = c0 endif if (present(fiso_evap) ) then allocate(l_fiso_evap(size(fiso_evap))) l_fiso_evap = fiso_evap else allocate(l_fiso_evap(1)) l_fiso_evap = c0 endif l_HDO_ocn = c0 if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn l_H2_16O_ocn = c0 if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn l_H2_18O_ocn = c0 if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn l_fswthru_vdr = c0 if (present(fswthru_vdr) ) l_fswthru_vdr = fswthru_vdr l_fswthru_vdf = c0 if (present(fswthru_vdf) ) l_fswthru_vdf = fswthru_vdf l_fswthru_idr = c0 if (present(fswthru_idr) ) l_fswthru_idr = fswthru_idr l_fswthru_idf = c0 if (present(fswthru_idf) ) l_fswthru_idf = fswthru_idf allocate(l_fswthrun_vdr(ncat)) allocate(l_fswthrun_vdf(ncat)) allocate(l_fswthrun_idr(ncat)) allocate(l_fswthrun_idf(ncat)) l_fswthrun_vdr = c0 if (present(fswthrun_vdr) ) l_fswthrun_vdr = fswthrun_vdr l_fswthrun_vdf = c0 if (present(fswthrun_vdf) ) l_fswthrun_vdf = fswthrun_vdf l_fswthrun_idr = c0 if (present(fswthrun_idr) ) l_fswthrun_idr = fswthrun_idr l_fswthrun_idf = c0 if (present(fswthrun_idf) ) l_fswthrun_idf = fswthrun_idf !----------------------------------------------------------------- ! Adjust frzmlt to account for ice-ocean heat fluxes since last ! call to coupler. ! Compute lateral and bottom heat fluxes. !----------------------------------------------------------------- call frzmlt_bottom_lateral (dt, ncat, & nilyr, nslyr, & aice, frzmlt, & vicen, vsnon, & zqin, zqsn, & sst, Tf, & ustar_min, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & rside, Cdn_ocn, & fside) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Update the neutral drag coefficients to account for form drag ! Oceanic and atmospheric drag coefficients !----------------------------------------------------------------- if (formdrag) then call neutral_drag_coeffs (apnd , & hpnd , ipnd , & alvl , vlvl , & aice , vice, & vsno , aicen , & vicen , & 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) if (icepack_warnings_aborted(subname)) return endif do n = 1, ncat meltsn (n) = c0 melttn (n) = c0 meltbn (n) = c0 congeln(n) = c0 snoicen(n) = c0 dsnown (n) = c0 Trefn = c0 Qrefn = c0 Qrefn_iso(:) = c0 fiso_ocnn(:) = c0 fiso_evapn(:) = c0 Urefn = c0 lhcoef = c0 shcoef = c0 worka = c0 workb = c0 if (aicen_init(n) > puny) then if (calc_Tsfc .or. calc_strair) then !----------------------------------------------------------------- ! Atmosphere boundary layer calculation; compute coefficients ! for sensible and latent heat fluxes. ! ! NOTE: The wind stress is computed here for later use if ! calc_strair = .true. Otherwise, the wind stress ! components are set to the data values. !----------------------------------------------------------------- call icepack_atm_boundary('ice', & Tsfc(n), potT, & uatm, vatm, & wind, zlvl, & Qa, rhoa, & strairxn, strairyn, & Trefn, Qrefn, & worka, workb, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & Qa_iso=l_Qa_iso, & Qref_iso=Qrefn_iso, & uvel=uvel, vvel=vvel, & Uref=Urefn, zlvs=zlvs) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc or calc_strair if (.not.(calc_strair)) then #ifndef CICE_IN_NEMO ! Set to data values (on T points) strairxn = strax strairyn = stray #else ! NEMO wind stress is supplied on u grid, multipied ! by ice concentration and set directly in evp, so ! strairxT/yT = 0. Zero u-components here for safety. strairxn = c0 strairyn = c0 #endif endif !----------------------------------------------------------------- ! Update ice age ! This is further adjusted for freezing in the thermodynamics. ! Melting does not alter the ice age. !----------------------------------------------------------------- if (tr_iage) call increment_age (dt, iage(n)) if (icepack_warnings_aborted(subname)) return if (tr_FY) call update_FYarea (dt, & lmask_n, lmask_s, & yday, FY(n)) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Vertical thermodynamics: Heat conduction, growth and melting. !----------------------------------------------------------------- if (.not.(calc_Tsfc)) then ! If not calculating surface temperature and fluxes, set ! surface fluxes (flatn, fsurfn, and fcondtopn) to be used ! in thickness_changes ! hadgem routine sets fluxes to default values in ice-only mode call set_sfcflux(aicen (n), & flatn_f (n), fsensn_f (n), & fcondtopn_f(n), & fsurfn_f (n), & flatn (n), fsensn (n), & fsurfn (n), & fcondtopn (n)) if (icepack_warnings_aborted(subname)) return endif call thermo_vertical(nilyr, nslyr, & dt, aicen (n), & vicen (n), vsnon (n), & Tsfc (n), zSin (:,n), & zqin (:,n), zqsn (:,n), & apnd (n), hpnd (n), & tr_pond_topo, & flw, potT, & Qa, rhoa, & fsnow, fpond, & fbot, Tbot, & Tsnice, sss, & lhcoef, shcoef, & fswsfcn (n), fswintn (n), & Sswabsn(:,n), Iswabsn(:,n), & fsurfn (n), fcondtopn(n), & fcondbotn(n), & fsensn (n), flatn (n), & flwoutn, evapn, & evapsn, evapin, & freshn, fsaltn, & fhocnn, & melttn (n), meltsn (n), & meltbn (n), & congeln (n), snoicen (n), & mlt_onset, frz_onset, & yday, dsnown (n), & prescribed_ice) if (icepack_warnings_aborted(subname)) then call icepack_warnings_add(subname//' ice: Vertical thermo error: ') return endif !----------------------------------------------------------------- ! Total absorbed shortwave radiation !----------------------------------------------------------------- fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n) !----------------------------------------------------------------- ! Aerosol update !----------------------------------------------------------------- if (tr_aero) then call update_aerosol (dt, & nilyr, nslyr, n_aero, & melttn (n), meltsn (n), & meltbn (n), congeln (n), & snoicen (n), fsnow, & aerosno(:,:,n), aeroice(:,:,n), & aicen_init (n), vicen_init (n), & vsnon_init (n), & vicen (n), vsnon (n), & aicen (n), & faero_atm , faero_ocn) if (icepack_warnings_aborted(subname)) return endif if (tr_iso) then call update_isotope (dt = dt, & nilyr = nilyr, nslyr = nslyr, & meltt = melttn(n),melts = meltsn(n), & meltb = meltbn(n),congel=congeln(n), & snoice=snoicen(n),evap=evapn, & fsnow=fsnow, Tsfc=Tsfc(n), & Qref_iso=Qrefn_iso(:), & isosno=l_isosno(:,n),isoice=l_isoice(:,n), & aice_old=aicen_init(n),vice_old=vicen_init(n), & vsno_old=vsnon_init(n), & vicen=vicen(n),vsnon=vsnon(n), & aicen=aicen(n), & fiso_atm=l_fiso_atm(:), & fiso_evapn=fiso_evapn(:), & fiso_ocnn=fiso_ocnn(:), & HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, & H2_18O_ocn=l_H2_18O_ocn) if (icepack_warnings_aborted(subname)) return endif endif ! aicen_init !----------------------------------------------------------------- ! Melt ponds ! If using tr_pond_cesm, the full calculation is performed here. ! If using tr_pond_topo, the rest of the calculation is done after ! the surface fluxes are merged, below. !----------------------------------------------------------------- !call ice_timer_start(timer_ponds) if (tr_pond) then if (tr_pond_cesm) then rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n) call compute_ponds_cesm(dt, hi_min, & pndaspect, rfrac, & melttn(n), meltsn(n), & frain, & aicen (n), vicen (n), & Tsfc (n), & apnd (n), hpnd (n)) if (icepack_warnings_aborted(subname)) return elseif (tr_pond_lvl) then rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n) call compute_ponds_lvl(dt, nilyr, & ktherm, & hi_min, & dpscale, frzpnd, & pndaspect, rfrac, & melttn(n), meltsn(n), & frain, Tair, & fsurfn(n), & dhsn (n), ffracn(n), & aicen (n), vicen (n), & vsnon (n), & zqin(:,n), zSin(:,n), & Tsfc (n), alvl (n), & apnd (n), hpnd (n), & ipnd (n)) if (icepack_warnings_aborted(subname)) return elseif (tr_pond_topo) then if (aicen_init(n) > puny) then ! collect liquid water in ponds ! assume salt still runs off rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n) pond = rfrac/rhofresh * (melttn(n)*rhoi & + meltsn(n)*rhos & + frain *dt) ! if pond does not exist, create new pond over full ice area ! otherwise increase pond depth without changing pond area if (apnd(n) < puny) then hpnd(n) = c0 apnd(n) = c1 endif hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n) fpond = fpond + pond * aicen(n) ! m endif ! aicen_init endif endif ! tr_pond !call ice_timer_stop(timer_ponds) !----------------------------------------------------------------- ! Increment area-weighted fluxes. !----------------------------------------------------------------- if (aicen_init(n) > puny) & call merge_fluxes (aicen=aicen_init(n), & flw=flw, & strairxn=strairxn, strairyn=strairyn,& Cdn_atm_ratio_n=Cdn_atm_ratio_n, & fsurfn=fsurfn(n), fcondtopn=fcondtopn(n),& fcondbotn=fcondbotn(n), & fsensn=fsensn(n), flatn=flatn(n), & fswabsn=fswabsn, flwoutn=flwoutn, & evapn=evapn, & evapsn=evapsn, evapin=evapin, & Trefn=Trefn, Qrefn=Qrefn, & freshn=freshn, fsaltn=fsaltn, & fhocnn=fhocnn, & fswthrun=fswthrun(n), & fswthrun_vdr=l_fswthrun_vdr(n), & fswthrun_vdf=l_fswthrun_vdf(n), & fswthrun_idr=l_fswthrun_idr(n), & fswthrun_idf=l_fswthrun_idf(n), & strairxT=strairxT, strairyT=strairyT,& Cdn_atm_ratio=Cdn_atm_ratio, & fsurf=fsurf, fcondtop=fcondtop,& fcondbot=fcondbot, & fsens=fsens, flat=flat, & fswabs=fswabs, flwout=flwout, & evap=evap, & evaps=evaps, evapi=evapi, & Tref=Tref, Qref=Qref, & fresh=fresh, fsalt=fsalt, & fhocn=fhocn, & fswthru=fswthru, & fswthru_vdr=l_fswthru_vdr, & fswthru_vdf=l_fswthru_vdf, & fswthru_idr=l_fswthru_idr, & fswthru_idf=l_fswthru_idf, & melttn=melttn (n), meltsn=meltsn(n), & meltbn=meltbn (n), congeln=congeln(n),& snoicen=snoicen(n), & meltt=meltt, melts=melts, & meltb=meltb, congel=congel, & snoice=snoice, & Uref=Uref, Urefn=Urefn, & Qref_iso=l_Qref_iso, & Qrefn_iso=Qrefn_iso, & fiso_ocn=l_fiso_ocn, & fiso_ocnn=fiso_ocnn, & fiso_evap=l_fiso_evap, & fiso_evapn=fiso_evapn) if (icepack_warnings_aborted(subname)) return enddo ! ncat if (present(isosno) ) isosno = l_isosno if (present(isoice) ) isoice = l_isoice if (present(Qa_iso) ) Qa_iso = l_Qa_iso if (present(Qref_iso) ) Qref_iso = l_Qref_iso if (present(fiso_atm) ) fiso_atm = l_fiso_atm if (present(fiso_ocn) ) fiso_ocn = l_fiso_ocn if (present(fiso_evap)) fiso_evap= l_fiso_evap if (present(fswthrun_vdr)) fswthrun_vdr= l_fswthrun_vdr if (present(fswthrun_vdf)) fswthrun_vdf= l_fswthrun_vdf if (present(fswthrun_idr)) fswthrun_idr= l_fswthrun_idr if (present(fswthrun_idf)) fswthrun_idf= l_fswthrun_idf if (present(fswthru_vdr)) fswthru_vdr= l_fswthru_vdr if (present(fswthru_vdf)) fswthru_vdf= l_fswthru_vdf if (present(fswthru_idr)) fswthru_idr= l_fswthru_idr if (present(fswthru_idf)) fswthru_idf= l_fswthru_idf deallocate(l_isosno) deallocate(l_isoice) deallocate(l_Qa_iso) deallocate(l_Qref_iso) deallocate(l_fiso_atm) deallocate(l_fiso_ocn) deallocate(l_fiso_evap) deallocate(l_fswthrun_vdr) deallocate(l_fswthrun_vdf) deallocate(l_fswthrun_idr) deallocate(l_fswthrun_idf) !----------------------------------------------------------------- ! Calculate ponds from the topographic scheme !----------------------------------------------------------------- !call ice_timer_start(timer_ponds) if (tr_pond_topo) then call compute_ponds_topo(dt, ncat, nilyr, & ktherm, heat_capacity, & aice, aicen, & vice, vicen, & vsno, vsnon, & meltt, & fsurf, fpond, & Tsfc, Tf, & zqin, zSin, & apnd, hpnd, ipnd ) if (icepack_warnings_aborted(subname)) return endif !call ice_timer_stop(timer_ponds) end subroutine icepack_step_therm1 !======================================================================= end module icepack_therm_vertical !=======================================================================