mo_grotz resetted

mo_grotz.f90

@@ -227,6 +227,8 @@ CONTAINS
                     CALL snow_precip (m_snow, H_abs_snow, thick_snow, dt, liquid_precip, T2m)
                 ELSE IF(MAX(liquid_precip, solid_precip)>0.0_wp .AND. N_Active==1) THEN !Precip in water layer
                     CALL snow_precip_0 (H_abs(1), S_abs(1), m(1), T(1), dt, liquid_precip, T2m)
+                ELSE IF(MIN(liquid_precip, solid_precip)<0.0_wp .AND. thick_snow>0._wp .AND. N_Active>1) THEN      !Erosion of snow layer
+                    CALL snow_erosion (m_snow, H_abs_snow, thick_snow, dt, liquid_precip, T2m)
                 END IF
 
             ELSE IF (precip_flag==0) THEN
@@ -235,6 +237,8 @@ CONTAINS
                     test = thick_snow
                 ELSE IF (MAX(liquid_precip, solid_precip)>0.0_wp .AND. N_Active==1) THEN !Precip in water layer
                     CALL snow_precip_0 (H_abs(1), S_abs(1), m(1), T(1), dt, liquid_precip, T2m, solid_precip)
+                ELSE IF (MAX(liquid_precip, solid_precip)>0.0_wp .AND. N_Active>1) THEN !Erosion of snowlayer
+                    CALL snow_erosion (m_snow, H_abs_snow, thick_snow, dt, liquid_precip, T2m, solid_precip)
                 END IF
             END IF
 
@@ -330,12 +334,15 @@ CONTAINS
                     grav_drain = grav_drain / time_out
                 end if
 
+                grav_heat_flux_down = grav_heat_flux_down / time_out
+                grav_heat_flux_up = grav_heat_flux_up / time_out
+
                 !Calling standard output
                 CALL output(Nlayer, T, psi_s, psi_l, thick, S_bu, ray, format_T, format_psi, &
                         format_thick, format_snow, freeboard, thick_snow, T_snow, psi_l_snow, psi_s_snow, &
                         energy_stored, freshwater, total_resist, thickness, bulk_salin, &
                         grav_drain, grav_salt, grav_temp, T2m, T_top, perm, format_perm, flush_v, flush_h, psi_g, &
-                        & melt_thick_output, format_melt)
+                        & melt_thick_output, grav_heat_flux_down, grav_heat_flux_up, format_melt)
 
                 IF (bgc_flag==2) THEN
                     CALL output_bgc(Nlayer, N_active, bgc_bottom, N_bgc, bgc_abs, psi_l, thick, m, format_bgc)
@@ -355,6 +362,7 @@ CONTAINS
                         ',  T_snow:', T_snow, &
                         ',  T2m:', T2m
 
+
                 ! if you uncomment snow gt zero and flush, you need the following write statement:
                 ! WRITE(*,'(A10,I2,A15,F4.3,A12,F5.3,A14,F7.3,A30,F3.1,A14,F5.4,A16,L1,A10,F7.3,A7,F7.3,A10,A3)')&
                 flush_question = 'no!'
@@ -366,7 +374,9 @@ CONTAINS
                         'S_bu_bottom: ', s_bu_bottom, &
                         'T_top: ', T_top, &
                         'fl_sw: ', fl_sw, &
-                        'fl_rest; ', fl_rest
+                        'fl_rest; ', fl_rest, &
+						', Habs_top:', H_abs(2), &
+						', Habs_bot:', H_abs(N_active)
 
 
 
@@ -375,7 +385,8 @@ CONTAINS
                 grav_drain = 0.0_wp
                 grav_salt = 0.0_wp
                 grav_temp = 0.0_wp
-
+                grav_heat_flux_down = 0._wp
+                grav_heat_flux_up = 0._wp
                 melt_thick_output(:) = 0._wp
 
                 n_time_out = 0
@@ -452,11 +463,11 @@ CONTAINS
                 IF (bgc_flag == 2) THEN
                     CALL fl_grav_drain (S_br, S_bu, psi_l, psi_s, thick, S_abs, H_abs, T, m, dt, Nlayer, N_active, &
                             &ray, T_bottom, S_bu_bottom, grav_drain, grav_temp, grav_salt, &
-                            grav_heat_flag, harmonic_flag, fl_brine_bgc)
+                            grav_heat_flag, harmonic_flag, grav_heat_flux_down, grav_heat_flux_up, fl_brine_bgc)
                 ELSE
                     CALL fl_grav_drain (S_br, S_bu, psi_l, psi_s, thick, S_abs, H_abs, T, m, dt, Nlayer, N_active, &
                             &ray, T_bottom, S_bu_bottom, grav_drain, grav_temp, grav_salt, &
-                            grav_heat_flag, harmonic_flag)
+                            grav_heat_flag, harmonic_flag, grav_heat_flux_down, grav_heat_flux_up)
                 END IF
 
             ELSE IF (grav_flag==3 .AND. N_active>1) THEN
@@ -794,6 +805,8 @@ CONTAINS
         CLOSE(48)
         CLOSE(49)
         CLOSE(50)
+        CLOSE(51)
+        CLOSE(52)
         CLOSE(66)
 
         !Deallocating arrays