diff --git a/README.md b/README.md index 4bd0e0ab5..84a06b056 100644 --- a/README.md +++ b/README.md @@ -26,6 +26,8 @@ where "branch name" can also be a version name ## More information +In progress documentation: https://duvivier.github.io/CICE/ + Detailed and searchable online documentation of CICE can be found [here](https://cice-consortium.github.io/CICE/). In this documentation, a [“Quick Start”](https://cice-consortium.github.io/CICE/cice_2_quick_start.html) subsection is available with instructions for running the model. A [“Testing”](https://cice-consortium.github.io/CICE/cice_7_testing.html) subsection with instructions for setting up standard tests (e.g. regression, restart) is also available. diff --git a/doc/source/.DS_Store b/doc/source/.DS_Store new file mode 100644 index 000000000..82c0810f5 Binary files /dev/null and b/doc/source/.DS_Store differ diff --git a/doc/source/cice_1_introduction.rst b/doc/source/cice_1_introduction.rst index 57200276e..1be56648a 100644 --- a/doc/source/cice_1_introduction.rst +++ b/doc/source/cice_1_introduction.rst @@ -36,9 +36,164 @@ File and directory names are in **boldface**. A comprehensive :ref:`index`, including glossary of symbols with many of their values, appears at the end of this guide. -====================== -Major updates in V5.1 -====================== +================= +Quick Start guide +================= + +~~~~~~~~~~~~~ +Get the model +~~~~~~~~~~~~~ + +Checkout the model from the CICE-Consortium repository, + + github.com/CICE-Consortium + +For more details about how to work in github with CICE, a document can be +found `here `_. + +~~~~~~~~~~~~~~~~~ +Running the model +~~~~~~~~~~~~~~~~~ + +> cd consortium + +> ./create.case -c ~/mycase1 -g gx3 -m thunder -s diag1,thread -p 8x1 + +> cd ~/mycase1 + +> ./cice.build + +> ./cice.submit/Users/duvivier/Documents/Research/github/CICE-Consortium/CICE/doc/source/all_orig/cice_2_quick_start.rst + +~~~~~~~~~~~~ +More Details +~~~~~~~~~~~~ + +create.case generates a case, use "create.case -h" for help with the tool. + -c is the case name and location (required) + + -m is the machine name (required). Currently, there are working ports for NCAR yellowstone and cheyenne, AFRL thunder, NavyDSRC gordon and conrad, and LANL’s wolf machines. + + -g is the resolution (default is gx3) + + -p is the task x thread/task values (default is 4x1) + + -s are comma separated optional env or namelist settings (default is "null") + + -t is the test name and location (cannot be used with -c). + + -bd is used to specify the location of the baseline datasets (only used with -t) + + -bg is used to specify the cice version name for generating baseline datasets (only used with -t) + + -bc is used to specify the cice versoin name for comparison. I.e., the version name for the baseline dataset (only used with -t) + + -testid is used to specify a test ID (used only with -t or -ts) + + -ts is used to generate all test cases for a given test suite. + + +Several files are placed in the case directory + + - env.${machine} defines the environment + + - cice.settings defines many variables associated with building and running the model + + - makdep.c is a tool that will automatically generate the make dependencies + + - Macros.${machine} defines the Makefile Macros + + - Makefile is the makefile used to build the model + + - cice.build is a script that build the model + + - ice_in is the namelist file + + - cice.run is a batch run script + + - cice.submit is a simple script that submits the cice.run script + +Once the case is created, all scripts and namelist are fully resolved. Users can edit any +of the files in the case directory manually to change the model configuration. The file +dependency is indicated in the above list. For instance, if any of the files before +cice.build in the list are edited, cice.build should be rerun. + +The casescripts directory holds scripts used to create the case and can largely be ignored. + +In general, when cice.build is executed, the model will build from scratch due to the large +dependence on cpps. To change this behavior, edit the env variable ICE_CLEANBUILD in +cice.settings. + +The cice.submit script just submits the cice.run script. You can use cice.submit or just +submit the cice.run script on the command line. + +The model will run in the directory defined by the env variable CICE_RUNDIR in cice.settings. +Build and run logs will be copied into the case logs directory when complete. + +To port, an env.machine and Macros.machine file have to be added to scripts/machines and the cice.run.setup.csh file needs to be modified. + - cd to consortium/scripts/machines + - Copy an existing env and Macros file to new names for your new machine + - Edit the env and Macros file + - cd to consortium/scripts + - Edit the cice.run.setup.csh script to add a section for your machine for the batch settings and for the job launch settings + - Download and untar the 1997 dataset to the location defined by ICE_MACHINE_INPUTDATA in the env file + - Create a file in your home directory called .cice_proj and add your preferred account name to the first line. + - You can now create a case and test. If there are problems, you can manually edit the env, Macros, and cice.run files in the case directory until things are working properly. Then you can copy the env and Macros files back to consortium/scripts/machines. You will have to manually modify the cice.run.setup.csh script if there any changes needed there. + +~~~~~~~~~~~~ +Forcing data +~~~~~~~~~~~~ + +The code is currently configured to run in standalone mode on a 3 degree grid using +atmospheric data from 1997, available as detailed on the `wiki `_. +These data files are designed only for testing the code, not for use in production +runs or as observational data. Please do not publish results based on these data +sets. Module cicecore/dynamics/cicedynB/ice_forcing.F90 can be modified to change the +forcing data. + +As currently configured, the model runs on 4 processors. MPI is used for message passing +between processors, and OpenMP threading is available. The grid provided here is too +small for the code to scale well beyond about 8 processors. A 1 degree grid is provided also, +and details about this grid can be found on the `wiki `_. + +~~~~~~~~~~~~~~~~ +Online resources +~~~~~~~~~~~~~~~~ + +**DO WE WANT TO KEEP THESE?** + +primary wiki page: + + +FAQ: + + +instructions for code developers: + + +ongoing or planned development projects: + + +list of users and publications: + + +Please send references to your publications using the CICE model to ... + + +Please report any bugs to +Elizabeth Hunke (eclare@lanl.gov) + +Good luck! + + +============= +Major updates +============= + +~~~~~~~~~ +CICE V5.1 +~~~~~~~~~ + - include ice velocity in atm-ice coupling updates (e.g. stress) for high-frequency coupling - allow a variable coefficient for the ice-ocean heat flux - several new namelist options improve flexibility, especially for coupled model configurations: @@ -52,9 +207,10 @@ Major updates in V5.1 - new diagnostics and updated documentation - various bug fixes -====================== -Major updates in V5.0 -====================== +~~~~~~~~~ +CICE V5.0 +~~~~~~~~~ + - A method for prognosing sea ice salinity, including improved snow-ice formation - Two new explicit melt pond parameterizations (topo and level-ice) - Sea ice biogeochemistry @@ -74,3 +230,56 @@ Major updates in V5.0 - CPP options for categories, layers and tracers - Corrected bugs, particularly for nonstandard configurations. +====================== +Acknowledgements +====================== +This work has been supported under the Department of Energy’s Climate, +Ocean and Sea Ice Modeling project through the Computer Hardware Applied +Mathematics and Model Physics (CHAMMP) program, Climate Change +Prediction Program (CCPP), Improving the Characterization of Clouds, +Aerosols and the Cryosphere in Climate Models (Cloud-Cryo) program and +Scientific Discovery through Advanced Computing (SCIDAC) program, with +additional support from the T-3 Fluid Dynamics and Solid Mechanics Group +at Los Alamos National Laboratory. Special thanks are due to the +following people: + +- members of the CESM Polar Climate Working Group, including David + Bailey, Alice DuVivier, Cecilia Bitz, Bruce Briegleb, Tony Craig, + Marika Holland, John Dennis, Julie Schramm, Bonnie Light and Phil Jones. + +- Andrew Roberts of the Naval Postgraduate School, + +- David Hebert and Olivier Lecomte for their melt pond work, + +- Jonathan Gregory of the University of Reading and the U.K. MetOffice + for supplying tripole T-fold code and documentation, + +- Alison McLaren, Ann Keen and others working with the Hadley Centre + GCM for testing non-standard model configurations and providing their + code to us, + +- Daniel Feltham and his research group for several new + parameterizations and documentation, + +- Sylvain Bouillon for the revised EVP approach, + +- the many researchers who tested beta versions of CICE 5 and waited + patiently for the official release. + +====================== +Copyright +====================== +© Copyright 2013, LANS LLC. All rights reserved. Unless otherwise +indicated, this information has been authored by an employee or +employees of the Los Alamos National Security, LLC (LANS), operator of +the Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396 +with the U.S. Department of Energy. The U.S. Government has rights to +use, reproduce, and distribute this information. The public may copy and +use this information without charge, provided that this Notice and any +statement of authorship are reproduced on all copies. Neither the +Government nor LANS makes any warranty, express or implied, or assumes +any liability or responsibility for the use of this information. +Beginning with version 4.0, the CICE code carries Los Alamos Software +Release number LA-CC-06-012. + + diff --git a/doc/source/cice_2_quick_start.rst b/doc/source/cice_2_quick_start.rst deleted file mode 100644 index a44908152..000000000 --- a/doc/source/cice_2_quick_start.rst +++ /dev/null @@ -1,149 +0,0 @@ -Quick Start guide -============================================ - -Get the model -------------- - -Checkout the model from the CICE-Consortium repository, - - github.com/CICE-Consortium - -For more details about how to work in github with CICE, a document can be -found `here `_. - - -Running the model ------------------ - -> cd consortium - -> ./create.case -c ~/mycase1 -g gx3 -m thunder -s diag1,thread -p 8x1 - -> cd ~/mycase1 - -> ./cice.build - -> ./cice.submit - - -More Details: -------------- - -create.case generates a case, use "create.case -h" for help with the tool. - -c is the case name and location (required) - - -m is the machine name (required). Currently, there are working ports for NCAR yellowstone and cheyenne, AFRL thunder, NavyDSRC gordon and conrad, and LANL’s wolf machines. - - -g is the resolution (default is gx3) - - -p is the task x thread/task values (default is 4x1) - - -s are comma separated optional env or namelist settings (default is "null") - - -t is the test name and location (cannot be used with -c). - - -bd is used to specify the location of the baseline datasets (only used with -t) - - -bg is used to specify the cice version name for generating baseline datasets (only used with -t) - - -bc is used to specify the cice versoin name for comparison. I.e., the version name for the baseline dataset (only used with -t) - - -testid is used to specify a test ID (used only with -t or -ts) - - -ts is used to generate all test cases for a given test suite. - - -Several files are placed in the case directory - - - env.${machine} defines the environment - - - cice.settings defines many variables associated with building and running the model - - - makdep.c is a tool that will automatically generate the make dependencies - - - Macros.${machine} defines the Makefile Macros - - - Makefile is the makefile used to build the model - - - cice.build is a script that build the model - - - ice_in is the namelist file - - - cice.run is a batch run script - - - cice.submit is a simple script that submits the cice.run script - -Once the case is created, all scripts and namelist are fully resolved. Users can edit any -of the files in the case directory manually to change the model configuration. The file -dependency is indicated in the above list. For instance, if any of the files before -cice.build in the list are edited, cice.build should be rerun. - -The casescripts directory holds scripts used to create the case and can largely be ignored. - -In general, when cice.build is executed, the model will build from scratch due to the large -dependence on cpps. To change this behavior, edit the env variable ICE_CLEANBUILD in -cice.settings. - -The cice.submit script just submits the cice.run script. You can use cice.submit or just -submit the cice.run script on the command line. - -The model will run in the directory defined by the env variable CICE_RUNDIR in cice.settings. -Build and run logs will be copied into the case logs directory when complete. - -To port, an env.machine and Macros.machine file have to be added to scripts/machines and the cice.run.setup.csh file needs to be modified. - - cd to consortium/scripts/machines - - Copy an existing env and Macros file to new names for your new machine - - Edit the env and Macros file - - cd to consortium/scripts - - Edit the cice.run.setup.csh script to add a section for your machine for the batch settings and for the job launch settings - - Download and untar the 1997 dataset to the location defined by ICE_MACHINE_INPUTDATA in the env file - - Create a file in your home directory called .cice_proj and add your preferred account name to the first line. - - You can now create a case and test. If there are problems, you can manually edit the env, Macros, and cice.run files in the case directory until things are working properly. Then you can copy the env and Macros files back to consortium/scripts/machines. You will have to manually modify the cice.run.setup.csh script if there any changes needed there. - - -Forcing data ------------- - -The code is currently configured to run in standalone mode on a 3 degree grid using -atmospheric data from 1997, available as detailed on the `wiki `_. -These data files are designed only for testing the code, not for use in production -runs or as observational data. Please do not publish results based on these data -sets. Module cicecore/dynamics/cicedynB/ice_forcing.F90 can be modified to change the -forcing data. - -As currently configured, the model runs on 4 processors. MPI is used for message passing -between processors, and OpenMP threading is available. The grid provided here is too -small for the code to scale well beyond about 8 processors. A 1 degree grid is provided also, -and details about this grid can be found on the `wiki `_. - - -Online resources ----------------- - -**DO WE WANT TO KEEP THESE?** - -primary wiki page: - - -FAQ: - - -instructions for code developers: - - -ongoing or planned development projects: - - -list of users and publications: - - -Please send references to your publications using the CICE model to ... - - -Authors -------- - -Please report any bugs to -Elizabeth Hunke (eclare@lanl.gov) - -Good luck! diff --git a/doc/source/cice_4_model_components.rst b/doc/source/cice_2_science_guide.rst similarity index 86% rename from doc/source/cice_4_model_components.rst rename to doc/source/cice_2_science_guide.rst index bfed09704..5fb222326 100644 --- a/doc/source/cice_4_model_components.rst +++ b/doc/source/cice_2_science_guide.rst @@ -1,6 +1,564 @@ +Science Guide +================ + +.. _coupl: + +-------------------------------------------- +Coupling with other climate model components +-------------------------------------------- + +The sea ice model exchanges information with the other model components +via a flux coupler. CICE has been coupled into numerous climate models +with a variety of coupling techniques. This document is oriented +primarily toward the CESM Flux Coupler :cite:`KL02` +from NCAR, the first major climate model to incorporate CICE. The flux +coupler was originally intended to gather state variables from the +component models, compute fluxes at the model interfaces, and return +these fluxes to the component models for use in the next integration +period, maintaining conservation of momentum, heat, and fresh water. +However, several of these fluxes are now computed in the ice model +itself and provided to the flux coupler for distribution to the other +components, for two reasons. First, some of the fluxes depend strongly +on the state of the ice, and vice versa, implying that an implicit, +simultaneous determination of the ice state and the surface fluxes is +necessary for consistency and stability. Second, given the various ice +types in a single grid cell, it is more efficient for the ice model to +determine the net ice characteristics of the grid cell and provide the +resulting fluxes, rather than passing several values of the state +variables for each cell. These considerations are explained in more +detail below. + +The fluxes and state variables passed between the sea ice model and the +CESM flux coupler are listed in :ref:`tab-flux-cpl`. By convention, +directional fluxes are positive downward. In CESM, the sea ice model may +exchange coupling fluxes using a different grid than the computational +grid. This functionality is activated using the namelist variable +``gridcpl_file``. Another namelist variable ``highfreq``, allows the +high-frequency coupling procedure implemented in the Regional Arctic +System Model (RASM). In particular, the relative atmosphere-ice velocity +(:math:`\vec{U}_a-\vec{u}`) is used instead of the full atmospheric +velocity for computing turbulent fluxes in the atmospheric boundary +layer. + +:ref:`tab-flux-cpl`: *Data exchanged between the CESM flux coupler and the sea ice model* + +.. _tab-flux-cpl: + +.. table:: Table 1 + + =========================== ====================================== ======================================================================================= + Variable Description Interaction with flux coupler + =========================== ====================================== ======================================================================================= + :math:`z_o` Atmosphere level height From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`\vec{U}_a` Wind velocity From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`Q_a` Specific humidity From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`\rho_a` Air density From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`\Theta_a` Air potential temperature From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`T_a` Air temperature From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`F_{sw\downarrow}` Incoming shortwave radiation From *atmosphere model* via flux coupler **to** *sea ice model* + (4 bands) + + :math:`F_{L\downarrow}` Incoming longwave radiation From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`F_{rain}` Rainfall rate From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`F_{snow}` Snowfall rate From *atmosphere model* via flux coupler **to** *sea ice model* + + :math:`F_{frzmlt}` Freezing/melting potential From *ocean model* via flux coupler **to** *sea ice model* + + :math:`T_w` Sea surface temperature From *ocean model* via flux coupler **to** *sea ice model* + + :math:`S` Sea surface salinity From *ocean model* via flux coupler **to** *sea ice model* + + :math:`\nabla H_o` Sea surface slope From *ocean model* via flux coupler **to** *sea ice model* + + :math:`\vec{U}_w` Surface ocean currents From *ocean model* via flux coupler **to** *sea ice model* + + :math:`\vec{\tau}_a` Wind stress From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`F_s` Sensible heat flux From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`F_l` Latent heat flux From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`F_{L\uparrow}` Outgoing longwave radiation From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`F_{evap}` Evaporated water From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`\alpha` Surface albedo (4 bands) From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`T_{sfc}` Surface temperature From *sea ice model* via flux coupler **to** *atmosphere model* + + :math:`F_{sw\Downarrow}` Penetrating shortwave radiation From *sea ice model* via flux coupler **to** *ocean model* + + :math:`F_{water}` Fresh water flux From *sea ice model* via flux coupler **to** *ocean model* + + :math:`F_{hocn}` Net heat flux to ocean From *sea ice model* via flux coupler **to** *ocean model* + + :math:`F_{salt}` Salt flux From *sea ice model* via flux coupler **to** *ocean model* + + :math:`\vec{\tau}_w` Ice-ocean stress From *sea ice model* via flux coupler **to** *ocean model* + + :math:`F_{bio}` Biogeochemical fluxes From *sea ice model* via flux coupler **to** *ocean model* + + :math:`a_{i}` Ice fraction From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* + + :math:`T^{ref}_{a}` 2m reference temperature (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* + + :math:`Q^{ref}_{a}` 2m reference humidity (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* + + :math:`F_{swabs}` Absorbed shortwave (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* + =========================== ====================================== ======================================================================================= + +The ice fraction :math:`a_i` (aice) is the total fractional ice +coverage of a grid cell. That is, in each cell, + +.. math:: + \begin{array}{cl} + a_{i}=0 & \mbox{if there is no ice} \\ + a_{i}=1 & \mbox{if there is no open water} \\ + 0 0 + +where :math:`\cos Z` is the cosine of the solar zenith angle. + +.. _ocean: + +~~~~~ +Ocean +~~~~~ + +New sea ice forms when the ocean temperature drops below its freezing +temperature. In the Bitz and Lipscomb thermodynamics, +:cite:`BL99` :math:`T_f=-\mu S`, where :math:`S` is the +seawater salinity and :math:`\mu=0.054 \ ^\circ`/ppt is the ratio of the +freezing temperature of brine to its salinity (linear liquidus +approximation). For the mushy thermodynamics, :math:`T_f` is given by a +piecewise linear liquidus relation. The ocean model calculates the new +ice formation; if the freezing/melting potential +:math:`F_{frzmlt}` is positive, its value represents a certain +amount of frazil ice that has formed in one or more layers of the ocean +and floated to the surface. (The ocean model assumes that the amount of +new ice implied by the freezing potential actually forms.) + +If :math:`F_{frzmlt}` is negative, it is used to heat already +existing ice from below. In particular, the sea surface temperature and +salinity are used to compute an oceanic heat flux :math:`F_w` +(:math:`\left|F_w\right| \leq \left|F_{frzmlt}\right|`) which +is applied at the bottom of the ice. The portion of the melting +potential actually used to melt ice is returned to the coupler in +:math:`F_{hocn}`. The ocean model adjusts its own heat budget +with this quantity, assuming that the rest of the flux remained in the +ocean. + +In addition to runoff from rain and melted snow, the fresh water flux +:math:`F_{water}` includes ice melt water from the top surface +and water frozen (a negative flux) or melted at the bottom surface of +the ice. This flux is computed as the net change of fresh water in the +ice and snow volume over the coupling time step, excluding frazil ice +formation and newly accumulated snow. Setting the namelist option +update\_ocn\_f to true causes frazil ice to be included in the fresh +water and salt fluxes. + +There is a flux of salt into the ocean under melting conditions, and a +(negative) flux when sea water is freezing. However, melting sea ice +ultimately freshens the top ocean layer, since the ocean is much more +saline than the ice. The ice model passes the net flux of salt +:math:`F_{salt}` to the flux coupler, based on the net change +in salt for ice in all categories. In the present configuration, +ice\_ref\_salinity is used for computing the salt flux, although the ice +salinity used in the thermodynamic calculation has differing values in +the ice layers. + +A fraction of the incoming shortwave :math:`F_{sw\Downarrow}` +penetrates the snow and ice layers and passes into the ocean, as +described in Section :ref:`sfc-forcing`. + +Many ice models compute the sea surface slope :math:`\nabla H_\circ` +from geostrophic ocean currents provided by an ocean model or other data +source. In our case, the sea surface height :math:`H_\circ` is a +prognostic variable in POP—the flux coupler can provide the surface +slope directly, rather than inferring it from the currents. (The option +of computing it from the currents is provided in subroutine +*evp\_prep*.) The sea ice model uses the surface layer currents +:math:`\vec{U}_w` to determine the stress between the ocean and the ice, +and subsequently the ice velocity :math:`\vec{u}`. This stress, relative +to the ice, + +.. math:: + \begin{aligned} + \vec{\tau}_w&=&c_w\rho_w\left|{\vec{U}_w-\vec{u}}\right|\left[\left(\vec{U}_w-\vec{u}\right)\cos\theta + +\hat{k}\times\left(\vec{U}_w-\vec{u}\right)\sin\theta\right] \end{aligned} + +is then passed to the flux coupler (relative to the ocean) for use by +the ocean model. Here, :math:`\theta` is the turning angle between +geostrophic and surface currents, :math:`c_w` is the ocean drag +coefficient, :math:`\rho_w` is the density of seawater, and +:math:`\hat{k}` is the vertical unit vector. The turning angle is +necessary if the top ocean model layers are not able to resolve the +Ekman spiral in the boundary layer. If the top layer is sufficiently +thin compared to the typical depth of the Ekman spiral, then +:math:`\theta=0` is a good approximation. Here we assume that the top +layer is thin enough. + +For CICE run in stand-alone mode (i.e., uncoupled), a thermodynamic slab +ocean mixed-layer parameterization is available in **ice\_ocean.F90**. +The turbulent fluxes are computed above the water surface using the same +parameterizations as for sea ice, but with parameters appropriate for +the ocean. The surface flux balance takes into account the turbulent +fluxes, oceanic heat fluxes from below the mixed layer, and shortwave +and longwave radiation, including that passing through the sea ice into +the ocean. If the resulting sea surface temperature falls below the +salinity-dependent freezing point, then new ice (frazil) forms. +Otherwise, heat is made available for melting the ice. + +.. _formdrag: + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Variable exchange coefficients +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In the default CICE setup, atmospheric and oceanic neutral drag +coefficients (:math:`c_u` and :math:`c_w`) are assumed constant in time +and space. These constants are chosen to reflect friction associated +with an effective sea ice surface roughness at the ice–atmosphere and +ice–ocean interfaces. Sea ice (in both Arctic and Antarctic) contains +pressure ridges as well as floe and melt pond edges that act as discrete +obstructions to the flow of air or water past the ice, and are a source +of form drag. Following :cite:`TFSFFKLB14` and based on +recent theoretical developments :cite:`LGHA12,LLCL11`, the +neutral drag coefficients can now be estimated from properties of the +ice cover such as ice concentration, vertical extent and area of the +ridges, freeboard and floe draft, and size of floes and melt ponds. The +new parameterization allows the drag coefficients to be coupled to the +sea ice state and therefore to evolve spatially and temporally. This +parameterization is contained in the subroutine *neutral\_drag\_coeffs* +and is accessed by setting `formdrag` = true in the namelist. + +Following :cite:`TFSFFKLB14`, consider the general case of +fluid flow obstructed by N randomly oriented obstacles of height +:math:`H` and transverse length :math:`L_y`, distributed on a domain +surface area :math:`S_T`. Under the assumption of a logarithmic fluid +velocity profile, the general formulation of the form drag coefficient +can be expressed as + +.. math:: + C_d=\frac{N c S_c^2 \gamma L_y H}{2 S_T}\left[\frac{\ln(H/z_0)}{\ln(z_{ref}/z_0)}\right]^2, + :label: formdrag + +where :math:`z_0` is a roughness length parameter at the top or bottom +surface of the ice, :math:`\gamma` is a geometric factor, :math:`c` is +the resistance coefficient of a single obstacle, and :math:`S_c` is a +sheltering function that takes into account the shielding effect of the +obstacle, + +.. math:: + S_{c}=\left(1-\exp(-s_l D/H)\right)^{1/2}, + :label: shelter + +with :math:`D` the distance between two obstacles and :math:`s_l` an +attenuation parameter. + +As in the original drag formulation in CICE (sections :ref:`atmo` and +:ref:`ocean`), :math:`c_u` and :math:`c_w` along with the transfer +coefficients for sensible heat, :math:`c_{\theta}`, and latent heat, +:math:`c_{q}`, are initialized to a situation corresponding to neutral +atmosphere–ice and ocean–ice boundary layers. The corresponding neutral +exchange coefficients are then replaced by coefficients that explicitly +account for form drag, expressed in terms of various contributions as + +.. math:: + \tt{Cdn\_atm} = \tt{Cdn\_atm\_rdg} + \tt{Cdn\_atm\_floe} + \tt{Cdn\_atm\_skin} + \tt{Cdn\_atm\_pond} , + :label: Cda + +.. math:: + \tt{Cdn\_ocn} = \tt{Cdn\_ocn\_rdg} + \tt{Cdn\_ocn\_floe} + \tt{Cdn\_ocn\_skin}. + :label: Cdw + +The contributions to form drag from ridges (and keels underneath the +ice), floe edges and melt pond edges can be expressed using the general +formulation of equation :eq:`formdrag` (see :cite:`TFSFFKLB14` for +details). Individual terms in equation :eq:`Cdw` are fully described in +:cite:`TFSFFKLB14`. Following :cite:`Arya75` +the skin drag coefficient is parametrized as + +.. math:: + { \tt{Cdn\_(atm/ocn)\_skin}}=a_{i} \left(1-m_{(s/k)} \frac{H_{(s/k)}}{D_{(s/k)}}\right)c_{s(s/k)}, \mbox{ if $\displaystyle\frac{H_{(s/k)}}{D_{(s/k)}}\ge\frac{1}{m_{(s/k)}}$,} + :label: skindrag + +where :math:`m_s` (:math:`m_k`) is a sheltering parameter that depends +on the average sail (keel) height, :math:`H_s` (:math:`H_k`), but is +often assumed constant, :math:`D_s` (:math:`D_k`) is the average +distance between sails (keels), and :math:`c_{ss}` (:math:`c_{sk}`) is +the unobstructed atmospheric (oceanic) skin drag that would be attained +in the absence of sails (keels) and with complete ice coverage, +:math:`a_{ice}=1`. + +Calculation of equations :eq:`formdrag` – :eq:`skindrag` requires that small-scale geometrical +properties of the ice cover be related to average grid cell quantities +already computed in the sea ice model. These intermediate quantities are +briefly presented here and described in more detail in +:cite:`TFSFFKLB14`. The sail height is given by + +.. math:: + H_{s} = \displaystyle 2\frac{v_{rdg}}{a_{rdg}}\left(\frac{\alpha\tan \alpha_{k} R_d+\beta \tan \alpha_{s} R_h}{\phi_r\tan \alpha_{k} R_d+\phi_k \tan \alpha_{s} R_h^2}\right), + :label: Hs + +and the distance between sails\ + +.. math:: + D_{s} = \displaystyle 2 H_s\frac{a_{i}}{a_{rdg}} \left(\frac{\alpha}{\tan \alpha_s}+\frac{\beta}{\tan \alpha_k}\frac{R_h}{R_d}\right), + :label: Ds + +where :math:`0<\alpha<1` and :math:`0<\beta<1` are weight functions, +:math:`\alpha_{s}` and :math:`\alpha_{k}` are the sail and keel slope, +:math:`\phi_s` and :math:`\phi_k` are constant porosities for the sails +and keels, and we assume constant ratios for the average keel depth and +sail height (:math:`H_k/H_s=R_h`) and for the average distances between +keels and between sails (:math:`D_k/D_s=R_d`). With the assumption of +hydrostatic equilibrium, the effective ice plus snow freeboard is +:math:`H_{f}=\bar{h_i}(1-\rho_i/\rho_w)+\bar{h_s}(1-\rho_s/\rho_w)`, +where :math:`\rho_i`, :math:`\rho_w` and :math:`\rho_s` are +respectively the densities of sea ice, water and snow, :math:`\bar{h_i}` +is the mean ice thickness and :math:`\bar{h_s}` is the mean snow +thickness (means taken over the ice covered regions). For the melt pond +edge elevation we assume that the melt pond surface is at the same level +as the ocean surface surrounding the floes +:cite:`FF07,FFT10,FSFH12` and use the simplification +:math:`H_p = H_f`. Finally to estimate the typical floe size +:math:`L_A`, distance between floes, :math:`D_F`, and melt pond size, +:math:`L_P` we use the parameterizations of :cite:`LGHA12` +to relate these quantities to the ice and pond concentrations. All of +these intermediate quantities are available as history output, along +with `Cdn\_atm`, `Cdn\_ocn` and the ratio `Cdn\_atm\_ratio\_n` between the +total atmospheric drag and the atmospheric neutral drag coefficient. + +We assume that the total neutral drag coefficients are thickness +category independent, but through their dependance on the diagnostic +variables described above, they vary both spatially and temporally. The +total drag coefficients and heat transfer coefficients will also depend +on the type of stratification of the atmosphere and the ocean, and we +use the parameterization described in section :ref:`atmo` that accounts +for both stable and unstable atmosphere–ice boundary layers. In contrast +to the neutral drag coefficients the stability effect of the atmospheric +boundary layer is calculated separately for each ice thickness category. + +The transfer coefficient for oceanic heat flux to the bottom of the ice +may be varied based on form drag considerations by setting the namelist +variable `fbot\_xfer\_type` to `Cdn\_ocn`; this is recommended when using +the form drag parameterization. Its default value of the transfer +coefficient is 0.006 (`fbot\_xfer\_type = ’constant’`). + + +---------------- Model components -================ +---------------- The Arctic and Antarctic sea ice packs are mixtures of open water, thin first-year ice, thicker multiyear ice, and thick pressure ridges. The @@ -166,8 +724,9 @@ computed at the end of the last timestep are scaled for the new forcing. .. _tracers: +~~~~~~~ Tracers -------- +~~~~~~~ The basic conservation equations for ice area fraction :math:`a_{in}`, ice volume :math:`v_{in}`, and snow volume :math:`v_{sn}` for each @@ -233,8 +792,9 @@ guidance on adding tracers. .. _pondtr: +******************************************************* Tracers that depend on other tracers (e.g., melt ponds) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +******************************************************* Tracers may be defined that depend on other tracers. Melt pond tracers provide an example (these equations pertain to cesm and topo tracers; @@ -335,8 +895,9 @@ section :ref:`ponds`. .. _ice-age: +******* Ice age -~~~~~~~ +******* The age of the ice, :math:`\tau_{age}`, is treated as an ice-volume tracer (`trcr\_depend` = 1). It is initialized at 0 when ice @@ -359,8 +920,9 @@ is discussed in :cite:`ABTH11`. .. _ice-bgc: +*********************** Sea ice biogeochemistry -~~~~~~~~~~~~~~~~~~~~~~~ +*********************** Ice algal photosynthesis leads to carbon fixation and pigment buildup throughout much of the pack ice in springtime, including warm layers in @@ -587,8 +1149,9 @@ developed for future release in CICE. .. _aero: +******** Aerosols -~~~~~~~~ +******** Aerosols may be deposited on the ice and gradually work their way through it until the ice melts and they are passed into the ocean. They @@ -629,8 +1192,9 @@ oceanic fluxes, for each species. .. _brine-ht: +************ Brine height -~~~~~~~~~~~~ +************ The brine height, :math:`h_b`, is the distance from the ice–ocean interface to the brine surface. When `tr\_brine` is set true in @@ -773,8 +1337,9 @@ where the sums are taken over thickness categories. .. _horiz-trans: +~~~~~~~~~~~~~~~~~~~~ Horizontal transport --------------------- +~~~~~~~~~~~~~~~~~~~~ We wish to solve the continuity or transport equation (Equation :eq:`transport-ai`) for the fractional ice area in each @@ -863,8 +1428,9 @@ below. .. _reconstruct: +************************************* Reconstructing area and tracer fields -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +************************************* First, using the known values of the state variables, the ice area and tracer fields are reconstructed in each grid cell as linear functions of @@ -1054,8 +1620,9 @@ for generality. .. _loc-dep-triangles: +**************************** Locating departure triangles -~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +**************************** The method for locating departure triangles is discussed in detail by :cite:`DB00`. The basic idea is illustrated in @@ -1259,8 +1826,9 @@ trajectory. .. _integ-flux: +****************** Integrating fields -~~~~~~~~~~~~~~~~~~ +****************** Next, we integrate the reconstructed fields over the departure triangles to find the total area, volume, and energy transported across each cell @@ -1328,8 +1896,9 @@ repeatedly. .. _updating-state-var: +************************ Updating state variables -~~~~~~~~~~~~~~~~~~~~~~~~ +************************ Finally, we compute new values of the state variables in each ice category and grid cell. The new fractional ice areas @@ -1366,1103 +1935,1115 @@ values, with non-negative weights :math:`a` and :math:`ah`. Thus the new-time values must lie between the maximum and minimum of the old-time values. -.. _itd-trans: +.. _dynam: -Transport in thickness space ----------------------------- +~~~~~~~~ +Dynamics +~~~~~~~~ -Next we solve the equation for ice transport in thickness space due to -thermodynamic growth and melt, +There are now different rheologies available in the CICE code. The +elastic-viscous-plastic (EVP) model represents a modification of the +standard viscous-plastic (VP) model for sea ice dynamics +:cite:`Hibler79`. The elastic-anisotropic-plastic (EAP) model, +on the other hand, explicitly accounts for the observed sub-continuum +anisotropy of the sea ice cover :cite:`WF06,WS09`. If +`kdyn` = 1 in the namelist then the EVP rheology is used (module +**ice\_dyn\_evp.F90**), while `kdyn` = 2 is associated with the EAP +rheology (**ice\_dyn\_eap.F90**). At times scales associated with the +wind forcing, the EVP model reduces to the VP model while the EAP model +reduces to the anisotropic rheology described in detail in +:cite:`WF06,TFW13`. At shorter time scales the +adjustment process takes place in both models by a numerically more +efficient elastic wave mechanism. While retaining the essential physics, +this elastic wave modification leads to a fully explicit numerical +scheme which greatly improves the model’s computational efficiency. -.. math:: - \frac{\partial g}{\partial t} + \frac{\partial}{\partial h} (f g) = 0, - :label: itd-transport +The EVP sea ice dynamics model is thoroughly documented in +:cite:`HD97`, :cite:`Hunke01`, +:cite:`HD02` and :cite:`HD03` and the EAP +dynamics in :cite:`TFW13`. Simulation results and +performance of the EVP and EAP models have been compared with the VP +model and with each other in realistic simulations of the Arctic +respectively in :cite:`HZ99` and +:cite:`TFW13`. Here we summarize the equations and +direct the reader to the above references for details. The numerical +implementation in this code release is that of :cite:`HD02` +and :cite:`HD03`, with revisions to the numerical solver as +in :cite:`BFLM13`. The implementation of the EAP sea ice +dynamics into CICE is described in detail in +:cite:`TFW13`. -which is obtained from Equation :eq:`transport-g` by neglecting the first and -third terms on the right-hand side. We use the remapping method of -:cite:`Lipscomb01`, in which thickness categories are -represented as Lagrangian grid cells whose boundaries are projected -forward in time. The thickness distribution function :math:`g` is -approximated as a linear function of :math:`h` in each displaced -category and is then remapped onto the original thickness categories. -This method is numerically smooth and is not too diffusive. It can be -viewed as a 1D simplification of the 2D incremental remapping scheme -described above. +.. _momentum: -We first compute the displacement of category boundaries in thickness -space. Assume that at time :math:`m` the ice areas :math:`a_n^m` and -mean ice thicknesses :math:`h_n^m` are known for each thickness -category. (For now we omit the subscript :math:`i` that distinguishes -ice from snow.) We use a thermodynamic model (Section :ref:`thermo`) -to compute the new mean thicknesses :math:`h_n^{m+1}` at time -:math:`m+1`. The time step must be small enough that trajectories do not -cross; i.e., :math:`h_n^{m+1} < h_{n+1}^{m+1}` for each pair of adjacent -categories. The growth rate at :math:`h = h_n` is given by -:math:`f_n = (h_n^{m+1} - h_n^m) / \Delta t`. By linear interpolation we -estimate the growth rate :math:`F_n` at the upper category boundary -:math:`H_n`: +******** +Momentum +******** -.. math:: - F_n = f_n + \frac{f_{n+1}-f_n}{h_{n+1}-h_n} \, (H_n - h_n). +The force balance per unit area in the ice pack is given by a +two-dimensional momentum equation :cite:`Hibler79`, obtained +by integrating the 3D equation through the thickness of the ice in the +vertical direction: -If :math:`a_n` or :math:`a_{n+1} = 0`, :math:`F_n` is set to the growth -rate in the nonzero category, and if :math:`a_n = a_{n+1} = 0`, we set -:math:`F_n = 0`. The temporary displaced boundaries are given by +.. math:: + m{\partial {\bf u}\over\partial t} = \nabla\cdot{\bf \sigma} + + \vec{\tau}_a+\vec{\tau}_w - \hat{k}\times mf{\bf u} - mg\nabla H_\circ, + :label: vpmom -.. math:: - H_n^* = H_n + F_n \, \Delta t, \ n = 1 \ {\rm to} \ N-1 +where :math:`m` is the combined mass of ice and snow per unit area and +:math:`\vec{\tau}_a` and :math:`\vec{\tau}_w` are wind and ocean +stresses, respectively. The strength of the ice is represented by the +internal stress tensor :math:`\sigma_{ij}`, and the other two terms on +the right hand side are stresses due to Coriolis effects and the sea +surface slope. The parameterization for the wind and ice–ocean stress +terms must contain the ice concentration as a multiplicative factor to +be consistent with the formal theory of free drift in low ice +concentration regions. A careful explanation of the issue and its +continuum solution is provided in :cite:`HD03` and +:cite:CGHM04`. -The boundaries must not be displaced by more than one category to the -left or right; that is, we require :math:`H_{n-1} < H_n^* < H_{n+1}`. -Without this requirement we would need to do a general remapping rather -than an incremental remapping, at the cost of added complexity. +The momentum equation is discretized in time as follows, for the classic +EVP approach. First, for clarity, the two components of Equation :eq:`vpmom` are -Next we construct :math:`g(h)` in the displaced thickness categories. -The ice areas in the displaced categories are :math:`a_n^{m+1} = a_n^m`, -since area is conserved following the motion in thickness space (i.e., -during vertical ice growth or melting). The new ice volumes are -:math:`v_n^{m+1} = (a_n h_n)^{m+1} = a_n^m h_n^{m+1}`. For conciseness, -define :math:`H_L = H_{n-1}^*` and :math:`H_R = H_{n}^*` and drop the -time index :math:`m+1`. We wish to construct a continuous function -:math:`g(h)` within each category such that the total area and volume at -time :math:`m+1` are :math:`a_n` and :math:`v_n`, respectively: +.. math:: + \begin{aligned} + m{\partial u\over\partial t} &=& {\partial\sigma_{1j}\over\partial x_j} + \tau_{ax} + + a_i c_w \rho_w + \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\cos\theta - \left(V_w-v\right)\sin\theta\right] + +mfv - mg{\partial H_\circ\over\partial x}, \\ + m{\partial v\over\partial t} &=& {\partial\sigma_{2j}\over\partial x_j} + \tau_{ay} + + a_i c_w \rho_w + \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\sin\theta - \left(V_w-v\right)\cos\theta\right] + -mfu - mg{\partial H_\circ\over\partial y}. \end{aligned} + +In the code, +:math:`{\tt vrel}=a_i c_w \rho_w\left|{\bf U}_w - {\bf u}^k\right|`, +where :math:`k` denotes the subcycling step. The following equations +illustrate the time discretization and define some of the other +variables used in the code. .. math:: - \int_{H_L}^{H_R} g \, dh = a_n, - :label: area-cons + \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta\right)}_{\tt cca} u^{k+1} + - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} + = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} + + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} + + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t_e}u^k, + :label: umom .. math:: - \int_{H_L}^{H_R} h \, g \, dh = v_n. - :label: volume-cons + \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} + + \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta\right)}_{\tt cca}v^{k+1} + = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} + + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} + + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t_e}v^k, + :label: vmom -The simplest polynomial that can satisfy both equations is a line. It -is convenient to change coordinates, writing -:math:`g(\eta) = g_1 \eta + g_0`, where :math:`\eta = h - H_L` and the -coefficients :math:`g_0` and :math:`g_1` are to be determined. Then -Equations :eq:`area-cons` and :eq:`volume-cons` can be written as +and vrel\ :math:`\cdot`\ waterx(y) = taux(y). -.. math:: - g_1 \frac{\eta_R^2}{2} + g_0 \eta_R = a_n, - -.. math:: - g_1 \frac{\eta_R^3}{3} + g_0 \frac{\eta_R^2}{2} = a_n \eta_n, - -where :math:`\eta_R = H_R - H_L` and :math:`\eta_n = h_n - H_L`. These -equations have the solution +We solve this system of equations analytically for :math:`u^{k+1}` and +:math:`v^{k+1}`. Define .. math:: - g_0 = \frac{6 a_n}{\eta_R^2} \left(\frac{2 \eta_R}{3} - \eta_n\right), - :label: g0 + \hat{u} = F_u + \tau_{ax} - mg{\partial H_\circ\over\partial x} + {\tt vrel} \left(U_w\cos\theta - V_w\sin\theta\right) + {m\over\Delta t_e}u^k + :label: cevpuhat .. math:: - g_1 = \frac{12 a_n}{\eta_R^3} \left(\eta_n - \frac{\eta_R}{2}\right). - :label: g1 + \hat{v} = F_v + \tau_{ay} - mg{\partial H_\circ\over\partial y} + {\tt vrel} \left(U_w\sin\theta + V_w\cos\theta\right) + {m\over\Delta t_e}v^k, + :label: cevpvhat -Since :math:`g` is linear, its maximum and minimum values lie at the -boundaries, :math:`\eta = 0` and :math:`\eta_R`: +where :math:`{\bf F} = \nabla\cdot\sigma^{k+1}`. Then .. math:: - g(0)=\frac{6 a_n}{\eta_R^2} \, \left(\frac{2 \eta_R}{3} - \eta_n\right) = g_0, - :label: gmin - -.. math:: - g(\eta_R) = \frac{6 a_n}{\eta_R^2} \, \left(\eta_n - \frac{\eta_R}{3}\right). - :label: gmax + \begin{aligned} + \left({m\over\Delta t_e} +{\tt vrel}\cos\theta\right)u^{k+1} - \left(mf + {\tt vrel}\sin\theta\right) v^{k+1} &=& \hat{u} \\ + \left(mf + {\tt vrel}\sin\theta\right) u^{k+1} + \left({m\over\Delta t_e} +{\tt vrel}\cos\theta\right)v^{k+1} &=& \hat{v}.\end{aligned} -Equation :eq:`gmin` implies that :math:`g(0) < 0` when -:math:`\eta_n > 2 \eta_R/3`, i.e., when :math:`h_n` lies in the right -third of the thickness range :math:`(H_L, H_R)`. Similarly, Equation :eq:`gmax` -implies that :math:`g(\eta_R) < 0` when :math:`\eta_n < \eta_R/3`, i.e., -when :math:`h_n` is in the left third of the range. Since negative -values of :math:`g` are unphysical, a different solution is needed when -:math:`h_n` lies outside the central third of the thickness range. If -:math:`h_n` is in the left third of the range, we define a cutoff -thickness, :math:`H_C = 3 h_n - 2 H_L`, and set :math:`g = 0` between -:math:`H_C` and :math:`H_R`. Equations :eq:`g0` and :eq:`g1` are then -valid with :math:`\eta_R` redefined as :math:`H_C - H_L`. And if -:math:`h_n` is in the right third of the range, we define -:math:`H_C = 3 h_n - 2 H_R` and set :math:`g = 0` between :math:`H_L` -and :math:`H_C`. In this case, :eq:`g0` and :eq:`g1` apply with -:math:`\eta_R = H_R - H_C` and :math:`\eta_n = h_n - H_C`. +Solving simultaneously for :math:`u^{k+1}` and :math:`v^{k+1}`, -:ref:`fig-gplot` illustrates the linear reconstruction of :math:`g` -for the simple cases :math:`H_L = 0`, :math:`H_R = 1`, :math:`a_n = 1`, -and :math:`h_n =` 0.2, 0.4, 0.6, and 0.8. Note that :math:`g` slopes -downward (:math:`g_1 < 0`) when :math:`h_n` is less than the midpoint -thickness, :math:`(H_L + H_R)/2 = 1/2`, and upward when :math:`h_n` -exceeds the midpoint thickness. For :math:`h_n = 0.2` and 0.8, -:math:`g = 0` over part of the range. +.. math:: + \begin{aligned} + u^{k+1} = {a \hat{u} + b \hat{v} \over a^2 + b^2} \\ + v^{k+1} = {a \hat{v} - b \hat{u} \over a^2 + b^2}, \end{aligned} -.. _fig-gplot: +where -.. figure:: ./figures/gplot.png - :align: center - :scale: 20% +.. math:: + a = {m\over\Delta t_e} + {\tt vrel}\cos\theta \\ + :label: cevpa - Figure 4 +.. math:: + b = mf + {\tt vrel}\sin\theta. + :label: cevpb -:ref:`fig-gplot` : Linear approximation of the thickness distribution -function :math:`g(h)` for an ice category with left boundary -:math:`H_L = 0`, right boundary :math:`H_R = 1`, fractional area -:math:`a_n = 1`, and mean ice thickness :math:`h_n = 0.2, 0.4, 0.6,` and :math:`0.8`. +When the subcycling is finished for each (thermodynamic) time step, the +ice–ocean stress must be constructed from `taux(y)` and the terms +containing `vrel` on the left hand side of the equations. -Finally, we remap the thickness distribution to the original boundaries -by transferring area and volume between categories. We compute the ice -area :math:`\Delta a_n` and volume :math:`\Delta v_n` between each -original boundary :math:`H_n` and displaced boundary :math:`H_n^*`. If -:math:`H_n^* > H_n`, ice moves from category :math:`n` to :math:`n+1`. -The area and volume transferred are +The Hibler-Bryan form for the ice-ocean stress :cite:`HB87` +is included in **ice\_dyn\_shared.F90** but is currently commented out, +pending further testing. -.. math:: - \Delta a_n = \int_{H_n}^{H_n^*} g \, dh, - :label: move-area +.. _internal-stress: -.. math:: - \Delta v_n = \int_{H_n}^{H_n^*} h \, g \, dh. - :label: move-volume +*************** +Internal stress +*************** -If :math:`H_n^* < H_N`, ice area and volume are transferred from -category :math:`n+1` to :math:`n` using Equations :eq:`move-area` and -:eq:`move-volume` with the limits of integration reversed. To evaluate -the integrals we change coordinates from :math:`h` to -:math:`\eta = h - H_L`, where :math:`H_L` is the left limit of the range -over which :math:`g > 0`, and write :math:`g(\eta)` using Equations :eq:`g0` and -:eq:`g1`. In this way we obtain the new areas :math:`a_n` and volumes -:math:`v_n` between the original boundaries :math:`H_{n-1}` and -:math:`H_n` in each category. The new thicknesses, -:math:`h_n = v_n/a_n`, are guaranteed to lie in the range -:math:`(H_{n-1}, H_n)`. If :math:`g = 0` in the part of a category that -is remapped to a neighboring category, no ice is transferred. +For convenience we formulate the stress tensor :math:`\bf \sigma` in +terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, +:math:`\sigma_2=\sigma_{11}-\sigma_{22}`, and introduce the +divergence, :math:`D_D`, and the horizontal tension and shearing +strain rates, :math:`D_T` and :math:`D_S` respectively. -Other conserved quantities are transferred in proportion to the ice -volume :math:`\Delta v_{in}`. For example, the transferred ice energy in -layer :math:`k` is -:math:`\Delta e_{ink} = e_{ink} (\Delta v_{in} / v_{in})`. +*Elastic-Viscous-Plastic* -The left and right boundaries of the domain require special treatment. -If ice is growing in open water at a rate :math:`F_0`, the left boundary -:math:`H_0` is shifted to the right by :math:`F_0 \Delta t` before -:math:`g` is constructed in category 1, then reset to zero after the -remapping is complete. New ice is then added to the grid cell, -conserving area, volume, and energy. If ice cannot grow in open water -(because the ocean is too warm or the net surface energy flux is -downward), :math:`H_0` is fixed at zero, and the growth rate at the left -boundary is estimated as :math:`F_0 = f_1`. If :math:`F_0 < 0`, all ice -thinner than :math:`\Delta h_0 = -F_0 \Delta t` is assumed to have -melted, and the ice area in category 1 is reduced accordingly. The area -of new open water is +In the EVP model the internal stress tensor is determined from a +regularized version of the VP constitutive law, -.. math:: - \Delta a_0 = \int_{0}^{\Delta h_0} g \, dh. +.. math:: + {1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} + + {P\over 2\zeta} = D_D, \\ + :label: sig1 -The right boundary :math:`H_N` is not fixed but varies with -:math:`h_N`, the mean ice thickness in the thickest category. Given -:math:`h_N`, we set :math:`H_N = 3 h_N - 2 H_{N-1}`, which ensures that -:math:`g(h) > 0` for :math:`H_{N-1} < h < H_N` and :math:`g(h) = 0` for -:math:`h \geq H_N`. No ice crosses the right boundary. If the ice growth -or melt rates in a given grid cell are too large, the thickness -remapping scheme will not work. Instead, the thickness categories in -that grid cell are treated as delta functions following -:cite:`BHWE01`, and categories outside their prescribed -boundaries are merged with neighboring categories as needed. For time -steps of less than a day and category thickness ranges of 10 cm or more, -this simplification is needed rarely, if ever. +.. math:: + {1\over E}{\partial\sigma_2\over\partial t} + {\sigma_2\over 2\eta} = D_T, + :label: sig2 -The linear remapping algorithm for thickness is not monotonic for -tracers, although significant errors rarely occur. Usually they appear -as snow temperatures (enthalpy) outside the physical range of values in -very small snow volumes. In this case we transfer the snow and its heat -and tracer contents to the ocean. +.. math:: + {1\over E}{\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over + 2\eta} = {1\over 2}D_S, + :label: sig12 -.. _mech-red: +where -Mechanical redistribution -------------------------- +.. math:: + D_D = \dot{\epsilon}_{11} + \dot{\epsilon}_{22}, -The last term on the right-hand side of Equation :eq:`transport-g` -is :math:`\psi`, which describes the redistribution -of ice in thickness space due to ridging and other mechanical processes. -The mechanical redistribution scheme in CICE is based on -:cite:`TRMC75`, :cite:`Rothrock75`, -:cite:`Hibler80`, :cite:`FH95`, and -:cite:`LHMJ07`. This scheme converts thinner ice to thicker -ice and is applied after horizontal transport. When the ice is -converging, enough ice ridges to ensure that the ice area does not -exceed the grid cell area. +.. math:: + D_T = \dot{\epsilon}_{11} - \dot{\epsilon}_{22}, -First we specify the participation function: the thickness distribution -:math:`a_P(h) = b(h) \, g(h)` of the ice participating in ridging. (We -use “ridging” as shorthand for all forms of mechanical redistribution, -including rafting.) The weighting function :math:`b(h)` favors ridging -of thin ice and closing of open water in preference to ridging of -thicker ice. There are two options for the form of :math:`b(h)`. If -`krdg\_partic` = 0 in the namelist, we follow :cite:`TRMC75` -and set +.. math:: + D_S = 2\dot{\epsilon}_{12}, .. math:: - b(h) = \left\{\begin{array}{ll} - \frac{2}{G^*}(1-\frac{G(h)}{G^*}) & \mbox{if $G(h) G^*`, then :math:`a_{Pn} = 0`. If -the open water fraction :math:`a_0 > G^*`, no ice can ridge, because -“ridging” simply reduces the area of open water. As in -:cite:`TRMC75` we set :math:`G^* = 0.15`. +.. math:: + \Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2}, -If the spatial resolution is too fine for a given time step -:math:`\Delta t`, the weighting function Equation :eq:`partic-old-contin` can -promote numerical instability. For :math:`\Delta t = \mbox{1 hour}`, -resolutions finer than :math:`\Delta x \sim \mbox{10 km}` are typically -unstable. The instability results from feedback between the ridging -scheme and the dynamics via the ice strength. If the strength changes -significantly on time scales less than :math:`\Delta t`, the -viscous-plastic solution of the momentum equation is inaccurate and -sometimes oscillatory. As a result, the fields of ice area, thickness, -velocity, strength, divergence, and shear can become noisy and -unphysical. +and :math:`P` is a function of the ice thickness and concentration, +described in Section :ref:`mech-red`. The dynamics component +employs a “replacement pressure” (see :cite:`GHA98`, for +example), which serves to prevent residual ice motion due to spatial +variations of :math:`P` when the rates of strain are exactly zero. -A more stable weighting function was suggested by -:cite:`LHMJ07`: +Viscosities are updated during the subcycling, so that the entire +dynamics component is subcycled within the time step, and the elastic +parameter :math:`E` is defined in terms of a damping timescale :math:`T` +for elastic waves, :math:`\Delta t_e < T < \Delta t`, as -.. math:: - b(h) = \frac{\exp[-G(h)/a^*]} - {a^*[1-\exp(-1/a^*)]} - :label: partic-new-contin +.. math:: + E = {\zeta\over T}, -When integrated between category boundaries, Equation :eq:`partic-new-contin` -implies +where :math:`T=E_\circ\Delta t` and :math:`E_\circ` (eyc) is a tunable +parameter less than one. The stress equations :eq:`sig1`–:eq:`sig12` +become .. math:: - a_{Pn} = \frac {\exp(-G_{n-1}/a^*) - \exp(-G_{n}/a^*)} - {1 - \exp(-1/a^*)} - :label: partic-new-discrete + \begin{aligned} + {\partial\sigma_1\over\partial t} + {\sigma_1\over 2T} + + {P\over 2T} &=& {P\over 2T\Delta} D_D, \\ + {\partial\sigma_2\over\partial t} + {e^2\sigma_2\over 2T} &=& {P\over + 2T\Delta} D_T,\\ + {\partial\sigma_{12}\over\partial t} + {e^2\sigma_{12}\over 2T} &=& + {P\over 4T\Delta}D_S.\end{aligned} -This weighting function is used if `krdg\_partic` = 1 in the namelist. -From Equation :eq:`partic-new-contin`, the mean value of :math:`G` for ice -participating in ridging is :math:`a^*`, as compared to :math:`G^*/3` -for Equation :eq:`partic-old-contin`. For typical ice thickness distributions, -setting :math:`a^* = 0.05` with `krdg\_partic` = 1 gives participation -fractions similar to those given by :math:`G^* = 0.15` with `krdg\_partic` -= 0. See :cite:`LHMJ07` for a detailed comparison of these -two participation functions. - -Thin ice is converted to thick, ridged ice in a way that reduces the -total ice area while conserving ice volume and internal energy. There -are two namelist options for redistributing ice among thickness -categories. If `krdg\_redist` = 0, ridging ice of thickness :math:`h_n` -forms ridges whose area is distributed uniformly between -:math:`H_{\min} = 2 h_n` and :math:`H_{\max} = 2 \sqrt{H^* h_n}`, as in -:cite:`Hibler80`. The default value of :math:`H^*` is 25 m, as -in earlier versions of CICE. Observations suggest that -:math:`H^* = 50` m gives a better fit to first-year ridges -:cite:`AMI04`, although the lower value may be appropriate -for multiyear ridges :cite:`FH95`. The ratio of the mean -ridge thickness to the thickness of ridging ice is -:math:`k_n = (H_{\min} + H_{\max}) / (2 h_n)`. If the area of category -:math:`n` is reduced by ridging at the rate :math:`r_n`, the area of -thicker categories grows simultaneously at the rate :math:`r_n/k_n`. -Thus the *net* rate of area loss due to ridging of ice in category -:math:`n` is :math:`r_n(1-1/k_n)`. - -The ridged ice area and volume are apportioned among categories in the -thickness range :math:`(H_{\min}, H_{\max})`. The fraction of the new -ridge area in category :math:`m` is - -.. math:: - f_m^{\mathrm{area}} = \frac{H_R - H_L} - {H_{\max} - H_{\min}}, - :label: ridge-area-old +All coefficients on the left-hand side are constant except for +:math:`P`, which changes only on the longer time step :math:`\Delta t`. +This modification compensates for the decreased efficiency of including +the viscosity terms in the subcycling. (Note that the viscosities do not +appear explicitly.) Choices of the parameters used to define :math:`E`, +:math:`T` and :math:`\Delta t_e` are discussed in +Sections :ref:`revp` and :ref:`parameters`. -where :math:`H_L = \max(H_{m-1},H_{\min})` and -:math:`H_R= \min(H_m,H_{\max})`. The fraction of the ridge volume going -to category :math:`m` is +The bilinear discretization used for the stress terms +:math:`\partial\sigma_{ij}/\partial x_j` in the momentum equation is +now used, which enabled the discrete equations to be derived from the +continuous equations written in curvilinear coordinates. In this +manner, metric terms associated with the curvature of the grid are +incorporated into the discretization explicitly. Details pertaining to +the spatial discretization are found in :cite:`HD02`. -.. math:: - f_m^{\mathrm{vol}} = \frac{(H_R)^2 - (H_L)^2} - {(H_{\max})^2 - (H_{\min})^2}. - :label: ridge-volume-old +*Elastic-Anisotropic-Plastic* -This uniform redistribution function tends to produce too little ice in -the 3–5 m range and too much ice thicker than 10 m -:cite:`AMI04`. Observations show that the ITD of ridges is -better approximated by a negative exponential. Setting `krdg\_redist` = 1 -gives ridges with an exponential ITD :cite:`LHMJ07`: +In the EAP model the internal stress tensor is related to the +geometrical properties and orientation of underlying virtual diamond +shaped floes (see :ref:`fig-EAP`). In contrast to the isotropic EVP +rheology, the anisotropic plastic yield curve within the EAP rheology +depends on the relative orientation of the diamond shaped floes (unit +vector :math:`\mathbf r` in :ref:`fig-EAP`), with respect to the +principal direction of the deformation rate (not shown). Local +anisotropy of the sea ice cover is accounted for by an additional +prognostic variable, the structure tensor :math:`\mathbf{A}` defined +by -.. math:: - g_R(h) \propto \exp[-(h - H_{\min})/\lambda] - :label: redist-new +.. math:: + {\mathbf A}=\int_{\mathbb{S}}\vartheta(\mathbf r)\mathbf r\mathbf r d\mathbf r\label{structuretensor}. -for :math:`h \ge H_{\min}`, with :math:`g_R(h) = 0` for -:math:`h < H_{\min}`. Here, :math:`\lambda` is an empirical *e*-folding -scale and :math:`H_{\min}=2h_n` (where :math:`h_n` is the thickness of -ridging ice). We assume that :math:`\lambda = \mu h_n^{1/2}`, where -:math:`\mu` (mu\_rdg) is a tunable parameter with units . Thus the mean -ridge thickness increases in proportion to :math:`h_n^{1/2}`, as in -:cite:`Hibler80`. The value :math:`\mu = 4.0`  gives -:math:`\lambda` in the range 1–4 m for most ridged ice. Ice strengths -with :math:`\mu = 4.0`  and `krdg\_redist` = 1 are roughly comparable to -the strengths with :math:`H^* = 50` m and `krdg\_redist` = 0. +where :math:`\mathbb{S}` is a unit-radius circle; **A** is a unit +trace, 2\ :math:`\times`\ 2 matrix. From now on we shall describe the +orientational distribution of floes using the structure tensor. For +simplicity we take the probability density function +:math:`\vartheta(\mathbf r )` to be Gaussian, +:math:`\vartheta(z)=\omega_{1}\exp(-\omega_{2}z^{2})`, where :math:`z` +is the ice floe inclination with respect to the axis :math:`x_{1}` of +preferential alignment of ice floes (see :ref:`fig-EAP`), +:math:`\vartheta(z)` is periodic with period :math:`\pi`, and the +positive coefficients :math:`\omega_{1}` and :math:`\omega_{2}` are +calculated to ensure normalization of :math:`\vartheta(z)`, i.e. +:math:`\int_{0}^{2\pi}\vartheta(z)dz=1`. The ratio of the principal +components of :math:`\mathbf{A}`, :math:`A_{1}/A_{2}`, are derived +from the phenomenological evolution equation for the structure tensor +:math:`\mathbf A`, -From Equation :eq:`redist-new` it can be shown that the fractional area going -to category :math:`m` as a result of ridging is +.. math:: + \frac{D\mathbf{A}}{D t}=\mathbf{F}_{iso}(\mathbf{A})+\mathbf{F}_{frac}(\mathbf{A},\boldsymbol\sigma), + :label: evolutionA -.. math:: - f_m^{\mathrm{area}} = \exp[-(H_{m-1} - H_{\min}) / \lambda] - - \exp[-(H_m - H_{\min}) / \lambda]. - :label: ridge-area-new +where :math:`t` is the time, and :math:`D/Dt` is the co-rotational +time derivative accounting for advection and rigid body rotation +(:math:`D\mathbf A/Dt = d\mathbf A/dt -\mathbf W \cdot \mathbf A -\mathbf A \cdot \mathbf W^{T}`) +with :math:`\mathbf W` being the vorticity tensor. +:math:`\mathbf F_{iso}` is a function that accounts for a variety of +processes (thermal cracking, melting, freezing together of floes) that +contribute to a more isotropic nature to the ice cover. +:math:`\mathbf F_{frac}` is a function determining the ice floe +re-orientation due to fracture, and explicitly depends upon sea ice +stress (but not its magnitude). Following :cite:`WF06`, +based on laboratory experiments by :cite:`Schulson01` we +consider four failure mechanisms for the Arctic sea ice cover. These +are determined by the ratio of the principal values of the sea ice +stress :math:`\sigma_{1}` and :math:`\sigma_{2}`: (i) under biaxial +tension, fractures form across the perpendicular principal axes and +therefore counteract any apparent redistribution of the floe +orientation; (ii) if only one of the principal stresses is +compressive, failure occurs through axial splitting along the +compression direction; (iii) under biaxial compression with a low +confinement ratio, (:math:`\sigma_{1}/\sigma_{2} 1`, the ridging is -repeated with a value of :math:`R_{\mathrm{net}}` sufficient to yield -:math:`a_i = 1`. +.. math:: + \boldsymbol\sigma^{EAP}(h)=P_{r}(h)\int_{\mathbb{S}}\vartheta(\mathbf r)\left[\boldsymbol\sigma_{r}^{b}(\mathbf r)+ k \boldsymbol\sigma_{s}^{b}(\mathbf r)\right]d\mathbf r + :label: stressaverage -Two tracers for tracking the ridged ice area and volume are available. -The actual tracers are for level (undeformed) ice area (`alvl`) and volume -(`vlvl`), which are easier to implement for a couple of reasons: (1) ice -ridged in a given thickness category is spread out among the rest of the -categories, making it more difficult (and expensive) to track than the -level ice remaining behind in the original category; (2) previously -ridged ice may ridge again, so that simply adding a volume of freshly -ridged ice to the volume of previously ridged ice in a grid cell may be -inappropriate. Although the code currently only tracks level ice -internally, both level ice and ridged ice are offered as history output. -They are simply related: +where we have introduced the friction parameter :math:`k=P_{s}/P_{r}` +and where we identify the ridging ice strength :math:`P_{r}(h)` with the +strength :math:`P` described in section 1 and used within the EVP +framework. -.. math:: - \begin{aligned} - a_{lvl} + a_{rdg} &=& a_i, \\ - v_{lvl} + v_{rdg} &=& v_i.\end{aligned} +As is the case for the EVP rheology, elasticity is included in the EAP +description not to describe any physical effect, but to make use of the +efficient, explicit numerical algorithm used to solve the full sea ice +momentum balance. We use the analogous EAP stress equations, -Level ice area fraction and volume increase with new ice formation and -decrease steadily via ridging processes. Without the formation of new -ice, level ice asymptotes to zero because we assume that both level ice -and ridged ice ridge, in proportion to their fractional areas in a grid -cell (in the spirit of the ridging calculation itself which does not -prefer level ice over previously ridged ice). +.. math:: + \frac{\partial \sigma_{1}}{\partial t}+\frac{\sigma_1}{2T} = \frac{\sigma^{EAP}_{1}}{2T} \mbox{,} + :label: EAPsigma1 -The ice strength :math:`P` may be computed in either of two ways. If the -namelist parameter kstrength = 0, we use the strength formula from -:cite:`Hibler79`: +.. math:: + \frac{\partial \sigma_{2}}{\partial t}+\frac{\sigma_2}{2T} = \frac{\sigma^{EAP}_{2}}{2T} \mbox{,} + :label: EAPsigma2 .. math:: - P = P^* h \exp[-C(1-a_i)], - :label: hib-strength + \frac{\partial \sigma_{12}}{\partial t}+\frac{\sigma_{12}}{2T} = \frac{\sigma^{EAP}_{12}}{2T} \mbox{,} + :label: EAPsigma12 -where :math:`P^* = 27,500 \, \mathrm {N/m}` and :math:`C = 20` are -empirical constants, and :math:`h` is the mean ice thickness. -Alternatively, setting kstrength = 1 gives an ice strength closely -related to the ridging scheme. Following -:cite:`Rothrock75`, the strength is assumed proportional -to the change in ice potential energy :math:`\Delta E_P` per unit area -of compressive deformation. Given uniform ridge ITDs (krdg\_redist = 0), -we have +where the anisotropic stress :math:`\boldsymbol\sigma^{EAP}` is defined +in a look-up table for the current values of strain rate and structure +tensor. The look-up table is constructed by computing the stress +(normalized by the strength) from Equations :eq:`EAPsigma1`–:eq:`EAPsigma12` +for discrete values of the largest eigenvalue of the structure tensor, +:math:`\frac{1}{2}\le A_{1}\le 1`, the angle :math:`0\le\theta\le2\pi`, +and the angle :math:`-\pi/2\le y\le\pi/2` between the major principal +axis of the strain rate tensor and the structure tensor +:cite:`TFW13`. The updated stress, after the elastic +relaxation, is then passed to the momentum equation and the sea ice +velocities are updated in the usual manner within the subcycling loop of +the EVP rheology. The structure tensor evolution equations are solved +implicitly at the same frequency, :math:`\Delta t_{e}`, as the ice +velocities and internal stresses. Finally, to be coherent with our new +rheology we compute the area loss rate due to ridging as +:math:`\vert\dot{\boldsymbol\epsilon}\vert\alpha_{r}(\theta)`, with +:math:`\alpha_r(\theta)` and :math:`\alpha_s(\theta)` given by +:cite:`WF04`, .. math:: - P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} - \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} - \left( \frac{(H_n^{\max})^3 - (H_n^{\min})^3} - {3(H_n^{\max}-H_n^{\min})} \right) \right], - :label: roth-strength0 + \begin{aligned} + \alpha_{r}(\theta)=\frac{\sigma^{r}_{ij}\dot\epsilon_{ij}}{P_{r} \vert\dot{\boldsymbol\epsilon}\vert } , \qquad \alpha_{s}(\theta)=\frac{\sigma^{s}_{ij}\dot\epsilon_{ij}}{P_{s} \vert\dot{\boldsymbol\epsilon}\vert }.\label{alphas}\end{aligned} -where :math:`C_P = (g/2)(\rho_i/\rho_w)(\rho_w-\rho_i)`, -:math:`\beta =R_{\mathrm{tot}}/R_{\mathrm{net}} > 1` -from Equation :eq:`Rtot-Rnet`, and :math:`C_f` is an empirical parameter that -accounts for frictional energy dissipation. Following -:cite:`FH95`, we set :math:`C_f = 17`. The first term in -the summation is the potential energy of ridging ice, and the second, -larger term is the potential energy of the resulting ridges. The factor -of :math:`\beta` is included because :math:`a_{Pn}` is normalized with -respect to the total area of ice ridging, not the net area removed. -Recall that more than one unit area of ice must be ridged to reduce the -net ice area by one unit. For exponential ridge ITDs (`krdg\_redist` = 1), -the ridge potential energy is modified: +Both ridging rate and sea ice strength are computed in the outer loop +of the dynamics. -.. math:: - P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} - \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} - \left( H_{\min}^2 + 2H_{\min}\lambda + 2 \lambda^2 \right) \right] % CHECK BRACES - :label: roth-strength1 +.. _revp: -The energy-based ice strength given by Equations :eq:`roth-strength0` or -:eq:`roth-strength1` is more physically realistic than the strength -given by Equation :eq:`hib-strength`. However, use of Equation :eq:`hib-strength` is -less likely to allow numerical instability at a given resolution and -time step. See :cite:`LHMJ07` for more details. +**************** +Revised approach +**************** -.. _dynam: +A modification of the standard elastic-viscous-plastic (EVP) approach +for sea ice dynamics has been proposed by :cite:`BFLM13`, +that generalizes the EVP elastic modulus :math:`E` and the time +stepping approach for both momentum and stress to use an +under-relaxation technique. In general terms, the momentum and stress +equations become -Dynamics --------- +.. math:: + \begin{aligned} + {\bf u}^{k+1} &=& {\bf u}^k + \left(\breve{{\bf u}}^k - {\bf u}^{k+1}\right){1\over\beta} \\ + \sigma^{k+1} &=& \sigma^k + \left(\breve{\sigma}^k - \sigma^{k+1}\right){1\over\alpha} \end{aligned} -There are now different rheologies available in the CICE code. The -elastic-viscous-plastic (EVP) model represents a modification of the -standard viscous-plastic (VP) model for sea ice dynamics -:cite:`Hibler79`. The elastic-anisotropic-plastic (EAP) model, -on the other hand, explicitly accounts for the observed sub-continuum -anisotropy of the sea ice cover :cite:`WF06,WS09`. If -`kdyn` = 1 in the namelist then the EVP rheology is used (module -**ice\_dyn\_evp.F90**), while `kdyn` = 2 is associated with the EAP -rheology (**ice\_dyn\_eap.F90**). At times scales associated with the -wind forcing, the EVP model reduces to the VP model while the EAP model -reduces to the anisotropic rheology described in detail in -:cite:`WF06,TFW13`. At shorter time scales the -adjustment process takes place in both models by a numerically more -efficient elastic wave mechanism. While retaining the essential physics, -this elastic wave modification leads to a fully explicit numerical -scheme which greatly improves the model’s computational efficiency. +where :math:`\breve{{\bf u}}` and :math:`\breve{\sigma}` represent +the converged VP solution and :math:`\alpha, \beta < 1`. -The EVP sea ice dynamics model is thoroughly documented in -:cite:`HD97`, :cite:`Hunke01`, -:cite:`HD02` and :cite:`HD03` and the EAP -dynamics in :cite:`TFW13`. Simulation results and -performance of the EVP and EAP models have been compared with the VP -model and with each other in realistic simulations of the Arctic -respectively in :cite:`HZ99` and -:cite:`TFW13`. Here we summarize the equations and -direct the reader to the above references for details. The numerical -implementation in this code release is that of :cite:`HD02` -and :cite:`HD03`, with revisions to the numerical solver as -in :cite:`BFLM13`. The implementation of the EAP sea ice -dynamics into CICE is described in detail in -:cite:`TFW13`. +*Momentum* -.. _momentum: +The momentum equations become -Momentum -~~~~~~~~ +.. math:: + \begin{aligned} + \beta{m\over\Delta t} \left(u^{k+1}-u^k\right) &=& \overline{u} + {\tt vrel}\left(-u^{k+1}\cos\theta + v^{k+1}\sin\theta\right) + mfv^{k+1} - {m\over \Delta t} u^{k+1} \\ + \beta{m\over\Delta t} \left(v^{k+1}-v^k\right) &=& \overline{v} - {\tt vrel}\left(u^{k+1}\sin\theta + v^{k+1}\cos\theta\right) - mfu^{k+1} - {m\over \Delta t} v^{k+1} \end{aligned} -The force balance per unit area in the ice pack is given by a -two-dimensional momentum equation :cite:`Hibler79`, obtained -by integrating the 3D equation through the thickness of the ice in the -vertical direction: +where .. math:: - m{\partial {\bf u}\over\partial t} = \nabla\cdot{\bf \sigma} - + \vec{\tau}_a+\vec{\tau}_w - \hat{k}\times mf{\bf u} - mg\nabla H_\circ, - :label: vpmom + \overline{u} = F_u + \tau_{ax} - mg{\partial H_\circ\over\partial x} + {\tt vrel} \left(U_w\cos\theta - V_w\sin\theta\right) + {m\over\Delta t}u^\circ + :label: revpuhat -where :math:`m` is the combined mass of ice and snow per unit area and -:math:`\vec{\tau}_a` and :math:`\vec{\tau}_w` are wind and ocean -stresses, respectively. The strength of the ice is represented by the -internal stress tensor :math:`\sigma_{ij}`, and the other two terms on -the right hand side are stresses due to Coriolis effects and the sea -surface slope. The parameterization for the wind and ice–ocean stress -terms must contain the ice concentration as a multiplicative factor to -be consistent with the formal theory of free drift in low ice -concentration regions. A careful explanation of the issue and its -continuum solution is provided in :cite:`HD03` and -:cite:CGHM04`. +.. math:: + \overline{v} = F_v + \tau_{ay} - mg{\partial H_\circ\over\partial y} + {\tt vrel} \left(U_w\sin\theta + V_w\cos\theta\right) + {m\over\Delta t}v^\circ, + :label: revpvhat -The momentum equation is discretized in time as follows, for the classic -EVP approach. First, for clarity, the two components of Equation :eq:`vpmom` are +:math:`{\bf u}^\circ` is the initial value of velocity at the +beginning of the subcycling (:math:`k=0`), and we use +:math:`{\bf u}^{k+1}` for the ice–ocean stress and Coriolis terms. +Equations :eq:`revpuhat` and :eq:`revpvhat` differ from +Equations :eq:`cevpuhat` and :eq:`cevpvhat` only in the last term. + +Solving simultaneously for :math:`{\bf u}^{k+1}` as before, we have .. math:: \begin{aligned} - m{\partial u\over\partial t} &=& {\partial\sigma_{1j}\over\partial x_j} + \tau_{ax} + - a_i c_w \rho_w - \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\cos\theta - \left(V_w-v\right)\sin\theta\right] - +mfv - mg{\partial H_\circ\over\partial x}, \\ - m{\partial v\over\partial t} &=& {\partial\sigma_{2j}\over\partial x_j} + \tau_{ay} + - a_i c_w \rho_w - \left|{\bf U}_w - {\bf u}\right| \left[\left(U_w-u\right)\sin\theta - \left(V_w-v\right)\cos\theta\right] - -mfu - mg{\partial H_\circ\over\partial y}. \end{aligned} + u^{k+1} = {\tilde{a} \tilde{u} + b \tilde{v} \over \tilde{a}^2 + b^2} \\ + v^{k+1} = {\tilde{a} \tilde{v} - b \tilde{u} \over \tilde{a}^2 + b^2}, \end{aligned} -In the code, -:math:`{\tt vrel}=a_i c_w \rho_w\left|{\bf U}_w - {\bf u}^k\right|`, -where :math:`k` denotes the subcycling step. The following equations -illustrate the time discretization and define some of the other -variables used in the code. +where .. math:: - \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta\right)}_{\tt cca} u^{k+1} - - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb}v^{k+1} - = \underbrace{{\partial\sigma_{1j}^{k+1}\over\partial x_j}}_{\tt strintx} - + \underbrace{\tau_{ax} - mg{\partial H_\circ\over\partial x} }_{\tt forcex} - + {\tt vrel}\underbrace{\left(U_w\cos\theta-V_w\sin\theta\right)}_{\tt waterx} + {m\over\Delta t_e}u^k, - :label: umom + \tilde{a} = \left(1+\beta\right){m\over\Delta t} + {\tt vrel}\cos\theta \\ .. math:: - \underbrace{\left(mf+{\tt vrel}\sin\theta\right)}_{\tt ccb} u^{k+1} - + \underbrace{\left({m\over\Delta t_e}+{\tt vrel} \cos\theta\right)}_{\tt cca}v^{k+1} - = \underbrace{{\partial\sigma_{2j}^{k+1}\over\partial x_j}}_{\tt strinty} - + \underbrace{\tau_{ay} - mg{\partial H_\circ\over\partial y} }_{\tt forcey} - + {\tt vrel}\underbrace{\left(U_w\sin\theta+V_w\cos\theta\right)}_{\tt watery} + {m\over\Delta t_e}v^k, - :label: vmom + \tilde{\bf u} = \overline{\bf u} + \beta {m\over\Delta t}{\bf u}^k, + :label: tildeu -and vrel\ :math:`\cdot`\ waterx(y) = taux(y). +and :math:`b` is the same as in Equation :eq:`cevpb`. -We solve this system of equations analytically for :math:`u^{k+1}` and -:math:`v^{k+1}`. Define +*Stress* -.. math:: - \hat{u} = F_u + \tau_{ax} - mg{\partial H_\circ\over\partial x} + {\tt vrel} \left(U_w\cos\theta - V_w\sin\theta\right) + {m\over\Delta t_e}u^k - :label: cevpuhat +In CICE’s classic approach, the update to :math:`\sigma_1` at subcycle +step :math:`k+1` is .. math:: - \hat{v} = F_v + \tau_{ay} - mg{\partial H_\circ\over\partial y} + {\tt vrel} \left(U_w\sin\theta + V_w\cos\theta\right) + {m\over\Delta t_e}v^k, - :label: cevpvhat + \sigma_1^{k+1} + = \left(\sigma_1^{k} + {P\over\Delta}{\Delta t_e\over 2T} \left(\dot{\epsilon} - \Delta\right)\right) * \left(1 + {\Delta t_e\over 2T}\right) + :label: sig1time -where :math:`{\bf F} = \nabla\cdot\sigma^{k+1}`. Then +If we set -.. math:: - \begin{aligned} - \left({m\over\Delta t_e} +{\tt vrel}\cos\theta\right)u^{k+1} - \left(mf + {\tt vrel}\sin\theta\right) v^{k+1} &=& \hat{u} \\ - \left(mf + {\tt vrel}\sin\theta\right) u^{k+1} + \left({m\over\Delta t_e} +{\tt vrel}\cos\theta\right)v^{k+1} &=& \hat{v}.\end{aligned} +.. math:: + \alpha_1 = {2T\over \Delta t_e}, -Solving simultaneously for :math:`u^{k+1}` and :math:`v^{k+1}`, +then Equation :eq:`sig1time` becomes -.. math:: - \begin{aligned} - u^{k+1} = {a \hat{u} + b \hat{v} \over a^2 + b^2} \\ - v^{k+1} = {a \hat{v} - b \hat{u} \over a^2 + b^2}, \end{aligned} +.. math:: + \sigma_1^{k+1}\left(1+\alpha_1\right) = \alpha_1\sigma_1^k + {P\over\Delta} \left(\dot{\epsilon} - \Delta\right). -where +This is equivalent to Eq. (23) in :cite:`BFLM13`, but +using :math:`\sigma` at the current subcycle :math:`k+1` in the last +term on the right-hand side. Likewise, setting -.. math:: - a = {m\over\Delta t_e} + {\tt vrel}\cos\theta \\ - :label: cevpa +.. math:: + \alpha_2 = {2T\over e^2\Delta t_e} = {\alpha_1\over e^2} + +produces equations equivalent to Eq. (23) in +:cite:`BFLM13` for :math:`\sigma_2` and +:math:`\sigma_{12}`. Therefore the only change needed in the stress +code is to use :math:`\alpha_1` and :math:`\alpha_2` instead of +:math:`2T / \Delta t_e` and :math:`2T /e^2 \Delta t_e`. + +However, :cite:`BFLM13` introduce another change to the EVP +stress equations by altering the form of Young’s modulus in the elastic +term: the coefficient of :math:`\partial\sigma_1/\partial t` is +:math:`1/E`, but it is :math:`e^2/E` in the :math:`\sigma_2` and +:math:`\sigma_{12}` equations. This change does not affect the VP +equations to which the EVP equations should converge, but it does affect +the transient path taken during the subcycling. Since EVP subcycling is +finite, the numerical solutions obtained using this method differ from +the original EVP code. + +To implement this second change, we need define only +:math:`\alpha_1 = {2T/\Delta t_e}` as above and incorporate the factor +of :math:`e^2` from :math:`\alpha_2` into the equations for +:math:`\sigma_2` and :math:`\sigma_{12}`: .. math:: - b = mf + {\tt vrel}\sin\theta. - :label: cevpb + \begin{aligned} + \sigma_1^{k+1}\left(1+\alpha_1\right) &=&\sigma_1^k + {\alpha_1}{P\over\Delta} D_D, \\ + \sigma_2^{k+1}\left(1+\alpha_1\right) &=&\sigma_2^k + {\alpha_1\over e^2}{P\over\Delta} D_T, \\ + \sigma_{12}^{k+1}\left(1+\alpha_1\right) &=&\sigma_{12}^k + {\alpha_1\over 2e^2}{P\over\Delta} D_S.\end{aligned} -When the subcycling is finished for each (thermodynamic) time step, the -ice–ocean stress must be constructed from `taux(y)` and the terms -containing `vrel` on the left hand side of the equations. +To minimize code changes and unify the two approaches, we define and +apply :math:`1/\alpha_1` and :math:`\beta` in the classic EVP code, and +modify the elastic stress term. These under-relaxation parameters +control the rate at which the iteration converges. Thus for classic EVP +we set -The Hibler-Bryan form for the ice-ocean stress :cite:`HB87` -is included in **ice\_dyn\_shared.F90** but is currently commented out, -pending further testing. +.. math:: + \begin{aligned} + {\tt arlx1i} &=& {1\over\alpha_1} = {\Delta t_e\over 2T} \\ + {\tt brlx} &=& \beta = {\Delta t\over\Delta t_e}. \end{aligned} -.. _internal-stress: +Then -Internal stress -~~~~~~~~~~~~~~~ +.. math:: + \begin{aligned} + {\tt denom1} &=& {1\over{1+{\tt arlx1i}}} = {1\over{1+1/\alpha_1}} = {1\over{1+\Delta t_e/ 2T}} \\ + {\tt c1} &=& {P\over\Delta}\,{\tt arlx1i} = {P\over\Delta}{\Delta t_e\over 2T} \\ + {\tt c0} &=& {{\tt c1}\over e^2} = {P\over\Delta}{\Delta t_e\over 2Te^2} .\end{aligned} -For convenience we formulate the stress tensor :math:`\bf \sigma` in -terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, -:math:`\sigma_2=\sigma_{11}-\sigma_{22}`, and introduce the -divergence, :math:`D_D`, and the horizontal tension and shearing -strain rates, :math:`D_T` and :math:`D_S` respectively. +The stress equations for `stressp` (:math:`\sigma_1`) are unchanged; the +modified equations for `stressm` (:math:`\sigma_2`) and `stress12` +(:math:`\sigma_{12}`) take the form -*Elastic-Viscous-Plastic* +.. math:: + \begin{aligned} + {\tt stressm} &=& {\tt stressm + c0}\,D_T \,{\tt denom1}\\ + {\tt stress12} &=& {\tt stress12 + 0.5\,c0}\,D_S \,{\tt denom1}.\end{aligned} -In the EVP model the internal stress tensor is determined from a -regularized version of the VP constitutive law, +For classic EVP, -.. math:: - {1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} - + {P\over 2\zeta} = D_D, \\ - :label: sig1 +.. math:: + {\tt cca} = a = {\tt brlx}\,{m\over\Delta t} + {\tt vrel}\cos\theta ={m\over\Delta t_e} + {\tt vrel}\cos\theta. -.. math:: - {1\over E}{\partial\sigma_2\over\partial t} + {\sigma_2\over 2\eta} = D_T, - :label: sig2 +For revised EVP, arlx1i and brlx are defined separately from +:math:`\Delta t`, :math:`\Delta t_e`, :math:`T` and :math:`e`, and -.. math:: - {1\over E}{\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over - 2\eta} = {1\over 2}D_S, - :label: sig12 +.. math:: + {\tt cca} = \tilde{a} = \left(1+ {\tt brlx}\right){m\over\Delta t} + {\tt vrel}\cos\theta= \left(1+\beta\right){m\over\Delta t} + {\tt vrel}\cos\theta. -where +:math:`\tilde{\bf u}` must also be defined for revised EVP as in +Equation :eq:`tildeu`. The extra terms in :math:`\tilde{a}` and +:math:`\tilde{\bf u}` are multiplied by a flag (revp) that equals 1 for +revised EVP and 0 for classic EVP. Revised EVP is activated by setting +the namelist parameter `revised\_evp` = true. Note that in the current +implementation, only the modified version of the elastic term is +available for either the classic (`revised\_evp` = false) or the revised +EVP method. A final difference is that the revised approach initializes +the stresses to 0 at the beginning of each time step, while the classic +EVP approach uses the previous time step value. -.. math:: - D_D = \dot{\epsilon}_{11} + \dot{\epsilon}_{22}, -.. math:: - D_T = \dot{\epsilon}_{11} - \dot{\epsilon}_{22}, +~~~~~~~~~~~~~~~~~~~~ +Thickness changes +~~~~~~~~~~~~~~~~~~~~ -.. math:: - D_S = 2\dot{\epsilon}_{12}, +.. _itd-trans: -.. math:: - \dot{\epsilon}_{ij} = {1\over 2}\left({{\partial u_i}\over{\partial x_j}} + {{\partial u_j}\over{\partial x_i}}\right), +**************************** +Transport in thickness space +**************************** -.. math:: - \zeta = {P\over 2\Delta}, +Next we solve the equation for ice transport in thickness space due to +thermodynamic growth and melt, .. math:: - \eta = {P\over {2\Delta e^2}}, + \frac{\partial g}{\partial t} + \frac{\partial}{\partial h} (f g) = 0, + :label: itd-transport -.. math:: - \Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2}, +which is obtained from Equation :eq:`transport-g` by neglecting the first and +third terms on the right-hand side. We use the remapping method of +:cite:`Lipscomb01`, in which thickness categories are +represented as Lagrangian grid cells whose boundaries are projected +forward in time. The thickness distribution function :math:`g` is +approximated as a linear function of :math:`h` in each displaced +category and is then remapped onto the original thickness categories. +This method is numerically smooth and is not too diffusive. It can be +viewed as a 1D simplification of the 2D incremental remapping scheme +described above. -and :math:`P` is a function of the ice thickness and concentration, -described in Section :ref:`mech-red`. The dynamics component -employs a “replacement pressure” (see :cite:`GHA98`, for -example), which serves to prevent residual ice motion due to spatial -variations of :math:`P` when the rates of strain are exactly zero. +We first compute the displacement of category boundaries in thickness +space. Assume that at time :math:`m` the ice areas :math:`a_n^m` and +mean ice thicknesses :math:`h_n^m` are known for each thickness +category. (For now we omit the subscript :math:`i` that distinguishes +ice from snow.) We use a thermodynamic model (Section :ref:`thermo`) +to compute the new mean thicknesses :math:`h_n^{m+1}` at time +:math:`m+1`. The time step must be small enough that trajectories do not +cross; i.e., :math:`h_n^{m+1} < h_{n+1}^{m+1}` for each pair of adjacent +categories. The growth rate at :math:`h = h_n` is given by +:math:`f_n = (h_n^{m+1} - h_n^m) / \Delta t`. By linear interpolation we +estimate the growth rate :math:`F_n` at the upper category boundary +:math:`H_n`: -Viscosities are updated during the subcycling, so that the entire -dynamics component is subcycled within the time step, and the elastic -parameter :math:`E` is defined in terms of a damping timescale :math:`T` -for elastic waves, :math:`\Delta t_e < T < \Delta t`, as +.. math:: + F_n = f_n + \frac{f_{n+1}-f_n}{h_{n+1}-h_n} \, (H_n - h_n). + +If :math:`a_n` or :math:`a_{n+1} = 0`, :math:`F_n` is set to the growth +rate in the nonzero category, and if :math:`a_n = a_{n+1} = 0`, we set +:math:`F_n = 0`. The temporary displaced boundaries are given by .. math:: - E = {\zeta\over T}, + H_n^* = H_n + F_n \, \Delta t, \ n = 1 \ {\rm to} \ N-1 -where :math:`T=E_\circ\Delta t` and :math:`E_\circ` (eyc) is a tunable -parameter less than one. The stress equations :eq:`sig1`–:eq:`sig12` -become +The boundaries must not be displaced by more than one category to the +left or right; that is, we require :math:`H_{n-1} < H_n^* < H_{n+1}`. +Without this requirement we would need to do a general remapping rather +than an incremental remapping, at the cost of added complexity. -.. math:: - \begin{aligned} - {\partial\sigma_1\over\partial t} + {\sigma_1\over 2T} - + {P\over 2T} &=& {P\over 2T\Delta} D_D, \\ - {\partial\sigma_2\over\partial t} + {e^2\sigma_2\over 2T} &=& {P\over - 2T\Delta} D_T,\\ - {\partial\sigma_{12}\over\partial t} + {e^2\sigma_{12}\over 2T} &=& - {P\over 4T\Delta}D_S.\end{aligned} +Next we construct :math:`g(h)` in the displaced thickness categories. +The ice areas in the displaced categories are :math:`a_n^{m+1} = a_n^m`, +since area is conserved following the motion in thickness space (i.e., +during vertical ice growth or melting). The new ice volumes are +:math:`v_n^{m+1} = (a_n h_n)^{m+1} = a_n^m h_n^{m+1}`. For conciseness, +define :math:`H_L = H_{n-1}^*` and :math:`H_R = H_{n}^*` and drop the +time index :math:`m+1`. We wish to construct a continuous function +:math:`g(h)` within each category such that the total area and volume at +time :math:`m+1` are :math:`a_n` and :math:`v_n`, respectively: -All coefficients on the left-hand side are constant except for -:math:`P`, which changes only on the longer time step :math:`\Delta t`. -This modification compensates for the decreased efficiency of including -the viscosity terms in the subcycling. (Note that the viscosities do not -appear explicitly.) Choices of the parameters used to define :math:`E`, -:math:`T` and :math:`\Delta t_e` are discussed in -Sections :ref:`revp` and :ref:`parameters`. +.. math:: + \int_{H_L}^{H_R} g \, dh = a_n, + :label: area-cons -The bilinear discretization used for the stress terms -:math:`\partial\sigma_{ij}/\partial x_j` in the momentum equation is -now used, which enabled the discrete equations to be derived from the -continuous equations written in curvilinear coordinates. In this -manner, metric terms associated with the curvature of the grid are -incorporated into the discretization explicitly. Details pertaining to -the spatial discretization are found in :cite:`HD02`. +.. math:: + \int_{H_L}^{H_R} h \, g \, dh = v_n. + :label: volume-cons -*Elastic-Anisotropic-Plastic* +The simplest polynomial that can satisfy both equations is a line. It +is convenient to change coordinates, writing +:math:`g(\eta) = g_1 \eta + g_0`, where :math:`\eta = h - H_L` and the +coefficients :math:`g_0` and :math:`g_1` are to be determined. Then +Equations :eq:`area-cons` and :eq:`volume-cons` can be written as -In the EAP model the internal stress tensor is related to the -geometrical properties and orientation of underlying virtual diamond -shaped floes (see :ref:`fig-EAP`). In contrast to the isotropic EVP -rheology, the anisotropic plastic yield curve within the EAP rheology -depends on the relative orientation of the diamond shaped floes (unit -vector :math:`\mathbf r` in :ref:`fig-EAP`), with respect to the -principal direction of the deformation rate (not shown). Local -anisotropy of the sea ice cover is accounted for by an additional -prognostic variable, the structure tensor :math:`\mathbf{A}` defined -by +.. math:: + g_1 \frac{\eta_R^2}{2} + g_0 \eta_R = a_n, .. math:: - {\mathbf A}=\int_{\mathbb{S}}\vartheta(\mathbf r)\mathbf r\mathbf r d\mathbf r\label{structuretensor}. + g_1 \frac{\eta_R^3}{3} + g_0 \frac{\eta_R^2}{2} = a_n \eta_n, -where :math:`\mathbb{S}` is a unit-radius circle; **A** is a unit -trace, 2\ :math:`\times`\ 2 matrix. From now on we shall describe the -orientational distribution of floes using the structure tensor. For -simplicity we take the probability density function -:math:`\vartheta(\mathbf r )` to be Gaussian, -:math:`\vartheta(z)=\omega_{1}\exp(-\omega_{2}z^{2})`, where :math:`z` -is the ice floe inclination with respect to the axis :math:`x_{1}` of -preferential alignment of ice floes (see :ref:`fig-EAP`), -:math:`\vartheta(z)` is periodic with period :math:`\pi`, and the -positive coefficients :math:`\omega_{1}` and :math:`\omega_{2}` are -calculated to ensure normalization of :math:`\vartheta(z)`, i.e. -:math:`\int_{0}^{2\pi}\vartheta(z)dz=1`. The ratio of the principal -components of :math:`\mathbf{A}`, :math:`A_{1}/A_{2}`, are derived -from the phenomenological evolution equation for the structure tensor -:math:`\mathbf A`, +where :math:`\eta_R = H_R - H_L` and :math:`\eta_n = h_n - H_L`. These +equations have the solution -.. math:: - \frac{D\mathbf{A}}{D t}=\mathbf{F}_{iso}(\mathbf{A})+\mathbf{F}_{frac}(\mathbf{A},\boldsymbol\sigma), - :label: evolutionA +.. math:: + g_0 = \frac{6 a_n}{\eta_R^2} \left(\frac{2 \eta_R}{3} - \eta_n\right), + :label: g0 -where :math:`t` is the time, and :math:`D/Dt` is the co-rotational -time derivative accounting for advection and rigid body rotation -(:math:`D\mathbf A/Dt = d\mathbf A/dt -\mathbf W \cdot \mathbf A -\mathbf A \cdot \mathbf W^{T}`) -with :math:`\mathbf W` being the vorticity tensor. -:math:`\mathbf F_{iso}` is a function that accounts for a variety of -processes (thermal cracking, melting, freezing together of floes) that -contribute to a more isotropic nature to the ice cover. -:math:`\mathbf F_{frac}` is a function determining the ice floe -re-orientation due to fracture, and explicitly depends upon sea ice -stress (but not its magnitude). Following :cite:`WF06`, -based on laboratory experiments by :cite:`Schulson01` we -consider four failure mechanisms for the Arctic sea ice cover. These -are determined by the ratio of the principal values of the sea ice -stress :math:`\sigma_{1}` and :math:`\sigma_{2}`: (i) under biaxial -tension, fractures form across the perpendicular principal axes and -therefore counteract any apparent redistribution of the floe -orientation; (ii) if only one of the principal stresses is -compressive, failure occurs through axial splitting along the -compression direction; (iii) under biaxial compression with a low -confinement ratio, (:math:`\sigma_{1}/\sigma_{2} 2 \eta_R/3`, i.e., when :math:`h_n` lies in the right +third of the thickness range :math:`(H_L, H_R)`. Similarly, Equation :eq:`gmax` +implies that :math:`g(\eta_R) < 0` when :math:`\eta_n < \eta_R/3`, i.e., +when :math:`h_n` is in the left third of the range. Since negative +values of :math:`g` are unphysical, a different solution is needed when +:math:`h_n` lies outside the central third of the thickness range. If +:math:`h_n` is in the left third of the range, we define a cutoff +thickness, :math:`H_C = 3 h_n - 2 H_L`, and set :math:`g = 0` between +:math:`H_C` and :math:`H_R`. Equations :eq:`g0` and :eq:`g1` are then +valid with :math:`\eta_R` redefined as :math:`H_C - H_L`. And if +:math:`h_n` is in the right third of the range, we define +:math:`H_C = 3 h_n - 2 H_R` and set :math:`g = 0` between :math:`H_L` +and :math:`H_C`. In this case, :eq:`g0` and :eq:`g1` apply with +:math:`\eta_R = H_R - H_C` and :math:`\eta_n = h_n - H_C`. -:ref:`fig-EAP` : Geometry of interlocking diamond-shaped floes (taken from -:cite:`WF06`). :math:`\phi` is half of the acute angle -of the diamonds. :math:`L` is the edge length. -:math:`\boldsymbol n_{1}`, :math:`\boldsymbol n_{2}` and -:math:`\boldsymbol\tau_{1}`, :math:`\boldsymbol\tau_{2}` are -respectively the normal and tangential unit vectors along the diamond edges. -:math:`\mathbf v=L\boldsymbol\tau_{2}\cdot\dot{\boldsymbol\epsilon}` -is the relative velocity between the two floes connected by the -vector :math:`L \boldsymbol \tau_{2}`. :math:`\mathbf r` is the unit -vector along the main diagonal of the diamond. Note that the diamonds -illustrated here represent one possible realisation of all possible -orientations. The angle :math:`z` represents the rotation of the -diamonds’ main axis relative to their preferential orientation along -the axis :math:`x_1`. +:ref:`fig-gplot` illustrates the linear reconstruction of :math:`g` +for the simple cases :math:`H_L = 0`, :math:`H_R = 1`, :math:`a_n = 1`, +and :math:`h_n =` 0.2, 0.4, 0.6, and 0.8. Note that :math:`g` slopes +downward (:math:`g_1 < 0`) when :math:`h_n` is less than the midpoint +thickness, :math:`(H_L + H_R)/2 = 1/2`, and upward when :math:`h_n` +exceeds the midpoint thickness. For :math:`h_n = 0.2` and 0.8, +:math:`g = 0` over part of the range. -The new anisotropic rheology requires solving the evolution -Equation :eq:`evolutionA` for the structure tensor in addition to the momentum -and stress equations. The evolution equation for :math:`\mathbf{A}` is -solved within the EVP subcycling loop, and consistently with the -momentum and stress evolution equations, we neglect the advection term -for the structure tensor. Equation :eq:`evolutionA` then reduces to the system -of two equations: +.. _fig-gplot: -.. math:: - \begin{aligned} - \frac{\partial A_{11}}{\partial t}&=&-k_{t}\left(A_{11}-\frac{1}{2}\right)+M_{11} \mbox{,} \\ - \frac{\partial A_{12}}{\partial t}&=&-k_{t} A_{12}+M_{12} \mbox{,}\end{aligned} +.. figure:: ./figures/gplot.png + :align: center + :scale: 20% -where the first terms on the right hand side correspond to the -isotropic contribution, :math:`F_{iso}`, and :math:`M_{11}` and -:math:`M_{12}` are the components of the term :math:`F_{frac}` in -Equation :eq:`evolutionA` that are given in :cite:`WF06` and -:cite:`TFW13`. These evolution equations are -discretized semi-implicitly in time. The degree of anisotropy is -measured by the largest eigenvalue (:math:`A_{1}`) of this tensor -(:math:`A_{2}=1-A_{1}`). :math:`A_{1}=1` corresponds to perfectly -aligned floes and :math:`A_{1}=0.5` to a uniform distribution of floe -orientation. Note that while we have specified the aspect ratio of the -diamond floes, through prescribing :math:`\phi`, we make no assumption -about the size of the diamonds so that formally the theory is scale -invariant. + Figure 4 -As described in greater detail in :cite:`WF06`, the -internal ice stress for a single orientation of the ice floes can be -calculated explicitly and decomposed, for an average ice thickness -:math:`h`, into its ridging (r) and sliding (s) contributions +:ref:`fig-gplot` : Linear approximation of the thickness distribution +function :math:`g(h)` for an ice category with left boundary +:math:`H_L = 0`, right boundary :math:`H_R = 1`, fractional area +:math:`a_n = 1`, and mean ice thickness :math:`h_n = 0.2, 0.4, 0.6,` and :math:`0.8`. + +Finally, we remap the thickness distribution to the original boundaries +by transferring area and volume between categories. We compute the ice +area :math:`\Delta a_n` and volume :math:`\Delta v_n` between each +original boundary :math:`H_n` and displaced boundary :math:`H_n^*`. If +:math:`H_n^* > H_n`, ice moves from category :math:`n` to :math:`n+1`. +The area and volume transferred are .. math:: - \boldsymbol \sigma^{b}(\mathbf r,h)=P_{r}(h) \boldsymbol \sigma_{r}^{b}(\mathbf r)+P_{s}(h) \boldsymbol \sigma_{s}^{b}(\mathbf r), - :label: stress1 + \Delta a_n = \int_{H_n}^{H_n^*} g \, dh, + :label: move-area -where :math:`P_{r}` and :math:`P_{s}` are the ridging and sliding -strengths and the ridging and sliding stresses are functions of the -angle :math:`\theta= \arctan(\dot\epsilon_{II}/\dot\epsilon_{I})`, the -angle :math:`y` between the major principal axis of the strain rate -tensor (not shown) and the structure tensor (:math:`x_1` axis in -:ref:`fig-EAP`, and the angle :math:`z` defined in :ref:`fig-EAP`. In -the stress expressions above the underlying floes are assumed parallel, -but in a continuum-scale sea ice region the floes can possess different -orientations in different places and we take the mean sea ice stress -over a collection of floes to be given by the average +.. math:: + \Delta v_n = \int_{H_n}^{H_n^*} h \, g \, dh. + :label: move-volume + +If :math:`H_n^* < H_N`, ice area and volume are transferred from +category :math:`n+1` to :math:`n` using Equations :eq:`move-area` and +:eq:`move-volume` with the limits of integration reversed. To evaluate +the integrals we change coordinates from :math:`h` to +:math:`\eta = h - H_L`, where :math:`H_L` is the left limit of the range +over which :math:`g > 0`, and write :math:`g(\eta)` using Equations :eq:`g0` and +:eq:`g1`. In this way we obtain the new areas :math:`a_n` and volumes +:math:`v_n` between the original boundaries :math:`H_{n-1}` and +:math:`H_n` in each category. The new thicknesses, +:math:`h_n = v_n/a_n`, are guaranteed to lie in the range +:math:`(H_{n-1}, H_n)`. If :math:`g = 0` in the part of a category that +is remapped to a neighboring category, no ice is transferred. + +Other conserved quantities are transferred in proportion to the ice +volume :math:`\Delta v_{in}`. For example, the transferred ice energy in +layer :math:`k` is +:math:`\Delta e_{ink} = e_{ink} (\Delta v_{in} / v_{in})`. + +The left and right boundaries of the domain require special treatment. +If ice is growing in open water at a rate :math:`F_0`, the left boundary +:math:`H_0` is shifted to the right by :math:`F_0 \Delta t` before +:math:`g` is constructed in category 1, then reset to zero after the +remapping is complete. New ice is then added to the grid cell, +conserving area, volume, and energy. If ice cannot grow in open water +(because the ocean is too warm or the net surface energy flux is +downward), :math:`H_0` is fixed at zero, and the growth rate at the left +boundary is estimated as :math:`F_0 = f_1`. If :math:`F_0 < 0`, all ice +thinner than :math:`\Delta h_0 = -F_0 \Delta t` is assumed to have +melted, and the ice area in category 1 is reduced accordingly. The area +of new open water is .. math:: - \boldsymbol\sigma^{EAP}(h)=P_{r}(h)\int_{\mathbb{S}}\vartheta(\mathbf r)\left[\boldsymbol\sigma_{r}^{b}(\mathbf r)+ k \boldsymbol\sigma_{s}^{b}(\mathbf r)\right]d\mathbf r - :label: stressaverage + \Delta a_0 = \int_{0}^{\Delta h_0} g \, dh. -where we have introduced the friction parameter :math:`k=P_{s}/P_{r}` -and where we identify the ridging ice strength :math:`P_{r}(h)` with the -strength :math:`P` described in section 1 and used within the EVP -framework. +The right boundary :math:`H_N` is not fixed but varies with +:math:`h_N`, the mean ice thickness in the thickest category. Given +:math:`h_N`, we set :math:`H_N = 3 h_N - 2 H_{N-1}`, which ensures that +:math:`g(h) > 0` for :math:`H_{N-1} < h < H_N` and :math:`g(h) = 0` for +:math:`h \geq H_N`. No ice crosses the right boundary. If the ice growth +or melt rates in a given grid cell are too large, the thickness +remapping scheme will not work. Instead, the thickness categories in +that grid cell are treated as delta functions following +:cite:`BHWE01`, and categories outside their prescribed +boundaries are merged with neighboring categories as needed. For time +steps of less than a day and category thickness ranges of 10 cm or more, +this simplification is needed rarely, if ever. -As is the case for the EVP rheology, elasticity is included in the EAP -description not to describe any physical effect, but to make use of the -efficient, explicit numerical algorithm used to solve the full sea ice -momentum balance. We use the analogous EAP stress equations, +The linear remapping algorithm for thickness is not monotonic for +tracers, although significant errors rarely occur. Usually they appear +as snow temperatures (enthalpy) outside the physical range of values in +very small snow volumes. In this case we transfer the snow and its heat +and tracer contents to the ocean. -.. math:: - \frac{\partial \sigma_{1}}{\partial t}+\frac{\sigma_1}{2T} = \frac{\sigma^{EAP}_{1}}{2T} \mbox{,} - :label: EAPsigma1 +.. _mech-red: -.. math:: - \frac{\partial \sigma_{2}}{\partial t}+\frac{\sigma_2}{2T} = \frac{\sigma^{EAP}_{2}}{2T} \mbox{,} - :label: EAPsigma2 +************************* +Mechanical redistribution +************************* + +The last term on the right-hand side of Equation :eq:`transport-g` +is :math:`\psi`, which describes the redistribution +of ice in thickness space due to ridging and other mechanical processes. +The mechanical redistribution scheme in CICE is based on +:cite:`TRMC75`, :cite:`Rothrock75`, +:cite:`Hibler80`, :cite:`FH95`, and +:cite:`LHMJ07`. This scheme converts thinner ice to thicker +ice and is applied after horizontal transport. When the ice is +converging, enough ice ridges to ensure that the ice area does not +exceed the grid cell area. + +First we specify the participation function: the thickness distribution +:math:`a_P(h) = b(h) \, g(h)` of the ice participating in ridging. (We +use “ridging” as shorthand for all forms of mechanical redistribution, +including rafting.) The weighting function :math:`b(h)` favors ridging +of thin ice and closing of open water in preference to ridging of +thicker ice. There are two options for the form of :math:`b(h)`. If +`krdg\_partic` = 0 in the namelist, we follow :cite:`TRMC75` +and set .. math:: - \frac{\partial \sigma_{12}}{\partial t}+\frac{\sigma_{12}}{2T} = \frac{\sigma^{EAP}_{12}}{2T} \mbox{,} - :label: EAPsigma12 + b(h) = \left\{\begin{array}{ll} + \frac{2}{G^*}(1-\frac{G(h)}{G^*}) & \mbox{if $G(h) G^*`, then :math:`a_{Pn} = 0`. If +the open water fraction :math:`a_0 > G^*`, no ice can ridge, because +“ridging” simply reduces the area of open water. As in +:cite:`TRMC75` we set :math:`G^* = 0.15`. -.. _revp: +If the spatial resolution is too fine for a given time step +:math:`\Delta t`, the weighting function Equation :eq:`partic-old-contin` can +promote numerical instability. For :math:`\Delta t = \mbox{1 hour}`, +resolutions finer than :math:`\Delta x \sim \mbox{10 km}` are typically +unstable. The instability results from feedback between the ridging +scheme and the dynamics via the ice strength. If the strength changes +significantly on time scales less than :math:`\Delta t`, the +viscous-plastic solution of the momentum equation is inaccurate and +sometimes oscillatory. As a result, the fields of ice area, thickness, +velocity, strength, divergence, and shear can become noisy and +unphysical. -Revised approach -~~~~~~~~~~~~~~~~ +A more stable weighting function was suggested by +:cite:`LHMJ07`: -A modification of the standard elastic-viscous-plastic (EVP) approach -for sea ice dynamics has been proposed by :cite:`BFLM13`, -that generalizes the EVP elastic modulus :math:`E` and the time -stepping approach for both momentum and stress to use an -under-relaxation technique. In general terms, the momentum and stress -equations become +.. math:: + b(h) = \frac{\exp[-G(h)/a^*]} + {a^*[1-\exp(-1/a^*)]} + :label: partic-new-contin + +When integrated between category boundaries, Equation :eq:`partic-new-contin` +implies .. math:: - \begin{aligned} - {\bf u}^{k+1} &=& {\bf u}^k + \left(\breve{{\bf u}}^k - {\bf u}^{k+1}\right){1\over\beta} \\ - \sigma^{k+1} &=& \sigma^k + \left(\breve{\sigma}^k - \sigma^{k+1}\right){1\over\alpha} \end{aligned} + a_{Pn} = \frac {\exp(-G_{n-1}/a^*) - \exp(-G_{n}/a^*)} + {1 - \exp(-1/a^*)} + :label: partic-new-discrete -where :math:`\breve{{\bf u}}` and :math:`\breve{\sigma}` represent -the converged VP solution and :math:`\alpha, \beta < 1`. +This weighting function is used if `krdg\_partic` = 1 in the namelist. +From Equation :eq:`partic-new-contin`, the mean value of :math:`G` for ice +participating in ridging is :math:`a^*`, as compared to :math:`G^*/3` +for Equation :eq:`partic-old-contin`. For typical ice thickness distributions, +setting :math:`a^* = 0.05` with `krdg\_partic` = 1 gives participation +fractions similar to those given by :math:`G^* = 0.15` with `krdg\_partic` += 0. See :cite:`LHMJ07` for a detailed comparison of these +two participation functions. -*Momentum* +Thin ice is converted to thick, ridged ice in a way that reduces the +total ice area while conserving ice volume and internal energy. There +are two namelist options for redistributing ice among thickness +categories. If `krdg\_redist` = 0, ridging ice of thickness :math:`h_n` +forms ridges whose area is distributed uniformly between +:math:`H_{\min} = 2 h_n` and :math:`H_{\max} = 2 \sqrt{H^* h_n}`, as in +:cite:`Hibler80`. The default value of :math:`H^*` is 25 m, as +in earlier versions of CICE. Observations suggest that +:math:`H^* = 50` m gives a better fit to first-year ridges +:cite:`AMI04`, although the lower value may be appropriate +for multiyear ridges :cite:`FH95`. The ratio of the mean +ridge thickness to the thickness of ridging ice is +:math:`k_n = (H_{\min} + H_{\max}) / (2 h_n)`. If the area of category +:math:`n` is reduced by ridging at the rate :math:`r_n`, the area of +thicker categories grows simultaneously at the rate :math:`r_n/k_n`. +Thus the *net* rate of area loss due to ridging of ice in category +:math:`n` is :math:`r_n(1-1/k_n)`. -The momentum equations become +The ridged ice area and volume are apportioned among categories in the +thickness range :math:`(H_{\min}, H_{\max})`. The fraction of the new +ridge area in category :math:`m` is .. math:: - \begin{aligned} - \beta{m\over\Delta t} \left(u^{k+1}-u^k\right) &=& \overline{u} + {\tt vrel}\left(-u^{k+1}\cos\theta + v^{k+1}\sin\theta\right) + mfv^{k+1} - {m\over \Delta t} u^{k+1} \\ - \beta{m\over\Delta t} \left(v^{k+1}-v^k\right) &=& \overline{v} - {\tt vrel}\left(u^{k+1}\sin\theta + v^{k+1}\cos\theta\right) - mfu^{k+1} - {m\over \Delta t} v^{k+1} \end{aligned} + f_m^{\mathrm{area}} = \frac{H_R - H_L} + {H_{\max} - H_{\min}}, + :label: ridge-area-old -where +where :math:`H_L = \max(H_{m-1},H_{\min})` and +:math:`H_R= \min(H_m,H_{\max})`. The fraction of the ridge volume going +to category :math:`m` is .. math:: - \overline{u} = F_u + \tau_{ax} - mg{\partial H_\circ\over\partial x} + {\tt vrel} \left(U_w\cos\theta - V_w\sin\theta\right) + {m\over\Delta t}u^\circ - :label: revpuhat + f_m^{\mathrm{vol}} = \frac{(H_R)^2 - (H_L)^2} + {(H_{\max})^2 - (H_{\min})^2}. + :label: ridge-volume-old + +This uniform redistribution function tends to produce too little ice in +the 3–5 m range and too much ice thicker than 10 m +:cite:`AMI04`. Observations show that the ITD of ridges is +better approximated by a negative exponential. Setting `krdg\_redist` = 1 +gives ridges with an exponential ITD :cite:`LHMJ07`: .. math:: - \overline{v} = F_v + \tau_{ay} - mg{\partial H_\circ\over\partial y} + {\tt vrel} \left(U_w\sin\theta + V_w\cos\theta\right) + {m\over\Delta t}v^\circ, - :label: revpvhat + g_R(h) \propto \exp[-(h - H_{\min})/\lambda] + :label: redist-new -:math:`{\bf u}^\circ` is the initial value of velocity at the -beginning of the subcycling (:math:`k=0`), and we use -:math:`{\bf u}^{k+1}` for the ice–ocean stress and Coriolis terms. -Equations :eq:`revpuhat` and :eq:`revpvhat` differ from -Equations :eq:`cevpuhat` and :eq:`cevpvhat` only in the last term. +for :math:`h \ge H_{\min}`, with :math:`g_R(h) = 0` for +:math:`h < H_{\min}`. Here, :math:`\lambda` is an empirical *e*-folding +scale and :math:`H_{\min}=2h_n` (where :math:`h_n` is the thickness of +ridging ice). We assume that :math:`\lambda = \mu h_n^{1/2}`, where +:math:`\mu` (mu\_rdg) is a tunable parameter with units . Thus the mean +ridge thickness increases in proportion to :math:`h_n^{1/2}`, as in +:cite:`Hibler80`. The value :math:`\mu = 4.0`  gives +:math:`\lambda` in the range 1–4 m for most ridged ice. Ice strengths +with :math:`\mu = 4.0`  and `krdg\_redist` = 1 are roughly comparable to +the strengths with :math:`H^* = 50` m and `krdg\_redist` = 0. -Solving simultaneously for :math:`{\bf u}^{k+1}` as before, we have +From Equation :eq:`redist-new` it can be shown that the fractional area going +to category :math:`m` as a result of ridging is .. math:: - \begin{aligned} - u^{k+1} = {\tilde{a} \tilde{u} + b \tilde{v} \over \tilde{a}^2 + b^2} \\ - v^{k+1} = {\tilde{a} \tilde{v} - b \tilde{u} \over \tilde{a}^2 + b^2}, \end{aligned} - -where + f_m^{\mathrm{area}} = \exp[-(H_{m-1} - H_{\min}) / \lambda] + - \exp[-(H_m - H_{\min}) / \lambda]. + :label: ridge-area-new -.. math:: - \tilde{a} = \left(1+\beta\right){m\over\Delta t} + {\tt vrel}\cos\theta \\ +The fractional volume going to category :math:`m` is .. math:: - \tilde{\bf u} = \overline{\bf u} + \beta {m\over\Delta t}{\bf u}^k, - :label: tildeu + f_m^{\mathrm{vol}} = \frac{(H_{m-1}+\lambda) \exp[-(H_{m-1}-H_{\min})/\lambda] + - (H_m + \lambda) \exp[-(H_m - H_{\min}) / \lambda]} + {H_{min} + \lambda}. + :label: ridge-volume-new -and :math:`b` is the same as in Equation :eq:`cevpb`. +Equations :eq:`ridge-area-new` and :eq:`ridge-volume-new` replace +Equations :eq:`ridge-area-old` and :eq:`ridge-volume-old` when `krdg\_redist` += 1. -*Stress* +Internal ice energy is transferred between categories in proportion to +ice volume. Snow volume and internal energy are transferred in the same +way, except that a fraction of the snow may be deposited in the ocean +instead of added to the new ridge. -In CICE’s classic approach, the update to :math:`\sigma_1` at subcycle -step :math:`k+1` is +The net area removed by ridging and closing is a function of the strain +rates. Let :math:`R_{\mathrm{net}}` be the net rate of area loss for the +ice pack (i.e., the rate of open water area closing, plus the net rate +of ice area loss due to ridging). Following :cite:`FH95`, +:math:`R_{\mathrm{net}}` is given by .. math:: - \sigma_1^{k+1} - = \left(\sigma_1^{k} + {P\over\Delta}{\Delta t_e\over 2T} \left(\dot{\epsilon} - \Delta\right)\right) * \left(1 + {\Delta t_e\over 2T}\right) - :label: sig1time - -If we set - -.. math:: - \alpha_1 = {2T\over \Delta t_e}, - -then Equation :eq:`sig1time` becomes + R_{\mathrm{net}} = \frac{C_s}{2} + (\Delta - |D_D|) - \min(D_D,0), + :label: Rnet -.. math:: - \sigma_1^{k+1}\left(1+\alpha_1\right) = \alpha_1\sigma_1^k + {P\over\Delta} \left(\dot{\epsilon} - \Delta\right). +where :math:`C_s` is the fraction of shear dissipation energy that +contributes to ridge-building, :math:`D_D` is the divergence, and +:math:`\Delta` is a function of the divergence and shear. These strain +rates are computed by the dynamics scheme. The default value of +:math:`C_s` is 0.25. -This is equivalent to Eq. (23) in :cite:`BFLM13`, but -using :math:`\sigma` at the current subcycle :math:`k+1` in the last -term on the right-hand side. Likewise, setting +Next, define :math:`R_{\mathrm{tot}} = \sum_{n=0}^N r_n`. This rate is +related to :math:`R_{\mathrm{net}}` by -.. math:: - \alpha_2 = {2T\over e^2\Delta t_e} = {\alpha_1\over e^2} +.. math:: + R_{\mathrm{net}} = + \left[ a_{P0} + \sum_{n=1}^N a_{Pn}\left(1-{1\over k_n}\right)\right] + R_{\mathrm{tot}}. + :label: Rtot-Rnet -produces equations equivalent to Eq. (23) in -:cite:`BFLM13` for :math:`\sigma_2` and -:math:`\sigma_{12}`. Therefore the only change needed in the stress -code is to use :math:`\alpha_1` and :math:`\alpha_2` instead of -:math:`2T / \Delta t_e` and :math:`2T /e^2 \Delta t_e`. +Given :math:`R_{\mathrm{net}}` from Equation :eq:`Rnet`, we +use Equation :eq:`Rtot-Rnet` to compute :math:`R_{\mathrm{tot}}`. Then the area +ridged in category :math:`n` is given by :math:`a_{rn} = r_n \Delta t`, +where :math:`r_n = a_{Pn} R_{\mathrm{tot}}`. The area of new ridges is +:math:`a_{rn} / k_n`, and the volume of new ridges is :math:`a_{rn} h_n` +(since volume is conserved during ridging). We remove the ridging ice +from category :math:`n` and use Equations :eq:`ridge-area-old` +and :eq:`ridge-volume-old`: (or :eq:`ridge-area-new` and +:eq:`ridge-volume-new`) to redistribute the ice among thicker +categories. -However, :cite:`BFLM13` introduce another change to the EVP -stress equations by altering the form of Young’s modulus in the elastic -term: the coefficient of :math:`\partial\sigma_1/\partial t` is -:math:`1/E`, but it is :math:`e^2/E` in the :math:`\sigma_2` and -:math:`\sigma_{12}` equations. This change does not affect the VP -equations to which the EVP equations should converge, but it does affect -the transient path taken during the subcycling. Since EVP subcycling is -finite, the numerical solutions obtained using this method differ from -the original EVP code. +Occasionally the ridging rate in thickness category :math:`n` may be +large enough to ridge the entire area :math:`a_n` during a time interval +less than :math:`\Delta t`. In this case :math:`R_{\mathrm{tot}}` is +reduced to the value that exactly ridges an area :math:`a_n` during +:math:`\Delta t`. After each ridging iteration, the total fractional ice +area :math:`a_i` is computed. If :math:`a_i > 1`, the ridging is +repeated with a value of :math:`R_{\mathrm{net}}` sufficient to yield +:math:`a_i = 1`. -To implement this second change, we need define only -:math:`\alpha_1 = {2T/\Delta t_e}` as above and incorporate the factor -of :math:`e^2` from :math:`\alpha_2` into the equations for -:math:`\sigma_2` and :math:`\sigma_{12}`: +Two tracers for tracking the ridged ice area and volume are available. +The actual tracers are for level (undeformed) ice area (`alvl`) and volume +(`vlvl`), which are easier to implement for a couple of reasons: (1) ice +ridged in a given thickness category is spread out among the rest of the +categories, making it more difficult (and expensive) to track than the +level ice remaining behind in the original category; (2) previously +ridged ice may ridge again, so that simply adding a volume of freshly +ridged ice to the volume of previously ridged ice in a grid cell may be +inappropriate. Although the code currently only tracks level ice +internally, both level ice and ridged ice are offered as history output. +They are simply related: .. math:: \begin{aligned} - \sigma_1^{k+1}\left(1+\alpha_1\right) &=&\sigma_1^k + {\alpha_1}{P\over\Delta} D_D, \\ - \sigma_2^{k+1}\left(1+\alpha_1\right) &=&\sigma_2^k + {\alpha_1\over e^2}{P\over\Delta} D_T, \\ - \sigma_{12}^{k+1}\left(1+\alpha_1\right) &=&\sigma_{12}^k + {\alpha_1\over 2e^2}{P\over\Delta} D_S.\end{aligned} - -To minimize code changes and unify the two approaches, we define and -apply :math:`1/\alpha_1` and :math:`\beta` in the classic EVP code, and -modify the elastic stress term. These under-relaxation parameters -control the rate at which the iteration converges. Thus for classic EVP -we set + a_{lvl} + a_{rdg} &=& a_i, \\ + v_{lvl} + v_{rdg} &=& v_i.\end{aligned} -.. math:: - \begin{aligned} - {\tt arlx1i} &=& {1\over\alpha_1} = {\Delta t_e\over 2T} \\ - {\tt brlx} &=& \beta = {\Delta t\over\Delta t_e}. \end{aligned} +Level ice area fraction and volume increase with new ice formation and +decrease steadily via ridging processes. Without the formation of new +ice, level ice asymptotes to zero because we assume that both level ice +and ridged ice ridge, in proportion to their fractional areas in a grid +cell (in the spirit of the ridging calculation itself which does not +prefer level ice over previously ridged ice). -Then +The ice strength :math:`P` may be computed in either of two ways. If the +namelist parameter kstrength = 0, we use the strength formula from +:cite:`Hibler79`: .. math:: - \begin{aligned} - {\tt denom1} &=& {1\over{1+{\tt arlx1i}}} = {1\over{1+1/\alpha_1}} = {1\over{1+\Delta t_e/ 2T}} \\ - {\tt c1} &=& {P\over\Delta}\,{\tt arlx1i} = {P\over\Delta}{\Delta t_e\over 2T} \\ - {\tt c0} &=& {{\tt c1}\over e^2} = {P\over\Delta}{\Delta t_e\over 2Te^2} .\end{aligned} + P = P^* h \exp[-C(1-a_i)], + :label: hib-strength -The stress equations for `stressp` (:math:`\sigma_1`) are unchanged; the -modified equations for `stressm` (:math:`\sigma_2`) and `stress12` -(:math:`\sigma_{12}`) take the form +where :math:`P^* = 27,500 \, \mathrm {N/m}` and :math:`C = 20` are +empirical constants, and :math:`h` is the mean ice thickness. +Alternatively, setting kstrength = 1 gives an ice strength closely +related to the ridging scheme. Following +:cite:`Rothrock75`, the strength is assumed proportional +to the change in ice potential energy :math:`\Delta E_P` per unit area +of compressive deformation. Given uniform ridge ITDs (krdg\_redist = 0), +we have .. math:: - \begin{aligned} - {\tt stressm} &=& {\tt stressm + c0}\,D_T \,{\tt denom1}\\ - {\tt stress12} &=& {\tt stress12 + 0.5\,c0}\,D_S \,{\tt denom1}.\end{aligned} - -For classic EVP, - -.. math:: - {\tt cca} = a = {\tt brlx}\,{m\over\Delta t} + {\tt vrel}\cos\theta ={m\over\Delta t_e} + {\tt vrel}\cos\theta. + P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} + \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} + \left( \frac{(H_n^{\max})^3 - (H_n^{\min})^3} + {3(H_n^{\max}-H_n^{\min})} \right) \right], + :label: roth-strength0 -For revised EVP, arlx1i and brlx are defined separately from -:math:`\Delta t`, :math:`\Delta t_e`, :math:`T` and :math:`e`, and +where :math:`C_P = (g/2)(\rho_i/\rho_w)(\rho_w-\rho_i)`, +:math:`\beta =R_{\mathrm{tot}}/R_{\mathrm{net}} > 1` +from Equation :eq:`Rtot-Rnet`, and :math:`C_f` is an empirical parameter that +accounts for frictional energy dissipation. Following +:cite:`FH95`, we set :math:`C_f = 17`. The first term in +the summation is the potential energy of ridging ice, and the second, +larger term is the potential energy of the resulting ridges. The factor +of :math:`\beta` is included because :math:`a_{Pn}` is normalized with +respect to the total area of ice ridging, not the net area removed. +Recall that more than one unit area of ice must be ridged to reduce the +net ice area by one unit. For exponential ridge ITDs (`krdg\_redist` = 1), +the ridge potential energy is modified: -.. math:: - {\tt cca} = \tilde{a} = \left(1+ {\tt brlx}\right){m\over\Delta t} + {\tt vrel}\cos\theta= \left(1+\beta\right){m\over\Delta t} + {\tt vrel}\cos\theta. +.. math:: + P = C_f \, C_p \, \beta \sum_{n=1}^{N_C} + \left[ -a_{Pn} \, h_n^2 + \frac{a_{Pn}}{k_n} + \left( H_{\min}^2 + 2H_{\min}\lambda + 2 \lambda^2 \right) \right] % CHECK BRACES + :label: roth-strength1 -:math:`\tilde{\bf u}` must also be defined for revised EVP as in -Equation :eq:`tildeu`. The extra terms in :math:`\tilde{a}` and -:math:`\tilde{\bf u}` are multiplied by a flag (revp) that equals 1 for -revised EVP and 0 for classic EVP. Revised EVP is activated by setting -the namelist parameter `revised\_evp` = true. Note that in the current -implementation, only the modified version of the elastic term is -available for either the classic (`revised\_evp` = false) or the revised -EVP method. A final difference is that the revised approach initializes -the stresses to 0 at the beginning of each time step, while the classic -EVP approach uses the previous time step value. +The energy-based ice strength given by Equations :eq:`roth-strength0` or +:eq:`roth-strength1` is more physically realistic than the strength +given by Equation :eq:`hib-strength`. However, use of Equation :eq:`hib-strength` is +less likely to allow numerical instability at a given resolution and +time step. See :cite:`LHMJ07` for more details. .. _thermo: +************** Thermodynamics --------------- +************** The current CICE version includes three thermodynamics options, the “zero-layer" thermodynamics of :cite:`Semtner76` @@ -2511,8 +3092,9 @@ surface temperatures. .. _ponds: +`````````` Melt ponds -~~~~~~~~~~ +`````````` Three explicit melt pond parameterizations are available in CICE, and all must use the delta-Eddington radiation scheme, described below. The @@ -3116,8 +3698,9 @@ the newly formed ice area :math:`\Delta a_i = \Delta a_{lvl}`. .. _sfc-forcing: +````````````````````````````````````` Thermodynamic surface forcing balance -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +````````````````````````````````````` The net surface energy flux from the atmosphere to the ice (with all fluxes defined as positive downward) is @@ -3292,8 +3875,9 @@ Section :ref:`thermo-growth`). .. _thermo-temp: +```````````````` New temperatures -~~~~~~~~~~~~~~~~ +```````````````` **Zero-layer thermodynamics** (`ktherm` = 0) An option for zero-layer thermodynamics :cite:`Semtner76` is @@ -4101,8 +4685,9 @@ may reduce ice formation by a small amount afterwards. .. _thermo-growth: +`````````````````` Growth and melting -~~~~~~~~~~~~~~~~~~ +`````````````````` Melting at the top surface is given by diff --git a/doc/source/cice_3_coupling.rst b/doc/source/cice_3_coupling.rst deleted file mode 100644 index 58a590404..000000000 --- a/doc/source/cice_3_coupling.rst +++ /dev/null @@ -1,549 +0,0 @@ -.. _coupl: - -Coupling with other climate model components -============================================ - -The sea ice model exchanges information with the other model components -via a flux coupler. CICE has been coupled into numerous climate models -with a variety of coupling techniques. This document is oriented -primarily toward the CESM Flux Coupler :cite:`KL02` -from NCAR, the first major climate model to incorporate CICE. The flux -coupler was originally intended to gather state variables from the -component models, compute fluxes at the model interfaces, and return -these fluxes to the component models for use in the next integration -period, maintaining conservation of momentum, heat, and fresh water. -However, several of these fluxes are now computed in the ice model -itself and provided to the flux coupler for distribution to the other -components, for two reasons. First, some of the fluxes depend strongly -on the state of the ice, and vice versa, implying that an implicit, -simultaneous determination of the ice state and the surface fluxes is -necessary for consistency and stability. Second, given the various ice -types in a single grid cell, it is more efficient for the ice model to -determine the net ice characteristics of the grid cell and provide the -resulting fluxes, rather than passing several values of the state -variables for each cell. These considerations are explained in more -detail below. - -The fluxes and state variables passed between the sea ice model and the -CESM flux coupler are listed in :ref:`tab-flux-cpl`. By convention, -directional fluxes are positive downward. In CESM, the sea ice model may -exchange coupling fluxes using a different grid than the computational -grid. This functionality is activated using the namelist variable -``gridcpl_file``. Another namelist variable ``highfreq``, allows the -high-frequency coupling procedure implemented in the Regional Arctic -System Model (RASM). In particular, the relative atmosphere-ice velocity -(:math:`\vec{U}_a-\vec{u}`) is used instead of the full atmospheric -velocity for computing turbulent fluxes in the atmospheric boundary -layer. - -:ref:`tab-flux-cpl`: *Data exchanged between the CESM flux coupler and the sea ice model* - -.. _tab-flux-cpl: - -.. table:: Table 1 - - =========================== ====================================== ======================================================================================= - Variable Description Interaction with flux coupler - =========================== ====================================== ======================================================================================= - :math:`z_o` Atmosphere level height From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`\vec{U}_a` Wind velocity From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`Q_a` Specific humidity From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`\rho_a` Air density From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`\Theta_a` Air potential temperature From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`T_a` Air temperature From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`F_{sw\downarrow}` Incoming shortwave radiation From *atmosphere model* via flux coupler **to** *sea ice model* - (4 bands) - - :math:`F_{L\downarrow}` Incoming longwave radiation From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`F_{rain}` Rainfall rate From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`F_{snow}` Snowfall rate From *atmosphere model* via flux coupler **to** *sea ice model* - - :math:`F_{frzmlt}` Freezing/melting potential From *ocean model* via flux coupler **to** *sea ice model* - - :math:`T_w` Sea surface temperature From *ocean model* via flux coupler **to** *sea ice model* - - :math:`S` Sea surface salinity From *ocean model* via flux coupler **to** *sea ice model* - - :math:`\nabla H_o` Sea surface slope From *ocean model* via flux coupler **to** *sea ice model* - - :math:`\vec{U}_w` Surface ocean currents From *ocean model* via flux coupler **to** *sea ice model* - - :math:`\vec{\tau}_a` Wind stress From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`F_s` Sensible heat flux From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`F_l` Latent heat flux From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`F_{L\uparrow}` Outgoing longwave radiation From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`F_{evap}` Evaporated water From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`\alpha` Surface albedo (4 bands) From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`T_{sfc}` Surface temperature From *sea ice model* via flux coupler **to** *atmosphere model* - - :math:`F_{sw\Downarrow}` Penetrating shortwave radiation From *sea ice model* via flux coupler **to** *ocean model* - - :math:`F_{water}` Fresh water flux From *sea ice model* via flux coupler **to** *ocean model* - - :math:`F_{hocn}` Net heat flux to ocean From *sea ice model* via flux coupler **to** *ocean model* - - :math:`F_{salt}` Salt flux From *sea ice model* via flux coupler **to** *ocean model* - - :math:`\vec{\tau}_w` Ice-ocean stress From *sea ice model* via flux coupler **to** *ocean model* - - :math:`F_{bio}` Biogeochemical fluxes From *sea ice model* via flux coupler **to** *ocean model* - - :math:`a_{i}` Ice fraction From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* - - :math:`T^{ref}_{a}` 2m reference temperature (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* - - :math:`Q^{ref}_{a}` 2m reference humidity (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* - - :math:`F_{swabs}` Absorbed shortwave (diagnostic) From *sea ice model* via flux coupler **to** both *ocean and atmosphere models* - =========================== ====================================== ======================================================================================= - -The ice fraction :math:`a_i` (aice) is the total fractional ice -coverage of a grid cell. That is, in each cell, - -.. math:: - \begin{array}{cl} - a_{i}=0 & \mbox{if there is no ice} \\ - a_{i}=1 & \mbox{if there is no open water} \\ - 0 0 - -where :math:`\cos Z` is the cosine of the solar zenith angle. - -.. _ocean: - -Ocean ------ - -New sea ice forms when the ocean temperature drops below its freezing -temperature. In the Bitz and Lipscomb thermodynamics, -:cite:`BL99` :math:`T_f=-\mu S`, where :math:`S` is the -seawater salinity and :math:`\mu=0.054 \ ^\circ`/ppt is the ratio of the -freezing temperature of brine to its salinity (linear liquidus -approximation). For the mushy thermodynamics, :math:`T_f` is given by a -piecewise linear liquidus relation. The ocean model calculates the new -ice formation; if the freezing/melting potential -:math:`F_{frzmlt}` is positive, its value represents a certain -amount of frazil ice that has formed in one or more layers of the ocean -and floated to the surface. (The ocean model assumes that the amount of -new ice implied by the freezing potential actually forms.) - -If :math:`F_{frzmlt}` is negative, it is used to heat already -existing ice from below. In particular, the sea surface temperature and -salinity are used to compute an oceanic heat flux :math:`F_w` -(:math:`\left|F_w\right| \leq \left|F_{frzmlt}\right|`) which -is applied at the bottom of the ice. The portion of the melting -potential actually used to melt ice is returned to the coupler in -:math:`F_{hocn}`. The ocean model adjusts its own heat budget -with this quantity, assuming that the rest of the flux remained in the -ocean. - -In addition to runoff from rain and melted snow, the fresh water flux -:math:`F_{water}` includes ice melt water from the top surface -and water frozen (a negative flux) or melted at the bottom surface of -the ice. This flux is computed as the net change of fresh water in the -ice and snow volume over the coupling time step, excluding frazil ice -formation and newly accumulated snow. Setting the namelist option -update\_ocn\_f to true causes frazil ice to be included in the fresh -water and salt fluxes. - -There is a flux of salt into the ocean under melting conditions, and a -(negative) flux when sea water is freezing. However, melting sea ice -ultimately freshens the top ocean layer, since the ocean is much more -saline than the ice. The ice model passes the net flux of salt -:math:`F_{salt}` to the flux coupler, based on the net change -in salt for ice in all categories. In the present configuration, -ice\_ref\_salinity is used for computing the salt flux, although the ice -salinity used in the thermodynamic calculation has differing values in -the ice layers. - -A fraction of the incoming shortwave :math:`F_{sw\Downarrow}` -penetrates the snow and ice layers and passes into the ocean, as -described in Section :ref:`sfc-forcing`. - -Many ice models compute the sea surface slope :math:`\nabla H_\circ` -from geostrophic ocean currents provided by an ocean model or other data -source. In our case, the sea surface height :math:`H_\circ` is a -prognostic variable in POP—the flux coupler can provide the surface -slope directly, rather than inferring it from the currents. (The option -of computing it from the currents is provided in subroutine -*evp\_prep*.) The sea ice model uses the surface layer currents -:math:`\vec{U}_w` to determine the stress between the ocean and the ice, -and subsequently the ice velocity :math:`\vec{u}`. This stress, relative -to the ice, - -.. math:: - \begin{aligned} - \vec{\tau}_w&=&c_w\rho_w\left|{\vec{U}_w-\vec{u}}\right|\left[\left(\vec{U}_w-\vec{u}\right)\cos\theta - +\hat{k}\times\left(\vec{U}_w-\vec{u}\right)\sin\theta\right] \end{aligned} - -is then passed to the flux coupler (relative to the ocean) for use by -the ocean model. Here, :math:`\theta` is the turning angle between -geostrophic and surface currents, :math:`c_w` is the ocean drag -coefficient, :math:`\rho_w` is the density of seawater, and -:math:`\hat{k}` is the vertical unit vector. The turning angle is -necessary if the top ocean model layers are not able to resolve the -Ekman spiral in the boundary layer. If the top layer is sufficiently -thin compared to the typical depth of the Ekman spiral, then -:math:`\theta=0` is a good approximation. Here we assume that the top -layer is thin enough. - -For CICE run in stand-alone mode (i.e., uncoupled), a thermodynamic slab -ocean mixed-layer parameterization is available in **ice\_ocean.F90**. -The turbulent fluxes are computed above the water surface using the same -parameterizations as for sea ice, but with parameters appropriate for -the ocean. The surface flux balance takes into account the turbulent -fluxes, oceanic heat fluxes from below the mixed layer, and shortwave -and longwave radiation, including that passing through the sea ice into -the ocean. If the resulting sea surface temperature falls below the -salinity-dependent freezing point, then new ice (frazil) forms. -Otherwise, heat is made available for melting the ice. - -.. _formdrag: - -Variable exchange coefficients ------------------------------- - -In the default CICE setup, atmospheric and oceanic neutral drag -coefficients (:math:`c_u` and :math:`c_w`) are assumed constant in time -and space. These constants are chosen to reflect friction associated -with an effective sea ice surface roughness at the ice–atmosphere and -ice–ocean interfaces. Sea ice (in both Arctic and Antarctic) contains -pressure ridges as well as floe and melt pond edges that act as discrete -obstructions to the flow of air or water past the ice, and are a source -of form drag. Following :cite:`TFSFFKLB14` and based on -recent theoretical developments :cite:`LGHA12,LLCL11`, the -neutral drag coefficients can now be estimated from properties of the -ice cover such as ice concentration, vertical extent and area of the -ridges, freeboard and floe draft, and size of floes and melt ponds. The -new parameterization allows the drag coefficients to be coupled to the -sea ice state and therefore to evolve spatially and temporally. This -parameterization is contained in the subroutine *neutral\_drag\_coeffs* -and is accessed by setting `formdrag` = true in the namelist. - -Following :cite:`TFSFFKLB14`, consider the general case of -fluid flow obstructed by N randomly oriented obstacles of height -:math:`H` and transverse length :math:`L_y`, distributed on a domain -surface area :math:`S_T`. Under the assumption of a logarithmic fluid -velocity profile, the general formulation of the form drag coefficient -can be expressed as - -.. math:: - C_d=\frac{N c S_c^2 \gamma L_y H}{2 S_T}\left[\frac{\ln(H/z_0)}{\ln(z_{ref}/z_0)}\right]^2, - :label: formdrag - -where :math:`z_0` is a roughness length parameter at the top or bottom -surface of the ice, :math:`\gamma` is a geometric factor, :math:`c` is -the resistance coefficient of a single obstacle, and :math:`S_c` is a -sheltering function that takes into account the shielding effect of the -obstacle, - -.. math:: - S_{c}=\left(1-\exp(-s_l D/H)\right)^{1/2}, - :label: shelter - -with :math:`D` the distance between two obstacles and :math:`s_l` an -attenuation parameter. - -As in the original drag formulation in CICE (sections :ref:`atmo` and -:ref:`ocean`), :math:`c_u` and :math:`c_w` along with the transfer -coefficients for sensible heat, :math:`c_{\theta}`, and latent heat, -:math:`c_{q}`, are initialized to a situation corresponding to neutral -atmosphere–ice and ocean–ice boundary layers. The corresponding neutral -exchange coefficients are then replaced by coefficients that explicitly -account for form drag, expressed in terms of various contributions as - -.. math:: - \tt{Cdn\_atm} = \tt{Cdn\_atm\_rdg} + \tt{Cdn\_atm\_floe} + \tt{Cdn\_atm\_skin} + \tt{Cdn\_atm\_pond} , - :label: Cda - -.. math:: - \tt{Cdn\_ocn} = \tt{Cdn\_ocn\_rdg} + \tt{Cdn\_ocn\_floe} + \tt{Cdn\_ocn\_skin}. - :label: Cdw - -The contributions to form drag from ridges (and keels underneath the -ice), floe edges and melt pond edges can be expressed using the general -formulation of equation :eq:`formdrag` (see :cite:`TFSFFKLB14` for -details). Individual terms in equation :eq:`Cdw` are fully described in -:cite:`TFSFFKLB14`. Following :cite:`Arya75` -the skin drag coefficient is parametrized as - -.. math:: - { \tt{Cdn\_(atm/ocn)\_skin}}=a_{i} \left(1-m_{(s/k)} \frac{H_{(s/k)}}{D_{(s/k)}}\right)c_{s(s/k)}, \mbox{ if $\displaystyle\frac{H_{(s/k)}}{D_{(s/k)}}\ge\frac{1}{m_{(s/k)}}$,} - :label: skindrag - -where :math:`m_s` (:math:`m_k`) is a sheltering parameter that depends -on the average sail (keel) height, :math:`H_s` (:math:`H_k`), but is -often assumed constant, :math:`D_s` (:math:`D_k`) is the average -distance between sails (keels), and :math:`c_{ss}` (:math:`c_{sk}`) is -the unobstructed atmospheric (oceanic) skin drag that would be attained -in the absence of sails (keels) and with complete ice coverage, -:math:`a_{ice}=1`. - -Calculation of equations :eq:`formdrag` – :eq:`skindrag` requires that small-scale geometrical -properties of the ice cover be related to average grid cell quantities -already computed in the sea ice model. These intermediate quantities are -briefly presented here and described in more detail in -:cite:`TFSFFKLB14`. The sail height is given by - -.. math:: - H_{s} = \displaystyle 2\frac{v_{rdg}}{a_{rdg}}\left(\frac{\alpha\tan \alpha_{k} R_d+\beta \tan \alpha_{s} R_h}{\phi_r\tan \alpha_{k} R_d+\phi_k \tan \alpha_{s} R_h^2}\right), - :label: Hs - -and the distance between sails\ - -.. math:: - D_{s} = \displaystyle 2 H_s\frac{a_{i}}{a_{rdg}} \left(\frac{\alpha}{\tan \alpha_s}+\frac{\beta}{\tan \alpha_k}\frac{R_h}{R_d}\right), - :label: Ds - -where :math:`0<\alpha<1` and :math:`0<\beta<1` are weight functions, -:math:`\alpha_{s}` and :math:`\alpha_{k}` are the sail and keel slope, -:math:`\phi_s` and :math:`\phi_k` are constant porosities for the sails -and keels, and we assume constant ratios for the average keel depth and -sail height (:math:`H_k/H_s=R_h`) and for the average distances between -keels and between sails (:math:`D_k/D_s=R_d`). With the assumption of -hydrostatic equilibrium, the effective ice plus snow freeboard is -:math:`H_{f}=\bar{h_i}(1-\rho_i/\rho_w)+\bar{h_s}(1-\rho_s/\rho_w)`, -where :math:`\rho_i`, :math:`\rho_w` and :math:`\rho_s` are -respectively the densities of sea ice, water and snow, :math:`\bar{h_i}` -is the mean ice thickness and :math:`\bar{h_s}` is the mean snow -thickness (means taken over the ice covered regions). For the melt pond -edge elevation we assume that the melt pond surface is at the same level -as the ocean surface surrounding the floes -:cite:`FF07,FFT10,FSFH12` and use the simplification -:math:`H_p = H_f`. Finally to estimate the typical floe size -:math:`L_A`, distance between floes, :math:`D_F`, and melt pond size, -:math:`L_P` we use the parameterizations of :cite:`LGHA12` -to relate these quantities to the ice and pond concentrations. All of -these intermediate quantities are available as history output, along -with `Cdn\_atm`, `Cdn\_ocn` and the ratio `Cdn\_atm\_ratio\_n` between the -total atmospheric drag and the atmospheric neutral drag coefficient. - -We assume that the total neutral drag coefficients are thickness -category independent, but through their dependance on the diagnostic -variables described above, they vary both spatially and temporally. The -total drag coefficients and heat transfer coefficients will also depend -on the type of stratification of the atmosphere and the ocean, and we -use the parameterization described in section :ref:`atmo` that accounts -for both stable and unstable atmosphere–ice boundary layers. In contrast -to the neutral drag coefficients the stability effect of the atmospheric -boundary layer is calculated separately for each ice thickness category. - -The transfer coefficient for oceanic heat flux to the bottom of the ice -may be varied based on form drag considerations by setting the namelist -variable `fbot\_xfer\_type` to `Cdn\_ocn`; this is recommended when using -the form drag parameterization. Its default value of the transfer -coefficient is 0.006 (`fbot\_xfer\_type = ’constant’`). diff --git a/doc/source/cice_5_numerical_implementation.rst b/doc/source/cice_3_user_guide.rst similarity index 61% rename from doc/source/cice_5_numerical_implementation.rst rename to doc/source/cice_3_user_guide.rst index 9da32c54f..1c074a515 100644 --- a/doc/source/cice_5_numerical_implementation.rst +++ b/doc/source/cice_3_user_guide.rst @@ -1,6 +1,9 @@ +User Guide +========== +------------------------ Numerical implementation -======================== +------------------------ CICE is written in FORTRAN90 and runs on platforms using UNIX, LINUX, and other operating systems. The code is parallelized via grid @@ -20,8 +23,9 @@ another system without using MPI. .. _dirstructure: +~~~~~~~~~~~~~~~~~~~ Directory structure -------------------- +~~~~~~~~~~~~~~~~~~~ The present code distribution includes make files, several scripts and some input files. The main directory is **cice/**, and a run directory @@ -342,8 +346,9 @@ execution or “run” directory created when the code is compiled using the **run\_ice** batch run script file from **cice/input\_templates/** +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Grid, boundary conditions and masks ------------------------------------ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The spatial discretization is specialized for a generalized orthogonal B-grid as in :cite:`Murray96` or @@ -382,9 +387,9 @@ In CESM, the sea ice model may exchange coupling fluxes using a different grid than the computational grid. This functionality is activated using the namelist variable `gridcpl\_file`. - +*********************** Grid domains and blocks -~~~~~~~~~~~~~~~~~~~~~~~ +*********************** In general, the global gridded domain is `nx\_global` :math:`\times`\ `ny\_global`, while the subdomains used in the @@ -449,8 +454,9 @@ Alternatively, a new variable is provided in the history files, `blkmask`, which labels the blocks in the grid decomposition according to `blkmask` = `my\_task` + `iblk/100`. +************* Tripole grids -~~~~~~~~~~~~~ +************* The tripole grid is a device for constructing a global grid with a normal south pole and southern boundary condition, which avoids placing @@ -511,8 +517,9 @@ excluded. .. _bio-grid: +******** Bio-grid -~~~~~~~~ +******** The bio-grid is a vertical grid used for solving the brine height variable :math:`h_b`. In the future, it will also be used for @@ -533,8 +540,9 @@ spaced at :math:`1/n_b` intervals beginning with `bgrid(2)` :math:` = equidistant with the same spacing, but physically coincide with points midway between those of `bgrid`. +******************** Column configuration -~~~~~~~~~~~~~~~~~~~~ +******************** A column modeling capability is available. Because of the boundary conditions and other spatial assumptions in the model, this is not a @@ -552,8 +560,9 @@ History variables available for column output are ice and snow temperature, `Tinz` and `Tsnz`. These variables also include thickness category as a fourth dimension. +******************* Boundary conditions -~~~~~~~~~~~~~~~~~~~ +******************* Much of the infrastructure used in CICE, including the boundary routines, is adopted from POP. The boundary routines perform boundary @@ -593,8 +602,9 @@ the tripole grid, both within the dynamics calculation and for restarts. This has not been implemented yet for tripoleT grids, pending further testing. +***** Masks -~~~~~ +***** A land mask hm (:math:`M_h`) is specified in the cell centers, with 0 representing land and 1 representing ocean cells. A corresponding mask @@ -631,10 +641,15 @@ or southern hemispheres, respectively. Special constants (`spval` and `spval\_dbl`, each equal to :math:`10^{30}`) are used to indicate land points in the history files and diagnostics. +~~~~~~~~~~~~~~~~~~~ +Test configurations +~~~~~~~~~~~~~~~~~~~ + .. _init: +~~~~~~~~~~~~~~~~~~~~~~~~~~~ Initialization and coupling ---------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~ The ice model’s parameters and variables are initialized in several steps. Many constants and physical parameters are set in @@ -707,8 +722,9 @@ reset to ‘none.’* .. _parameters: +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Choosing an appropriate time step ---------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The time step is chosen based on stability of the transport component (both horizontal and in thickness space) and on resolution of the @@ -802,13 +818,15 @@ temperature :math:`T_{sfc}` is computed internally. The numerical constraint on the thermodynamics time step is associated with the transport scheme rather than the thermodynamic solver. +~~~~~~~~~~~~ Model output ------------- +~~~~~~~~~~~~ .. _history: +************* History files -~~~~~~~~~~~~~ +************* Model output data is averaged over the period(s) given by `histfreq` and `histfreq\_n`, and written to binary or  files prepended by `history\_file` @@ -901,8 +919,9 @@ another that is multiplied by :math:`a_i`, representing an average over the grid cell area. Our naming convention attaches the suffix “\_ai" to the grid-cell-mean variable names. +**************** Diagnostic files -~~~~~~~~~~~~~~~~ +**************** Like `histfreq`, the parameter `diagfreq` can be used to regulate how often output is written to a log file. The log file unit to which diagnostic @@ -980,8 +999,9 @@ The timers use *MPI\_WTIME* for parallel runs and the F90 intrinsic | 16 | BGC | biogeochemistry | +--------------+-------------+----------------------------------------------------+ +************* Restart files -~~~~~~~~~~~~~ +************* CICE now provides restart data in binary unformatted or  formats, via the `IO\_TYPE` flag in **comp\_ice** and namelist variable @@ -1040,6 +1060,7 @@ and format, provided that the same number of ice layers and basic physics packages will be used for the new runs. See Section :ref:`restarttrouble` for details. +-------------------- Execution procedures -------------------- @@ -1188,8 +1209,26 @@ filename can be assigned to ice\_ic in **ice\_in**. Consult Section :ref:`init` for more details. Restarts are exact for MPI or single processor runs. +~~~~~~~ +Scripts +~~~~~~~ + +~~~~~~~~~~~ +Directories +~~~~~~~~~~~ + +~~~~~~~~~~~~~~~~~~~ +Local modifications +~~~~~~~~~~~~~~~~~~~ + +~~~~~~~~~~~~ +Forcing data +~~~~~~~~~~~~ + + .. _performance: +----------- Performance ----------- @@ -1340,11 +1379,13 @@ Delta-Eddington radiation scheme is required for all melt pond schemes and the aerosol tracers, and the level-ice pond parameterization additionally requires the level-ice tracers. +------------- Adding things ------------- .. _addtimer: +~~~~~~ Timers ~~~~~~ @@ -1364,6 +1405,7 @@ desired, by including the block ID in the timer calls. .. _addhist: +~~~~~~~~~~~~~~ History fields ~~~~~~~~~~~~~~ @@ -1394,6 +1436,7 @@ section :ref:`history`. .. _addtrcr: +~~~~~~~ Tracers ~~~~~~~ @@ -1485,3 +1528,943 @@ a pattern. #. If strict conservation is necessary, add diagnostics as noted for topo ponds in Section :ref:`ponds`. +--------------- +Troubleshooting +--------------- + +Check the FAQ: http://oceans11.lanl.gov/drupal/CICE/FAQ. + +.. _setup: + +~~~~~~~~~~~~~ +Initial setup +~~~~~~~~~~~~~ + +The script **comp\_ice** is configured so that the files **grid**, +**kmt**, **ice\_in**, **run\_ice**, **iced\_gx3\_v5.0** and +**ice.restart\_file** are NOT overwritten after the first setup. If you +wish to make changes to the original files in **input\_templates/** +rather than those in the run directory, either remove the files from the +run directory before executing **comp\_ice** or edit the script. + +The code may abort during the setup phase for any number of reasons, and +often the buffer containing the diagnostic output fails to print before +the executable exits. The quickest way to get the diagnostic information +is to run the code in an interactive shell with just the command `cice` +for serial runs or “`mpirun -np N cice`” for MPI runs, where N is the +appropriate number of processors (or a command appropriate for your +computer’s software). + +If the code fails to compile or run, or if the model configuration is +changed, try the following: + +- create **Macros.\***. **Makefile.\*** and **run\_ice.\*** files for + your particular platform, if they do not already exist (type ‘uname + -s’ at the prompt and compare the result with the file suffixes; we + rename `UNICOS/mp` as `UNICOS` for simplicity). + +- modify the `INCLUDE` directory path and other settings for your system + in the scripts, **Macros.\*** and **Makefile.\*** files. + +- alter directory paths, file names and the execution command as needed + in **run\_ice** and **ice\_in**. + +- ensure that `nprocs` in **ice\_in** is equal to `NTASK` in **comp\_ice**. + +- ensure that the block size `NXBLOCK`, `NYBLOCK` in **comp\_ice** is + compatible with the processor\_shape and other domain options in + **ice\_in** + +- if using the rake or space-filling curve algorithms for block + distribution (`distribution\_type` in **ice\_in**) the code will abort + if `MXBLCKS` is not large enough. The correct value is provided in the + diagnostic output. + +- if starting from a restart file, ensure that kcatbound is the same as + that used to create the file (`kcatbound` = 0 for the files included in + this code distribution). Other configuration parameters, such as + `NICELYR`, must also be consistent between runs. + +- for stand-alone runs, check that `-Dcoupled` is *not* set in the + **Macros.\*** file. + +- for coupled runs, check that `-Dcoupled` and other + coupled-model-specific (e.g., CESM, popcice or hadgem) preprocessing + options are set in the **Macros.\*** file. + +- edit the grid size and other parameters in **comp\_ice**. + +- remove the **compile/** directory completely and recompile. + +.. _restarttrouble: + +~~~~~~~~ +Restarts +~~~~~~~~ + +CICE version 5 introduces a new model configuration that makes +restarting from older simulations difficult. In particular, the number +of ice categories, the category boundaries, and the number of vertical +layers within each category must be the same in the restart file and in +the run restarting from that file. Moreover, significant differences in +the physics, such as the salinity profile, may cause the code to fail +upon restart. Therefore, new model configurations may need to be started +using `runtype` = ‘initial’. Binary restart files that were provided with +CICE v4.1 were made using the BL99 thermodynamics with 4 layers and 5 +thickness categories (`kcatbound` = 0) and therefore can not be used for +the default CICE v5 configuration (7 layers). In addition, CICE’s +default restart file format is now  instead of binary. + +Restarting a run using `runtype` = ‘continue’ requires restart data for +all tracers used in the new run. If tracer restart data is not +available, use `runtype` = ‘initial’, setting `ice\_ic` to the name of the +core restart file and setting to true the namelist restart flags for +each tracer that is available. The unavailable tracers will be +initialized to their default settings. + +On tripole grids, use `restart\_ext` = true when using either binary or +regular (non-PIO) netcdf. + +Provided that the same number of ice layers (default: 4) will be used +for the new runs, it is possible to convert v4.1 restart files to the +new file structure and then to  format. If the same physical +parameterizations are used, the code should be able to execute from +these files. However if different physics is used (for instance, mushy +thermo instead of BL99), the code may still fail. To convert a v4.1 +restart file: + +#. Edit the code **input\_templates/convert\_restarts.f90** for your + model configuration and path names. Compile and run this code to + create a binary restart file that can be read using v5. Copy the + resulting file to the **restart/** subdirectory in your working + directory. + +#. In your working directory, turn off all tracer restart flags in + **ice\_in** and set the following: + + - runtype = ‘initial’ + + - ice\_ic = ‘./restart/[your binary file name]’ + + - restart = .true. + + - use\_restart\_time = .true. + +#. In **CICE\_InitMod.F90**, comment out the call to + restartfile(ice\_ic) and uncomment the call to + restartfile\_v4(ice\_ic) immediately below it. This will read the + v4.1 binary file and write a v5  file containing the same + information. + +If restart files are taking a long time to be written serially (i.e., +not using PIO), see the next section. + +~~~~~~~~~~~~~~ +Slow execution +~~~~~~~~~~~~~~ + +On some architectures, underflows (:math:`10^{-300}` for example) are +not flushed to zero automatically. Usually a compiler flag is available +to do this, but if not, try uncommenting the block of code at the end of +subroutine *stress* in **ice\_dyn\_evp.F90** or **ice\_dyn\_eap.F90**. +You will take a hit for the extra computations, but it will not be as +bad as running with the underflows. + +In some configurations, multiple calls to scatter or gather global +variables may overfill MPI’s buffers, causing the code to slow down +(particularly when writing large output files such as restarts). To +remedy this problem, set `BARRIERS yes` in **comp\_ice**. This +synchronizes MPI messages, keeping the buffers in check. + +~~~~~~~~~~~~~~~ +Debugging hints +~~~~~~~~~~~~~~~ + +Several utilities are available that can be helpful when debugging the +code. Not all of these will work everywhere in the code, due to possible +conflicts in module dependencies. + +*debug\_ice* (**CICE.F90**) + A wrapper for *print\_state* that is easily called from numerous + points during the timestepping loop (see + **CICE\_RunMod.F90\_debug**, which can be substituted for + **CICE\_RunMod.F90**). + +*print\_state* (**ice\_diagnostics.F90**) + Print the ice state and forcing fields for a given grid cell. + +`dbug` = true (**ice\_in**) + Print numerous diagnostic quantities. + +`print\_global` (**ice\_in**) + If true, compute and print numerous global sums for energy and mass + balance analysis. This option can significantly degrade code + efficiency. + +`print\_points` (**ice\_in**) + If true, print numerous diagnostic quantities for two grid cells, + one near the north pole and one in the Weddell Sea. This utility + also provides the local grid indices and block and processor numbers + (`ip`, `jp`, `iblkp`, `mtask`) for these points, which can be used in + conjunction with `check\_step`, to call *print\_state*. These flags + are set in **ice\_diagnostics.F90**. This option can be fairly slow, + due to gathering data from processors. + +*global\_minval, global\_maxval, global\_sum* (**ice\_global\_reductions.F90**) + Compute and print the minimum and maximum values for an individual + real array, or its global sum. + +~~~~~~~~~~ +Known bugs +~~~~~~~~~~ + +#. Fluxes sent to the CESM coupler may have incorrect values in grid + cells that change from an ice-free state to having ice during the + given time step, or vice versa, due to scaling by the ice area. The + authors of the CESM flux coupler insist on the area scaling so that + the ice and land models are treated consistently in the coupler (but + note that the land area does not suddenly become zero in a grid cell, + as does the ice area). + +#. With the old CCSM radiative scheme (`shortwave` = ‘default’ or + ‘ccsm3’), a sizable fraction (more than 10%) of the total shortwave + radiation is absorbed at the surface but should be penetrating into + the ice interior instead. This is due to use of the aggregated, + effective albedo rather than the bare ice albedo when + `snowpatch` :math:`< 1`. + +#. The date-of-onset diagnostic variables, `melt\_onset` and `frz\_onset`, + are not included in the core restart file, and therefore may be + incorrect for the current year if the run is restarted after Jan 1. + Also, these variables were implemented with the Arctic in mind and + may be incorrect for the Antarctic. + +#. The single-processor *system\_clock* time may give erratic results on + some architectures. + +#. History files that contain time averaged data (`hist\_avg` = true in + **ice\_in**) will be incorrect if restarting from midway through an + averaging period. + +#. In stand-alone runs, restarts from the end of `ycycle` will not be + exact. + +#. Using the same frequency twice in `histfreq` will have unexpected + consequences and causes the code to abort. + +#. Latitude and longitude fields in the history output may be wrong when + using padding. + +~~~~~~~~~~~~~~~~~~~~~~~~~ +Interpretation of albedos +~~~~~~~~~~~~~~~~~~~~~~~~~ + +The snow-and-ice albedo, `albsni`, and diagnostic albedos `albice`, `albsno`, +and `albpnd` are merged over categories but not scaled (divided) by the +total ice area. (This is a change from CICE v4.1 for `albsni`.) The latter +three history variables represent completely bare or completely snow- or +melt-pond-covered ice; that is, they do not take into account the snow +or melt pond fraction (`albsni` does, as does the code itself during +thermodyamic computations). This is to facilitate comparison with +typical values in measurements or other albedo parameterizations. The +melt pond albedo `albpnd` is only computed for the Delta-Eddington +shortwave case. + +With the Delta-Eddington parameterization, the albedo depends on the +cosine of the zenith angle (:math:`\cos\varphi`, `coszen`) and is zero if +the sun is below the horizon (:math:`\cos\varphi < 0`). Therefore +time-averaged albedo fields would be low if a diurnal solar cycle is +used, because zero values would be included in the average for half of +each 24-hour period. To rectify this, a separate counter is used for the +averaging that is incremented only when :math:`\cos\varphi > 0`. The +albedos will still be zero in the dark, polar winter hemisphere. + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Proliferating subprocess parameterizations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +With the addition of several alternative parameterizations for sea ice +processes, a number of subprocesses now appear in multiple parts of the +code with differing descriptions. For instance, sea ice porosity and +permeability, along with associated flushing and flooding, are +calculated separately for mushy thermodynamics, topo and level-ice melt +ponds, and for the brine height tracer, each employing its own +equations. Likewise, the BL99 and mushy thermodynamics compute freeboard +and snow–ice formation differently, and the topo and level-ice melt pond +schemes both allow fresh ice to grow atop melt ponds, using slightly +different formulations for Stefan freezing. These various process +parameterizations will be compared and their subprocess descriptions +possibly unified in the future. + +------------ +Testing CICE +------------ + +Version 6, August 2017 +This documents how to use the testing features developed for the +CICE Consortium CICE sea ice model. + +.. _basic: + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Individual tests and test suites +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The CICE scripts support both setup of individual tests as well as test suites. Individual +tests are run from the command line like + + > create.case -t smoke -m wolf -g gx3 -p 8x2 -s diag1,run5day -testid myid + +where -m designates a specific machine. Test suites are multiple tests that are specified in +an input file and are started on the command line like + + > create.case -ts base_suite -m wolf -testid myid + +create.case with -t or -ts require a testid to uniquely name test directories. The format +of the case directory name for a test will always be +${machine}_${test}_${grid}_${pes}_${soptions}.${testid} + +To build and run a test, the process is the same as a case, + cd into the test directory, + + run cice.build + + run cice.submit + +The test results will be generated in a local file called "test_output". + +When running a test suite, the create.case command line automatically generates all the tests +under a directory names ${test_suite}.${testid}. It then automatically builds and submits all +tests. When the tests are complete, run the results.csh script to see the results from all the +tests. + +Tests are defined under configuration/scripts/tests. The tests currently supported are: + smoke - Runs the model for default length. The length and options can + be set with the -s commmand line option. The test passes if the + model completes successfully. + restart - Runs the model for 10 days, writing a restart file at day 5 and + again at day 10. Runs the model a second time starting from the + day 5 restart and writing a restart at day 10 of the model run. + The test passes if both the 10 day and 5 day restart run complete and + if the restart files at day 10 from both runs are bit-for-bit identical. + +Please run './create.case -h' for additional details. + +.. _additional: + +~~~~~~~~~~~~~~~~~~~~~~~~~~ +Additional testing options +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +There are several additional options on the create.case command line for testing that +provide the ability to regression test and compare tests to each other. + + -bd defines a baseline directory where tests can be stored for regression testing + + -bg defines a version name that where the current tests can be saved for regression testing + + -bc defines a version name that the current tests should be compared to for regression testing + + -td provides a way to compare tests with each other + +To use -bg, + > create.case -ts base_suite -m wolf -testid v1 -bg version1 -bd $SCRATCH/CICE_BASELINES + will copy all the results from the test suite to $SCRATCH/CICE_BASELINES/version1. + +To use -bc, + > create.case -ts base_suite -m wolf -testid v2 -bc version1 -bd $SCRATCH/CICE_BASELINES + will compare all the results from this test suite to results saved before in $SCRATCH/CICE_BASELINES/version1. + +-bc and -bg can be combined, + >create.case -ts base_suite -m wolf -testid v2 -bg version2 -bc version1 -bd $SCRATCH/CICE_BASELINES + will save the current results to $SCRATCH/CICE_BASELINES/version2 and compare the current results to + results save before in $SCRATCH/CICE_BASELINES/version1. + +-bg, -bc, and -bd are used for regression testing. There is a default -bd on each machine. + +-td allows a user to compare one test result to another. For instance, + > create.case -t smoke -m wolf -g gx3 -p 8x2 -s run5day -testid t01 + > create.case -t smoke -m wolf -g gx3 -p 4x2 -s run5day -testid t01 -td smoke_gx3_8x2_run5day + +An additional check will be done for the second test (because of the -td argument), and it will compare +the output from the first test "smoke_gx3_8x2_run5day" to the output from it's test "smoke_gx3_4x2_run5day" +and generate a result for that. It's important that the first test complete before the second test is +done. Also, the -td option works only if the testid and the machine are the same for the baseline +run and the current run. + +.. _format: + +~~~~~~~~~~~~~~~~~ +Test suite format +~~~~~~~~~~~~~~~~~ + +The format for the test suite file is relatively simple. It is a text file with white space delimited +columns like, + +.. _tab-test: + +.. csv-table:: Table 7 + :header: "#Test", "Grid", "PEs", "Sets", "BFB-compare" + :widths: 7, 7, 7, 15, 15 + + "smoke", "gx3", "8x2", "diag1,run5day", "" + "smoke", "gx3", "8x2", "diag24,run1year,medium", "" + "smoke", "gx3", "4x1", "debug,diag1,run5day", "" + "smoke", "gx3", "8x2", "debug,diag1,run5day", "" + "smoke", "gx3", "4x2", "diag1,run5day", "smoke_gx3_8x2_diag1_run5day" + "smoke", "gx3", "4x1", "diag1,run5day,thread", "smoke_gx3_8x2_diag1_run5day" + "smoke", "gx3", "4x1", "diag1,run5day", "smoke_gx3_4x1_diag1_run5day_thread" + "restart", "gx3", "8x1", "", "" + "restart", "gx3", "4x2", "debug", "" + + +The first column is the test name, the second the grid, the third the pe count, the fourth column is +the -s options and the fifth column is the -td argument. The fourth and fifth columns are optional. +The argument to -ts defines which filename to choose and that argument can contain a path. create.case +will also look for the filename in configuration/scripts/tests where some preset test suites are defined. + +~~~~~~~~~~~~~~~~~~~~~~~~~~ +Example Tests (Quickstart) +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +********************************************** +To generate a baseline dataset for a test case +********************************************** + +./create.case -t smoke -m wolf -bg cicev6.0.0 -testid t00 + +cd wolf_smoke_gx3_4x1.t00 + +./cice.build + +./cice.submit + +# After job finishes, check output + +cat test_output + +**************************************************** +To run a test case and compare to a baseline dataset +**************************************************** + +./create.case -t smoke -m wolf -bc cicev6.0.0 -testid t01 + +cd wolf_smoke_gx3_4x1.t01 + +./cice.build + +./cice.submit + +# After job finishes, check output + +cat test_output + +********************************************* +To run a test suite to generate baseline data +********************************************* + +./create.case -m wolf -ts base_suite -testid t02 -bg cicev6.0.0bs + +cd base_suite.t02 + +# Once all jobs finish, concatenate all output + +./results.csh # All tests results will be stored in results.log + +# To plot a timeseries of "total ice extent", "total ice area", and "total ice volume" + +./timeseries.csh + +ls \*.png + +*********************************************** +To run a test suite to compare to baseline data +*********************************************** + +./create.case -m wolf -ts base_suite -testid t03 -bc cicev6.0.0bs + +cd base_suite.t03 + +# Once all jobs finish, concatenate all output + +./results.csh # All tests results will be stored in results.log + +# To plot a timeseries of "total ice extent", "total ice area", and "total ice volume" + +./timeseries.csh + +ls \*.png + +************************** +To compare to another test +************************** +`First:` + +./create.case -m wolf -t smoke -testid t01 -p 8x2 + +cd wolf_smoke_gx3_8x2.t01 + +./cice.build + +./cice.submit + +# After job finishes, check output + +cat test_output + +`Then, do the comparison:` + +./create.case -m wolf -t smoke -testid t01 -td smoke_gx3_8x2 -s thread -p 4x1 + +cd wolf_smoke_gx3_4x1_thread.t01 + +./cice.build + +./cice.submit + +# After job finishes, check output + +cat test_output + +****************** +Additional Details +****************** + +- In general, the baseline generation, baseline compare, and test diff are independent. +- Use the '-bd' flag to specify the location where you want the baseline dataset + to be written. Without specifying '-bd', the baseline dataset will be written + to the default baseline directory found in the env. file (ICE_MACHINE_BASELINE). +- If '-bd' is not passed, the scripts will look for baseline datasets in the default + baseline directory found in the env. file (ICE_MACHINE_BASELINE). + If the '-bd' option is passed, the scripts will look for baseline datasets in the + location passed to the -bd argument. +- To generate a baseline dataset for a specific version (for regression testing), + use '-bg '. The scripts will then place the baseline dataset + in $ICE_MACHINE_BASELINE// +- The '-testid' flag allows users to specify a testing id that will be added to the + end of the case directory. For example, "./create.case -m wolf -t smoke -testid t12 -p 4x1" + creates the directory wolf_smoke_gx3_4x1.t12. This flag is REQUIRED if using -t or -ts. + +~~~~~~~~~~~~~~~~~~~~ +Code compliance test +~~~~~~~~~~~~~~~~~~~~ + +Additions and changes to CICE and Icepack are expected to be bit-for-bit +unless there is a strong justification for non-reproducibility, such as +a bug-fix or approved scientific alteration to existing code. However, +situations do arise when additions to CICE or Icepack are not +bit-for-bit, but are also not expected to change the science of CICE and +Icepack. In that instant, further evidence is required in the initial +testing phase to support the premise that code changes have not altered +the science of the model. To support this testing, a :math:`t`-test is +being/has been implemented in the CICE testing infrastructure. + +****** +Method +****** + +Welch’s two-sided :math:`t`-test is used to help determine whether or +not two different simulations that should be identical are significantly +different for any grid cell of the CICE gx-1/3 domain for grid-cell +averaged sea ice thickness, :math:`h`, ice concentration, :math:`c`, and +pack velocity components :math:`\pmb{u}=\pmb{u}(u,v)`. In this +circumstance, we seek to determine whether or not the null hypothesis, +:math:`H_0`, is true. The null hypothesis is: Two simulations that are +not bit-for-bit identical are ostensibly the same at every model grid +point. The test begins from the standpoint that two CICE simulations +*should* be bit-for-bit but are suspected of only being different at the +level of computational innaccuracy. Therefore, we seek to limit a +:math:`t`-test Type II error, where a test would erroneously confirm the +null hypothesis, :math:`H_0`. To that end, we choose to test the +hypothesis that grid-point means from CICE simulation ‘:math:`a`’ are +different from CICE simulation ‘:math:`b`’ at a relatively low +confidence interval. Formally, we test the hypothesis +:math:`H_0:\bar{x}_a=\bar{x}_b`, :math:`H_1:\bar{x}_a\neq\bar{x}_b` for +each of the aforementioned variables at every model grid point using a +two-sided t-test with a 68, 80 and 95% confidence interval. Here, +:math:`\bar{x}{=}\tfrac{1}{n}\sum_{i=1}^n x_i` is the time series mean +of :math:`n` samples :math:`x_i` representing :math:`h`, :math:`c`, +:math:`u` or :math:`v`, and daily samples are used from 5-year +stand-alone CICE simulations. More frequent output is unnecessary, +because each of :math:`h`, :math:`c`, :math:`u` and :math:`v` typically +have a high degree of auto-correlation in sea ice models. + +Due to the strong auto-correlation in geo-located sea ice time series, +we calculate a two-sided :math:`t`-statistic to compare +:math:`\bar{x}_a` and :math:`\bar{x}_b`, given their respective standard +deviations, :math:`\sigma_a` and :math:`\sigma_b`, and effective sample +sizes, :math:`n'_a` and :math:`n'_b`, following +:cite:`vSZ99` : + +.. math:: + t=\frac{\bar{x}_a - \bar{x}_b}{\sqrt{\frac{\sigma^2_a}{n'_a}+\frac{\sigma^2_b}{n'_b}}}. + :label: t-distribution + +The null hypothesis :math:`H_0:\bar{x}_a=\bar{x}_b` is true when + +.. math:: + -t_{crit}({1{-}\alpha/2},N)``", "string", "frequency units for writing ```` to history", "" + "", "``y``", "write history every ``histfreq_n`` years", "" + "", "``m``", "write history every ``histfreq_n`` months", "" + "", "``d``", "write history every ``histfreq_n`` days", "" + "", "``h``", "write history every ``histfreq_n`` hours", "" + "", "``1``", "write history every time step", "" + "", "``x``", "do not write ```` to history", "" + "", "``md``", "*e.g.,* write both monthly and daily files", "" + "``f__ai``", "", "grid cell average of ```` (:math:`\times a_i`)", "" + diff --git a/doc/source/cice_10_index.rst b/doc/source/cice_4_index.rst similarity index 100% rename from doc/source/cice_10_index.rst rename to doc/source/cice_4_index.rst diff --git a/doc/source/cice_6_troubleshooting.rst b/doc/source/cice_6_troubleshooting.rst deleted file mode 100644 index 7e5dfc787..000000000 --- a/doc/source/cice_6_troubleshooting.rst +++ /dev/null @@ -1,263 +0,0 @@ -.. role:: math(raw) - :format: html latex -.. - -Troubleshooting -================ - -Check the FAQ: http://oceans11.lanl.gov/drupal/CICE/FAQ. - -.. _setup: - -Initial setup -------------- - -The script **comp\_ice** is configured so that the files **grid**, -**kmt**, **ice\_in**, **run\_ice**, **iced\_gx3\_v5.0** and -**ice.restart\_file** are NOT overwritten after the first setup. If you -wish to make changes to the original files in **input\_templates/** -rather than those in the run directory, either remove the files from the -run directory before executing **comp\_ice** or edit the script. - -The code may abort during the setup phase for any number of reasons, and -often the buffer containing the diagnostic output fails to print before -the executable exits. The quickest way to get the diagnostic information -is to run the code in an interactive shell with just the command `cice` -for serial runs or “`mpirun -np N cice`” for MPI runs, where N is the -appropriate number of processors (or a command appropriate for your -computer’s software). - -If the code fails to compile or run, or if the model configuration is -changed, try the following: - -- create **Macros.\***. **Makefile.\*** and **run\_ice.\*** files for - your particular platform, if they do not already exist (type ‘uname - -s’ at the prompt and compare the result with the file suffixes; we - rename `UNICOS/mp` as `UNICOS` for simplicity). - -- modify the `INCLUDE` directory path and other settings for your system - in the scripts, **Macros.\*** and **Makefile.\*** files. - -- alter directory paths, file names and the execution command as needed - in **run\_ice** and **ice\_in**. - -- ensure that `nprocs` in **ice\_in** is equal to `NTASK` in **comp\_ice**. - -- ensure that the block size `NXBLOCK`, `NYBLOCK` in **comp\_ice** is - compatible with the processor\_shape and other domain options in - **ice\_in** - -- if using the rake or space-filling curve algorithms for block - distribution (`distribution\_type` in **ice\_in**) the code will abort - if `MXBLCKS` is not large enough. The correct value is provided in the - diagnostic output. - -- if starting from a restart file, ensure that kcatbound is the same as - that used to create the file (`kcatbound` = 0 for the files included in - this code distribution). Other configuration parameters, such as - `NICELYR`, must also be consistent between runs. - -- for stand-alone runs, check that `-Dcoupled` is *not* set in the - **Macros.\*** file. - -- for coupled runs, check that `-Dcoupled` and other - coupled-model-specific (e.g., CESM, popcice or hadgem) preprocessing - options are set in the **Macros.\*** file. - -- edit the grid size and other parameters in **comp\_ice**. - -- remove the **compile/** directory completely and recompile. - -.. _restarttrouble: - -Restarts --------- - -CICE version 5 introduces a new model configuration that makes -restarting from older simulations difficult. In particular, the number -of ice categories, the category boundaries, and the number of vertical -layers within each category must be the same in the restart file and in -the run restarting from that file. Moreover, significant differences in -the physics, such as the salinity profile, may cause the code to fail -upon restart. Therefore, new model configurations may need to be started -using `runtype` = ‘initial’. Binary restart files that were provided with -CICE v4.1 were made using the BL99 thermodynamics with 4 layers and 5 -thickness categories (`kcatbound` = 0) and therefore can not be used for -the default CICE v5 configuration (7 layers). In addition, CICE’s -default restart file format is now  instead of binary. - -Restarting a run using `runtype` = ‘continue’ requires restart data for -all tracers used in the new run. If tracer restart data is not -available, use `runtype` = ‘initial’, setting `ice\_ic` to the name of the -core restart file and setting to true the namelist restart flags for -each tracer that is available. The unavailable tracers will be -initialized to their default settings. - -On tripole grids, use `restart\_ext` = true when using either binary or -regular (non-PIO) netcdf. - -Provided that the same number of ice layers (default: 4) will be used -for the new runs, it is possible to convert v4.1 restart files to the -new file structure and then to  format. If the same physical -parameterizations are used, the code should be able to execute from -these files. However if different physics is used (for instance, mushy -thermo instead of BL99), the code may still fail. To convert a v4.1 -restart file: - -#. Edit the code **input\_templates/convert\_restarts.f90** for your - model configuration and path names. Compile and run this code to - create a binary restart file that can be read using v5. Copy the - resulting file to the **restart/** subdirectory in your working - directory. - -#. In your working directory, turn off all tracer restart flags in - **ice\_in** and set the following: - - - runtype = ‘initial’ - - - ice\_ic = ‘./restart/[your binary file name]’ - - - restart = .true. - - - use\_restart\_time = .true. - -#. In **CICE\_InitMod.F90**, comment out the call to - restartfile(ice\_ic) and uncomment the call to - restartfile\_v4(ice\_ic) immediately below it. This will read the - v4.1 binary file and write a v5  file containing the same - information. - -If restart files are taking a long time to be written serially (i.e., -not using PIO), see the next section. - -Slow execution --------------- - -On some architectures, underflows (:math:`10^{-300}` for example) are -not flushed to zero automatically. Usually a compiler flag is available -to do this, but if not, try uncommenting the block of code at the end of -subroutine *stress* in **ice\_dyn\_evp.F90** or **ice\_dyn\_eap.F90**. -You will take a hit for the extra computations, but it will not be as -bad as running with the underflows. - -In some configurations, multiple calls to scatter or gather global -variables may overfill MPI’s buffers, causing the code to slow down -(particularly when writing large output files such as restarts). To -remedy this problem, set `BARRIERS yes` in **comp\_ice**. This -synchronizes MPI messages, keeping the buffers in check. - -Debugging hints ---------------- - -Several utilities are available that can be helpful when debugging the -code. Not all of these will work everywhere in the code, due to possible -conflicts in module dependencies. - -*debug\_ice* (**CICE.F90**) - A wrapper for *print\_state* that is easily called from numerous - points during the timestepping loop (see - **CICE\_RunMod.F90\_debug**, which can be substituted for - **CICE\_RunMod.F90**). - -*print\_state* (**ice\_diagnostics.F90**) - Print the ice state and forcing fields for a given grid cell. - -`dbug` = true (**ice\_in**) - Print numerous diagnostic quantities. - -`print\_global` (**ice\_in**) - If true, compute and print numerous global sums for energy and mass - balance analysis. This option can significantly degrade code - efficiency. - -`print\_points` (**ice\_in**) - If true, print numerous diagnostic quantities for two grid cells, - one near the north pole and one in the Weddell Sea. This utility - also provides the local grid indices and block and processor numbers - (`ip`, `jp`, `iblkp`, `mtask`) for these points, which can be used in - conjunction with `check\_step`, to call *print\_state*. These flags - are set in **ice\_diagnostics.F90**. This option can be fairly slow, - due to gathering data from processors. - -*global\_minval, global\_maxval, global\_sum* (**ice\_global\_reductions.F90**) - Compute and print the minimum and maximum values for an individual - real array, or its global sum. - -Known bugs ----------- - -#. Fluxes sent to the CESM coupler may have incorrect values in grid - cells that change from an ice-free state to having ice during the - given time step, or vice versa, due to scaling by the ice area. The - authors of the CESM flux coupler insist on the area scaling so that - the ice and land models are treated consistently in the coupler (but - note that the land area does not suddenly become zero in a grid cell, - as does the ice area). - -#. With the old CCSM radiative scheme (`shortwave` = ‘default’ or - ‘ccsm3’), a sizable fraction (more than 10%) of the total shortwave - radiation is absorbed at the surface but should be penetrating into - the ice interior instead. This is due to use of the aggregated, - effective albedo rather than the bare ice albedo when - `snowpatch` :math:`< 1`. - -#. The date-of-onset diagnostic variables, `melt\_onset` and `frz\_onset`, - are not included in the core restart file, and therefore may be - incorrect for the current year if the run is restarted after Jan 1. - Also, these variables were implemented with the Arctic in mind and - may be incorrect for the Antarctic. - -#. The single-processor *system\_clock* time may give erratic results on - some architectures. - -#. History files that contain time averaged data (`hist\_avg` = true in - **ice\_in**) will be incorrect if restarting from midway through an - averaging period. - -#. In stand-alone runs, restarts from the end of `ycycle` will not be - exact. - -#. Using the same frequency twice in `histfreq` will have unexpected - consequences and causes the code to abort. - -#. Latitude and longitude fields in the history output may be wrong when - using padding. - -Interpretation of albedos -------------------------- - -The snow-and-ice albedo, `albsni`, and diagnostic albedos `albice`, `albsno`, -and `albpnd` are merged over categories but not scaled (divided) by the -total ice area. (This is a change from CICE v4.1 for `albsni`.) The latter -three history variables represent completely bare or completely snow- or -melt-pond-covered ice; that is, they do not take into account the snow -or melt pond fraction (`albsni` does, as does the code itself during -thermodyamic computations). This is to facilitate comparison with -typical values in measurements or other albedo parameterizations. The -melt pond albedo `albpnd` is only computed for the Delta-Eddington -shortwave case. - -With the Delta-Eddington parameterization, the albedo depends on the -cosine of the zenith angle (:math:`\cos\varphi`, `coszen`) and is zero if -the sun is below the horizon (:math:`\cos\varphi < 0`). Therefore -time-averaged albedo fields would be low if a diurnal solar cycle is -used, because zero values would be included in the average for half of -each 24-hour period. To rectify this, a separate counter is used for the -averaging that is incremented only when :math:`\cos\varphi > 0`. The -albedos will still be zero in the dark, polar winter hemisphere. - -Proliferating subprocess parameterizations ------------------------------------------- - -With the addition of several alternative parameterizations for sea ice -processes, a number of subprocesses now appear in multiple parts of the -code with differing descriptions. For instance, sea ice porosity and -permeability, along with associated flushing and flooding, are -calculated separately for mushy thermodynamics, topo and level-ice melt -ponds, and for the brine height tracer, each employing its own -equations. Likewise, the BL99 and mushy thermodynamics compute freeboard -and snow–ice formation differently, and the topo and level-ice melt pond -schemes both allow fresh ice to grow atop melt ponds, using slightly -different formulations for Stefan freezing. These various process -parameterizations will be compared and their subprocess descriptions -possibly unified in the future. diff --git a/doc/source/cice_7_testing.rst b/doc/source/cice_7_testing.rst deleted file mode 100644 index 079345887..000000000 --- a/doc/source/cice_7_testing.rst +++ /dev/null @@ -1,242 +0,0 @@ -.. - -Testing CICE -================ - -Version 6, August 2017 -This documents how to use the testing features developed for the -CICE Consortium CICE sea ice model. - -.. _basic: - -Individual tests and test suites --------------------------------- - -The CICE scripts support both setup of individual tests as well as test suites. Individual -tests are run from the command line like - - > create.case -t smoke -m wolf -g gx3 -p 8x2 -s diag1,run5day -testid myid - -where -m designates a specific machine. Test suites are multiple tests that are specified in -an input file and are started on the command line like - - > create.case -ts base_suite -m wolf -testid myid - -create.case with -t or -ts require a testid to uniquely name test directories. The format -of the case directory name for a test will always be -${machine}_${test}_${grid}_${pes}_${soptions}.${testid} - -To build and run a test, the process is the same as a case, - cd into the test directory, - - run cice.build - - run cice.submit - -The test results will be generated in a local file called "test_output". - -When running a test suite, the create.case command line automatically generates all the tests -under a directory names ${test_suite}.${testid}. It then automatically builds and submits all -tests. When the tests are complete, run the results.csh script to see the results from all the -tests. - -Tests are defined under configuration/scripts/tests. The tests currently supported are: - smoke - Runs the model for default length. The length and options can - be set with the -s commmand line option. The test passes if the - model completes successfully. - restart - Runs the model for 10 days, writing a restart file at day 5 and - again at day 10. Runs the model a second time starting from the - day 5 restart and writing a restart at day 10 of the model run. - The test passes if both the 10 day and 5 day restart run complete and - if the restart files at day 10 from both runs are bit-for-bit identical. - -Please run './create.case -h' for additional details. - -.. _additional: - -Additional testing options --------------------------- - -There are several additional options on the create.case command line for testing that -provide the ability to regression test and compare tests to each other. - - -bd defines a baseline directory where tests can be stored for regression testing - - -bg defines a version name that where the current tests can be saved for regression testing - - -bc defines a version name that the current tests should be compared to for regression testing - - -td provides a way to compare tests with each other - -To use -bg, - > create.case -ts base_suite -m wolf -testid v1 -bg version1 -bd $SCRATCH/CICE_BASELINES -will copy all the results from the test suite to $SCRATCH/CICE_BASELINES/version1. - -To use -bc, - > create.case -ts base_suite -m wolf -testid v2 -bc version1 -bd $SCRATCH/CICE_BASELINES -will compare all the results from this test suite to results saved before in $SCRATCH/CICE_BASELINES/version1. - --bc and -bg can be combined, - >create.case -ts base_suite -m wolf -testid v2 -bg version2 -bc version1 -bd $SCRATCH/CICE_BASELINES -will save the current results to $SCRATCH/CICE_BASELINES/version2 and compare the current results to -results save before in $SCRATCH/CICE_BASELINES/version1. - --bg, -bc, and -bd are used for regression testing. There is a default -bd on each machine. - --td allows a user to compare one test result to another. For instance, - > create.case -t smoke -m wolf -g gx3 -p 8x2 -s run5day -testid t01 - > create.case -t smoke -m wolf -g gx3 -p 4x2 -s run5day -testid t01 -td smoke_gx3_8x2_run5day -An additional check will be done for the second test (because of the -td argument), and it will compare -the output from the first test "smoke_gx3_8x2_run5day" to the output from it's test "smoke_gx3_4x2_run5day" -and generate a result for that. It's important that the first test complete before the second test is -done. Also, the -td option works only if the testid and the machine are the same for the baseline -run and the current run. - -.. _format: - -Test suite format ------------------ - -The format for the test suite file is relatively simple. It is a text file with white space delimited -columns like, - -.. _tab-test: - -.. csv-table:: Table 7 - :header: "#Test", "Grid", "PEs", "Sets", "BFB-compare" - :widths: 7, 7, 7, 15, 15 - - "smoke", "gx3", "8x2", "diag1,run5day", "" - "smoke", "gx3", "8x2", "diag24,run1year,medium", "" - "smoke", "gx3", "4x1", "debug,diag1,run5day", "" - "smoke", "gx3", "8x2", "debug,diag1,run5day", "" - "smoke", "gx3", "4x2", "diag1,run5day", "smoke_gx3_8x2_diag1_run5day" - "smoke", "gx3", "4x1", "diag1,run5day,thread", "smoke_gx3_8x2_diag1_run5day" - "smoke", "gx3", "4x1", "diag1,run5day", "smoke_gx3_4x1_diag1_run5day_thread" - "restart", "gx3", "8x1", "", "" - "restart", "gx3", "4x2", "debug", "" - - -The first column is the test name, the second the grid, the third the pe count, the fourth column is -the -s options and the fifth column is the -td argument. The fourth and fifth columns are optional. -The argument to -ts defines which filename to choose and that argument can contain a path. create.case -will also look for the filename in configuration/scripts/tests where some preset test suites are defined. - -Example Tests (Quickstart) --------------------------- - -To generate a baseline dataset for a test case -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -./create.case -t smoke -m wolf -bg cicev6.0.0 -testid t00 - -cd wolf_smoke_gx3_4x1.t00 - -./cice.build - -./cice.submit - -# After job finishes, check output - -cat test_output - - -To run a test case and compare to a baseline dataset -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -./create.case -t smoke -m wolf -bc cicev6.0.0 -testid t01 - -cd wolf_smoke_gx3_4x1.t01 - -./cice.build - -./cice.submit - -# After job finishes, check output - -cat test_output - - -To run a test suite to generate baseline data -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -./create.case -m wolf -ts base_suite -testid t02 -bg cicev6.0.0bs - -cd base_suite.t02 - -# Once all jobs finish, concatenate all output - -./results.csh # All tests results will be stored in results.log - -# To plot a timeseries of "total ice extent", "total ice area", and "total ice volume" - -./timeseries.csh - -ls *.png - - -To run a test suite to compare to baseline data -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -./create.case -m wolf -ts base_suite -testid t03 -bc cicev6.0.0bs - -cd base_suite.t03 - -# Once all jobs finish, concatenate all output - -./results.csh # All tests results will be stored in results.log - -# To plot a timeseries of "total ice extent", "total ice area", and "total ice volume" - -./timeseries.csh - -ls *.png - - -To compare to another test -~~~~~~~~~~~~~~~~~~~~~~~~~~ -`First:` - -./create.case -m wolf -t smoke -testid t01 -p 8x2 - -cd wolf_smoke_gx3_8x2.t01 - -./cice.build - -./cice.submit - -# After job finishes, check output - -cat test_output - -`Then, do the comparison:` - -./create.case -m wolf -t smoke -testid t01 -td smoke_gx3_8x2 -s thread -p 4x1 - -cd wolf_smoke_gx3_4x1_thread.t01 - -./cice.build - -./cice.submit - -# After job finishes, check output - -cat test_output - - -Additional Details ------------------- -- In general, the baseline generation, baseline compare, and test diff are independent. -- Use the '-bd' flag to specify the location where you want the baseline dataset - to be written. Without specifying '-bd', the baseline dataset will be written - to the default baseline directory found in the env. file (ICE_MACHINE_BASELINE). -- If '-bd' is not passed, the scripts will look for baseline datasets in the default - baseline directory found in the env. file (ICE_MACHINE_BASELINE). - If the '-bd' option is passed, the scripts will look for baseline datasets in the - location passed to the -bd argument. -- To generate a baseline dataset for a specific version (for regression testing), - use '-bg '. The scripts will then place the baseline dataset - in $ICE_MACHINE_BASELINE// -- The '-testid' flag allows users to specify a testing id that will be added to the - end of the case directory. For example, "./create.case -m wolf -t smoke -testid t12 -p 4x1" - creates the directory wolf_smoke_gx3_4x1.t12. This flag is REQUIRED if using -t or -ts. diff --git a/doc/source/cice_8_copyright.rst b/doc/source/cice_8_copyright.rst deleted file mode 100644 index 56abd2af8..000000000 --- a/doc/source/cice_8_copyright.rst +++ /dev/null @@ -1,54 +0,0 @@ -Acknowledgments and Copyright -============================= - -====================== -Acknowledgements -====================== -This work has been supported under the Department of Energy’s Climate, -Ocean and Sea Ice Modeling project through the Computer Hardware Applied -Mathematics and Model Physics (CHAMMP) program, Climate Change -Prediction Program (CCPP), Improving the Characterization of Clouds, -Aerosols and the Cryosphere in Climate Models (Cloud-Cryo) program and -Scientific Discovery through Advanced Computing (SCIDAC) program, with -additional support from the T-3 Fluid Dynamics and Solid Mechanics Group -at Los Alamos National Laboratory. Special thanks are due to the -following people: - -- members of the CESM Polar Climate Working Group, including David - Bailey, Cecilia Bitz, Bruce Briegleb, Tony Craig, Marika Holland, - John Dennis, Julie Schramm, Bonnie Light and Phil Jones, - -- Andrew Roberts of the Naval Postgraduate School, - -- David Hebert and Olivier Lecomte for their melt pond work, - -- Jonathan Gregory of the University of Reading and the U.K. MetOffice - for supplying tripole T-fold code and documentation, - -- Alison McLaren, Ann Keen and others working with the Hadley Centre - GCM for testing non-standard model configurations and providing their - code to us, - -- Daniel Feltham and his research group for several new - parameterizations and documentation, - -- Sylvain Bouillon for the revised EVP approach, - -- the many researchers who tested beta versions of CICE 5 and waited - patiently for the official release. - -====================== -Copyright -====================== -© Copyright 2013, LANS LLC. All rights reserved. Unless otherwise -indicated, this information has been authored by an employee or -employees of the Los Alamos National Security, LLC (LANS), operator of -the Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396 -with the U.S. Department of Energy. The U.S. Government has rights to -use, reproduce, and distribute this information. The public may copy and -use this information without charge, provided that this Notice and any -statement of authorship are reproduced on all copies. Neither the -Government nor LANS makes any warranty, express or implied, or assumes -any liability or responsibility for the use of this information. -Beginning with version 4.0, the CICE code carries Los Alamos Software -Release number LA-CC-06-012. diff --git a/doc/source/cice_9_namelist_opts.rst b/doc/source/cice_9_namelist_opts.rst deleted file mode 100644 index 6039a43a8..000000000 --- a/doc/source/cice_9_namelist_opts.rst +++ /dev/null @@ -1,278 +0,0 @@ -.. _tabnamelist: - -Table of namelist options -========================= - -============================== -Comprehensive Namelist Options -============================== - -.. _tab-namelist: - -.. csv-table:: Table 8 - :header: "variable", "options/format", "description", "recommended value" - :widths: 15, 15, 30, 15 - - "*setup_nml*", "", "", "" - "", "", "*Time, Diagnostics*", "" - "``days_per_year``", "``360`` or ``365``", "number of days in a model year", "365" - "``use_leap_years``", "true/false", "if true, include leap days", "" - "``year_init``", "yyyy", "the initial year, if not using restart", "" - "``istep0``", "integer", "initial time step number", "0" - "``dt``", "seconds", "thermodynamics time step length", "3600." - "``npt``", "integer", "total number of time steps to take", "" - "``ndtd``", "integer", "number of dynamics/advection/ridging/steps per thermo timestep", "1" - "", "", "*Initialization/Restarting*", "" - "``runtype``", "``initial``", "start from ``ice_ic``", "" - "", "``continue``", "restart using ``pointer_file``", "" - "``ice_ic``", "``default``", "latitude and sst dependent", "default" - "", "``none``", "no ice", "" - "", "path/file", "restart file name", "" - "``restart``", "true/false", "initialize using restart file", "``.true.``" - "``use_restart_time``", "true/false", "set initial date using restart file", "``.true.``" - "``restart_format``", "nc", "read/write  restart files (use with PIO)", "" - "", "bin", "read/write binary restart files", "" - "``lcdf64``", "true/false", "if true, use 64-bit  format", "" - "``restart_dir``", "path/", "path to restart directory", "" - "``restart_ext``", "true/false", "read/write halo cells in restart files", "" - "``restart_file``", "filename prefix", "output file for restart dump", "‘iced’" - "``pointer_file``", "pointer filename", "contains restart filename", "" - "``dumpfreq``", "``y``", "write restart every ``dumpfreq_n`` years", "y" - "", "``m``", "write restart every ``dumpfreq_n`` months", "" - "", "``d``", "write restart every ``dumpfreq_n`` days", "" - "``dumpfreq_n``", "integer", "frequency restart data is written", "1" - "``dump_last``", "true/false", "if true, write restart on last time step of simulation", "" - "", "", "*Model Output*", "" - "``bfbflag``", "true/false", "for bit-for-bit diagnostic output", "" - "``diagfreq``", "integer", "frequency of diagnostic output in ``dt``", "24" - "", "*e.g.*, 10", "once every 10 time steps", "" - "``diag_type``", "``stdout``", "write diagnostic output to stdout", "" - "", "``file``", "write diagnostic output to file", "" - "``diag_file``", "filename", "diagnostic output file (script may reset)", "" - "``print_global``", "true/false", "print diagnostic data, global sums", "``.false.``" - "``print_points``", "true/false", "print diagnostic data for two grid points", "``.false.``" - "``latpnt``", "real", "latitude of (2) diagnostic points", "" - "``lonpnt``", "real", "longitude of (2) diagnostic points", "" - "``dbug``", "true/false", "if true, write extra diagnostics", "``.false.``" - "``histfreq``", "string array", "defines output frequencies", "" - "", "``y``", "write history every ``histfreq_n`` years", "" - "", "``m``", "write history every ``histfreq_n`` months", "" - "", "``d``", "write history every ``histfreq_n`` days", "" - "", "``h``", "write history every ``histfreq_n`` hours", "" - "", "``1``", "write history every time step", "" - "", "``x``", "unused frequency stream (not written)", "" - "``histfreq_n``", "integer array", "frequency history output is written", "" - "", "0", "do not write to history", "" - "``hist_avg``", "true", "write time-averaged data", "``.true.``" - "", "false", "write snapshots of data", "" - "``history\_dir``", "path/", "path to history output directory", "" - "``history\_file``", "filename prefix", "output file for history", "‘iceh’" - "``write\_ic``", "true/false", "write initial condition", "" - "``incond\_dir``", "path/", "path to initial condition directory", "" - "``incond\_file``", "filename prefix", "output file for initial condition", "‘iceh’" - "``runid``", "string", "label for run (currently CESM only)", "" - "", "", "", "" - "*grid_nml*", "", "", "" - "", "", "*Grid*", "" - "``grid_format``", "``nc``", "read  grid and kmt files", "‘bin’" - "", "``bin``", "read direct access, binary file", "" - "``grid_type``", "``rectangular``", "defined in *rectgrid*", "" - "", "``displaced_pole``", "read from file in *popgrid*", "" - "", "``tripole``", "read from file in *popgrid*", "" - "", "``regional``", "read from file in *popgrid*", "" - "``grid_file``", "filename", "name of grid file to be read", "‘grid’" - "``kmt_file``", "filename", "name of land mask file to be read", "‘kmt’" - "``gridcpl_file``", "filename", "input file for coupling grid info", "" - "``kcatbound``", "``0``", "original category boundary formula", "0" - "", "``1``", "new formula with round numbers", "" - "", "``2``", "WMO standard categories", "" - "", "``-1``", "one category", "" - "", "", "", "" - "*domain_nml*", "", "", "" - "", "", "*Domain*", "" - "``nprocs``", "integer", "number of processors to use", "" - "``processor_shape``", "``slenderX1``", "1 processor in the y direction (tall, thin)", "" - "", "``slenderX2``", "2 processors in the y direction (thin)", "" - "", "``square-ice``", "more processors in x than y, :math:`\sim` square", "" - "", "``square-pop``", "more processors in y than x, :math:`\sim` square", "" - "``distribution_type``", "``cartesian``", "distribute blocks in 2D Cartesian array", "" - "", "``roundrobin``", "1 block per proc until blocks are used", "" - "", "``sectcart``", "blocks distributed to domain quadrants", "" - "", "``sectrobin``", "several blocks per proc until used", "" - "", "``rake``", "redistribute blocks among neighbors", "" - "", "``spacecurve``", "distribute blocks via space-filling curves", "" - "``distribution_weight``", "``block``", "full block size sets ``work_per_block``", "" - "", "``latitude``", "latitude/ocean sets ``work_per_block``", "" - "``ew_boundary_type``", "``cyclic``", "periodic boundary conditions in x-direction", "" - "", "``open``", "Dirichlet boundary conditions in x", "" - "``ns_boundary_type``", "``cyclic``", "periodic boundary conditions in y-direction", "" - "", "``open``", "Dirichlet boundary conditions in y", "" - "", "``tripole``", "U-fold tripole boundary conditions in y", "" - "", "``tripoleT``", "T-fold tripole boundary conditions in y", "" - "``maskhalo_dyn``", "true/false", "mask unused halo cells for dynamics", "" - "``maskhalo_remap``", "true/false", "mask unused halo cells for transport", "" - "``maskhalo_bound``", "true/false", "mask unused halo cells for boundary updates", "" - "", "", "", "" - "*tracer_nml*", "", "", "" - "", "", "*Tracers*", "" - "``tr_iage``", "true/false", "ice age", "" - "``restart_age``", "true/false", "restart tracer values from file", "" - "``tr_FY``", "true/false", "first-year ice area", "" - "``restart_FY``", "true/false", "restart tracer values from file", "" - "``tr_lvl``", "true/false", "level ice area and volume", "" - "``restart_lvl``", "true/false", "restart tracer values from file", "" - "``tr_pond_cesm``", "true/false", "CESM melt ponds", "" - "``restart_pond_cesm``", "true/false", "restart tracer values from file", "" - "``tr_pond_topo``", "true/false", "topo melt ponds", "" - "``restart_pond_topo``", "true/false", "restart tracer values from file", "" - "``tr_pond_lvl``", "true/false", "level-ice melt ponds", "" - "``restart_pond_lvl``", "true/false", "restart tracer values from file", "" - "``tr_aero``", "true/false", "aerosols", "" - "``restart_aero``", "true/false", "restart tracer values from file", "" - "*thermo_nml*", "", "", "" - "", "", "*Thermodynamics*", "" - "``kitd``", "``0``", "delta function ITD approximation", "1" - "", "``1``", "linear remapping ITD approximation", "" - "``ktherm``", "``0``", "zero-layer thermodynamic model", "" - "", "``1``", "Bitz and Lipscomb thermodynamic model", "" - "", "``2``", "mushy-layer thermodynamic model", "" - "``conduct``", "``MU71``", "conductivity :cite:`MU71`", "" - "", "``bubbly``", "conductivity :cite:`PETB07`", "" - "``a_rapid_mode``", "real", "brine channel diameter", "0.5x10\:math:`^{-3}` m" - "``Rac_rapid_mode``", "real", "critical Rayleigh number", "10" - "``aspect_rapid_mode``", "real", "brine convection aspect ratio", "1" - "``dSdt_slow_mode``", "real", "drainage strength parameter", "-1.5x10\:math:`^{-7}` m/s/K" - "``phi_c_slow_mode``", ":math:`0<\phi_c < 1`", "critical liquid fraction", "0.05" - "``phi_i_mushy``", ":math:`0<\phi_i < 1`", "solid fraction at lower boundary", "0.85" - "", "", "", "" - "*dynamics_nml*", "", "", "" - "", "", "*Dynamics*", "" - "``kdyn``", "``0``", "dynamics OFF", "1" - "", "``1``", "EVP dynamics", "" - "", "``2``", "EAP dynamics", "" - "``revised_evp``", "true/false", "use revised EVP formulation", "" - "``ndte``", "integer", "number of EVP subcycles", "120" - "``advection``", "``remap``", "linear remapping advection", "‘remap’" - "", "``upwind``", "donor cell advection", "" - "``kstrength``", "``0``", "ice strength formulation :cite:`Hibler79`", "1" - "", "``1``", "ice strength formulation :cite:`Rothrock75`", "" - "``krdg_partic``", "``0``", "old ridging participation function", "1" - "", "``1``", "new ridging participation function", "" - "``krdg_redist``", "``0``", "old ridging redistribution function", "1" - "", "``1``", "new ridging redistribution function", "" - "``mu_rdg``", "real", "e-folding scale of ridged ice", "" - "``Cf``", "real", "ratio of ridging work to PE change in ridging", "17." - "", "", "", "" - "*shortwave_nml*", "", "", "" - "", "", "*Shortwave*", "" - "``shortwave``", "``default``", "NCAR CCSM3 distribution method", "" - "", "``dEdd``", "Delta-Eddington method", "" - "``albedo_type``", "``default``", "NCAR CCSM3 albedos", "‘default’" - "", "``constant``", "four constant albedos", "" - "``albicev``", ":math:`0<\alpha <1`", "visible ice albedo for thicker ice", "" - "``albicei``", ":math:`0<\alpha <1`", "near infrared ice albedo for thicker ice", "" - "``albsnowv``", ":math:`0<\alpha <1`", "visible, cold snow albedo", "" - "``albsnowi``", ":math:`0<\alpha <1`", "near infrared, cold snow albedo", "" - "``ahmax``", "real", "albedo is constant above this thickness", "0.3 m" - "``R_ice``", "real", "tuning parameter for sea ice albedo from Delta-Eddington shortwave", "" - "``R_pnd``", "real", "... for ponded sea ice albedo …", "" - "``R_snw``", "real", "... for snow (broadband albedo) …", "" - "``dT_mlt``", "real", ":math:`\Delta` temperature per :math:`\Delta` snow grain radius", "" - "``rsnw_mlt``", "real", "maximum melting snow grain radius", "" - "``kalg``", "real", "absorption coefficient for algae", "" - "", "", "", "" - "*ponds_nml*", "", "", "" - "", "", "*Melt Ponds*", "" - "``hp1``", "real", "critical ice lid thickness for topo ponds", "0.01 m" - "``hs0``", "real", "snow depth of transition to bare sea ice", "0.03 m" - "``hs1``", "real", "snow depth of transition to pond ice", "0.03 m" - "``dpscale``", "real", "time scale for flushing in permeable ice", ":math:`1\times 10^{-3}`" - "``frzpnd``", "``hlid``", "Stefan refreezing with pond ice thickness", "‘hlid’" - "", "``cesm``", "CESM refreezing empirical formula", "" - "``rfracmin``", ":math:`0 \le r_{min} \le 1`", "minimum melt water added to ponds", "0.15" - "``rfracmax``", ":math:`0 \le r_{max} \le 1`", "maximum melt water added to ponds", "1.0" - "``pndaspect``", "real", "aspect ratio of pond changes (depth:area)", "0.8" - "", "", "", "" - "*zbgc_nml*", "", "", "" - "", "", "*Biogeochemistry*", "" - "``tr_brine``", "true/false", "brine height tracer", "" - "``restart_hbrine``", "true/false", "restart tracer values from file", "" - "``skl_bgc``", "true/false", "biogeochemistry", "" - "``bgc_flux_type``", "``Jin2006``", "ice–ocean flux velocity of :cite:`JDWSTWLG06`", "" - "", "``constant``", "constant ice–ocean flux velocity", "" - "``restart_bgc``", "true/false", "restart tracer values from file", "" - "``restore_bgc``", "true/false", "restore nitrate/silicate to data", "" - "``bgc_data_dir``", "path/", "data directory for bgc", "" - "``sil_data_type``", "``default``", "default forcing value for silicate", "" - "", "``clim``", "silicate forcing from ocean climatology :cite:`GLBA06`", "" - "``nit_data_type``", "``default``", "default forcing value for nitrate", "" - "", "``clim``", "nitrate forcing from ocean climatology :cite:`GLBA06`", "" - "", "``sss``", "nitrate forcing equals salinity", "" - "``tr_bgc_C_sk``", "true/false", "algal carbon tracer", "" - "``tr_bgc_chl_sk``", "true/false", "algal chlorophyll tracer", "" - "``tr_bgc_Am_sk``", "true/false", "ammonium tracer", "" - "``tr_bgc_Sil_sk``", "true/false", "silicate tracer", "" - "``tr_bgc_DMSPp_sk``", "true/false", "particulate DMSP tracer", "" - "``tr_bgc_DMSPd_sk``", "true/false", "dissolved DMSP tracer", "" - "``tr_bgc_DMS_sk``", "true/false", "DMS tracer", "" - "``phi_snow``", "real", "snow porosity for brine height tracer", "" - "", "", "", "" - "*forcing_nml*", "", "", "" - "", "", "*Forcing*", "" - "``formdrag``", "true/false", "calculate form drag", "" - "``atmbndy``", "``default``", "stability-based boundary layer", "‘default’" - "", "``constant``", "bulk transfer coefficients", "" - "``fyear_init``", "yyyy", "first year of atmospheric forcing data", "" - "``ycycle``", "integer", "number of years in forcing data cycle", "" - "``atm_data_format``", "``nc``", "read  atmo forcing files", "" - "", "``bin``", "read direct access, binary files", "" - "``atm_data_type``", "``default``", "constant values defined in the code", "" - "", "``LYq``", "AOMIP/Large-Yeager forcing data", "" - "", "``monthly``", "monthly forcing data", "" - "", "``ncar``", "NCAR bulk forcing data", "" - "", "``oned``", "column forcing data", "" - "``atm_data_dir``", "path/", "path to atmospheric forcing data directory", "" - "``calc_strair``", "true", "calculate wind stress and speed", "" - "", "false", "read wind stress and speed from files", "" - "``highfreq``", "true/false", "high-frequency atmo coupling", "" - "``natmiter``", "integer", "number of atmo boundary layer iterations", "" - "``calc_Tsfc``", "true/false", "calculate surface temperature", "``.true.``" - "``precip_units``", "``mks``", "liquid precipitation data units", "" - "", "``mm_per_month``", "", "" - "", "``mm_per_sec``", "(same as MKS units)", "" - "``tfrz_option``", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8\degC`)", "" - "", "``linear_salt``", "linear function of salinity (ktherm=1)", "" - "", "``mushy_layer``", "matches mushy-layer thermo (ktherm=2)", "" - "``ustar_min``", "real", "minimum value of ocean friction velocity", "0.0005 m/s" - "``fbot_xfer_type``", "``constant``", "constant ocean heat transfer coefficient", "" - "", "``Cdn\_ocn``", "variable ocean heat transfer coefficient", "" - "``update_ocn_f``", "true", "include frazil water/salt fluxes in ocn fluxes", "" - "", "false", "do not include (when coupling with POP)", "" - "``l_mpond_fresh``", "true", "retain (topo) pond water until ponds drain", "" - "", "false", "release (topo) pond water immediately to ocean", "" - "``oceanmixed_ice``", "true/false", "active ocean mixed layer calculation", "``.true.`` (if uncoupled)" - "``ocn_data_format``", "``nc``", "read  ocean forcing files", "" - "", "``bin``", "read direct access, binary files", "" - "``sss_data_type``", "``default``", "constant values defined in the code", "" - "", "``clim``", "climatological data", "" - "", "``near``", "POP ocean forcing data", "" - "``sst_data_type``", "``default``", "constant values defined in the code", "" - "", "``clim``", "climatological data", "" - "", "``ncar``", "POP ocean forcing data", "" - "``ocn_data_dir``", "path/", "path to oceanic forcing data directory", "" - "``oceanmixed_file``", "filename", "data file containing ocean forcing data", "" - "``restore_sst``", "true/false", "restore sst to data", "" - "``trestore``", "integer", "sst restoring time scale (days)", "" - "``restore_ice``", "true/false", "restore ice state along lateral boundaries", "" - "", "", "", "" - "*icefields_tracer_nml*", "", "", "" - "", "", "*History Fields*", "" - "``f_``", "string", "frequency units for writing ```` to history", "" - "", "``y``", "write history every ``histfreq_n`` years", "" - "", "``m``", "write history every ``histfreq_n`` months", "" - "", "``d``", "write history every ``histfreq_n`` days", "" - "", "``h``", "write history every ``histfreq_n`` hours", "" - "", "``1``", "write history every time step", "" - "", "``x``", "do not write ```` to history", "" - "", "``md``", "*e.g.,* write both monthly and daily files", "" - "``f__ai``", "", "grid cell average of ```` (:math:`\times a_i`)", "" diff --git a/doc/source/conf.py b/doc/source/conf.py index 4751a71c7..a9c7291d2 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -116,12 +116,12 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. -html_theme = 'bizstyle' +html_theme = 'classic' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. -#html_theme_options = {} +html_theme_options = {"stickysidebar": "true"} # Add any paths that contain custom themes here, relative to this directory. #html_theme_path = [] diff --git a/doc/source/index.rst b/doc/source/index.rst index 6719d6d25..3c73abbeb 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -11,19 +11,13 @@ Table of Contents: ------------------ .. toctree:: - :maxdepth: 3 + :maxdepth: 5 :numbered: cice_1_introduction.rst - cice_2_quick_start.rst - cice_3_coupling.rst - cice_4_model_components.rst - cice_5_numerical_implementation.rst - cice_6_troubleshooting.rst - cice_7_testing.rst - cice_8_copyright.rst - cice_9_namelist_opts.rst - cice_10_index.rst + cice_2_science_guide.rst + cice_3_user_guide.rst + cice_4_index.rst zreferences.rst Useful tools diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 84d9a9d1a..e572f7d9f 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -207,7 +207,7 @@ @Article{RM88 author = "A. Rosati and K. Miyakoda", title = "{A general circulation model for upper ocean simulation}", journal = JPO, - year = "1999", + year = "1988", volume = {18}, pages = {1601-1626}, url = {https://doi.org/10.1175/1520-0485(1988)018<1601:AGCMFU>2.0.CO;2} @@ -369,6 +369,14 @@ @Article{JAM99 pages = {7785-7806}, url = {http://dx.doi.org/10.1029/1999JC900011} } +@Book{vSZ99 + author = "H. von Storch and F.W. Zwiers", + title = "{Statistical Analysis in Climate Research}", + publisher = "Cambridge University Press", + note = "Cambridge, UK", + year = "1999", + pages = {484 pp}, +} @Article{DB00 author = "J.K. Dukowicz and J.R. Baumgardner", title = "{Incremental remapping as a transport/advection algorithm}", @@ -553,6 +561,14 @@ @Article{WF06 pages = {22-32}, url = {http://dx.doi.org/10.1016/j.jnnfm.2006.05.001} } +@Book{Wilks06 + author = "D.S. Wilks", + title = "{Statistical methods in the atmospheric sciences}", + publisher = "Academic Press", + note = "2nd ed.", + year = "2006", + pages = {627 pp}, +} @Manual{BL07 author = "B.P. Briegleb and B. Light", title = "{A Delta-Eddington multiple scattering parameterization for solar radiation in the sea ice component of the Community Climate System Model}",