diff --git a/mo_grav_drain.f90 b/mo_grav_drain.f90 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..6f074b3c21b1c781460f0fb635f964ddbda939fb 100644 --- a/mo_grav_drain.f90 +++ b/mo_grav_drain.f90 @@ -0,0 +1,281 @@ +!> +!! Computes the Salt fluxes caused by gravity drainage. +!! +!! +!! +!! +!! @author Philipp Griewank +!! +!! +!! COPYRIGHT +!! +!! This file is part of SAMSIM. +!! +!! SAMSIM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! +!! SAMSIM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License along with SAMSIM. If not, see <http://www.gnu.org/licenses/>. +!! +! +!! +!! +!! @par Revision History +!! Injected with life by Philipp Griewank, IMPRS (<2010-08-27>) +!! +!! + +MODULE mo_grav_drain + + USE mo_parameters + USE mo_thermo_functions, ONLY:func_S_br + USE mo_mass + IMPLICIT NONE + + PUBLIC :: fl_grav_drain,fl_grav_drain_simple + +CONTAINS + + + !> + !! Calculates fluxes caused by gravity drainage. + !! + !! If the Rayleigh number of a layer is higher then the critical value, brine leaves the layer by a theoretical brine channel. + !! The discharged brine flows downward through all layers directly into the underlying ocean. + !! To preserve mass the same amount of water flows upwards through all lower layers. + !! In contrast to the downward flux the upward flux is assumed to be in thermal equilibrium thus moving salt and heat to each layer. + !! The upward flux is a standard upwind advection. + !! The downward flux of a layer over the timestep is = x*(Ray-Ray_crit)*dt*thick. + !! _____________________________________ + !! | ->__|__<- + !! |______________| v |__________________ + !! | ^ ->| |<- ^ + !! |______________| |__________________ + !! | ^ ^ ->|||||<- ^ ^ + !! |______________|vvv|__________________ + !! ^ ^ ^ ^ ^ ^ + !! This superb ascii art is supposed to show how the assumed fluxes flow downward through a brine channel and upwards through the layer. + !! x and r are passed on to enable easy optimization. + !! The effect of the upward moving brine is calculated in mass_transfer. + !! + !! IMPORTANT: The height assumptions are special. The bottom of the ice edge is assumed to be at psi_s(N_active)/psi_s_min *thick_0 + !! + !! The first approach assumed that brine drainage occurred between two layers but performed poorly. + !! + !! If grav_heat_flag is set to 2 the amount of heat transported out of the ice will be compensated in the lowest layer + !! + !! @par Revision History + !! created by Philipp Griewank, IMPRS (2010-08-27) \n + !! Completely revised to assume brine channels by Philipp Griewank , IMPRS (2010-11-05) \n + !! Mass_transfer is used to advect H and S by Philipp Griewank, IMPRS (2010-11-05) \n + !! Added condition S_br(k)>S_br(k+1) by Philipp Griewank. IMPRS (2011-04-29) \n + !! 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) + 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 + REAL(wp), INTENT(inout) :: grav_drain,grav_temp,grav_salt + REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: S_abs,H_abs,m + REAL(wp), DIMENSION(Nlayer-1), INTENT(out) :: ray !< Rayleigh number + REAL(wp), DIMENSION(N_active) :: fl_up !< Upward brine flux [kg] 1 is the flux from 2 to 1, N_active flux from ocean to N_active + 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) :: 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 + + ray_mini = ray_crit!Arises from the optimization of the r + heat_loss = 0._wp + perm = 0.0_wp + perm(N_active) = 9999999._wp !lowest layer is considered fully liquid + ray = 0._wp + fl_up = 0._wp + fl_down = 0._wp + harmonic_perm = 0.0_wp + + !Permeability is determined + DO k=1,N_active + perm(k)=10.0_wp**(-17.0_wp) * (1000.0_wp*ABS(psi_l(k)))**3.10_wp + END DO + + !Harmonic mean permeability is determined + IF (harmonic_flag == 2) THEN + DO k=1,N_active-1 + test1 = minval(perm(k:N_active-1)) + IF (test1 .lt. 10._wp**(-14._wp)) THEN + harmonic_perm(k)=0.0_wp + ELSE + DO kk=k,N_active-1 + harmonic_perm(k)=harmonic_perm(k)+thick(kk)/perm(kk) + END DO + !bottom layer is included in a linear fashion + harmonic_perm(k)=harmonic_perm(k)+(thick(N_active)*psi_s(N_active)/psi_s_min)/perm(N_active) + harmonic_perm(k)=(SUM(thick(k:N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min)/harmonic_perm(k) + END IF + END DO + END IF + + !___________Raleigh number is calculated (see Notz 2005)__________________________ + DO k=1,N_active-1 + d_S_br = S_br(k)-S_br(N_active) + height = SUM(thick((k+1):N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min !Height adjustment + IF (harmonic_flag == 1) THEN + ray(k) = grav*rho_l*bbeta*d_S_br*height* MINVAL(perm(k:N_active)) + ELSE IF (harmonic_flag == 2) THEN + ray(k) = grav*rho_l*bbeta*d_S_br*height* harmonic_perm(k) + END IF + ray(k) = ray(k)/(kappa_l*mu) + ray(k) = MAX(ray(k),0._wp) + END DO + + + + + grav_salt = grav_salt + SUM(S_abs(:)) + + + DO k = 1,N_active-1 + IF (ray(k)>ray_mini .AND. psi_s(k)>0.001_wp .AND. S_abs(k)/m(k)>0.1_wp .AND. S_br(k)>S_br(k+1)) THEN !if psi_s equal zero then flushing is assumed + flux = x_grav*(ray(k)-ray_mini) *dt*thick(k) + flux = MIN(flux,psi_l(k)*rho_l*thick(k)) !should limit flux to brine volume of layer + S_abs(k) = S_abs(k) -flux*S_br(k) + IF (S_abs(k)<0.0_wp) THEN + PRINT*,'Negative Salinity due to gravity drainige overdrive in layer',k + PRINT*,S_abs(k),S_abs(k)+flux*S_br(k),flux,psi_l(k)*rho_l*thick(k) + STOP 21234 + END IF + + grav_temp = grav_temp + flux*T(k) + + H_abs(k) = H_abs(k) -flux*c_l*T(k) + heat_loss = heat_loss +flux*c_l*T(k) + + + fl_down(k) = flux + FORALL(kk=k:N_active) + 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 + + + END IF + END DO + + grav_salt = grav_salt - SUM(S_abs(:)) + + !Influence of summed up upward transport resulting from the downward drainage + fl_m(1)=0.0_wp + fl_m(2:N_active+1)=fl_up(1:N_active) + + 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 + fl_brine_bgc(k+1,k) = fl_brine_bgc(k+1,k) + fl_up(k) + END DO + + END IF + + CALL mass_transfer (Nlayer,N_active,T,H_abs,S_abs,S_bu,T_bottom,S_bu_bottom,fl_m) + + grav_drain = grav_drain+fl_m(N_active+1) + + !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 + END IF + + + IF(MINVAL(S_abs)<0.0_wp) THEN + PRINT*,'negative salinity after gravity drainige, aborted',S_abs + STOP 1337 + END IF + END SUBROUTINE fl_grav_drain + + + + + + !> + !! Calculates salinity to imitate the effects gravity drainage. + !! + !! Based on the assumption that super critical Rayleigh numbers are quickly reduced below the critical Rayleigh number. + !! Proposed as a very simplified parametrisation of gravity drainage. + !! Includes no fluxes of any kind, instead bulk salinity is simply reduced when ever the Rayleigh number is above the critical values. + !! The parametrization begins from the bottom layers and moves upward. + !! + !! @par Revision History + !! created by Philipp Griewank, IMPRS (2012-01-01) + + SUBROUTINE fl_grav_drain_simple (psi_s,psi_l,thick,S_abs,S_br,Nlayer,N_active,ray, & + grav_drain,harmonic_flag) + INTEGER, INTENT(in) :: Nlayer, N_active, harmonic_flag + REAL(wp), DIMENSION(Nlayer), INTENT(in) :: thick,psi_l,psi_s,S_br + REAL(wp), INTENT(inout) :: grav_drain + REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: S_abs + REAL(wp), DIMENSION(Nlayer-1), INTENT(out) :: ray !< Rayleigh number + REAL(wp), DIMENSION(Nlayer) :: perm !< Permeability + REAL(wp), DIMENSION(Nlayer) :: harmonic_perm !< Harmonic permeability + REAL(wp) :: d_S_br,height,ray_mini,temp + INTEGER :: k,kk + + ray_mini = ray_crit + perm = 0.0_wp + perm(N_active) = 9999999._wp + ray = 0._wp + + !Permeability is determined + DO k = 1,N_active + perm(k)=10.0_wp**(-17.0_wp) * (1000.0_wp*ABS(psi_l(k)))**3.10_wp + END DO + + !Harmonic mean permeability is determined + IF (harmonic_flag == 2) THEN + DO k=1,N_active-1 + temp = minval(perm(k:N_active-1)) + IF (temp .lt. 10._wp**(-14._wp)) THEN + harmonic_perm(k)=0.0_wp + ELSE + DO kk=k,N_active-1 + harmonic_perm(k)=harmonic_perm(k)+thick(kk)/perm(kk) + END DO + !bottom layer is included in a linear fashion + harmonic_perm(k)=harmonic_perm(k)+(thick(N_active)*psi_s(N_active)/psi_s_min)/perm(N_active) + harmonic_perm(k)=(SUM(thick(k:N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min)/harmonic_perm(k) + END IF + END DO + END IF + + !___________Raleigh number is calculated (see Notz 2005)__________________________ + DO k=1,N_active-1 + d_S_br = S_br(k)-S_br(N_active) + height = SUM(thick((k+1):N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min !Height adjustment + IF (harmonic_flag == 1) THEN + ray(k) = grav*rho_l*bbeta*d_S_br*height* MINVAL(perm(k:N_active)) + ELSE IF (harmonic_flag == 2) THEN + ray(k) = grav*rho_l*bbeta*d_S_br*height* harmonic_perm(k) + END IF + ray(k) = ray(k)/(kappa_l*mu) + ray(k) = MAX(ray(k),0._wp) + END DO + + + !##########Where the magic happens################################# + DO k = N_active-1,1,-1 + IF (ray(k) > ray_mini ) THEN + S_abs(k) = S_abs(k)*0.99 + END IF + END DO + grav_drain = 0._wp + END SUBROUTINE fl_grav_drain_simple + +END MODULE mo_grav_drain + diff --git a/mo_grotz.f90 b/mo_grotz.f90 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..76242c4c1b879ae8749ca7c7a730d4e1a2a4e439 100644 --- a/mo_grotz.f90 +++ b/mo_grotz.f90 @@ -0,0 +1,803 @@ +!> +!! The most important module of SAMSIM +!! +!! The module mo_grotz contains the most important subroutine grotz (Named after GRiewank nOTZ). +!! Mo_grotz is called by SAMSIM.f90. +!! Subroutine grotz contains the time loop, as well as the initialization, and calls all other branches of the model. +!! This model was developed from scratch by Philipp Griewank during and after his PhD at Max Planck Institute of Meteorology from 2010-2014. +!! The code is intended to be understandable and most subroutines, modules, functions, parameters, and global variables have doxygen compatible descriptions. +!! In addition to the doxygen generated description, some python plotscripts are available to plot model output. +!! +!! +!! +!! +!! @author Philipp Griewank +!! +!! +!! COPYRIGHT +!! +!! This file is part of SAMSIM. +!! +!! SAMSIM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! +!! SAMSIM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License along with SAMSIM. If not, see <http://www.gnu.org/licenses/>. +! +!! +!! +!! +!! @par Revision History +!! Started by Philipp Griewank 2012-08-28 +!! + +MODULE mo_grotz + +CONTAINS + + !> + !! + !! Main subroutine of SAMSIM, a 1D thermodynamic seaice model. + !! A semi-adaptive grid is used which is managed by mo_layer_dynamics. + !! + !! The basic rundown of the time loop is: + !! 1. Calculate the current ice/snow state and forcing, as well as gravity drainage and flooding + !! 2. Apply all the fluxes, recalculate ice state + !! 3. Flushing and layer dynamics + !! + !! Here is the full rundown of what happens in mo_grotz: + !! + !! - Initialization: all fields are initialized for the given config.json file, and the output is formatted + !! - Input and Forcing read in. + !! TIME LOOP BEGINS: + !! - Calculate the total ice properties, total freshwater, thermal resistivity, energy, bulk salinity + !! - Determine snow and rain rates + !! - Calculate snow thermodynamics + !! - Calculate inner ice thermodynamic fluxes + !! - Calculate brine flux from expulsion + !! - Raw output written out if debug_flag is set to 2 + !! - Standard output written + !! - Flooding parametrized + !! - Lowest layer mixing with underlying water + !! - Gravity drainage parametrized + !! - Calcuating and applying the heat fluxes + !! - After heatfluxes are applied new liquidus thermal equilibrium is calculated + !! - Flushing is parametrized + !! - Chemistry advection calculated + !! - Layer Dynamics + !! TIME LOOP ENDS + !! -Final output, files closed, and fields deallocated + !! + !! + !! IMPORTANT: + !! To get the correct freshwater amount make sure the freshwater is calculated using a salinity value to compare against. + !! + !! + !! Common errors leading to termination are: too small timestep, bad programming + !! + !! @par Revision History + !! 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) + SUBROUTINE grotz () + + USE mo_parameters + USE mo_thermo_functions + USE mo_data + USE mo_init + USE mo_layer_dynamics + USE mo_mass + USE mo_grav_drain + USE mo_output + USE mo_flush + USE mo_flood + USE mo_snow + USE mo_functions + USE mo_heat_fluxes + + IMPLICIT NONE + + CHARACTER*12000 :: description !< String to describes simulation which is output into dat_settings + + !Bastard variables, are used for various dirty deeds + INTEGER :: jj + REAL(wp) :: temp, temp2, temp4, temp5, temp_2017_H, temp_2017_m !when a real is needed !< Niels, 2017 add: temp_2017_m + temp5 = 0. + + + + + !########################################################################################## + !Initialization + !########################################################################################## + CALL init() + + ! Get description + OPEN(unit = 99, file = 'Run_specifics/description.txt') + READ(99, *) description + CLOSE(99) + + CALL output_begin(Nlayer, debug_flag, format_T, format_psi, format_thick, format_snow, format_T2m_top, & + & format_perm, format_melt) + IF (bgc_flag==2) THEN + CALL output_begin_bgc(Nlayer, N_bgc, format_bgc) + END IF + CALL output_settings(description, N_top, N_bottom, Nlayer, fl_q_bottom, T_bottom, S_bu_bottom, thick_0, time_out, & + & time_total, dt, boundflux_flag, albedo_flag, grav_flag, flush_flag, flood_flag, & + & grav_heat_flag, flush_heat_flag, & + harmonic_flag, prescribe_flag, salt_flag, turb_flag, bottom_flag, tank_flag, precip_flag, bgc_flag, & + & N_bgc, k_snow_flush) + + + !########################################################################################## + !Time Loop + !########################################################################################## + DO i = 1, i_time + + ! IF (n_time_out == i_time_out .OR. i==1) THEN + ! WRITE(*,*), '0: Tsnow:', T_snow, ' so_prec:', solid_precip, ' T_Top:', T_top + ! END IF + + !########################################################################################## + !Vital signs. Calculates stored energy, freshwater column, thermal resistivity, ice thickness, and bulk salinity + !########################################################################################## + !stored energy + energy_stored = H_abs_snow + SUM(H_abs(1:N_active)) - T_bottom * SUM(m(1:N_active)) * c_l + + !freshwater + freshwater = SUM(m(1:N_active)) / rho_l + freshwater = freshwater * (1._wp - SUM(S_abs(1:N_active)) / SUM(m(1:N_active)) / ref_salinity) + freshwater = freshwater + m_snow / rho_l + + !Total resistivity only includes a fraction of the lowest layer. + total_resist = 0._wp + DO jj = 1, N_active - 1 + total_resist = total_resist + thick(jj) / (psi_l(jj) * k_l + psi_s(jj) * k_s) + END DO + total_resist = total_resist + thick(N_active) * psi_s(N_active) / psi_s_min * & + & (psi_s_min * k_s + 1._wp - psi_s_min * k_l) + IF (thick_snow>thick_min / 110._wp) THEN + total_resist = total_resist + thick_snow / func_k_snow(m_snow, thick_snow) + END IF + + !Thickness is the thickness of all layers plus a fraction of the lowest layer + IF (N_active>1) THEN + thickness = SUM(thick(1:N_active - 1)) + ELSE + thickness = 0.0 + END IF + thickness = thickness + thick(N_active) * psi_s(N_active) / psi_s_min + + !Bulk salinity of the ice + IF (N_active>1) THEN + bulk_salin = SUM(S_abs(1:N_active - 1)) + S_abs(N_active) * psi_s(N_active) / psi_s_min + bulk_salin = bulk_salin / (SUM(m(1:N_active - 1)) + m(N_active) * psi_s(N_active) / psi_s_min) + ELSE + bulk_salin = S_abs(1) / m(1) + END IF + + + + !########################################################################################## + !Selection and linear interpolation of boudary values + !########################################################################################## + + IF (time>time_input(time_counter)) THEN ! + time_counter = time_counter + 1 + END IF + + ! check if timestep for boundary values equals timestep of simulation + ! If so, use boundary value at same timestep + IF (abs(time - time_input(time_counter)) < tolerance) THEN + solid_precip = precip_s_input(time_counter) + liquid_precip = precip_l_input(time_counter) + T2m = T2m_input(time_counter) + T_bottom = T_bottom_input(time_counter) + fl_q_bottom = fl_q_bottom_input(time_counter) + IF (boundflux_flag == 1) THEN + T_top = T_top_input(time_counter) + END IF + ! in tank experiments, bulk salinity at bottom varies since total amount of salt must be conserved + IF (tank_flag .NE. 2) THEN + S_bu_bottom = S_bu_bottom_input(time_counter) + END IF + + ! If timestep is different, interpolate linearly between boundary values + ELSE + temp = (time - time_input(time_counter - 1)) / (time_input(time_counter) - time_input(time_counter - 1)) + solid_precip = (1._wp - temp) * precip_s_input(time_counter - 1) + temp * precip_s_input(time_counter) + liquid_precip = (1._wp - temp) * precip_l_input(time_counter - 1) + temp * precip_l_input(time_counter) + T2m = (1._wp - temp) * T2m_input(time_counter - 1) + temp * T2m_input(time_counter) + T_bottom = (1._wp - temp) * T_bottom_input(time_counter - 1) + temp * T_bottom_input(time_counter) + fl_q_bottom = (1._wp - temp) * fl_q_bottom_input(time_counter - 1) + temp * fl_q_bottom_input(time_counter) + IF (boundflux_flag == 1) THEN + T_top = (1._wp - temp) * T_top_input(time_counter - 1) + temp * T_top_input(time_counter) + END IF + IF (tank_flag .NE. 2) THEN + S_bu_bottom = (1._wp - temp) * S_bu_bottom_input(time_counter - 1) + temp * S_bu_bottom_input(time_counter) + END IF + END IF + + + !########################################################################################## + !Snow fall + !########################################################################################## + IF (precip_flag==1) THEN + IF (MAX(liquid_precip, solid_precip)>0.0_wp .AND. N_Active>1) THEN !Precip on/in snow layer + 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) + END IF + + ELSE IF (precip_flag==0) THEN + IF (MAX(liquid_precip, solid_precip)>0.0_wp .AND. N_Active>1) THEN !Precip on/in snowlayer + CALL snow_precip (m_snow, H_abs_snow, thick_snow, dt, liquid_precip, T2m, solid_precip) + 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) + END IF + END IF + + + !########################################################################################## + !Snow thermodynamics and volume adjustment + !########################################################################################## + + !< Niels, 2017 add: adjusted to snow_flush_flag==1 + IF (thick_snow>0.0_wp) THEN + IF (snow_flush_flag == 0) THEN + CALL snow_thermo (psi_l_snow, psi_s_snow, psi_g_snow, thick_snow, S_abs_snow, & + & H_abs_snow, m_snow, T_snow, m(1), thick(1), H_abs(1)) + melt_thick_snow = 0._wp + ELSE IF (snow_flush_flag == 1) THEN + melt_thick_snow = 0._wp + CALL snow_thermo_meltwater (psi_l_snow, psi_s_snow, psi_g_snow, thick_snow, S_abs_snow, & + & H_abs_snow, m_snow, T_snow, m(1), thick(1), H_abs(1), melt_thick_snow) + END IF + ELSE + thick_snow = 0.0_wp + m_snow = 0.0_wp + psi_s_snow = 0.0_wp + psi_l_snow = 0.0_wp + psi_g_snow = 0.0_wp + H_abs_snow = 0.0_wp + S_abs_snow = 0.0_wp + melt_thick_snow = 0.0_wp !< Niels, 2017 add: melt_thick_snow + END IF + + + !########################################################################################## + !Inner layer thermodynamics and Expulsion + !########################################################################################## + T_test = T_bottom + DO k = N_active, 1, -1 + S_bu(k) = S_abs(k) / m(k) + H(k) = H_abs(k) / m(k) + CALL getT(H(k), S_bu(k), T_test, T(k), phi(k), k) + T_test = T(k) + S_br(k) = func_S_br(T(k), S_bu(k)) + !Expulsion + CALL Expulsion(phi(k), thick(k), m(k), psi_s(k), psi_l(k), psi_g(k), V_ex(k)) + END DO + + !########################################################################################## + !Brine flux due to Expulsion + !########################################################################################## + CALL expulsion_flux (thick, V_ex, Nlayer, N_active, psi_g, fl_m, m) + 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) + IF (bgc_flag ==2) THEN + DO k = 1, N_active + fl_brine_bgc(k, k + 1) = -fl_m(k + 1) + END DO + END IF + END IF + + + !########################################################################################## + !Raw Output, if debug flag is set to 2 then all layer data is output + !each timestep + !########################################################################################## + IF (debug_flag==2) THEN + CALL output_raw(Nlayer, N_active, time, T, thick, S_bu, psi_s, psi_l, psi_g) + CALL output_raw_snow(time, T_snow, thick_snow, S_abs_snow, m_snow, psi_s_snow, psi_l_snow, psi_g_snow) + END IF + + DO k = N_active, 1, -1 + S_bu(k) = S_abs(k) / m(k) + END DO + + !########################################################################################## + !Standard output at each time_out + !########################################################################################## + IF (n_time_out == i_time_out .OR. i==1) THEN + + !Various things are calculated or recalculated before they are output + IF (N_active>1) THEN + freeboard = func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag) + ELSE + freeboard = 0._wp + END IF + + IF (grav_flag == 2) then + IF (abs(grav_drain) < tolerance) THEN + grav_temp = 0._wp + ELSE + grav_temp = grav_temp / grav_drain + END IF + grav_salt = grav_salt / time_out + grav_drain = grav_drain / time_out + end if + + !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) + + IF (bgc_flag==2) THEN + CALL output_bgc(Nlayer, N_active, bgc_bottom, N_bgc, bgc_abs, psi_l, thick, m, format_bgc) + END IF + + + !Printed to console + WRITE(*, '(A10,I2,A15,F6.3,A12,F5.3,A14,F7.3,A30,F3.1,A14,F7.4,A10,F7.3,A7,F8.3)')& + 'progress: ', INT(100._wp * (time + dt) / time_total), & + '%, thickness: ', thickness, & + 'm, albedo: ', func_albedo(thick_snow, T_snow, psi_l(1), thick_min, albedo_flag), & + ', surface T: ', T_top, & + 'C, thermal stability (<0.5): ', k_s * dt / rho_s / c_s / MIN(thick(1), thick_0)**2._wp, & + ', snow_thick:', thick_snow, & + !', snow gt zero: ', (thick_snow>0._wp),& !< uncomment if you wish to know if there is a very thin snow layer on the ice + !', Flush?:', flush_question,& !< tells you if flushing occurs in this iteration + ', 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!' + + ! Write other boundary data to console + WRITE(*,*)& + 'fl_q_bottom: ', fl_q_bottom, & + 'T_bottom: ', T_bottom, & + 'S_bu_bottom: ', s_bu_bottom, & + 'T_top: ', T_top, & + 'fl_sw: ', fl_sw, & + 'fl_rest; ', fl_rest + + + + + !Resetting gravity drainage things and time_out counter + grav_drain = 0.0_wp + grav_salt = 0.0_wp + grav_temp = 0.0_wp + + melt_thick_output(:) = 0._wp + + n_time_out = 0 + ELSE + n_time_out = n_time_out + 1 + END IF + + + + !########################################################################################## + !When the lowest layer contains gas, it is replaced with ocean water + !########################################################################################## + IF (psi_g(N_active)>0.0_wp) THEN + temp2 = psi_g(N_active) * thick(N_active) * rho_l + m(N_active) = m(N_active) + temp2 + S_abs(N_active) = S_abs(N_active) + temp2 * S_bu_bottom + H_abs(N_active) = H_abs(N_active) + temp2 * c_l * T_bottom + END IF + + + !########################################################################################## + !Snow top layer coupling for thin snow + !When snow is thinner than thick_min it is considered to be in thermal + !equilibrium with the top ice layer + !########################################################################################## + IF (m_snow>0.0_wp .AND. thick_snow<thick_min) THEN + CALL snow_coupling (H_abs_snow, phi_s, T_snow, H_abs(1), H(1), phi(1), T(1), m_snow, S_abs_snow, m(1), & + & S_bu(1)) + END IF + + + + + !########################################################################################## + !Flooding + !########################################################################################## + IF (N_active>1 .AND. flood_flag>1) THEN + freeboard = func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag) + IF (freeboard<0.0_wp) THEN + IF (flood_flag==2) THEN + IF (bgc_flag == 2) THEN + CALL flood (freeboard, psi_s, psi_l, S_abs, H_abs, m, T, thick, dt, Nlayer, N_active, & + &T_bottom, S_bu_bottom, H_abs_snow, m_snow, thick_snow, psi_g_snow, debug_flag,& + & fl_brine_bgc) + ELSE + CALL flood (freeboard, psi_s, psi_l, S_abs, H_abs, m, T, thick, dt, Nlayer, N_active, & + &T_bottom, S_bu_bottom, H_abs_snow, m_snow, thick_snow, psi_g_snow, debug_flag) + + END IF + ELSE IF (flood_flag==3 .AND. freeboard<neg_free) THEN + CALL flood_simple (freeboard, S_abs, H_abs, m, thick, T_bottom, S_bu_bottom, H_abs_snow, & + &m_snow, thick_snow, psi_g_snow, Nlayer, N_active, debug_flag) + END IF + END IF + END IF + + !########################################################################################## + !Turbulence, mixes lowest layer with underlying water + !########################################################################################## + IF (turb_flag==2) THEN + IF (bgc_flag == 2) then + Call sub_turb_flux(T_bottom, S_bu_bottom, T(N_active), S_abs(N_active), m(N_active), dt, & + & N_bgc, bgc_bottom, bgc_abs(N_active, :)) + ELSE + Call sub_turb_flux(T_bottom, S_bu_bottom, T(N_active), S_abs(N_active), m(N_active), dt, N_bgc) + end if + END IF + + + !########################################################################################## + !Gravity drainage: grav_flag 2 complex, 3 simple + !########################################################################################## + IF (grav_flag==2 .AND. N_active>1) THEN + 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) + 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) + END IF + + ELSE IF (grav_flag==3 .AND. N_active>1) THEN + CALL fl_grav_drain_simple (psi_s, psi_l, thick, S_abs, S_br, Nlayer, N_active, ray, & + grav_drain, harmonic_flag) + end if + + !########################################################################################## + !Prescribes salinity profile for prescribe_flag = =2 + !########################################################################################## + IF (prescribe_flag==2) THEN + + !Linear advance over the lowest 15 cm from S_bu_bottom to 4 ppt + k = N_active + DO WHILE (k>1 .AND. SUM(thick(k:N_active))<0.15_wp) + S_bu(k) = S_bu_bottom - SUM(thick(k:N_active)) / 0.15_wp * (S_bu_bottom - 4._wp) + k = k - 1 + END DO + DO WHILE (k>1 .AND. SUM(thick(k:N_active))>=0.15_wp) + S_bu(k) = 4._wp - 4._wp * (SUM(thick(k:N_active)) - 0.15_wp) / (SUM(thick(1:N_active)) - 0.15_wp) + k = k - 1 + S_bu(1) = 0._wp + END DO + S_bu(N_active) = S_bu_bottom + S_abs = S_bu * m + END IF + + + + !########################################################################################## + !Calculating water concentrations if tank_flag == 2 + !########################################################################################## + IF(tank_flag==2) then + S_bu_bottom = (S_total - SUM(S_abs(:))) / (m_total - SUM(m)) + IF(bgc_flag==2) then + bgc_bottom = (bgc_total(1) - SUM(bgc_abs(:, 1))) / (m_total - SUM(m)) + end if + end if + + !########################################################################################## + !Calculating and applying heatfluxes + !########################################################################################## + + call sub_heat_fluxes() + + + + !########################################################################################## + !T and phi are recalculated before applying layer dynamics and flushing + !########################################################################################## + + T_test = T_bottom + DO k = N_active, 1, -1 + S_bu(k) = S_abs(k) / m(k) + H(k) = H_abs(k) / m(k) + CALL getT(H(k), S_bu(k), T_test, T(k), phi(k), k) + T_test = T(k) !The temperature of the layer i is used to initialize the iteration for layer i-1 + END DO + + temp_2017_H = H_abs(1) + H_abs_snow + melt_thick_snow * rho_l * c_l * T_snow !< Niels, 2017 add: checking energy and mass conservation later + temp_2017_m = m(1) + m_snow + melt_thick_snow * rho_l !< Niels, 2017 + + melt_thick_snow_old = melt_thick_snow !< Niels, 2017 add: keep meltwater during recalculation of thermodynamics + IF (thick_snow>0.0_wp) THEN + IF (snow_flush_flag == 0) THEN + CALL snow_thermo (psi_l_snow, psi_s_snow, psi_g_snow, thick_snow, S_abs_snow, H_abs_snow, m_snow, & + & T_snow, m(1), thick(1), H_abs(1)) + melt_thick_snow = 0._wp + ELSE IF (snow_flush_flag == 1) THEN + melt_thick_snow = 0._wp + CALL snow_thermo_meltwater (psi_l_snow, psi_s_snow, psi_g_snow, thick_snow, S_abs_snow, H_abs_snow,& + & m_snow, T_snow, m(1), thick(1), H_abs(1), melt_thick_snow) + END IF + + ELSE + thick_snow = 0.0_wp + m_snow = 0.0_wp + psi_s_snow = 0.0_wp + psi_l_snow = 0.0_wp + psi_g_snow = 0.0_wp + H_abs_snow = 0.0_wp + S_abs_snow = 0.0_wp + melt_thick_snow = 0.0_wp !< Niels, 2017 + END IF + melt_thick_snow = melt_thick_snow_old + melt_thick_snow !< Niels, 2017 add: keep meltwater during recalculation of thermodynamics + + + !########################################################################################## + !Flushing preparations + !########################################################################################## + !Subroutines are called to calculate the melt thickness for flush_flag greater 3,4,5,6 for boundflux_flag 2 and 3 + IF (N_active>1 .AND. flush_flag>2) THEN + IF (boundflux_flag ==2) THEN + T_freeze = func_T_freeze(S_abs(1) / m(1), salt_flag) + melt_thick = 0.0_wp + IF (func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag)> & + & 0.0000000000001_wp) THEN + IF (psi_s(1)<psi_s_top_min .OR. T_top>=T_freeze) THEN + + CALL sub_melt_thick(psi_l(1), psi_s(1), psi_g(1), T(1), T_freeze, T_top, fl_Q(1), & + & thick_snow, dt, melt_thick, thick(1), thick_min) + IF (thick_snow>=thick_min / 100._wp .AND. melt_thick>0.00000000001_wp .AND. & + & abs(melt_thick_snow) < tolerance) THEN + !< Niels, 2017 add: .and. melt_thick_snow == 0. + CALL sub_melt_snow(melt_thick, thick(1), thick_snow, H_abs(1), H_abs_snow, m(1), & + & m_snow, psi_g_snow) + END IF + END IF + END IF + END IF + + IF (boundflux_flag==3) THEN + T_freeze = func_T_freeze(S_abs(1) / m(1), salt_flag) + melt_thick = 0.0_wp + !< Niels, 2017 add: freeboard_snow_flag, import parameter for lab experiments, see in definition + IF (func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag)> & + & 0.0000000000001_wp) THEN + IF (psi_s(1)<psi_s_top_min .OR. T2m>=T_freeze) THEN + CALL sub_melt_thick(psi_l(1), psi_s(1), psi_g(1), T(1), T_freeze, T2m, fl_Q(1), thick_snow,& + & dt, melt_thick, thick(1), & + &thick_min) + melt_thick = max(melt_thick, 0._wp) !< Niels, 2017 add: just for the case + IF (thick_snow>=thick_min / 100._wp .AND. melt_thick>0.00000000001_wp .AND. & + & abs(melt_thick_snow) < tolerance) THEN + CALL sub_melt_snow(melt_thick, thick(1), thick_snow, H_abs(1), H_abs_snow, m(1), & + & m_snow, psi_g_snow) + END IF + END IF + END IF + END IF + END IF + + + !########################################################################################## + !Flushing + !########################################################################################## + freeboard = func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag) + + melt_thick_output(1) = melt_thick_output(1) + melt_thick !< Niels, 2017 add: for later evaluation: accumulated melt_thick + melt_thick_output(2) = melt_thick_output(2) + melt_thick_snow !< Niels, 2017 add: accumulated melt_thick_snow + + melt_thick = melt_thick + melt_thick_snow !< Niels, 2017 add: Ice melt plus snow melt + + IF (melt_thick_snow > 0._wp) THEN + !< Niels, 2017 add: Add snow melt water to upper ice layer to keep balances + H_abs(1) = H_abs(1) + melt_thick_snow * rho_l * c_l * T_snow + S_abs(1) = S_abs(1) + melt_thick_snow * rho_l * func_S_br(T_snow, S_abs_snow / m_snow) + thick(1) = thick(1) + melt_thick_snow + m(1) = m(1) + melt_thick_snow * rho_l + S_bu(1) = S_abs(1) / m(1) + H(1) = H_abs(1) / m(1) + END IF + + !< Niels, 2017 add: energy and mass balance test + IF (temp_2017_H - (H_abs(1) + H_abs_snow) > 0.00000001_wp) THEN + print*, i, ' enthalpy problem during snow melt:', temp_2017_H, H_abs(1) + H_abs_snow + ELSE IF (temp_2017_m - (m(1) + m_snow) > 0.00000001_wp) THEN + print*, i, 'mass problem during snow melt:', temp_2017_m, m(1) + m_snow + END IF + + ! Start flushing routines + + flush_v_old(:) = flush_v(:) !< Niels, 2017 add: used for accumulated flush output + flush_h_old(:) = flush_h(:) + + flush_v(:) = 0._wp + flush_h(:) = 0._wp + + IF (N_active>1 .and. freeboard>0.001_wp) THEN + IF (flush_flag==4) THEN + IF (melt_thick>0.000000000001_wp .AND. N_active>2) THEN + H_abs(1) = H_abs(1) - melt_thick * rho_l * c_l * T(1) + S_abs(1) = S_abs(1) * (1._wp - (melt_thick * rho_l) / m(1)) + thick(1) = thick(1) - melt_thick + m(1) = m(1) - melt_thick * rho_l + IF (S_abs(1)<0.0_wp) THEN + PRINT*, 'sorry bro, but you got a problem', melt_thick, psi_l * thick(1) + END IF + END IF + + ELSE IF (flush_flag==5) THEN + IF (melt_thick>0.000000000001_wp .AND. N_active>2 .AND. freeboard>0.0_wp) THEN + freeboard = func_freeboard(N_active, Nlayer, psi_s, psi_g, m, thick, m_snow, freeboard_snow_flag) + flush_question = 'yes' !< Niels, 2017 add: used for stdout + IF (bgc_flag == 2) THEN + CALL flush3(freeboard, psi_l, thick, thick_0, S_abs, H_abs, m, T, dt, Nlayer, N_active, & + &T_bottom, S_bu_bottom, melt_thick, debug_flag, flush_heat_flag, melt_err, perm, & + & flush_v, flush_h, psi_g, rho_l, snow_flush_flag, fl_brine_bgc) !< Niels, 2017 add: increased number of arguments + ELSE + CALL flush3(freeboard, psi_l, thick, thick_0, S_abs, H_abs, m, T, dt, Nlayer, N_active, & + &T_bottom, S_bu_bottom, melt_thick, debug_flag, flush_heat_flag, melt_err, perm, & + & flush_v, flush_h, psi_g, rho_l, snow_flush_flag) !< Niels, 2017 add: increased number of arguments + END IF + END IF + ELSE IF (flush_flag==6) THEN + IF (melt_thick>0.000000000001_wp .AND. N_active>2 .AND. thick_snow<thick_0) THEN + CALL flush4 (psi_l, thick, T, S_abs, H_abs, m, Nlayer, N_active, melt_thick, debug_flag) + END IF + END IF + END IF + + flush_v(:) = flush_v(:) + flush_v_old(:) !< Niels, 2017 add: for output + flush_h(:) = flush_h(:) + flush_h_old(:) + + !########################################################################################## + !BGC advection + !########################################################################################## + IF (bgc_flag == 2) THEN + + CALL bgc_advection (Nlayer, N_active, N_bgc, fl_brine_bgc, bgc_abs, psi_l, thick, bgc_bottom) + fl_brine_bgc = 0._wp + + END IF + + temp4 = sum(H_abs(:)) + + !########################################################################################## + !Layer Dynamics + !########################################################################################## + IF (N_active>1) THEN + !########################################################################################## + !Bottom & top growth or melt + !########################################################################################## + IF (phi(N_active)> psi_s_min .OR. phi(N_active - 1)<=psi_s_min / 2._wp .OR. thick(1) / thick_0>1.5_wp & + &.OR. thick(1) / thick_0<0.5_wp) THEN + IF (bgc_flag == 2) THEN + CALL layer_dynamics(phi, N_active, Nlayer, N_bottom, N_middle, N_top, m, S_abs, H_abs, thick, & + & thick_0, T_bottom, S_bu_bottom, bottom_flag, debug_flag, melt_thick_output(3), & + & N_bgc, bgc_abs, bgc_bottom) + ELSE + CALL layer_dynamics(phi, N_active, Nlayer, N_bottom, N_middle, N_top, m, S_abs, H_abs, thick, & + & thick_0, T_bottom, S_bu_bottom, bottom_flag, debug_flag, melt_thick_output(3), N_bgc) + END IF + END IF + + + !If the lowest layer was deactivated, its remaining values are scrubbed. + IF (N_active<Nlayer .AND. abs(thick(MIN(N_active + 1, Nlayer))) < tolerance) THEN + T(N_active + 1) = T_bottom + S_bu(N_active + 1) = S_bu_bottom + H(N_active + 1) = 0.0_wp + psi_l(N_active + 1) = 1.0_wp + psi_s(N_active + 1) = 0.0_wp + IF (bgc_flag==2) THEN + bgc_abs(N_active + 1, :) = 0._wp + END IF + + END IF + + ELSE !If Ice appears in the only active surface layer + IF (phi(1).GT.psi_s_min) THEN + IF (bgc_flag == 2) THEN + CALL layer_dynamics(phi, N_active, Nlayer, N_bottom, N_middle, N_top, m, S_abs, H_abs, thick, & + & thick_0, T_bottom, S_bu_bottom, bottom_flag, debug_flag, melt_thick_output(3), & + & N_bgc, bgc_abs, bgc_bottom) + ELSE + CALL layer_dynamics(phi, N_active, Nlayer, N_bottom, N_middle, N_top, m, S_abs, H_abs, thick, & + & thick_0, T_bottom, S_bu_bottom, bottom_flag, debug_flag, melt_thick_output(3), N_bgc) + END IF + END IF + END IF + + temp5 = temp5 + (sum(H_abs(:)) - temp4)!/dt-fl_q(N_active+1)+fl_q(1) + + !########################################################################################## + !timestep + !########################################################################################## + time = time + dt + + + !########################################################################################## + !Health check + !########################################################################################## + IF (MINVAL(psi_s(1:N_active))<0.0_wp) THEN + PRINT*, 'negative solid fraction, aborted' + PRINT*, 'wtfrho', m(1), thick(1) + STOP 1337 + ELSE IF(MINVAL(S_abs(1:N_active))<0.0_wp) THEN + PRINT*, 'negative salinity, aborted' + PRINT*, 'wtfsalt', S_abs + FORALL (k = 1:N_active) + S_abs(k) = MAX(S_abs(k), 0._wp) + END FORALL + !STOP 1337 + END IF + + + + + + + + + + !########################################################################################## + !END OF Time LOOP + !########################################################################################## + + END DO + + !########################################################################################## + !Final Output + !########################################################################################## + WRITE(*, *)'Run completed, total ice thickness at end of run:', SUM(thick(1:N_active - 1)), & + & ' melt_err= ', melt_err + + + + !Closing output files + IF (debug_flag==2) THEN + DO k = 1, Nlayer + CLOSE(k + 100) + END DO + END IF + IF (bgc_flag==2) THEN + DO k = 1, N_bgc + CLOSE(2 * k + 400) + CLOSE(2 * k + 401) + END DO + END IF + CLOSE(30) + CLOSE(31) + CLOSE(32) + CLOSE(33) + CLOSE(34) + CLOSE(35) + CLOSE(40) + CLOSE(41) + CLOSE(42) + CLOSE(43) + CLOSE(44) + CLOSE(45) + CLOSE(46) !< Niels, 2017 add:46-50 for some output files + CLOSE(47) + CLOSE(48) + CLOSE(49) + CLOSE(50) + CLOSE(66) + + !Deallocating arrays + CALL sub_deallocate + END SUBROUTINE grotz + +END MODULE mo_grotz