Merge pull request #34 from eclare108213/prescribedice1

Icepack interface access

components/mpas-seaice/driver/ice_comp_mct.F

@@ -48,6 +48,8 @@ module ice_comp_mct
    use seaice_analysis_driver
    use seaice_column, only: seaice_column_reinitialize_fluxes, & !colpkg
                             seaice_column_coupling_prep
+   use seaice_icepack, only: seaice_icepack_reinitialize_fluxes, &
+                             seaice_icepack_coupling_prep
    use seaice_constants, only: coupleAlarmID, &
                                seaiceFreshWaterFreezingPoint, &
                                seaiceDensityIce, &
@@ -69,6 +71,7 @@ module ice_comp_mct
    use seaice_error, only: seaice_check_critical_error
 
    use ice_colpkg, only: colpkg_sea_freezing_temperature
+   use icepack_intfc, only: icepack_sea_freezing_temperature
 
 !
 ! !PUBLIC MEMBER FUNCTIONS:
@@ -745,8 +748,12 @@ end subroutine xml_stream_get_attributes
     call mpas_reset_clock_alarm(domain % clock, coupleAlarmID, ierr=ierr)
 
     ! Coupling prep
-    call seaice_column_coupling_prep(domain)
-
+    call MPAS_pool_get_config(domain % configs, "config_column_physics_type", tempCharConfig)
+    if (trim(tempCharConfig) == "icepack") then
+      call seaice_icepack_coupling_prep(domain)
+    else if (trim(tempCharConfig) == "column_package") then
+      call seaice_column_coupling_prep(domain)
+    endif ! config_column_physics_type
 
 !-----------------------------------------------------------------------
 !
@@ -1053,7 +1060,9 @@ subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i)!{{{
       logical :: streamActive, debugOn
       logical, pointer :: config_write_output_on_startup
       logical, save :: first=.true.
-      character (len=StrKIND), pointer :: config_restart_timestamp_name
+      character (len=StrKIND), pointer :: &
+            config_restart_timestamp_name, &
+            config_column_physics_type
       real(kind=RKIND), pointer :: &
             dayOfNextShortwaveCalculation ! needed for CESM like coupled simulations
 
@@ -1069,14 +1078,19 @@ subroutine ice_run_mct( EClock, cdata_i, x2i_i, i2x_i)!{{{
       if (debugOn) call mpas_log_write("=== Beginning ice_run_mct ===")
 
       call mpas_pool_get_config(domain % configs, 'config_restart_timestamp_name', config_restart_timestamp_name)
+      call MPAS_pool_get_config(domain % configs, "config_column_physics_type", config_column_physics_type)
 
       ! Setup log information.
       call shr_file_getLogUnit (shrlogunit)
       call shr_file_getLogLevel(shrloglev)
       call shr_file_setLogUnit (iceLogUnit)
 
       ! reinitialize fluxes
-      call seaice_column_reinitialize_fluxes(domain)
+      if (trim(config_column_physics_type) == "icepack") then
+         call seaice_icepack_reinitialize_fluxes(domain)
+      else if (trim(config_column_physics_type) == "column_package") then
+         call seaice_column_reinitialize_fluxes(domain)
+      endif ! config_column_physics_type
 
       ! Import state from coupler
       call ice_import_mct(x2i_i, ierr)
@@ -1803,6 +1817,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       config_use_column_biogeochemistry
 
    character(len=strKIND), pointer ::   &
+      config_column_physics_type,       &
       config_thermodynamics_type,       &
       config_ocean_surface_type
 
@@ -1913,6 +1928,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
    do while(associated(block_ptr))
 
       configs => block_ptr % configs
+      call mpas_pool_get_config(configs, "config_column_physics_type", config_column_physics_type)
       call mpas_pool_get_config(configs, "config_thermodynamics_type", config_thermodynamics_type)
       call mpas_pool_get_config(configs, "config_ocean_surface_type", config_ocean_surface_type)
       call mpas_pool_get_config(configs, "config_use_aerosols", config_use_aerosols)
@@ -1988,7 +2004,11 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
         seaSurfaceTemperature(i)       = x2i_i % rAttr(index_x2i_So_t, n)
         seaSurfaceSalinity(i)          = x2i_i % rAttr(index_x2i_So_s, n)
 
-        seaFreezingTemperature(i) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(i))
+        if (trim(config_column_physics_type) == "icepack") then
+           seaFreezingTemperature(i) = icepack_sea_freezing_temperature(seaSurfaceSalinity(i))
+        else if (trim(config_column_physics_type) == "column_package") then
+           seaFreezingTemperature(i) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(i))
+        endif ! config_column_physics_type
 
         uOceanVelocity(i)              = x2i_i % rAttr(index_x2i_So_u, n)
         vOceanVelocity(i)              = x2i_i % rAttr(index_x2i_So_v, n)
@@ -2012,8 +2032,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
            ! coupling step as a freshwater and salt flux. This step is required to balance mass
            ! and heat with the ocean.
 
-           call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i), &
-                                   config_thermodynamics_type)
+           call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i))
 
            frazilMassAdjust(i) = frazilMassFlux-frazilMassFluxRev
 
@@ -2754,8 +2773,7 @@ end subroutine basal_pressure!}}}
 ! !IROUTINE: frazil_mass
 !
 ! !INTERFACE
-   subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, &
-                                 config_thermodynamics_type)
+   subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity)
 !
 ! !DESCRIPTION:
 ! Calculate frazil mass based on on the sea surface salinity, and frazil heat flux
@@ -2769,14 +2787,20 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, &
          liquidus_temperature_mush, &
          enthalpy_mush
 
+      use icepack_intfc, only:  &
+         icepack_liquidus_temperature, &
+         icepack_enthalpy_mush
+
 ! !INPUT PARAMETERS:
       real (kind=RKIND),      intent(in)   :: freezingPotential
       real (kind=RKIND),      intent(in)   :: seaSurfaceSalinity
-      character(len=strKIND), intent(in)   :: config_thermodynamics_type
 
 ! !OUTPUT PARAMETERS:
       real (kind=RKIND),      intent(out)  :: frazilMassFlux
 
+      character(len=strKIND), pointer   :: config_thermodynamics_type
+      character(len=strKIND), pointer   :: config_column_physics_type
+
 !EOP
 !BOC
 !-----------------------------------------------------------------------
@@ -2790,6 +2814,8 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, &
       qi0new,          &
       vi0new
 
+   call MPAS_pool_get_config(domain % configs, "config_thermodynamics_type", config_thermodynamics_type)             
+   call MPAS_pool_get_config(domain % configs, "config_column_physics_type", config_column_physics_type)             
 
    if (freezingPotential > 0.0_RKIND) then
 
@@ -2799,8 +2825,15 @@ subroutine frazil_mass(freezingPotential, frazilMassFlux, seaSurfaceSalinity, &
          else
              Si0new = seaSurfaceSalinity**2 / (4.0_RKIND*seaiceFrazilSalinityReduction)
          endif
-         Ti = liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity)
-         qi0new = enthalpy_mush(Ti, Si0new)
+
+         if (trim(config_column_physics_type) == "icepack") then
+            Ti = icepack_liquidus_temperature(Si0new/seaiceFrazilIcePorosity)
+            qi0new = icepack_enthalpy_mush(Ti, Si0new)
+         else if (trim(config_column_physics_type) == "column_package") then
+            Ti = liquidus_temperature_mush(Si0new/seaiceFrazilIcePorosity)
+            qi0new = enthalpy_mush(Ti, Si0new)
+         endif ! config_column_physics_type
+
       else
          qi0new = -seaiceDensityIce*seaiceLatentHeatMelting
       endif    ! ktherm

components/mpas-seaice/src/analysis_members/mpas_seaice_temperatures.F

@@ -222,6 +222,10 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{
            colpkg_ice_temperature, &
            colpkg_snow_temperature
 
+      use icepack_intfc, only: &
+           icepack_ice_temperature, &
+           icepack_snow_temperature
+
       !-----------------------------------------------------------------
       !
       ! input variables
@@ -261,6 +265,9 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{
            tracersPool, &
            temperaturesAMPool
 
+      character(len=strKIND), pointer :: &
+           config_column_physics_type
+
       real(kind=RKIND), dimension(:,:,:), pointer :: &
            iceAreaCategory, &
            snowVolumeCategory, &
@@ -287,6 +294,8 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{
       block => domain % blocklist
       do while (associated(block))
 
+         call MPAS_pool_get_config(block % configs, "config_column_physics_type", config_column_physics_type)
+
          call MPAS_pool_get_subpool(block % structs, 'temperaturesAM', temperaturesAMPool)
          call MPAS_pool_get_subpool(block % structs, 'tracers', tracersPool)
 
@@ -309,6 +318,39 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{
          snowTemperature = 0.0_RKIND
 
          ! compute temperatures
+
+         if (trim(config_column_physics_type) == "icepack") then
+
+         do iCell = 1, nCellsSolve
+            do iCategory = 1, nCategories
+
+               ! check if ice present
+               if (iceAreaCategory(1,iCategory,iCell) > 1e-11_RKIND) then
+
+                  ! ice layers
+                  do iIceLayer = 1, nIceLayers
+                     iceTemperature(iIceLayer, iCategory, iCell) = &
+                          icepack_ice_temperature(iceEnthalpy(iIceLayer,iCategory,iCell), &
+                                                  iceSalinity(iIceLayer,iCategory,iCell))
+                  enddo ! iIceLayer
+
+                  ! snow layers
+                  if (snowVolumeCategory(1,iCategory,iCell) > 1e-11_RKIND) then
+
+                     do iSnowLayer = 1, nSnowLayers
+                        snowTemperature(iSnowLayer, iCategory, iCell) = &
+                             icepack_snow_temperature(snowEnthalpy(iSnowLayer,iCategory,iCell))
+                     enddo ! iIceLayer
+
+                  endif ! snowVolumeCategory
+
+               endif ! iceAreaCategory
+
+            enddo ! iCategory
+         enddo ! iCell
+
+         else if (trim(config_column_physics_type) == "column_package") then
+
          do iCell = 1, nCellsSolve
             do iCategory = 1, nCategories
 
@@ -337,6 +379,8 @@ subroutine seaice_compute_temperatures(domain, instance, timeLevel, err)!{{{
             enddo ! iCategory
          enddo ! iCell
 
+         endif ! config_column_physics_type
+
          block => block % next
       enddo
 

components/mpas-seaice/src/shared/mpas_seaice_constants.F

@@ -50,13 +50,7 @@ module seaice_constants
        seaiceLatentHeatVaporization    = SHR_CONST_LATVAP ,&! latent heat, vaporization freshwater (J/kg)
        seaiceLatentHeatMelting         = SHR_CONST_LATICE ,&! latent heat of melting of fresh ice (J/kg)
        seaiceReferenceSalinity         = SHR_CONST_ICE_REF_SAL,&! ice reference salinity (ppt)
-       seaiceSnowPatchiness            = 0.005_RKIND      ,&! parameter for fractional snow area (m)
-
-#ifdef RASM_MODS
-       seaiceIceOceanDragCoefficient   = 0.00962_RKIND      ! ice-ocn drag coefficient for RASM as temporary measure
-#else
-       seaiceIceOceanDragCoefficient   = 0.00536_RKIND      ! ice-ocn drag coefficient
-#endif
+       seaiceSnowPatchiness            = 0.005_RKIND        ! parameter for fractional snow area (m)
 
 !         R_gC2molC = SHR_CONST_MWC    ,&! molar mass of carbon
 
@@ -92,9 +86,7 @@ module seaice_constants
        seaiceLatentHeatMelting         = seaiceLatentHeatSublimation & ! latent heat of melting of fresh ice (J/kg)
                                        - seaiceLatentHeatVaporization, &
        seaiceReferenceSalinity         = 4._RKIND      ,&! ice reference salinity (ppt)
-       seaiceSnowPatchiness            = 0.02_RKIND    ,&! parameter for fractional snow area (m)
-
-       seaiceIceOceanDragCoefficient   = 0.00536_RKIND   ! ice-ocn drag coefficient
+       seaiceSnowPatchiness            = 0.02_RKIND      ! parameter for fractional snow area (m)
 
 !         R_gC2molC = 12.0107_RKIND, & ! molar mass of carbon
 #endif
@@ -144,7 +136,8 @@ module seaice_constants
 
   ! dynamics constants
   real(kind=RKIND), public :: &
-       seaiceIceStrengthConstantHiblerP   = 2.75e4_RKIND  ,&! P* constant in Hibler strength formula
+       seaiceIceOceanDragCoefficient      = 0.00536_RKIND, &! ice-ocn drag coefficient
+       seaiceIceStrengthConstantHiblerP   = 2.75e4_RKIND , &! P* constant in Hibler strength formula
        seaiceIceStrengthConstantHiblerC   = 20._RKIND       ! C* constant in Hibler strength formula
 
   ! minimum sea ice area

components/mpas-seaice/src/shared/mpas_seaice_forcing.F

@@ -2495,8 +2495,8 @@ end subroutine oceanic_forcing
 
   subroutine prepare_oceanic_coupling_variables_ncar(block, firstTimeStep)
 
-    use seaice_icepack, only: &
-         seaice_icepack_sea_freezing_temperature
+    use icepack_intfc, only: &
+         icepack_sea_freezing_temperature
     use seaice_column, only: &
          seaice_column_sea_freezing_temperature
 
@@ -2552,7 +2552,7 @@ subroutine prepare_oceanic_coupling_variables_ncar(block, firstTimeStep)
           oceanMixedLayerDepth(iCell) = max(oceanMixedLayerDepth(iCell), 0.0_RKIND)
 
           ! sea freezing temperature
-          seaFreezingTemperature(iCell) = seaice_icepack_sea_freezing_temperature(seaSurfaceSalinity(iCell))
+          seaFreezingTemperature(iCell) = icepack_sea_freezing_temperature(seaSurfaceSalinity(iCell))
 
        enddo ! iCell
     else if (trim(config_column_physics_type) == "column_package") then