From 29516032d74526b6cdbd84506a181b34f469bf01 Mon Sep 17 00:00:00 2001
From: JakobDeutloff <50237419+JakobDeutloff@users.noreply.github.com>
Date: Mon, 11 Apr 2022 12:37:40 +0200
Subject: [PATCH] mo_grav_drain resetted

---
 mo_grav_drain.f90 | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/mo_grav_drain.f90 b/mo_grav_drain.f90
index 6f074b3..b68f0a9 100644
--- a/mo_grav_drain.f90
+++ b/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
@@ -83,13 +84,14 @@ CONTAINS
     REAL(wp), DIMENSION(N_active)                :: fl_down               !< Downward brine flux [kg] 1 is the flux from 1 to N_active, N_active is from N_active to the ocean
     REAL(wp), DIMENSION(Nlayer)                  :: perm                  !< Permeability
     REAL(wp), DIMENSION(Nlayer)                  :: harmonic_perm         !< Harmonic mean permeability
-    REAL(wp), DIMENSION(Nlayer+1)                :: fl_m  
+    REAL(wp), DIMENSION(Nlayer+1)                :: fl_m
     REAL(wp)                                     :: flux                  !< Downward brine flux [kg]
     REAL(wp)                                     :: test1
     REAL(wp)                                     :: d_S_br,height,ray_mini
     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), 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
@@ -160,7 +162,7 @@ CONTAINS
 
           fl_down(k) = flux
           FORALL(kk=k:N_active)
-             fl_up(kk) = fl_up(kk)+flux 
+             fl_up(kk) = fl_up(kk)+flux
           END FORALL
 
           fl_up(k) = MIN(fl_up(k),psi_l(k)*rho_l*thick(k)) ! should limit flux to brine volume of layer
@@ -178,7 +180,7 @@ CONTAINS
     IF( PRESENT( fl_brine_bgc ) ) THEN
        fl_brine_bgc(1:(N_active-1),N_active+1) = fl_brine_bgc(1:(N_active-1),N_active) + fl_down(1:(N_active-1))
        !fl_brine_bgc(N_active,N_active+1)     = fl_brine_bgc(N_active,N_active+1)     + sum(fl_down(:))
-       DO k = 1,N_active 
+       DO k = 1,N_active
           fl_brine_bgc(k+1,k) = fl_brine_bgc(k+1,k) + fl_up(k)
        END DO
 
@@ -190,9 +192,10 @@ CONTAINS
 
     !Heatflux correction
     IF (grav_heat_flag==2) THEN
-       H_abs(N_active) = H_abs(N_active) +heat_loss-fl_up(N_active)*c_l*T_bottom       
+       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
-- 
GitLab