Merge pull request #37 from eclare108213/physics_cleanup

Physics cleanup: pi, emissivity, code deprecation, more icepack calls
eclare108213 · Sep 30, 2023 · 50b3913 · 50b3913
2 parents 6b7682d + 35980d6
commit 50b3913
Show file tree

Hide file tree

Showing 8 changed files with 481 additions and 461 deletions.
diff --git a/components/mpas-seaice/src/Makefile.icepack b/components/mpas-seaice/src/Makefile.icepack
@@ -13,7 +13,6 @@ OBJS =  icepack_aerosol.o \
 	icepack_itd.o \
 	icepack_kinds.o \
 	icepack_mechred.o \
-	icepack_meltpond_cesm.o \
 	icepack_meltpond_lvl.o \
 	icepack_meltpond_topo.o \
 	icepack_mushy_physics.o \
@@ -23,7 +22,6 @@ OBJS =  icepack_aerosol.o \
 	icepack_shortwave.o \
         icepack_shortwave_data.o \
 	icepack_snow.o \
-	icepack_therm_0layer.o \
 	icepack_therm_bl99.o \
 	icepack_therm_itd.o \
 	icepack_therm_mushy.o \
@@ -33,8 +31,7 @@ OBJS =  icepack_aerosol.o \
 	icepack_warnings.o \
 	icepack_wavefracspec.o \
 	icepack_zbgc.o \
-	icepack_zbgc_shared.o \
-	icepack_zsalinity.o
+	icepack_zbgc_shared.o
 
 
 all: $(OBJS)
@@ -110,9 +107,6 @@ icepack_mechred.o:icepack_kinds.o
 icepack_mechred.o:icepack_parameters.o
 icepack_mechred.o:icepack_tracers.o
 icepack_mechred.o:icepack_warnings.o
-icepack_meltpond_cesm.o:icepack_kinds.o
-icepack_meltpond_cesm.o:icepack_parameters.o
-icepack_meltpond_cesm.o:icepack_warnings.o
 icepack_meltpond_lvl.o:icepack_kinds.o
 icepack_meltpond_lvl.o:icepack_parameters.o
 icepack_meltpond_lvl.o:icepack_therm_shared.o
@@ -146,10 +140,6 @@ icepack_snow.o:icepack_kinds.o
 icepack_snow.o:icepack_parameters.o
 icepack_snow.o:icepack_therm_shared.o
 icepack_snow.o:icepack_warnings.o
-icepack_therm_0layer.o:icepack_kinds.o
-icepack_therm_0layer.o:icepack_parameters.o
-icepack_therm_0layer.o:icepack_therm_bl99.o
-icepack_therm_0layer.o:icepack_warnings.o
 icepack_therm_bl99.o:icepack_kinds.o
 icepack_therm_bl99.o:icepack_parameters.o
 icepack_therm_bl99.o:icepack_therm_shared.o
@@ -181,13 +171,11 @@ icepack_therm_vertical.o:icepack_firstyear.o
 icepack_therm_vertical.o:icepack_flux.o
 icepack_therm_vertical.o:icepack_isotope.o
 icepack_therm_vertical.o:icepack_kinds.o
-icepack_therm_vertical.o:icepack_meltpond_cesm.o
 icepack_therm_vertical.o:icepack_meltpond_lvl.o
 icepack_therm_vertical.o:icepack_meltpond_topo.o
 icepack_therm_vertical.o:icepack_mushy_physics.o
 icepack_therm_vertical.o:icepack_parameters.o
 icepack_therm_vertical.o:icepack_snow.o
-icepack_therm_vertical.o:icepack_therm_0layer.o
 icepack_therm_vertical.o:icepack_therm_bl99.o
 icepack_therm_vertical.o:icepack_therm_mushy.o
 icepack_therm_vertical.o:icepack_therm_shared.o
@@ -211,19 +199,10 @@ icepack_zbgc.o:icepack_therm_shared.o
 icepack_zbgc.o:icepack_tracers.o
 icepack_zbgc.o:icepack_warnings.o
 icepack_zbgc.o:icepack_zbgc_shared.o
-icepack_zbgc.o:icepack_zsalinity.o
 icepack_zbgc_shared.o:icepack_kinds.o
 icepack_zbgc_shared.o:icepack_parameters.o
 icepack_zbgc_shared.o:icepack_tracers.o
 icepack_zbgc_shared.o:icepack_warnings.o
-icepack_zsalinity.o:icepack_brine.o
-icepack_zsalinity.o:icepack_kinds.o
-icepack_zsalinity.o:icepack_parameters.o
-icepack_zsalinity.o:icepack_therm_shared.o
-icepack_zsalinity.o:icepack_tracers.o
-icepack_zsalinity.o:icepack_warnings.o
-icepack_zsalinity.o:icepack_zbgc_shared.o
-
 
 .F90.o:
 	$(RM) $@ $*.mod

diff --git a/components/mpas-seaice/src/analysis_members/Makefile b/components/mpas-seaice/src/analysis_members/Makefile
@@ -21,10 +21,10 @@ FW = ../../../mpas-framework/src
 	$(RM) $@ $*.mod
 ifeq "$(GEN_F90)" "true"
 	$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
-	$(FC) $(FFLAGS) -c $*.f90 -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../column -I../shared -I../model_forward $(FCINCLUDES)
+	$(FC) $(FFLAGS) -c $*.f90 -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../icepack/columnphysics -I../column -I../shared -I../model_forward $(FCINCLUDES)
 else
-	$(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../column -I../shared -I../model_forward $(CPPINCLUDES) $(FCINCLUDES)
+	$(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../icepack/columnphysics -I../column -I../shared -I../model_forward $(CPPINCLUDES) $(FCINCLUDES)
 endif
 
 .c.o:
-	$(CC) $(CPPFLAGS) $(CFLAGS) $(CINCLUDES) -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../column -I../shared -I../model_forward -c $<
+	$(CC) $(CPPFLAGS) $(CFLAGS) $(CINCLUDES) -I$(FW)/external/esmf_time_f90 -I$(FW)/framework -I$(FW)/operators -I../icepack/columnphysics -I../column -I../shared -I../model_forward -c $<
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_column.F b/components/mpas-seaice/src/shared/mpas_seaice_column.F
@@ -9929,7 +9929,7 @@ subroutine check_column_package_configs(domain)
          config_calc_surface_temperature, &
          config_use_form_drag, &
          config_use_level_ice, &
-         config_use_cesm_meltponds, &
+         config_use_cesm_meltponds, &  ! deprecated
          config_use_level_meltponds, &
          config_use_topo_meltponds, &
          config_use_vertical_zsalinity, &
@@ -9989,7 +9989,7 @@ subroutine check_column_package_configs(domain)
     call MPAS_pool_get_config(domain % configs, "config_snow_to_ice_transition_depth", config_snow_to_ice_transition_depth)
     call MPAS_pool_get_config(domain % configs, "config_use_form_drag", config_use_form_drag)
     call MPAS_pool_get_config(domain % configs, "config_use_level_ice", config_use_level_ice)
-    call MPAS_pool_get_config(domain % configs, "config_use_cesm_meltponds", config_use_cesm_meltponds)
+    call MPAS_pool_get_config(domain % configs, "config_use_cesm_meltponds", config_use_cesm_meltponds) ! deprecated
     call MPAS_pool_get_config(domain % configs, "config_use_level_meltponds", config_use_level_meltponds)
     call MPAS_pool_get_config(domain % configs, "config_use_topo_meltponds", config_use_topo_meltponds)
     call MPAS_pool_get_config(domain % configs, "config_ocean_heat_transfer_type", config_ocean_heat_transfer_type)
@@ -10026,10 +10026,15 @@ subroutine check_column_package_configs(domain)
     !-----------------------------------------------------------------------
 
     ! check config_thermodynamics_type value
-    if (.not. (trim(config_thermodynamics_type) == "zero layer" .or. &
-               trim(config_thermodynamics_type) == "BL99" .or. &
+    if (trim(config_thermodynamics_type) == "zero layer") then
+       call mpas_log_write(&
+            "check_column_package_configs: config_thermodynamics_type) = zero layer but 0-layer thermo is being deprecated", &
+            messageType=MPAS_LOG_WARN)
+    endif
+
+    if (.not. (trim(config_thermodynamics_type) == "BL99" .or. &
                trim(config_thermodynamics_type) == "mushy")) then
-       call config_error("config_thermodynamics_type", config_thermodynamics_type, "'zero layer', 'BL99' or 'mushy'")
+       call config_error("config_thermodynamics_type", config_thermodynamics_type, "'BL99' or 'mushy'")
     endif
 
     ! check config_heat_conductivity_type value
@@ -10111,6 +10116,13 @@ subroutine check_column_package_configs(domain)
             'Check for inconsistencies in restart file: config_nIceLayers /= nIceLayers', &
             messageType=MPAS_LOG_CRIT)
 
+    ! deprecate cesm ponds
+    if (config_use_cesm_meltponds) then
+       call mpas_log_write(&
+            "check_column_package_configs: config_use_cesm_meltponds = .true. but cesm ponds are being deprecated", &
+            messageType=MPAS_LOG_CRIT)
+    endif
+
     !-----------------------------------------------------------------------
     ! Check combinations
     !-----------------------------------------------------------------------
@@ -10256,6 +10268,13 @@ subroutine check_column_package_configs(domain)
             messageType=MPAS_LOG_CRIT)
     endif
 
+    ! deprecate zsalinity
+    if (config_use_vertical_zsalinity) then
+       call mpas_log_write(&
+            "check_column_package_configs: vertical zSalinity has been deprecated", &
+            messageType=MPAS_LOG_CRIT)
+    endif
+
     ! check biogeochemistry flags:
     if (.not. config_use_column_biogeochemistry .and. (config_use_brine .or. config_use_vertical_zsalinity .or. &
         config_use_vertical_biochemistry .or. config_use_shortwave_bioabsorption .or. config_use_vertical_tracers .or. &

diff --git a/components/mpas-seaice/src/shared/mpas_seaice_constants.F b/components/mpas-seaice/src/shared/mpas_seaice_constants.F
@@ -3,18 +3,21 @@
 !  seaice_constants
 !
 !> \brief
-!> \author Adrian K. Turner, LANL
-!> \date 2013-2014
+!> \author Adrian K. Turner and Elizabeth Hunke, LANL
+!> \date 2013-2014, 2023
 !> \details
 !> Values in this module have been copied from the two ice_constants_colpkg.F90 
-!  modules in the column physics package (colpkg, in column/).  They will remain 
-!  duplicates until the column physic package is deprecated.
+!> modules in the column physics package (in column/).  They will remain 
+!> duplicates until the column physics package is deprecated.
+!
 !-----------------------------------------------------------------------
 
 module seaice_constants
 
   use mpas_derived_types
+#ifdef CCSMCOUPLED
   use shr_const_mod
+#endif
 
   private
   save
@@ -57,7 +60,7 @@ module seaice_constants
 #else
       real (kind=RKIND), parameter, public :: &
   ! fundamental constants
-       pi                              = 3.14159265358979323846_RKIND
+       pi                              = 3.14159265358979323846_RKIND, &
        seaiceSecondsPerDay             = 24.0_RKIND * 3600.0_RKIND
 
       real (kind=RKIND), public :: &
@@ -83,8 +86,9 @@ module seaice_constants
        seaiceSeaWaterSpecificHeat      = 4218._RKIND   ,&! specific heat of ocn    (J/kg/K)
                                                          ! freshwater value needed for enthalpy
        seaiceLatentHeatVaporization    = 2.501e6_RKIND ,&! latent heat, vaporization freshwater (J/kg)
-       seaiceLatentHeatMelting         = seaiceLatentHeatSublimation & ! latent heat of melting of fresh ice (J/kg)
-                                       - seaiceLatentHeatVaporization, &
+!       seaiceLatentHeatMelting         = seaiceLatentHeatSublimation & ! latent heat of melting of fresh ice (J/kg)
+!                                       - seaiceLatentHeatVaporization, &
+       seaiceLatentHeatMelting         = 3.34e5_RKIND  ,&! latent heat of melting of fresh ice (J/kg)
        seaiceReferenceSalinity         = 4._RKIND      ,&! ice reference salinity (ppt)
        seaiceSnowPatchiness            = 0.02_RKIND      ! parameter for fractional snow area (m)
 
@@ -93,23 +97,17 @@ module seaice_constants
 
   ! fundamental constants
   real (kind=RKIND), parameter, public :: &
-       pii = 3.141592653589793_RKIND, & ! echmod - why is this different from pi above?
-       seaiceDegreesToRadians = pii / 180.0_RKIND, &
-       seaiceRadiansToDegrees = 180.0_RKIND / pii, &
-!echmod - BFB but slower???
-!       pii = pi, &
-!       seaiceDegreesToRadians = pi / 180.0_RKIND, &
-!       seaiceRadiansToDegrees = 180.0_RKIND / pi, &
-       seaiceSecondsPerYear = 24.0_RKIND * 3600.0_RKIND * 365.0_RKIND, & !echmod - incorrect for leap years
-       seaiceDaysPerSecond = 1.0_RKIND/seaiceSecondsPerDay
+       pii                                = pi, &    !echmod - replace pii with seaicePi elsewhere
+       seaiceDegreesToRadians             = pi / 180.0_RKIND, &
+       seaiceRadiansToDegrees             = 180.0_RKIND / pi, &
+       seaiceSecondsPerYear               = 24.0_RKIND * 3600.0_RKIND * 365.0_RKIND, & !echmod - incorrect for leap years
+       seaiceDaysPerSecond                = 1.0_RKIND/seaiceSecondsPerDay
 
   real (kind=RKIND), public :: &
-       seaicePi
-!echmod - BFB but slower???
-!       seaicePi                           = pi ! pi
+       seaicePi                           = pi ! pi
 
   character (len=*), public, parameter :: &
-       coupleAlarmID = 'coupling'
+       coupleAlarmID                      = 'coupling'
 
   real(kind=RKIND), public :: &
        seaicePuny                         = 1.0e-11_RKIND   ! a small value

diff --git a/components/mpas-seaice/src/shared/mpas_seaice_forcing.F b/components/mpas-seaice/src/shared/mpas_seaice_forcing.F
@@ -2599,6 +2599,9 @@ subroutine prepare_oceanic_coupling_variables_ISPOL(block, firstTimeStep)
     use ice_colpkg, only: &
          colpkg_sea_freezing_temperature
 
+    use icepack_intfc, only: &
+         icepack_sea_freezing_temperature
+
     type (block_type), pointer :: block
 
     logical, intent(in) :: &
@@ -2615,7 +2618,8 @@ subroutine prepare_oceanic_coupling_variables_ISPOL(block, firstTimeStep)
          seaFreezingTemperature
 
     character(len=strKIND), pointer :: &
-         config_sea_freezing_temperature_type
+         config_sea_freezing_temperature_type, &
+         config_column_physics_type
 
     logical, pointer :: &
          config_do_restart
@@ -2641,16 +2645,29 @@ subroutine prepare_oceanic_coupling_variables_ISPOL(block, firstTimeStep)
     call MPAS_pool_get_array(ocean_coupling, "oceanMixedLayerDepth", oceanMixedLayerDepth)
     call MPAS_pool_get_array(ocean_coupling, "seaFreezingTemperature", seaFreezingTemperature)
 
-    do iCell = 1, nCellsSolve
+    if (trim(config_column_physics_type) == "icepack") then
+       do iCell = 1, nCellsSolve
+
+          ! ensure physical realism
+          seaSurfaceSalinity(iCell)   = max(seaSurfaceSalinity(iCell), 0.0_RKIND)
+          oceanMixedLayerDepth(iCell) = max(oceanMixedLayerDepth(iCell), 0.0_RKIND)
+
+          ! sea freezing temperature
+          seaFreezingTemperature(iCell) = icepack_sea_freezing_temperature(seaSurfaceSalinity(iCell))
 
-       ! ensure physical realism
-       seaSurfaceSalinity(iCell)   = max(seaSurfaceSalinity(iCell), 0.0_RKIND)
-       oceanMixedLayerDepth(iCell) = max(oceanMixedLayerDepth(iCell), 0.0_RKIND)
+       enddo ! iCell
+    else if (trim(config_column_physics_type) == "column_package") then
+       do iCell = 1, nCellsSolve
 
-       ! sea freezing temperature
-       seaFreezingTemperature(iCell) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(iCell))
+          ! ensure physical realism
+          seaSurfaceSalinity(iCell)   = max(seaSurfaceSalinity(iCell), 0.0_RKIND)
+          oceanMixedLayerDepth(iCell) = max(oceanMixedLayerDepth(iCell), 0.0_RKIND)
 
-    enddo ! iCell
+          ! sea freezing temperature
+          seaFreezingTemperature(iCell) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(iCell))
+
+       enddo ! iCell
+    endif ! config_column_physics_type
 
     ! only update sea surface temperature on first non-restart timestep
     if (firstTimeStep .and. .not. config_do_restart) then