Update 3 files

- /mo_grotz.f90
- /mo_layer_dynamics.f90
- /mo_mass.f90

mo_grotz.f90

@@ -287,7 +287,7 @@ CONTAINS
             !##########################################################################################
             !Brine flux due to Expulsion
             !##########################################################################################
-            CALL expulsion_flux (thick, V_ex, Nlayer, N_active, psi_g, fl_m, m)
+            CALL expulsion_flux (thick, V_ex, Nlayer, N_active, psi_g, fl_m, m, H_abs)
             IF (i .NE. 1) THEN !is disabled on first timestep to avoid expulsion due to unbalanced initial conditions
                 !This ensures that the salinity is not effected if the initial mass of the layers is to high
                 CALL mass_transfer  (Nlayer, N_active, T, H_abs, S_abs, S_bu, T_bottom, S_bu_bottom, fl_m)

mo_layer_dynamics.f90

@@ -185,7 +185,7 @@ CONTAINS
   !!
   !! @par Revision History
   !! Pasted by Philipp Griewank, IMPRS (2011-05-10)
-  !!
+  !! Edited by Niels Fuchs, UHH (2024-05-03), changed thick_0 to thick(i) to work better with ice core initialization
 
 
   SUBROUTINE top_melt (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
@@ -254,10 +254,10 @@ CONTAINS
        !Removing bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top.
        !Is also activated if thick(N_top+1)=thick_0
        DO k = N_top,N_active-1
-          m(k)          = rho (k+1)        *thick_0
-          S_abs(k)      = S_bu(k+1)        *rho(k+1)*thick_0
-          H_abs(k)      = H (k+1)          *rho(k+1)*thick_0
-          bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick_0
+          m(k)          = rho (k+1)        *thick(k+1) ! thick_0 changed by NF, 2024
+          S_abs(k)      = S_bu(k+1)        *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
+          H_abs(k)      = H (k+1)          *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
+          bgc_temp(k,:) = bgc_bulk (k+1,:) *rho(k+1)*thick(k+1) ! thick_0 changed by NF, 2024
        END DO
        m(N_active)          = 0.0_wp 
        S_abs(N_active)      = 0.0_wp  
@@ -600,7 +600,7 @@ CONTAINS
   !!
   !! @par Revision History
   !! Started by Philipp Griewank, IMPRS (2011-05-10>)
-  !!
+  !! edited by Niels Fuchs, UHH (2024-05-03), changed according to changes in function top_melt
 
 
   SUBROUTINE top_grow (Nlayer,N_active,N_middle,N_top,thick_0,m,S_abs,H_abs,thick,N_bgc,bgc_abs)
@@ -663,10 +663,10 @@ CONTAINS
        !Adding new bottom layer and moving all other layers if N_active<Nlayer .and. N_active>N_top
     ELSE IF (N_active>N_top .AND. N_active<Nlayer) THEN
        DO k = N_top+1,N_active
-          m(k)         = rho (k-1)         *thick_0
-          S_abs(k)     = S_bu(k-1)         *rho(k-1)*thick_0 
-          H_abs(k)     = H (k-1)           *rho(k-1)*thick_0 
-          bgc_temp(k,:) = bgc_bulk(k-1,:)  *rho(k-1)*thick_0
+          m(k)         = rho (k-1)         *thick(k-1) !thick_0
+          S_abs(k)     = S_bu(k-1)         *rho(k-1)*thick(k-1) !thick_0 
+          H_abs(k)     = H (k-1)           *rho(k-1)*thick(k-1) !thick_0 
+          bgc_temp(k,:) = bgc_bulk(k-1,:)  *rho(k-1)*thick(k-1) !thick_0
        END DO
        N_active = N_active+1
        m(N_active)         = rho (N_active-1)         *thick_0 

mo_mass.f90

@@ -112,15 +112,17 @@ CONTAINS
   !! @par Revision History
   !! Brought to life by Philipp Griewank, IMPRS (2010-08-24) \n
   !! Simplified by Philipp Griewank, IMPRS (2010-11-27)
-  SUBROUTINE expulsion_flux (thick,V_ex,Nlayer,N_active,psi_g,fl_m,m)
+  !! edited by Niels Fuchs, UHH (2024-05-03), added brine mass removal from H_abs
+  SUBROUTINE expulsion_flux (thick,V_ex,Nlayer,N_active,psi_g,fl_m,m, H_abs)
 
     INTEGER,                        INTENT(in)    :: Nlayer, N_active
     REAL(wp), DIMENSION(Nlayer),    INTENT(in)    :: V_ex,thick
-    REAL(wp), DIMENSION(Nlayer),    INTENT(inout) :: psi_g,m
+    REAL(wp), DIMENSION(Nlayer),    INTENT(inout) :: psi_g,m, H_abs
     REAL(wp), DIMENSION(Nlayer+1),  INTENT(out)   :: fl_m
 
     INTEGER::k
 
+    H = H_abs(:)/m(:)
     fl_m(1:Nlayer+1)  =  0._wp
     fl_m(2)           = -V_ex(1)*rho_l
     DO k = 2,N_active
@@ -134,6 +136,7 @@ CONTAINS
 
     DO k = 1,N_active
        m(k) =  m(k) +fl_m(k+1)-fl_m(k)
+       H_abs(k) = H_abs(k) + (fl_m(k+1)-fl_m(k))*H(k) ! exp. brine mass must be removed from H_abs, temperature loss is calculated in mass transfer subs.
     END DO
 
   END SUBROUTINE expulsion_flux