My Project
Functions/Subroutines | Variables
ice_therm_vertical Module Reference

Functions/Subroutines

subroutine init_thermo_vertical
 
subroutine thermo_vertical
 
subroutine frzmlt_bottom_lateral (Tbot, fbot, rside)
 
subroutine init_thermo_vars (strxn, stryn, Trefn, Qrefn, fsensn, flatn, fswabsn, flwoutn, evapn, freshn, fsaltn, fhnetn, fswthrun, fsurf, fcondtop, fcondbot, fswint, einit, efinal, mvap)
 
subroutine init_vertical_profile (ni, icells, indxi, indxj, hin, hlyr, hsn, hin_init, hsn_init, qin, Tin, qsn, Tsn, Tsf, einit)
 
subroutine temperature_changes (ni, icells, indxi, indxj, hlyr, hsn, qin, Tin, qsn, Tsn, Tsf, Tbot, fbot, fsensn, flatn, fswabsn, flwoutn, fswthrun, fhnetn, fsurf, fcondtop, fcondbot, fswint, einit)
 
subroutine conductivity (icells, indxi, indxj, hlyr, hsn, Tin, Tbot, khi, khs)
 
subroutine absorbed_solar (ni, icells, indxi, indxj, hlyr, hsn, fswsfc, fswint, fswthrun, Iabs)
 
subroutine surface_fluxes (isolve, indxii, indxjj, Tsf, fswsfc, flwoutn, fsensn, flatn, fsurf, dflwout_dT, dfsens_dT, dflat_dT, dfsurf_dT)
 
subroutine tridiag_solver (isolve, indxii, indxjj, nmat, sbdiag, diag, spdiag, rhs, xout)
 
subroutine thickness_changes (ni, icells, indxi, indxj, hin, hlyr, hsn, qin, qsn, mvap, Tbot, fbot, flatn, fhnetn, fsurf, fcondtop, fcondbot, efinal)
 
subroutine conservation_check_vthermo (ni, icells, indxi, indxj, fsurf, flatn, fhnetn, fswint, einit, efinal)
 
subroutine add_new_snow (ni, icells, indxi, indxj, hsn, hsn_new, qsn, Tsf)
 
subroutine update_state_vthermo (ni, icells, indxi, indxj, hin, hsn, qin, qsn, Tsf)
 

Variables

real(kind=dbl_kind), parameter ferrmax = 1.0e-3_dbl_kind
 
real(kind=dbl_kind), parameter hsnomin = 1.0e-6_dbl_kind
 
real(kind=dbl_kind), dimension(nilyr+1) salin
 
real(kind=dbl_kind), dimension(nilyr+1) tmlt
 
real(kind=dbl_kind) ustar_scale
 
real(kind=dbl_kind), dimension(:,:,:), allocatable, save hicen_old
 
real(kind=dbl_kind), dimension(:,:), allocatable, save rside
 
logical(kind=log_kind) l_brine
 

Function/Subroutine Documentation

◆ absorbed_solar()

subroutine ice_therm_vertical::absorbed_solar ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  hlyr,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  hsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout)  fswsfc,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout)  fswint,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout)  fswthrun,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(inout)  Iabs 
)

Definition at line 1996 of file ice_therm_vertical.f90.

1996 !
1997 ! !USES:
1998 !
1999  use ice_albedo
2000 !
2001 ! !INPUT/OUTPUT PARAMETERS:
2002 !
2003  integer (kind=int_kind), intent(in) :: &
2004  ni & ! thickness category index
2005  , icells ! number of cells with aicen > puny
2006 
2007  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2008  intent(in) :: &
2009  indxi, indxj ! compressed indices for cells with aicen > puny
2010 
2011  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: &
2012  hlyr & ! ice layer thickness
2013  , hsn ! snow thickness (m)
2014 
2015  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout) :: &
2016  fswsfc & ! SW absorbed at ice/snow surface (W m-2)
2017  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
2018  , fswthrun ! SW through ice to ocean (W m-2)
2019 
2020  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
2021  intent(inout) :: &
2022  iabs ! SW absorbed in particular layer (W m-2)
2023 !
2024 !EOP
2025 !
2026  real (kind=dbl_kind), parameter :: &
2027  i0vis = 0.70_dbl_kind ! fraction of penetrating solar rad (visible)
2028 
2029  integer (kind=int_kind) :: &
2030  i, j & ! horizontal indices
2031  , ij & ! horizontal index, combines i and j loops
2032  , k ! ice layer index
2033 
2034  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2035  fswpen & ! SW penetrating beneath surface (W m-2)
2036  , trantop & ! transmitted frac of penetrating SW at layer top
2037  , tranbot ! transmitted frac of penetrating SW at layer bot
2038 
2039  real (kind=dbl_kind) :: &
2040  swabs & ! net SW down at surface (W m-2)
2041  , swabsv & ! swabs in vis (wvlngth < 700nm) (W/m^2)
2042  , swabsi & ! swabs in nir (wvlngth > 700nm) (W/m^2)
2043  , frsnow ! fractional snow coverage
2044 
2045 
2046  fswpen(:,:) = c0i
2047  trantop(:,:) = c0i
2048  tranbot(:,:) = c0i
2049 
2050  do ij = 1, icells
2051  i = indxi(ij)
2052  j = indxj(ij)
2053 
2054  !-----------------------------------------------------------------
2055  ! Fractional snow cover
2056  !-----------------------------------------------------------------
2057 
2058  if (hsn(i,j) > puny) then
2059  frsnow = hsn(i,j) / (hsn(i,j) + snowpatch)
2060  else
2061  frsnow = c0i
2062  endif
2063 
2064  !-----------------------------------------------------------------
2065  ! Shortwave flux absorbed at surface, absorbed internally,
2066  ! and penetrating to mixed layer.
2067  ! This parameterization assumes that all IR is absorbed at the
2068  ! surface; only visible is absorbed in the ice interior or
2069  ! transmitted to the ocean.
2070  !-----------------------------------------------------------------
2071 
2072  swabsv = swvdr(i,j)*(c1i-alvdrn(i,j,ni)) &
2073  + swvdf(i,j)*(c1i-alvdfn(i,j,ni))
2074  swabsi = swidr(i,j)*(c1i-alidrn(i,j,ni)) &
2075  + swidf(i,j)*(c1i-alidfn(i,j,ni))
2076  swabs = swabsv + swabsi
2077 
2078  fswpen(i,j) = swabsv * (c1i-frsnow) * i0vis
2079 !c & + swabsi * (c1i-frsnow) * i0nir ! i0nir = 0
2080  fswsfc(i,j) = swabs - fswpen(i,j)
2081 
2082  trantop(i,j) = c1i ! transmittance at top of ice
2083 
2084  enddo ! ij
2085 
2086  !-----------------------------------------------------------------
2087  ! penetrating SW absorbed in each ice layer
2088  !-----------------------------------------------------------------
2089 
2090  do k = 1, nilyr
2091 !DIR$ CONCURRENT !Cray
2092 !cdir nodep !NEC
2093 !ocl novrec !Fujitsu
2094  do ij = 1, icells
2095  i = indxi(ij)
2096  j = indxj(ij)
2097 
2098  tranbot(i,j) = exp(-kappav*hlyr(i,j)*real(k,kind=dbl_kind))
2099  iabs(i,j,k) = fswpen(i,j) * (trantop(i,j)-tranbot(i,j))
2100 
2101  ! bottom of layer k = top of layer k+1
2102  trantop(i,j) = tranbot(i,j)
2103 
2104  enddo ! ij
2105  enddo ! k
2106 
2107 !DIR$ CONCURRENT !Cray
2108 !cdir nodep !NEC
2109 !ocl novrec !Fujitsu
2110  do ij = 1, icells
2111  i = indxi(ij)
2112  j = indxj(ij)
2113 
2114  ! SW penetrating thru ice into ocean
2115  fswthrun(i,j) = fswpen(i,j) * tranbot(i,j)
2116 
2117  ! SW absorbed in ice interior
2118  fswint(i,j) = fswpen(i,j) - fswthrun(i,j)
2119 
2120  enddo ! ij
2121 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_albedo::alidrn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alidrn
Definition: ice_albedo.f90:76
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_flux::swvdf
real(kind=dbl_kind), dimension(:,:), allocatable, save swvdf
Definition: ice_flux.f90:91
ice_albedo::snowpatch
real(kind=dbl_kind), parameter snowpatch
Definition: ice_albedo.f90:59
ice_albedo::alidfn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alidfn
Definition: ice_albedo.f90:76
ice_albedo::alvdrn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alvdrn
Definition: ice_albedo.f90:76
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_flux::swvdr
real(kind=dbl_kind), dimension(:,:), allocatable, save swvdr
Definition: ice_flux.f90:91
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_flux::swidf
real(kind=dbl_kind), dimension(:,:), allocatable, save swidf
Definition: ice_flux.f90:91
ice_albedo
Definition: ice_albedo.f90:37
ice_flux::swidr
real(kind=dbl_kind), dimension(:,:), allocatable, save swidr
Definition: ice_flux.f90:91
ice_constants::kappav
real(kind=dbl_kind), parameter kappav
Definition: ice_constants.f90:118
ice_albedo::alvdfn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alvdfn
Definition: ice_albedo.f90:76
Here is the caller graph for this function:

◆ add_new_snow()

subroutine ice_therm_vertical::add_new_snow ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  hsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out)  hsn_new,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  qsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  Tsf 
)

Definition at line 2952 of file ice_therm_vertical.f90.

2952 !
2953 ! !USES:
2954 !
2955 ! !INPUT/OUTPUT PARAMETERS:
2956 !
2957  integer (kind=int_kind), intent(in) :: &
2958  ni & ! thickness category index
2959  , icells ! number of cells with aicen > puny
2960 
2961  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2962  intent(in) :: &
2963  indxi, indxj ! compressed indices for cells with aicen > puny
2964 
2965  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in):: &
2966  tsf ! ice/snow surface temperature, Tsfcn
2967 
2968  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2969  hsn & ! snow thickness (m)
2970  , qsn ! snow enthalpy
2971 
2972  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
2973  hsn_new ! thickness of new snow (m)
2974 !
2975 !EOP
2976 !
2977  real (kind=dbl_kind) :: &
2978  qsnew & ! enthalpy of new snow (J kg-1)
2979  , hstot ! snow thickness including new snow (J kg-1)
2980 
2981  integer (kind=int_kind) :: &
2982  i, j & ! horizontal indices
2983  , ij ! horizontal index, combines i and j loops
2984 
2985 
2986  hsn_new(:,:) = c0i
2987 
2988 !DIR$ CONCURRENT !Cray
2989 !cdir nodep !NEC
2990 !ocl novrec !Fujitsu
2991  do ij = 1, icells
2992  i = indxi(ij)
2993  j = indxj(ij)
2994 
2995  !----------------------------------------------------------------
2996  ! NOTE: If heat flux diagnostics are to work, new snow should
2997  ! have T = 0 (i.e. q = -rhos*Lfresh) and should not be
2998  ! converted to rain.
2999  !----------------------------------------------------------------
3000 
3001  if (fsnow(i,j) > c0i) then
3002 
3003 ! hsn_new(i,j) = fsnow(i,j)/rhos * dt
3004  hsn_new(i,j) = fsnow(i,j)/rhos * dtice
3005  qsnew = -rhos*lfresh
3006  hstot = hsn(i,j) + hsn_new(i,j)
3007  if (hstot > puny) then
3008  qsn(i,j) = (hsn(i,j) *qsn(i,j) &
3009  + hsn_new(i,j)*qsnew ) / hstot
3010  ! avoid roundoff errors if hsn=0
3011  qsn(i,j) = min(qsn(i,j), -rhos*lfresh)
3012  hsn(i,j) = hstot
3013  endif
3014  endif
3015 
3016  enddo ! ij
3017 
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_flux::fsnow
real(kind=dbl_kind), dimension(:,:), allocatable, save fsnow
Definition: ice_flux.f90:91
Here is the caller graph for this function:

◆ conductivity()

subroutine ice_therm_vertical::conductivity ( integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  hlyr,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  hsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(in)  Tin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  Tbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+1), intent(out)  khi,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out)  khs 
)

Definition at line 1887 of file ice_therm_vertical.f90.

1887 !
1888 ! !USES:
1889 !
1890 ! !INPUT/OUTPUT PARAMETERS:
1891 !
1892  integer (kind=int_kind), intent(in) :: &
1893  icells ! number of cells with aicen > puny
1894 
1895  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
1896  intent(in) :: &
1897  indxi, indxj ! compressed indices for cells with aicen > puny
1898 
1899  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
1900  hlyr & ! ice layer thickness
1901  , hsn & ! snow layer thickness
1902  , tbot ! ice bottom surface temperature (deg C)
1903 
1904  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
1905  intent(in) :: &
1906  tin ! internal ice layer temperatures
1907 
1908  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+1), &
1909  intent(out) :: &
1910  khi ! ki / hlyr
1911 
1912  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
1913  khs ! ksno / hsn
1914 !
1915 !EOP
1916 !
1917  real (kind=dbl_kind), parameter :: &
1918  betak = 0.13_dbl_kind &! constant in formula for k (W m-1 ppt-1)
1919  , kimin = 0.10_dbl_kind ! min conductivity of saline ice (W m-1 deg-1)
1920 
1921  integer (kind=int_kind) :: &
1922  i, j & ! horizontal indices
1923  , ij & ! horizontal index, combines i and j loops
1924  , k ! ice layer index
1925 
1926  real (kind=dbl_kind) :: &
1927  ki ! thermal cond of sea ice (W m-1 deg-1)
1928 
1929  khi(:,:,:) = c0i
1930  khs(:,:) = c0i
1931 
1932 !DIR$ CONCURRENT !Cray
1933 !cdir nodep !NEC
1934 !ocl novrec !Fujitsu
1935  do ij = 1, icells
1936  i = indxi(ij)
1937  j = indxj(ij)
1938 
1939  ! top surface of top ice layer
1940  ki = kice + betak*salin(1) / min(-puny,tin(i,j,1))
1941  ki = max(ki, kimin)
1942  khi(i,j,1) = ki / hlyr(i,j)
1943 
1944  ! bottom surface of bottom ice layer
1945  ki = kice + betak*salin(nilyr+1) / min(-puny,tbot(i,j))
1946  ki = max(ki, kimin)
1947  khi(i,j,nilyr+1) = ki / hlyr(i,j)
1948 
1949  ! if snow is present: top snow surface, snow-ice interface
1950  if (hsn(i,j) > hsnomin) then
1951  khs(i,j) = ksno / hsn(i,j)
1952  khi(i,j,1) = c2i*khi(i,j,1)*khs(i,j) / &
1953  (khi(i,j,1) + khs(i,j))
1954  endif
1955 
1956  enddo ! ij
1957 
1958  ! interior ice interfaces
1959  do k = 2, nilyr
1960 !DIR$ CONCURRENT !Cray
1961 !cdir nodep !NEC
1962 !ocl novrec !Fujitsu
1963  do ij = 1, icells
1964  i = indxi(ij)
1965  j = indxj(ij)
1966 
1967  ki = kice + betak*p5*(salin(k-1)+salin(k)) / &
1968  min(-puny, p5*(tin(i,j,k-1)+tin(i,j,k)))
1969  ki = max(ki, kimin)
1970  khi(i,j,k) = ki / hlyr(i,j)
1971 
1972  enddo ! ij
1973  enddo ! k
1974 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_constants::ksno
real(kind=dbl_kind), parameter ksno
Definition: ice_constants.f90:118
ice_constants::kice
real(kind=dbl_kind), parameter kice
Definition: ice_constants.f90:118
Here is the caller graph for this function:

◆ conservation_check_vthermo()

subroutine ice_therm_vertical::conservation_check_vthermo ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  fsurf,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  flatn,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  fhnetn,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  fswint,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  einit,
real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in)  efinal 
)

Definition at line 2837 of file ice_therm_vertical.f90.

2837 !
2838 ! !USES:
2839 !
2840 ! use ice_exit
2841 !
2842 ! !INPUT/OUTPUT PARAMETERS:
2843 !
2844  integer (kind=int_kind), intent(in) :: &
2845  ni & ! thickness category index (diagnostic only)
2846  , icells ! number of cells with aicen > puny
2847 
2848  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2849  intent(in) :: &
2850  indxi, indxj ! compressed indices for cells with aicen > puny
2851 
2852  real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in) :: &
2853  fsurf & ! net flux to top surface, not including fcondtop
2854  , flatn & ! surface downward latent heat (W m-2)
2855  , fhnetn & ! fbot, corrected for any surplus energy
2856  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
2857  , einit & ! initial energy of melting (J m-2)
2858  , efinal ! final energy of melting (J m-2)
2859 !
2860 !EOP
2861 !
2862  integer (kind=int_kind) :: &
2863  i, j & ! horizontal indices
2864  , ij ! horizontal index, combines i and j loops
2865 
2866  real (kind=dbl_kind) :: &
2867  einp & ! energy input during timestep (J m-2)
2868  , ferr ! energy conservation error (W m-2)
2869 
2870  logical (kind=log_kind) :: & ! for vector-friendly error checks
2871  ferr_flag ! flag for energy error, ferr > ferrmax
2872 
2873  ferr_flag = .false. ! ferr <= ferrmax, initialization
2874 
2875 !DIR$ CONCURRENT !Cray
2876 !cdir nodep !NEC
2877 !ocl novrec !Fujitsu
2878  do ij = 1, icells
2879  i = indxi(ij)
2880  j = indxj(ij)
2881 
2882  einp = (fsurf(i,j) - flatn(i,j) + fswint(i,j) - fhnetn(i,j)) &
2883  * dtice
2884 ! * dt
2885 ! ferr = abs(efinal(i,j)-einit(i,j)-einp) / dt
2886  ferr = abs(efinal(i,j)-einit(i,j)-einp) / dtice
2887  if (ferr > ferrmax) then
2888  ferr_flag = .true.
2889  endif
2890  enddo
2891  !----------------------------------------------------------------
2892  ! If energy is not conserved, print diagnostics and exit.
2893  !----------------------------------------------------------------
2894  if (ferr_flag .and. dbg_set(dbg_sbrio)) then
2895  do ij = 1, icells
2896  i = indxi(ij)
2897  j = indxj(ij)
2898 
2899  !-----------------------------------------------------------------
2900  ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent
2901  ! heat is not included in the energy input, since de is the energy
2902  ! change in the system ice + vapor, and the latent heat lost by
2903  ! the ice is equal to that gained by the vapor.
2904  !-----------------------------------------------------------------
2905 
2906  einp = (fsurf(i,j) - flatn(i,j) + fswint(i,j) - &
2907  fhnetn(i,j)) * dtice
2908 ! fhnetn(i,j)) * dt
2909 ! ferr = abs(efinal(i,j)-einit(i,j)-einp) / dt
2910  ferr = abs(efinal(i,j)-einit(i,j)-einp) / dtice
2911  if (ferr > ferrmax) then
2912 ! write(nu_diag,*) 'Energy error, cat',ni, &
2913 ! 'my_task,i,j',my_task,i,j
2914 ! write(nu_diag,*)'wind(I,j)',wind(i,j),shcoef(i,j), potT(i,j)
2915 ! write(nu_diag,*) 'Flux error (W/m^2) =', ferr
2916 ! write(nu_diag,*) 'Energy error (J) =', ferr*dt
2917 ! write(nu_diag,*) 'Energy error (J) =', ferr*dtice
2918 ! write(nu_diag,*) 'Initial energy =', einit(i,j)
2919 ! write(nu_diag,*) 'Final energy =', efinal(i,j)
2920 ! write(nu_diag,*) 'efinal - einit =', &
2921 ! efinal(i,j)-einit(i,j)
2922 ! write(nu_diag,*) 'Input energy =', einp
2923 ! stoplabel = 'ice state at stop'
2924 ! call print_state(stoplabel,i,j)
2925 ! call abort_ice &
2926 ! ('vertical thermo: energy conservation error')
2927  endif
2928  enddo
2929  endif
2930 
mod_time::abs
Definition: mod_time.f90:64
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
Here is the call graph for this function:
Here is the caller graph for this function:

◆ frzmlt_bottom_lateral()

subroutine ice_therm_vertical::frzmlt_bottom_lateral ( real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  Tbot,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fbot,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  rside 
)

Definition at line 477 of file ice_therm_vertical.f90.

477 !
478 ! !USES:
479 !
480 ! !INPUT/OUTPUT PARAMETERS:
481 !
482  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
483  tbot & ! ice bottom surface temperature (deg C)
484  , fbot & ! heat flux to ice bottom (W/m^2)
485  , rside ! fraction of ice that melts laterally
486 !
487 !EOP
488 !
489  integer (kind=int_kind) :: &
490  i, j & ! horizontal indices
491  , ni & ! thickness category index
492  , k & ! layer index
493  , ij & ! horizontal index, combines i and j loops
494  , icells ! number of cells with ice melting
495 
496  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
497  indxi, indxj ! compressed indices for cells with ice melting
498 
499  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
500  etot & ! total energy in column
501  , fside ! lateral heat flux (W/m^2)
502 
503  real (kind=dbl_kind) :: &
504  deltat &! SST - Tbot >= 0
505  , ustar &! skin friction velocity for fbot (m/s)
506  , wlat &! lateral melt rate (m/s)
507  , xtmp ! temporarty variable
508 
509  ! Parameters for bottom melting
510 
511  ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut
512 
513  real (kind=dbl_kind), parameter :: &
514  cpchr = -cp_ocn*rhow*0.006_dbl_kind &
515  , ustar_min = 5.e-3_dbl_kind
516 
517  ! Parameters for lateral melting
518 
519  real (kind=dbl_kind), parameter :: &
520  floediam = 300.0_dbl_kind & ! effective floe diameter (m)
521  , alpha = 0.66_dbl_kind & ! constant from Steele (unitless)
522  , m1 = 1.6e-6_dbl_kind & ! constant from Maykut & Perovich
523  ! (m/s/deg^(-m2))
524  , m2 = 1.36_dbl_kind ! constant from Maykut & Perovich
525  ! (unitless)
526 
527  !-----------------------------------------------------------------
528  ! Set bottom ice surface temperature and initialize fluxes.
529  !-----------------------------------------------------------------
530 
531 
532  do j = jlo, jhi
533  do i = ilo, ihi
534  tbot(i,j) = tf(i,j)
535  fbot(i,j) = c0i
536  rside(i,j) = c0i
537  fside(i,j) = c0i
538  etot(i,j) = c0i
539  enddo ! i
540  enddo ! j
541 
542  !-----------------------------------------------------------------
543  ! Identify grid cells where ice can melt.
544  !-----------------------------------------------------------------
545 
546  icells = 0
547  do j = jlo, jhi
548  do i = ilo, ihi
549  if (aice(i,j) > puny .and. frzmlt(i,j) < c0i) then ! ice can melt
550  icells = icells + 1
551  indxi(icells) = i
552  indxj(icells) = j
553  endif
554  enddo ! i
555  enddo ! j
556 
557  do ij = 1, icells ! cells where ice can melt
558  i = indxi(ij)
559  j = indxj(ij)
560 
561  !-----------------------------------------------------------------
562  ! Use boundary layer theory for fbot.
563  ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703.
564  !-----------------------------------------------------------------
565 
566  deltat = max((sst(i,j)-tbot(i,j)),c0i)
567 
568  ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
569  ustar = sqrt(sqrt(strocnxt(i,j)**2+strocnyt(i,j)**2)/rhow)
570  ustar = max(ustar,ustar_min*ustar_scale)
571 
572  fbot(i,j) = cpchr * deltat * ustar ! < 0
573 
574  fbot(i,j) = max(fbot(i,j), frzmlt(i,j)) ! frzmlt < fbot < 0
575 
576 !!! uncomment to use all frzmlt for standalone runs
577 ! fbot(i,j) = min (c0i, frzmlt(i,j))
578 
579  !-----------------------------------------------------------------
580  ! Compute rside. See these references:
581  ! Maykut and Perovich (1987): JGR, 92, 7032-7044
582  ! Steele (1992): JGR, 97, 17,729-17,738
583  !-----------------------------------------------------------------
584 
585  wlat = m1 * deltat**m2 ! Maykut & Perovich
586 ! rside(i,j) = wlat*dt*pi/(alpha*floediam) ! Steele
587  rside(i,j) = wlat*dtice*pi/(alpha*floediam) ! Steele
588 
589 
590  rside(i,j) = max(c0i,min(rside(i,j),c1i))
591 
592  enddo ! ij
593 
594  !-----------------------------------------------------------------
595  ! Compute heat flux associated with this value of rside.
596  !-----------------------------------------------------------------
597 
598  do ni = 1, ncat
599 
600  ! melting energy/unit area in each column, etot < 0
601  do ij = 1, icells
602  i = indxi(ij)
603  j = indxj(ij)
604  etot(i,j) = esnon(i,j,ni)
605  enddo ! ij
606 
607  do k = 1, nilyr
608 !DIR$ CONCURRENT !Cray
609 !cdir nodep !NEC
610 !ocl novrec !Fujitsu
611  do ij = 1, icells
612  i = indxi(ij)
613  j = indxj(ij)
614  etot(i,j) = etot(i,j) + eicen(i,j,ilyr1(ni)+k-1)
615  enddo ! ij
616  enddo ! k
617 
618 !DIR$ CONCURRENT !Cray
619 !cdir nodep !NEC
620 !ocl novrec !Fujitsu
621  do ij = 1, icells
622  i = indxi(ij)
623  j = indxj(ij)
624  ! lateral heat flux
625 ! fside(i,j) = fside(i,j) + rside(i,j)*etot(i,j)/dt ! fside < 0
626  fside(i,j) = fside(i,j) + rside(i,j)*etot(i,j)/dtice ! fside < 0
627  enddo ! ij
628 
629  enddo ! n
630 
631  !-----------------------------------------------------------------
632  ! Limit bottom and lateral heat fluxes if necessary.
633  !-----------------------------------------------------------------
634 
635 !DIR$ CONCURRENT !Cray
636 !cdir nodep !NEC
637 !ocl novrec !Fujitsu
638 
639 
640  do ij = 1, icells
641  i = indxi(ij)
642  j = indxj(ij)
643 
644  xtmp = frzmlt(i,j)/(fbot(i,j) + fside(i,j) + puny)
645  xtmp = min(xtmp, c1i)
646  fbot(i,j) = fbot(i,j) * xtmp
647  rside(i,j) = rside(i,j) * xtmp
648  enddo ! ij
649 
all_vars::alpha
real(sp), dimension(:), allocatable, target alpha
Definition: mod_main.f90:1336
ice_flux::tf
real(kind=dbl_kind), dimension(:,:), allocatable, save tf
Definition: ice_flux.f90:91
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_state::esnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save esnon
Definition: ice_state.f90:97
control::pi
real(dp), parameter pi
Definition: mod_main.f90:880
ice_flux::strocnxt
real(kind=dbl_kind), dimension(:,:), allocatable, save strocnxt
Definition: ice_flux.f90:55
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_flux::sst
real(kind=dbl_kind), dimension(:,:), allocatable, save sst
Definition: ice_flux.f90:91
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_constants::rhow
real(kind=dbl_kind), parameter rhow
Definition: ice_constants.f90:85
ice_state::aice
real(kind=dbl_kind), dimension(:,:), allocatable, target, save aice
Definition: ice_state.f90:82
ice_state::eicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save eicen
Definition: ice_state.f90:108
ice_model_size::ncat
integer(kind=int_kind), parameter ncat
Definition: ice_model_size.f90:48
ice_flux::frzmlt
real(kind=dbl_kind), dimension(:,:), allocatable, save frzmlt
Definition: ice_flux.f90:91
ice_flux::strocnyt
real(kind=dbl_kind), dimension(:,:), allocatable, save strocnyt
Definition: ice_flux.f90:55
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_therm_vertical::rside
real(kind=dbl_kind), dimension(:,:), allocatable, save rside
Definition: ice_therm_vertical.f90:88
ice_constants::cp_ocn
real(kind=dbl_kind), parameter cp_ocn
Definition: ice_constants.f90:85
Here is the caller graph for this function:

◆ init_thermo_vars()

subroutine ice_therm_vertical::init_thermo_vars ( real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  strxn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  stryn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  Trefn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  Qrefn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fsensn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  flatn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fswabsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  flwoutn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  evapn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  freshn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fsaltn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fhnetn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fswthrun,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fsurf,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fcondtop,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fcondbot,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  fswint,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  einit,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  efinal,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  mvap 
)

Definition at line 675 of file ice_therm_vertical.f90.

675 !
676 ! !USES:
677 !
678 ! !INPUT/OUTPUT PARAMETERS:
679 !
680  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
681 
682  ! fields aggregated for coupler
683  strxn & ! air/ice zonal strss (N/m^2)
684  , stryn & ! air/ice merdnl strss (N/m^2)
685  , trefn & ! air tmp rfrnc level (K)
686  , qrefn & ! air sp hum rfrnc level (kg/kg)
687  , fsensn & ! surface downward sensible heat (W m-2)
688  , flatn & ! surface downward latent heat (W m-2)
689  , fswabsn & ! SW absorbed by ice (W m-2)
690  , flwoutn & ! upward LW at surface (W m-2)
691  , evapn & ! flux of vapor, atmos to ice (kg m-2 s-1)
692  , freshn & ! flux of water, ice to ocean (kg m-2 s-1)
693  , fsaltn & ! flux of salt, ice to ocean (kg m-2 s-1)
694  , fhnetn & ! fbot, corrected for surplus energy (W m-2)
695  , fswthrun & ! SW through ice to ocean (W m-2)
696 
697  ! other fields passed among subroutines
698  , fsurf & ! net flux to top surface, not including fcondtop
699  , fcondtop & ! downward cond flux at top surface (W m-2)
700  , fcondbot & ! downward cond flux at bottom surface (W m-2)
701  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
702  , einit & ! initial energy of melting (J m-2)
703  , efinal & ! final energy of melting (J m-2)
704  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
705 !
706 !EOP
707 !
708  integer (kind=int_kind) :: &
709  i, j ! horizontal indices
710 
711  do j = jlo, jhi
712  do i = ilo, ihi
713 
714  strxn(i,j) = c0i
715  stryn(i,j) = c0i
716  trefn(i,j) = c0i
717  qrefn(i,j) = c0i
718  fsensn(i,j) = c0i
719  flatn(i,j) = c0i
720  fswabsn(i,j) = c0i
721  flwoutn(i,j) = c0i
722  evapn(i,j) = c0i
723  freshn(i,j) = c0i
724  fsaltn(i,j) = c0i
725  fhnetn(i,j) = c0i
726  fswthrun(i,j) = c0i
727 
728  fsurf(i,j) = c0i
729  fcondtop(i,j) = c0i
730  fcondbot(i,j) = c0i
731  fswint(i,j) = c0i
732  einit(i,j) = c0i
733  efinal(i,j) = c0i
734  mvap(i,j) = c0i
735 
736  enddo ! i
737  enddo ! j
738 
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
Here is the caller graph for this function:

◆ init_thermo_vertical()

subroutine ice_therm_vertical::init_thermo_vertical ( )

Definition at line 117 of file ice_therm_vertical.f90.

117 !
118 ! !USES:
119 !
120  use ice_itd
121 !
122 ! !INPUT/OUTPUT PARAMETERS:
123 !
124 !EOP
125 !
126  real (kind=dbl_kind), parameter ::&
127  nsal = 0.407_dbl_kind&
128  , msal = 0.573_dbl_kind&
129  , saltmax = 3.2_dbl_kind &! max salinity at ice base (ppt)
130  , min_salin = 0.1_dbl_kind ! threshold for brine pocket treatment
131 
132  integer (kind=int_kind) :: k ! ice layer index
133  real (kind=dbl_kind) :: zn ! normalized ice thickness
134 
135  !-----------------------------------------------------------------
136  ! Determine l_brine based on saltmax.
137  ! Thermodynamic solver will not converge if l_brine is true and
138  ! saltmax is close to zero.
139  !-----------------------------------------------------------------
140 
141  if (saltmax > min_salin) then
142  l_brine = .true.
143  else
144  l_brine = .false.
145  endif
146 
147  ! salinity and melting temperature profile
148  if (l_brine) then
149  do k = 1, nilyr
150  zn = (real(k,kind=dbl_kind)-p5) / real(nilyr,kind=dbl_kind)
151  salin(k)=(saltmax/c2i)*(c1i-cos(pi*zn**(nsal/(msal+zn))))
152 ! salin(k)=saltmax ! for isosaline ice
153  enddo
154  salin(nilyr+1) = saltmax
155  do k = 1, nilyr+1
156  tmlt(k) = -salin(k)*depresst
157  enddo
158  else
159  do k = 1, nilyr+1
160  salin(k) = c0i
161  tmlt(k) = c0i
162  enddo
163  endif
164  ustar_scale = c1i ! for nonzero currents
165 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_constants::depresst
real(kind=dbl_kind), parameter depresst
Definition: ice_constants.f90:118
control::pi
real(dp), parameter pi
Definition: mod_main.f90:880
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_itd
Definition: ice_itd.f90:45

◆ init_vertical_profile()

subroutine ice_therm_vertical::init_vertical_profile ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hin,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hlyr,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hin_init,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hsn_init,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(out)  qin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(out)  Tin,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  qsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  Tsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  Tsf,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout)  einit 
)

Definition at line 763 of file ice_therm_vertical.f90.

763 !
764 ! !USES:
765 !
766 ! use ice_exit
767 !
768 ! !INPUT/OUTPUT PARAMETERS:
769 !
770  integer (kind=int_kind), intent(in) :: &
771  ni & ! thickness category index
772  , icells ! number of cells with aicen > puny
773 
774  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
775  intent(in) :: &
776  indxi, indxj ! compressed indices for cells with aicen > puny
777 
778  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
779  hin &! ice thickness (m)
780  , hsn &! snow thickness (m)
781  , hlyr &! ice layer thickness
782  , hin_init &! initial value of hin
783  , hsn_init &! initial value of hsn
784  , qsn &! snow enthalpy
785  , tsn &! snow temperature
786  , tsf ! ice/snow surface temperature, Tsfcn
787 
788  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
789  intent(out) :: &
790  qin &! ice layer enthalpy (J m-3)
791  , tin ! internal ice layer temperatures
792 
793  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout) :: &
794  einit ! initial energy of melting (J m-2)
795 !
796 !EOP
797 !
798  real (kind=dbl_kind), parameter :: &
799  tmin = -100._dbl_kind ! min allowed internal temperature (deg C)
800 
801  integer (kind=int_kind) :: &
802  i, j & ! horizontal indices
803  , ij & ! horizontal index, combines i and j loops
804  , k ! ice layer index
805 
806  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
807  tmax ! maximum allowed snow temperature (deg C)
808 
809  real (kind=dbl_kind) :: &
810  aa1, bb1, cc1 ! terms in quadratic formula
811 
812  logical (kind=log_kind) :: & ! for vector-friendly error checks
813  tsno_high & ! flag for Tsn > Tmax
814  , tice_high & ! flag for Tin > Tmlt
815  , tsno_low & ! flag for Tsn < Tmin
816  , tice_low ! flag for Tin < Tmin
817 
818  !-----------------------------------------------------------------
819  ! Initialize arrays
820  !-----------------------------------------------------------------
821 
822  do j = jlo, jhi
823  do i = ilo, ihi
824  hin(i,j) = c0i
825  hsn(i,j) = c0i
826  hlyr(i,j) = c0i
827  hin_init(i,j) = c0i
828  hsn_init(i,j) = c0i
829  qsn(i,j) = c0i
830  tsn(i,j) = c0i
831  tsf(i,j) = c0i
832  tmax(i,j) = c0i
833  enddo
834  enddo
835 
836  do k = 1, nilyr
837  do j = jlo, jhi
838  do i = ilo, ihi
839  qin(i,j,k) = c0i
840  tin(i,j,k) = c0i
841  enddo
842  enddo
843  enddo
844 
845  !-----------------------------------------------------------------
846  ! Initialize flags
847  !-----------------------------------------------------------------
848 
849  tsno_high = .false.
850  tice_high = .false.
851  tsno_low = .false.
852  tice_low = .false.
853 
854  !-----------------------------------------------------------------
855  ! Load arrays for vertical thermo calculation.
856  !-----------------------------------------------------------------
857 !DIR$ CONCURRENT !Cray
858 !cdir nodep !NEC
859 !ocl novrec !Fujitsu
860  do ij = 1, icells
861  i = indxi(ij)
862  j = indxj(ij)
863 
864  !-----------------------------------------------------------------
865  ! Surface temperature, ice and snow thickness, snow enthalpy
866  !
867  ! Tmax based on the idea that dT ~ dq / (rhos*cp_ice)
868  ! dq ~ q dv / v
869  ! dv ~ eps11
870  ! where 'd' denotes an error due to roundoff.
871  !-----------------------------------------------------------------
872 
873  tsf(i,j) = tsfcn(i,j,ni)
874  hin(i,j) = vicen(i,j,ni) / aicen(i,j,ni)
875  hlyr(i,j) = hin(i,j) / real(nilyr,kind=dbl_kind)
876  hsn(i,j) = vsnon(i,j,ni) / aicen(i,j,ni)
877  hin_init(i,j) = hin(i,j)
878  hsn_init(i,j) = hsn(i,j)
879 
880  if (hsn(i,j) > hsnomin) then
881  qsn(i,j) = esnon(i,j,ni) / vsnon(i,j,ni) ! qsn, esnon < 0
882  tmax(i,j) = -qsn(i,j)*eps11 / (rhos*cp_ice*vsnon(i,j,ni))
883  else
884  qsn(i,j) = -rhos * lfresh
885  tmax(i,j) = puny
886  endif
887 
888  !-----------------------------------------------------------------
889  ! Compute snow temperatures from enthalpies.
890  !-----------------------------------------------------------------
891  tsn(i,j) = (lfresh + qsn(i,j)/rhos)/cp_ice ! qsn <= -rhos*Lfresh
892 
893  !-----------------------------------------------------------------
894  ! Check for Tsn > Tmax (allowing for roundoff error)
895  ! and Tsn < Tmin.
896  !-----------------------------------------------------------------
897  if (tsn(i,j) > tmax(i,j)) then
898  tsno_high = .true.
899  elseif (tsn(i,j) < tmin) then
900  tsno_low = .true.
901  endif
902 
903  enddo
904 
905  !-----------------------------------------------------------------
906  ! If Tsn is out of bounds, print diagnostics and exit.
907  !-----------------------------------------------------------------
908  if (tsno_high) then
909  do ij = 1, icells
910  i = indxi(ij)
911  j = indxj(ij)
912 
913  if (tsn(i,j) > tmax(i,j)) then ! allowing for roundoff error
914 ! write(nu_diag,*) ' '
915 ! write(nu_diag,*) 'Starting thermo, Tsn > Tmax, cat',ni
916 ! write(nu_diag,*) 'Tsn=',Tsn(i,j)
917 ! write(nu_diag,*) 'Tmax=',Tmax(i,j)
918 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
919 ! istep1, my_task, i, j
920 ! write(nu_diag,*) 'qsn',qsn(i,j)
921 ! stoplabel = 'ice state at stop'
922 ! call print_state(stoplabel,i,j)
923 ! call abort_ice('vertical thermo: Tsn must be <= 0')
924  endif
925 
926  enddo ! ij
927  endif ! tsno_high
928 
929  if (tsno_low) then
930  do ij = 1, icells
931  i = indxi(ij)
932  j = indxj(ij)
933 
934 ! if (Tsn(i,j) < Tmin) then ! allowing for roundoff error
935 ! write(nu_diag,*) ' '
936 ! write(nu_diag,*) 'Starting thermo, Tsn < Tmin, cat',ni
937 ! write(nu_diag,*) 'Tsn=', Tsn(i,j)
938 ! write(nu_diag,*) 'Tmin=', Tmin
939 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
940 ! istep1, my_task, i, j
941 ! write(nu_diag,*) 'qsn', qsn(i,j)
942 ! stoplabel = 'ice state at stop'
943 ! call print_state(stoplabel,i,j)
944 ! call abort_ice('vertical thermo: Tsn must be <= 0')
945 ! endif
946 
947  enddo ! ij
948  endif ! tsno_low
949 
950  !-----------------------------------------------------------------
951  ! initial energy per unit area of ice/snow, relative to 0 C
952  ! incremented later for ice
953  !-----------------------------------------------------------------
954 
955 !DIR$ CONCURRENT !Cray
956 !cdir nodep !NEC
957 !ocl novrec !Fujitsu
958  do ij = 1, icells
959  i = indxi(ij)
960  j = indxj(ij)
961 
962  if (tsn(i,j) > c0i) then ! correct roundoff error
963  tsn(i,j) = c0i
964  qsn(i,j) = -rhos*lfresh
965  endif
966 
967  einit(i,j) = hsn(i,j)*qsn(i,j)
968 
969  enddo ! ij
970 
971  do k = 1, nilyr
972 !DIR$ CONCURRENT !Cray
973 !cdir nodep !NEC
974 !ocl novrec !Fujitsu
975  do ij = 1, icells
976  i = indxi(ij)
977  j = indxj(ij)
978 
979  !-----------------------------------------------------------------
980  ! Compute ice enthalpy
981  !-----------------------------------------------------------------
982 
983  qin(i,j,k) = eicen(i,j,ilyr1(ni)+k-1) & ! qin < 0
984  * real(nilyr,kind=dbl_kind)/vicen(i,j,ni)
985 
986  !-----------------------------------------------------------------
987  ! Compute ice temperatures from enthalpies using quadratic formula
988  !-----------------------------------------------------------------
989 
990  if (l_brine) then
991  aa1 = cp_ice
992  bb1 = (cp_ocn-cp_ice)*tmlt(k) - qin(i,j,k)/rhoi - lfresh
993  cc1 = lfresh * tmlt(k)
994  tin(i,j,k) = (-bb1 - sqrt(bb1*bb1 - c4i*aa1*cc1)) / &
995  (c2i*aa1)
996  tmax(i,j) = tmlt(k)
997 
998  else ! fresh ice
999  tin(i,j,k) = (lfresh + qin(i,j,k)/rhoi) / cp_ice
1000  tmax(i,j) = -qin(i,j,k)*eps11/(rhos*cp_ice*vicen(i,j,ni))
1001  ! as above for snow
1002  endif
1003 
1004  IF (tin(i,j,k) < tmin)THEN
1005  !(Tin=Tmlt(k))
1006  eicen(i,j,ilyr1(ni)+k-1) = &
1007  (rhoi *cp_ocn*tmlt(k)) &
1008  * vicen(i,j,ni)/real(nilyr,kind=dbl_kind)
1009  qin(i,j,k) = eicen(i,j,ilyr1(ni)+k-1) & ! qin < 0
1010  * real(nilyr,kind=dbl_kind)/vicen(i,j,ni)
1011  aa1 = cp_ice
1012  bb1 = (cp_ocn-cp_ice)*tmlt(k) - qin(i,j,k)/rhoi - lfresh
1013  cc1 = lfresh * tmlt(k)
1014  tin(i,j,k) = (-bb1 - sqrt(bb1*bb1 - c4i*aa1*cc1)) / &
1015  (c2i*aa1)
1016  end if
1017 
1018 
1019 
1020  !-----------------------------------------------------------------
1021  ! Check for Tin > Tmax and Tin < Tmin
1022  !-----------------------------------------------------------------
1023 
1024  if (tin(i,j,k) > tmax(i,j)) then
1025  tice_high = .true.
1026  elseif (tin(i,j,k) < tmin) then
1027  tice_low = .true.
1028  endif
1029 
1030  enddo ! ij
1031 
1032  !-----------------------------------------------------------------
1033  ! If Tin is out of bounds, print diagnostics and exit.
1034  !-----------------------------------------------------------------
1035 
1036  if (tice_high) then
1037  do ij = 1, icells
1038  i = indxi(ij)
1039  j = indxj(ij)
1040 
1041  if (tin(i,j,k) > tmax(i,j)) then
1042 ! write(nu_diag,*) ' '
1043 ! write(nu_diag,*) 'Starting thermo, T > Tmax, cat',&
1044 ! ni,', layer',k
1045 ! write(nu_diag,*) 'Tin=',Tin(i,j,k),', Tmax=',Tmax(i,j)
1046 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
1047 ! istep1, my_task, i, j
1048 ! write(nu_diag,*) 'qin',qin(i,j,k)
1049 ! stoplabel = 'ice state at stop'
1050 ! call print_state(stoplabel,i,j)
1051 ! call abort_ice('vertical thermo: Tin must be <= Tmax')
1052  endif
1053 
1054  enddo ! ij
1055  endif ! tice_high
1056 
1057  if (tice_low) then
1058  do ij = 1, icells
1059  i = indxi(ij)
1060  j = indxj(ij)
1061 
1062 ! if (Tin(i,j,k) < Tmin) then
1063 ! write(nu_diag,*) ' '
1064 ! write(nu_diag,*) 'Starting thermo T < Tmin, cat', &
1065 ! ni,', layer',k
1066 ! write(nu_diag,*) 'Tin =', Tin(i,j,k)
1067 ! write(nu_diag,*) 'Tmin =', Tmin
1068 !! write(nu_diag,*) 'istep1, my_task, i, j:', &
1069 ! istep1, my_task, i, j
1070 ! stoplabel = 'ice state at stop'
1071 ! call print_state(stoplabel,i,j)
1072 ! call abort_ice('vertical_thermo: Tin must be > Tmin')
1073 ! endif
1074  enddo ! ij
1075  endif ! tice_low
1076 
1077  !-----------------------------------------------------------------
1078  ! initial energy per unit area of ice/snow, relative to 0 C
1079  !-----------------------------------------------------------------
1080 
1081 !DIR$ CONCURRENT !Cray
1082 !cdir nodep !NEC
1083 !ocl novrec !Fujitsu
1084  do ij = 1, icells
1085  i = indxi(ij)
1086  j = indxj(ij)
1087 
1088  if (tin(i,j,k) > c0i) then ! correct roundoff error
1089  tin(i,j,k) = c0i
1090  qin(i,j,k) = -rhoi*lfresh
1091  endif
1092 
1093  einit(i,j) = einit(i,j) + hlyr(i,j)*qin(i,j,k)
1094 
1095  enddo ! ij
1096  enddo ! k
1097 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_state::esnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save esnon
Definition: ice_state.f90:97
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_constants::c4i
real(kind=dbl_kind), parameter c4i
Definition: ice_constants.f90:150
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_constants::cp_ice
real(kind=dbl_kind), parameter cp_ice
Definition: ice_constants.f90:85
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_state::tsfcn
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save tsfcn
Definition: ice_state.f90:97
ice_state::vicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vicen
Definition: ice_state.f90:97
ice_state::eicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save eicen
Definition: ice_state.f90:108
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
ice_constants::eps11
real(kind=dbl_kind), parameter eps11
Definition: ice_constants.f90:150
ice_itd::ilyr1
integer(kind=int_kind), dimension(ncat) ilyr1
Definition: ice_itd.f90:62
ice_state::vsnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vsnon
Definition: ice_state.f90:97
ice_constants::cp_ocn
real(kind=dbl_kind), parameter cp_ocn
Definition: ice_constants.f90:85
Here is the caller graph for this function:

◆ surface_fluxes()

subroutine ice_therm_vertical::surface_fluxes ( integer (kind=int_kind), intent(in)  isolve,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxii,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxjj,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  Tsf,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  fswsfc,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  flwoutn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fsensn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  flatn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fsurf,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  dflwout_dT,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  dfsens_dT,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  dflat_dT,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  dfsurf_dT 
)

Definition at line 2146 of file ice_therm_vertical.f90.

2146 !
2147 ! !USES:
2148 !
2149 ! !INPUT/OUTPUT PARAMETERS:
2150 !
2151  integer (kind=int_kind), intent(in) :: &
2152  isolve ! number of cells with temps not converged
2153 
2154  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2155  intent(in) :: &
2156  indxii, indxjj ! compressed indices for cells not converged
2157 
2158  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
2159  tsf & ! ice/snow surface temperature, Tsfcn
2160  , fswsfc ! SW absorbed at ice/snow surface (W m-2)
2161 
2162  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2163  fsensn & ! surface downward sensible heat (W m-2)
2164  , flatn & ! surface downward latent heat (W m-2)
2165  , flwoutn & ! upward LW at surface (W m-2)
2166  , fsurf & ! net flux to top surface, not including fcondtop
2167  , dfsens_dt & ! deriv of fsens wrt Tsf (W m-2 deg-1)
2168  , dflat_dt & ! deriv of flat wrt Tsf (W m-2 deg-1)
2169  , dflwout_dt & ! deriv of flwout wrt Tsf (W m-2 deg-1)
2170  , dfsurf_dt ! derivative of fsurf wrt Tsf
2171 !
2172 !EOP
2173 !
2174  integer (kind=int_kind) :: &
2175  i, j & ! horizontal indices
2176  , ij & ! horizontal index, combines i and j loops
2177  , k ! ice layer index
2178 
2179  real (kind=dbl_kind) :: &
2180  tsfk & ! ice/snow surface temperature (K)
2181  , qsfc & ! saturated surface specific humidity (kg/kg)
2182  , dqsfcdt & ! derivative of Qsfc wrt surface temperature
2183  , qsat & ! the saturation humidity of air (kg/m^3)
2184  , flwdabs & ! downward longwave absorbed heat flx (W/m^2)
2185  , tmpvar ! 1/TsfK
2186 
2187 !DIR$ CONCURRENT !Cray
2188 !cdir nodep !NEC
2189 !ocl novrec !Fujitsu
2190  do ij = 1, isolve
2191  i = indxii(ij) ! NOTE: not indxi and indxj
2192  j = indxjj(ij)
2193 
2194 
2195  ! ice surface temperature in Kelvin
2196  tsfk = tsf(i,j) + tffresh
2197  tmpvar = c1i/tsfk
2198 
2199  ! saturation humidity
2200  qsat = qqqice * exp(-tttice*tmpvar)
2201  qsfc = qsat / rhoa(i,j)
2202  dqsfcdt = tttice * tmpvar*tmpvar * qsfc
2203 
2204  ! longwave radiative flux
2205  flwdabs = emissivity * flw(i,j)
2206  flwoutn(i,j) = -emissivity * stefan_boltzmann * tsfk**4
2207 
2208  ! downward latent and sensible heat fluxes
2209  fsensn(i,j) = shcoef(i,j) * (pott(i,j) - tsfk)
2210  flatn(i,j) = lhcoef(i,j) * (qa(i,j) - qsfc)
2211 
2212  ! derivatives wrt surface temp
2213  dflwout_dt(i,j) = - emissivity*stefan_boltzmann * c4i*tsfk**3
2214  dfsens_dt(i,j) = - shcoef(i,j)
2215  dflat_dt(i,j) = - lhcoef(i,j) * dqsfcdt
2216 
2217  fsurf(i,j) = fswsfc(i,j) + flwdabs + flwoutn(i,j) &
2218  + fsensn(i,j) + flatn(i,j)
2219  dfsurf_dt(i,j) = dflwout_dt(i,j) &
2220  + dfsens_dt(i,j) + dflat_dt(i,j)
2221  enddo
2222 
ice_flux::lhcoef
real(kind=dbl_kind), dimension(:,:), allocatable, save lhcoef
Definition: ice_flux.f90:164
ice_constants::tttice
real(kind=dbl_kind), parameter tttice
Definition: ice_constants.f90:134
ice_constants::qqqice
real(kind=dbl_kind), parameter qqqice
Definition: ice_constants.f90:134
ice_flux::rhoa
real(kind=dbl_kind), dimension(:,:), allocatable, save rhoa
Definition: ice_flux.f90:91
ice_flux::pott
real(kind=dbl_kind), dimension(:,:), allocatable, save pott
Definition: ice_flux.f90:91
ice_constants::stefan_boltzmann
real(kind=dbl_kind), parameter stefan_boltzmann
Definition: ice_constants.f90:85
ice_constants::c4i
real(kind=dbl_kind), parameter c4i
Definition: ice_constants.f90:150
ice_flux::flw
real(kind=dbl_kind), dimension(:,:), allocatable, save flw
Definition: ice_flux.f90:91
ice_constants::emissivity
real(kind=dbl_kind), parameter emissivity
Definition: ice_constants.f90:118
ice_constants::tffresh
real(kind=dbl_kind), parameter tffresh
Definition: ice_constants.f90:85
ice_flux::shcoef
real(kind=dbl_kind), dimension(:,:), allocatable, save shcoef
Definition: ice_flux.f90:164
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_flux::qa
real(kind=dbl_kind), dimension(:,:), allocatable, save qa
Definition: ice_flux.f90:91
Here is the caller graph for this function:

◆ temperature_changes()

subroutine ice_therm_vertical::temperature_changes ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  hlyr,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  hsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(inout)  qin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(inout)  Tin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  qsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  Tsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  Tsf,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  Tbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fsensn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  flatn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fswabsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  flwoutn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fswthrun,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fhnetn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fsurf,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fcondtop,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fcondbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fswint,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  einit 
)

Definition at line 1128 of file ice_therm_vertical.f90.

1128 !
1129 ! !USES:
1130 !
1131 ! use ice_exit
1132 !
1133 ! !INPUT/OUTPUT PARAMETERS:
1134 !
1135  integer (kind=int_kind), intent(in) :: &
1136  ni & ! thickness category index
1137  , icells ! number of cells with aicen > puny
1138 
1139  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
1140  intent(in) :: &
1141  indxi, indxj ! compressed indices for cells with aicen > puny
1142 
1143  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
1144  hlyr & ! ice layer thickness
1145  , hsn & ! snow thickness (m)
1146  , tbot & ! ice bottom surface temperature (deg C)
1147  , fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
1148  , qsn & ! snow enthalpy
1149  , tsn & ! internal snow temperature
1150  , tsf & ! ice/snow surface temperature, Tsfcn
1151  , fsensn & ! surface downward sensible heat (W m-2)
1152  , flatn & ! surface downward latent heat (W m-2)
1153  , fswabsn & ! shortwave absorbed by ice (W m-2)
1154  , flwoutn & ! upward LW at surface (W m-2)
1155  , fswthrun & ! SW through ice to ocean (W m-2)
1156  , fhnetn & ! fbot, corrected for any surplus energy
1157  , fsurf & ! net flux to top surface, not including fcondtop
1158  , fcondtop & ! downward cond flux at top surface (W m-2)
1159  , fcondbot & ! downward cond flux at bottom surface (W m-2)
1160  , fswint ! SW absorbed in ice interior, below surface (W m-2)
1161 
1162  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
1163  intent(inout) :: &
1164  qin & ! ice layer enthalpy (J m-3)
1165  , tin ! internal ice layer temperatures
1166 
1167  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in):: &
1168  einit ! initial energy of melting (J m-2)
1169 !
1170 !EOP
1171 !
1172  integer (kind=int_kind), parameter :: &
1173  nitermax = 50 ! max number of iterations in temperature solver
1174 
1175  real (kind=dbl_kind), parameter :: &
1176  alph = c3i & ! constant used to get 2nd order accurate fluxes
1177  , bet = -p333 & ! constant used to get 2nd order accurate fluxes
1178  , tsf_errmax = 5.e-4 ! max allowed error in Tsf
1179  ! recommend Tsf_errmax < 0.01 K
1180 
1181  integer (kind=int_kind) :: &
1182  i, j & ! horizontal indices
1183  , ij & ! horizontal index, combines i and j loops
1184  , k & ! ice layer index
1185  , nmat & ! matrix dimension
1186  , niter & ! iteration counter in temperature solver
1187  , isolve ! number of cells with temps not converged
1188 
1189  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
1190  indxii, indxjj ! compressed indices for cells not converged
1191 
1192  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
1193  tsn_init & ! Tsn at beginning of time step
1194  , tsn_start & ! Tsn at start of iteration
1195  , tsf_start & ! Tsf at start of iteration
1196  , dtsf & ! Tsf - Tsf_start
1197  , dtsf_prev & ! dTsf from previous iteration
1198  , dfsurf_dt & ! derivative of fsurf wrt Tsf
1199  , dfsens_dt & ! deriv of fsens wrt Tsf (W m-2 deg-1)
1200  , dflat_dt & ! deriv of flat wrt Tsf (W m-2 deg-1)
1201  , dflwout_dt & ! deriv of flwout wrt Tsf (W m-2 deg-1)
1202  , fswsfc ! SW absorbed at ice/snow surface (W m-2)
1203 
1204  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
1205  khs & ! ksno / hsn
1206  , etas & ! dt / (rhos * cp_ice * hsn)
1207  , dt_rhoi_hlyr & ! dt/(rhoi*hlyr)
1208  , avg_tin & ! = 1. if Tin averaged w/Tin_start, else = 0.
1209  , enew ! new energy of melting after temp change (J m-2)
1210 
1211  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
1212  tin_init & ! Tin at beginning of time step
1213  , tin_start & ! Tin at start of iteration
1214  , iabs & ! SW absorbed in particular layer (W m-2)
1215  , etai ! dt / (rho * cp * h) for ice layers
1216 
1217  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+1) :: &
1218  khi ! ki / hlyr
1219 
1220  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+2) :: &
1221  diag & ! diagonal matrix elements
1222  , sbdiag & ! sub-diagonal matrix elements
1223  , spdiag & ! super-diagonal matrix elements
1224  , rhs & ! rhs of tri-diagonal matrix equation
1225  , tmat ! matrix output temperatures
1226 
1227  real (kind=dbl_kind) :: &
1228  ci & ! specific heat of sea ice (J kg-1 deg-1)
1229  , spdiag2 & ! term to right of superdiag on top row
1230  , avg_tsf & ! = 1. if Tsf averaged w/Tsf_start, else = 0.
1231  , avg_tsn & ! = 1. if Tsn averaged w/Tsn_start, else = 0.
1232  , ferr & ! energy conservation error (W m-2)
1233  , w1,w2,w3,w4,w5,w6,w7 ! temporary variables
1234 
1235  logical (kind=log_kind), dimension (ilo:ihi,jlo:jhi) :: &
1236  converged ! = true when local solution has converged
1237 
1238  logical (kind=log_kind) :: &
1239  all_converged ! = true when all cells have converged
1240 
1241  !-----------------------------------------------------------------
1242  ! Initialize
1243  !-----------------------------------------------------------------
1244 
1245  all_converged = .false.
1246 
1247 !DIR$ CONCURRENT !Cray
1248 !cdir nodep !NEC
1249 !ocl novrec !Fujitsu
1250  do ij = 1, icells
1251 
1252  i = indxi(ij)
1253  j = indxj(ij)
1254 
1255  converged(i,j) = .false.
1256  khs(i,j) = c0i
1257  dtsf(i,j) = c0i
1258  dtsf_prev(i,j) = c0i
1259  tsf_start(i,j) = c0i
1260  enew(i,j) = c0i
1261 
1262  dfsurf_dt(i,j) = c0i
1263  dfsens_dt(i,j) = c0i
1264  dflat_dt(i,j) = c0i
1265  dflwout_dt(i,j) = c0i
1266 
1267  fhnetn(i,j) = fbot(i,j) ! ocean energy used by the ice, <= 0
1268 
1269  fswsfc(i,j) = c0i
1270 ! dt_rhoi_hlyr(i,j) = dt / (rhoi*hlyr(i,j)) ! used by matrix solver
1271  dt_rhoi_hlyr(i,j) = dtice / (rhoi*hlyr(i,j)) ! used by matrix solver
1272 
1273  if (hsn(i,j) > hsnomin) then
1274 ! etas(i,j) = dt / (rhos*cp_ice*hsn(i,j)) ! used by matrix solver
1275  etas(i,j) = dtice / (rhos*cp_ice*hsn(i,j)) ! used by matrix solver
1276  else
1277  etas(i,j) = c0i
1278  endif
1279 
1280  tsn_init(i,j) = tsn(i,j) ! beginning of time step
1281  tsn_start(i,j) = tsn(i,j) ! beginning of iteration
1282 
1283  enddo ! ij
1284 
1285  do k = 1, nilyr
1286 !DIR$ CONCURRENT !Cray
1287 !cdir nodep !NEC
1288 !ocl novrec !Fujitsu
1289  do ij = 1, icells
1290  i = indxi(ij)
1291  j = indxj(ij)
1292 
1293  iabs(i,j,k) = c0i
1294  tin_init(i,j,k) = tin(i,j,k) ! beginning of time step
1295  tin_start(i,j,k) = tin(i,j,k) ! beginning of iteration
1296  enddo ! ij
1297  enddo ! k
1298 
1299  !-----------------------------------------------------------------
1300  ! Compute thermal conductivity at interfaces (held fixed during
1301  ! the subsequent iteration).
1302  !-----------------------------------------------------------------
1303  call conductivity &
1304  (icells, indxi, indxj, &
1305  hlyr, hsn, tin, tbot, khi, khs)
1306 
1307  !-----------------------------------------------------------------
1308  ! Compute solar radiation absorbed in ice and penetrating to ocean
1309  !-----------------------------------------------------------------
1310  call absorbed_solar &
1311  (ni, icells, indxi, indxj, &
1312  hlyr, hsn, fswsfc, fswint, fswthrun, iabs)
1313 
1314  !-----------------------------------------------------------------
1315  ! Solve for new temperatures.
1316  ! Iterate until temperatures converge with minimal energy error.
1317  !-----------------------------------------------------------------
1318 
1319  do niter = 1, nitermax
1320 
1321  if (all_converged) then ! thermo calculation is done
1322  exit
1323  else ! identify cells not yet converged
1324  isolve = 0
1325  do ij = 1, icells
1326  i = indxi(ij)
1327  j = indxj(ij)
1328  if (.not.converged(i,j)) then
1329  isolve = isolve + 1
1330  indxii(isolve) = i
1331  indxjj(isolve) = j
1332  endif
1333  enddo ! ij
1334  endif
1335 
1336  !-----------------------------------------------------------------
1337  ! Update radiative and turbulent fluxes and their derivatives
1338  ! with respect to Tsf.
1339  !-----------------------------------------------------------------
1340 
1341  call surface_fluxes &
1342  (isolve, indxii, indxjj, &
1343  tsf, fswsfc, &
1344  flwoutn, fsensn, flatn, fsurf, &
1345  dflwout_dt, dfsens_dt, dflat_dt, dfsurf_dt)
1346 
1347 !DIR$ CONCURRENT !Cray
1348 !cdir nodep !NEC
1349 !ocl novrec !Fujitsu
1350  do ij = 1, isolve
1351  i = indxii(ij) ! NOTE: not indxi and indxj
1352  j = indxjj(ij)
1353 
1354  !-----------------------------------------------------------------
1355  ! Compute conductive flux at top surface, fcondtop.
1356  ! If fsurf < fcondtop and Tsf = 0, then reset Tsf to slightly less
1357  ! than zero (but not less than -puny).
1358  !-----------------------------------------------------------------
1359 
1360  if (hsn(i,j) > hsnomin) then
1361  fcondtop(i,j) = c2i * khs(i,j) * (tsf(i,j) - tsn(i,j))
1362  else
1363  fcondtop(i,j) = khi(i,j,1) &
1364  * (alph*(tsf(i,j) - tin(i,j,1)) &
1365  + bet*(tsf(i,j) - tin(i,j,2)))
1366  endif
1367  if (fsurf(i,j) < fcondtop(i,j)) &
1368  tsf(i,j) = min(tsf(i,j), -puny)
1369 
1370  !-----------------------------------------------------------------
1371  ! Save surface temperature at start of iteration
1372  !-----------------------------------------------------------------
1373  tsf_start(i,j) = tsf(i,j)
1374 
1375  enddo ! ij
1376 
1377  !-----------------------------------------------------------------
1378  ! Compute specific heat of each ice layer
1379  !-----------------------------------------------------------------
1380 
1381  do k = 1, nilyr
1382  do ij = 1, isolve
1383  i = indxii(ij)
1384  j = indxjj(ij)
1385 
1386  if (l_brine) then
1387  ci = cp_ice - lfresh*tmlt(k) / &
1388  (tin_init(i,j,k)*tin(i,j,k))
1389  else
1390  ci = cp_ice
1391  endif
1392  etai(i,j,k) = dt_rhoi_hlyr(i,j) / ci
1393 
1394  enddo ! ij
1395  enddo ! k
1396 
1397 !DIR$ CONCURRENT !Cray
1398 !cdir nodep !NEC
1399 !ocl novrec !Fujitsu
1400  do ij = 1, isolve
1401  i = indxii(ij)
1402  j = indxjj(ij)
1403 
1404  !-----------------------------------------------------------------
1405  ! Compute matrix elements
1406  !
1407  ! Four possible cases to solve:
1408  ! (1) Cold surface (Tsf < 0), snow present
1409  ! (2) Melting surface (Tsf = 0), snow present
1410  ! (3) Cold surface (Tsf < 0), no snow
1411  ! (4) Melting surface (Tsf = 0), no snow
1412  !-----------------------------------------------------------------
1413 
1414  !-----------------------------------------------------------------
1415  ! first row
1416  ! Tsf for case 1; dummy equation for cases 2, 3 and 4
1417  !-----------------------------------------------------------------
1418 
1419  sbdiag(i,j,1) = c0i
1420  spdiag2 = c0i
1421 
1422  if (hsn(i,j) > hsnomin) then
1423  if (tsf(i,j) <= -puny) then
1424  spdiag(i,j,1) = c2i*khs(i,j)
1425  diag(i,j,1) = dfsurf_dt(i,j) - c2i*khs(i,j)
1426  rhs(i,j,1) = dfsurf_dt(i,j)*tsf(i,j) - fsurf(i,j)
1427  else
1428  spdiag(i,j,1) = c0i
1429  diag(i,j,1) = c1i
1430  rhs(i,j,1) = c0i
1431  endif
1432  else ! hsn <= hsnomin
1433  if (tsf(i,j) <= -puny) then
1434  spdiag(i,j,1) = c0i
1435  diag(i,j,1) = c1i
1436  rhs(i,j,1) = c0i
1437  else
1438  spdiag(i,j,1) = c0i
1439  diag(i,j,1) = c1i
1440  rhs(i,j,1) = c0i
1441  endif
1442  endif
1443 
1444  !-----------------------------------------------------------------
1445  ! second row
1446  ! Tsn for cases 1 and 2; Tsf for case 3; dummy for case 4
1447  !-----------------------------------------------------------------
1448 
1449  if (hsn(i,j) > hsnomin) then
1450  if (tsf(i,j) <= -puny) then
1451  sbdiag(i,j,2) = -etas(i,j) * c2i * khs(i,j)
1452  spdiag(i,j,2) = -etas(i,j) * khi(i,j,1)
1453  spdiag2 = c0i
1454  diag(i,j,2) = c1i &
1455  + etas(i,j) * (c2i*khs(i,j) + khi(i,j,1))
1456  rhs(i,j,2) = tsn_init(i,j)
1457  else
1458  sbdiag(i,j,2) = c0i
1459  spdiag(i,j,2) = -etas(i,j) * khi(i,j,1)
1460  spdiag2 = c0i
1461  diag(i,j,2) = c1i &
1462  + etas(i,j) * (c2i*khs(i,j) + khi(i,j,1))
1463  rhs(i,j,2) = tsn_init(i,j) &
1464  + etas(i,j)*c2i*khs(i,j)*tsf(i,j)
1465  endif
1466  else ! hsn <= hsnomin
1467  if (tsf(i,j) <= -puny) then
1468  sbdiag(i,j,2) = c0i
1469  spdiag(i,j,2) = alph * khi(i,j,1)
1470  spdiag2 = bet * khi(i,j,1)
1471  diag(i,j,2) = dfsurf_dt(i,j) - alph*khi(i,j,1) &
1472  - bet*khi(i,j,1)
1473  rhs(i,j,2) = dfsurf_dt(i,j)*tsf(i,j) - fsurf(i,j)
1474  else
1475  sbdiag(i,j,2) = c0i
1476  spdiag(i,j,2) = c0i
1477  spdiag2 = c0i
1478  diag(i,j,2) = c1i
1479  rhs(i,j,2) = c0i
1480  endif
1481  endif
1482 
1483  !-----------------------------------------------------------------
1484  ! third row, Tin(i,j,1)
1485  !
1486  ! For each internal ice layer compute the specific heat ci,
1487  ! a function of the starting temperature and of the latest
1488  ! guess for the final temperature.
1489  !-----------------------------------------------------------------
1490 
1491  if (hsn(i,j) > hsnomin) then
1492  if (tsf(i,j) <= -puny) then
1493  sbdiag(i,j,3) = -etai(i,j,1) * khi(i,j,1)
1494  spdiag(i,j,3) = -etai(i,j,1) * khi(i,j,2)
1495  diag(i,j,3) = c1i &
1496  + etai(i,j,1)*(khi(i,j,1) + khi(i,j,2))
1497  rhs(i,j,3) = tin_init(i,j,1) &
1498  + etai(i,j,1)*iabs(i,j,1)
1499  else
1500  sbdiag(i,j,3) = -etai(i,j,1) * khi(i,j,1)
1501  spdiag(i,j,3) = -etai(i,j,1) * khi(i,j,2)
1502  diag(i,j,3) = c1i &
1503  + etai(i,j,1)*(khi(i,j,1) + khi(i,j,2))
1504  rhs(i,j,3) = tin_init(i,j,1) &
1505  + etai(i,j,1)*iabs(i,j,1)
1506  endif
1507  else ! hsn <= hsnomin
1508  if (tsf(i,j) <= -puny) then
1509  sbdiag(i,j,3) = -(alph+bet) * etai(i,j,1) * khi(i,j,1)
1510  spdiag(i,j,3) = -etai(i,j,1) &
1511  * (khi(i,j,2) - bet*khi(i,j,1))
1512  diag(i,j,3) = c1i + etai(i,j,1) &
1513  * (khi(i,j,2) + alph*khi(i,j,1))
1514  rhs(i,j,3) = tin_init(i,j,1) &
1515  + etai(i,j,1)*iabs(i,j,1)
1516  else
1517  sbdiag(i,j,3) = c0i
1518  spdiag(i,j,3) = -etai(i,j,1) &
1519  * (khi(i,j,2) - bet*khi(i,j,1))
1520  diag(i,j,3) = c1i + etai(i,j,1) &
1521  * (khi(i,j,2) + alph*khi(i,j,1))
1522  rhs(i,j,3) = tin_init(i,j,1) &
1523  + etai(i,j,1)*iabs(i,j,1)
1524 !!! & + (alph+bet)*etai(i,j,k)*khi(i,j,1)*Tsf(i,j)
1525  ! commented line needed if Tsf were in Kelvin
1526  endif
1527  endif
1528 
1529  if (hsn(i,j) <= hsnomin .and. tsf(i,j) <= -puny) then ! case 3
1530  ! remove spdiag2 in row 2 by adding rows 2 and 3
1531  diag(i,j,2) = spdiag(i,j,3) * diag(i,j,2) &
1532  - spdiag2 * sbdiag(i,j,3)
1533  spdiag(i,j,2) = spdiag(i,j,3) * spdiag(i,j,2) &
1534  - spdiag2 * diag(i,j,3)
1535  rhs(i,j,2) = spdiag(i,j,3) * rhs(i,j,2) &
1536  - spdiag2 * rhs(i,j,3)
1537  endif
1538 
1539  !-----------------------------------------------------------------
1540  ! Bottom row, Tin(i,j,nilyr)
1541  !-----------------------------------------------------------------
1542 
1543  k = nilyr
1544  sbdiag(i,j,k+2) = -etai(i,j,k) &
1545  * (khi(i,j,k) - bet*khi(i,j,k+1))
1546  spdiag(i,j,k+2) = c0i
1547  diag(i,j,k+2) = c1i + etai(i,j,k) &
1548  * (khi(i,j,k) + alph*khi(i,j,k+1))
1549  rhs(i,j,k+2) = tin_init(i,j,k) + etai(i,j,k)*iabs(i,j,k) &
1550  + (alph+bet)*etai(i,j,k) &
1551  *khi(i,j,k+1)*tbot(i,j)
1552 
1553  enddo ! ij
1554 
1555  !-----------------------------------------------------------------
1556  ! Ice interior
1557  !-----------------------------------------------------------------
1558  do k = 2, nilyr-1
1559  do ij = 1, isolve
1560  i = indxii(ij)
1561  j = indxjj(ij)
1562 
1563  sbdiag(i,j,k+2) = -etai(i,j,k) * khi(i,j,k)
1564  spdiag(i,j,k+2) = -etai(i,j,k) * khi(i,j,k+1)
1565  diag(i,j,k+2) = c1i - sbdiag(i,j,k+2) - spdiag(i,j,k+2)
1566  rhs(i,j,k+2) = tin_init(i,j,k) &
1567  + etai(i,j,k) * iabs(i,j,k)
1568 
1569  enddo ! ij
1570  enddo ! k
1571 
1572  !-----------------------------------------------------------------
1573  ! Solve tridiagonal matrix to obtain the new temperatures.
1574  !-----------------------------------------------------------------
1575 
1576  nmat = nilyr + 2
1577 
1578  call tridiag_solver &
1579  (isolve, indxii, indxjj, &
1580  nmat, sbdiag, diag, spdiag, rhs, tmat)
1581 
1582  !-----------------------------------------------------------------
1583  ! Reload temperatures from matrix solution.
1584  !-----------------------------------------------------------------
1585 
1586  do ij = 1, isolve
1587  i = indxii(ij)
1588  j = indxjj(ij)
1589 
1590  if (hsn(i,j) > hsnomin) then
1591  if (tsf(i,j) <= -puny) then
1592  tsf(i,j) = tmat(i,j,1)
1593  tsn(i,j) = tmat(i,j,2)
1594  else
1595  tsf(i,j) = c0i
1596  tsn(i,j) = tmat(i,j,2)
1597  endif
1598  else ! hsn <= hsnomin
1599  if (tsf(i,j) <= -puny) then
1600  tsf(i,j) = tmat(i,j,2)
1601  else
1602  tsf(i,j) = c0i
1603  endif
1604  endif
1605 
1606  enddo
1607 
1608  do k = 1, nilyr
1609  do ij = 1, isolve
1610  i = indxii(ij)
1611  j = indxjj(ij)
1612 
1613  tin(i,j,k) = tmat(i,j,k+2)
1614  enddo
1615  enddo
1616 
1617  !-----------------------------------------------------------------
1618  ! Determine whether the computation has converged to an acceptable
1619  ! solution. Five conditions must be satisfied:
1620  !
1621  ! (1) Tsf <= 0 C.
1622  ! (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
1623  ! dTsf(niter-1) have magnitudes greater than puny, then
1624  ! dTsf(niter)/dTsf(niter-1) cannot be a negative number
1625  ! with magnitude greater than 0.5.
1626  ! (3) abs(dTsf) < Tsf_errmax
1627  ! (4) If Tsf = 0 C, then the downward turbulent/radiative
1628  ! flux, fsurf, must be greater than or equal to the downward
1629  ! conductive flux, fcondtop.
1630  ! (5) The net energy added to the ice per unit time must equal
1631  ! the net change in internal ice energy per unit time,
1632  ! within the prescribed error ferrmax.
1633  !
1634  ! For briny ice (the standard case), Tsn and Tin are limited
1635  ! to prevent them from exceeding their melting temperatures.
1636  ! (Note that the specific heat formula for briny ice assumes
1637  ! that T < Tmlt.)
1638  ! For fresh ice there is no limiting, since there are cases
1639  ! when the only convergent solution has Tsn > 0 and/or Tin > 0.
1640  ! Above-zero temperatures are then reset to zero (with melting
1641  ! to conserve energy) in the thickness_changes subroutine.
1642  !-----------------------------------------------------------------
1643 
1644  ! initialize global convergence flag
1645  all_converged = .true.
1646 
1647 !DIR$ CONCURRENT !Cray
1648 !cdir nodep !NEC
1649 !ocl novrec !Fujitsu
1650  do ij = 1, isolve
1651  i = indxii(ij)
1652  j = indxjj(ij)
1653 
1654  if (l_brine) tsn(i,j) = min(tsn(i,j), c0i)
1655 
1656  !-----------------------------------------------------------------
1657  ! Initialize convergence flag (true until proven false), dTsf,
1658  ! and temperature-averaging coefficients.
1659  ! Average only if test 1 or 2 fails.
1660  !-----------------------------------------------------------------
1661 
1662  converged(i,j) = .true.
1663  dtsf(i,j) = tsf(i,j) - tsf_start(i,j)
1664  avg_tsf = c0i
1665  avg_tsn = c0i
1666  avg_tin(i,j) = c0i
1667 
1668  !-----------------------------------------------------------------
1669  ! Condition 1: check for Tsf > 0
1670  ! If Tsf > 0, set Tsf = 0, then average Tsn and Tin to force
1671  ! internal temps below their melting temps.
1672  !-----------------------------------------------------------------
1673 
1674  if (tsf(i,j) > puny) then
1675  tsf(i,j) = c0i
1676  dtsf(i,j) = -tsf_start(i,j)
1677  if (l_brine) then ! average with starting temp
1678  avg_tsn = c1i
1679  avg_tin(i,j) = c1i
1680  endif
1681  converged(i,j) = .false.
1682  all_converged = .false.
1683 
1684  !-----------------------------------------------------------------
1685  ! Condition 2: check for oscillating Tsf
1686  ! If oscillating, average all temps to increase rate of convergence.
1687  !-----------------------------------------------------------------
1688 
1689  elseif (niter > 1 &! condition (2)
1690  .and. tsf_start(i,j) <= -puny &
1691  .and. abs(dtsf(i,j)) > puny &
1692  .and. abs(dtsf_prev(i,j)) > puny &
1693  .and. -dtsf(i,j)/(dtsf_prev(i,j)+puny*puny) > p5) then
1694 
1695  if (l_brine) then ! average with starting temp
1696  avg_tsf = c1i
1697  avg_tsn = c1i
1698  avg_tin(i,j) = c1i
1699  endif
1700  dtsf(i,j) = p5 * dtsf(i,j)
1701  converged(i,j) = .false.
1702  all_converged = .false.
1703  endif
1704 
1705  !-----------------------------------------------------------------
1706  ! If condition 1 or 2 failed, average new surface/snow
1707  ! temperatures with their starting values.
1708  ! (No change if both tests passed)
1709  !-----------------------------------------------------------------
1710  tsf(i,j) = tsf(i,j) &
1711  + avg_tsf * p5 * (tsf_start(i,j) - tsf(i,j))
1712  tsn(i,j) = tsn(i,j) &
1713  + avg_tsn * p5 * (tsn_start(i,j) - tsn(i,j))
1714 
1715  !-----------------------------------------------------------------
1716  ! Compute qsn and increment new energy.
1717  !-----------------------------------------------------------------
1718  qsn(i,j) = -rhos * (lfresh - cp_ice*tsn(i,j))
1719  enew(i,j) = hsn(i,j) * qsn(i,j)
1720 
1721  tsn_start(i,j) = tsn(i,j) ! for next iteration
1722 
1723  enddo ! ij
1724 
1725  do k = 1, nilyr
1726 !DIR$ CONCURRENT !Cray
1727 !cdir nodep !NEC
1728 !ocl novrec !Fujitsu
1729  do ij = 1, isolve
1730  i = indxii(ij)
1731  j = indxjj(ij)
1732 
1733  if (l_brine) tin(i,j,k) = min(tin(i,j,k), tmlt(k))
1734 
1735  !-----------------------------------------------------------------
1736  ! If condition 1 or 2 failed, average new ice layer
1737  ! temperatures with their starting values.
1738  !-----------------------------------------------------------------
1739  tin(i,j,k) = tin(i,j,k) &
1740  + avg_tin(i,j)*p5*(tin_start(i,j,k)-tin(i,j,k))
1741 
1742  !-----------------------------------------------------------------
1743  ! Compute qin and increment new energy.
1744  !-----------------------------------------------------------------
1745  qin(i,j,k) = -rhoi * (cp_ice*(tmlt(k)-tin(i,j,k)) &
1746  + lfresh*(c1i-tmlt(k)/tin(i,j,k)) &
1747  - cp_ocn*tmlt(k))
1748  enew(i,j) = enew(i,j) + hlyr(i,j)*qin(i,j,k)
1749 
1750  tin_start(i,j,k) = tin(i,j,k) ! for next iteration
1751 
1752  enddo ! ij
1753  enddo ! k
1754 
1755 !DIR$ CONCURRENT !Cray
1756 !cdir nodep !NEC
1757 !ocl novrec !Fujitsu
1758  do ij = 1, isolve
1759  i = indxii(ij)
1760  j = indxjj(ij)
1761 
1762  !-----------------------------------------------------------------
1763  ! Condition 3: check for large change in Tsf
1764  !-----------------------------------------------------------------
1765 
1766  if (abs(dtsf(i,j)) > tsf_errmax) then
1767  converged(i,j) = .false.
1768  all_converged = .false.
1769  endif
1770 
1771  !-----------------------------------------------------------------
1772  ! Condition 4: check for fsurf < fcondtop with Tsf > 0
1773  !-----------------------------------------------------------------
1774 
1775  fsurf(i,j) = fsurf(i,j) + dtsf(i,j)*dfsurf_dt(i,j)
1776  if (hsn(i,j) > hsnomin) then
1777  fcondtop(i,j) = c2i * khs(i,j) * (tsf(i,j)-tsn(i,j))
1778  else
1779  fcondtop(i,j) = khi(i,j,1) &
1780  * (alph*(tsf(i,j)-tin(i,j,1)) &
1781  + bet*(tsf(i,j)-tin(i,j,2)))
1782  endif
1783 
1784  if (tsf(i,j) > -puny .and. fsurf(i,j) < fcondtop(i,j)) then
1785  converged(i,j) = .false.
1786  all_converged = .false.
1787  endif
1788 
1789  !-----------------------------------------------------------------
1790  ! Condition 5: check for energy conservation error
1791  ! Change in internal ice energy should equal net energy input.
1792  !-----------------------------------------------------------------
1793 
1794  fcondbot(i,j) = khi(i,j,nilyr+1) * &
1795  (alph*(tin(i,j,nilyr) - tbot(i,j)) &
1796  + bet*(tin(i,j,nilyr-1) - tbot(i,j)))
1797 
1798 ! ferr = abs( (enew(i,j)-einit(i,j))/dt &
1799  ferr = abs( (enew(i,j)-einit(i,j))/dtice &
1800  - (fcondtop(i,j) - fcondbot(i,j) + fswint(i,j)) )
1801 
1802  ! factor of 0.9 allows for roundoff errors later
1803  if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5)
1804  converged(i,j) = .false.
1805  all_converged = .false.
1806  endif
1807 
1808  dtsf_prev(i,j) = dtsf(i,j)
1809 
1810  enddo ! ij
1811 
1812  enddo ! temperature iteration
1813 
1814  if (.not.all_converged .and.dbg_set(dbg_sbrio) ) then
1815  do ij = 1, icells
1816  i = indxi(ij)
1817  j = indxj(ij)
1818 
1819  !-----------------------------------------------------------------
1820  ! Check for convergence failures.
1821  !-----------------------------------------------------------------
1822  if (.not.converged(i,j)) then
1823 ! write(nu_diag,*) 'Thermo iteration does not converge,', &
1824 ! 'istep1, my_task, i, j, n:', &
1825 ! istep1, my_task, i, j, ni
1826 ! write(nu_diag,*) 'Ice thickness:', hlyr(i,j)*nilyr
1827 ! write(nu_diag,*) 'Snow thickness:', hsn(i,j)
1828 ! write(nu_diag,*) 'dTsf, Tsf_errmax:',dTsf(i,j),Tsf_errmax
1829 ! write(nu_diag,*) 'Tsf:', Tsf(i,j)
1830 ! write(nu_diag,*) 'fsurf:', fsurf(i,j)
1831 ! write(nu_diag,*) 'fcondtop, fcondbot, fswint', &
1832 ! fcondtop(i,j), fcondbot(i,j), fswint(i,j)
1833 ! write(nu_diag,*) 'Flux conservation error =', ferr
1834 ! write(nu_diag,*) 'Initial snow and ice temperatures:'
1835 ! write(nu_diag,*) Tsn_init(i,j), &
1836 ! (Tin_init(i,j,k),k=1,nilyr)
1837 ! write(nu_diag,*) 'Final temperatures:'
1838 ! write(nu_diag,*) Tsn(i,j), (Tin(i,j,k),k=1,nilyr)
1839 ! stoplabel = 'ice state at stop'
1840 ! call print_state(stoplabel,i,j)
1841 ! call abort_ice('vertical thermo: convergence error')
1842  endif
1843  enddo ! ij
1844  endif ! all_converged
1845 
1846 !DIR$ CONCURRENT !Cray
1847 !cdir nodep !NEC
1848 !ocl novrec !Fujitsu
1849  do ij = 1, icells
1850  i = indxi(ij)
1851  j = indxj(ij)
1852 
1853  ! update fluxes that depend on Tsf
1854  flwoutn(i,j) = flwoutn(i,j) + dtsf(i,j) * dflwout_dt(i,j)
1855  fsensn(i,j) = fsensn(i,j) + dtsf(i,j) * dfsens_dt(i,j)
1856  flatn(i,j) = flatn(i,j) + dtsf(i,j) * dflat_dt(i,j)
1857 
1858  ! absorbed shortwave flux for coupler
1859  fswabsn(i,j) = fswsfc(i,j) + fswint(i,j) + fswthrun(i,j)
1860 
1861  enddo ! ij
1862 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
mod_time::abs
Definition: mod_time.f90:64
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::c3i
real(kind=dbl_kind), parameter c3i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_constants::cp_ice
real(kind=dbl_kind), parameter cp_ice
Definition: ice_constants.f90:85
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
ice_constants::p333
real(kind=dbl_kind), parameter p333
Definition: ice_constants.f90:150
ice_constants::cp_ocn
real(kind=dbl_kind), parameter cp_ocn
Definition: ice_constants.f90:85
Here is the call graph for this function:
Here is the caller graph for this function:

◆ thermo_vertical()

subroutine ice_therm_vertical::thermo_vertical ( )

Definition at line 188 of file ice_therm_vertical.f90.

188 !
189 ! !USES:
190 !
191  use ice_atmo
192 ! use ice_timers
193 ! use ice_ocean
194 ! ggao
195  use ice_work, only: worka, workb
196 
197 !
198 ! !INPUT/OUTPUT PARAMETERS:
199 !
200 !EOP
201 !
202  integer (kind=int_kind) :: &
203  i, j &! horizontal indices
204  , ij &! horizontal index, combines i and j loops
205  , ni &! thickness category index
206  , k &! ice layer index
207  , icells ! number of cells with aicen > puny
208 
209  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
210  indxi, indxj ! compressed indices for cells with aicen > puny
211 
212  real (kind=dbl_kind) :: &
213  dhi & ! change in ice thickness
214  , dhs ! change in snow thickness
215 
216 ! 2D coupler variables (computed for each category, then aggregated)
217 
218  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
219  strxn & ! air/ice zonal strss, (N/m^2)
220  , stryn & ! air/ice merdnl strss, (N/m^2)
221  , fsensn & ! surface downward sensible heat (W/m^2)
222  , flatn & ! surface downward latent heat (W/m^2)
223  , fswabsn & ! shortwave absorbed by ice (W/m^2)
224  , flwoutn & ! upward LW at surface (W/m^2)
225  , evapn & ! flux of vapor, atmos to ice (kg m-2 s-1)
226  , trefn & ! air tmp reference level (K)
227  , qrefn & ! air sp hum reference level (kg/kg)
228  , freshn & ! flux of water, ice to ocean (kg/m^2/s)
229  , fsaltn & ! flux of salt, ice to ocean (kg/m^2/s)
230  , fhnetn & ! fbot corrected for leftover energy (W/m^2)
231  , fswthrun ! SW through ice to ocean (W/m^2)
232 
233 ! 2D state variables (thickness, temperature, enthalpy)
234 
235  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
236  hin & ! ice thickness (m)
237  , hsn & ! snow thickness (m)
238  , hlyr & ! ice layer thickness
239  , hin_init & ! initial value of hin
240  , hsn_init & ! initial value of hsn
241  , hsn_new & ! thickness of new snow (m)
242  , qsn & ! snow enthalpy, qsn < 0 (J m-3)
243  , tsn & ! internal snow temperature (deg C)
244  , tsf ! ice/snow top surface temp, same as Tsfcn (deg C)
245 
246  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
247  hnew & ! new ice layer thicknesses (m)
248  , qin & ! ice layer enthalpy, qin < 0 (J m-3)
249  , tin ! internal ice layer temperatures
250 
251 ! other 2D flux and energy variables
252 
253  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
254  tbot & ! ice bottom surface temperature (deg C)
255  , fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
256  , fsurf & ! net flux to top surface, not including fcondtop
257  , fcondtop & ! downward cond flux at top surface (W m-2)
258  , fcondbot & ! downward cond flux at bottom surface (W m-2)
259  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
260  , einit & ! initial energy of melting (J m-2)
261  , efinal & ! final energy of melting (J m-2)
262  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
263 
264 ! call ice_timer_start(4) ! column model
265 ! call ice_timer_start(5) ! thermodynamics
266 
267  !-----------------------------------------------------------------
268  ! Initialize coupler variables sent to the atmosphere.
269  !
270  ! We cannot initialize the ocean coupler fluxes (fresh, fsalt,
271  ! fhnet, and fswthru) because these may have changed during the
272  ! previous time step after the call to_coupler. The ocean
273  ! coupler fluxes are initialized by init_flux_ocn after
274  ! the call to_coupler.
275  !-----------------------------------------------------------------
276  call init_flux_atm
277 
278  !-----------------------------------------------------------------
279  ! Initialize diagnostic variables for history files.
280  !-----------------------------------------------------------------
281  call init_diagnostics
282 
283  !-----------------------------------------------------------------
284  ! Adjust frzmlt to account for ice-ocean heat fluxes since last
285  ! call to coupler.
286  ! Compute lateral and bottom heat fluxes.
287  !-----------------------------------------------------------------
288  call frzmlt_bottom_lateral (tbot, fbot, rside)
289 
290  !-----------------------------------------------------------------
291  ! Vertical thermodynamics: growth rates, updated ice state, and
292  ! coupler variables
293  !-----------------------------------------------------------------
294 
295  do ni = 1, ncat
296 
297  !-----------------------------------------------------------------
298  ! Save initial ice thickness (used later for ITD remapping)
299  !-----------------------------------------------------------------
300 
301  do j = jlo, jhi
302  do i = ilo, ihi
303  if (aicen(i,j,ni) > puny) then
304  hicen_old(i,j,ni) = vicen(i,j,ni) / aicen(i,j,ni)
305  else
306  hicen_old(i,j,ni) = c0i
307  endif
308  enddo ! i
309  enddo ! j
310 
311 
312  !-----------------------------------------------------------------
313  ! Initialize the single-category versions of coupler fluxes,
314  ! along with other thermodynamic arrays.
315  !-----------------------------------------------------------------
316 
317  call init_thermo_vars (strxn, stryn, trefn, qrefn, &
318  fsensn, flatn, fswabsn, flwoutn, &
319  evapn, freshn, fsaltn, fhnetn, &
320  fswthrun, fsurf, fcondtop, fcondbot, &
321  fswint, einit, efinal, mvap)
322 
323  !-----------------------------------------------------------------
324  ! prep for air-to-ice heat, momentum, radiative and water fluxes
325  !-----------------------------------------------------------------
326 
327  call atmo_boundary_layer (ni, 'ice', tsfcn(ilo:ihi,jlo:jhi,ni),&
328  strxn, stryn, trefn, qrefn, worka, workb)
329 
330  !-----------------------------------------------------------------
331  ! Identify cells with nonzero ice area
332  !-----------------------------------------------------------------
333 
334  icells = 0
335  do j = jlo, jhi
336  do i = ilo, ihi
337  if (aicen(i,j,ni) > puny) then
338  icells = icells + 1
339  indxi(icells) = i
340  indxj(icells) = j
341  endif
342  enddo ! i
343  enddo ! j
344 
345  !-----------------------------------------------------------------
346  ! Compute variables needed for vertical thermo calculation
347  !-----------------------------------------------------------------
348 
349  call init_vertical_profile &
350  (ni, icells, indxi, indxj, &
351  hin, hlyr, hsn, hin_init, hsn_init, &
352  qin, tin, qsn, tsn, tsf, einit)
353 
354  !-----------------------------------------------------------------
355  ! Compute new surface temperature and internal ice and snow
356  ! temperatures
357  !-----------------------------------------------------------------
358 
359  call temperature_changes &
360  (ni, icells, indxi, indxj, &
361  hlyr, hsn, qin, tin, qsn, tsn, &
362  tsf, tbot, fbot, fsensn, flatn, fswabsn, &
363  flwoutn, fswthrun, fhnetn, fsurf, fcondtop, fcondbot, &
364  fswint, einit)
365 
366  !-----------------------------------------------------------------
367  ! Compute growth and/or melting at the top and bottom surfaces.
368  ! Repartition ice into equal-thickness layers, conserving energy.
369  !-----------------------------------------------------------------
370 
371  call thickness_changes &
372  (ni, icells, indxi, indxj, &
373  hin, hlyr, hsn, qin, qsn, &
374  mvap, tbot, fbot, flatn, fhnetn, &
375  fsurf, fcondtop, fcondbot, efinal)
376 
377  !-----------------------------------------------------------------
378  ! Check for energy conservation by comparing the change in energy
379  ! to the net energy input
380  !-----------------------------------------------------------------
381 
382  call conservation_check_vthermo &
383  (ni, icells, indxi, indxj, &
384  fsurf, flatn, fhnetn, fswint, einit, efinal)
385 
386  !-----------------------------------------------------------------
387  ! Let it snow
388  !-----------------------------------------------------------------
389 
390  call add_new_snow &
391  (ni, icells, indxi, indxj, &
392  hsn, hsn_new, qsn, tsf)
393 
394  !-----------------------------------------------------------------
395  ! Compute fluxes of water and salt from ice to ocean
396  !-----------------------------------------------------------------
397 
398  do ij = 1, icells
399  i = indxi(ij)
400  j = indxj(ij)
401 
402  dhi = hin(i,j) - hin_init(i,j)
403  dhs = hsn(i,j) - hsn_init(i,j)
404 
405 ! evapn(i,j) = mvap(i,j)/dt ! mvap < 0 => sublimation,
406  evapn(i,j) = mvap(i,j)/dtice ! mvap < 0 => sublimation,
407  ! mvap > 0 => condensation
408  freshn(i,j) = evapn(i,j) - &
409  (rhoi*dhi + rhos*(dhs-hsn_new(i,j))) / dtice
410  fsaltn(i,j) = -rhoi*dhi*ice_ref_salinity*p001/dtice
411 ! (rhoi*dhi + rhos*(dhs-hsn_new(i,j))) / dt
412 ! fsaltn(i,j) = -rhoi*dhi*ice_ref_salinity*p001/dt
413 
414  enddo ! ij
415 
416  !-----------------------------------------------------------------
417  ! Increment area-weighted fluxes.
418  ! Note: Must not change aicen before calling this subroutine.
419  !-----------------------------------------------------------------
420 
421  call merge_fluxes (ni, &
422  strxn, stryn, fsensn, flatn, fswabsn, &
423  flwoutn, evapn, trefn, qrefn, &
424  freshn, fsaltn, fhnetn, fswthrun)
425 
426  !-----------------------------------------------------------------
427  ! Given the vertical thermo state variables (hin, hsn, qin,
428  ! qsn, compute the new ice state variables (vicen, vsnon,
429  ! eicen, esnon).
430  !-----------------------------------------------------------------
431 
432  call update_state_vthermo &
433  (ni, icells, indxi, indxj, &
434  hin, hsn, qin, qsn, tsf)
435 
436  enddo ! ncat
437 
438  !-----------------------------------------------------------------
439  ! Update mixed layer with heat from ice
440  !-----------------------------------------------------------------
441 
442 ! if (oceanmixed_ice) then
443 ! do j=jlo,jhi
444 ! do i=ilo,ihi
445 ! if (hmix(i,j) > c0i) sst(i,j) = sst(i,j) &
446 ! + (fhnet(i,j)+fswthru(i,j))*dt/(cp_ocn*rhow*hmix(i,j))
447 ! enddo
448 ! enddo
449 ! endif
450 
451 ! call ice_timer_stop(5) ! thermodynamics
452 ! call ice_timer_stop(4) ! column model
453 ! ggao
454 
ice_therm_vertical::hicen_old
real(kind=dbl_kind), dimension(:,:,:), allocatable, save hicen_old
Definition: ice_therm_vertical.f90:84
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::ice_ref_salinity
real(kind=dbl_kind), parameter ice_ref_salinity
Definition: ice_constants.f90:85
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_flux::init_flux_atm
subroutine init_flux_atm
Definition: ice_flux.f90:182
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_state::tsfcn
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save tsfcn
Definition: ice_state.f90:97
ice_work
Definition: ice_work.f90:38
ice_work::worka
real(kind=dbl_kind), dimension(:,:), allocatable worka
Definition: ice_work.f90:61
ice_constants::p001
real(kind=dbl_kind), parameter p001
Definition: ice_constants.f90:150
ice_atmo
Definition: ice_atmo.f90:36
ice_flux::merge_fluxes
subroutine merge_fluxes(ni, strxn, stryn, fsensn, flatn, fswabsn, flwoutn, evapn, Trefn, Qrefn, freshn, fsaltn, fhnetn, fswthrun)
Definition: ice_flux.f90:339
ice_state::vicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vicen
Definition: ice_state.f90:97
ice_atmo::atmo_boundary_layer
subroutine atmo_boundary_layer(ni, sfctype, Tsf, strx, stry, Trf, Qrf, delt, delq)
Definition: ice_atmo.f90:62
ice_model_size::ncat
integer(kind=int_kind), parameter ncat
Definition: ice_model_size.f90:48
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_flux::init_diagnostics
subroutine init_diagnostics
Definition: ice_flux.f90:286
ice_work::workb
real(kind=dbl_kind), dimension(:,:), allocatable workb
Definition: ice_work.f90:61
ice_therm_vertical::rside
real(kind=dbl_kind), dimension(:,:), allocatable, save rside
Definition: ice_therm_vertical.f90:88
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
Here is the call graph for this function:

◆ thickness_changes()

subroutine ice_therm_vertical::thickness_changes ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out)  hin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  hlyr,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  hsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(inout)  qin,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  qsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  mvap,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  Tbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  fbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  flatn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  fhnetn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  fsurf,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  fcondtop,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  fcondbot,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout)  efinal 
)

Definition at line 2349 of file ice_therm_vertical.f90.

2349 !
2350 ! !USES:
2351 !
2352 ! !INPUT/OUTPUT PARAMETERS:
2353 !
2354  integer (kind=int_kind), intent(in) :: &
2355  ni & ! thickness category index
2356  , icells ! number of cells with aicen > puny
2357 
2358  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2359  intent(in) :: &
2360  indxi, indxj ! compressed indices for cells with aicen > puny
2361 
2362  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
2363  fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
2364  , tbot & ! ice bottom surface temperature (deg C)
2365  , flatn & ! surface downward latent heat (W m-2)
2366  , fsurf & ! net flux to top surface, not including fcondtop
2367  , fcondtop & ! downward cond flux at top surface (W m-2)
2368  , fcondbot ! downward cond flux at bottom surface (W m-2)
2369 
2370  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
2371  intent(inout) :: &
2372  qin ! ice layer enthalpy (J m-3)
2373 
2374  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2375  fhnetn & ! fbot, corrected for any surplus energy
2376  , hlyr & ! ice layer thickness
2377  , hsn & ! snow thickness (m)
2378  , qsn & ! snow enthalpy
2379  , efinal & ! final energy of melting (J m-2)
2380  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
2381 
2382  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
2383  hin ! total ice thickness
2384 !
2385 !EOP
2386 !
2387  real (kind=dbl_kind), parameter :: &
2388  qbotmax = -p5*rhoi*lfresh ! max enthalpy of ice growing at bottom
2389 
2390  integer (kind=int_kind) :: &
2391  i, j & ! horizontal indices
2392  , ij & ! horizontal index, combines i and j loops
2393  , k, kold ! ice layer indices
2394 
2395  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2396  esub &! energy for sublimation, > 0 (J m-2)
2397  , etop_mlt &! energy for top melting, > 0 (J m-2)
2398  , ebot_mlt &! energy for bottom melting, > 0 (J m-2)
2399  , rhlyr ! reciprocal ice layer thickness
2400 
2401  real (kind=dbl_kind) :: &
2402  dhi & ! change in ice thickness
2403  , dhs & ! change in snow thickness
2404  , ti & ! ice temperature
2405  , ts & ! snow temperature
2406  , econ & ! energy for condensation, < 0 (J m-2)
2407  , ebot_gro & ! energy for bottom growth, < 0 (J m-2)
2408  , qbot & ! enthalpy of ice growing at bottom surface (J m-3)
2409  , qsub & ! energy/unit volume to sublimate ice/snow (J m-3)
2410  , hqtot & ! sum of h*q for two layers
2411  , w1 & ! temporary variable
2412  , hovlp ! overlap between old and new layers (m)
2413 
2414  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
2415  hnew & ! new layer thickness
2416  , hq ! h * q for a layer
2417 
2418  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,0:nilyr) :: &
2419  zold & ! depth of layer boundaries (m)
2420  , znew ! adjusted depths, with equal hlyr (m)
2421 
2422  !-----------------------------------------------------------------
2423  ! Initialize ice layer thicknesses
2424  !-----------------------------------------------------------------
2425 
2426  do k = 1, nilyr
2427  do ij = 1, icells
2428  i = indxi(ij)
2429  j = indxj(ij)
2430  hnew(i,j,k) = hlyr(i,j)
2431  enddo ! if
2432  enddo
2433 
2434  !-----------------------------------------------------------------
2435  ! For l_brine = false (fresh ice), check for temperatures > 0.
2436  ! Melt ice or snow as needed to bring temperatures back to 0.
2437  ! For l_brine = true, this should not be necessary.
2438  !-----------------------------------------------------------------
2439 
2440  if (.not. l_brine) then
2441 !DIR$ CONCURRENT !Cray
2442 !cdir nodep !NEC
2443 !ocl novrec !Fujitsu
2444  do ij = 1, icells
2445  i = indxi(ij)
2446  j = indxj(ij)
2447 
2448  ts = (lfresh + qsn(i,j)/rhos) / cp_ice
2449  if (ts > c0i) then
2450  dhs = cp_ice*ts*hsn(i,j) / lfresh
2451  hsn(i,j) = hsn(i,j) - dhs
2452  qsn(i,j) = -rhos*lfresh
2453  endif
2454  enddo
2455 
2456  do k = 1, nilyr
2457 !DIR$ CONCURRENT !Cray
2458 !cdir nodep !NEC
2459 !ocl novrec !Fujitsu
2460  do ij = 1, icells
2461  i = indxi(ij)
2462  j = indxj(ij)
2463 
2464  ti = (lfresh + qin(i,j,k)/rhoi) / cp_ice
2465  if (ti > c0i) then
2466  dhi = cp_ice*ti*hnew(i,j,k) / lfresh
2467  hnew(i,j,k) = hnew(i,j,k) - dhi
2468  qin(i,j,k) = -rhoi*lfresh
2469  endif
2470  enddo ! ij
2471  enddo ! k
2472 
2473  endif ! .not. l_brine
2474 
2475  !-----------------------------------------------------------------
2476  ! Sublimation/condensation of ice/snow
2477  !-----------------------------------------------------------------
2478 
2479 !DIR$ CONCURRENT !Cray
2480 !cdir nodep !NEC
2481 !ocl novrec !Fujitsu
2482  do ij = 1, icells
2483  i = indxi(ij)
2484  j = indxj(ij)
2485 
2486 ! w1 = -flatn(i,j) * dt
2487  w1 = -flatn(i,j) * dtice
2488  esub(i,j) = max(w1, c0i) ! energy for sublimation, > 0
2489  econ = min(w1, c0i) ! energy for condensation, < 0
2490 
2491  !--------------------------------------------------------------
2492  ! Condensation (mvap > 0)
2493  !--------------------------------------------------------------
2494 
2495  if (hsn(i,j) > puny) then ! add snow with enthalpy qsn
2496  dhs = econ / (qsn(i,j) - rhos*lvap) ! econ < 0, dhs > 0
2497  hsn(i,j) = hsn(i,j) + dhs
2498  mvap(i,j) = mvap(i,j) + dhs*rhos
2499  else ! add ice with enthalpy qin(i,j,1)
2500  dhi = econ / (qin(i,j,1) - rhoi*lvap) ! econ < 0, dhi > 0
2501  hnew(i,j,1) = hnew(i,j,1) + dhi
2502  mvap(i,j) = mvap(i,j) + dhi*rhoi
2503  endif
2504 
2505  !--------------------------------------------------------------
2506  ! Sublimation of snow (mvap < 0)
2507  !--------------------------------------------------------------
2508  qsub = qsn(i,j) - rhos*lvap ! qsub < 0
2509  dhs = max(-hsn(i,j), esub(i,j)/qsub) ! esub > 0, dhs < 0
2510  hsn(i,j) = hsn(i,j) + dhs
2511  esub(i,j) = esub(i,j) - dhs*qsub
2512  esub(i,j) = max(esub(i,j), c0i) ! in case of roundoff error
2513  mvap(i,j) = mvap(i,j) + dhs*rhos
2514  enddo ! ij
2515 
2516  !--------------------------------------------------------------
2517  ! Sublimation of ice (mvap < 0)
2518  !--------------------------------------------------------------
2519 
2520  do k = 1, nilyr
2521 !DIR$ CONCURRENT !Cray
2522 !cdir nodep !NEC
2523 !ocl novrec !Fujitsu
2524  do ij = 1, icells
2525  i = indxi(ij)
2526  j = indxj(ij)
2527  qsub = qin(i,j,k) - rhoi*lvap ! qsub < 0
2528  dhi = max(-hnew(i,j,k), esub(i,j)/qsub) ! esub < 0, dhi < 0
2529  hnew(i,j,k) = hnew(i,j,k) + dhi
2530  esub(i,j) = esub(i,j) - dhi*qsub
2531  esub(i,j) = max(esub(i,j), c0i)
2532  mvap(i,j) = mvap(i,j) + dhi*rhoi
2533  enddo ! ij
2534  enddo ! k
2535 
2536  !-----------------------------------------------------------------
2537  ! Top melt
2538  ! (There is no top growth because there is no liquid water to
2539  ! freeze at the top.)
2540  !-----------------------------------------------------------------
2541 
2542  !--------------------------------------------------------------
2543  ! Melt snow
2544  !--------------------------------------------------------------
2545 !DIR$ CONCURRENT !Cray
2546 !cdir nodep !NEC
2547 !ocl novrec !Fujitsu
2548  do ij = 1, icells
2549  i = indxi(ij)
2550  j = indxj(ij)
2551 
2552 ! w1 = (fsurf(i,j) - fcondtop(i,j)) * dt
2553  w1 = (fsurf(i,j) - fcondtop(i,j)) * dtice
2554  etop_mlt(i,j) = max(w1, c0i) ! etop_mlt > 0
2555  dhs = max(-hsn(i,j), etop_mlt(i,j)/qsn(i,j)) ! qsn < 0, dhs < 0
2556  hsn(i,j) = hsn(i,j) + dhs
2557  etop_mlt(i,j) = etop_mlt(i,j) - dhs*qsn(i,j)
2558  etop_mlt(i,j) = max(etop_mlt(i,j), c0i) ! in case of roundoff error
2559 
2560  ! history diagnostics
2561  if (dhs < -puny .and. mlt_onset(i,j) < puny) &
2562  mlt_onset(i,j) = yday
2563 
2564  enddo ! ij
2565 
2566  !--------------------------------------------------------------
2567  ! Melt ice
2568  !--------------------------------------------------------------
2569 
2570  do k = 1, nilyr
2571 !DIR$ CONCURRENT !Cray
2572 !cdir nodep !NEC
2573 !ocl novrec !Fujitsu
2574  do ij = 1, icells
2575  i = indxi(ij)
2576  j = indxj(ij)
2577 
2578  dhi = max(-hnew(i,j,k), etop_mlt(i,j)/qin(i,j,k))
2579  hnew(i,j,k) = hnew(i,j,k) + dhi ! qin < 0, dhi < 0
2580  etop_mlt(i,j) = etop_mlt(i,j) - dhi*qin(i,j,k)
2581  etop_mlt(i,j) = max(etop_mlt(i,j), c0i)
2582 
2583  ! history diagnostics
2584  if (dhi < -puny .and. mlt_onset(i,j) < puny) &
2585  mlt_onset(i,j) = yday
2586  meltt(i,j) = meltt(i,j) - dhi*aicen(i,j,ni)
2587 
2588  enddo ! ij
2589  enddo ! k
2590 
2591  !----------------------------------------------------------------
2592  ! Bottom growth and melt
2593  !----------------------------------------------------------------
2594 
2595 !DIR$ CONCURRENT !Cray
2596 !cdir nodep !NEC
2597 !ocl novrec !Fujitsu
2598  do ij = 1, icells
2599  i = indxi(ij)
2600  j = indxj(ij)
2601 
2602 ! w1 = (fcondbot(i,j) - fbot(i,j)) * dt
2603  w1 = (fcondbot(i,j) - fbot(i,j)) * dtice
2604  ebot_mlt(i,j) = max(w1, c0i) ! ebot_mlt > 0
2605  ebot_gro = min(w1, c0i) ! ebot_gro < 0
2606 
2607  !--------------------------------------------------------------
2608  ! Grow ice
2609  !--------------------------------------------------------------
2610 
2611  ! enthalpy of new ice growing at bottom surface
2612  qbot = -rhoi * (cp_ice * (tmlt(nilyr+1)-tbot(i,j)) &
2613  + lfresh * (c1i-tmlt(nilyr+1)/tbot(i,j)) &
2614  - cp_ocn * tmlt(nilyr+1))
2615  qbot = min(qbot, qbotmax) ! in case Tbot is close to Tmlt
2616  dhi = ebot_gro / qbot ! dhi > 0
2617 
2618  hqtot = hnew(i,j,nilyr)*qin(i,j,nilyr) + dhi*qbot
2619  hnew(i,j,nilyr) = hnew(i,j,nilyr) + dhi
2620 
2621  if (hnew(i,j,nilyr) > puny) &
2622  qin(i,j,nilyr) = hqtot / hnew(i,j,nilyr)
2623 
2624  ! history diagnostics
2625  congel(i,j) = congel(i,j) + dhi*aicen(i,j,ni)
2626  if (dhi > puny .and. frz_onset(i,j) < puny) &
2627  frz_onset(i,j) = yday
2628 
2629  enddo ! ij
2630 
2631 
2632  !--------------------------------------------------------------
2633  ! Melt ice
2634  !--------------------------------------------------------------
2635 
2636  do k = nilyr, 1, -1
2637 !DIR$ CONCURRENT !Cray
2638 !cdir nodep !NEC
2639 !ocl novrec !Fujitsu
2640  do ij = 1, icells
2641  i = indxi(ij)
2642  j = indxj(ij)
2643 
2644  dhi = max(-hnew(i,j,k), ebot_mlt(i,j)/qin(i,j,k))
2645  hnew(i,j,k) = hnew(i,j,k) + dhi ! qin < 0, dhi < 0
2646  ebot_mlt(i,j) = ebot_mlt(i,j) - dhi*qin(i,j,k)
2647  ebot_mlt(i,j) = max(ebot_mlt(i,j), c0i)
2648 
2649  ! history diagnostics
2650  meltb(i,j) = meltb(i,j) - dhi*aicen(i,j,ni)
2651 
2652  enddo ! ij
2653  enddo ! k
2654 
2655  !--------------------------------------------------------------
2656  ! Melt snow (only if all the ice has melted)
2657  !--------------------------------------------------------------
2658 
2659 !DIR$ CONCURRENT !Cray
2660 !cdir nodep !NEC
2661 !ocl novrec !Fujitsu
2662  do ij = 1, icells
2663  i = indxi(ij)
2664  j = indxj(ij)
2665 
2666  dhs = max(-hsn(i,j), ebot_mlt(i,j)/qsn(i,j))
2667  hsn(i,j) = hsn(i,j) + dhs ! qsn < 0, dhs < 0
2668  ebot_mlt(i,j) = ebot_mlt(i,j) - dhs*qsn(i,j)
2669  ebot_mlt(i,j) = max(ebot_mlt(i,j), c0i)
2670 
2671  ! snow-to-ice conversion (very rare)
2672  if (hsn(i,j) > puny .and. ebot_mlt(i,j) > puny*lfresh) then
2673  hnew(i,j,1) = hsn(i,j) * rhos/rhoi ! reduce thickness
2674  qin(i,j,1) = qsn(i,j) * rhoi/rhos ! increase q to conserve energy
2675  hsn(i,j) = c0i
2676  endif
2677 
2678  enddo ! ij
2679 
2680  !-----------------------------------------------------------------
2681  ! Give the ocean any energy left over after melting
2682  !-----------------------------------------------------------------
2683 
2684 !DIR$ CONCURRENT !Cray
2685 !cdir nodep !NEC
2686 !ocl novrec !Fujitsu
2687  do ij = 1, icells
2688  i = indxi(ij)
2689  j = indxj(ij)
2690 
2691  fhnetn(i,j) = fhnetn(i,j) &
2692  + (esub(i,j) + etop_mlt(i,j) + ebot_mlt(i,j))/dtice
2693 ! + (esub(i,j) + etop_mlt(i,j) + ebot_mlt(i,j))/dt
2694 
2695  enddo ! ij
2696 
2697 !---!-------------------------------------------------------------------
2698 !---! Repartition the ice into equal-thickness layers, conserving energy.
2699 !---!-------------------------------------------------------------------
2700 
2701  !-----------------------------------------------------------------
2702  ! Initialize the new ice thickness.
2703  ! Initialize the final energy: snow energy plus energy of
2704  ! sublimated/condensed ice.
2705  !-----------------------------------------------------------------
2706 
2707  do ij = 1, icells
2708  i = indxi(ij)
2709  j = indxj(ij)
2710  hin(i,j) = c0i
2711  efinal(i,j) = hsn(i,j)*qsn(i,j) - mvap(i,j)*lvap
2712  enddo
2713 
2714  !-----------------------------------------------------------------
2715  ! Compute the new ice thickness.
2716  ! Initialize h*q for new layers.
2717  !-----------------------------------------------------------------
2718 
2719  do k = 1, nilyr
2720 !DIR$ CONCURRENT !Cray
2721 !cdir nodep !NEC
2722 !ocl novrec !Fujitsu
2723  do ij = 1, icells
2724  i = indxi(ij)
2725  j = indxj(ij)
2726  hin(i,j) = hin(i,j) + hnew(i,j,k)
2727  hq(i,j,k) = c0i
2728  enddo ! ij
2729  enddo ! k
2730 
2731  !-----------------------------------------------------------------
2732  ! Make sure hin > puny.
2733  ! Compute depths zold of old layers (unequal thicknesses).
2734  ! Compute depths znew of new layers (all the same thickness).
2735  !-----------------------------------------------------------------
2736 
2737 !DIR$ CONCURRENT !Cray
2738 !cdir nodep !NEC
2739 !ocl novrec !Fujitsu
2740  do ij = 1, icells
2741  i = indxi(ij)
2742  j = indxj(ij)
2743 
2744  if (hin(i,j) > puny) then
2745  hlyr(i,j) = hin(i,j) / real(nilyr,kind=dbl_kind)
2746  rhlyr(i,j) = c1i / hlyr(i,j)
2747  else
2748  hin(i,j) = c0i
2749  hlyr(i,j) = c0i
2750  rhlyr(i,j) = c0i
2751  endif
2752 
2753  zold(i,j,0) = c0i
2754  znew(i,j,0) = c0i
2755 
2756  zold(i,j,nilyr) = hin(i,j)
2757  znew(i,j,nilyr) = hin(i,j)
2758 
2759  enddo ! ij
2760 
2761  do k = 1, nilyr-1
2762  do ij = 1, icells
2763  i = indxi(ij)
2764  j = indxj(ij)
2765 
2766  zold(i,j,k) = zold(i,j,k-1) + hnew(i,j,k) ! old unequal layers
2767  znew(i,j,k) = znew(i,j,k-1) + hlyr(i,j) ! new equal layers
2768 
2769  end do ! ij
2770  enddo ! k
2771 
2772  !-----------------------------------------------------------------
2773  ! Compute h*q for new layer (k) given overlap with old layers (kold)
2774  !-----------------------------------------------------------------
2775 
2776  do k = 1, nilyr
2777  do kold = 1, nilyr
2778 !DIR$ CONCURRENT !Cray
2779 !cdir nodep !NEC
2780 !ocl novrec !Fujitsu
2781  do ij = 1, icells
2782  i = indxi(ij)
2783  j = indxj(ij)
2784 
2785  hovlp = min(zold(i,j,kold), znew(i,j,k)) &
2786  - max(zold(i,j,kold-1), znew(i,j,k-1))
2787  hovlp = max(hovlp, c0i)
2788 
2789  hq(i,j,k) = hq(i,j,k) + hovlp*qin(i,j,kold)
2790 
2791  enddo ! ij
2792  enddo ! kold
2793  enddo ! k
2794 
2795  !-----------------------------------------------------------------
2796  ! Compute new values of qin and increment ice-snow energy.
2797  ! Note: Have to finish previous loop before updating qin
2798  ! in this loop.
2799  !-----------------------------------------------------------------
2800 
2801  do k = 1, nilyr
2802 !DIR$ CONCURRENT !Cray
2803 !cdir nodep !NEC
2804 !ocl novrec !Fujitsu
2805  do ij = 1, icells
2806  i = indxi(ij)
2807  j = indxj(ij)
2808 
2809  qin(i,j,k) = hq(i,j,k) * rhlyr(i,j)
2810  efinal(i,j) = efinal(i,j) + hlyr(i,j)*qin(i,j,k)
2811 
2812  enddo ! ij
2813  enddo !
2814 
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::lvap
real(kind=dbl_kind), parameter lvap
Definition: ice_constants.f90:85
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_flux::congel
real(kind=dbl_kind), dimension(:,:), allocatable, save congel
Definition: ice_flux.f90:140
ice_constants::cp_ice
real(kind=dbl_kind), parameter cp_ice
Definition: ice_constants.f90:85
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_flux::meltt
real(kind=dbl_kind), dimension(:,:), allocatable, save meltt
Definition: ice_flux.f90:140
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_flux::meltb
real(kind=dbl_kind), dimension(:,:), allocatable, save meltb
Definition: ice_flux.f90:140
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_flux::frz_onset
real(kind=dbl_kind), dimension(:,:), allocatable, save frz_onset
Definition: ice_flux.f90:140
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
ice_flux::mlt_onset
real(kind=dbl_kind), dimension(:,:), allocatable, save mlt_onset
Definition: ice_flux.f90:140
ice_constants::cp_ocn
real(kind=dbl_kind), parameter cp_ocn
Definition: ice_constants.f90:85
Here is the caller graph for this function:

◆ tridiag_solver()

subroutine ice_therm_vertical::tridiag_solver ( integer (kind=int_kind), intent(in)  isolve,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxii,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxjj,
integer (kind=int_kind), intent(in)  nmat,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), intent(in)  sbdiag,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), intent(in)  diag,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), intent(in)  spdiag,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), intent(in)  rhs,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), intent(out)  xout 
)

Definition at line 2245 of file ice_therm_vertical.f90.

2245 !
2246 ! !USES:
2247 !
2248 ! !INPUT/OUTPUT PARAMETERS:
2249 !
2250  integer (kind=int_kind), intent(in) :: &
2251  isolve ! number of cells with temps not converged
2252 
2253  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)),&
2254  intent(in) :: &
2255  indxii, indxjj ! compressed indices for cells not converged
2256 
2257  integer (kind=int_kind), intent(in) :: &
2258  nmat ! matrix dimension
2259 
2260  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), &
2261  intent(in) :: &
2262  diag & ! diagonal matrix elements
2263  , sbdiag & ! sub-diagonal matrix elements
2264  , spdiag & ! super-diagonal matrix elements
2265  , rhs ! rhs of tri-diagonal matrix eqn.
2266 
2267  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), &
2268  intent(out) :: &
2269  xout ! solution vector
2270 !
2271 !EOP
2272 !
2273  integer (kind=int_kind) :: &
2274  i, j & ! horizontal indices
2275  , ij & ! horizontal index, combines i and j loops
2276  , k ! row counter
2277 
2278  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2279  wbeta ! temporary matrix variable
2280 
2281  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat) :: &
2282  wgamma ! temporary matrix variable
2283 
2284 
2285 !DIR$ CONCURRENT !Cray
2286 !cdir nodep !NEC
2287 !ocl novrec !Fujitsu
2288  do ij = 1, isolve
2289  i = indxii(ij)
2290  j = indxjj(ij)
2291 
2292  wbeta(i,j) = diag(i,j,1)
2293  xout(i,j,1) = rhs(i,j,1) / wbeta(i,j)
2294 
2295  enddo ! ij
2296 
2297  do k = 2, nmat
2298 !DIR$ CONCURRENT !Cray
2299 !cdir nodep !NEC
2300 !ocl novrec !Fujitsu
2301  do ij = 1, isolve
2302  i = indxii(ij)
2303  j = indxjj(ij)
2304 
2305  wgamma(i,j,k) = spdiag(i,j,k-1) / wbeta(i,j)
2306  wbeta(i,j) = diag(i,j,k) - sbdiag(i,j,k)*wgamma(i,j,k)
2307  xout(i,j,k) = (rhs(i,j,k) - sbdiag(i,j,k)*xout(i,j,k-1)) &
2308  / wbeta(i,j)
2309 
2310  enddo ! ij
2311  enddo ! k
2312 
2313  do k = nmat-1, 1, -1
2314 !DIR$ CONCURRENT !Cray
2315 !cdir nodep !NEC
2316 !ocl novrec !Fujitsu
2317  do ij = 1, isolve
2318  i = indxii(ij)
2319  j = indxjj(ij)
2320 
2321  xout(i,j,k) = xout(i,j,k) - wgamma(i,j,k+1)*xout(i,j,k+1)
2322 
2323  enddo ! ij
2324  enddo ! k
2325 
Here is the caller graph for this function:

◆ update_state_vthermo()

subroutine ice_therm_vertical::update_state_vthermo ( integer (kind=int_kind), intent(in)  ni,
integer (kind=int_kind), intent(in)  icells,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxi,
integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), intent(in)  indxj,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  hin,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  hsn,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), intent(in)  qin,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  qsn,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  Tsf 
)

Definition at line 3042 of file ice_therm_vertical.f90.

3042 !
3043 ! !USES:
3044 !
3045 ! !INPUT/OUTPUT PARAMETERS:
3046 !
3047  integer (kind=int_kind), intent(in) :: &
3048  ni & ! thickness category index
3049  , icells ! number of cells with aicen > puny
3050 
3051  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
3052  intent(in) :: &
3053  indxi, indxj ! compressed indices for cells with aicen > puny
3054 
3055 
3056  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: &
3057  hin & ! ice thickness (m)
3058  , hsn & ! snow thickness (m)
3059  , qsn & ! snow enthalpy
3060  , tsf ! ice/snow surface temperature, Tsfcn
3061 
3062  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
3063  intent(in) :: &
3064  qin ! ice layer enthalpy (J m-3)
3065 !
3066 !EOP
3067 !
3068  integer (kind=int_kind) :: &
3069  i, j & ! horizontal indices
3070  , ij & ! horizontal index, combines i and j loops
3071  , k ! ice layer index
3072 
3073 !DIR$ CONCURRENT !Cray
3074 !cdir nodep !NEC
3075 !ocl novrec !Fujitsu
3076  do ij = 1, icells
3077  i = indxi(ij)
3078  j = indxj(ij)
3079 
3080  if (hin(i,j) > c0i) then
3081  ! aicen is already up to date
3082  vicen(i,j,ni) = aicen(i,j,ni) * hin(i,j)
3083  vsnon(i,j,ni) = aicen(i,j,ni) * hsn(i,j)
3084  esnon(i,j,ni) = qsn(i,j) * vsnon(i,j,ni)
3085  tsfcn(i,j,ni) = tsf(i,j)
3086  else ! (hin(i,j) == c0i)
3087  aicen(i,j,ni) = c0i
3088  vicen(i,j,ni) = c0i
3089  vsnon(i,j,ni) = c0i
3090  esnon(i,j,ni) = c0i
3091  tsfcn(i,j,ni) = tf(i,j)
3092  endif
3093 
3094  enddo ! ij
3095 
3096  do k = 1, nilyr
3097 !DIR$ CONCURRENT !Cray
3098 !cdir nodep !NEC
3099 !ocl novrec !Fujitsu
3100  do ij = 1, icells
3101  i = indxi(ij)
3102  j = indxj(ij)
3103 
3104  if (hin(i,j) > c0i) then
3105  eicen(i,j,ilyr1(ni)+k-1) = &
3106  qin(i,j,k) * vicen(i,j,ni)/real(nilyr,kind=dbl_kind)
3107  else
3108  eicen(i,j,ilyr1(ni)+k-1) = c0i
3109  endif
3110 
3111  enddo ! ij
3112  enddo ! k
3113 
ice_flux::tf
real(kind=dbl_kind), dimension(:,:), allocatable, save tf
Definition: ice_flux.f90:91
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_state::esnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save esnon
Definition: ice_state.f90:97
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_state::tsfcn
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save tsfcn
Definition: ice_state.f90:97
ice_state::vicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vicen
Definition: ice_state.f90:97
ice_state::eicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save eicen
Definition: ice_state.f90:108
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_itd::ilyr1
integer(kind=int_kind), dimension(ncat) ilyr1
Definition: ice_itd.f90:62
ice_state::vsnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vsnon
Definition: ice_state.f90:97
Here is the caller graph for this function:

Variable Documentation

◆ ferrmax

real (kind=dbl_kind), parameter ice_therm_vertical::ferrmax = 1.0e-3_dbl_kind

Definition at line 73 of file ice_therm_vertical.f90.

73  real (kind=dbl_kind), parameter :: &
74  ferrmax = 1.0e-3_dbl_kind & ! max allowed energy flux error (W m-2)
75  ! rec&ommend ferrmax < 0.01 W m-2
76  , hsnomin = 1.0e-6_dbl_kind ! min thickness for which Tsno computed (m)

◆ hicen_old

real (kind=dbl_kind), dimension (:,:,:), allocatable, save ice_therm_vertical::hicen_old

Definition at line 84 of file ice_therm_vertical.f90.

84  real (kind=dbl_kind), dimension (:,:,:),allocatable,save ::&
85  hicen_old ! old ice thick ness (m), used by ice_therm_itd
ice_therm_vertical::hicen_old
real(kind=dbl_kind), dimension(:,:,:), allocatable, save hicen_old
Definition: ice_therm_vertical.f90:84

◆ hsnomin

real (kind=dbl_kind), parameter ice_therm_vertical::hsnomin = 1.0e-6_dbl_kind

Definition at line 73 of file ice_therm_vertical.f90.

◆ l_brine

logical (kind=log_kind) ice_therm_vertical::l_brine

Definition at line 91 of file ice_therm_vertical.f90.

91  logical (kind=log_kind) :: &
92  l_brine ! if true, treat brine pocket effects

◆ rside

real (kind=dbl_kind), dimension (:,:), allocatable, save ice_therm_vertical::rside

Definition at line 88 of file ice_therm_vertical.f90.

88  real (kind=dbl_kind), dimension (:,:) ,allocatable,save ::&
89  rside ! fraction of i ce that melts laterally
ice_therm_vertical::rside
real(kind=dbl_kind), dimension(:,:), allocatable, save rside
Definition: ice_therm_vertical.f90:88

◆ salin

real (kind=dbl_kind), dimension(nilyr+1) ice_therm_vertical::salin

Definition at line 78 of file ice_therm_vertical.f90.

78  real (kind=dbl_kind) :: &
79  salin(nilyr+1) &! salinity (ppt)
80  , tmlt(nilyr+1) &! melting temp, -mu * salinity
81  , ustar_scale ! scaling for i ce-ocean heat flux
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48

◆ tmlt

real (kind=dbl_kind), dimension(nilyr+1) ice_therm_vertical::tmlt

Definition at line 78 of file ice_therm_vertical.f90.

◆ ustar_scale

real (kind=dbl_kind) ice_therm_vertical::ustar_scale

Definition at line 78 of file ice_therm_vertical.f90.

Generated by   doxygen 1.8.15