anchor.diff:  
           end if
           
  -        !now%t2m = t2m_new 
  -        ! add remainder of melted snow/ice and refreezing to surface temperature
  -        !where (now%mask == 2)         ! ice
  -        !    now%tsurf = now%tsurf  - (now%melted_snow + now%melted_ice - now%refr)*par%tsticsub*rhow*clm/par%ceff
  -        !elsewhere                     ! land and sea
  -        !    now%tsurf = now%tsurf - (now%melted_snow - now%refr)*par%tsticsub*rhow*clm/par%ceff
  -        !end where
  -
           ! write prognostic and diagnostic output
           now%subl = now%subl/rhow
   
  
positive.code:             now%hsnow = 0.0_dp
          else where
              ! update snow height
              now%hsnow = dmax1(0.0_dp, now%hsnow + now%smb_snow*par%tsticsub)
          end where
          
          ! Relax snow height to maximum (eg, 5 m)
          snow_to_ice   = dmax1(0.d0,now%hsnow-hsmax)
          now%hsnow   = now%hsnow - snow_to_ice
          now%smb_ice = snow_to_ice/par%tsticsub - now%melted_ice + now%refr   ! Use to force ice sheet model
          now%hice    = now%hice + now%smb_ice*par%tsticsub ! update new ice budget: remove or add ice
  
          if (.not. bnd%smb) then
              where (now%mask == 2)
                  now%smb = now%smb_snow + now%smb_ice - snow_to_ice/par%tsticsub
              elsewhere
                  now%smb = now%smb_snow + dmax1(0.0_dp,now%smb_ice - snow_to_ice/par%tsticsub)
              end where
          end if
  
  
negative.code:         allocate(now%lwu(npts))
          allocate(now%runoff(npts))
          allocate(now%evap(npts))
          allocate(now%qmr(npts))
  
          ! forcing fields
          allocate(now%sf(npts))
  