diff --git a/components/mpas-seaice/bld/build-namelist b/components/mpas-seaice/bld/build-namelist index b23329e53f6e..a79e2ee7ee06 100755 --- a/components/mpas-seaice/bld/build-namelist +++ b/components/mpas-seaice/bld/build-namelist @@ -567,6 +567,7 @@ add_default($nl, 'config_forcing_cycle_start'); add_default($nl, 'config_forcing_cycle_duration'); add_default($nl, 'config_forcing_precipitation_units'); add_default($nl, 'config_forcing_sst_type'); +add_default($nl, 'config_forcing_bgc_type'); add_default($nl, 'config_update_ocean_fluxes'); add_default($nl, 'config_include_pond_freshwater_feedback'); diff --git a/components/mpas-seaice/bld/build-namelist-section b/components/mpas-seaice/bld/build-namelist-section index 9a3251fbdfe8..5256c85c6c2d 100644 --- a/components/mpas-seaice/bld/build-namelist-section +++ b/components/mpas-seaice/bld/build-namelist-section @@ -109,6 +109,7 @@ add_default($nl, 'config_forcing_cycle_start'); add_default($nl, 'config_forcing_cycle_duration'); add_default($nl, 'config_forcing_precipitation_units'); add_default($nl, 'config_forcing_sst_type'); +add_default($nl, 'config_forcing_bgc_type'); add_default($nl, 'config_update_ocean_fluxes'); add_default($nl, 'config_include_pond_freshwater_feedback'); diff --git a/components/mpas-seaice/bld/namelist_files/namelist_defaults_mpassi.xml b/components/mpas-seaice/bld/namelist_files/namelist_defaults_mpassi.xml index f1c04cef9324..83a228341776 100644 --- a/components/mpas-seaice/bld/namelist_files/namelist_defaults_mpassi.xml +++ b/components/mpas-seaice/bld/namelist_files/namelist_defaults_mpassi.xml @@ -107,6 +107,7 @@ '2-00-00_00:00:00' 'mm_per_sec' 'ncar' +'ISPOL' true false diff --git a/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml b/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml index f1240794b100..4c87838a0a38 100644 --- a/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml +++ b/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml @@ -319,7 +319,7 @@ Default: Defined in namelist_defaults.xml category="initialize" group="initialize"> Initial condition type for sea ice state. -Valid values: 'restart', 'uniform', 'circle', 'square' or 'uniform_interior' +Valid values: 'restart', 'uniform', 'circle', 'square', 'uniform_interior' or 'uniform_1D' Default: Defined in namelist_defaults.xml @@ -461,7 +461,7 @@ Default: Defined in namelist_defaults.xml category="forcing" group="forcing"> Atmospheric forcing type. -Valid values: 'CORE' +Valid values: 'CORE' or 'ISPOL' Default: Defined in namelist_defaults.xml @@ -501,7 +501,15 @@ Default: Defined in namelist_defaults.xml category="forcing" group="forcing"> Sea surface temperature ocean forcing type. -Valid values: 'ncar' +Valid values: 'ncar' or 'ISPOL' +Default: Defined in namelist_defaults.xml + + + +Ocean biogeochemistry forcing type. + +Valid values: 'ISPOL' Default: Defined in namelist_defaults.xml diff --git a/components/mpas-seaice/cime_config/buildnml b/components/mpas-seaice/cime_config/buildnml index b7429ff23f1f..e6622283c953 100755 --- a/components/mpas-seaice/cime_config/buildnml +++ b/components/mpas-seaice/cime_config/buildnml @@ -518,6 +518,13 @@ def buildnml(case, caseroot, compname): lines.append(' filename_interval="none"') lines.append(' input_interval="initial_only" />') lines.append('') + lines.append('') + lines.append('') lines.append('') diff --git a/components/mpas-seaice/src/Registry.xml b/components/mpas-seaice/src/Registry.xml index 8a48e739e27f..a74ca1688574 100644 --- a/components/mpas-seaice/src/Registry.xml +++ b/components/mpas-seaice/src/Registry.xml @@ -450,7 +450,7 @@ /> + + + + + + + + + + + + + + + + + + + + + + + + + + +▶ + + + diff --git a/components/mpas-seaice/src/shared/mpas_seaice_constants.F b/components/mpas-seaice/src/shared/mpas_seaice_constants.F index 252015c31c20..d96edeff891b 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_constants.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_constants.F @@ -103,5 +103,21 @@ module seaice_constants skeletalLayerThickness = sk_l , & gramsCarbonPerMolCarbon = R_gC2molC ! g carbon per mol carbon + ! ocean biogeochemistry ISPOL values + real(kind=RKIND), parameter, public :: & + oceanAmmoniumISPOL = 1.0_RKIND, & ! mmol N m-3 + oceanDMSISPOL = 0.1_RKIND, & ! mmol S m-3 + oceanDMSPISPOL = 0.1_RKIND, & ! mmol S m-3 + oceanDiatomsISPOL = 1.0_RKIND, & ! mmol N m-3 + oceanSmallAlgaeISPOL = 0.0057_RKIND, & ! mmol N m-3 + oceanPhaeocystisISPOL = 0.0027_RKIND, & ! mmol N m-3 + oceanPolysaccharidsISPOL = 16.2_RKIND, & ! mmol C m-3 + oceanLipidsISPOL = 9.0_RKIND, & ! mmol C m-3 + oceanProteinsCarbonISPOL = 9.0_RKIND, & ! mmol C m-3 + oceanDICISPOL = 1.0_RKIND, & ! mmol C m-3 + oceanProteinsISPOL = 12.9_RKIND, & ! mmol N m-3 + oceanDissolvedIronISPOL = 0.4_RKIND, & ! mmol Fe m-3 + oceanParticulateIronISPOL = 2.0_RKIND,& ! mmol Fe m-3 + oceanHumicsISPOL = 1.0_RKIND ! mmol C m-3 end module seaice_constants diff --git a/components/mpas-seaice/src/shared/mpas_seaice_forcing.F b/components/mpas-seaice/src/shared/mpas_seaice_forcing.F index c3dc2453f394..e5c82692d75c 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_forcing.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_forcing.F @@ -73,10 +73,12 @@ subroutine seaice_forcing_init(domain, clock) logical, pointer :: & config_use_forcing, & - config_use_data_icebergs + config_use_data_icebergs, & + config_use_column_biogeochemistry call MPAS_pool_get_config(domain % configs, "config_use_forcing", config_use_forcing) call MPAS_pool_get_config(domain % configs, "config_use_data_icebergs", config_use_data_icebergs) + call MPAS_pool_get_config(domain % configs, "config_use_column_biogeochemistry", config_use_column_biogeochemistry) if (config_use_forcing) then @@ -86,6 +88,10 @@ subroutine seaice_forcing_init(domain, clock) ! init the ocean forcing call init_oceanic_forcing(domain) + ! init the ocean bgc forcing + if (config_use_column_biogeochemistry) & + call init_oceanic_bgc_forcing(domain) + endif ! init the data iceberg forcing @@ -120,6 +126,8 @@ subroutine init_atmospheric_forcing(domain) select case (trim(config_atmospheric_forcing_type)) case ("CORE") call init_atmospheric_forcing_CORE(domain) + case ("ISPOL") + call init_atmospheric_forcing_ISPOL(domain) case ("hourly_velocities") call init_atmospheric_forcing_hourly_velocities(domain) case default @@ -288,6 +296,183 @@ subroutine init_atmospheric_forcing_CORE(domain) end subroutine init_atmospheric_forcing_CORE +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! init_atmospheric_forcing_ISPOL +! +!> \brief Initialize the forcing objects +!> \author Nicole Jeffery, LANL +!> \date Nov, 2022 +!> \details +!> This routine calls the MPAS_forcing module subroutines that initializes +!> a 1D forcing dataset. The data corresponds to the Weddell Sea (67.9S, 54W) +!> from June 16, 2004 to December 31, 2004, and is from the 2004 +!> Ice Station Polarstern (ISPOL) field experiment. +! +!----------------------------------------------------------------------- + + subroutine init_atmospheric_forcing_ISPOL(domain) + + type(domain_type) :: domain + + real(kind=RKIND), pointer :: & + config_dt + + character(len=strKIND), pointer :: & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration + + logical, pointer :: & + config_do_restart + + character(len=strKIND) :: & + forcingIntervalSixHourly, & + forcingReferenceTimeSixHourly, & + forcingIntervalDaily, & + forcingReferenceTimeDaily + + ! get atmospheric forcing configuration options + call MPAS_pool_get_config(domain % configs, "config_forcing_start_time", config_forcing_start_time) + call MPAS_pool_get_config(domain % configs, "config_dt", config_dt) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_start", config_forcing_cycle_start) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_duration", config_forcing_cycle_duration) + call MPAS_pool_get_config(domain % configs, "config_do_restart", config_do_restart) + + ! create the six hourly forcing group + call MPAS_forcing_init_group(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + domain, & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration, & + config_do_restart, & + .false.) + + forcingIntervalSixHourly = "06:00:00" + forcingReferenceTimeSixHourly = "2004-01-01_00:00:00" + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + "shortwaveDown", & + "ISPOLSixHourlyForcing", & + "atmos_forcing", & + "shortwaveDown", & + "linear", & + forcingReferenceTimeSixHourly, & + forcingIntervalSixHourly, & + "next") + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + "longwaveDown", & + "ISPOLSixHourlyForcing", & + "atmos_coupling", & + "longwaveDown", & + "linear", & + forcingReferenceTimeSixHourly, & + forcingIntervalSixHourly, & + "next") + + call MPAS_forcing_init_field_data(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + domain % streamManager, & + config_do_restart, & + .false.) + + ! create the daily forcing group + call MPAS_forcing_init_group(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + domain, & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration, & + config_do_restart) + + forcingIntervalDaily = "00-00-01_00:00:00" + forcingReferenceTimeDaily = "2004-01-01_00:00:00" + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + "airTemperature", & + "ISPOLDailyForcing", & + "atmos_coupling", & + "airTemperature", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily, & + "next") + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + "airSpecificHumidity", & + "ISPOLDailyForcing", & + "atmos_coupling", & + "airSpecificHumidity", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily, & + "next") + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + "uAirVelocity", & + "ISPOLDailyForcing", & + "atmos_coupling", & + "uAirVelocity", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily, & + "next") + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + "vAirVelocity", & + "ISPOLDailyForcing", & + "atmos_coupling", & + "vAirVelocity", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily, & + "next") + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + "rainfallRate", & + "ISPOLDailyForcing", & + "atmos_coupling", & + "rainfallRate", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily, & + "next") + + call MPAS_forcing_init_field_data(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + domain % streamManager, & + config_do_restart, & + .false.) + + end subroutine init_atmospheric_forcing_ISPOL + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! init_atmospheric_forcing_hourly_velocities @@ -501,6 +686,20 @@ subroutine atmospheric_forcing(& streamManager, & config_dt) + else if (trim(config_atmospheric_forcing_type) == "ISPOL") then + + call MPAS_forcing_get_forcing(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + streamManager, & + config_dt) + + call MPAS_forcing_get_forcing(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_daily", & + streamManager, & + config_dt) + endif block => domain % blocklist @@ -510,6 +709,8 @@ subroutine atmospheric_forcing(& select case (trim(config_atmospheric_forcing_type)) case ("CORE") call prepare_atmospheric_coupling_variables_CORE(block) + case ("ISPOL") + call prepare_atmospheric_coupling_variables_ISPOL(block) end select ! perform post coupling operations @@ -697,112 +898,260 @@ end subroutine prepare_atmospheric_coupling_variables_CORE !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! post_atmospheric_coupling +! prepare_atmospheric_coupling_variables_ISPOL ! !> \brief -!> \author Adrian K. Turner, LANL +!> \author Nicole Jeffery, LANL !> \date !> \details -!> +!> Defines the cell fields corresponding to the 1D forcing fields. ! !----------------------------------------------------------------------- - subroutine post_atmospheric_coupling(block) + subroutine prepare_atmospheric_coupling_variables_ISPOL(block) - use seaice_mesh, only: & - seaice_latlon_vector_rotation_forward + use seaice_constants, only: & + seaiceFreshWaterFreezingPoint, & + pii type (block_type), pointer :: block type (mpas_pool_type), pointer :: & mesh, & atmosCoupling, & - atmosForcing + atmosForcing, & + oceanCoupling, & + tracers_aggregate real(kind=RKIND), dimension(:), pointer :: & + airLevelHeight, & + airPotentialTemperature, & + airTemperature, & + airSpecificHumidity, & + airDensity, & + shortwaveDown, & shortwaveVisibleDirectDown, & shortwaveVisibleDiffuseDown, & shortwaveIRDirectDown, & shortwaveIRDiffuseDown, & - uAirVelocity, & - vAirVelocity, & - windSpeed, & - shortwaveDown, & - latCell, & + rainfallRate, & + snowfallRate, & + cloudFraction, & lonCell, & - xCell, & - yCell, & - zCell + latCell, & + iceAreaCell - real(kind=RKIND), pointer :: & - sphere_radius + type (MPAS_time_type) :: & + currentForcingTime - logical, pointer :: & - config_rotate_cartesian_grid, & - on_a_sphere + real(kind=RKIND) :: & + secondsToday integer, pointer :: & nCellsSolve integer :: & + dayOfYear, & iCell call MPAS_pool_get_subpool(block % structs, "mesh", mesh) call MPAS_pool_get_subpool(block % structs, "atmos_coupling", atmosCoupling) + call MPAS_pool_get_subpool(block % structs, "ocean_coupling", oceanCoupling) + call MPAS_pool_get_subpool(block % structs, "tracers_aggregate", tracers_aggregate) call MPAS_pool_get_subpool(block % structs, "atmos_forcing", atmosForcing) - call MPAS_pool_get_config(block % configs, "config_rotate_cartesian_grid", config_rotate_cartesian_grid) - call MPAS_pool_get_config(mesh, "on_a_sphere", on_a_sphere) - - call MPAS_pool_get_config(mesh, "sphere_radius", sphere_radius) call MPAS_pool_get_dimension(mesh, "nCellsSolve", nCellsSolve) - call MPAS_pool_get_array(mesh, "latCell", latCell) call MPAS_pool_get_array(mesh, "lonCell", lonCell) - call MPAS_pool_get_array(mesh, "xCell", xCell) - call MPAS_pool_get_array(mesh, "yCell", yCell) - call MPAS_pool_get_array(mesh, "zCell", zCell) + call MPAS_pool_get_array(mesh, "latCell", latCell) + call MPAS_pool_get_array(atmosCoupling, "airLevelHeight", airLevelHeight) + call MPAS_pool_get_array(atmosCoupling, "airPotentialTemperature", airPotentialTemperature) + call MPAS_pool_get_array(atmosCoupling, "airTemperature", airTemperature) + call MPAS_pool_get_array(atmosCoupling, "airSpecificHumidity", airSpecificHumidity) + call MPAS_pool_get_array(atmosCoupling, "airDensity", airDensity) call MPAS_pool_get_array(atmosCoupling, "shortwaveVisibleDirectDown", shortwaveVisibleDirectDown) call MPAS_pool_get_array(atmosCoupling, "shortwaveVisibleDiffuseDown", shortwaveVisibleDiffuseDown) call MPAS_pool_get_array(atmosCoupling, "shortwaveIRDirectDown", shortwaveIRDirectDown) call MPAS_pool_get_array(atmosCoupling, "shortwaveIRDiffuseDown", shortwaveIRDiffuseDown) - call MPAS_pool_get_array(atmosCoupling, "uAirVelocity", uAirVelocity) - call MPAS_pool_get_array(atmosCoupling, "vAirVelocity", vAirVelocity) + call MPAS_pool_get_array(atmosCoupling, "rainfallRate", rainfallRate) + call MPAS_pool_get_array(atmosCoupling, "snowfallRate", snowfallRate) - call MPAS_pool_get_array(atmosForcing, "windSpeed", windSpeed) + call MPAS_pool_get_array(atmosForcing, "cloudFraction", cloudFraction) call MPAS_pool_get_array(atmosForcing, "shortwaveDown", shortwaveDown) - do iCell = 1, nCellsSolve - - ! wind speed - windSpeed(iCell) = sqrt(uAirVelocity(iCell)**2 + vAirVelocity(iCell)**2) - - enddo ! iCell + call MPAS_pool_get_array(tracers_aggregate, "iceAreaCell", iceAreaCell) - if (on_a_sphere .and. config_rotate_cartesian_grid) then + ! get the current time + call MPAS_forcing_get_forcing_time(& + seaiceForcingGroups, & + "seaice_atmospheric_forcing_sixhrly", & + currentForcingTime) - do iCell = 1, nCellsSolve + ! get the number of seconds so far today + call get_seconds_today(& + currentForcingTime, & + secondsToday, & + dayOfYear) - ! rotate velocities from geographical to local - call seaice_latlon_vector_rotation_forward(& - uAirVelocity(iCell), & - vAirVelocity(iCell), & - uAirVelocity(iCell), & - vAirVelocity(iCell), & - latCell(iCell), & - lonCell(iCell), & - xCell(iCell), & - yCell(iCell), & - zCell(iCell), & - sphere_radius, & - config_rotate_cartesian_grid) + do iCell = 1, nCellsSolve - enddo ! iCell + ! not sure if we need this + cloudFraction(iCell) = 0.5_RKIND - endif ! + ! limit air temperature values where ice is present + if (iceAreaCell(iCell) > 0.1_RKIND) then + airTemperature(iCell) = min(airTemperature(iCell), seaiceFreshWaterFreezingPoint + 0.1_RKIND) + endif - do iCell = 1, nCellsSolve + ! prevent supersaturated humidity + call limit_specific_humidity(& + airTemperature(iCell), & + airSpecificHumidity(iCell)) + + shortwaveVisibleDirectDown(iCell) = shortwaveDown(iCell) * fracShortwaveVisibleDirect + shortwaveVisibleDiffuseDown(iCell) = shortwaveDown(iCell) * fracShortwaveVisibleDiffuse + shortwaveIRDirectDown(iCell) = shortwaveDown(iCell) * fracShortwaveIRDirectDown + shortwaveIRDiffuseDown(iCell) = shortwaveDown(iCell) * fracShortwaveIRDiffuseDown + + ! ensure physically realistic values + shortwaveDown(iCell) = max(shortwaveDown(iCell),0.0_RKIND) + rainfallRate(iCell) = max(rainfallRate(iCell),0.0_RKIND) + airSpecificHumidity(iCell) = max(airSpecificHumidity(iCell),0.0_RKIND) + + ! atmospheric level height + airLevelHeight(iCell) = 10.0_RKIND ! for velocities, 2-m for air T and specific humidity + + ! air potential temperature + airPotentialTemperature(iCell) = airTemperature(iCell) + + ! air density + airDensity(iCell) = 1.3_RKIND + + ! precipitation + + ! divide total precipitation between rain and snow + snowfallRate(iCell) = 0.0_RKIND + + if (airTemperature(iCell) < seaiceFreshWaterFreezingPoint) then + + snowfallRate(iCell) = rainfallRate(iCell) + rainfallRate(iCell) = 0.0_RKIND + + endif + + enddo ! iCell + + end subroutine prepare_atmospheric_coupling_variables_ISPOL + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! post_atmospheric_coupling +! +!> \brief +!> \author Adrian K. Turner, LANL +!> \date +!> \details +!> +! +!----------------------------------------------------------------------- + + subroutine post_atmospheric_coupling(block) + + use seaice_mesh, only: & + seaice_latlon_vector_rotation_forward + + type (block_type), pointer :: block + + type (mpas_pool_type), pointer :: & + mesh, & + atmosCoupling, & + atmosForcing + + real(kind=RKIND), dimension(:), pointer :: & + shortwaveVisibleDirectDown, & + shortwaveVisibleDiffuseDown, & + shortwaveIRDirectDown, & + shortwaveIRDiffuseDown, & + uAirVelocity, & + vAirVelocity, & + windSpeed, & + shortwaveDown, & + latCell, & + lonCell, & + xCell, & + yCell, & + zCell + + real(kind=RKIND), pointer :: & + sphere_radius + + logical, pointer :: & + config_rotate_cartesian_grid, & + on_a_sphere + + integer, pointer :: & + nCellsSolve + + integer :: & + iCell + + call MPAS_pool_get_subpool(block % structs, "mesh", mesh) + call MPAS_pool_get_subpool(block % structs, "atmos_coupling", atmosCoupling) + call MPAS_pool_get_subpool(block % structs, "atmos_forcing", atmosForcing) + + call MPAS_pool_get_config(block % configs, "config_rotate_cartesian_grid", config_rotate_cartesian_grid) + call MPAS_pool_get_config(mesh, "on_a_sphere", on_a_sphere) + + call MPAS_pool_get_config(mesh, "sphere_radius", sphere_radius) + call MPAS_pool_get_dimension(mesh, "nCellsSolve", nCellsSolve) + + call MPAS_pool_get_array(mesh, "latCell", latCell) + call MPAS_pool_get_array(mesh, "lonCell", lonCell) + call MPAS_pool_get_array(mesh, "xCell", xCell) + call MPAS_pool_get_array(mesh, "yCell", yCell) + call MPAS_pool_get_array(mesh, "zCell", zCell) + + call MPAS_pool_get_array(atmosCoupling, "shortwaveVisibleDirectDown", shortwaveVisibleDirectDown) + call MPAS_pool_get_array(atmosCoupling, "shortwaveVisibleDiffuseDown", shortwaveVisibleDiffuseDown) + call MPAS_pool_get_array(atmosCoupling, "shortwaveIRDirectDown", shortwaveIRDirectDown) + call MPAS_pool_get_array(atmosCoupling, "shortwaveIRDiffuseDown", shortwaveIRDiffuseDown) + call MPAS_pool_get_array(atmosCoupling, "uAirVelocity", uAirVelocity) + call MPAS_pool_get_array(atmosCoupling, "vAirVelocity", vAirVelocity) + + call MPAS_pool_get_array(atmosForcing, "windSpeed", windSpeed) + call MPAS_pool_get_array(atmosForcing, "shortwaveDown", shortwaveDown) + + do iCell = 1, nCellsSolve + + ! wind speed + windSpeed(iCell) = sqrt(uAirVelocity(iCell)**2 + vAirVelocity(iCell)**2) + + enddo ! iCell + + if (on_a_sphere .and. config_rotate_cartesian_grid) then + + do iCell = 1, nCellsSolve + + ! rotate velocities from geographical to local + call seaice_latlon_vector_rotation_forward(& + uAirVelocity(iCell), & + vAirVelocity(iCell), & + uAirVelocity(iCell), & + vAirVelocity(iCell), & + latCell(iCell), & + lonCell(iCell), & + xCell(iCell), & + yCell(iCell), & + zCell(iCell), & + sphere_radius, & + config_rotate_cartesian_grid) + + enddo ! iCell + + endif ! + + do iCell = 1, nCellsSolve ! shortwave - comment out for bit for bit tests shortwaveDown(iCell) = & @@ -1286,6 +1635,8 @@ subroutine init_oceanic_forcing(domain) select case (trim(config_forcing_sst_type)) case ("ncar") call init_oceanic_forcing_ncar(domain) + case ("ISPOL") + call init_oceanic_forcing_ISPOL(domain) case ("hourly_velocities") call init_oceanic_forcing_hourly_velocities(domain) case ("som_varying_thickness") @@ -1296,6 +1647,34 @@ subroutine init_oceanic_forcing(domain) end subroutine init_oceanic_forcing +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! init_ocean_bgc_forcing +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date 3rd March 2023 +!> \details +!> Uses an ocean biogeochemistry forcing file +! for bgc enabled sea ice only test cases. +!----------------------------------------------------------------------- + + subroutine init_oceanic_bgc_forcing(domain) + + type (domain_type) :: domain + + character(len=strKIND), pointer :: & + config_forcing_bgc_type + + call MPAS_pool_get_config(domain % configs, "config_forcing_bgc_type", config_forcing_bgc_type) + + select case (trim(config_forcing_bgc_type)) + case ("ISPOL") + call init_oceanic_bgc_forcing_ISPOL(domain) + end select + + end subroutine init_oceanic_bgc_forcing + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! init_oceanic_forcing_ncar @@ -1465,6 +1844,262 @@ subroutine init_oceanic_forcing_ncar(domain) end subroutine init_oceanic_forcing_ncar +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! init_oceanic_forcing_ISPOL +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date October 2022 +!> \details +!> The physical data corresponds to the Weddell Sea (67.9S, 54W) location from a +!> POP 1-degree (gx1v3) field derived from an ocean-ice simulation +!> forced with Ice Station Polarstern (ISPOL) 2004 atmospheric data. +!> Biogeochemical data is World Ocean Atlas (NCEP 2013) surface nutrient data. +! +!----------------------------------------------------------------------- + + subroutine init_oceanic_forcing_ISPOL(domain) + + type (domain_type) :: domain + + logical, pointer :: & + config_do_restart + + character(len=strKIND), pointer :: & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration + + character(len=strKIND) :: & + forcingIntervalMonthly, & + forcingReferenceTimeMonthly, & + forcingIntervaldaily, & + forcingReferenceTimedaily + + ! get oceanic forcing configuration options + call MPAS_pool_get_config(domain % configs, "config_do_restart", config_do_restart) + + ! get ispol forcing configuration options + call MPAS_pool_get_config(domain % configs, "config_forcing_start_time", config_forcing_start_time) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_start", config_forcing_cycle_start) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_duration", config_forcing_cycle_duration) + + ! create the sea surface temperature forcing group + call MPAS_forcing_init_group(& + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + domain, & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration, & + config_do_restart) + + forcingIntervalMonthly = "00-01-00_00:00:00" + forcingReferenceTimeMonthly = "0001-01-15_00:00:00" + + ! sea surface temperature + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "seaSurfaceTemperature", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "seaSurfaceTemperature", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! sea surface salinity + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "seaSurfaceSalinity", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "seaSurfaceSalinity", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! u ocean velocity + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "uOceanVelocity", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "uOceanVelocity", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! v ocean velocity + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "vOceanVelocity", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "vOceanVelocity", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! ocean mixed layer depth + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "oceanMixedLayerDepth", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "oceanMixedLayerDepth", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! ocean sea surface tilt + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "seaSurfaceTiltU", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "seaSurfaceTiltU", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "seaSurfaceTiltV", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "seaSurfaceTiltV", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + ! deep ocean heat flux convergence + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + "oceanHeatFluxConvergence", & + "ISPOLMonthlyForcing", & + "ocean_coupling", & + "oceanHeatFluxConvergence", & + "linear", & + forcingReferenceTimeMonthly, & + forcingIntervalMonthly) + + call MPAS_forcing_init_field_data(& + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + domain % streamManager, & + config_do_restart, & + .false.) + + end subroutine init_oceanic_forcing_ISPOL + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! init_oceanic_forcing_ISPOL +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date March 2023 +!> \details +!> surface nutrient data corresponds to the Weddell Sea (67.9S, 54W) location from +!> the World Ocean Atlas (NCEP 2013). +! +!----------------------------------------------------------------------- + + subroutine init_oceanic_bgc_forcing_ISPOL(domain) + + type (domain_type) :: domain + + logical, pointer :: & + config_do_restart + + character(len=strKIND), pointer :: & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration + + character(len=strKIND) :: & + forcingIntervalMonthly, & + forcingReferenceTimeMonthly, & + forcingIntervaldaily, & + forcingReferenceTimedaily + + ! get oceanic forcing configuration options + call MPAS_pool_get_config(domain % configs, "config_do_restart", config_do_restart) + + ! get ispol forcing configuration options + call MPAS_pool_get_config(domain % configs, "config_forcing_start_time", config_forcing_start_time) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_start", config_forcing_cycle_start) + call MPAS_pool_get_config(domain % configs, "config_forcing_cycle_duration", config_forcing_cycle_duration) + + + ! ocean surface biogeochemistry fields + + ! create the sea surface biogeochemistry forcing group + call MPAS_forcing_init_group(& + seaiceForcingGroups, & + "seaice_oceanic_biogeochemistry_forcing_daily", & + domain, & + config_forcing_start_time, & + config_forcing_cycle_start, & + config_forcing_cycle_duration, & + config_do_restart) + + forcingIntervalDaily = "00-00-01_00:00:00" + forcingReferenceTimeDaily = "2004-01-01_00:00:00" + + ! sea surface nitrate + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_biogeochemistry_forcing_daily", & + "oceanNitrateConc1D", & + "ISPOLDailyWOAForcing", & + "1D_ocean_biogeochemistry_forcing", & + "oceanNitrateConc1D", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily) + + ! sea surface silicate + call MPAS_forcing_init_field(& + domain % streamManager, & + seaiceForcingGroups, & + "seaice_oceanic_biogeochemistry_forcing_daily", & + "oceanSilicateConc1D", & + "ISPOLDailyWOAForcing", & + "1D_ocean_biogeochemistry_forcing", & + "oceanSilicateConc1D", & + "linear", & + forcingReferenceTimeDaily, & + forcingIntervalDaily) + + call MPAS_forcing_init_field_data(& + seaiceForcingGroups, & + "seaice_oceanic_biogeochemistry_forcing_daily", & + domain % streamManager, & + config_do_restart, & + .false.) + + end subroutine init_oceanic_bgc_forcing_ISPOL + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! init_oceanic_forcing_hourly_velocities @@ -1758,15 +2393,19 @@ subroutine oceanic_forcing(& config_dt character(len=strKIND), pointer :: & - config_forcing_sst_type + config_forcing_sst_type, & + config_forcing_bgc_type logical, pointer :: & - config_do_restart + config_do_restart, & + config_use_column_biogeochemistry ! configurations call mpas_pool_get_config(domain % configs, 'config_dt', config_dt) call mpas_pool_get_config(domain % configs, 'config_forcing_sst_type', config_forcing_sst_type) + call mpas_pool_get_config(domain % configs, 'config_forcing_bgc_type', config_forcing_bgc_type) call mpas_pool_get_config(domain % configs, 'config_do_restart', config_do_restart) + call mpas_pool_get_config(domain % configs, 'config_use_column_biogeochemistry', config_use_column_biogeochemistry) ! use the forcing layer to get data if (trim(config_forcing_sst_type) == 'ncar' .or. & @@ -1797,17 +2436,43 @@ subroutine oceanic_forcing(& streamManager, & config_dt) + else if (trim(config_forcing_sst_type) == 'ISPOL') then + + call MPAS_forcing_get_forcing(& + seaiceForcingGroups, & + "seaice_oceanic_forcing_monthly", & + streamManager, & + config_dt) + endif + if (config_use_column_biogeochemistry) then + if (trim(config_forcing_bgc_type) == 'ISPOL') then + + call MPAS_forcing_get_forcing(& + seaiceForcingGroups, & + "seaice_oceanic_biogeochemistry_forcing_daily", & + streamManager, & + config_dt) + endif endif block => domain % blocklist do while (associated(block)) - if (trim(config_forcing_sst_type) == 'ncar' .or. & - trim(config_forcing_sst_type) == 'som_varying_thickness') then - + ! convert the input forcing variables to the coupling variables + select case (trim(config_forcing_sst_type)) + case ("ncar","som_varying_thickness") call prepare_oceanic_coupling_variables_ncar(block, firstTimeStep) + case ("ISPOL") + call prepare_oceanic_coupling_variables_ISPOL(block, firstTimeStep) + end select - endif + if (config_use_column_biogeochemistry) then + ! convert the input forcing variables to the coupling variables + select case (trim(config_forcing_bgc_type)) + case ("ISPOL") + call prepare_oceanic_coupling_bgc_variables_ISPOL(block, firstTimeStep) + end select + end if call post_oceanic_coupling(block) @@ -1900,6 +2565,205 @@ subroutine prepare_oceanic_coupling_variables_ncar(block, firstTimeStep) end subroutine prepare_oceanic_coupling_variables_ncar +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! prepare_oceanic_coupling_variables_ISPOL +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date October, 2022 +!> \details +!> Prepares the ocean coupling variables from the ISPOL 1D forcing +! +!----------------------------------------------------------------------- + + subroutine prepare_oceanic_coupling_variables_ISPOL(block, firstTimeStep) + + use ice_colpkg, only: & + colpkg_sea_freezing_temperature + + type (block_type), pointer :: block + + logical, intent(in) :: & + firstTimeStep + + type (mpas_pool_type), pointer :: & + mesh, & + ocean_coupling + + real(kind=RKIND), dimension(:), pointer :: & + seaSurfaceTemperature, & + seaSurfaceSalinity, & + oceanMixedLayerDepth, & + seaFreezingTemperature + + character(len=strKIND), pointer :: & + config_sea_freezing_temperature_type + + logical, pointer :: & + config_do_restart + + integer, pointer :: & + nCellsSolve, & + nCells + + integer :: & + iCell + + call MPAS_pool_get_config(block % configs, "config_do_restart", config_do_restart) + call MPAS_pool_get_config(block % configs, "config_sea_freezing_temperature_type", config_sea_freezing_temperature_type) + + call MPAS_pool_get_subpool(block % structs, "mesh", mesh) + call MPAS_pool_get_subpool(block % structs, "ocean_coupling", ocean_coupling) + + call MPAS_pool_get_dimension(mesh, "nCellsSolve", nCellsSolve) + call MPAS_pool_get_dimension(mesh, "nCells", nCells) + + call MPAS_pool_get_array(ocean_coupling, "seaSurfaceTemperature", seaSurfaceTemperature) + call MPAS_pool_get_array(ocean_coupling, "seaSurfaceSalinity", seaSurfaceSalinity) + call MPAS_pool_get_array(ocean_coupling, "oceanMixedLayerDepth", oceanMixedLayerDepth) + call MPAS_pool_get_array(ocean_coupling, "seaFreezingTemperature", seaFreezingTemperature) + + 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) = colpkg_sea_freezing_temperature(seaSurfaceSalinity(iCell)) + + enddo ! iCell + + ! only update sea surface temperature on first non-restart timestep + if (firstTimeStep .and. .not. config_do_restart) then + + do iCell = 1, nCellsSolve + + ! sea surface temperature + seaSurfaceTemperature(iCell) = max(seaSurfaceTemperature(iCell), seaFreezingTemperature(iCell)) + + enddo ! iCell + + endif + + end subroutine prepare_oceanic_coupling_variables_ISPOL + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! prepare_oceanic_coupling_variables_ISPOL +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date March, 2023 +!> \details +!> Defines the ocean bgc cell fields corresponding to the 1D forcing fields. +! +!----------------------------------------------------------------------- + + subroutine prepare_oceanic_coupling_bgc_variables_ISPOL(block, firstTimeStep) + + use seaice_constants, only: & + oceanAmmoniumISPOL, oceanDMSISPOL, oceanDMSPISPOL, oceanDiatomsISPOL, & + oceanSmallAlgaeISPOL, oceanPhaeocystisISPOL, oceanPolysaccharidsISPOL, & + oceanLipidsISPOL, oceanDICISPOL, oceanProteinsISPOL, oceanDissolvedIronISPOL, & + oceanParticulateIronISPOL, oceanHumicsISPOL, oceanProteinsCarbonISPOL + + type (block_type), pointer :: block + + logical, intent(in) :: & + firstTimeStep + + type (mpas_pool_type), pointer :: & + mesh, & + oceanBGCForcing, & + biogeochemistry + + real(kind=RKIND), dimension(:,:), pointer :: & + oceanDONConc, & + oceanDICConc, & + oceanDOCConc, & + oceanAlgaeConc, & + oceanParticulateIronConc, & + oceanDissolvedIronConc + + real(kind=RKIND), dimension(:), pointer :: & + oceanNitrateConc, & + oceanSilicateConc, & + oceanAmmoniumConc, & + oceanDMSConc, & + oceanDMSPConc, & + oceanHumicsConc + + real(kind=RKIND), pointer :: & + oceanNitrateConc1D, & + oceanSilicateConc1D + + logical, pointer :: & + config_do_restart, & + config_use_column_biogeochemistry + + integer, pointer :: & + nCellsSolve, & + nCells + + integer :: & + iCell + + call MPAS_pool_get_config(block % configs, "config_do_restart", config_do_restart) + call MPAS_pool_get_config(block % configs, "config_use_column_biogeochemistry", config_use_column_biogeochemistry) + + call MPAS_pool_get_subpool(block % structs, "mesh", mesh) + + call MPAS_pool_get_dimension(mesh, "nCellsSolve", nCellsSolve) + call MPAS_pool_get_dimension(mesh, "nCells", nCells) + + call MPAS_pool_get_subpool(block % structs, "biogeochemistry",biogeochemistry) + call MPAS_pool_get_subpool(block % structs, "1D_ocean_biogeochemistry_forcing", & + oceanBGCForcing) + + call MPAS_pool_get_array(oceanBGCForcing, "oceanNitrateConc1D", oceanNitrateConc1D) + call MPAS_pool_get_array(oceanBGCForcing, "oceanSilicateConc1D", oceanSilicateConc1D) + call MPAS_pool_get_array(biogeochemistry, "oceanSilicateConc", oceanSilicateConc) + call MPAS_pool_get_array(biogeochemistry, "oceanNitrateConc", oceanNitrateConc) + call MPAS_pool_get_array(biogeochemistry, "oceanAlgaeConc", oceanAlgaeConc) + call MPAS_pool_get_array(biogeochemistry, "oceanAmmoniumConc", oceanAmmoniumConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDMSConc", oceanDMSConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDMSPConc", oceanDMSPConc) + call MPAS_pool_get_array(biogeochemistry, "oceanHumicsConc", oceanHumicsConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDONConc", oceanDONConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDICConc", oceanDICConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDOCConc", oceanDOCConc) + call MPAS_pool_get_array(biogeochemistry, "oceanParticulateIronConc", oceanParticulateIronConc) + call MPAS_pool_get_array(biogeochemistry, "oceanDissolvedIronConc", oceanDissolvedIronConc) + + do iCell = 1, nCellsSolve + ! define cell quantities + oceanNitrateConc(iCell) = oceanNitrateConc1D + oceanSilicateConc(iCell) = oceanSilicateConc1D + end do + + if (firstTimeStep) then + do iCell = 1, nCellsSolve + oceanAmmoniumConc(iCell) = oceanAmmoniumISPOL + oceanDMSConc(iCell) = oceanDMSISPOL + oceanDMSPConc(iCell) = oceanDMSPISPOL + oceanAlgaeConc(:,iCell) = oceanDiatomsISPOL + oceanAlgaeConc(2,iCell) = oceanSmallAlgaeISPOL + oceanAlgaeConc(3,iCell) = oceanPhaeocystisISPOL + oceanDOCConc(:,iCell) = oceanPolysaccharidsISPOL + oceanDOCConc(2,iCell) = oceanLipidsISPOL + oceanDOCConc(3,iCell) = oceanProteinsCarbonISPOL + oceanDONConc(:,iCell) = oceanProteinsISPOL + oceanDICConc(:,iCell) = oceanDICISPOL + oceanHumicsConc(iCell) = oceanHumicsISPOL + oceanParticulateIronConc(:,iCell) = oceanParticulateIronISPOL + oceanDissolvedIronConc(:,iCell) = oceanDissolvedIronISPOL + enddo ! iCell + endif + + end subroutine prepare_oceanic_coupling_bgc_variables_ISPOL + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! post_oceanic_coupling diff --git a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F index f182abbcaeca..77fb3a7f2c84 100644 --- a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F +++ b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F @@ -317,6 +317,12 @@ subroutine init_ice_state(& block, & domain % configs) + else if (trim(config_initial_condition_type) == "uniform_1D") then + + call init_ice_state_uniform_1D(& + block, & + domain % configs) + else if (trim(config_initial_condition_type) == "special") then call init_ice_state_special(& @@ -564,6 +570,170 @@ subroutine init_ice_state_uniform_ice(& end subroutine init_ice_state_uniform_ice!}}} +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! init_ice_state_uniform_1D +! +!> \brief +!> \author Nicole Jeffery, LANL +!> \date November, 2022 +!> \details +!> Initializes sea ice in the first category only with constant +!> values and latitude bounds defined in the namelist. +! +!----------------------------------------------------------------------- + + subroutine init_ice_state_uniform_1D(& + block, & + configs)!{{{ + + use seaice_constants, only: & + seaiceDegreesToRadians + + use ice_colpkg, only: & + colpkg_init_trcr, & + colpkg_enthalpy_snow + + type(block_type), intent(inout) :: & + block !< Input/Output: + + type (MPAS_pool_type), pointer, intent(in) :: & + configs !< Input: + + type (MPAS_pool_type), pointer :: & + mesh, & + tracers, & + tracers_aggregate, & + atmos_coupling, & + ocean_coupling, & + initial + + integer, pointer :: & + nCellsSolve, & + nSnowLayers, & + nIceLayers, & + nCategories + + real(kind=RKIND), dimension(:), pointer :: & + latCell, & + airTemperature, & + seaSurfaceTemperature, & + seaFreezingTemperature + + real(kind=RKIND), dimension(:,:), pointer :: & + initialSalinityProfile, & + initialMeltingTemperatureProfile + + real(kind=RKIND), dimension(:,:,:), pointer :: & + iceAreaCategory, & + iceVolumeCategory, & + snowVolumeCategory, & + surfaceTemperature, & + snowEnthalpy, & + iceEnthalpy, & + iceSalinity + + real(kind=RKIND), dimension(:), pointer :: & + iceAreaCell, & + surfaceTemperatureCell + + real(kind=RKIND), pointer :: & + config_initial_ice_area, & + config_initial_ice_volume, & + config_initial_snow_volume, & + config_initial_latitude_north, & + config_initial_latitude_south + + integer :: & + iCell, & + iSnowLayer, & + iIceLayer, & + iCategory + + call MPAS_pool_get_config(configs, "config_initial_ice_area", config_initial_ice_area) + call MPAS_pool_get_config(configs, "config_initial_ice_volume", config_initial_ice_volume) + call MPAS_pool_get_config(configs, "config_initial_snow_volume", config_initial_snow_volume) + call MPAS_pool_get_config(configs, "config_initial_latitude_north", config_initial_latitude_north) + call MPAS_pool_get_config(configs, "config_initial_latitude_south", config_initial_latitude_south) + + call MPAS_pool_get_subpool(block % structs, "mesh", mesh) + call MPAS_pool_get_subpool(block % structs, "tracers", tracers) + call MPAS_pool_get_subpool(block % structs, "tracers_aggregate", tracers_aggregate) + call MPAS_pool_get_subpool(block % structs, "ocean_coupling", ocean_coupling) + call MPAS_pool_get_subpool(block % structs, "atmos_coupling", atmos_coupling) + call MPAS_pool_get_subpool(block % structs, "initial", initial) + + call MPAS_pool_get_dimension(mesh, "nCells", nCellsSolve) + call MPAS_pool_get_dimension(mesh, "nSnowLayers", nSnowLayers) + call MPAS_pool_get_dimension(mesh, "nIceLayers", nIceLayers) + call MPAS_pool_get_dimension(mesh, "nCategories", nCategories) + + call MPAS_pool_get_array(mesh, "latCell", latCell) + + call MPAS_pool_get_array(tracers, "iceAreaCategory", iceAreaCategory, 1) + call MPAS_pool_get_array(tracers, "iceVolumeCategory", iceVolumeCategory, 1) + call MPAS_pool_get_array(tracers, "snowVolumeCategory", snowVolumeCategory, 1) + call MPAS_pool_get_array(tracers, "surfaceTemperature", surfaceTemperature, 1) + + call MPAS_pool_get_array(atmos_coupling, "airTemperature", airTemperature) + + call MPAS_pool_get_array(ocean_coupling, "seaSurfaceTemperature", seaSurfaceTemperature) + call MPAS_pool_get_array(ocean_coupling, "seaFreezingTemperature", seaFreezingTemperature) + + call MPAS_pool_get_array(tracers_aggregate, "iceAreaCell", iceAreaCell) + call MPAS_pool_get_array(tracers_aggregate, "surfaceTemperatureCell", surfaceTemperatureCell) + call MPAS_pool_get_array(tracers, "snowEnthalpy", snowEnthalpy, 1) + call MPAS_pool_get_array(tracers, "iceEnthalpy", iceEnthalpy, 1) + call MPAS_pool_get_array(tracers, "iceSalinity", iceSalinity, 1) + + call MPAS_pool_get_array(initial, "initialSalinityProfile", initialSalinityProfile) + call MPAS_pool_get_array(initial, "initialMeltingTemperatureProfile", initialMeltingTemperatureProfile) + + do iCell = 1, nCellsSolve + + iceAreaCategory(1,:,iCell) = 0.0_RKIND + iceVolumeCategory(1,:,iCell) = 0.0_RKIND + snowVolumeCategory(1,:,iCell) = 0.0_RKIND + surfaceTemperature(1,:,iCell) = seaFreezingTemperature(iCell) + iceEnthalpy(:,:,iCell) = 0.0_RKIND + + do iSnowLayer = 1, nSnowLayers + snowEnthalpy(iSnowLayer,:,iCell) = colpkg_enthalpy_snow(0.0_RKIND) + end do + + if (latCell(iCell) > config_initial_latitude_north * seaiceDegreesToRadians .or. & + latCell(iCell) < config_initial_latitude_south * seaiceDegreesToRadians) then + + ! has ice + iceAreaCategory(1,1,iCell) = config_initial_ice_area + iceVolumeCategory(1,1,iCell) = config_initial_ice_volume + snowVolumeCategory(1,1,iCell) = config_initial_snow_volume + surfaceTemperature(1,1,iCell) = -1.0_RKIND + + call colpkg_init_trcr(& + airTemperature(iCell), & + seaFreezingTemperature(iCell), & + initialSalinityProfile(:,iCell), & + initialMeltingTemperatureProfile(:,iCell), & + surfaceTemperature(1,1,iCell), & + nIceLayers, & + nSnowLayers, & + iceEnthalpy(:,1,iCell), & + snowEnthalpy(:,1,iCell)) + endif + + iceAreaCell(iCell) = sum(iceAreaCategory(1,:,iCell)) + surfaceTemperatureCell(iCell) = -20.15_RKIND + + do iCategory = 1, nCategories + do iIceLayer = 1, nIceLayers + iceSalinity(iIceLayer,iCategory,iCell) = initialSalinityProfile(iIceLayer,iCell) + enddo ! iIceLayer + enddo ! iCategory + enddo ! iCell + + end subroutine init_ice_state_uniform_1D!}}} + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! init_ice_cice_default