Merge branch 'debug_summertime_expulsion' into 'master'

Changes in layer dynamics and enthalpy calculations

See merge request !3

mo_data.f90

@@ -21,6 +21,7 @@
 !! @par Revision History
 !! Initialized by Philipp Griewank, IMPRS (2010-07-14) \n
 !! Add several variables by Niels Fuchs, MPIMET (2017-03-01)
+!! edited by Niels Fuchs, UHH (2024-05-03), corrected Enthalpy declaration (before specific and total enthalpy were confused)
 
 !!
 MODULE mo_data
@@ -31,8 +32,8 @@ MODULE mo_data
 
   
   !----Arrays
-  REAL(wp), DIMENSION(:), ALLOCATABLE   :: H                 !<  Enthalpy [J]
-  REAL(wp), DIMENSION(:), ALLOCATABLE   :: H_abs             !<  specific Enthalpy [J/kg]
+  REAL(wp), DIMENSION(:), ALLOCATABLE   :: H                 !<  specific Enthalpy [J/kg]
+  REAL(wp), DIMENSION(:), ALLOCATABLE   :: H_abs             !<  total Enthalpy [J]
   REAL(wp), DIMENSION(:), ALLOCATABLE   :: Q                 !<  Heat in layer [J]
   REAL(wp), DIMENSION(:), ALLOCATABLE   :: fl_Q              !<  Heat flux between layers [J/s]
   REAL(wp), DIMENSION(:), ALLOCATABLE   :: T                 !<  Temperature [C]

mo_grotz.f90

@@ -79,6 +79,7 @@ CONTAINS
     !! Basic thermodynamics and layer_dynamics for fixed boundaries seem stable, backup made. by griewank (2010-08-10) \n
     !! Add some more outputs, changed routine names and arguments with respect to newly introduces flags by Niels Fuchs, MPIMET (2017-03-01) \n
     !! Added a bit of description with the run down of what happends by Philipp Griewank, Uni K (2018-08-08)
+    !! Edited by Niels Fuchs, UHH (2024-05-03), correct layer enthalpy balance in brine expulsion routine
     SUBROUTINE grotz ()
 
         USE mo_parameters
@@ -287,7 +288,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(k) in middle layers to work better with ice core initialization, since middle layer thickness can vary
 
 
   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