Home

CICE: the Los Alamos Sea Ice Model Documentation and

1. length of eastern edge Ay of T cell length of northern edge Ax of T cell length of southern edge Ax of T cell length of western edge of Ay T cell global domain location for each grid cell fraction of penetrating visible solar radiation block on which to write debugging data Cartesian i j position of block e choice of initial conditions see Table 4 reference salinity for ice ocean exchanges number of grid cells with specified property for vectorization 1ce surface TOUSHNESS aia ada senden 63 9 80616 m s 0 15 0 05 m 0 01 m 20 m 0 2 m 0 005 m 1 x107 m 1 x10 4m 0 03 m 53538585R 5 0 70 4 psu lO m 64 icetmask iceumask idate idateO ierr 1G hi i j lo ilyrl ilyrn incond_dir incond_file int_kind integral_order ip jp istep istep0 istepl Iswabs K kappan kappav kcatbound kdyn kg_to_g kice kimin kitd kmt file krdg partic krdg redist kseaice ksno kstrength L 1_brine 1_conservation_check 1 fixed_area latpt latt u _bounds Lfresh lhcoef Imask_n s local_id log kind Index of primary variables and parameters ice extent mask T cell ice extent mask U cell the date at the end of the current time step yyyymmdd initial date general use error flag last ij index of physical domain local first iG index of physical domain local index of the top layer in each cat for ei
2. 0 e eee eee eee day of year that surface melt begins the month number previous month number m per s to cm per day conversion nn local processor number that writes debugging data e e folding scale of ridged ice task ID for the current processor number of blocks on current processor total number of blocks in decomposition total number of blocks in x y direction number of ice categories 0 066 cee cece eee eee e number of subcycles ne e number of dynamics advection steps under thermo flag for beginning new day flag for beginning new hour flag for beginning new month flag for beginning new year number of rows of ghost cells surrounding each subdomain number of groups of flux triangles in remapping northern latitude of artificial mask edge number of ice layers in each category e total number of processors e total number of time steps dt 65 degrees E degrees E 2 835x 10 I kg 2 501 x 10 J kg 0 01 kg m 100 1 6x 1078 m s deg 1 36 1x1076 8 64x 10 30 S 66 ns_boundary_type nslyr nspint nt_ trcr ntilyr ntrace ntrcr ntslyr nu_diag nu_dump nu_dump age nu_dump_pond nu_forcing nu_grid nu_hdr nu_history nu_kmt nu_nml nu_restart nu_restart_age nu_restart_pond nu_rst_pointer nx y _block nx y global nyr O oceanmixed_file oceanmixed_ice ocn_data_dir ocn_data_format omega opening
3. distribution original round WMO kcatbound 0 1 2 No 5 5 5 6 7 category lower bound m 1 0 00 0 00 0 00 0 00 0 00 2 0 64 0 60 0 30 0 15 0 10 3 1 39 1 40 0 70 0 30 0 15 4 2 47 2 40 1 20 0 70 0 30 5 4 57 3 60 2 00 1 20 0 70 6 2 00 1 20 7 2 00 Table 2 Lower boundary values for thickness categories in meters for the three distribution options kcatbound In the WMO case the distribution used depends on the number of categories used ice lies in each thickness range Thus the basic problem in sea ice modeling is to describe the evolution of the ice thickness distribution ITD in time and space The fundamental equation solved by CICE is 43 7 V gu 2 ito 2 where u is the horizontal ice velocity V 4 37 f is the rate of thermodynamic ice growth y is a ridging redistribution function and g is the ice thickness distribution function We define g x h t dh as the fractional area covered by ice in the thickness range h h dh at a given time and location Equation 2 is solved by partitioning the ice pack in each grid cell into discrete thickness categories The number of categories can be set by the user with a default value Nc 5 Five categories plus open water are generally sufficient to simulate the annual cycles of ice thickness ice strength and surface fluxes 3 24 Each category n has lower thickness bound H _ and upper bound H The lower bound of the thinnest ice
4. Ayy where az Qq 0a 0x and ay aq Oa Oy are the limited gradients in the x and y directions respec tively and the components of T T f4 xdA A and y J ydA A are evaluated using the triangle integration formulas described in Section 3 1 3 These means along with higher order means such as 22 zy and y are computed once and stored Next consider the ice and snow thickness and enthalpy fields Thickness is analogous to the tracer concentration T in 7 but there is no analog in 7 to the enthalpy The reconstructed ice or snow thickness h r and enthalpy q r must satisfy il ahdA A ahqdA A II al gt D 13 II el gt gt D 14 where h h is the thickness at the center of ice area and q q f is the enthalpy at the center of ice or snow volume Equations 13 and 14 are satisfied when h r and q r are given by h r h an Vh r f 15 q r aq Vq x 16 where ay and a are limiting coefficients The center of ice area r and the center of ice or snow volume f are given by 1 A dA f oraa 1 f ahrdA ahA Ja Evaluating the integrals we find that the components of f are BL AcE aga ayry y a Rr A OpEY ayy a lt Q gt and the components of f are T coe C3TY Car c522y cory ah 8 y Y CTG cgy eur y cry coy ah Kay y Horizontal transport 13
5. Tf Tffresh Tfrzpt time time _forc Timelt tinyarea TLAT TLON tmask tmass Tmin Tmlt Tocnfrz triage tr_pond trer trer_depend Tref trestore tripole tripoleT Tsf_errmax Tsfc n Tsmelt TTTice TTTocn U uarea uarear uatm ULAT ULON umask ice ocean stress in the x y direction U cell ice ocean stress x y dir T cell surface stress due to sea surface slope incoming shortwave radiation visible near IR direct diffuse air temperature at 10 m area of T cell area of northern hemisphere T cells l tarea area of southern hemisphere T cells absolute day number freezing temperature freezing temp of fresh ice oenen e type of freezing temperature constant or linear_S total elapsed time time of last forcing update melting temperature of ice top surface puny tarea latitude of cell center longitude of cell center land boundary mask thickness T cell total mass of ice and snow minimum allowed internal temperature melting temperature of ice temperature of constant freezing point parameterization e if true use ice age tracer e if true use explicit melt ponds ice tracers tracer dependency on basic state variables 2m atmospheric reference temperature e sst restoring time scale if true block lies along tripole boundary if true tripole boundary is T fold if false U fold max allowed T fc error thermodynamics
6. and is provided here for direct communication compatibility with POP The latitude Performance 47 CICE4 0 timings 320x384 grid 240 dt January X NPX x 1 ome NPX x 2 Sa Me square ice 100 x thy me RA eee square pop al n ke E o RF 4 o Xx n Ps A ee cartesian E x rake weight block size 7 rake weight latitude nocn i 1 fi 1 1 1 1 fi L i 10 100 number of processors 220 A 4 lis x 4 E 160 x 96 B 200 180 al if rocessors 2 x 4 4 E bien size 20 x 24 x 4 S 160 4x2 zal 8 C O 80 x 192 ken 4 8 x 1 10x12 El 140 E 40 x 384 i X L 2 8x 1 Xx 7 pa 20 x 24 E 10 x 12 J7 120 w el E 20 x 24 X 7 1000 oy q 4 o y yl 10 100 1000 number of blocks Figure 7 a CICE 4 0 timings for 8 16 32 and 64 processors of the Linux cluster coyote after 240 1 hr time steps on a 320x384 displaced pole grid b 8 processor timings for various processor shapes indicated by the upper pair of numbers eg 8x1 slenderX 1 and block sizes 48 Numerical implementation option weights the blocks based on latitude and the number of ocean grid cells they contain The rake distribution type is initialized as a standard Cartesian distribution Using the work per block estimates blocks are raked onto neighboring processors as needed to improve load balancing charact
7. cm to meters conversion enn value for constant albedo parameterization value for constant albedo parameterization e conductivity parameterization basal ice growth cosine of the turning angle in water cosine of the zenith angle proportionality constant for potential energy specific heat of air eee eee eee specific heat of fresh ice specific heat of sea water 0000 specific heat of water vapor diffuse fresnel reflectivity above diffuse fresnel reflectivity below fraction of shear energy contributing to ridging constant in Hibler ice strength formula ice area tendency due to dynamics transport ice area tendency due to thermodynamics see ice_shortwave F90 see ice_shortwave F90 see ice_shortwave F90 rate of fractional area loss by ridging ice rate of fractional area gain by new ridges 1080 256 0 01 0 70 0 81 m 1 kg m s 1005 0 J kg K 2106 J kg K 4218 J kg K 1 81x 10 J kg K 0 063 0 455 0 25 20 1 s 1 s 0 075 0 100 0 150 1 s 1 s Index of primary variables and parameters daymo daycal days_per_year dbl_kind dbug Delta depressT diag file diag type diagfreq distrb_info distribution type distribution_weight divu divu_adv dragio dragw dt dT_mlt dte dtei dump file dumpfreq dumpfreq_n dxt dx
8. coupling see flux coupler currents ocean 4 7 damping timescale 25 Delta Eddington see radiation density atmosphere 4 6 29 ice or snow 29 33 ocean 7 29 diagnostics 6 36 43 51 distribution block see blocks thickness see thickness distribution dynamics elastic viscous plastic see elastic viscous plastic ridging see ridging transport see transport ECMWE 57 elastic 71 viscous plastic dynamics 2 9 24 26 36 40 42 44 waves 24 25 42 energy see enthalpy enthalpy 8 11 13 26 32 33 evaporation 4 33 EVP see elastic viscous plastic dynamics flux coupler 3 7 9 34 41 44 50 51 fraction ice see ice fraction frazil 7 freeboard 34 freezing potential 4 7 fresh water flux 4 5 7 grid 14 25 35 36 3740 41 Hadley Centre 2 5 52 halo 38 height reference 4 6 sea surface 7 history 36 37 42 43 45 49 51 humidity reference 4 6 specific 4 6 28 ice see individual variables age 9 10 36 49 fraction 3 5 8 11 17 19 21 24 51 growth 8 18 20 32 34 level 23 36 49 ice ocean stress 4 5 7 24 initial condition 41 43 internal stress 24 26 37 43 LANL 2 52 53 latent heat 4 6 27 28 33 lateral melt 29 34 leads see open water longwave see radiation longwave Los Alamos National Laboratory see LANL Los Alamos National Security LLC 53 masks 38 40 mechanical distribution see ridging melt pond 5 27 37 49 72
9. 10 km Practical limits may be somewhat less depending on the strength of the atmospheric winds Transport in thickness space imposes a similar restraint on the time step given by the ice growth melt rate and the smallest range of thickness among the categories At lt min AH 2max f where AH is the distance between category boundaries and f is the thermodynamic growth rate For the 5 category ice thickness distribution used as the default in this distribution this is not a stringent limitation At lt 19 4 hr assuming max f 40 cm day The dynamics component is subcycled ndte N times per dynamics time step so that the elastic waves essentially disappear before the next time step The subcycling time step Ate is thus dte dt dyn ndte A second parameter E eyc must be selected which defines the elastic wave damping timescale T described in Section 3 4 as eyc dt_dyn The forcing terms are not updated during the subcycling Given the small step dt e at which the EVP dynamics model is subcycled the elastic parameter F is also limited by stability constraints as discussed in 17 Linear stability analysis for the dynamics component shows that the numerical method is stable as long as the subcycling time step At sufficiently resolves the damping timescale T For the stability analysis we had to make several simplifications of the problem hence the location of the boundary between stable and unstable regions is merely an estima
10. 4255 AAN O d a bata Sak oa a ls Ge he dese sd 4 3 Initialization and coupling oaa ee ee 4 4 Choosing an appropriate time step een 41 4 3 Model output 4 2 35 a aa os Be ba a Bt ne AE ee ee Se eel 42 4 51 History le iria Harde Bh Sod ae Se arr Yao fa as ee o e ada 42 4 3 2 Diagnostic Alesin nn a wt ae ennen etn A ete 43 4 5 3 Restart files arg eers e A be aa eed a A 44 4 6 Execution procedures eeen 44 4 7 hPerlormances ua natto avant al Ae deg ea a eek AE ng hv PR nde el 46 48 Adding things 5e sor ran en meededen a a ie ede ol a Qa es a in de ee as be 48 ABT Nn 48 4 8 2 History Helds sn 20 4 A te ee ew rel iede id Sow ae a 49 4 93 Bracerse un oan ek Gee AE ee Ae vene rl dn ea te an eht el ee 49 5 Troubleshooting 50 Sl Unitialsettip nnn ener en eer ete ae eid Ge en te ot Ere ede 50 52 SSIOWEKECUUON saas aco ks wet Er Oe eh A AED ee 51 3 3 Debugging hints 2 5 00 a oon uhh arab bh ass daha dP Aledo dd God ae sala cd 51 5 4 Known bugs nnn aa en been we A ke A a EEE a Grek 51 5 5 Multi dimensional output ee 52 5 6 Interpretation of albedos eee 52 Acknowledgments and Copyright 52 Table of namelist options 54 Index of primary variables and parameters 59 General Index 71 Bibliography 74 1 Introduction The Los Alamos sea ice model CICE is the result of an effort to develop a computationally efficient sea ice component for a fully coupled atmosphere ice ocea
11. M M Holland and B P Briegleb Impact of sea ice melt ponds and solar radiation param eterization in the Community Climate System Model Manuscript in preparation 2008 C M Bitz M M Holland A J Weaver and M Eby Simulating the ice thickness distribution in a coupled climate model J Geophys Res Oceans 106 2441 2463 2001 C M Bitz and W H Lipscomb An energy conserving thermodynamic sea ice model for climate study J Geophys Res Oceans 104 15669 15677 1999 B P Briegleb and B Light A Delta Eddington multiple scattering parameterization for solar radiation in the sea ice component of the Community Climate System Model NCAR Tech Note NCAR TN 4724 STR National Center for Atmospheric Research 2007 W M Connolley J M Gregory E C Hunke and A J McLaren On the consistent scaling of terms in the sea ice dynamics equation J Phys Oceanogr 34 1776 1780 2004 J K Dukowicz and J R Baumgardner Incremental remapping as a transport advection algorithm J Comput Phys 160 318 335 2000 J K Dukowicz R D Smith and R C Malone A reformulation and implementation of the Bryan Cox Semtner ocean model on the connection machine J Atmos Oceanic Technol 10 195 208 1993 J K Dukowicz R D Smith and R C Malone Implicit free surface method for the Bryan Cox Semtner ocean model J Geophys Res Oceans 99 7991 8014 1994 E E Ebert J L Schramm and J A Curry Dispositio
12. P p001 pol p027 p05 p055 p111 p15 p166 p2 p222 Index of primary variables and parameters e type of north south boundary condition number of snow layers in each category number of solar spectral intervals tracer index sum of number of ice layers in all categories number of fields being transported number of tracers transported in remapping sum of number of snow layers in all categories unit number for diagnostics output file unit number for dump file for restarting unit number for age dump file for restarting unit number for melt pond dump file for restarting unit number for forcing data file unit number for grid file unit number for binary history header file unit number for history file unit number for land mask file unit number for namelist input file unit number for restart input file unit number for age restart input file unit number for melt pond restart input file unit number for pointer to latest restart file total number of gridpoints on block in x y direction number of physical gridpoints in x y direction global domain year number e data file containing ocean forcing data e if true use internal ocean mixed layer e directory for ocean forcing data e format of ocean forcing files angular velocity of Earth ee rate of ice opening due to divergence and shear 1 1000 1 100 1 36 1 20 1 18 1 10 1 9 15 100 1 6 1 5 2 9 7 292 x10 rad s l s Index of primary variables and p
13. TBE gt 36 2 aPn Ge Gn Gi Sa Ge where ap is the ratio of the ice area ridging or open water area closing in category n to the total area ridging and closing and Gn is the total fractional ice area in categories O to n Equation 36 applies to categories with Gn lt G If Gn 1 lt G lt Gn then 36 is valid with G replacing Gn and if Gn 1 gt G then apn 0 If the open water fraction ag gt G no ice can ridge because ridging simply reduces the area of open water As in 43 we set G 0 15 If the spatial resolution is too fine for a given time step At the weighting function 35 can promote numerical instability For At 1 hour resolutions finer than Ax 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 At 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 A more stable weighting function was suggested by 26 exp G h a b h 37 WE p 1 07 qa When integrated between category boundaries 37 implies as exp G 1 a exp G a 38 1 exp 1 a This weighting function is used if krdg_partic 1 in the namelist From 37 the mean value of G for ice p
14. Tracer monotonicity is ensured because T fAahdA faadA JaahqdA E AA JxahdA where h and q are the new time thickness and enthalpy given by integrating the old time ice area volume and energy over a Lagrangian departure region with area A That is the new time thickness and enthalpy are weighted averages over old time values with non negative weights a and ah Thus the new time values must lie between the maximum and minimum of the old time values 3 2 Transport in thickness space Next we solve the equation for ice transport in thickness space due to thermodynamic growth and melt g 0 DE T ant 0 26 which is obtained from 2 by neglecting the first and third terms on the right hand side We use the remap ping method of 24 in which thickness categories are represented as Lagrangian grid cells whose bound aries are projected forward in time The thickness distribution function g is approximated as a linear function of 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 We first compute the displacement of category boundaries in thickness space Assume that at time m the ice areas a and mean ice thicknesses hij are known for each thickness category For now we omit the subscript that distinguishes ice from snow W
15. model suffers a linearization error associated with the viscosities which are lagged over the time step 15 This led to principal stress states that were widely scattered outside the elliptical yield curve in both models 20 We have addressed this problem by updating the viscosities during the subcycling so that the entire dynamics component is subcycled within the time step Taken alone this change would require an increased number operations to compute the viscosities However the new dynamics component is roughly as efficient as the earlier version due to a change in the definition of the elastic parameter E E is now defined in terms of a damping timescale T for elastic waves At lt T lt At as where T E At and Eo eyc is a tunable parameter less than one as before The stress equations 50 52 become Oo O1 P P m5 Dp o Tar or 2TA P Oos elos P pd NES y at OT MAA 0012 e o12 P SRP ECB oe Ot 2T ATA All coefficients on the left hand side are constant except for P which changes only on the longer time step At 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 E T and At are discussed in Section 4 4 A different discretization of the stress terms 00 0x in the momentum equation is now used which enabled the discrete equations to be derived from the continuous
16. restarts from the end of ycycle will not be exact 7 Using the same frequency twice in hist freq will have unexpected consequences and causes the code to abort 5 5 Multi dimensional output The column configuration of the model has not been thoroughly tested it is being provided for a few users who are beginning to look more closely at thermodynamic and biogeochemical processes within the sea ice column and problems with the configuration are likely to arise as the code is exercised An issue that has already arisen in this context is the format of the multi dimensional history variables Even though each history file includes only a single time slice or average all netcdf history variables include the time dimen sion for use with external post processing software such as NCO When the time dimension is included for the four dimensional variables x y z and categories strictly speaking time makes them 5D the Ferret software package misinterprets the data In order to look at these variables using Ferret uncomment the lines indicated by the tag ferret in ice_history F90 and comment out the lines that they replace 5 6 Interpretation of albedos Interpretation of the albedos in the history output file is a little tricky The snow and ice albedo albsni is merged across categories and then scaled divided by the total ice area as with other variables sent to the CCSM coupler The diagnostic albedos albice albsno and albpnd are m
17. setting kstrength 1 gives an ice strength closely related to the ridging scheme Following 35 the strength is assumed proportional to the change in ice potential energy A Ep per unit area of compressive deformation Given uniform ridge ITDs krdg_redist 0 we have Ne max 3 miny3 n H gt Ap P apn h2 2 Hn um 47 oe apn hy E zine Hem 47 where Cp g 2 pi Pw Pw pi B Rtot Rnet gt 1 from 45 and Cf is an empirical parameter that accounts for frictional energy dissipation Following 11 we set Cy 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 is included because app 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 Nc P C1O p apn hh 2 Hain 2Hmin 222 48 n 1 n The energy based ice strength given by 47 or 48 is more physically realistic than the strength given by 46 However use of 46 is less likely to allow numerical instability at a given resolution and time step See 26 for more details 24 Model components 3 4 Dynamics The elastic viscous plastic EVP model represents a modification of the standard viscous plastic VP
18. temperature of ice snow top surface in category n melting temperature of snow top surface for saturated specific humidity over ice 06 for saturated specific humidity over ocean area of U cell 1 uarea wind velocity in the x direction latitude of U cell centers longitude of U cell centers land boundary mask velocity U cell 69 N m N m N m W m s s 0 C ime radians radians kg m 100 C 1 8 C K days 5 x10 fdeg C 0 C 5897 8 K 5107 4 K m2 m2 m s radians radians 70 umax_stab umin uocn update_ocn_f ustar_min uvel uvm V vatm vice n vicen_init vocn vonkar vsno n vvel W warmice warmsno wind work_g1 23 work_gi4 work gi8 work_gr work gr3 work 1 2 worka bcd write_history write_ic write_restart Y ycycle yday year_init Z zlvl zref zTrf zvir Index of primary variables and parameters ice speed threashold diagnostics min wind speed for turbulent fluxes ocean current in the x direction e if true include frazil ice fluxes in ocean flux fields e minimum friction velocity under ice x component of ice velocity land boundary mask velocity U cell wind velocity in the y direction volume per unit area of ice in category n ice volume at beginning of timestep ocean current in the y direction von Karman constant volume per unit area of snow in category n y co
19. zp PERA by y PER z A 2 2 2 2 17 where Lis an edge length and the indices N S E W denote compass directions Equation 17 is equivalent to the divergence computed in the EVP dynamics Section 3 4 In general the fluxes in this expression are not equal to those implied by the above scheme for locating departure regions For some applications it may be desirable to prescribe the divergence by prescribing the area of the departure region for each edge This can be done in CICE 4 0 by setting 1_fixed_area true in ice transport driver F90 and passing the prescribed departure areas edgearea_e and edgearea_n into the remapping routine An extra triangle is then constructed for each departure region to ensure that the total area is equal to the prescribed value This idea was suggested and first implemented by Mats Bentsen of the Nansen Environmental and Remote Sensing Center Norway who applied an earlier version of the CICE remapping scheme to an ocean model The implementation in CICE 4 0 is somewhat more general allowing for departure regions lying on both sides of a cell edge The extra triangle is constrained to lie in one but not both of the grid cells that share the edge Since this option has yet to be fully tested in CICE the current default is 1_fixed_area false We made one other change in the scheme of 7 for locating triangles In their paper departure points are defined by projecting cell corner velocities directly back
20. 2003 20 E C Hunke and Y Zhang A comparison of sea ice dynamics models at high resolution Mon Wea Rev 127 396 408 1999 21 R E Jordan E L Andreas and A P Makshtas Heat budget of snow covered sea ice at North Pole 4 J Geophys Res Oceans 104 7785 7806 1999 22 B G Kauffman and W G Large The CCSM coupler version 5 0 1 Technical note National Center for Atmospheric Research August 2002 http www ccsm ucar edu models 23 W H Lipscomb Modeling the Thickness Distribution of Arctic Sea Ice PhD thesis Dept of Atmo spheric Sciences Univ of Washington Seattle 1998 24 W H Lipscomb Remapping the thickness distribution in sea ice models J Geophys Res Oceans 106 13 989 14 000 2001 25 W H Lipscomb and E C Hunke Modeling sea ice transport using incremental remapping Mon Wea Rev 132 1341 1354 2004 26 W H Lipscomb E C Hunke W Maslowski and J Jakacki Improving ridging schemes for high resolution sea ice models J Geophys Res Oceans 112 C03S91 doi 10 1029 2005JC003355 2007 27 W Maslowski D Marble W Walczowski U Schauer J L Clement and A J Semtner On climato logical mass heat and salt transports through the Barents Sea and Fram Strait from a pan Arctic cou pled ice ocean model simulation J Geophys Res 109 C03032 doi 10 1029 2001JC001039 2004 28 G A Maykut Large scale heat exchange and ice production in the central Arctic
21. 3 1 1 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 x and y For each field we compute the value at the cell center i e at the origin of a 2D Cartesian coordinate system defined for that grid cell along with gradients in the x and y directions The gradients are limited to preserve monotonicity When integrated over a grid cell the reconstructed fields must have mean values equal to the known state variables denoted by a for fractional area h for thickness and for enthalpy The mean values are not in general equal to the values at the cell center For example the mean ice area must equal the value at the centroid which may not lie at the cell center Consider first the fractional ice area the analog to fluid density p in 7 For each thickness category we construct a field a r whose mean is a where r x y is the position vector relative to the cell center That is we require n dad 11 A where A f4 dA is the grid cell area Equation 11 is satisfied if a r has the form a r aa Va r T 12 12 Model components where Va is a centered estimate of the area gradient within the cell a is a limiting coefficient that enforces monotonicity and r is the cell centroid an teh 5 ira It follows from 12 that the ice area at the cell center r 0 is Ac Q ArT
22. 4 umu 0 5 Pik cing 0 6 Peek 4 Y gt equa 0 7 In addition there are one or more equations describing tracer transport whose values are contained in the trern array These equations typically have one of the following three forms Gina y anTau 0 8 Win yon 0 9 Menta vont 0 10 Equation 8 describes the transport of surface temperature whereas 9 and 10 describe the transport of passive tracers such as volume weighted ice age and snow age Each tracer field is given an integer index txrcr_depend which has the value 0 1 or 2 depending on whether the appropriate conservation equation is 8 9 or 10 respectively The total number of tracers is Nir gt 1 In the default configuration there are two tracers surface temperature and volume weighted ice age Tracers for melt pond volume level ice area and level ice volume used to diagnose ridged ice area and volume are also available Users may add any number of additional tracers that are transported conservatively provided that t rcr_depend is defined appropriately See Section 4 8 3 for guidance on adding tracers The age of the ice is treated as an ice area tracer trcr_depend 0 It is initialized at O when ice forms as frazil and the ice ages the length of the timestep during each timestep Freezing directly onto the bottom of the ice does not affect the age nor does melting Mechanical redistribution processes and advection alter the age of
23. 5 1 Many ice models compute the sea surface slope VH from geostrophic ocean currents provided by an ocean model or other data source In our case the sea surface height Hs 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 to determine the stress between the ocean and the ice and subsequently the ice velocity w This stress relative to the ice Tij cupu Uno l Gu a cos0 k x Cw sin is then passed to the flux coupler relative to the ocean for use by the ocean model Here 6 is the turning angle between geostrophic and surface currents Cw is the ocean drag coefficient pw is the density of seawa ter dragw Cwpw and k is the vertical unit vector The turning angle is necessary 1f 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 6 0 is a good approximation Here we assume that the top layer is thin enough 3 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 thermodynamic and dynamic properties of the ice pack depend on how much 8 Model components
24. J Geophys Res Oceans 87 7971 7984 1982 29 G A Maykut and M G McPhee Solar heating of the Arctic mixed layer J Geophys Res Oceans 100 24691 24703 1995 30 G A Maykut and D K Perovich The role of shortwave radiation in the summer decay of a sea ice cover J Geophys Res 92 C7 7032 7044 1987 31 G A Maykut and N Untersteiner Some results from a time dependent thermodynamic model of sea ice J Geophys Res 76 1550 1575 1971 32 R J Murray Explicit generation of orthogonal grids for ocean models J Comput Phys 126 251 273 1996 33 N Ono Specific heat and heat of fusion of sea ice In H Oura editor Physics of Snow and Ice volume 1 pages 599 610 Institute of Low Temperature Science Hokkaido Japan 1967 34 D J Pringle H Eicken H J Trodahl and L G E Backstrom Thermal conductivity of landfast Antarctic and Arctic sea ice J Geophys Res 112 C04017 doi 10 1029 2006JC003641 2007 16 REFERENCES 35 D A Rothrock The energetics of the plastic deformation of pack ice by ridging J Geophys Res 80 4514 4519 1975 36 W Schwarzacher Pack ice studies in the Arctic Ocean J Geophys Res 64 2357 2367 1959 37 A J Semtner A model for the thermodynamic growth of sea ice in numerical investigations of climate J Phys Oceanogr 6 379 389 1976 38 R D Smith J K Dukowicz and R C Malone Parallel ocean general circulation modeling Phys
25. Table 6 If NTASK 1 then the serial code is used otherwise the code in mpi is used Note that the value of NTASK in comp_ice must equal the value of nprocs in ice in Generally the value of MKBLCKS computed by comp ice is sufficient but sometimes it will need to be set explicitly as discussed in Section 4 7 The scripts define a number of environment variables mostly as directories that you will need to edit for your own environment SSYSTEM USERDIR which on machines at Oak Ridge National Laboratory points automatically to scratch space is intended to be a disk where the run directory resides The reproducible option DITTO makes diagnostics bit for bit when varying the number of pro cessors The simulation results are bit for bit regardless because they do not require global sums or max mins as do the diagnostics This was done mainly by increasing the precision for the global reduc tion calculations except for regular double precision r8 calculations involving MPI MPI can not handle MPI_REAL16 on some architectures Instead these cases perform sums or max min calculations across the global block structure so that the results are bit for bit as long as the block distribution is the same the number of processors can be different CICE namelist variables available for changes after compile time appear in ice log with values read from the file ice_in their definitions are given in Section 5 6 For example to run for a differ
26. Tap lt 0 C 2 The change in T f since the previous iteration is less than a prescribed limit AT max 3 Fo gt Fe Uf Fo lt Fet ice would be growing at the top surface which is not allowed 4 The rate at which energy is added to the ice by the external fluxes equals the rate at which the internal ice energy is changing to within a prescribed limit A Fmax We also check the convergence rate of Tsp If Tsp is oscillating and failing to converge we average temper atures from successive iterations to improve convergence When all these conditions are satisfied usually within two to four iterations for AT max gt 0 01 C and AF nar 0 01 W m the calculation is complete 3 5 3 Growth and melting We first derive expressions for the enthalpy q The enthalpy of snow or fresh ice is given by qs T ps eoT Lo Sea ice enthalpy is more complex because of brine pockets whose salinity varies inversely with temperature The specific heat of sea ice given by 61 includes not only the energy needed to warm or cool ice but also the energy used to freeze or melt ice adjacent to brine pockets Equation 61 can be integrated to give the energy e required to raise the temperature of a unit mass of sea ice of salinity S from T to T 1 1 BTT eo T T Lows 7 If we let T Tm uS the temperature at which the ice is completely melted we have Tm e T Tri co Tm T Lo 1 The
27. a system of equations for the new temperatures For Tr lt 0 C we require Fo Eu 66 where Fe is the conductive flux from the top surface to the ice interior and both fluxes are evaluated at time m 1 Although Fo is a nonlinear function of T s we can make the linear approximation dFy men A Gat sf where T s is the surface temperature from the most recent iteration and Fy and dFo dTf are functions of T of We initialize T sf T sf and update it with each iteration Thus we can rewrite 66 as dF i Fo K TRT E Sir dT Rearranging terms we obtain dFo i m 1 ammtl1 _ dFo j x Km el TE 67 the first equation in the set of equations 59 The temperature change in ice snow layer k is Gd At Ahk PkCk KT TET Ke TET AD 68 where To T y in the equation for layer 1 In tridiagonal matrix form 68 becomes kK eTR 1 Ke Keo Ty Kra Ted Te mele 69 where nk At prckAh In the equation for the bottom ice layer the temperature at the ice ocean interface is held fixed at T s the freezing temperature of the mixed layer thus the last term on the LHS is 32 Model components known and is moved to the RHS If T 0 C then there is no surface flux equation In this case the first equation in 59 is similar to 69 but with the first term on the LHS moved to the RHS These equations are modified if T p and F are computed within the atmospheric m
28. 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 ncat 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 eg CCSM 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 5 2 Slow execution 51 5 2 Slow execution On some architectures underflows 107300 for example are not flushed to zero automatically Usually a c
29. below sea level when Pw pi hi Ps h h gt 0 In this case we raise the snow base to sea level by converting some snow to ice bh zph Pw tras BA Pw In rare cases this process can increase the ice thickness substantially For this reason snow ice conversions are postponed until after the remapping in thickness space Section 3 2 which assumes that ice growth during a single time step is fairly small Lateral melting is accomplished by multiplying the state variables by 1 rs ge Where TP side is the fraction of ice melted laterally and adjusting the ice energy and fluxes as appropriate 4 Numerical implementation CICE is written in FORTRAN90 and runs on platforms using UNIX LINUX and other operating sys tems The code is parallelized via grid decomposition with MPI and includes some optimizations for vector architectures A second external layer of parallelization involves message passing between CICE and the flux cou pler which may be running on different machines in a distributed system The parallelization scheme for CICE was designed so that MPI could be used for the coupling along with either MPI or no parallelization internally The internal parallelization method is set at compile time with the NTASK definition in the com pile script Message passing between the ice model and the CCSM flux coupler is accomplished with MPI regardless of the type of internal parallelization used for CICE although
30. block The physical portion of a subdomain is indexed as ilo ihi jlo jhi with nghost ghost or halo cells outside the domain used for boundary conditions These parameters are illustrated in Figure 5 in one dimension The routines global_scatter and global_gather distribute infor mation from the global domain to the local domains and back respectively If MPI is not being used for grid decomposition in the ice model these routines simply adjust the indexing on the global domain to the single local domain index coordinates Although we recommend that the user choose the local domains so that the global domain is evenly divided if this is not possible then the furthest east and or north blocks will contain nonphysical points padding These points are excluded from the computation domain and have little effect on model performance The user chooses a block size BLCKXxBLCKY and the number of processors NTASK in comp_ice Parameters in the domain_nml namelist in ice_in determine how the blocks are distributed across the proces sors and how the processors are distributed across the grid domain Recommended combinations of these parameters for best performance are given in Section 4 7 The script comp_in computes the maximum num ber of blocks on each processor for typical Cartesian distributions but for non Cartesian cases MKBLCKS may need to be set in the script The code will print this information to the log file before aborting an
31. compute an oceanic heat flux Fy Fil lt Ffremi 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 Fhocn 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 Fwater includes ice meltwater 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 Fox 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 Fswy penetrates the snow and ice layers and passes into the ocean as described in Section 3
32. densities of sea ice and pure ice Whereas the parameterization in 62 asymptotes to a constant conductivity of 2 03 W m t K with decreasing T K in 63 continues to increase with colder temperatures The equation for temperature changes in snow is analogous to 60 with ps 330 kg m3 cs co and K 0 30 W m deg replacing the corresponding ice values If shortwave default then the penetrating solar radiation is equal to zero for snow covered ice since most of the incoming sunlight is absorbed near the top surface If shortwave dEdd however then Ipen is nonzero in snow layers We now convert 60 to finite difference form The resulting equations are second order accurate in space except possibly at material boundaries and first order accurate in time Before writing the equations in full we give finite difference expressions for some of the terms First consider the terms on the left hand side of 60 We write the time derivatives as OT T Tm OE At where T is the temperature of either ice or snow and m is a time index The specific heat of ice layer k is approximated as Lop siz ml Ek Tik which ensures that energy is conserved during a change in temperature This can be shown by using 61 to integrate c dT from T to TOA the result is Cik TF Ty where ci is given by 64 The specific heat is a nonlinear function of rhage the unknown new temperature We can retain a set of linear equations
33. equations written in curvilinear coordi nates In this manner metric terms associated with the curvature of the grid were incorporated into the discretization explicitly We no longer use a triangle discretization in which the strain rates and stresses were constant over each of four subtriangles within each grid cell and instead use a bilinear approximation for the velocities and stresses Details pertaining to the spatial discretization are found in 18 The momentum equation is discretized in time as follows First for clarity the two components of 49 are Ou 29 001 e i H mz DE Tax QiCwPw Uw ul Uw u cos 0 Vi v sin 0 mfu mg a Hee ey Tay iCwPw Uw ul U u sind Vw v cos 6 mfu mee at Oz Oy In the code vrel aicwPw Uw u where k denotes the subcycling step The following equations illustrate the time discretization and define some of the other variables used in the code da OH vrel cos 9 u t1_ mf vrelsin 6 vtt Ht ran mg vrel Uy cos 9 Vy sin 0 ti At AAA Ox Ox A AS AG ccb han amd waterx cca strintx forcex k 1 OH mf vrel sin 6 utta 2 vrel cos 9 yt 2 Tay Mg vrel Uw sin 6 Vy cos 0 qa IA At Ox Oy NS AE en ee ccb _ OS vatery cea strinty forcey 26 Model components and vrel waterx y taux y We solve this system of equations analytically for u and v When th
34. extra diagnostics false Table 7 Namelist options continued next page Namelist options 55 variable options format description recommended value Model Output continued histfreq string array defines output frequencies histfreq n hist_avg history dir history file history_format write_ic incond_dir incond_file runid grid_nml grid_format grid_type grid_file kmt_file kcatbound domain_nml nprocs processor_shape distribution type distribution weight R5T0O3R x integer array 0 true false path filename prefix ne bin true false path filename prefix string ne bin rectangular displaced_pole tripole panarctic filename filename 0 1 2 integer slenderX1 slenderX2 square ice square pop cartesian rake spacecurve block latitude write history every hist freq_n years write history every hist freq_n months write history every hist freq_n days write history every hist freq_n hours write history every time step unused frequency stream i e do not write frequency history output is written do not write to history write time averaged data write snapshots of data path to history output directory output file for history write netCDF history file write direct access binary file write initial condition path to initial condition directory output file for initial condition label for run currently CCSM only Grid read netC
35. ice in any given grid cell in a conservative manner following changes in ice area The sea ice age tracer is validated in 16 By default ice and snow are assumed to have constant densities so that volume conservation is equiv alent to mass conservation Variable density ice and snow layers can be transported conservatively by defining tracers corresponding to ice and snow density as explained in the introductory comments in ice_transport_remap F90 Prognostic equations for ice and or snow density may be included in future model versions but have not yet been implemented Two transport schemes are available upwind and the incremental remapping scheme of 7 as modified for sea ice by 25 The remapping scheme has several desirable features e It conserves the quantity being transported area volume or energy e It is non oscillatory that is it does not create spurious ripples in the transported fields e It preserves tracer monotonicity That is it does not create new extrema in the thickness and enthalpy fields the values at time m 1 are bounded by the values at time m Horizontal transport 11 e It is second order accurate in space and therefore is much less diffusive than first order schemes e g upwind The accuracy may be reduced locally to first order to preserve monotonicity e It is efficient for large numbers of categories or tracers Much of the work is geometrical and is performed only once per grid cell instead of be
36. ice is transferred Other conserved quantities are transferred in proportion to the ice volume Avin We now use the subscripts 2 and s to distinguish ice from snow For example the transferred snow volume is Avs Usn Avin vin and the transferred ice energy in layer k is Aeing Cinkl AVin Vin The left and right boundaries of the domain require special treatment If ice is growing in open water at a rate Fo the left boundary Ho is shifted to the right by Fo At before 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 Ho is fixed at zero and the growth rate at the left boundary is estimated as Fo f If Fo lt 0 all ice thinner than Ahy FoAt is assumed to have melted and the ice area in category 1 is reduced accordingly The area of new open water is Aho Aao gdh 0 The right boundary Hy is not fixed but varies with hy the mean ice thickness in the thickest category Given hy we set Hy 3hy 2Hy 1 which ensures that g h gt 0 for Hy_1 lt h lt Hy and g h 0 for h gt Hy 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 functi
37. melting potential 4 7 33 meltwater 5 7 27 33 mixed layer 26 33 37 momentum equation 24 monotonicity 10 13 18 MPI 34 36 38 40 41 44 45 50 namelist 36 54 58 National Center for Atmospheric Research see NCAR NCAR 2 3 56 57 Oak Ridge National Laboratory 45 ocean 7 currents see currents ocean density see density ocean heat 4 7 29 mixed layer see mixed layer salinity see salinity ocean stress see ice ocean stress surface height see height sea surface surface slope see slope sea surface temperature see temperature ocean open water 3 5 7 8 20 21 22 Parallel Ocean Program see POP parallelization 34 36 45 POP 2 7 38 40 41 57 radiation Delta Eddington 9 27 longwave 4 5 27 28 shortwave 4 5 7 9 27 28 29 51 rain 4 5 7 reference height see height reference humidity see humidity reference temperature see temperature reference regional 40 remapping incremental 9 9 20 37 41 linear see transport thickness replacement pressure 24 reproducible 45 restart 35 37 41 44 45 52 restoring 37 ridging 2 8 9 21 23 37 44 49 salinity ice 7 26 30 32 33 INDEX ocean 4 7 salt see salinity sensible heat 4 6 27 28 shortwave see radiation shortwave slope sea surface 4 5 7 24 snow 2 4 5 7 8 10 12 20 22 24 26 34 51 solar see radiation shortwave space filling curve 37 specific humidity see humidity specific s
38. model for sea ice dynamics 13 It reduces to the VP model at time scales associated with the wind forcing while at shorter time scales the adjustment process takes place 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 The EVP sea ice dynamics model is thoroughly documented in 17 15 18 and 19 Simulation results and performance of this model have been compared with the VP model in realistic simulations of the Arctic 20 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 18 and 19 The force balance per unit area in the ice pack is given by a two dimensional momentum equation 13 obtained by integrating the 3D equation through the thickness of the ice in the vertical direction 0 a mo V o Fat Fo k x mfu mgVHo 49 where m is the combined mass of ice and snow per unit area and 7 and 7 are wind and ocean stresses respectively The strength of the ice is represented by the internal stress tensor o 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
39. module has been generalized to allow output at different frequencies Five output frequen cies 1 h d m y are available simultaneously during a run The same variable can be output at different frequencies say daily and monthly via its namelist flag var which is now a character string corre sponding to hist freq or x for none Grid variable flags are still logicals since they are written to all files no matter what the frequency is If there are no namelist flags with a given hist freq value or if an element of hist freq_n is 0 then no file will be written at that frequency The output period can be discerned from the filenames Model output 43 For example in namelist histireg Tit Pals Tid mir ey histfregn 1 6 0 1 1 fhi 1 fhs p f_Tsfc d f_aice m fmeltb mh fiage x Here hi will be written to a file on every timestep hs will be written once every 6 hours aice once a month meltb once a month AND once every 6 hours and Ts fc and iage will not be written From an efficiency standpoint it is best to set unused frequencies in histfreq to x Having output at all 5 frequencies takes nearly 5 times as long as for a single frequency If you only want monthly output the most efficient setting is hist freq m x x x x The code counts the number of desired streams nst reams based on hist freq The history variable names must be unique f
40. of ice_ tracer F90 in case other users do not want to use that tracer Four optional tracers are available in the code ice age melt pond volume and level ice area and volume from which ridged ice quantities are derived Age is a volume weighted tracer while melt pond volume is an area weighted tracer In the absence of sources and sinks the total mass of a volume weighted tracer such as soot kg is conserved under transport in horizontal and thickness space the mass in a given grid cell will change whereas the soot concentration kg m is unchanged following the motion and in particular the concentration is unchanged when there is surface or basal melting The proper units for a volume weighted mass tracer in the tracer array are kg m To add a tracer follow these steps the ice age tracer can be used as a pattern 1 ice_domain size F90 increase max_ntrcr 2 ice_state F90 declare nt_ tracer 3 ice_ tracer F90 create initialization physics restart routines 4 ice fileunits F90 add new dump and restart file units 5 ice_init F90 e add new module and tr_ t racer to list of used modules and variables e add logical namelist variable tr tracer e initialize namelist variable e broadcast namelist variable e print namelist variable to diagnostic output file e increment number of tracers in use based on namelist input nt rcr e define tracer types t rcr_depend 0 for ice area tracers 1 for ice volume 2 for snow vol
41. of each thickness category r 0 1 specifies the fraction of available liquid water captured by the ponds and r2 0 01 Here Tp is a reference temperature equal to 2 C Pond depth is assumed to be a linear function of the pond fraction hp 0 8ap and is limited by the category ice thickness hp lt 0 9h Ponds are allowed only on ice at least 10 cm thick The namelist parameters Rice and R_pnd adjust the albedo of bare or ponded ice by the product of the namelist value and one standard deviation For example if R ice 0 1 the albedo increases by 0 10 Similarly setting R snw 0 1 decreases the snow grain radius by 0 10 thus increasing the albedo See 5 for details the melt pond and Delta Eddington parameterizations are further explained and validated in a manuscript currently in preparation 2 Shortwave radiation CCSM3 In the parameterization used in the previous version of the Community Climate System Model CCSM3 the albedo depends on the temperature and thickness of ice and snow and on the spectral distribution of the incoming solar radiation Albedo parameters have been chosen to fit observations from the SHEBA field ex periment For Tp lt 1 C and h gt ahmax the bare ice albedo is 0 78 for visible wavelengths lt 700 nm and 0 36 for near IR wavelengths gt 700 nm As h decreases from ahmax to zero the ice albedo decreases smoothly using an arctangent function to the ocean albedo 0 06 The ice albedo in b
42. of the three terms on the right set to zero in each stage We compute horizontal transport using the incremental remapping scheme of 7 as adapted for sea ice by 25 this scheme is discussed in Section 3 1 Ice is transported in thickness space using the remapping scheme of 24 as described in Section 3 2 The mechanical redistribution scheme based on 43 35 14 11 and 26 is outlined in Section 3 3 To solve the horizontal transport and ridging equations we need the ice velocity u and to compute transport in thickness space we must know the the ice growth rate f in each thickness category We use the elastic viscous plastic EVP ice dynamics scheme of 17 as modified by 6 15 18 and 19 to find the velocity as described in Section 3 4 Finally we use the thermodynamic model of 4 discussed in Section 3 5 to compute f The order in which these computations are performed in the code itself was chosen so that quantities sent to the coupler are consistent with each other and as up to date as possible In earlier versions of CICE fluxes and state variables were sent to the coupler during the timestep after initial thermodynamic surface fluxes had been computed but before the rest of the thermodynamics and dynamics calculations This was done to improve load balancing in a concurrent coupling scheme model components ran simultaneously rather than sequentially CCSM no longer requires this and so the coupling is now done at
43. one dimensional 20 cell global domain decomposed into four local subdomains Each local domain has one ghost cell on each side and the physical portion of the local domains are labeled ilo ihi The parameter nx block is the total number of cells in the local domain including ghost cells and the same numbering system is applied to each of the four subdomains on the northeast corner of the corresponding T cells and have velocity in the center of each The velocity components are aligned along grid lines The user has several choices of grid routines popgrid reads grid lengths and other parameters for a nonuniform grid including tripole grids and rectgrid creates a regular rectangular grid including that used for the column configuration The panarctic grid is a regional configuration 27 The input files global_gx3 grid and global_gx3 kmt contain the 3 POP grid and land mask global_gx1 grid and global_gx1 kmt contain the 1 grid and land mask These are binary unformatted direct access files produced on an SGI Big Endian If you are using an incompatible Little Endian architecture choose rectangular instead of displaced pole in ice in or follow procedures as for coyote OS SITE machine Linux LANL coyote There are netCDF versions of the gx3 grid files available 4 2 1 Grid domains and blocks In general the global gridded domain is nx_global xny_global while the subdomains used in the block distribution are nx blockxny
44. output file ice_in namelist input data from cice input_templates gx3 hist iceh timeID ne output history file kmt land mask file from cice input_templates gx3 restart restart directory iced_gx3_v4 0_kcatbound0 initial condition from cice input_templates gx3 ice restart file restart pointer from cice input templates gx3 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 32 or 39 The ice and snow area volume and energy are given at the center of the cell velocity is defined at the corners and the internal ice stress tensor takes four different values within a grid cell bilinear approximations are used for the stress tensor and the ice velocity across the cell as described in 18 This tends to avoid the grid decoupling problems associated with the B grid EVP is available on the C grid through the MITgcm code distribution http mitgem org cgi bin viewcvs cgi MITgcm pkg seaice Since ice thickness and thermodynamic variables such as temperature are given in the center of each cell the grid cells are referred to as T cells We also occasionally refer to U cells which are centered 38 Numerical implementation Global Physical Domain 123 45 6 7 8 9 101112131415 1617181920 imt_global Local Domain nghost 1 nx_block o ilo ihi y y Figure 5 Grid parameters for a sample
45. physical measurements reflected absorbed and transmitted shortwave radiation apparent optical properties are then computed for each snow and ice layer in a self consistent manner Absorptive effects of inclusions in the ice snow matrix such as dust and algae can also be included along with radiative treatment of melt ponds and other changes in physical properties for example granu larization associated with snow aging The present code includes a simple parameterization of surface melt ponds The Delta Eddington formulation is described in detail in 5 Since publication of this technical paper a number of improvements have been made to the Delta Eddington scheme including a surface scat tering layer and internal shortwave absorption for snow generalization for multiple snow layers and more than four layers of ice and updated IOP values The melt pond parameterization is designed to be used in conjunction with the Delta Eddington short wave scheme Melt pond volume is carried on each ice thickness category as an area tracer Defined simply as the product of pond area ap and depth hp the melt pond volume vp grows through addition of ice or snow melt water or rain water and shrinks when the ice surface temperature becomes cold At pond growth v op t r1 an dhs En Brun Pw p p w W max Tp Ts fc 0 Tp pond contraction vp t At vpexp where dh and dh represent ice and snow melt at the top surface
46. restoring is not needed if the boundaries are only an ice sink the ice simply flows out but it is needed for an ice source The namelist variable rest ore_ice turns this functionality on and off and currently uses the restoring timescale t restore also used for restoring ocean sea surface temperature in stand alone ice runs This implementation is only intended to provide the hooks for a more sophisticated treatment the rectangular grid option can be used to test this configuration Much of the infrastructure used in CICE including the boundary routines is adopted from POP The boundary routines perform boundary communications among processors when MPI is in use and among blocks whenever there is more than one block per processor 4 2 5 Masks A land mask hm Mp is specified in the cell centers with O representing land and 1 representing ocean cells A corresponding mask uvm Mu for velocity and other corner quantities is given by Muli j min Mh l L i j 1 9 i j 1 1 974 1 The logical masks tmask and umask which correspond to the real masks hm and uvm respectively are useful in conditional statements In addition to the land masks two other masks are implemented in evp_prep in order to reduce the dynamics component s work on a global grid At each time step the logical masks ice_tmask and ice_umask are determined from the current ice extent such that they have the value true wherever ice exists The
47. to be fresh and the midpoint salinity S in each ice layer is given by S _l 1 Ga ik max cos Tz E 56 where z k 1 2 Ni Smax 3 2 psu and a 0 407 and b 0 573 are determined from a least squares fit to the salinity profile observed in multiyear sea ice by 36 This profile varies from S 0 at the top surface z 0 to S Smax at the bottom surface z 1 and is similar to that used by 31 Equation 56 is fairly accurate for ice that has drained at the top surface due to summer melting It is not a good approximation for cold first year ice which has a more vertically uniform salinity because it has not yet drained However the effects of salinity on heat capacity are small for temperatures well below freezing so the salinity error does not lead to significant temperature errors Each ice layer has an enthalpy q x defined as the negative of the energy required to melt a unit volume of ice and raise its temperature to 0 C Because of internal melting and freezing in brine pockets the ice enthalpy depends on the brine pocket volume and is a function of temperature and salinity Since the salinity is prescribed there is a one to one relationship between temperature and enthalpy We can also define a snow enthalpy qs which depends on temperature alone Expressions for the enthalpy are derived in Section 3 5 3 Given surface forcing at the atmosphere ice and ice ocean interfaces along with the ice and snow thicknesse
48. to the northeast corner The east cell triangles and selecting conditions are identical except for a rotation through 90 degrees This scheme was originally designed for rectangular grids Grid cells in CICE actually lie on the surface of a sphere and must be projected onto a plane The projection used in CICE 4 0 maps each grid cell to a square with sides of unit length Departure triangles across a given cell edge are computed in a coordinate system whose origin lies at the midpoint of the edge and whose vertices are at 0 5 0 and 0 5 0 Inter section points are computed assuming Cartesian geometry with cell edges meeting at right angles Let CL and CR denote the left and right vertices which are joined by line CLR Similarly let DL and DR denote the departure points which are joined by line DLR Also let IL and IR denote the intersection points 0 ya and 0 yp respectively and let IC xe 0 denote the intersection of CLR and DLR It can be shown that Ya Yb and x are given by LCL YDM YDL TDMYDL TDLYDM Ya EDM CDL j ICR YDR YDM TDMYDR TDRYDM b TDR TDM A EDR LDL fe fpr YDL YDR YDL Horizontal transport 15 Triangle Triangle Selecting logical group label condition 1 NW Ya gt Oand y gt 0 and z lt 0 NW1 Ya lt Oand yy gt 0 and z lt 0 W Ya lt Oand y lt 0 and z lt 0 W2 Ya gt Oand y lt 0 and z lt 0 2 NE yo gt O
49. where C1 ache C2 Gelet agr he C3 achy ayhe Ca Ache C5 achy ayhe ce ayhy From 15 and 16 the thickness and enthalpy at the cell center are given by he h het hi de 4 Quit qy where hz hy qx and qy are the limited gradients of thickness and enthalpy The surface temperature is treated the same way as ice or snow thickness but it has no associated enthalpy Tracers obeying conserva tion equations of the form 9 and 10 are treated in analogy to ice and snow enthalpy respectively We preserve monotonicity by van Leer limiting If i 7 denotes the mean value of some field in grid cell i j we first compute centered gradients of in the x and y directions then check whether these gradients give values of p within cell i 7 that lie outside the range of in the cell and its eight neighbors Let Pmax and Omin be the maximum and minimum values of over the cell and its neighbors and let max and din be the maximum and minimum values of the reconstructed within the cell Since the reconstruction is linear max and min are located at cell corners If max gt Omar or Omin lt Ditton we multiply the unlimited gradient by a min amax Amin Where Omax Pmax 0 max En d Amin Omin P min nt o Otherwise the gradient need not be limited Earlier versions of CICE through 3 14 computed gradients in physical space In version 4 0 gradients are computed in a scaled sp
50. where q 1 16378 x 107 kg m q2 5897 8 K TH is the surface temperature in Kelvin and pa is the surface air density The net downward heat flux from the ice to the ocean is given by 29 Foot Pi Ci Erie Tw LE Ty 58 where pw is the density of seawater c is the specific heat of seawater cy 0 006 is a heat transfer coefficient ux Y Twl pw is the friction velocity and Tw is the sea surface temperature A minimum value of u is available we recommend us min 5 X 1074 m s but the optimal value may depend on the ocean forcing used and can be as low as 0 Foot is limited by the total amount of heat available from the ocean Ffrzm Additional heat Fyide is used to melt the ice laterally following 30 and 40 Frog and the fraction of ice melting laterally are scaled so that Fit Fside Ffremit in the case that Ffrzmit lt 0 melting 3 5 2 New temperatures An option for zero layer thermodynamics 37 is available in this version of CICE by setting the namelist pa rameter heat_capacity to false and changing the number of ice layers ni lyr in ice_domain_size F90 to 1 In the zero layer case the ice is fresh and the thermodynamic calculations are much simpler than in the default configuration which we describe here Given the temperatures a T and Tij at time m we solve a set of finite difference equations to obtain the new temperatures at time m 1 Each temperature is coupled to the temperatures of the layers im
51. 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 19 and 6 For convenience we formulate the stress tensor in terms of 0 011 029 09 011 022 and introduce the divergence Dp and the horizontal tension and shearing strain rates Dr and Ds respectively The internal stress tensor is determined from a regularized version of the VP constitutive law 1 Oor O1 P H Dp E ot X 7 D 50 1 Oo 02 E 0t Ry a Ov 1 0012 012 1 H Ds 2 E ot 2n pre 6 where Dp 11 22 53 Dr n 22 54 Ds 2 12 55 et ee 1 Oui 0g da 2 Ox Ox i P ZA al P 1 AZ 1 1 2 E 2 2 2 A Dj DF D3 and P is a function of the ice thickness and concentration described in Section 3 3 The dynamics compo nent employs a replacement pressure see 12 for example which serves to prevent residual ice motion due to spatial variations of P when the rates of strain are exactly zero Dynamics 25 Several changes have been made to the EVP model since the original release In the earlier version the viscosities were held fixed while the stress and momentum equations were subcycled with the smaller time step dte The reason for implementing the EVP model in this way was to reproduce the results of the original VP model as closely as possible When solved with time steps of several hours or more the VP
52. wrapping the code with start and stop commands for that timer Printing of the timer output is done simultaneously for all timers To add a timer first declare it timer_ tmr at the top of ice_timers F90 we recommend doing this in both the mpi and serial directories then add a call to get_ice_timer in the subroutine init_ice_timers In the module containing the code to be timed call ice_timer_start timer_ tmr atthe beginning of the section to be timed and a similar call to ice timer_stop at the end Be careful not to have one command outside of a loop and the other command inside Timers can be run for individual blocks if desired by including the block ID in the timer calls Adding things 49 4 8 2 History fields To add a variable to be printed in the history output search for example in ice_history F90 1 add a frequency flag for the new field 2 add the flag to the namelist here and also in ice in 3 add an index number 4 broadcast the flag 5 add a call to define hist field 6 add a call to accum_hist_field To add an output frequency for an existing variable see section 4 5 1 4 8 3 Tracers For backward compatibility of restart files each new tracer has its own restart file generated in its own mod ule ice tracer F90 which also contains as much of the additional tracer code as possible We recommend that the logical namelist variable tr_ tracer be used for all calls involving the new tracer outside
53. 28 Hr The simplest polynomial that can satisfy both equations is a line It is convenient to change coordinates writing g n 919 go where y h H and the coefficients go and g are to be determined Then 27 and 28 can be written as nh g GONR An 2 3 2 nt go Anh where nr Hr Hr and nn hn Hr These equations have the solution 6a 2 R N m 29 NR 3 12a R n m ze 30 TR 2 Since g is linear its maximum and minimum values lie at the boundaries y 0 and nr 6a 2 R 0 55 m go 31 NR 3 6a R g nR m ze 32 NR 3 Equation 31 implies that g 0 lt 0 when 7 gt 27R 3 i e when hn lies in the right third of the thickness range Hr Hr Similarly 32 implies that g nr lt 0 when nn lt nr 3 i e when hn is in the left third of the range Since negative values of g are unphysical a different solution is needed when Ap lies outside the central third of the thickness range If h is in the left third of the range we define a cutoff thickness Ho 3h 2H and set g O between Hc and Hp Equations 29 and 30 are then valid with nr redefined as Ho Hr And if h is in the right third of the range we define Ho 3h 2Hp and set g 0 between Hr and He In this case 29 and 30 apply with nr Hr He and ny hn He Figure 3 illustrates the linear reconstruction of g for the simple cases Hz 0 Hr 1 an 1 and hn 0 2 0 4 0 6 an
54. CICE the Los Alamos Sea Ice Model Documentation and Software User s Manual Version 4 1 LA CC 06 012 Elizabeth C Hunke and William H Lipscomb T 3 Fluid Dynamics Group Los Alamos National Laboratory Los Alamos NM 87545 May 5 2010 Contents 1 Introduction 2 Coupling with other climate model components QA Atmosphere ann daar ke a wee a Sas A A a A DD OCEAN 5 oor BoD ar dea a kee tte ERLE wha dee ot ae A Ad A ide 3 Model components 3 1 Horizontal transport sn la Be A A Wk as BOE es ar eas 3 1 1 Reconstructing area and tracer fields 2 2 0 2 2 eee ee ee 3 1 2 Locating departure triangles 2 2 eee 3 1 3 Integrating fields viudos eek oe a be A eg El a de de 3 1 4 Updating state variables eee 3 2 Transport in thickness space 2 1 eee 3 3 Mechanical redistribution aaa 3 4 Dynamics ste O ect BOs 3 3 Thermodynamics Li i ewe ak rn de ee an ta ae Ge Be Soe aha 3 5 1 Thermodynamic surface forcing eee 3 32 New temperatures ce ca mera ee ee er ae ere a ee 3 5 3 Growth and melting a 4 Numerical implementation 4 1 Directory structure nnn ne aar a a annen raadden re hae ER EE BAREL 4 2 Grid boundary conditions and masks 2 o o e e 4 2 1 Grid domains and blocks o o e eee eee 4 22 Tmpolegrids nn verser ob Gran ASA ler En te ERE er Gor 4 2 3 Column configuration adai a d da a e i a e E E E A a 4 2 4 Boundary conditions deaa eeen
55. DF grid and kmt files read direct access binary file defined in rectgrid read from file in popgrid read from file in popgrid read from file in panarctic_grid name of grid file to be read name of land mask file to be read original category boundary formula new formula with round numbers WMO standard categories Domain number of processors to use 1 processor in the y direction tall thin 2 processors in the y direction thin more processors in x than y square more processors in y than x square distribute blocks in 2D Cartesian array redistribute blocks among neighbors distribute blocks via space filling curves full block size sets work_per_block latitude ocean sets work per block true iceh nc iceh bin displaced_pole grid kmt 0 Table 7 Namelist options continued 56 Namelist options variable options format description recommended value Domain continued ew_boundary_type cyclic periodic boundary conditions in x direction open Neumann boundary conditions in x closed Dirichlet boundary conditions in x ns_boundary _type cyclic periodic boundary conditions in y direction open Neumann boundary conditions in y closed Dirichlet boundary conditions in y tripole U fold tripole boundary conditions in y tripoleT T fold tripole boundary conditions in y tracer_nml Tracers tr_iage true false ice age restart_age true false rest
56. David Bailey Cecilia Bitz Bruce Briegleb Tony Craig Marika Holland John Dennis Julie Schramm Bonnie Light and Phil Jones e Jonathan Gregory of the University of Reading and the U K MetOffice for supplying tripole T fold code and documentation e Alison McLaren Ann Keen and others working with the Hadley Centre GCM for testing non standard model configurations and providing their code to us e the many researchers who tested beta versions of CICE 4 0 and waited patiently for the official release Acknowledgments and Copyright 53 Copyright 2010 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 54 Table of namelist options Namelist options variable options format description recommend
57. ace in which each grid cell has sides of unit length The origin is at the cell center and the four vertices are located at 0 5 0 5 0 5 0 5 0 5 0 5 and 0 5 0 5 In this coordinate system several of the above grid cell mean quantities vanish because they are odd functions of x and or y but they have been retained in the code for generality 3 1 2 Locating departure triangles The method for locating departure triangles is discussed in detail by 7 The basic idea is illustrated in Figure 1 which shows a shaded quadrilateral departure region whose contents are transported to the target or home grid cell labeled H The neighboring grid cells are labeled by compass directions NW N NE W and E The four vectors point along the velocity field at the cell corners and the departure region is formed by joining the starting points of these vectors Instead of integrating over the entire departure region it is convenient to compute fluxes across cell edges We identify departure regions for the north and east edges of each cell which are also the south and west edges of neighboring cells Consider the north edge of the home cell across which there are fluxes from the neighboring NW and N cells The contributing region from the NW cell is a triangle with vertices abc and that from the N cell is a quadrilateral that can be divided into two triangles with vertices acd and ade Focusing on triangle abc we first determine the coordinates of
58. alues can be altered at run time are handled in input_data and other initialization routines These variables are given default values in the code which may then be changed when the input file ice_in is read Other physical constants numerical parameters and variables are first set in initialization routines for each ice model component or module Then if the ice model is being restarted from a previous run some variables are read and reinitialized in restartfile Finally albedo and other quantities dependent on the initial ice state are set Some of these parameters will be described in more detail in Table 7 Three namelist variables control model initialization ice_ic runtype and restart as described in Table 4 It is possible to do an initial run from a file filename in two ways 1 set runtype ini tial restart true and ice_ic filename or 2 runtype continue and pointer_file Jrestart ice restart file where restart ice restart file contains the line restart filename The first op tion is convenient when repeatedly starting from a given file when subsequent restart files have been written MPI is initialized in init_communicate for both coupled and stand alone MPI runs The ice component communicates with a flux coupler or other climate components via external routines that handle the variables listed in Table 1 For stand alone runs routines in ice_forcing F90 read and interpolate data from files and are
59. and ya gt 0 and z2 gt 0 NE1 Yo lt 0 and y2 gt O and z2 gt 0 E yp lt O and y2 lt 0 and z2 gt 0 E2 yp gt O and y2 lt 0 and z2 gt 0 3 WI Ya lt 0 and y gt 0 and z lt 0 NW2 Ya gt Qand y lt 0 and z lt 0 El Yb lt 0 and yo gt 0 and x2 gt 0 NE2 yp gt 0 and y2 lt 0 and x2 gt 0 4 Hla Yayo gt 0 and Ya yo lt 0 Nla YaYo Oand ya yb gt 0 Hib Yayr lt 0 and y lt 0 Nlb YaYo lt 0 and y gt 0 5 H2a Yayo gt O and Ya yo lt 0 N2a Yayo Z 0 and Ya Yo gt 0 H2b YaYo lt 0 and ya lt 0 N2b YaYp lt 0 and ya gt 0 Table 3 Evaluation of contributions from the 20 triangles across the north cell edge The coordinates 71 2 Y1 Y2 Ya and y are defined in the text We define y y if x gt 0 else J Ya Similarly yo Ya if x2 lt 0 else Yo yp 16 Model components N soot SON NW iy N SNE sd 1 Le w aw Home Ding E W WSs A Home a b NW NE Lz pe a x BEA F gt eck Ww vy Home Ga E w 7 Home ee E i d1 d2 Figure 2 The 20 possible triangles that can contribute fluxes across the north edge of a grid cell Each departure triangle is defined by three of the seven points CL CR DL DR IL IR IC Given a 2D velocity field u the divergence V u in a given grid cell can be computed from the local velocities and written in terms of fluxes across each cell edge 1 are PERES ze ES
60. ansfer coefficient for sensible heat strain rate II component southern latitude of artificial mask edge e flag for shortwave parameterization default or dEdd principal stress components diagnostic sine of the turning angle in water 00 snow ice formation length scale for parameterizing nonuniform snow coverage special value single precision special value double precision sea surface slope in the x y direction sea surface salinity e source of surface salinity data sea surface temperature e source of surface temperature data shortwave radiation absorbed in snow layers Stefan Boltzmann constant 0 scene ee if 1 end program execution stress on ice by air in the x y direction centered in U cell stress on ice by air x y direction centered in T cell wind stress components from data ice strength pressure internal ice stress 012 internal ice stress 011 092 internal ice stress 011 092 divergence of internal ice stress x y kg m 1000 0 kg m 917 kg m 330 kg m 1026 kg m 100 x 10 m 1000 x 10 m 500 x 10 m 250 x 10 m psu 3 2 psu 86400 l s 30 N 0 02 m 193 1030 m m psu C W m 5 67x 108 W m K N m N m N m N m N m N m N m N m Index of primary variables and parameters strocnx y strocnx y T strtltx y swv n dr f T Tair tarea tarean tarear tareas tday
61. arameters p25 p333 p4 p5 p52083 p5625m p6 p666 p75 pi pi2 pih pointer_file potT precip_units print global print_points processor_shape prs_sig Pstar puny Q Qa qdp qqqice qqgocn Qref R R_ice R_pnd R_snw rl6_kind rad_to_deg radius rdg_conv rdg_shear real_kind refindx restart restart_age restart_dir restart_file restart_pond restore_ice 1 4 1 3 2 5 1 2 25 48 9 16 3 5 2 3 3 4 T 277 m 2 e input file for restarting atmospheric potential temperature e liquid precipitation data units e if true print global data 0 e if true print point data 02 eee e descriptor for processor aspect ratio replacement pressure ice strength parameter ne a small positive number ooo specific humidity at 10 m deep ocean heat flux for saturated specific humidity over ice for saturated specific humidity over ocean 2m atmospheric reference specific humidity e parameter for Delta Eddington ice albedo e parameter for Delta Eddington pond albedo e parameter for Delta Eddington snow albedo definition of quad precision degree radian conversion eee eee earth radius eeen convergence for ridging shear for ridging definition of single precision real refractive index of sea ice n e if true initialize using restart file instead of defaults e if true read age restart file
62. arameters effective for tuning sea ice thickness to the namelist e Implemented tracers for tracking the level and ridged ice area and volume e Generalized the history module to allow output at different frequencies simultaneously e Added multi dimensional variables including vertical and category dimensions to the history output e Reduced the work and memory required for multiple tracers e Added an alternative T fold tripole grid option e Added a column configuration e Corrected bugs particularly for nonstandard configurations Generally speaking subroutine names are given in italic and file names are boldface in this document Symbols used in the code are typewritten while corresponding symbols in this document are in the math font which granted looks a lot like italic A comprehensive index including a glossary of symbols with many of their values appears at the end The organization of this software distribution is described in Section 4 1 most files and subroutines referred to in this documentation are part of the CICE code found in cice source unless otherwise noted After many years CICE has finally become an acronym for Community Ice CodE We still pro nounce the name as sea ice but there has been a small grass roots movement underway to alter the model name s pronunciation from sea ice to what an Italian might say ch cha or chee chay Others choose to say sis English rhymes wi
63. art tracer values from file Ere vl true false level ice area and volume restart _lvl true false restart tracer values from file tr_pond true false melt ponds restart_pond true false restart tracer values from file ice_nml Physical Parameterizations kitd 0 delta function ITD approximation 1 linear remapping ITD approximation kdyn 0 dynamics OFF 1 EVP dynamics ndte integer number of EVP subcycles 120 kstrength 0 ice strength formulation 13 1 ice strength formulation 35 krdg_partic 0 old ridging participation function 1 new ridging participation function krdg_redist 0 old ridging redistribution function 1 new ridging redistribution function mu_rdg real e folding scale of ridged ice advection remap linear remapping advection remap upwind donor cell advection heat_capacity true salinity dependent thermodynamics true false zero layer thermodynamic model conduct MU71 conductivity 31 bubbly conductivity 34 shortwave default NCAR CCSM3 distribution method default dEdd Delta Eddington method albedo type default NCAR CCSM3 albedos default constant four constant albedos albicev 0 lt a lt l visible ice albedo for thicker ice albicei 0 lt a lt l near infrared ice albedo for thicker ice albsnowv 0O lt a lt l visible cold snow albedo albsnowi 0O lt a lt l near infrared cold snow albedo ahmax real albedo is constant above this thickness 0 3 m Table 7 Namelist options continued Namelist options 57 variable option
64. articipating in ridging is a as compared to G 3 for 35 For typical ice thickness distributions setting a 0 05 with krdg_partic 1 gives participation fractions similar to those given by G 0 15 with krdg_partic 0 See 26 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 vol ume and internal energy There are two namelist options for redistributing ice among thickness categories If krdg_redist 0 ridging ice of thickness hn forms ridges whose area is distributed uniformly between Hmin 2hn and Hmax 2V H hn as in 14 The default value of H is 25 m as in earlier versions of CICE Observations suggest that H 50 m gives a better fit to first year ridges 1 although the lower 22 Model components value may be appropriate for multiyear ridges 11 The ratio of the mean ridge thickness to the thickness of ridging ice is kn Hmin Hmax 2hn If the area of category n is reduced by ridging at the rate rp the area of thicker categories grows simultaneously at the rate r k Thus the net rate of area loss due to ridging of ice in category n is rn 1 1 k The ridged ice area and volume are apportioned among categories in the thickness range Hin Hmax The fraction of the new ridge area in category m is et nO hh 39 Amax Amin where Hy max Hm 1 Hmin and Hr min Hm Hmax The fraction of
65. atent heat flux is negative Le latent heat is transferred from the ice to the atmosphere snow or snow free ice sublimates at the top surface If the latent heat flux is positive vapor from the atmosphere is deposited at the surface as snow or ice The thickness change of the surface layer is given by Fo Tin Tf pLy q h FA 73 where p is the density of the surface material snow or ice and L 2 501 x 10 J kg is the latent heat of vaporization of liquid water at 0 C Note that pL is nearly an order of magnitude larger than typical values of q For positive latent heat fluxes the deposited snow or ice is assumed to have the same enthalpy as the existing surface layer After growth and melting the various ice layers no longer have equal thicknesses We therefore adjust the layer interfaces conserving energy so as to restore layers of equal thickness Ah h N This is done by computing the overlap nj of each new layer k with each old layer m Nkm Min Zm Zk max 2m 1 2k 1 34 Numerical implementation where Zm and zj are the vertical coordinates of the old and new layers respectively The enthalpies of the new layers are i dk Ah pe MkmUm m 1 At the end of the time step we check whether the snow is deep enough to lie partially below the surface of the ocean freeboard From Archimedes principle the base of the snow is at sea level when Pihi T pshs Pwhi Thus the snow base lies
66. ating system mpi modules that require MPI calls ice_boundary F90 boundary conditions ice_broadcast F90 routines for broadcasting data across processors ice_communicate F90 routines for communicating between processors ice_exit F90 aborts or exits the run ice gather scatter F90 gathers scatters data to from one processor from to all processors ice_global_reductions F90 global sums minvals maxvals etc across processors ice_timers F90 timing routines serial same modules as in mpi but without MPI calls source general CICE source code ice_age F90 handles most work associated with the age tracer ice atmo F90 stability based parameterization for calculation of turbulent ice atmosphere fluxes ice blocks F90 for decomposing global domain into blocks ice calendar F90 keeps track of what time it is ice_diagnostics F90 miscellaneous diagnostic and debugging routines ice_distribution F90 for distributing blocks across processors ice_domain F90 decompositions distributions and related parallel processing info ice_ domain size F90 domain and block sizes ice dyn _evp F90 elastic viscous plastic dynamics component ice fileunits F90 unit numbers for I O ice flux F90 fluxes needed produced by the model ice forcing F90 routines to read and interpolate forcing data for stand alone ice model runs ice grid F90 grid and land masks ice history F90 netCDF or binary output routines ice init F90 namelist and initializations ice itd F90 utilities f
67. between the sea ice model and the CCSM flux coupler are listed in Table 1 By convention directional fluxes are positive downward The ice fraction a aice is the total fractional ice coverage of a grid cell That is in each cell a 0 if there is no ice a 1 if there is no open water 0 lt a lt 1 if there is both ice and open water Typewritten equivalents used in the code are described in the index Coupling with other climate model components Atmosphere Ocean Provided by the flux coupler to the sea ice model Zo Atmosphere level height Pireomi Freezing melting potential Ua Wind velocity Tw Sea surface temperature Qa Specific humidity S Sea surface salinity Pa Air density VH Sea surface slope Oa Air potential temperature Ui Surface ocean currents Ta Air temperature Fy Shortwave radiation 4 bands Frl Incoming longwave radiation Frain Rainfall rate Fnow Snowfall rate Provided by the sea ice model to the flux coupler Ta Wind stress Fwy Penetrating shortwave F Sensible heat flux Fuater Fresh water flux F Latent heat flux Fhocn Net heat flux to ocean Fry Outgoing longwave Pan Salt flux Fevap Evaporated water Tij Ice ocean stress a Surface albedo 4 bands Tote Surface temperature a Ice fraction Tref 2m reference temperature diagnostic QT 2m reference humidity diagnostic Fswabs Absorbed shortwave diagnostic Table 1 Data exchanged between the CCSM flux coupler and the sea ice model 2 1 Atmo
68. binary files sss data type default constant values defined in the code clim climatological data ncar 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 Table 7 Namelist options continued 58 variable Namelist options options format description recommended value icefields_nml f_ var History Fields string frequency units for writing var to history write history every hist freq_n years write history every hist freq_n months write history every hist freq_n days write history every hist freq_n hours write history every time step do not write var to history md e g write both monthly and daily files Xx 503 Table 7 Namelist options continued Index of primary variables and parameters Index of primary variables and parameters 59 This index defines many of the symbols used frequently in the ice model code Values appearing in this list are fixed or recommended most namelist parameters are indicated e with their default values For other namelist options see Table 7 All quantities in the code are expressed in MKS units temperatures may tak
69. category Ho is set to zero The other boundaries are chosen with greater resolution for small h since the properties of the ice pack are especially sensitive to the amount of thin ice 28 The continuous function g h is replaced by the discrete variable ain defined as the fractional area covered by ice in the thickness range H _1 Hn We denote the fractional area of open water by ajo giving TAS Gin 1 by definition Category boundaries are computed in init_itd using one of several formulas summarized in Table 2 Setting the namelist variable kcat bound equal to 0 or 1 gives lower thickness boundaries for any number of thickness categories Nc Table 2 shows the boundary values for Nc 5 A third option specifies the boundaries based on the World Meteorological Organization classification the full WMO thickness distribution is used if Nc 7 if Nc 5 or 6 some of the thinner categories are combined The original formula kcatbound 0 is the default because it was used to create the restart files included with the code distribution Users may substitute their own preferred boundaries in init_itd In addition to the fractional ice area ain we define the following state variables for each category n e Vin the ice volume equal to the product of ain and the ice thickness hin e Usn the snow volume equal to the product of ain and the snow thickness Asn e cink the internal ice energy in layer k equal to the product of the ice layer vo
70. cen index of the bottom layer in each cat for eicen e directory to write snapshot of initial condition e prefix for initial condition file name definition of an integer veen eee eee selected real _kind 6 polynomial order of quadrature integrals in remapping 3 local processor coordinates on which to write debugging data local step counter for time loop e number of steps taken in previous run 0 total number of steps at current time step shortwave radiation absorbed in ice layers W m visible extinction coefficient in ice wavelength gt 700nm 17 6 m7 visible extinction coefficient in ice wavelength lt 700nm 14m e category boundary formula e type of dynamics 1 EVP 0 off 1 kg to g conversion factor noe 1000 thermal conductivity of fresh ice 0 eee eee eee 2 03 W m deg minimum conductivity of saline ice ooooocoommmo 0 10 W m deg e type of itd conversions 0 delta function 1 linear remap 1 e input file for land mask info e ridging participation function 0 eee eee eee 1 e ridging redistribution function 00 e eee eee 1 thermal conductivity of ice for zero layer thermodynamics 2 0 W m deg thermal conductivity of snow 0 30 W m deg e ice stength formulation 1 Rothrock 1975 0 Hibler 1979 1 flag for brine pocket effects if true check conservation when ridging flag for prescribing remapping fluxes e latitude of diagnostic poin
71. clean the compile directory and start fresh alter directories in the script clean_ice and execute the script In the run directory 1 Alter atm_data_dir and ocn_data_dir in the namelist file ice_in Execution procedures 45 variable options description RES gx3 gxl grid resolution BINTYPE MPI use MPI for internal parallelization NTASK integer total number of processors BLCKX integer number of grid cells on each block in the x direction BLCKY integer number of grid cells on each block in the y direction MXBLCKS integer maximum number of blocks per processor USE ESMF yes no for coupling with the Earth System Modeling Framework CAMICE yes no for single column CAM runs NETCDF yes no use no if netCDF library is unavailable DITTO yes no for reproducible diagnostics T Does not include ghost cells Table 6 Configuration options available in comp_ice 2 Alter the script run_ice for your system 3 Execute run_ice If this fails see Section 5 1 This procedure creates the output log file ice log ID and if npt is long enough compared with dumpfreg and hist freq dump files iced timeID and netCDF or binary history output files iceh_ timeID nc da Using the 3 grid the log file should be similar to ice log OS provided for the user s convenience These log files were created using MPI on 4 processors on the 3 grid Several options are available in comp ice for configuring the run shown in
72. d 00 eee e eee eee unknown location of field 00 field centered on west face cee eee eee angle field type iii eee ignore field type enn scalar field type eee unknown field type vector field type eeen latent heat flux effective floe diameter for lateral melt incoming longwave radiation outgoing longwave radiation Coriolis parameter mass in U cell rainfall rate frazil ice growth fresh water flux to ocean grid box mean fresh water flux fresh day of year that freezing begins freezing melting potential maximum magnitude of freezing melting potential net salt flux to ocean grid box mean salt flux to ocean fsalt sensible heat flux snowfall rate snow fraction that survives in ridging net surface heat flux excluding fcondtop incoming shortwave radiation absorbed shortwave radiation scaling factor to adjust ice quantities for updated data shortwave penetrating to ocean grid box mean shortwave penetrating to ocean fswthru current data year last data year e initial data year W m 1 s 1 x 1073 W m W m W m 1 W m 300 m W m W m kg s kg m s m kg m s kg m s W m 1000 W m kg m s kg m s W m kg m s 0 5 W m W m W m W m W m Index of primary variables and parameters G gravit grid_file grid format grid_type Gstar H halo_info heat_capacity hfrazi
73. d 0 8 Note that g slopes downward g lt 0 when hp is less than the midpoint thickness Hz Hr 2 1 2 and upward when h exceeds the midpoint thickness For h 0 2 and 0 8 g 0 over part of the range Finally we remap the thickness distribution to the original boundaries by transferring area and volume between categories We compute the ice area Aan and volume Av between each original boundary H and displaced boundary Af If H gt Hn ice moves from category n to n 1 The area and volume transferred are Hy Kees os wd 33 Ht Ne J hadik 34 An 20 Model components 0 0 2 0 4 0 6 0 8 1 Ice thickness Figure 3 Linear approximation of the thickness distribution function g h for an ice category with left boundary Hy 0 right boundary Hp 1 fractional area a 1 and mean ice thickness hn 0 2 0 4 0 6 and 0 8 If Af lt Hy ice area and volume are transferred from category n 1 to n using 33 and 34 with the limits of integration reversed To evaluate the integrals we change coordinates from h to y h Hz where Hr is the left limit of the range over which g gt 0 and write g 7 using 29 and 30 In this way we obtain the new areas an and volumes v between the original boundaries H _ and Hn in each category The new thicknesses hn Up dn are guaranteed to lie in the range H _1 Hn If g 0 in the part of a category that is remapped to a neighboring category no
74. d print the minimum and maximum values for an individual real array or its global sum 5 4 Known bugs 1 Fluxes sent to the CCSM 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 CCSM 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 2 With the standard CCSM radiative scheme shortwave default 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 lt 1 3 The date of onset diagnostic variables melt onset and frz_onset are not included in the 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 4 The single processor system_clock time may give erratic results on some architectures 5 History files that contain time averaged data hist avg true in ice_in will be incorrect if restarting from midway through an averaging period 52 Acknowledgments and Copyright 6 In stand alone runs
75. d the Grid boundary conditions and masks 39 user will need to adjust MXBLCKS in comp_ice and recompile The code will also print a warning if the maximum number of blocks is too large Although this is not fatal it does require excess memory A loop at the end of routine create_blocks in module ice_blocks F90 will print the locations for all of the blocks on the global grid if dbug is set to be true Likewise a similar loop at the end of routine create_local_block_ids in module ice_distribution F90 will print the processor and local block number for each block With this information the grid decomposition into processors and blocks can be ascertained The dbug flag must be manually set in the code in each case independently of the dbug flag in ice_in as there may be hundreds or thousands of blocks to print and this information should be needed only rarely This information is much easier to look at using a debugger such as Totalview 4 2 2 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 a physical boundary or grid singularity in the Arctic Ocean Instead of a single north pole it has two poles in the north both located on land with a line of grid points between them This line of points is called the fold and it is the top row of the physical grid One pole is at the left hand end of the top row and the ot
76. derX1 or slenderX2 often better for sea ice simulations on global grids where nearly all of the work is at the top and bottom of the grid with little to do in between and close to square domains which maximize the volume to surface ratio and therefore on processor compu tations to message passing if there were ice in every grid cell In cases where the number of processors is not a perfect square 4 9 16 the processor_shape namelist variable allows the user to choose how the processors are arranged Here again it is better in the sea ice model to have more processors in x than in y for example 8 processors arranged 4x2 square ice rather than 2x4 square pop The latter option is offered for direct communication compatibility with POP in which this is the default The distribution_type options allow standard Cartesian distribution of blocks redistribution via a rake algorithm for improved load balancing across processors and redistribution based on space filling curves The rake and space filling curve algorithms are primarily helpful when using squarish processor domains where some processors located near the equator would otherwise have little work to do Processor domains need not be rectangular however distribution_wght chooses how the work per block estimates are weighted The block option is the default in POP which uses a lot of array syntax requiring calculations over entire blocks whether or not land is present
77. e either Celsius or Kelvin units A a_min advection ahmax aice_extmin aice_init aice0 aice n albedo type albicei albicev albocn albsnowi albsnowv alpha alv n dr f alv n dr f _gbm ANGLE ANGLET apondn astar atm_data_dir atm_data_format atm_data_type atmbndy awtidf awtidr awtvdf awtvdr avgsiz minimum area concentration for computing velocity e type of advection algorithm used remap or upwind e thickness above which ice albedo is constant minimum value for ice extent diagnostic 004 concentration of ice at beginning of timestep fractional open water area total concentration of ice in grid cell in category n e type of albedo parameterization default or constant e near infrared ice albedo for thicker ice e visible ice albedo for thicker ice Ocean albedo roppa a hed ceed venne e near infrared cold snow albedo e visible cold snow albedo floe shape constant for lateral melt albedo visible near IR direct diffuse grid box mean value of alv n dr f for conversions between the POP grid and latitude longitude grids ANGLE converted to T cells area concentration of melt ponds e folding scale for participation function e directory for atmospheric forcing data e format of atmospheric forcing files e type of atmospheric forcing e atmo boundary layer parameterization default or constant weighti
78. e physical grid and points with T index i and nx_global i 2 are coincident Let the fold have T row index n on the global grid It is usual for the northernmost row of the physical domain to be a U row but in the case of the T fold the U row of index n is beyond the fold although it is not a ghost row it is not physically independent because it coincides with U row n 1 and it therefore has to be treated like a ghost row Points i on U row n coincides with nx_global i 1 on U row n 1 There are still ghost T and U rows n 1 to the north of U row n Ghost T row n 1 coincides with T row n 1 and ghost U row n 1 coincides with U row n 2 The tripole grid thus requires two special kinds of treatment for certain rows arranged by the halo update routines First within rows along the fold coincident points must always have the same value This is achieved by averaging them in pairs Second values for ghost rows and the quasi ghost U row on the T fold grid are reflected copies of the coincident physical rows Both operations involve the tripole buffer which is used to assemble the data for the affected rows Special treatment is also required in the scattering routine and when computing global sums one of each pair of coincident points has to be excluded 40 Numerical implementation 4 2 3 Column configuration A column modeling capability is available Because of the boundary conditions and other spatial a
79. e flux coupler The snowfall precipitation rate provided as liquid water equivalent and converted by the ice model to snow depth also contributes to the heat and water mass budgets of the ice layer Melt ponds generally form on the ice surface in the Arctic and refreeze later in the fall reducing the total amount of fresh water that reaches the ocean and altering the heat budget of the ice this version includes a simple melt pond parameterization Rain and all melted snow end up in the ocean Wind stress and transfer coefficients for the turbulent heat fluxes are computed in subroutine atmo_boundary_layer following 22 For clarity the equations are reproduced here in the present notation The wind stress and turbulent heat flux calculation accounts for both stable and unstable atmosphere ice boundary layers Define the stability ap 9 o T Q u 10 14 0 606Q 1 0 606 Q where is the von Karman constant g is gravitational acceleration and u O and Q are turbulent scales for velocity temperature and humidity respectively CS ve a 6 Coupling with other climate model components e co Oa Tafe 1 Q Cq Qa E The wind speed has a minimum value of 1 m s We have ignored ice motion in u and Tsfe and Qsf are the surface temperature and specific humidity respectively The latter is calculated by assuming a saturated surface as described in Section 3 5 1 The exchange coefficie
80. e path to restart dump files e restart file prefix e if true read melt pond restart file e if true restore ice state along lateral boundaries 67 F F N m 2 75x 104 N m 1x1071 kg kg W m 1 16378x 10 kg m 6 275724 x 10 kg m kg kg selected_real_kind 26 180 7 6 37x 10 m l s l s selected_real_kind 6 1 310 T 68 restore_sst rhoa rhofresh rhoi rhos rhow rnilyr rside rsnw fresh rsnw_melt rsnw_nonmelt rsnw_sig runtype S salin saltmax scale_factor sec secday shcoef shear shlat shortwave sigl 2 sinw snoice snowpatch spval spval_dbl ss_tltx y SSS sss_data_type sst sst_data_type Sswabs stefan boltzmann stop_now strairx y strairx y T strax y strength stress 12 stressm stressp strintx y Index of primary variables and parameters e restore sst to data air density density of fresh water enne density OF Ei density of SNOW oenen rr density of seawater eeen eene real nlyr fraction of ice that melts laterally freshly fallen snow grain radius melting snow grain radius 0 e eee eee eee nonmelting snow grain radius 000 000 standard deviation of snow grain radius e type of initialization used ice salinity max Salinity at ice base nnee scaling factor for shortwave radiation components seconds elasped into idate number of seconds in a day eee eee eee tr
81. e subcycling is finished for each thermodynamic time step the ice ocean stress must be con structed from taux y and the terms containing vrel on the left hand side of the equations This is done in subroutine evp finish 3 5 Thermodynamics The thermodynamic sea ice model is based on 31 and 4 and is described more fully in 23 For each thickness category the model computes changes in the ice and snow thickness and vertical temperature profile resulting from radiative turbulent and conductive heat fluxes The ice has a temperature dependent specific heat to simulate the effect of brine pocket melting and freezing Each thickness category n in each grid cell is treated as a horizontally uniform column with ice thickness hin Vin Qin and snow thickness hs Usn Qin Henceforth we omit the category index n Each column is divided into N ice layers of thickness Ah h N and N snow layers of thickness Ah hs Ns The surface temperature i e the temperature of ice or snow at the interface with the atmosphere is Tsp which cannot exceed 0 C The temperature at the midpoint of the snow layer is 75 and the midpoint ice layer temperatures are T where k ranges from 1 to N The temperature at the bottom of the ice is held at Tr the freezing temperature of the ocean mixed layer All temperatures are in degrees Celsius unless stated otherwise The vertical salinity profile is prescribed and is unchanging in time The snow is assumed
82. e use a thermodynamic model Section 3 5 to compute the new mean thicknesses h at time m 1 The time step must be small enough that trajectories do not cross i e AMPF lt peri for each pair of adjacent categories The growth rate at h hp is given by fn htt h At By linear interpolation we estimate the growth rate F at the upper category boundary Hp Hn hn pe fn En n f t Pati hn If an or Gn41 0 Fn is set to the growth rate in the nonzero category and if an an 1 0 we set Fn 0 The temporary displaced boundaries are given by H H F At n 1toN 1 The boundaries must not be displaced by more than one category to the left or right that is we require Hi lt H lt Hy 1 Without this requirement we would need to do a general remapping rather than an incremental remapping at the cost of added complexity Next we construct g h in the displaced thickness categories The ice areas in the displaced categories are a q since area is conserved following the motion in thickness space i e during vertical ice rowth or melting The new ice volumes are t anhn t a h For conciseness define n n n Transport in thickness space 19 Hzr H _ and Hp H and drop the time index m 1 We wish to construct a continuous function g h within each category such that the total area and volume at time m 1 are a and vn respectively Hr gdh an 27 Hr Hr haa ee
83. ed value setup_nml Time days per year 360or 365 number of days in a model year 365 vear 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 ndyn_dt integer number of dynamics advection ridging 1 steps per thermo timestep 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 restart_dir path path to restart directory 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 dumpfreg n days dumpfreq_n integer frequency restart data is written 1 Model Output diagfreg integer frequency of diagnostic output in dt 24 eg 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
84. edirect standard out as in run_ice otherwise it is written to the file given by diag_file In addition to the standard di agnostic output maximum area averaged thickness velocity average albedo total ice area and total ice and snow volumes the namelist options print_points and print_global cause additional diag nostic information to be computed and written print global outputs global sums that are useful for checking global conservation of mass and energy print_points writes data for two specific grid points Currently one point is near the North Pole and the other is in the Weddell Sea these may be changed in ice_diagnostics F90 Timers are declared and initialized in ice_timers F90 and the code to be timed is wrapped with calls to ice_timer_start and ice_timer_stop Finally ice_timer_print writes the results to the log file The optional stats argument true false prints additional statistics Calling ice_timer_print_all prints all of the timings at once rather than having to call each individually Currently the timers are set up as in Table 5 Section 4 8 1 contains instructions for adding timers The timings provided by these timers are not mutually exclusive For example the column timer 5 includes the timings from 6 7 8 and 9 and subroutine bound timer 14 is called from many different 44 Numerical implementation Timer Index Label 1 Total the entire run 2 Step total minus initialization and exit 3 Dynamic
85. ent length of time say three days set npt 72 in ice_in At present the user supplies the time step dt the number of dynamics advection ridging subcycles ndyn_dt and the number of evp subcycles ndte and dte is then calculated in subroutine init_evp The primary reason for doing it this way is to ensure that ndte is an integer To restart from a previous run set restart true inice in There are two ways of restarting from a given file The restart pointer file ice restart_file created by the previous run contains the name of the last written data file iced timeID Alternatively a filename can be assigned to ice_ic in ice_in Consult Section 4 3 for more details Restarts are exact for MPI or single processor runs 46 Numerical implementation Figure 6 Distribution of 256 blocks across 16 processors represented by colors on the gxl grid a cartesian slenderX1 b cartesian slenderX2 c cartesian square ice square pop is equivalent here d rake with block weighting e rake with latitude weighting f spacecurve Each block consists of 20x24 grid cells and white blocks consist entirely of land cells 4 7 Performance Namelist options domain nml provide considerable flexibility for finding the most efficient processor and block configuration Some of these choices are illustration in Figure 6 processor_shape chooses be tween tall thin processor domains s len
86. erged over categories but not scaled 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 cos p coszen and is zero if the sun is below the horizon cosy lt 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 cos y gt 0 The albedos will still be zero in the dark polar winter hemisphere Acknowledgments and Copyright This work has been supported through the Department of Energy Computer Hardware Applied Mathemat ics and Model Physics CHAMMP program Climate Change Prediction Program CCPP and Scientific Discovery through Advanced Computing SCIDAC program with additional support from the T 3 Fluid Dynamics Group at Los Alamos National Laboratory Special thanks are due to the following people e members of the CCSM Polar Climate Working Group including
87. eris tics among processors first in the x direction and then in y Space filling curves reduce a multi dimensional space 2D in our case to one dimension The curve is composed of a string of blocks that is snipped into sections again based on the work per processor and each piece is placed on a processor for optimal load balancing This option requires that the block size be chosen such that the number of blocks in the x direction equals the number of blocks in the y direction and that number must be factorable as 2 3 5P where n m p are integers For example a 16x16 array of blocks each containing 20x24 grid cells fills the gx1 grid n 4 m p 0 If either of these conditions is not met a Cartesian distribution is used instead The user provides the total number of processors and the block dimensions in the setup script comp_ice When moving toward smaller more numerous blocks there is a point where the code becomes less efficient blocks should not have fewer than about 20 grid cells in each direction Squarish blocks are probably better again to optimize the volume to surface ratio for communications In general the following rules of thumb seem to work best e For large numbers of processors processor_shape slenderX1 with one block per proces sor distribution type and distribution wght do not matter they default to cartesian slenderX2 is not quite as efficient for the same number of processors but is bette
88. h category are returned to the atmosphere model via the coupler Since the ice surface temperature is treated explicitly the effective conductivity may need to be limited to ensure stability As a result accuracy may be significantly reduced especially for thin ice or snow layers A more stable and accurate procedure would be to compute the temperature profiles for both the atmosphere and ice together with the surface fluxes in a single implicit calculation This was judged impractical however given that the atmosphere and sea ice models generally exist on different grids and or processor sets 2 1 Atmosphere The wind velocity specific humidity air density and potential temperature at the given level height z are used to compute transfer coefficients used in formulas for the surface wind stress and turbulent heat fluxes Ta Fs and Fj as described below Wind stress is arguably the primary forcing mechanism for the ice motion although the ice ocean stress Coriolis force and slope of the ocean surface are also important 41 The sensible and latent heat fluxes Fs and F along with shortwave and longwave radiation Fy FL and Fry are included in the flux balance that determines the ice or snow surface temperature when calc_Tsfc true As described in Section 3 5 these fluxes depend nonlinearly on the ice surface temperature T sf The balance equation is iterated until convergence and the resulting fluxes and Tsf are then passed to th
89. her is in the middle of the row The grid is constructed by folding the top row so that the left hand half and the right hand half of it coincide Two choices for constructing the tripole grid are available The one first introduced to CICE is called U fold which means that the poles and the grid cells between them are U cells on the grid Alternatively the poles and the cells between them can be grid T cells making a T fold Both of these options are also supported by the OPA NEMO ocean model which calls the U fold an f fold because it uses the Arakawa C grid in which U cells are on Trows The choice of tripole grid is given by the namelist variable ns boundary type tripole for the U fold and tripoleT for the T fold grid In the U fold tripole grid the poles have U index nx_globa1 2 and nx_global on the top U row of the physical grid and points with U index i and nx_global i are coincident Let the fold have U row index n on the global grid this will also be the T row index of the T row to the south of the fold There are ghost T and U rows to the north beyond the fold on the logical grid The point with index i along the ghost T row of index n 1 physically coincides with point nx global i 1 on the T row of index n The ghost U row of index n 1 physically coincides with the U row of index n 1 In the T fold tripole grid the poles have T index 1 and and nx_globa1 2 1 on the top T row of th
90. however by initially guessing Gavi Ti and then iterating the solution updating qua in 64 with each iteration until the solution converges Next consider the first term on the right hand side of 60 The first term describes heat diffusion and is discretized for a given ice or snow layer k as Cik Co 64 o OT 1 een q g m P ae x 7 NF T g a ES KE 1 TEAD 65 Thermodynamics 31 where Ah is the layer thickness and K is the effective conductivity at the upper boundary of layer k This discretization is centered and second order accurate in space except at the boundaries The flux terms on the right hand side RHS are treated implicitly i e they depend on the temperatures at the new time m 1 The resulting scheme is first order accurate in time and unconditionally stable The effective conductivity K at the interface of layers k 1 and k is defined as E Krab Kila which reduces to the appropriate values in the limits K gt Kp or vice versa and hy gt hpi or vice versa The effective conductivity at the top bottom interface of the ice snow column is given by K 2K Ah where K and Ah are the thermal conductivity and thickness of the top bottom layer The second term on the RHS of 60 is discretized as O Tk 1 Tk Ik De tien 2 I Ah Ah where 7 is the fraction of the penetrating solar radiation Jp that is transmitted through layer k without being absorbed We now construct
91. ica D 60 38 61 1992 39 R D Smith S Kortas and B Meltz Curvilinear coordinates for global ocean models Technical Report LA UR 95 1146 Los Alamos National Laboratory 1995 40 M Steele Sea ice melting and floe geometry in a simple ice ocean model J Geophys Res 97 C11 17729 17738 1992 41 M Steele J Zhang D Rothrock and H Stern The force balance of sea ice in a numerical model of the Arctic Ocean J Geophys Res Oceans 102 21061 21079 1997 42 A H Stroud Approximate Calculation of Multiple Integrals Prentice Hall Englewood Cliffs New Jersey 1971 431 pp 43 A S Thorndike D A Rothrock G A Maykut and R Colony The thickness distribution of sea ice J Geophys Res 80 4501 4513 1975 44 H J Trodahl S O F Wilkinson M J McGuiness and T G Haskell Thermal conductivity of sea ice dependence on temperature and depth Geophys Res Lett 28 1279 1282 2001 45 N Untersteiner Calculations of temperature regime and heat budget of sea ice in the Central Arctic J Geophys Res 69 4755 4766 1964
92. ient to yield a 1 Two tracers for tracking the ridged ice area and volume are available The actual tracers are for level undeformed ice area alv1 and volume v1v1 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 Alul T Ardg Qi Ulul T Urdg Vi Level ice area fraction and volume increase with new ice formation and decrease steadily via ridging pro cesses 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 The ice strength P may be computed in either of two ways If the namelist parameter kst rength 0 we use the strength formula from 13 P P hexp C 1 a 46 where P 27 500N m and C 20 are empirical constants and h is the mean ice thickness Alterna tively
93. ing repeated for each quantity being transported The time step is limited by the requirement that trajectories projected backward from grid cell corners are confined to the four surrounding cells this is what is meant by incremental remapping as opposed to general remapping This requirement leads to a CFL like condition max u At lt A Az T For highly divergent velocity fields the maximum time step must be reduced by a factor of two to ensure that trajectories do not cross However ice velocity fields in climate models usually have small divergences per time step relative to the grid size The remapping algorithm can be summarized as follows 1 Given mean values of the ice area and tracer fields in each grid cell construct linear approximations of these fields Limit the field gradients to preserve monotonicity 2 Given ice velocities at grid cell corners identify departure regions for the fluxes across each cell edge Divide these departure regions into triangles and compute the coordinates of the triangle vertices 3 Integrate the area and tracer fields over the departure triangles to obtain the area volume and energy transported across each cell edge 4 Given these transports update the state variables Since all scalar fields are transported by the same velocity field step 2 is done only once per time step The other three steps are repeated for each field in each thickness category These steps are described below
94. intended merely to provide guidance for the user to write his or her own routines Whether the code is to be run in stand alone or coupled mode is determined at compile time as described below 4 4 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 physical forcing CICE allows the dynamics advection and ridging portion of the code to be run with a shorter timestep At yn dt dyn than the thermodynamics timestep At dt In this case dt and the integer ndyn_dt are specified and dt_dyn dt ndyn_dt A conservative estimate of the horizontal transport time step bound or CFL condition under remapping yields min Az Ay Atdyn lt dyn 2 max u v Numerical estimates for this bound for several POP grids assuming max u v 0 5 m s are as follows grid label N pole singularity dimensions minyAzx Ay max Atayn gx3 Greenland 100x116 39x10 m 10 8hr gxl Greenland 320 x 384 18 x 10 m 5 0 hr p4 Canada 900 x 600 6 5x10 m 1 8hr 42 Numerical implementation As discussed in section 3 3 and 26 the maximum time step in practice is usually determined by the time scale for large changes in the ice strength which depends in part on wind strength Using the strength parameterization of 35 as in eq 47 limits the time step to 30 minutes for the old ridging scheme and to 2 hours for the new scheme assuming Az
95. io squared energy of melting of ice per unit area in category n emissivity of snow and ice eee eee eee eee a small number a small number a small number energy of melting of snow per unit area in category n evaporative water flux e if true use evp damping procedure 15 e type of east west boundary condition coefficient for calculating the parameter E 0 lt eyc lt 1 61 365 selected_real_kind 13 false l s 0 054 deg psu stdout l s l s 0 00536 kg m 3600 s 1 deg s 62 F fcondtop n _f fcor_blk ferrmax fhocn fhocn_gbm field_loc_center field_loc_Eface field_loc_NEcorner field_loc_Nface field_loc_noupdate field_loc_unknown field_loc_Wface field_type_angle field_type_noupdate field_type_scalar field type unknown field type vector flat floediam flw flwout fm frain frazil fresh fresh_gbm frz_onset frzmlt frzmlt max fsalt fsalt gbm fsens fsnow fsnowrdg fsurf n _f fsw fswabs fswfac fswthru fswthru_gbm fyear fyear_final fyear_init Index of primary variables and parameters conductive heat flux Coriolis parameter max allowed energy flux error thermodynamics net heat flux to ocean grid box mean net heat flux to ocean fhocn field centered on grid cell 0 c eee eee eee field centered on east face nn field on northeast corner ec eee eee eee eee field centered on north face ignore location of fiel
96. lmin hi_min hicen hin_max hist_avg histfreq histfreq_n history _dir history file history format hm hmix hour hp0 hpmin hpondn hs_min hsmin hs0 Hstar HTE HTN HTS HTW I 1 5 glob 10vis iblkp iG block ice_ic ice_ref_salinity icells iceruf gravitational acceleration oenen e input file for grid info e format of grid files e rectangular or displaced_pole or column piecewise linear ridging participation function parameter information for updating ghost cells e if true use salinity dependent thermodynamics minimum thickness of new frazil ice minimum ice thickness for thinnest ice category ice thickness in category n category thickness limits e if true write averaged data instead of snapshots e units of history output frequency y m w d or 1 e integer output frequency in histfreq units e path to history output files e history output file prefix e format of history files land boundary mask thickness T cell ocean mixed layer depth neen hour of the year pond depth at which transition to bare ice occurs minimum melt pond depth 0 cece eee eee eee melt pond depth minimum thickness for which Ts is computed minimum snow depth for Delta Eddington snow depth at which transition to ice occurs dEdd determines mean thickness of ridged ice
97. lume v N and the ice layer enthalpy ging Here N is the total number of ice layers with a default value N 4 and 3 1 Horizontal transport 9 dink 18 the negative of the energy needed to melt a unit volume of ice and raise its temperature to 0 C it is discussed in Section 3 5 NOTE In the current code e lt 0 and q lt 0 with e viqi e snk the internal snow energy in layer k equal to the product of the snow layer volume vs N and the snow layer enthalpy qsnk where Ns is the number of snow layers Similarly es lt 0 in the code Earlier versions of CICE had a single snow layer but multiple layers are now allowed The default value is N 1 e Tsm the surface temperature e Tage the volume weighted mean ice age Since the fractional area is unitless the volume variables have units of meters 1 e m of ice or snow per m of grid cell area and the energy variables have units of J m Both Ty and Tage are assigned to the tracer array trcrn which consists of N fields By default Nj 2 but users may create any number of additional tracer fields The three terms on the right hand side of 2 describe three kinds of sea ice transport 1 horizontal transport in x y space 2 transport in thickness space h due to thermodynamic growth and melting and 3 transport in thickness space h due to ridging and other mechanical processes We solve the equation by operator splitting in three stages with two
98. m T and Ty using the coefficients cu cg and cg Although the sea ice model does not use this quantity it is convenient for the ice model to perform this calculation The atmospheric reference temperature is returned to the flux coupler as a climate diagnostic The same is true for the reference humidity Qef Additional details about the latent and sensible heat fluxes and other quantities referred to here can be found in Section 3 5 1 2 2 Ocean 7 2 2 Ocean New sea ice forms when the ocean temperature drops below its freezing temperature Ty uS where S is the seawater salinity and u 0 054 psu is the ratio of the freezing temperature of brine to its salinity The ocean model performs this calculation if the freezing melting potential Ffrzm 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 In general this ice is added to the thinnest ice category The new ice is grown in the open water area of the grid cell to a specified minimum thickness if the open water area is nearly zero or if there is more new ice than will fit into the thinnest ice category then the new ice is spread over the entire cell If Ffrzm is negative it is used to heat already existing ice from below In particular the sea surface temperature and salinity are used to
99. mediately above and below by heat conduction terms that are treated implicitly For example the rate of change of Tij depends on the new temperatures in layers k 1 k and k 1 Thus we have a set of equations of the form Ax b 59 where A is a tridiagonal matrix x is a column vector whose components are the unknown new temperatures and b is another column vector Given A and b we can compute x with a standard tridiagonal solver There are four general cases 1 Tef lt 0 C snow present 2 Taf 0 C snow present 3 Ts lt 0 C snow absent and 4 Tsf 0 C snow absent For case 1 we have one equation the top row of the matrix for the new surface temperature N equations for the new snow temperatures and N equations for the new ice temperatures For cases 2 and 4 we omit the equation for the surface temperature which is held at 0 C and for cases 3 and 4 we omit the snow temperature equations Snow is considered absent if the snow depth is less than a user specified minimum value hs_min Very thin snow layers are still transported conservatively by the transport modules they are simply ignored by the thermodynamics The rate of temperature change in the ice interior is given by 31 OT 2 0 a nal ed where p 917 kg m is the sea ice density assumed to be uniform c T S is the specific heat of sea ice K T S is the thermal conductivity of sea ice pen is the flux of penetrating solar radiation a
100. meters esmf Earth System Modeling Framework driver www esmf ucar edu CICE F90 main program CICE ComponentMod F90 subroutinized version of CICE F90 for direct coupling CICE FinalMod F90 routines for finishing and exiting a run CICE InitMod F90 routines for initializing a run CICE RunMod F90 main driver routines for time stepping ice_constants F90 physical and numerical constants and parameters ice log OS SITE machine sample diagnostic output files input_templates input files that may be modified for other CICE configurations gx1 1 displaced pole grid files global_gx1 grid 1 displaced pole grid binary global_gx1 kmt 1 land mask binary ice restart file pointer for restart file name ice_in namelist input data data paths depend on particular system iced gx1 v4 0 kcatbound0 restart file used for initial condition gx3 3 displaced pole grid files global _gx3 grid 3 displaced pole grid binary 36 Numerical implementation global_gx3 kmt 3 land mask binary global_gx3 grid nc 3 displaced pole grid netCDF global_gx3 kmt ne 3 land mask netCDF ice restart file pointer for restart file name ice_in namelist input data data paths depend on particular system iced _gx3_v4 0_kcatbound0 restart file used for initial condition col column configuration files ice_in namelist input data data paths depend on particular system run_ice OS SITE machine sample script for running on the given oper
101. mponent of ice velocity value for constant albedo parameterization value for constant albedo parameterization wind speed allocatable 2D dbl_kind work array allocatable 2D integer work array allocatable 2D long integer work array allocatable 2D real_kind work array allocatable 3D real_kind work array nx_block ny_block max_blocks work array nx_block ny_block work array if true write history now e if true write initial conditions if 1 write restart now e number of years in forcing data cycle day of the year e the initial year atmospheric level height reference height for stability reference height for Tref Qref ennn gas constant water vapor gas constant air 1 1 m s 1 m s m s m s 0 68 0 77 m s 10 m 2 m 0 606 Index advection see transport albedo 4 9 27 37 41 51 52 AOMIP 57 area ice see ice fraction bilinear 25 37 blocks 36 38 40 45 46 48 50 boundary communication 40 condition 37 40 44 layer 5 6 7 thickness category 8 18 20 C grid 37 categories thickness see thickness distribution CCSM 2 3 9 27 51 52 CFL condition 11 18 32 4 column configuration 36 38 40 52 Community Climate System Model see CCSM concentration see ice fraction conductivity 29 30 30 conservation 3 10 18 20 21 30 33 43 conservation equation see transport continuity equation see transport Coriolis 5 24
102. n and ig 0 for near IR Radiation penetrating into the ice is attenuated according to Beer s Law I z lpexp K 2 57 where J z is the shortwave flux that reaches depth z beneath the surface without being absorbed and k is the bulk extinction coefficient for solar radiation in ice set to 1 4 m7 for visible wavelengths 10 A fraction exp h of the penetrating solar radiation passes through the ice to the ocean Fs Longwave radiation turbulent fluxes While incoming shortwave and longwave radiation are obtained from the atmosphere outgoing long wave radiation and the turbulent heat fluxes are derived quantities Outgoing longwave takes the standard 4 blackbody form Fry eo 1 where e 0 95 is the emissivity of snow or ice is the Stefan Boltzmann constant and T E is the surface temperature in Kelvin The longwave fluxes are partitioned such that eFL is absorbed at the surface the remaining 1 e Fr being returned to the atmosphere via Fzt The sensible heat flux is proportional to the difference between air potential temperature and the surface temperature of the snow or snow free ice F Co Oa T Cs and C below are nonlinear turbulent heat transfer coefficients described in Section 2 1 Similarly the latent heat flux is proportional to the difference between Qa and the surface saturation specific humidity Thermodynamics 29 sf Fi Ci Qa Qof Qs qi pa exp q2 T5
103. n land global climate model It was designed to be compatible with the Parallel Ocean Program POP an ocean circulation model developed at Los Alamos National Laboratory for use on massively parallel computers 38 8 9 The current version of the model has been enhanced greatly through collaborations with members of the Community Climate System Model CCSM Polar Climate Working Group based at the National Center for Atmospheric Research NCAR and researchers at the U K Met Office Hadley Centre CICE has several interacting components a thermodynamic model that computes local growth rates of snow and ice due to vertical conductive radiative and turbulent fluxes along with snowfall a model of ice dynamics which predicts the velocity field of the ice pack based on a model of the material strength of the ice a transport model that describes advection of the areal concentration ice volumes and other state variables and a ridging parameterization that transfers ice among thickness categories based on energetic balances and rates of strain Additional routines prepare and execute data exchanges with an external flux coupler which then passes the data to other climate model components such as POP This model release is CICE version 4 1 available from http climate lanl gov Models CICE It updates CICE 4 0 which was released in August 2008 e Added an alternative parameterization for ice conductivity suggested by 34 e Moved several p
104. n of solar radiation in sea ice and the upper ocean J Geophys Res Oceans 100 15 965 15 975 1995 G M Flato and W D Hibler Ridging and strength in modeling the thickness distribution of Arctic sea ice J Geophys Res Oceans 100 18611 18626 1995 C A Geiger W D Hibler and S F Ackley Large scale sea ice drift and deformation Comparison between models and observations in the western Weddell Sea during 1992 J Geophys Res Oceans 103 21893 21913 1998 W D Hibler A dynamic thermodynamic sea ice model J Phys Oceanogr 9 817 846 1979 W D Hibler Modeling a variable thickness sea ice cover Mon Wea Rev 108 1943 1973 1980 E C Hunke Viscous plastic sea ice dynamics with the EVP model Linearization issues J Comput Phys 170 18 38 2001 E C Hunke and C M Bitz Age characteristics in a multidecadal arctic sea ice simulation J Geophys Res 114 C08013 doi 10 1029 2008JC005186 2009 E C Hunke and J K Dukowicz An elastic viscous plastic model for sea ice dynamics J Phys Oceanogr 27 1849 1867 1997 REFERENCES 75 18 E C Hunke and J K Dukowicz The Elastic Viscous Plastic sea ice dynamics model in general orthogonal curvilinear coordinates on a sphere Effect of metric terms Mon Wea Rev 130 1848 1865 2002 19 E C Hunke and J K Dukowicz The sea ice momentum equation in the free drift regime Technical Report LA UR 03 2219 Los Alamos National Laboratory
105. ng factor for near ir diffuse albedo weighting factor for near ir direct albedo oo weighting factor for visible diffuse albedo weighting factor for visible direct albedo o oo number of cell averaged fields that can be written to history file 0 001 remap 0 3 m 0 15 0 06 0 66 radians radians 0 05 0 36218 0 63282 0 00182 0 00318 81 60 B bignum block block_id block_size_x y blockGlobalID blockLocalID blockLocation blocks_ice C c n calc_strair calc_Tsfc Cf char_len char_len_long check_step check_umax cldf cm_to_m coldice coldsnow conduct congel COSw coszen Cp cp air cp_ice cp_ocn cp wv cp063 cp455 Cs Cstar D daidtd daidtt dalb_mlt dalb_mlti dalb_mltv dardg1dt dardg2dt Index of primary variables and parameters Alar ire zeeen arn ate aas ne data type for blocks global block number number of cells along x y direction of block global block IDs local block IDs processor location of block local block IDs real n e if true calculate wind stress oo e if true calculate surface temperature ratio of ridging work to PE change in ridging length of character variable strings length of longer character variable strings time step on which to begin writing debugging data if true check for ice speed gt umax_stab cloud fraction
106. nts cu cg and cq are initialized as eee ee Inn Zer Zice and updated during a short iteration as they depend upon the turbulent scales Here z is a reference height of 10 m and Zice is the roughness length scale for the given sea ice category Y is constrained to have magnitude less than 10 Further defining y 1 16r and x gt 1 the integrated flux profiles for momentum and stability in the unstable Y lt 0 case are given by bm 21n 0 5 1 x Im 0 5 1 x2 2tanx Sl p 2m 051 22 In a departure from the parameterization used in 22 we use profiles for the stable case following 21 Um Ys 0 7Y 0 75 Y 14 3 exp 0 35Y 10 7 The coefficients are then updated as fe tome Cu a ST es Oa IR eye 1 co A Ys K E where A In zo zref The first iteration ends with new turbulent scales from equations 1 After five iterations the latent and sensible heat flux coefficients are computed along with the wind stress Ci Pa Loap F Tice U Cy Cs Patpu cg 1 5 pau Ua Ta gt gt Ua where Lap and Lice are latent heats of vaporization and fusion pa is the density of air and c is its specific heat Again following 21 we have added a constant to the sensible heat flux coefficient in order to allow some heat to pass between the atmosphere and the ice surface in stable calm conditions The atmospheric reference temperature T7 f is computed fro
107. nts lying halfway between the midpoint and the three vertices 7 use this formula to compute transports of the product p T which is analogous to ice volume Equation 23 does not work for ice and snow energies which are cubic functions products of area thickness and enthalpy Integrals of a cubic polynomial over a triangle can be evaluated using a four point formula 42 9 Dn Fy Ar ig xo A 48 foe 24 where x 3x9 2x 5 To evaluate functions at specific points we must compute many products of the form a x h x and a x h x q x where each term in the product is the sum of a cell center value and two displacement terms In the code the computation is sped up by storing some sums that are used repeatedly 3 1 4 Updating state variables Finally we compute new values of the state variables in each ice category and grid cell The new fractional ice areas a i j are given by 1 he We eroi Fy i 1 3 Fa i J Fa i j 1 Fa j din id amli j EEI Th Sen Dd as 18 Model components where Fur i j and Fan i 7 are the area transports across the east and north edges respectively of cell i j and A i 7 is the grid cell area All transports added to one cell are subtracted from a neighboring cell thus 25 conserves total ice area The new ice volumes and energies are computed analogously New thicknesses are given by the ratio of volume to area and enthalpies by the ratio of energy to volume
108. odel and passed to CICE calc_Tsfc false see Section 2 In this case there is no surface flux equation The top layer temperature is computed by an equation similar to 69 but with the first term on the LHS replaced by n Fet and moved to the RHS The main drawback of treating the surface temperature and fluxes explicitly is that the solution scheme is no longer unconditionally stable Instead the effective conductivity in the top layer must satisfy a diffusive CFL condition 5 x PC K lt Ar For thin layers and typical coupling intervals 1 hr K may need to be limited before being passed to the atmosphere via the coupler Otherwise the fluxes that are returned to CICE may result in oscillating highly inaccurate temperatures The effect of limiting is to treat the ice as a poor heat conductor As a result winter growth rates are reduced and the ice is likely to be too thin other things being equal The values of hs min and At must therefore be chosen with care If hs min is too small frequent limiting is required but if hs_min is too large snow will be ignored when its thermodynamic effects are significant Likewise infrequent coupling requires more limiting whereas frequent coupling is computationally expensive This completes the specification of the matrix equations for the four cases We compute the new tem peratures using a tridiagonal solver After each iteration we check to see whether the following conditions hold 1
109. ompiler flag is available to do this but 1f not try uncommenting the block of code at the end of subroutine stress in ice dyn evp F90 You will take a hit for the extra computations but it will not be as bad as running with the underflows 5 3 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 F90 debug which can be substituted for CICE 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 an
110. ons following 3 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 3 3 Mechanical redistribution 21 3 3 Mechanical redistribution The last term on the right hand side of 2 is y 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 43 35 14 11 and 26 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 ap 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 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 b h If krdg_partic 0 in the namelist we follow 43 and set b h Ad ah 35 0 otherwise where G h is the fractional area covered by ice thinner than h and G is an empirical constant Integrating ap h between category boundaries H _ and Hn we obtain the mean value of ap in category n Gn 1 gt
111. or managing ice thickness distribution ice_kinds_mod F90 basic definitions of reals integers etc ice_1v1 F90 handles most work associated with the level ice area and volume tracers 4 2 Grid boundary conditions and masks 37 ice mechred F90 mechanical redistribution component ridging ice_meltpond F90 melt pond parameterization ice ocean F90 mixed layer ocean model ice_orbital F90 orbital parameters for Delta Eddington shortwave parameterization ice_read_write F90 utilities for reading and writing files ice restart F90 read write core restart file ice_restoring F90 basic restoring for open boundary conditions ice_shortwave F90 shortwave and albedo parameterizations ice spacecurve F90 space filling curves distribution method ice_state F90 essential arrays to describe the state of the ice ice_step_mod F90 routines for time stepping the major code components ice therm itd F90 thermodynamic changes mostly related to ice thickness distribution ice therm vertical F90 vertical growth rates and fluxes ice_transport_driver F90 driver for horizontal advection ice_transport_remap F90 horizontal advection via incremental remapping ice_work F90 globally accessible work arrays rundir execution or run directory created when the code is compiled using the comp_ice script gx3 4 2 cice code executable compile directory containing object files etc grid horizontal grid file from cice input_templates gx3 ice log ID diagnostic
112. or netcdf so in cases where a variable is written at more than one frequency the variable name is appended with the frequency in files after the first one In the example above meltb is called me1tb in the monthly file for backward compatibility with the default configuration and me1tb_h in the 6 hourly file Using the same frequency twice in hist freq will have unexpected consequences and currently will cause the code to abort It is not possible at the moment to output averages once a month and also once every 3 months for example If write ic is set to T in ice in a snapshot of the same set of history fields at the start of the run will be written to the history directory in iceh_ic timeID nc da Nine history variables are hard coded for instantaneous output regardless of the averaging flag at the frequency given by their namelist flag The normalized principal components of internal ice stress are computed in principal_stress and writ ten to the history file This calculation is not necessary for the simulation principal stresses are merely computed for diagnostic purposes and included here for the user s convenience 4 5 2 Diagnostic files Like hist freq the parameters diagfreq and diagfreq_n can be used to regulate how often output is written to a log file The log file unit to which diagnostic output is written is set in ice fileunits F90 If diag type stdout then it is written to standard out or to ice log ID if you r
113. oth spectral bands decreases by 0 075 as Tp rises from 1 C to 0 C The albedo of cold snow Tyr lt 1 C is 0 98 for visible wavelengths and 0 70 for near IR wavelengths The visible snow albedo decreases by 0 10 and the near IR albedo by 0 15 as Tf increases from 1 C to 0 C The total albedo is an area weighted average of the ice and snow albedos where the fractional snow covered area is hs Fsnow TER Pa a hs hsnowpatch 28 Model components wei 1 C max snow 4 effective albedo O C max snow 4 we2 1 C bare ice Tee 0 C bare ice 4 gt 0 0 0 0 0 2 0 4 0 6 0 8 1 0 a ice thickness m Figure 4 Albedo as a function of ice thickness and temperature for the two extrema in snow depth for the default CCSM3 shortwave option Maximum snow depth is computed based on Archimedes Principle for the given ice thickness These curves represent the envelope of possible albedo values and Asnowpatch 0 02 m The envelope of albedo values is shown in Figure 4 This albedo formulation incorporates the effects of melt ponds implicitly the explicit melt pond parameterization is not used in this case The net absorbed shortwave flux is Fswabs X 1 aj Few jj Where the summation is over four radiative categories direct and diffuse visible direct and diffuse near infrared The flux penetrating into the ice is lo io Fswabs Where io 0 70 1 fsnow for visible radiatio
114. r than either square option and is sometimes necessary to keep the processor domain width sufficiently large e For small numbers of processors distribution_type rake or spacecurve with distribution_wght latitude and processor_shape square ice e The cross over point i e the number of processors where these choices result in about the same timings will vary depending on the grid and the processor architecture Figure 7 shows the model s scaling characteristics for up to 64 processors on the Linux cluster coyote These runs were made on the gx1 global grid 320x384 approximately 1 degree resolution from the restart file included with the code distribution On the gx1 grid the scaling is nearly linear up to 64 processors for slender x1 or x2 decompositions More than 64 processors do not have enough ice work to do on this grid and the timings do not scale as well The lower plot shows timings for 8 processor block distributions versus the total number of blocks illustrating the performance gains that are available by redistributing smaller blocks across processors Throughout the code i j loops have been combined into a single loop often over just ocean cells or those containing sea ice This was done to reduce unnecessary operations and to improve vector perfor mance 4 8 Adding things 4 8 1 Timers Timing any section of code or multiple sections consists of defining the timer and then
115. rmodynamics 33 Multiplying by p to change the units from J kg to J m and adding a term for the energy needed to raise the meltwater temperature to 0 C we obtain the sea ice enthalpy T q T S pi colTm T Lo 1 cu Tn 70 Note that 70 is a quadratic equation in T Given the layer enthalpies we can compute the temperatures using the quadratic formula b yb dac T 2a where a Co qi Cu E co Tm E L Pi C LoEm The other root is unphysical Melting at the top surface is given by qoh Fo Fa At if Fo gt Fa 0 otherwise CED where q is the enthalpy of the surface ice or snow layer recall that g lt 0 and dh is the change in thickness If the layer melts completely the remaining flux is used to melt the layers beneath Any energy left over when the ice and snow are gone is added to the ocean mixed layer Ice cannot grow at the top surface due to conductive fluxes however snow ice can form New snowfall is added at the end of the thermodynamic time step Growth and melting at the bottom ice surface are governed by qoh Tek Fo a Foot At 72 where Fot is given by 58 and Fi is the conductive heat flux at the bottom surface Kina Ah If ice is melting at the bottom surface q in 72 is the enthalpy of the bottom ice layer If ice is growing q is the enthalpy of new ice with temperature T and salinity Sar This ice is added to the bottom layer If the l
116. s EVP 4 Advectn horizontal transport 5 Column all vertical column processes 6 Thermo vertical thermodynamics 7 Shortwave SW radiation and albedo 8 Ridging mechanical redistribution 9 Cat Cony transport in thickness space 10 Coupling sending receiving coupler messages 11 ReadWrite reading writing files 12 Diags diagnostics log file 13 History history output 14 Bound boundary conditions and subdomain communications Table 5 CICE timers places in the code including the dynamics and advection routines The timers use MPI_WTIME for parallel runs and the F90 intrinsic system_clock for single processor runs 4 5 3 Restart files A binary unformatted file is created that contains all of the data that CICE needs for a full restart The file name begins with the character string dumpfile and the restart dump frequency is given by dumpfreq and dumpfreq_n The pointer to the filename from which the restart data is to be read for a continuation run is set in pointer_file 4 6 Execution procedures To compile and execute the code in the source directory 1 Download the forcing data used for testing from the CICE website http climate lanl gov Models CICE 2 Create Macros and run_ice files for your particular platform if they do not already exist type uname s at the prompt to get OS 3 Alter directories in the script comp ice 4 Run comp ice to set up the run directory and make the executable cice 5 To
117. s and temperatures enthalpies at time m the thermodynamic model advances these quantities to time m 1 The calculation proceeds in two steps First we solve a set of equations for the new temper atures as discussed in Section 3 5 2 Then we compute the melting if any of ice or snow at the top surface and the growth or melting of ice at the bottom surface as described in Section 3 5 3 We begin by describing the surface forcing parameterizations which are closely related to the ice and snow surface temperatures 3 5 1 Thermodynamic surface forcing The net energy flux from the atmosphere to the ice with all fluxes defined as positive downward is Fo F Fi Fry Fry 1 a 1 to Few Thermodynamics 27 where F is the sensible heat flux F is the latent heat flux Fz is the incoming longwave flux Fr is the outgoing longwave flux Fsw is the incoming shortwave flux a is the shortwave albedo and t is the fraction of absorbed shortwave flux that penetrates into the ice The albedo may be altered by the presence of melt ponds Shortwave radiation Delta Eddington Two methods for computing albedo and shortwave fluxes are available the default cesm3 method described below and a multiple scattering radiative transfer scheme that uses a Delta Eddington approach Inherent optical properties IOPs for snow and sea ice such as extinction coefficient and single scatter ing albedo are prescribed based on
118. s format description recommended value Physical Parameterizations continued 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 Forcing 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 netCDF atmo forcing files bin bin read direct access binary files atm data type default constant values defined in the code ecmwf ECMWF forcing data near NCAR bulk forcing data LYq AOMIP Large Yeager forcing data monthly monthly 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 calc_Tsfce true false calculate surface temperature true precip units mks liquid precipitation data units mm_per_month mm_per_sec same as MKS units TE FZpPE constant freezing temperature 1 8 C linear S linear function of salinity profile ustar_min real minimum value of ocean friction velocity 0 0005 m s update_ocn_f true include frazil water salt fluxes in ocn fluxes false do not include when coupling with POP oceanmixed_ice true false active ocean mixed layer calculation true if uncoupled ocn_data_format ne read netCDF ocean forcing files bin bin read direct access
119. sphere 5 where a is the sum of fractional ice areas for each category of ice The ice fraction is used by the flux coupler to merge fluxes from the ice model with fluxes from the other components For example the penetrating shortwave radiation flux weighted by a is combined with the net shortwave radiation flux through ice free leads weighted by 1 a to obtain the net shortwave flux into the ocean over the entire grid cell The flux coupler requires the fluxes to be divided by the total ice area so that the ice and land models are treated identically land also may occupy less than 100 of an atmospheric grid cell These fluxes are per unit ice area rather than per unit grid cell area In some coupled climate models for example recent versions of the U K Hadley Centre model the surface air temperature and fluxes are computed within the atmosphere model and are passed to CICE In this case the logical parameter calc_Tsfc in ice_therm_vertical is set to false The fields fsurfn the net surface heat flux from the atmosphere flatn the surface latent heat flux and fcondtopn the conductive flux at the top surface for each ice thickness category are copied or derived from the input coupler fluxes and are passed to the thermodynamic driver subroutine thermo_vertical At the end of the time step the surface temperature and effective conductivity 1 e thermal conductivity divided by thickness of the top ice snow layer in eac
120. ssump tions in the model this is not a single column but a small array of columns minimum grid size is 5x5 However the code is set up so that only the single central column is used all other columns are designated as land The column is located near Barrow 71 35N 156 5W Options for choosing the column configu ration are given in comp ice choose RES col and in the namelist file input templates col ice in Here istepo and the initial conditions are set such that the run begins September 1 with no ice The grid type is rectangular dynamics are turned off kdyn 0 and one processor is used History variables available for column output are ice and snow temperature Tinz and Tsnz These variables also include thickness category as a fourth dimension 4 2 4 Boundary conditions Open boundary conditions are the default in CICE the physical domain can still be closed using the land mask In our bipolar displaced pole grids one row of grid cells along the north and south boundaries is located on land and along east west domain boundaries not masked by land periodic conditions wrap the domain around the globe CICE can be run on regional grids with open boundary conditions Except for variables describing grid lengths non land ghost cells along the grid edge are filled using Neumann reflecting conditions as would be appropriate for an open boundary in mid ocean The ice state can be restored to data values along open boundaries also
121. t depth z and z is the vertical coordinate defined to be positive downward with z 0 at the top surface If shortwave default the penetrating radiation is given by Beer s Law Ipen 2 Jo exp 12 30 Model components where Ip is the penetrating solar flux at the top ice surface and 1 is an extinction coefficient If shortwave dEdd then solar absorption is computed by the Delta Eddington scheme The specific heat of sea ice is given to an excellent approximation by 33 Lops ee where cy 2106 J kg deg is the specific heat of fresh ice at 0 C Lo 3 34 x 10 J kg is the latent heat of fusion of fresh ice at 0 C and u 0 054 deg psu is the ratio between the freezing temperature and salinity of brine Following 45 and 31 the standard thermal conductivity conduct MU71 is given by Ki T S Ko zn 62 where Ko 2 03 W m deg is the conductivity of fresh ice and 4 0 13 W m psu is an empirical constant Experimental results 44 suggest that 62 may not be a good description of the thermal conductivity of sea ice In particular the measured conductivity does not markedly decrease as T approaches 0 C but does decrease near the top surface regardless of temperature An alternative parameterization based on the bubbly brine model of 34 for conductivity is available conduct bubbly ci T S co 61 K 2 11 0 011T 0 098 T 63 0 where p and po 917 kg m are
122. tability 3 5 6 21 36 41 42 state variables 3 4 8 17 18 37 strain rate 2 22 24 strength 2 23 stress ice ocean see ice ocean stress principal 25 43 tensor see internal stress wind see wind stress subcycling 25 42 sublimation see evaporation surface height see height sea surface temperature 26 34 59 atmospheric 4 freezing 7 ice 9 26 33 ocean 4 7 potential 4 5 reference 4 6 surface 5 6 9 13 thermodynamics 26 34 thickness distribution 2 8 9 11 18 23 26 36 37 42 ice or snow 7 9 11 13 18 24 26 33 34 37 space see transport thickness timers 36 43 48 51 tracer 9 tracers 10 11 13 36 49 50 65 66 transport 2 8 20 37 41 42 44 horizontal 9 18 thickness 18 20 tripole 38 39 turbulent fluxes latent heat see latent heat sensible heat see sensible heat wind stress see wind stress upwind 10 INDEX van Leer 13 velocity ice 2 6 9 11 13 16 24 26 37 38 volume ice or snow 8 10 12 17 19 22 23 37 water open see open water wind stress 4 5 6 24 velocity 4 6 42 WMO 8 zero layer model 29 73 74 REFERENCES References 1 2 Ly 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 T L Amundrud H Malling and R G Ingram Geometrical constraints on the evolution of ridged sea ice J Geophys Res 109 2004 C06005 doi 10 1029 2003JC002251 D Bailey
123. te In practice the ratio Ate T At 1 40 120 provides both stability and acceptable efficiency for time steps At on the order of 1 hour Note that only T and At figure into the stability of the dynamics component At does not Although the time step may not be tightly limited by stability considerations large time steps eg At 1 day given daily forcing do not produce accurate results in the dynamics component The reasons for this error are discussed in 17 see 20 for its practical effects The thermodynamics component is stable for any time step as long as the surface temperature T fe is computed internally 4 5 Model output 4 5 1 History files Model output data is averaged over the period s given by hist freq and histfregq_n and written to binary or netCDF files prepended by history file inice in That is if history file iceh then the filenames will have the form iceh timeID nc or iceh timeID da depending on the output file format history_ format chosen The netCDF history files are CF compliant header information for data contained in the netCDF files is displayed with the command ncdump h filename nc With binary files a separate header file is written with equivalent information Standard fields are output accord ing to settings in the icefields nml namelist in ice in The user may add or subtract variables not already available in the namelist by following the instructions in section 4 8 2 The history
124. tead of added to the new ridge The net area removed by ridging and closing is a function of the strain rates Let Rnet 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 11 Rnet is given by Cs Rnet z Dp min Dp 0 44 where Cs is the fraction of shear dissipation energy that contributes to ridge building Dp is the divergence and A is a function of the divergence and shear These strain rates are computed by the dynamics scheme The default value of C is 0 25 Next define Riot ys rn This rate is related to Rnet by N 1 ort Dern 1 E Rar 45 n 1 Rnet Mechanical redistribution 23 Given Rnet from 44 we use 45 to compute Rio Then the area ridged in category n is given by arn Tn At where rn apn Rtot The area of new ridges is arn kn and the volume of new ridges is arn hy since volume is conserved during ridging We remove the ridging ice from category n and use 39 and 40 or 42 and 43 to redistribute the ice among thicker categories Occasionally the ridging rate in thickness category n may be large enough to ridge the entire area an during a time interval less than At In this case Rot is reduced to the value that exactly ridges an area an during At After each ridging iteration the total fractional ice area a is computed If a gt 1 the ridging is repeated with a value of Rnet suffic
125. th ice ses French like cease or she 1 s o Shii aisu Japanese 2 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 CCSM Flux Coupler 22 from NCAR the first major climate model to incorporate CICE The flux coupler was originally intended to gather state variables from the component models com pute 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 1s 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
126. the end of the timestep simplifying the code a great deal However the Delta Eddington radiative scheme computes albedo and shortwave components simultaneously and in order to have the most up to date values available for the coupler at the end of the timestep the order of radiation calculations needed to be shifted Albedo and shortwave components are computed after the ice state has been modified by both thermodynamics and dynamics so that they are consistent with the ice area and thickness at the end of the step when sent to the coupler However they are computed using the downwelling shortwave from the beginning of the timestep Rather than recompute the albedo and shortwave components at the beginning of the next timestep using new values of the downwelling shortwave forcing the shortwave components computed at the end of the last timestep are scaled for the new forcing 3 1 Horizontal transport We wish to solve the continuity or transport equation Odin Ot V aimu 0 3 10 Model components for the fractional ice area in each thickness category n Equation 3 describes the conservation of ice area under horizontal transport It is obtained from 2 by discretizing g and neglecting the second and third terms on the right hand side which are treated separately Sections 3 2 and 3 3 There are similar conservation equations for ice volume snow volume ice energy and snow energy Sun Y vmu 0 4 Cue
127. the ice model may be coupled to another system without using MPI 4 1 Directory structure The present code distribution includes make files several scripts and some input files The main directory is cice and a run directory rundir is created upon initial execution of the script comp ice One year of atmospheric forcing data is also available from the code distribution web site see the README file for details cice README v4 0 basic information Directory structure 35 bld makefiles Macros OS SITE machine macro definitions for the given operating system used by Makefile OS Makefile OS primary makefile for the given operating system std works for most systems makedep c perl script that determines module dependencies clean_ice script that removes files from the compile directory comp_ice script that sets up the run directory and compiles the code csm_share modules based on shared code in CCSM shr_orb_mod F90 orbital parameterizations doc documentation cicedoc pdf this document PDF PDF documents of numerous publications related to CICE drivers institution specific modules cice4 official driver for CICE v 4 LANL CICE F90 main program CICE F90_debug debugging version of CICE F90 CICE FinalMod F90 routines for finishing and exiting a run CICE InitMod F90 routines for initializing a run CICE RunMod F90 main driver routines for time stepping ice constants F90 physical and numerical constants and para
128. the ridge volume going to category m is la a a 40 a Hmax Hmin l This uniform redistribution function tends to produce too little ice in the 3 5 m range and too much ice thicker than 10 m 1 Observations show that the ITD of ridges is better approximated by a negative exponential Setting krdg_redist 1 gives ridges with an exponential ITD 26 gr h Xx exp h Hmin 41 for h gt Hmin with gr h 0 for h lt Hmin Here A is an empirical e folding scale and Hmin 2an where hn is the thickness of ridging ice We assume that A uhl 2 where u mu_rdg is a tunable parameter with units m Thus the mean ridge thickness increases in proportion to hal E as in 14 The value y 4 0 ml gives A in the range 1 4 m for most ridged ice Ice strengths with y 4 0 m and krdg_redist 1 are roughly comparable to the strengths with A 50 m and krdg_redist 0 From 41 it can be shown that the fractional area going to category m as a result of ridging is ta 7 exp Hm 1 Hain exp Hm z Fis Al 42 The fractional volume going to category m is vol Hm 1 A exp Hin Hmin A Hm A expl Hm Hmin A 43 Equations 42 and 43 replace 39 and 40 when krdg_redist 1 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 ins
129. ts degrees N latitude of T U grid cell corners degrees N latent heat of melting of fresh ice Lsub Lvap J kg transfer coefficient for latent heat northern southern hemisphere mask local address of block in current distribution definition of a logical variable o o oooooooonmmmm m kind true Index of primary variables and parameters lonpt lont u _bounds Lsub Itripole_grid Lvap M m_min m_to_cm ml m2 m2_to_km2 master_task max_blocks max_ntrer maxraft mday meltb meltl melts meltt min salin mlt_onset month monthp mps_to_cmpdy mtask mu_rdg my_task N nblocks nblocks_tot nblocks_x y ncat ndte ndyn_dt new_day new_hour new_month new_year nghost ngroups nhlat nilyr nprocs npt e longitude of diagnostic points longitude of T U grid cell corners latent heat of sublimation for fresh water flag to signal use of tripole grid latent heat of vaporization for fresh water minimum mass for computing velocity meters to cm conversion ensen eenen constant for lateral melt rate 00 eee eee eee m to km conversion creek ie eG NEE MEERN task ID for the controlling processor maximum number of blocks per processor maximum number of tracers available maximum thickness of ice that rafts day of the month basal ice melt lateral ice melt snow melt top ice melt threshold for brine pockets
130. u dyn_dt dyt dyu dvidtd dvidtt dvirdgdt E ecci eice n emissivity eps11 eps13 eps16 esno n evap evp_damping ew_boundary_type eyc number of days in one month day number at end of month e number of days in one year definition of double precision 0 0 e eee eee eee e write extra diagnostics function of strain rates see Section 3 4 ratio of freezing temperature to salinity of brine e diagnostic output file alternative to standard out e where diagnostic output is written e how often diagnostic output is written 10 once per 10 dt block distribution information e cartesian or rake or spacecurve e weighting method used to compute work per block strain rate I component velocity divergence divergence associated with advection drag coefficient for water on ice drag coefficient for water on ice pw e thermodynamics time step see ice_shortwave F90 subcycling time step for EVP dynamics Ate 1 dte where dte is the EVP subcycling time step e output file for restart dump e dump frequency for restarts y m or d e restart output frequency width of T cell Az through the middle width of U cell Ax through the middle dynamics and transport time step Atgyn height of T cell Ay through the middle height of U cell Ay through the middle ice volume tendency due to dynamics transport ice volume tendency due to thermodynamics ice volume ridging rate yield curve minor major axis rat
131. ume 6 CICE InitMod F90 initialize tracer includes reading restart file 50 7 8 9 5 TROUBLESHOOTING CICE_RunMod F90 ice_step_mod F90 e call routine to write tracer restart file e call physics routines in ice_ tracer F90 ice_history F90 add history variables Section 4 8 2 ice_in add namelist variables to tracer_nml and icefields_nml 5 Troubleshooting 5 1 Initial setup The script comp ice is configured so that the files grid kmt ice in run ice iced_gx3_v4 0_kcatbound0 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 al ready exist type uname s at the prompt
132. vertices 14 Model components Figure 1 In incremental remapping conserved quantities are remapped from the shaded departure region a quadrilateral formed by connecting the backward trajectories from the four cell corners to the grid cell labeled H The region fluxed across the north edge of cell H consists of a triangle abc in the NW cell and a quadrilateral two triangles acd and ade in the N cell b and c relative to the cell corner vertex a using Euclidean geometry to find vertex c Then we translate the three vertices to a coordinate system centered in the NW cell This translation is needed in order to integrate fields Section 3 1 3 in the coordinate system where they have been reconstructed Section 3 1 1 Repeating this process for the north and east edges of each grid cell we compute the vertices of all the departure triangles associated with each cell edge Figure 2 reproduced from 7 shows all possible triangles that can contribute fluxes across the north edge of a grid cell There are 20 triangles which can be organized into five groups of four mutually exclu sive triangles as shown in Table 3 In this table x1 y1 and 2 y2 are the Cartesian coordinates of the departure points relative to the northwest and northeast cell corners respectively The departure points are joined by a straight line that intersects the west edge at 0 ya relative to the northwest corner and intersects the east edge at 0 y relative
133. ward That is Xp uAt 18 where xp is the location of the departure point relative to the cell corner and u is the velocity at the corner This approximation is only first order accurate Accuracy can be improved by estimating the velocity at the midpoint of the trajectory Horizontal transport 17 3 1 3 Integrating fields Next we integrate the reconstructed fields over the departure triangles to find the total area volume and energy transported across each cell edge Area transports are easy to compute since the area is linear in x and y Given a triangle with vertices xj zi yi 1 2 3 the triangle area is 1 Ar z ea 21 ys y1 ya y1 z3 zi 19 The integral F of any linear function f r over a triangle is given by Fa Arf xo 20 where xo Zo Yo is the triangle midpoint 1 3 AQ 21 3 i 1 To compute the area transport we evaluate the area at the midpoint a xo Ae AzXQ Ayyo 22 and multiply by Ar By convention northward and eastward transport is positive while southward and westward transport is negative Equation 20 cannot be used for volume transport because the reconstructed volumes are quadratic functions of position They are products of two linear functions area and thickness The integral of a quadratic polynomial over a triangle requires function evaluations at three points A 3 Fr gt df 23 i 1 where x xo x 2 are poi
134. y also include a border of cells around the ice pack for numerical purposes These masks are used in the dynamics component to prevent unnecessary calculations on grid points where there is no ice They are not used in the thermodynamics component so that ice may form in previously ice free cells Like the land masks hm and uvm the ice extent masks ice tmask and ice umask are for T cells and U cells respectively Two additional masks are created for the user s convenience 1mask_n and lmask s can be used to compute or write data only for the northern or southern hemispheres respectively Special constants spval and spval _db1 each equal to 10 are used to indicate land points in the history files and diagnostics 4 3 Initialization and coupling 41 ice_ic runtype restart initial false initial true continue true or false none no ice no ice restart using pointer file default SST latitude dependent SST latitude dependent restart using pointer_file filename no ice start from filename restart using pointer_file Table 4 Ice initial state resulting from combinations of ice_ic runtype and restart If false restart is reset to true restart is reset to false ice_ic is reset to none 4 3 Initialization and coupling The ice model s parameters and variables are initialized in several steps Many constants and physical parameters are set in ice_constants F90 Namelist variables Table 7 whose v

CICE: the Los Alamos Sea Ice Model Documentation and

Contents

Download Pdf Manuals

Related Search

Related Contents