mo_grav_drain resetted

mo_grav_drain.f90

@@ -72,7 +72,8 @@ CONTAINS
   !! Added harmonic mean for permeability by Philipp Griewank (2014-01-05)
 
   SUBROUTINE 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)
+       T_bottom,S_bu_bottom,grav_drain,grav_temp,grav_salt,grav_heat_flag,harmonic_flag, &
+       grav_heat_flux_down, grav_heat_flux_up, fl_brine_bgc)
     INTEGER,                       INTENT(in)    :: Nlayer, N_active,grav_heat_flag,harmonic_flag
     REAL(wp), DIMENSION(Nlayer),   INTENT(in)    :: S_br,S_bu,thick,T,psi_l,psi_s
     REAL(wp),                      INTENT(in)    :: dt,S_bu_bottom,T_bottom
@@ -90,6 +91,7 @@ CONTAINS
     REAL(wp)                                     :: heat_loss             !< Amount of heat transported from the ice [J]
     INTEGER                                      :: k,kk
     REAL(wp), DIMENSION(Nlayer+1,Nlayer+1), INTENT(inout),OPTIONAL :: fl_brine_bgc
+    REAL(wp),                      INTENT(out)   :: grav_heat_flux_down, grav_heat_flux_up
 
     ray_mini       = ray_crit!Arises from the optimization of the r
     heat_loss      = 0._wp
@@ -192,7 +194,8 @@ CONTAINS
     IF (grav_heat_flag==2) THEN
        H_abs(N_active) = H_abs(N_active) +heat_loss-fl_up(N_active)*c_l*T_bottom
     END IF
-
+    grav_heat_flux_down = grav_heat_flux_down + heat_loss
+    grav_heat_flux_up = grav_heat_flux_up - fl_up(N_active)*c_l*T_bottom
 
     IF(MINVAL(S_abs)<0.0_wp) THEN
        PRINT*,'negative salinity after gravity drainige, aborted',S_abs