@@ -547,7 +547,7 @@ subroutine StieglitzSnow_snowrt(tile_lon, tile_lat, &
547547
548548 do i= 1 ,N_snow
549549
550- call StieglitzSnow_calc_tpsnow(htsnn(i),wesn(i),tdum,fdum, ice1(i),tzero(i), ignore_pos_tpsnow= .true. )
550+ call StieglitzSnow_calc_tpsnow(htsnn(i),wesn(i),tdum,fdum, ice1(i),tzero(i), .true. , ignore_pos_tpsnow= .true. )
551551
552552 if (ice1(i)) then
553553 cl(i) = df(i)
@@ -625,7 +625,7 @@ subroutine StieglitzSnow_snowrt(tile_lon, tile_lat, &
625625 - fhsn(i)- df(i)* (dtc(i-1 )- dtc(i))
626626 HTSPRIME= HTSNN(I)+ AREASC* FLXNET* DTS
627627
628- call StieglitzSnow_calc_tpsnow( HTSPRIME, wesn(i), tdum, fnew, logdum, logdum, ignore_pos_tpsnow= .true. )
628+ call StieglitzSnow_calc_tpsnow( HTSPRIME, wesn(i), tdum, fnew, logdum, logdum, .true. , ignore_pos_tpsnow= .true. )
629629
630630 fnew= amax1(0 ., amin1(1 ., fnew))
631631
@@ -650,7 +650,7 @@ subroutine StieglitzSnow_snowrt(tile_lon, tile_lat, &
650650
651651 HTSPRIME= HTSNN(I)+ AREASC* FLXNET* DTS
652652
653- call StieglitzSnow_calc_tpsnow( HTSPRIME, wesn(i), tdum, fnew, logdum, logdum, ignore_pos_tpsnow= .true. )
653+ call StieglitzSnow_calc_tpsnow( HTSPRIME, wesn(i), tdum, fnew, logdum, logdum, .true. , ignore_pos_tpsnow= .true. )
654654
655655 fnew= amax1(0 ., amin1(1 ., fnew))
656656
@@ -875,23 +875,23 @@ subroutine StieglitzSnow_snowrt(tile_lon, tile_lat, &
875875
876876 if (dens(i) > StieglitzSnow_RHOMA) then
877877
878- if (tileType== MAPL_LANDICE) then ! restrict SWE adjustment to LANDICE tiles
878+ ! ! if (tileType==MAPL_LANDICE) then ! restrict SWE adjustment to LANDICE tiles
879879
880- ! excs = SWE in excess of max density given fixed snow depth
881-
882- excs(i) = (dens(i)- StieglitzSnow_RHOMA)* sndz(i) ! solid + liquid
883- wlossfrac= excs(i)/ wesn(i)
884- wesn(i) = wesn(i) - excs(i) ! remove EXCS from SWE
885- do k= 1 ,N_constit
886- rmelt(k)= rmelt(k)+ rconstit(i,k)* wlossfrac/ dts
887- rconstit(i,k)= rconstit(i,k)* (1 .- wlossfrac)
888- rconstit(i,k)= amax1(0 .,rconstit(i,k)) ! guard against truncation error
889- enddo
890- hnew = (StieglitzSnow_CPW* tpsn(i)- fices(i)* alhm)* wesn(i) ! adjust heat content accordingly
891- hcorr= hcorr+ (htsnn(i)- hnew)/ dts ! add excess heat content into residual accounting term
892- htsnn(i)= hnew
893-
894- end if
880+ ! excs = SWE in excess of max density given fixed snow depth
881+
882+ excs(i) = (dens(i)- StieglitzSnow_RHOMA)* sndz(i) ! solid + liquid
883+ wlossfrac= excs(i)/ wesn(i)
884+ wesn(i) = wesn(i) - excs(i) ! remove EXCS from SWE
885+ do k= 1 ,N_constit
886+ rmelt(k)= rmelt(k)+ rconstit(i,k)* wlossfrac/ dts
887+ rconstit(i,k)= rconstit(i,k)* (1 .- wlossfrac)
888+ rconstit(i,k)= amax1(0 .,rconstit(i,k)) ! guard against truncation error
889+ enddo
890+ hnew = (StieglitzSnow_CPW* tpsn(i)- fices(i)* alhm)* wesn(i) ! adjust heat content accordingly
891+ hcorr= hcorr+ (htsnn(i)- hnew)/ dts ! add excess heat content into residual accounting term
892+ htsnn(i)= hnew
893+
894+ ! ! end if
895895
896896 dens(i) = StieglitzSnow_RHOMA
897897 endif
@@ -901,12 +901,12 @@ subroutine StieglitzSnow_snowrt(tile_lon, tile_lat, &
901901
902902 wesndens = wesn - wesndens
903903
904- if (tileType== MAPL_LANDICE) then ! finish SWE adjustment for LANDICE tiles
905-
906- pre = pre + sum (excs* max (1 .- fices,0.0 ))/ dts
907- excs = excs * fices / dts
904+ ! ! if (tileType==MAPL_LANDICE) then ! finish SWE adjustment for LANDICE tiles
908905
909- end if
906+ pre = pre + sum (excs* max (1 .- fices,0.0 ))/ dts
907+ excs = excs * fices / dts
908+
909+ ! ! end if
910910
911911 snowd= sum (wesn)
912912 call StieglitzSnow_calc_asnow( snowd, areasc0 )
@@ -1169,7 +1169,7 @@ subroutine StieglitzSnow_relayer(N_snow, N_constit, tileType, targetthick, &
11691169 if (conserve_ice10_tzero_tmp) then
11701170
11711171 do i= 1 ,N_snow
1172- call StieglitzSnow_calc_tpsnow(htsnn(i),wesn(i),tdum,fdum,ice10(i),tzero0(i),ignore_pos_tpsnow= .true. )
1172+ call StieglitzSnow_calc_tpsnow(htsnn(i),wesn(i),tdum,fdum,ice10(i),tzero0(i),.true. , ignore_pos_tpsnow= .true. )
11731173 enddo
11741174
11751175 end if
@@ -1288,7 +1288,7 @@ subroutine StieglitzSnow_relayer(N_snow, N_constit, tileType, targetthick, &
12881288 !
12891289 ! --+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
12901290
1291- ! diagnose snow temperature and ice fraction
1291+ ! diagnose snow temperature and ice fraction (*_calc_tpsnow_vector always has use_threshold_fac=.false.)
12921292
12931293 call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices, rc= rc_tmp)
12941294
@@ -1302,6 +1302,9 @@ subroutine StieglitzSnow_relayer(N_snow, N_constit, tileType, targetthick, &
13021302
13031303 ! for each layer, check snow conditions (partially/fully frozen, temp at/below zero)
13041304 ! before and after relayer; in select cases, adjust snow heat content and temp
1305+ !
1306+ ! NOTE: logicals before relayer action were computed with "buffer" (use_threshold_fac=.true. )
1307+ ! reals after relayer action were computed without "buffer" (use_threshold_fac=.false.)
13051308
13061309 do i= 1 ,N_snow
13071310
@@ -1341,7 +1344,7 @@ end subroutine StieglitzSnow_relayer
13411344
13421345 ! **********************************************************************
13431346
1344- subroutine StieglitzSnow_calc_tpsnow_scalar ( h , w , t , f , ice1 , tzero , ignore_pos_tpsnow , rc )
1347+ subroutine StieglitzSnow_calc_tpsnow_scalar ( h , w , t , f , ice1 , tzero , use_threshold_fac , ignore_pos_tpsnow , rc )
13451348
13461349 ! diagnose snow temperature and frozen fraction from snow mass and snow heat content
13471350 !
@@ -1371,6 +1374,8 @@ subroutine StieglitzSnow_calc_tpsnow_scalar( h, w, t, f, ice1, tzero, ignore_pos
13711374
13721375 logical , intent (out ) :: ice1, tzero ! frozen fraction==1?, snow temp at 0 deg C?
13731376
1377+ logical , intent (in ) :: use_threshold_fac
1378+
13741379 logical , intent (in ), optional :: ignore_pos_tpsnow
13751380
13761381 integer , intent (out ), optional :: rc ! return code for detecting positive snow heat content and temperature
@@ -1407,21 +1412,21 @@ subroutine StieglitzSnow_calc_tpsnow_scalar( h, w, t, f, ice1, tzero, ignore_pos
14071412
14081413 ! -------------------------------------------------------------------
14091414
1410- ! if (use_threshold_fac) then
1411- !
1412- ! ! replicates original get_tf0d()
1413- !
1414- ! threshold1 = -1.00001*alhm
1415- ! threshold2 = -0.99999*alhm
1416- !
1417- ! else
1415+ if (use_threshold_fac) then
1416+
1417+ ! replicates original get_tf0d()
1418+
1419+ threshold1 = - 1.00001 * alhm
1420+ threshold2 = - 0.99999 * alhm
1421+
1422+ else
14181423
1419- ! replicates original get_tf_nd() / StieglitzSnow_calc_tpsnow[_vector]()
1420-
1421- threshold1 = - alhm
1422- threshold2 = - alhm
1423-
1424- ! end if
1424+ ! replicates original get_tf_nd() / StieglitzSnow_calc_tpsnow[_vector]()
1425+
1426+ threshold1 = - alhm
1427+ threshold2 = - alhm
1428+
1429+ end if
14251430
14261431 ! -------------------------------------------------------------------
14271432
@@ -1544,7 +1549,7 @@ subroutine StieglitzSnow_calc_tpsnow_vector( N, h, w, t, f, ignore_pos_tpsnow, r
15441549 do ii= 1 ,N
15451550
15461551 call StieglitzSnow_calc_tpsnow_scalar( h(ii), w(ii), t(ii), f(ii), ice1, tzero, &
1547- ignore_pos_tpsnow= ignore_pos_tpsnow_tmp, rc= rc_tmp )
1552+ .false. , ignore_pos_tpsnow= ignore_pos_tpsnow_tmp, rc= rc_tmp )
15481553
15491554 if (rc_tmp/= 0 ) then
15501555
0 commit comments