Index: issm/trunk/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- issm/trunk/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 27347)
+++ issm/trunk/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 28013)
@@ -240,5 +240,5 @@
 
 		// change in pure snow albedo due to soot loading
-		IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
+		IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S2,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
 
 		// determine the effective change due to finite depth and soot loading
Index: issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 27347)
+++ issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 28013)
@@ -67,5 +67,4 @@
 			}
 
-			smb[v]=smb[v]/1000*rho_water/rho_ice;      // SMB in m/y ice
 		}  //end of the loop over the vertices
 
@@ -187,6 +186,6 @@
 	femmodel->parameters->FindParam(&arlagcoefs,&M,&N,SmbARMAarlagcoefsEnum);             _assert_(M==numbasins); _assert_(N==arorder);
    femmodel->parameters->FindParam(&malagcoefs,&M,&N,SmbARMAmalagcoefsEnum);             _assert_(M==numbasins); _assert_(N==maorder);
-   femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);                 _assert_(M==numbasins); _assert_(N==numelevbins);
-   femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);                _assert_(M==numbasins); _assert_(N==numelevbins-1);
+   femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);                 _assert_(M==numbasins); _assert_(N==numelevbins*12);
+   femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);                _assert_(M==numbasins); _assert_(N==(numelevbins-1)*12);
    femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum);                _assert_(M==numbasins);
 
@@ -504,131 +503,9 @@
 
 }/*}}}*/
-void SmbDebrisMLx(FemModel* femmodel){/*{{{*/
-
-	//      The function is based on:
-	//      Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2
-	//      Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7
-	//      function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90
-
-	/*Intermediaries*/
-	// altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002)
-	IssmDouble LW=2.9;          // W/m^2 /100m                       2.9
-	IssmDouble SW=1.3;          // W/m^2 /100m                       1.3
-	IssmDouble HumidityG=0;     // % /100m         rough estimate
-	IssmDouble AirTemp=0.7;     // C /100m
-	IssmDouble WindSpeed=0.02;  // m/s /100m       rough estimate    0.2
-
-	// accumulation follows a linear increase above the ELA up to a plateau
-	IssmDouble AccG=0.1;                    // m w.e. /100m
-	IssmDouble AccMax=1.;                    // m w.e.
-	IssmDouble ReferenceElevation=2200.;     // m M&L
-	IssmDouble AblationDays=120.;            //
-
-	IssmDouble In=100.;                 // Wm^-2        incoming long wave
-	IssmDouble Q=500.;                  // Wm^-2        incoming short wave
-	IssmDouble K=0.585;                // Wm^-1K^-1    thermal conductivity          0.585
-	IssmDouble Qm=0.0012;              // kg m^-3      measured humiditiy level
-	IssmDouble Qh=0.006 ;              // kg m^-3      saturated humidity level
-	IssmDouble Tm=2.;                   // C            air temperature
-	IssmDouble Rhoaa=1.22;             // kgm^-3       air densitiy
-	IssmDouble Um=1.5;                 // ms^-1        measured wind speed
-	IssmDouble Xm=1.5;                 // ms^-1        measurement height
-        IssmDouble Xr=0.01;                // ms^-1        surface roughness             0.01
-        IssmDouble Alphad=0.07;            //              debris albedo                 0.07
-        IssmDouble Alphai=0.4;             //              ice ablbedo
-        IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
-        IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
-        IssmDouble Lm;//=3.34E+05;            // jkg^-1K^-1   latent heat of ice melt
-        IssmDouble Lv=2.50E+06;            // jkg^-1K^-1   latent heat of evaporation
-        IssmDouble Tf=273.;                 // K            water freeezing temperature
-        IssmDouble Eps=0.95;               //              thermal emissivity
-        IssmDouble Rhoi=900.;               // kgm^-3       ice density
-        IssmDouble Sigma=5.67E-08;         // Wm^-2K^-4    Stefan Boltzmann constant
-        IssmDouble Kstar=0.4;              //              von kármán constant
-        IssmDouble Gamma=180.;              // m^-1         wind speed attenuation        234
-	IssmDouble PhiD;//=0.005;              //              debris packing fraction       0.01
-	IssmDouble Humidity=0.2;           //              relative humidity
-
-	IssmDouble smb,yts,z,debris;
-	IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
-
-	/*Get material parameters and constants */
-	//femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as  benchmark was run with different densities for debris and FS
-	femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum);
-	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); 
-	PhiD=0.;
-	//if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
-
-	/* Loop over all the elements of this partition */
-	for(Object* & object : femmodel->elements->objects){
-		Element* element=xDynamicCast<Element*>(object);
-
-		/* Allocate all arrays */
-		int         numvertices=element->GetNumberOfVertices();
-		IssmDouble* surfacelist=xNew<IssmDouble>(numvertices);
-		IssmDouble* smb=xNew<IssmDouble>(numvertices);
-		IssmDouble* debriscover=xNew<IssmDouble>(numvertices);
-		element->GetInputListOnVertices(surfacelist,SurfaceEnum);
-
-		/* Get inputs */
-		element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
-
-		/*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
-		for(int v=0;v<numvertices;v++){
-
-                        /*Get vertex elevation */
-                        z=surfacelist[v];
-
-                        /* Get debris cover */
-                        debris=debriscover[v];
-
-                        /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
-                        IssmDouble n=debris/dk;
-                        IssmDouble nmax=1000;
-                        IssmDouble Alphaeff;
-                        if(n>nmax){
-                                Alphaeff=Alphad;
-                        } else {
-                                Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
-                        }*/
-
-                        // M&L
-                       IssmDouble Alphaeff=Alphad;
-
-                        /* compute smb */
-                        for (int ismb=0;ismb<2;ismb++){
-                                if(ismb==0){
-                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
-                                        debris=0.;
-                                }else{
-                                	// only in the meltregime debris develops
-                                        if(-MassBalanceCmDayDebris<0) debris=debriscover[v]; 
-                                }
-                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
-                                    (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
-                                    (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
-                                    ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
-                                    ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
-                                    K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
-                                    ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
-                                    Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
-                                    -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
-                        }
-
-                        /* account form ablation days, and convert to m/s */
-			MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts;
-
-			/*Update array accordingly*/
-			smb[v]=MassBalanceMYearDebris;
-		}
-
-		/*Add input to element and Free memory*/
-		element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
-		xDelete<IssmDouble>(surfacelist);
-		xDelete<IssmDouble>(smb);
-		xDelete<IssmDouble>(debriscover);
-	}
+void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
+        for(Object* & object : femmodel->elements->objects){
+                Element* element=xDynamicCast<Element*>(object);
+                element->SmbDebrisEvatt();
+        }
 }/*}}}*/
 void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
@@ -641,9 +518,10 @@
 }/*}}}*/
 #ifdef _HAVE_SEMIC_
-void SmbSemicx(FemModel* femmodel){/*{{{*/
-
-	for(Object* & object : femmodel->elements->objects){
-		Element* element=xDynamicCast<Element*>(object);
-		element->SmbSemic();
+void SmbSemicx(FemModel* femmodel,int ismethod){/*{{{*/
+
+	for(Object* & object : femmodel->elements->objects){
+		Element* element=xDynamicCast<Element*>(object);
+		if (ismethod == 1) element->SmbSemicTransient(); // Inwoo's version.
+		else element->SmbSemic(); // original SmbSEMIC
 	}
 
Index: issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 27347)
+++ issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 28013)
@@ -23,7 +23,7 @@
 void SmbMeltComponentsx(FemModel* femmodel);
 void SmbGradientsComponentsx(FemModel* femmodel);
-void SmbDebrisMLx(FemModel* femmodel);
+void SmbDebrisEvattx(FemModel* femmodel);
 /* SEMIC: */
-void SmbSemicx(FemModel* femmodel);
+void SmbSemicx(FemModel* femmodel, int ismethod);
 /*GEMB: */
 void       Gembx(FemModel* femmodel);
Index: issm/trunk/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90
===================================================================
--- issm/trunk/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90	(revision 28013)
+++ issm/trunk/src/c/modules/SurfaceMassBalancex/run_semic_transient.f90	(revision 28013)
@@ -0,0 +1,248 @@
+subroutine run_semic_transient(nx, ntime, nloop, sf_in, rf_in, swd_in, lwd_in, wind_in, &
+      sp_in, rhoa_in, qq_in, tt_in, tsurf_in, qmr_in, &
+      tstic, &
+      hcrit, rcrit, &
+      mask, hice, hsnow, &
+      albedo, albedo_snow, &
+      alb_scheme, alb_smax, alb_smin, albi, albl, &
+      Tamp, &
+      tmin, tmax, tmid, mcrit, w_crit, tau_a, tau_f, afac, verbose, &
+      tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out,  refr_out, alb_out, & 
+      alb_snow_out,hsnow_out,hice_out,qmr_out) !{{{
+
+   use utils
+   use surface_physics
+   implicit none
+
+   ! declare surface physics class
+   type(surface_physics_class) :: surface
+   ! declare forcing class
+   !type(forc_class) :: forc
+   ! declare validation class
+   !type(vali_class) :: vali	! validation not needed here
+
+   integer, parameter:: dp=kind(0.d0)  !< define precision (machine specific)
+   integer :: i, k, n
+   integer :: nnx, nny
+   integer :: year=0
+   integer :: day =1 !< not used value.
+   integer, intent(in) :: nx, ntime, nloop   ! number of grid / number of times
+   logical, intent(in) :: verbose            ! verbosity
+   logical :: debug=.false.
+
+   ! forcing data    
+   ! input argument format array size (nx, ntime)...
+   double precision, intent(in), dimension(nx):: sf_in    ! snow fall rate [m/s]
+   double precision, intent(in), dimension(nx):: rf_in    ! rain fall rate [m/s]
+   double precision, intent(in), dimension(nx):: swd_in   ! downwelling shortwave radiation [W/m2]
+   double precision, intent(in), dimension(nx):: lwd_in   ! downwelling longwave radiation [W/m2]
+   double precision, intent(in), dimension(nx):: wind_in  ! surface wind speed [m/s]
+   double precision, intent(in), dimension(nx):: sp_in    ! surface pressure [Pa]
+   double precision, intent(in), dimension(nx):: rhoa_in  ! air density [kg/m3]
+   double precision, intent(in), dimension(nx):: qq_in    ! air specific humidity [kg/kg]
+   double precision, intent(in), dimension(nx):: tt_in    ! air temperature [K]
+
+   ! input data
+   double precision, intent(in) :: tstic  ! time step from ISSM [sec].
+
+   ! output data
+   ! Ice surface Temperature [K]
+   double precision, intent(out), dimension(nx):: tsurf_out    
+   ! surface mass balance=(Accu-Melt) [m/s]
+   double precision, intent(out), dimension(nx):: smb_out     
+   double precision, intent(out), dimension(nx):: smbi_out     ! SMB ice  [water equivalent m/s]
+   double precision, intent(out), dimension(nx):: smbs_out     ! SMB snow [water equivalent m/s]
+   double precision, intent(out), dimension(nx):: saccu_out    ! accumulation [m/s]
+   double precision, intent(out), dimension(nx):: smelt_out    ! ablation [m/s]
+   double precision, intent(out), dimension(nx):: refr_out     ! freezing [m/s]
+   double precision, intent(out), dimension(nx):: alb_out      ! grid-averaged albedo [no unit] 
+   double precision, intent(out), dimension(nx):: alb_snow_out 
+   double precision, intent(out), dimension(nx):: hice_out    
+   double precision, intent(out), dimension(nx):: hsnow_out   
+   double precision, intent(out), dimension(nx):: qmr_out     
+
+   double precision :: total_time, start, finish
+
+   ! set parameters
+   !character (len=256) :: name         ! not used(?)
+   !character (len=256) :: boundary(30) ! not used(?)
+   !character (len=256), intent(in) :: alb_scheme  !< name of albedo scheme
+   integer, intent(in)          :: alb_scheme
+   !integer :: n_ksub    
+   !double precision             :: ceff         !< surface heat heat capacity of snow/ice [J/W m2]
+   double precision, intent(in), dimension(nx):: albedo
+   double precision, intent(in), dimension(nx):: albedo_snow !< spatial..
+   double precision, intent(in), dimension(nx):: hsnow
+   double precision, intent(in), dimension(nx):: hice
+   double precision, intent(in), dimension(nx):: tsurf_in    !< input temperature [K]
+   double precision, intent(in), dimension(nx):: qmr_in 
+   double precision, intent(in), dimension(nx):: mask
+
+   double precision, intent(in) :: albi
+   double precision, intent(in) :: albl
+   double precision, intent(in) :: alb_smax
+   double precision, intent(in) :: alb_smin
+   double precision, intent(in) :: hcrit !< critical snow height for which grid cell is 50% snow covered [m]
+   double precision, intent(in) :: rcrit !< critical snow height fro which refreezing fraction is 50% [m]
+   double precision, intent(in) :: Tamp
+   !double precision    :: csh
+   !double precision    :: clh
+   double precision, intent(in) :: tmin
+   double precision, intent(in) :: tmax
+   !double precision    :: tsticsub
+   ! parameters for isba albedo scheme.
+   double precision, intent(in) :: tau_a  !< critical liquide water concent for "isba" albedo scheme [kg/m2]
+   double precision, intent(in) :: tau_f
+   double precision, intent(in) :: w_crit
+   double precision, intent(in) :: mcrit
+   double precision, intent(in) :: afac !< param
+   double precision, intent(in) :: tmid !< param for "alex" albedo parameterization [K]
+
+   if (debug) then
+      print*,'   ntime: ', ntime
+      print*,'   nx   : ', nx
+   end if
+
+   ! set vector length
+   surface%par%nx = nx
+
+   ! FIXME should be user input
+   !boundary = "" "" ""
+   if (debug) then
+      print*, "run_semic_transient: initialize parameters."
+   end if
+   surface%par%tstic = tstic      !< time step [s]
+   surface%par%ceff= 2.0e6_dp     !< surface heat capacity of snow/ice [J/K/m2]
+   surface%par%csh = 2.0e-3_dp    !< turbulent heat exchange coefficient 
+   surface%par%clh = 5.0e-4_dp    !< turbulent heat exchange coefficient [no unit]
+   surface%par%alb_smax = alb_smax !0.79_dp !< max snow albedo
+   surface%par%alb_smin = alb_smin !0.6_dp  !< min snow albedo
+   surface%par%albi = albi ! 0.41_dp     !< albedo for ice
+   surface%par%albl = albl ! 0.07_dp     !< albedo for land
+   surface%par%tmin = tmin ! -999_dp
+   surface%par%tmax = tmax ! 273.15_dp
+   surface%par%hcrit = hcrit !0.028_dp   !< critical snow height for which grid cell is 50 % snow colvered [m]
+   surface%par%rcrit = rcrit !0.85_dp    !< refreezing fraction is 50% [m]
+   surface%par%amp   = Tamp !3.0_dp   !< amplitude of diurnal cycle [K]
+   if (alb_scheme == 0) then
+      surface%par%alb_scheme="none"
+   else if (alb_scheme == 1) then
+      surface%par%alb_scheme = "slater"
+   else if (alb_scheme == 2) then
+      surface%par%alb_scheme = "denby"
+   else if (alb_scheme == 3) then
+      surface%par%alb_scheme = "isba"
+   else
+      print*, "ERROR: current albedo scheme is not available."
+      call exit(1)
+   end if 
+   surface%par%tau_a  = tau_a  !0.008_dp
+   surface%par%tau_f  = tau_f  !0.24_dp
+   surface%par%w_crit = w_crit !15.0_dp ! default value
+   surface%par%mcrit  = mcrit  !6.0e-8_dp
+   surface%par%n_ksub = 3      ! sub ...
+   ! snow albedo of alex
+   surface%par%afac   = afac
+   surface%par%tmid   = tmid
+   
+   ! initialize sub-daily time step tsticsub
+   surface%par%tsticsub = surface%par%tstic / dble(surface%par%n_ksub)
+
+   ! allocate necessary arrays for surface_physics module
+   call surface_alloc(surface%now,surface%par%nx)
+
+   ! initialise prognostic variables
+   if (debug) then
+      print*,"run_semic_transient: initialize variables."
+   end if
+   ! these values will be updated through "surface_energy_and_mass_balance" function.
+   surface%now%mask    (:) = mask       (:) ! 2.0_dp  !loi_mask(:nx)
+   if (debug) then
+      print*,"run_semic_transient: initialize variables: mask"
+   end if
+   surface%now%hsnow   (:) = hsnow      (:) ! initial snow height...
+   surface%now%hice    (:) = hice       (:) ! initial ice height..
+   surface%now%tsurf   (:) = tsurf_in   (:) !< initial ice surface temperature
+   surface%now%alb     (:) = albedo     (:) !< initial albedo for energy balance.
+   surface%now%alb_snow(:) = albedo_snow(:) !< initial albedo for ISBA albedo method.
+   if (debug) then
+      print*,"run_semic_transient: initialize variables. DONE."
+   end if
+
+   if (debug) then
+      !print*, "====== global variable =========="
+      !print*, "nloop          :", nloop
+      !print*, "nx             :", surface%par%nx
+      !print*, "======  parameters ======"
+      !print*, "csh            :", surface%par%csh
+      !print*, "clh            :", surface%par%clh
+      !print*, "albeo scheme   :", surface%par%alb_scheme
+      !print*, "albeo ice      :", surface%par%albi
+      !print*, "tstic          :", surface%par%tstic
+      !print*, "tsticsub       :", surface%par%tsticsub
+      !print*, "n_ksub         :", surface%par%n_ksub
+      print*, "====== inputs ========="
+      print*, "hsnow          :", hsnow
+      print*, "======  state variables ======"
+      print*, "hsnow          :", surface%now%hsnow
+      print*, "hice           :", surface%now%hice
+      print*, "albeo          :", surface%now%alb
+      print*, "albeo snow     :", surface%now%alb_snow
+      print*, "mask           :", surface%now%mask
+      print*, "tsurf          :", surface%now%tsurf
+      print*, "sf             :", sf_in
+   end if
+
+   ! define boundary conditions (not used, here!)
+   call surface_boundary_define(surface%bnd,surface%par%boundary)
+   !call print_boundary_opt(surface%bnd)
+
+   ! input with single value
+   do k =1,nloop
+      do i =1,ntime
+         if (debug) then
+            print*,"run_semic_transient: forcing data: ntime = ", i
+         end if
+         surface%now%sf   = sf_in  !(:,i)
+         surface%now%rf   = rf_in  !(:,i)
+         surface%now%sp   = sp_in  !(:,i)
+         surface%now%lwd  = lwd_in !(:,i)
+         surface%now%swd  = swd_in !(:,i)
+         surface%now%wind = wind_in!(:,i)
+         surface%now%rhoa = rhoa_in!(:,i)
+         surface%now%t2m  = tt_in  !(:,i)
+         surface%now%qq   = qq_in  !(:,i)
+         ! qmr_res is used to "energy_balance" in semic.
+         surface%now%qmr_res = qmr_in
+
+         ! calculate prognostic and diagnsotic variables
+         call surface_energy_and_mass_balance(surface%now,surface%par,surface%bnd,i,year)
+         
+         if (debug) then
+            print*,"done..."
+         end if
+         if (k == nloop) then
+            tsurf_out         = surface%now%tsurf
+            ! melt - potential surface melt [m/s]
+            ! smb = SMB_ice + SMB_snow
+            ! smbi  - SMB_ice  (water equivalent m/sec)
+            ! smbs  - SMB_snow (water equivalent m/sec)
+            smb_out           =surface%now%smb      ! smb = smb_snow + smb_ice
+            smbi_out          =surface%now%smb_ice  ! Csi (snow>ice) - melted_ice + refrezon_rain.
+            smbs_out          =surface%now%smb_snow ! smb_snow = snowfall - sublimiation - melted_snow + refrozen_snow
+            saccu_out         =surface%now%acc      ! acc      = snowfall - sublimiation - refreezing 
+            smelt_out         =surface%now%melt     ! potential surface melt = melt_ice + melt_snow
+            refr_out          =surface%now%refr     ! refreezing values. [m/sec]
+            alb_out           =surface%now%alb
+            alb_snow_out      =surface%now%alb_snow
+            hsnow_out         =surface%now%hsnow
+            hice_out          =surface%now%hice
+            qmr_out           =surface%now%qmr_res
+         end if
+      end do
+   end do
+
+   ! de-allocate surface_physics arrays
+   call surface_dealloc(surface%now)
+
+end subroutine run_semic_transient ! }}}
