,," ,

"

'.

. .~

,'! ...

, 

" "',,4

...-"

, ,

r :,

-'

"

I,

"

. I..' '~,

,

"

\, .1

"

" , r ......

.. .

, ,,',' ;> 'I •.#

"

... !:

~, ' OJ' ...

, ,

- " , -

.,

, ,

.,

,',

.1".

• I ...

" "

"

, .

'lo.,'

I',

"

, .'

'1

.', .

,

I

" ,

,"', ;.~'

.'

'\'.

r ,~

"

,

,.' ,

... \'

.. ' \".r

,,',.,

','

.'

'\', .

"

\.

"

. ' _',.1

',,;,'

j

t··;,..

" 1., I~ 0

.\

'. \

\,'

.. .... ", ''-~,

• '. t

, ... '~",). .

.'

" ,

'.' .

.,

i'

"'I.,

, .

, .~.\'

Jo r. .....<f.... c'

.'

.,

"

; ,'.. .. 't- '.,,' '.

"'." .

,

,

'-.

"

. , .)' ,,~

"

l.'·...•

J"

0,> •

,. I

'C

_ .J

.'....

;..

•." .

. J ,

.. t .. ' •• 1

1.

., .

.. ~',~ -'\..,

: .

'. )

, ;:,~ .

"

, '.

'.,

.<,

'"

, , ,',

. i"~, :," ".k.

."'.",

,. .

,....

.~

.•. ''i,

'.'

, f.·'

' .. ":

I.'

,- .

, , ,"

:,'-:,

,. ,

"

" '"I

,~~ :..'

\ .

",-.. 'fI,.~. ",

'1

, I

.... .

~., .'

" .;: ",.~",~) "0

.. \

... ~:.. (

'.:....'~...."'~. ! " . .-'>,

.... ).' / y , c'" f

,~ ..

,'\,

:','

., ;,

'...(

.... '.

1 .

, '

"

\.

~)

'(',1;

"",

! '.'''1

.'".-

'·t,

, " .

.)

, '

'r

"

"

, \

.f'" ..: 'I'

'., ! ; ~. '-- , ~.:;. ,

; ,

",

..t' f. .. .-

,~,

.,'

......

'.~

',),,'

, ··r',

.', .

, ,.'

.': r

"

,~ \

, " . " .,"

. tJ.S~ Departineritof't1).ei~t~rl.or . Mitier;als.Man:agem·~ritService,

:',;Ail.,chorage, 'Al~ka' i . ."

. :~

.. • . ,0. :-"' ~. '. ;' • V" 1 .r,

"C6htractN6~' i4-35:0r~96~CT-3'6818 ,::-,,, ." '<," :\

"'.

:' r' <\

.~." '. ,

....:.:.:

, 1'"

,t·._

'.il

. i

,,'

',' ,

" .,' '-', ,Kafherirte;$:'Hedstr'oth:"~, . ...:. '. J..., ',' 1\,' '. • ·....v ~_.'-' '

'~risti tut,E! :0~. ~a;~ihe 'aqd:~mistaJ;Scien~~s.

. ",'J3..ut~ers"UJiiversity,.'. ... (. ' " .~ "

. '. '.

'.";' ,

"

, '. ~.

'.

,,'.:: 'J""

or,' •

i..' "

,.' '

'-':,

.. '.'

,,.

.,'):

\J

.,

'.

.... ,

'4- <:•

. '"

,. .~'

..: .. ."

_', ":~\,,,' " '. " ',',", ','" ":> ,: ';'" ~._~;. ". "".; ,', 'TechriicaI,M~nu.alfor a" C.oup~~~"S~~~Ic<¢10c'~an:·:

.~' ,',' -;~'~ ~ir~ul~tioD;;:M9del '(Versi(jJ;l' 2,) '~<, ',', • '

,~ \ '.~~.- - ~.':' ~,'.. .~' .. '" • -. ,.,: , •.; oj '<" -', '\·~.f:·} •. ~"

. '~.

',',: ,~

.,

, \' \;,

,','

, . '-.

...,.... '. '~

'/\1 .'

, . ( '1,·

·i

....

"

,.

r:,,- ..,/

.... "

.......

;: 1"', .'.'

',~ 'I

'.'

.. ,.,j

,j

"

",.,- " ... •• C'

, "f' 1"." , ~.(' .1

.,

"

, .

,"I.; T'

" . ~. ','

, ';;.

'.' "J

.,',

'r" ~.' .

"'''...'

.j- '."

~, ,.' "',

.'"

'\

.;-

:/'.- -i" .'l:~·~

J \ ... • .I . 1O~ .'

, ..'" ., ''''''.' , ....

.

~. "

;'",' /. 1,

.....

. "

/" .' • 'It'.

~, '1

"

• ,i. i';

.......

.'..'....

, ,,.. ... ,,.,. /... .

,.,.'t., • , , .-1'
I""~

-:! '

.'.

. , ,~

. )

, '

I" ,

':

" .r,~

... > f ' )

1

,1:

;-'. t

, ,

"".

",

,

"
OCS Study MMS 2000-047

Technical Manual for a Coupled Sea-Ice/Ocean

Circulation Model (Version 2)

Katherine S. Hedstrom

Institute of Marine and Coastal Sciences

Rutgers University

June 2000

This study was funded by the Alaska Outer Continental Shelf Region of the Minerals Management Service, U.S. Department ofthe Interior, Anchorage, Alaska, through Contract 14-35-01-

96-CT-30818 with Rutgers University, Institute of Marine and Coastal Sciences.

The opinions, findings, conclusions, or recommendations expressed in this report or product are

those of the authors and do not necessarily reflect the views of the U.S. Department of the Interior,

nor does mention of trade names or commercial products constitute endorsement or recommendation for use by the Federal Government.
".••••••••••••••••••••••••••••••••••••••••••

.,.
•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.'.•

•

•

•

•

'.

,.

•

•

•

•

•

•

•

•

•

_e•

~I.

This document was prepared with Jb.'fEX and xfig.

Acknowledgments

The SPEM model had a free-surface version written by Dale Haidvogel and myself, which was

never adequately well-behaved. We asked Tony Song to try to fix it and Tony made many changes,

including the introduction of the vertical s-coordinate. John Wilkin therefore named it SCRUM,

the S-Coordinate Rutgers University Model. Hernan Arango has since made even more changes,

cleaning up the numerics and the code, writing new NetCDF I/O routines, etc. Bob Chant has

also contributed to the SCRUM effort, providing the Smolarkiewicz advection scheme. Bernard

Barnier and Anne-Marie Treguier convinced us to try the vertical finite difference approach which

has proven to be numerically much more stable than either the original spectral scheme or the

finite elements that Tony introduced. John Wilkin, Aike Beckmann and Dale Haidvogel provided

the rotated mixing tensors for the horizontal viscosity/diffusion. Bill Large provided us with the

code for their planetary boundary layer option and Scott Durski put it into SPEM and provided

me with its description. I owe all these people and the rest of the SPEM/SCRUM community,

especially Hernan Arango for making his notes available and for the comments describing all the

subroutines and variables in the model.

Bill Hibler first came up with the viscous-plastic rheology we are using. Paul Budgell has

rewritten the dynamic sea-ice model, improving the solution procedure and making the waterstress term implicit in time. We are very grateful that he is allowing us to use his version of

the code. The sea-ice thermodynamics is derived from Sirpa Hakkinen's implementation of the

Mellor-Kantha scheme. She was kind enough to allow us to start with her code.

Thanks to the Usenet community for providing great tools like perl, patch, cpp, res, and

imake to aid in software development (and to make it more fun). Cathy Lascara talked me into

trying imake with SPEM which has been well worth the trouble.

Development and testing of the SPEM model has been funded by the Office of Naval Research

(SC-53789), the National Science Foundation (OCE 90-12754-01), the Minerals Management Service (14-35-0001-30675), and the National Aeronautics and Space Administration (GC-R-261348-

006-C).

Development and testing of the SCRUM model has been funded by the Minerals Management

Service (14-35-01-96-CT030818) and the Office of Naval Research (N00014-93-1-0758, N00014-95-

1-0457, and N00014-93-1-0197).

UNIX is a registered trademark of the Open Group.

Sun is a trademark of Sun Microsystems, Inc.

SGI is a trademark of Silicon Graphics, Inc.

This is Contribution #2000-07 of the Institute of Marine and Coastal Sciences, Rutgers University.
••••••••••••••••••••••••••••••••••••••••••••
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

~.

•

;.•

".•

•

'.•

,:. '.

•

•

•

•

•

•

•

•

"

Abstract

The S-Coordinate Rutgers University Model (SCRUM), authored by Dr. Hernan Arango et

al. of the Institute of Marine and Coastal Sciences at Rutgers University, is one approach to

regional and basin-scale ocean modeling. This user's manual for SCRUM describes the model

equations and algorithms, as well as additional user configurations necessary for specific applications. It also describes the sea~ice model that is based on Hibler's viscous-plastic dynamics

and the Mellor-Kantha thermodynamics.
••••••••••••••••••••••••

.,••••••••••••••••••

...
••••-••••••••

'.••

:.•

'.•••••

'.:.'.'.'.•••

:.•

I:."•••••••

'

Contents 1 Introduction

1.1 Acquiring the SCRUM code

1.2 The SCRUM email list

.

1.3 Future plans. . . .

1.4 Warnings and bugs ...

2 Ocean Model Formulation

2.1 Equations of motion ...

2.2 Vertical boundary conditions.

2.3 Horizontal boundary conditions

2.4s (stretched vertical) coordinate system

2.5 Horizontal curvilinear coordinates.

3 Numerical Solution Technique

3.1 Vertical and horizontal discretization.

3.1.1 Horizontal grid

3.1.2 Vertical grid

3.2 Masking of land areas

3.2.1 Velocity ....

3.2.2 Temperature, salinity and surface elevation

3.3 Conservation properties ....

3.4 Vertical viscosity and diffusion.

3.5 Depth~integrated equations

.

3.6 Time stepping: internal velocity modes, temperature, and salinity.

3.7 Advection schemes.

3.7.1 Smolarkiewicz.

3.7.2 Third-order Upwind

.

3.8 Determination of the vertical velocity and density fields

3.9 The pressure gradient terms . . .

3.10 Horizontal friction and diffusion.

3.10.1 Laplacian.

3.10.2 Biharmonic.

3.10.3 Rotated mixing tensors

3.11 Vertical mixing schemes . . . .

3.11.1 Brunt-ViiisiiHi frequency scheme

3.11.2 Pacanowski-Philander.

3.11.3 Mellor-Yamada.

3.11.4 The Large, McWilliams and Doney parameterization.

3.12 Open boundary conditions.

3.12.1 Gradient boundary condition

.

3.12.2 Radiation boundary condition.

4 Details of the Code

4.1 Main subroutines . . . . . .

4.2 Other subroutines and functions

4.3C preprocessor variables

4.4 Important parameters . . . . . .

1

12334556779

11

11

11

11

12

12

13

13

14

15

17

17

17

19

20

20

20

21

21

21

22

23

23

23

25

28

28

28

31

31

34

38

42
5 Support Programs for Initialization

5.1 Grid generation.

5.1.1 ezgrid ..

5.1.2 gridpak.

5.2 Masking.....

5.2.1 The serum_mask program

5.3 Objective Analysis . . . . . .

5.4 Forcing fields . . . . . . . . .

5.5 Initial and climatology fields.

6 Configuring SCRUM for

a Specific Application

6.1 Configuring SCRUM . . . . . . . . .

6.1.1 cppdefs.h and checkdefs.F

6.1.2 Model domain

6.1.3 x,y grid.

6.1.4 ~,17 grid.

6.1.5 Initial conditions

6.1.6 Equation ofstate.

6.1.7 Boundary conditions

6.1.8 Model forcing.

6.1.9 scrum.in.

6.1.10 User variables and subroutines

6.2 Upwelling/Downwelling Example

6.2.1 cppdefs.h ..

6.2.2 Model domain.

6.2.3 ana_grid . . . . . . . . .

6.2.4 Initial conditions and the equation of state

6.2.5 Boundary conditions

6.2.6 Model forcing.

6.2.7 scrum.in....

6.2.8 Output ....

6.3 North Atlantic example

6.3.1 cppdefs.h ..

6.3.2 Model domain

6.3.3 gridpak.....

6.3.4 Initial conditions

6.3.5 Boundary conditions

6.3.6 Forcing ...

6.3.7 Climatology.

6.3.8 scrum.in.

6.3.9 Output ...

7 Plotting Programs for Postprocessing 8 Ice Model Formulation

8.1 Model structure. . . . . . . . . . .

8.2 Horizontal curvilinear coordinates.

8.3 Numerical Scheme.

8.4 Horizontal boundary conditions . .

8.5 Thermodynamics..........

8.5.1 Ocean surface boundary conditions.

2

43

43

43

43

43

43

44

47

47

49

49

50

50

51

51

51

51

53

53

53

60

60

60

61

61

61

61

61

62

62

65

65

70

71

71

71

71

71

75

75

81

85

85

87

89

90

90

94

-.•••••••••••••••••••••••••••••••••••••••••••
•••••••••••••••

'.••

'.I.'.

:.'.•••':.:.

••'.',.'.••

'.•••

'.••••

'

8.5.2 Frazil ice formation. . . . . . . . . .

8.5.3 Differences from Mellor and Kantha

9 Description of the Ice Model and the Coupling

9.1 Ice model structure.

9.1.1 Thermodynamic subroutines

9.1.2 Initialization ..

.

9.1.3 Forcing fields ...

9.1.4 Other subroutines

9.2 Coupling strategy. . . .

9.3C preprocessor variables

A Model Timestep B The vertical s-coordinate C Horizontal curvilinear coordinates D Viscosity and Diffusion

D.1 Horizontal viscosity.

D.2 Horizontal Diffusion.

D.3 Vertical Viscosity and Diffusion

E Radiant heat fluxes

E.1 Shortwave radiation

E.2 Longwave radiation.

E.3 Sensible heat

E.4 Latent heat . . .

F The

C preprocessor

F.1 File inclusion ..

F.2 Macro substitution ..

F.3 Conditional inclusion.

F.4C comments . . . .

F.5 Potential problems

F.6 Modern Fortran.

G The patch program H Makefiles

H.1 imake .....

H.2 Your Makefile

I Perl scripts for Fortran

1.1 redo ..

1.2 findent

1.3 relabel

1.4 unenddo

1.5 ifspace

1.6 sfmakedepend.

3

95

96

97

97

97

99

99

100

100

100

103

105

107

109

109

109

109

111

111

111

111

111

113

113

113

114

115

115

116

117

119

119

. 120

123

123

124

124

124

125

125
List of Tables

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

3

11

11

12

16

32

33

44

47

52

66

5

6

86

92

94

96

112

106

121

69

72

73

74

79

91

91

98

101

4

The variables used in the description of the ocean model . . . . .. . . . . .

The variables used in the vertical boundary conditions for the ocean model

Variables used in the ice momentum equations

Variables used in the ice thermodynamics

Ocean surface variables .

Frazil ice variables . . . . . . . . . . . . . . . .

Variables used in computing the in,coming radiation and latent and sensible heat

Tree structure of anonymous ftp directory

Placement of variables on an Arakawa C grid . .

Placement of variables on staggered vertical grid

Masked region within the domain . . . . .

The split timestepping used in the model.

Flow chart of the model main program.

Flow chart of the initial subroutine. .

Small grid with masked regions . . . .

The serum_mask program in action .

The whole grid. . . . . .. . . . . . . .

The upwelling/downwelling bathymetry.

Surface velocities after one day, showing the flow to the left of the wind (southern

hemisphere). 67

Constant slices of the u, v, wand n fields at day 1. , . . .. 68

Constant slices of the T, S (tracer), kinetic energy and Ertel potential vorticity at

day 1 .

The North Atlantic grid. . . . . . . . . . . .

The raw bathymetry from etopo5. . . . . .

The smoothed North Atlantic bathymetry..

The annual mean surface elevation for year 10.

Diagram of the different locations where ice melting and freezing can occur.

Diagram of internal ice temperatures and fluxes. The hashed layer is the snow.

Flow chart for the sea-ice model. .

Flow chart for the coupled ice-ocean model. . . . . . . . . . . . . . . . . . . . .

The s-surfaces for the North Atlantic with (a) () = 0.0001 and b = 0, (b) () = 8 and

b = 0, (c) () = 8 and b = 1. (d) The actual values used in this domain were () = 5

and b = 0.4 .

24 Creating Makefiles

1

2

3

4

5

6

7

15

16

17

18

19

20

21

22

23

13

14

List of Figures

1

2

3

4

5

6

7

8

9

10

11

12
.'•

•

•

•

•

•

•

•

•

•

•

•

•

•

,.

•

:.•

•

•

•

•

•

•

:.•

:.

•

•

•

'.•

'..'•

•

'.'.•

•

'.

,.'

1 Introduction

This user's manual for the S-CoordiIiate Rutgers University Model (SCRUM) describes the model

equations and algorithms, as well as additional user configurations necessary for specific applications. Some initial tests of SCRUM are described in Song and Haidvogel [54], while others are

described in Haidvogel and Beckmann [19]. This manual also describes the sea-ice model that

we are using, derived from that of Hibler [22]. It was rewritten and coupled to SCRUM by Paul

Budgell.

The principle attributes of the model are as follows:

General

• Primitive equations with potential temperature, salinity, and an equation of state. '

• Hydrostatic and Boussinesq approximations.

• Optional third-order upwind advection scheme.

• Optional Smolarkiewicz advection scheme for tracers (potential temperature, salinity, etc.).

• Option for point sources and sinks.

Horizontal

• Orthogonal~curvilinear coordinates.

• Arakawa C grid.

• Closed basin, periodic, prescribed, radiation, and gradient open boundary conditions.

• Masking of land areas.

Vertical

• s (terrain-following) coordinate.

• Free surface.

• Tridiagonal solve with implicit treatment of vertical viscosity and diffusivity.

Ice

• Hibler viscous-plastic dynamics.

• Mellor-Kantha thermodynamics.

• Orthogonal-curvilinear coordinates.

• Arakawa C grid.

• Smolarkiewicz or third-order upwind advection of tracers.

• Optional ridging scheme.

Mixing options

• Horizontal Laplacian and biharmonic viscosity and diffusion along constant s, z or

density surfaces.

• Horizontal free-slip or no-slip boundaries.

• Vertical harmonic viscosity and diffusion with a spatially variable coefficient, with

options to compute the coefficients with Large et al. [28], Mellor-Yamada [36], or

Pacanowski-Philander [40] mixing schemes.

1
Implementation

• Dimensional in meter, kilogram, second (MKS) units.

• Fortran 77 and common extensions.

• Runs under UNIX, requires the C preprocessor.

• All input and output is done in NetCDF [47] (Network Common Data Format),

requires the NetCDF library.

• Optimized for vector computers.

• Pre- and post-processing graphics package available which uses the NCAR (National

Center for Atmospheric Research) graphics libraries.

Chapters 2 and 3 describe the model physics and numerical techniques and are an update of

the description in Song and Haidvogel [54]. Chapter 4 lists the model subroutines and functions.

Chapter 5 describes the support programs which are needed to provide SeRUM with data files.

As distributed, SCRUM is ready to run with a number of example problems. The process of

configuring SCRUM for a particular application is described in Chapter 6, including a discussion

of a few example applications. Chapter 7 describes Hernan Arango's plotting programs cnt, cent,

sec, and csec. Chapter 8 describes the ice equations while chapter 9 describes the ice subroutines

and the coupling procedure.

1.1 Acquiring the SeRUM code

The version of the model described in this document is a merger between SCRUM 4.0 and a sea-ice

model. SCRUM has since been replaces by ROMS as our ocean model of choice. If you would like

to obtain this code, please contact me (kate@imcs.rutgers.edu). Version 3.0 of SCRUM is available

from ftp://ahab.rutgers.edu/pub. When connected, you will be in the ftp directory. The directory

structure is shown in Fig. 1. Everyone has write permission in the pub/incoming directory. To

get to the SCRUM source code, type 'cd pub/serum/tars' and then 'get scrum3_f77.ta.r.gz'. Some

version 4.0 files are under pub/arango. It might be more convenient to access the version 3.0 files

through our web site:

http://marine.rutgers.edu/po/index.html

The UNIX convention is that a filename ending in .gz has been compressed with the gnu gzip

utility. The corresponding gunzip also comes with it. Likewise, a .tar ending designates a file

created with the Unix tar (tape archive) utility. The steps in unpacking these files are:

%gunzip gridpak.tar.gz

%tar xvi gridpak.tar

or

%gunzip < gridpak.tar.gz I tar xvi -

Note that both .tar and .gz files are binary and must be retrieved using binary mode in ftp.

If you are unable to acquire the code in this fashion feel free to contact Hernan Arango at:

Hernan G. Arango

Institute of Marine and Coastal Sciences

71 Dudley Road

New Brunswick, NJ 08903

(732)-932-3704

Internet: arango@imcs.rutgers.edu

2

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
3

1.2 The SeRUM email list

1.3 Future plans

• Alternate versions of the ice thermodynamics. Possible options are those of Ebert and

Curry[10] and Doug Martinson. On the ice dynamics, Bruno Tremblay is planning to convert

his model to a curvilinear version for parallel computers.

~mMk mexcdf4

mexcdf5

tools

serum ---....,..-- Examples-:·-------r--- Basin

incoming bin CanyonA

util em -'--, Patches CanyonB

gridpak config Grav

perl contrib Tasman

arango doc Upwelling

£77 --' Patches

forcing

grid

hydro

initial

lib

matlab

~:tcdf-f---------,Incedit

plot--c Data L tools

tars _L C Palettes

xcoast Patches

Figure 1: Tree structure of anonymous ftp directory

pub

• Lagrangian floats from SPEM. This is just a matter of time.

• ROMS is a multi-threaded parallel version of SCRUM for shared memory machines such as

those produced by SGI and Sun. Plans are to add an MPI option for distributed memory

systems. ROMS plans also include nesting, vertical splines, and data assimilation. For the

latest news on ROMS and the various support programs, please contact Hernan Arango

(arango@imcs.rutgers.edu).

Our group is continuing to explore new directions in ocean modeling. Our plans include:

• We have a 2-dimensional model which uses spectral finite elements in the horizontal [24].

The corresponding 3-dimensional model has been written and is being tested on more and

more realistic problems. Ask Mohamed Iskandarani (mohamed@imcs.rutgers.edu) for more

information.

We maintain an electronic mailing list of SeRUM users. If you would like to be added (or taken off)

the easiest thing to do is to send email to majordomo@imcs.rutgers.edu with a message containing

the word "help". Other possible messages include "subscribe serum" and "unsubscribe serum" ,

both with the obvious meanings. Hernan Arango, Dale Haidvogel and I use the list to send out

announcements of new model versions, patches to old versions, and news of SCRUM meetings.

There are even occasional announcements of job opportunities. The email list is a mail alias on

imcs.rutgers.edu such that all mail sent to scrum@imcs.rutgers.edu will go to everyone on the list.

We strongly recommend that you subscribe to this list if you use SCRUM (and learn to use patch!).

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

;.•

•

'.:,.".

:.•

•

•

•

•

•

'.•

:.•

i.

•

'.•

:.•

•

'.'.;.'
1.4 Warnings and bugs

SCRUM is not a large program by some standards, but it is still complex enough to require some

effort to use effectively. Section 6.1 attempts to describe what the user is responsible for-please

read it carefully.

More specific things to be wary of include:

• It is recommended that you use 64 bits of precision rather than 32 bits.

• The code must be run through the C preprocessor before it is compiled. This can occasionally

be dangerous, especially with the newer ANSI C versions of cpp. Potential problems are listed

in Appendix F. We have the source code for a version of cpp which is known to work at

ftp://ahab.rutgers.edu/pub/util/cpp.tar.gz.

• I have started declaring all variables as being of type BIGREAL and then defining them to

be real*8 (except on a Cray). This usually works for variables, but some compilers do not

like

real*8 function vmin(varl, var2)

If you run into one of these finicky compilers, you will just have to fix each JBIGREAL

function by hand.

• The SCRUM grid generation software was originally developed for us€~ with SPEM. Through

the years there have been several file formats for the SPEM grid file. Make sure that your

grid generation software is recent enough to create a NetCDF grid file as opposed to a binary

fort.3. The unformatted binary files created portability problems with different architectures

and with single vs. double precision.

• The vertical s coordinate was chosen as being a sensible way to handle variations in the water

depth. It has been used with success when the maximum and minimum depths differ by a

factor of twenty or less, and the value of the stretching parameter, THETA_5, is between

zero and five. It is also desirable to have the depth variations be well resolved by the horizontal

grid. For realistic problems we often fail to resolve the bathymetric slopes and WE~ then resort

to bathymetric smoothing. This in turn changes the shape of the basin and leads to its own

set of problems, such as altered sill depths. Also, the currents will react to the change in shelf

slope-you are now solving a different problem.

4

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
Equations (1) and (2) express the momentum balance in the x- and y-directions, respectively.

The time evolution of the potential temperature and salinity fields, T(x, y, z, t) and S(x, y, z, t), are

•

•

•

•

•

•

•

•

•

-

•

•

•

•

•

•

•

•

•

•

•

•

••

•

•

'.:.•

:.

,e

!:. '.

...

•

•

•

,:.'.:.;.'.'.'.'

2 Ocean Model Formulation

2.1 Equations of motion

The primitive equations in Cartesian coordinates can be written:

ov _ o¢

ot + v· \7v + fu:=: - oy + F v + V v

a;:; + v· \7T = Fr +Vr

~~ + v· \7S :=: Fs + V s

P = p(T,S,P)

o¢ -pg

OZ Po

OU ov. OW_ o

ox + oy + OZ - .

The variables are shown in Table 2.1.

Variable Description

Vu , Vv , Vr ,Vs diffusive terms

Fu , Fv , Fr, Fs 'forcing terms

f(x, y) Coriolis parameter

9 acceleration of gravity

h(x,y) bottom depth

V,'" horizontal viscosity and diffusivity

Km,Kr,Ks vertical viscosity and diffusivity

P total pressure P :::::: -Pogz

¢(x, y, z, t) dynamic pressure ¢ = (P/ Po)

Po + p(x, y, z, t) total in situ density

S(x, y, z, t) salinity

t time

T(x, y, z, t) potential temperature

u,v,w the (x, y, z) components of vector velocity v

x,y horizontal coordinates

z vertical coordinate

((x, y, t) the surface elevation

Table 1: The variables used in the description of the ocean model

5

(1)

(2)

(3)

(4)

(5)

(6)

(7)
6

Table 2: The variables used in the vertical boundary conditions for the ocean model

The vertical heat and salt flux may also be prescribed at the bottom, although they are usually set

to zero.

2.2 Vertical boundary conditions

The vertical boundary conditions can be prescribed as follows:

•

•

•

•

•

•

•

•

•

•

•

e

•

•

•

•

•

•

•

•

•

•

•

•

•

e.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

K m ~~ =T;(X, y, t)

K m ~~ = T%(X, y, t)

K aT - QT + 1 dQT (T T, ) . T oz - Poc? Poc? dT . - ref

K aS _ (E-P)S

S oz - Po

w- IK - at

K m ~~ .~ Tb(X, y, t)

K m ~~ = T~(X,y,t)

K oT_O T oz -

K as - 0 S oz -

-w +if· \7h = O.

top (z = ((x, y, t))

Variable Description

E-P evaporation minus precipitation

,1,,2 linear and quadratic bottom stress coefficients

QT surface heat flux x y T surface wind stress s , Ts x y T bottom stress b , Tb

Tref surface reference temperature

and bottom (z = -h(x, y))

Tb = (,on + ,2 Ju2 + v2 )u

r~ = (,1 +'Y2JU2 +v2 )v

The surface boundary condition variables are defined in Table 2.2. Since QT is a strong function

of the surface temperature, it is also prudent to include a correction term for the change in Q as

the surface temperature drifts away from the reference temperature that was used in computing

QT. On the variable bottom, z = -h(x, y), the horizontal velocity components are constrained to

accommodate a prescribed bottom stress which is a sum of linear and quadratic terms:

governed by the advective-diffusive equations (3) and (4). The equation ofstate is given by equation

(5). In the Boussinesq approximation, density variations are neglected in the momentum equations

except in their contribution to the buoyancy force in the vertical momentum equation (6). Under

the hydrostatic approximation, it is further assumed that the vertical pressure gradient balances

the buoyancy force. Lastly, equation (7) expresses the continuity equation for an incompressible

fluid. For the moment, the effects of forcing and dissipation will be represented by the schematic

terms :F and 'D, respectively. The horizontal and vertical mixing will be described more fully in

§3.1O.
where

and

i = t.

(8)

7

_. (~)fu -- - fv +v' \7u = --- - - - g- + Fu + 'Du at ax Po ax ax

az

Hz == as

As a trade-off for this geometric simplification, the dynamic equations become somewhat more

complicated. The resulting dynamic equations are, after dropping the carats:

See Appendix B for the form of s used here. In the stretched system, the vertical coordinate s

spans the range -1 ::; s ::; 0; we are therefore left with level upper (s = 0) and lower (s = -1)

bounding surfaces. The chain rules for this transformation are:

2.4 s (stretched vertical) coordinate system

From the point of view of the computational model, it is highly convenient to introduce a stretched

vertical coordinate system which essentially "flattens out" the variable bottom at z = -h(x, y).

Such "s" coordinate systems have long been used, with slight appropriate modification, in both

meteorology and oceanography (e.g., Phillips [43] and Freeman et al. [13]). To proceed, we make

the coordinate transformation:

y=y

s = s(x,y,z)

z = z(x,y,s)

x=x

2.3 Horizontal boundary conditions

As distributed, the model can easily qeconfigured for a periodic channel, a doubly periodic domain,

or a closed basin. Code is also included for open boundaries which mayor may not work for your

particular application. Appropriate boundary conditions are provided for u, v, T, 5, and (. At every

timestep the subroutines bcs2d and bcs3d are called to fill in the necessary boundary values.

The model domain is logically rectangular, but it is possible to mask out land areas on the'

boundary and in the interior. Boundary conditions on these masked regions are straightforward,

with a choice of no-slip or free-slip walls.

If biharmonic friction is used, a higher order boundary condition must also be provided. The

model currently has this built into the code where the biharmonic terms are calculated. The high

order boundary conditions used for u are tx (~~ ~) = 0 on the eastern and western boundaries

and ty (~~ ~) = 0 on the northern and southern boundaries. The boundary conditions for v, T,

and 5 are similar. These boundary conditions were chosen because they preserve the property of

no gain or loss of volume-integrated momentum, temperature, or salt.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

:.•

.'•

•

•

•

•

'.'.I:
Note the simplification of the boundary conditions on vertical velocity that arises from the s coordinate transformation.

8

if={u,v,n)

p = p(T,S,P)

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(9)

(14)

(13)

(12)

(11)

(10)

(lim.) &u = TX(X Y t) Hz &s s"

(lim.) &v =TY(X Y t) Hz &s S"

( KT) &T _ .9L + _l_dQ (T - T, ) Hz &s - PoCp PoCP dT ref

(

&) &5 _ (E-P)S

Hz &s - Po

n=o

(lim.) &u = TX(X Y t) Hz &s b"

( Km) &v = TY(X Y t) Hz &s b"

(!iT..) &T = 0 Hz &s

(!is.) &5 = 0

Hz &s

n=o.

top (s= 0)

if· V' = u~ + vi!- + n~. ox oy os

oHz + o(Hzu) + o(Hzv) + o(Hzn) = 0

ot ox oy os

a;:; + if ·V'T = FT + VT

~~ + if· V'S = Fs + Vs

o¢ =(-9HzP)

os Po

1 [ o( oz OZ] n(x y s t) = - w -(1 + s) - - u- - v- , , , Hz ot ox oy .

ov _. o¢ (gp) OZ 8( - +

m·fu + v· V'v = --- - - - g- + Fv + Vv

and bottom (s = -1)

OZ oz. oz

w = ot + u ox + v oy + nHz.

In the stretched coordinate system, the vertical boundary conditions become:

and

The vertical velocity in s coordinates is

where
2.5 Horizontal curvilinear coordinates

In many applications of interest (e.g., flow adjacent to a coastal boundary), the fluid may be confined

horizontally within an irregular region. In such problems, a horizontal coordinate system which

conforms to the irregular lateral boundaries is advantageous. It is often also true in many geophysical problems that the simulated flow fields have regions of enhanced structure (e.g., boundary

currents or fronts) which occupy a relatively small fraction of the physical/computational domain.

In these problems, added efficiency can be gained by placing more computational resolution in such

regions.

The requirement for a boundary-following coordinate system and for a laterally variable grid

resolution can both be met, for suitably smooth domains, by introducing an appropriate orthogonal

coordinate transformation in the horizontal. Let the new coordinates be ~(x, y) and 77(X, y), where

the relationship of horizontal arc length to the differential distance is given by:

(ds)~ = (~) d~ (15)

(ds)7J = (~) d77 (16)

Here, m(~, 77) and n(~, 77) are the scale factors which relate the differential distances (D.~, D.77) to

the actual (physical) arc lengths. Appendix C contains the curvilinear version of several common

vector quantities.

Denoting the velocity components in the new coordinate system by

(22)

(23)

(21)

(20)

(19)

(18)

(17)

(HzT) + (HzuT) + (HzVT) + (HzOT) = Hz (FT + 'Dr) ot mn o~ n 077 m os mn mn

9

P = p(T,S,P)

(HzS) + (HzuS) + (HzvS) + (HzOS) = Hz (Fs +'Ds) ot mn o~ n 077 m os mn mn

~(HzV)+~(HzUV)+~(HzV2)+~(HzVO) ot mn o)~ n 077 m os mn

+ {(~n) +v :~ (~) - u:77 (~)} Hzu =

_ (Hz) (01 +9Poz +90() + Hz (Fv+'Dv)

m 077 Po 077 077 mn

iJ·r,=v

the equations of motion (8)-(14) can be re-written (see, e.g., Arakawa and Lamb [2]) as:

and

'.'.:.•

•

•

'.'.•

•

•

•

•

•

•

•

•

•

I.•

•

•

•r.i.'.•

•

•

•

:.'.:.

••,,'.•

•

•

•

•

'.'.•

•
(}[z

)+ (HzU)+ (Hzv)+ (HzO) =: O. at mn a~ n aT} m as mn

Sincez isa linear function of (, equation (25) can be rewritten as:

~(--L)+ ~(HzU)+ (Hzv)+ (HzO) =: O. at mn a~ n aT} m as mn

All boundary conditions remain unchanged.

10

(24)

(25)

(26)

••••••••••••••••••••••••••••••••••••••••••••
'.•

•

'.•

•

•

•

•

•

,e•

•

•

•

•

•

•

I.•

ie'.•

'.,'.;.'.'I.".•

•

•

•

I.•

•

•

•

'.•

•

•

•

•

3 Numerical Solution Technique

3.1 Vertical and horizontal discretization

3.1.1 Horizontal grid

In the horizontal (~, "1), a traditional, centered, second-order finite-difference approximation is

adopted. In particular, the horizontal arrangement of variables is as shown in Fig. 2. This is

equivalent to the well known Arakawa "C" grid, which is well suited for problems with horizontal

resolution that is fine compared to the first radius of deformation (Arakawa and Lamb [2]).

.6.~

i

Vi,j+1 I

U' , (p, h,i, 0kj Ui+I,j t,)

- - 0 - - .6."1

+

v' ,I t,)

Figure 2: Placement of variables on an Arakawa C grid

3.1.2 Vertical grid

The vertical discretization also uses a second-order finite~difference approximation. Just as we use

a staggered horizontal grid, the model was found to be more well-behaved with a staggered vertical

grid. The vertical grid is shown in Fig. 3.

-----WN

• PN -- WNm

• PNm

•

•

W2

• P2 --- WI

• PI -----Wo

Figure 3: Placement of variables on staggered vertical grid

11
12

3.2.1 Velocity

Figure 4: Masked region within the domain

3.2 Masking of land areas

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

x .- u points

0- p points

0- v points

• -- 'l/J points

A B C

0 0 0 0

D E F - .~ ,

G H I J

,0 0 x 0 0

K L - - - -

M N

0 0 0 0

- - - -

At the end of every timestep, the values of many variables within the masked region are set to zero

by multiplying, by the mask for either the u, v or p points. This is appropriate for the v points E

and L in Fig. 4, since the flow in and out of the land should be zero. It is likewise appropriate for

the u point at I, but is not necessarily correct for point G. The only term in the u equation that

requires the u value at point G is the horizontal viscosity, which has a term of the form :T/v~~.

Since point G is used in this term by both points A and M,it is not sufficient to replace its value

with that of the image point for A. Instead, the term ~~ is computed and the values at points D

and K are replaced with the values appropriate for either free-slip or no-slip boundary conditions.

Likewise, the term :f,v~~ in the v equation must be corrected at the mask boundaries.

_ This is accomplished by having a fourth mask array defined at the 'l/J points, in which the values

depend on whether or not free-slip boundaries are desired. In the case of free-slip, the value of ~~

is simply set to zero at points D and K. For no-slip boundaries, we count on the values inside the

land (point G) having been zeroed out. For point 0, the image point at G should contain minus

the value of u at point A. The desired value of ~~ is therefore 2UA while instead we have simply

UA. In order to achieve the correct result, we multiply by a mask which contains the value 2 at

point D. It also contains a 2 at point K so that ~~ there will acquire the desired value of -2UM.

The corner point F is treated in the same way as points D and K. This can be changed in

geLmask if desired.

seRUM has the ability to work with interior land areas, although the computations occur over

the entire model domain. One grid cell is shown in Fig. 2 while several cells are shown in Fig. 4,

including two land cells. The process of defining which areas are to be masked is described in §5.2,

while this section describes how the masking ,affects the computation of the various terms in the

equations of motion.
(32)

(31)

(30)

(29)

(28)

(27)

p = p(S,T,P)

(hS) J: (uH;~S~) . J: (vH:'S'I) 0 (-S8 HZ

n ) = 1) F

at + U _c + UTJ ----n + 8 S + s

mn n" m'/ mn

13

Here o~, oTJ and 08 denote simple centered finite-difference approximations to a/oe a/ory and a/as

with the differences taken over the distances b.~, b.ry and b.s, respectively. -(-)~, r--T and -(-)8

3.3 Conservation properties

seRUM conserves the first moments of u, v, S, and T. This is accomplished by using the flux form

of the momentum and tracer equations. It is also necessary to be careful when averaging between

the velocity and tracer grids, for instance to obtain u at p points. The semi-discrete form of the

dynamic equations (19)-(22) is:

3.2.2 Temperature, salinity and surface elevation

The handling of masks by the temperature, salinity and surface elevation equations is similar to

that in the momentum equations, and is in fact simpler. Values of T, Sand ( inside the land

masks, such as point H in Fig. 4, are set to zero after every timestep. This point would be used by

the horizontal diffusion term for points B, J, and N. This is handled by setting the first derivative

terms at points E, I, and L to zero, to be consistent with a no-flux boundary condition. Note

that the equation of state must be able to handle T = S = 0 since this is the value inside masked

regions.
14

where Rrj> represents all of the forcing terms other than the vertical viscosity or diffusion. Since we

want the diffusion term to be evaluated partly at the current tirnestep n and partly at the next

timestep n + 1, we introduce the parameter>. and rewrite equation (34) as:

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(33)

(34)

(35)

(37)

(36)

a ( K a¢) - ----- as Hzmn as

. a(Hz¢) _ R (~ a¢)

at - mn rj> + as Hz as

a(Hz¢) =mnR + (1 _ >.)~ (K a¢n) + >.~ (~a¢n+l) . at rj> as Hz as as Hz as

The discrete form of equation (35) is:

where k is used as the vertical level index. This can be reorganized so that all the terms involving

¢n+1 are on the left and all the other terms are on the right. The equation for ¢~+1 will contain

terms involving the neighbors above and below (¢ni and ¢~+i) which leads to a set of coupled

equations with boundary conditions for the top and bottom. The general form of these equations

is:

where ¢ represents one of u, v, T or 5, and K is the corresponding vertical viscous or diffusive

coefficient. This is timestepped using a semi-implicit Crank-Nicholson scheme with a weighting of

0.5 on the old timestep and 0.5 on the new timestep. Specifically, the equation of motion for ¢ can

be written as:

represent averages taken over the distances~~, ~'T} and ~s. 12 indicates a second-order vertical

integral computed as a sum from level s to the surface at s = O.

This method of averaging was chosen because it internally conserves first moments in the model

domain, although it is still possible to exchange mass and energy through the open boundaries.

The method is similar to that used in Arakawa and Lamb [2]; however, their scheme also conserves

enstrophy.

The continuity equation will be discussed below in §3.8.

3.4 Vertical viscosity and diffusion

The 'Du , 'Dv , 'DT and 'Ds terms in equations (27)-(30) represent both horizontal and vertical mixing

processes. The horizontal options will be covered in §3.10. The model has several options for

computing the vertical coefficients; these will be described in §3.11. The vertical viscosity and

diffusion terms have the form:
15

where the boundary conditions are written into the coefficients for the end points. In this case the

coefficients become:

This is a standard tridiagonal system for which the solution procedure can be found in any standard

reference, such as Press et al. [44].

(50)

(49)

(48)

(43)

(44)

(45)

(46)

(47)

(42)

(41)

(40)

(38)

(39)

_ 1 fO A = D Hz Ads

-1

D = ((~, '1], t) + h(~, 'f/)

D(N) =

A(I)

A(2 : N)

B(I) =

B(2: Nm) -

B(N) =

C(1 : Nm) =

C(N)

D(I)

D(2: Nm)

(Dfi) + (DUU) + (DUV) _ Dfv ot mn o~ n O'f/ m mn

- [VV :~ (~) - UV:'f/ (~)] D =- (O~2 + 9~~ )

+ D (Fu +'Ph,,) + _1_ (r; -r;) (51)

mn mn

where (/>2 includes the ~~ term, 'Phu is the horizontal viscosity and the vertical viscosity only

contributes through the upper and lower boundary conditions. The corresponding vertical integral

is the total depth of the water column. The vertical integral of equation (19) is:

where the overbar indicates a vertically averaged quantity and

3.5 Depth-integrated equations

The depth average of a quantity A is given by:

•

•

•

.'•

•

•

•

•

•

••

•

•

•

•

•

•

•

•Ie

•

•

•

•

•

•

•

i.•

'.I.•

:.•

'.'.•

:.•

•

•

:.'.

•
16

Figure 5: The split timestepping used in the model.

We also need the vertical integral of equation (26). Using the vertical boundary conditions on n

we get:

When timestepping the model, we compute the right-hand-sides for equations (27) and (28) as

well as the right-hand-sides for equations (54) and (55). The vertical integral of tlle 3-D righthand-sides are obtained and then the 2-D right-hand-sides are subtracted. The resulting fields are

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

-I-I•

•

•

•

•

•

-

•

•

•

•

•

•

e

•

•

•

•

•

•

•

•

•

•

•

•

•

(53)

(52)

time

(i-) +~ (DU) +~ (D15)_O at mn a~ n aT) m -.

dtfast

f--dt--H (

n-1 n n+1 n+2

I I I I I I I I I I I I II I 1--i.1--l.1...;.L.....L-LI...L1-1-1..1-1.1-1.1--1~IIL>-Io'--i.I-i.d._--..

'--______ _ __-1

V

(avg

- - a (D15) +-a (DUV) -- +-\a (DVV) -- +-- Dfu at mn a~ n aT] m mn

+[uv~ (~) -uu~ (~)] D= _ D(a(!J2 +ga() a~ n aT) m m aT) aT)

+ -D (- ~) 1 ( TI) mn F v +1)h + - TTl -- T . v mn S b

of equation (20) is:

- - a (D15) +-a (Du15) -- +-a (D1515) -- +-- Dfu at mn a~ n aT] m mn

+[U15~ (~) - UU~ (~)] D= Rv - gD a( + D1)1j _ ._1TTl. (55) a~ n aT] m slow m aT] mn mn b

~(DU)+~(DUU)+~(DU15)_Df15 at mn a~ n aT] m mn

- [1515~ (~) - u15~ (~)] D= Ru - gD a( + D1)u; _ _ l_T~ (54) a~ n aT) m '- slow n a~ mn mn b

Some of the terms in equations (51) and (52) are updated on the short timestep while others

are not. The contributions from the slow terms are computed once per long timestep and stored.

If we call these terms RUSlow and Rvs1ow ' equations (51) and (52) become:

. The presence of a free surface introduces waves which propagate at a speed ofJ9h. These

waves usually impose a more severe timestep limit than any of the internal processes. We have

therefore chosen to solve the full equations by means of a split timestep. In other words, the depth

integrated equations (51), (52), and (53) are integrated using a short timestep and the values of

U and v are used to replace those found by integrating the full equations on a longer timestep. A .

diagram of the timestepping is shown in Fig. 5.
3.7.1 Smolarkiewicz

3.7 Advection schemes

Rather than computing the centered average for these terms, we use the "upwind" or "donor cell"

value for 'lj;. For instance:

Thus far, the advection scheme presented here is a centered second-order scheme. This scheme

is known to have some unfortunate properties in the presence of strong gradients, such as large

over- and under-shoots of tracers, leading to the need for large amounts of horizontal smoothing.

SeRUM also provides two alternative advection schemes with better behavior in many situations.

At present, the alternatives are only implemented in the full 3-D engine of the model.

(59)

(58)

(60)

(57)

(56)

17

-~

Fi~j,k = ~~ [max(O, Ui,j,k)'lj;i-l,j,k + min(O, Ui,j,k)'lj;i,j,k] .

Hz'lj; = -~F~ _ ~FTJ _ ~FS. at mn a~ a", as

where we have introduced the advective fluxes:

F~ = Hzu'lj;

n

FTJ = Hzv'lj;

m

F S = HzfJ'lj;.

mn

3.6 Time stepping: internal velocity modes, temperature, and salinity

The momentum equations (27) and(28) are advanced by computing all the terms except the vertical

viscosity and then using the implicit scheme described in §3.4 to find the new values for u and v.

The depth-averaged component is then removed and replaced by the IT and v computed as in §3.5.

A third-order Adams-Bashforth timestepping is used when computing the right-hand-side terms

(see Appendix A). The temperature and salinity equations (29) and (30) are also advanced as in

§3.4. There is also an option to advect the temperature and salinity using a Smolarkiewicz scheme;

this affects the time stepping as described in §3.7.1.

the slow forcings RUslOW and RVslow' .This was found to be the easiest way to retain the baroclinic

contributions of the non-linear terms such as uu - uu.

The model is timestepped from time n to time n +1 by using short timesteps on equations (54),

(55) and (53). Equation (53) is timestepped first, so that an estimate of the new D is available

for the time rate of change terms in equations (54) and (55). A third-order predictor-corrector

timestepping is used. In practice, we actually timestep all the way to time (n + 2) - dtfast and

then average the values of u, v and (. The averages are used to replace the values at time n + 1.

These time averages damp out certain instabilities which would otherwise grow to dominate the

solution.

Smolarkiewicz has written a series of papers ([50], [51], [52] and [53]) on an advection scheme known

'as MPDATA (Multidimensional Positive Definite Advection Transport Algorithm). MPDATA is

meant to be applied to positive-definite tracer fields only. It uses a first-order upwind scheme to

obtain a first estimate of the updated field, followed by one or more "antidiffusion" passes to obtain

the final updated field. Note that we add an offset to the temperature field before advecting with

this scheme, since ocean temperatures can be negative (or we could use the Kelvin scale instead).

Recall that the advection of a tracer 'lj; has an equation of the form

•

•

•

•

•

•

•

•

•

•

•

,e

•

~.

Ie:.'.I.''.•

•

•,'.,.

•

•

'.•',.•

•

'.'.'.,:.•

'.'.'.•

•

'.•
18

Then, we define the l3-ratios:

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

(64)

(63)

(62)

(61 )

(65)

(I3t)-l ('l/J0~k - 'l/Ji,j,k) = mt::.t [max(O,Ui,j,k)'l/Ji_l,j,k - min(O,ui+l,j,k)'l/Ji+l,j,k]

+ n t::.t [max(O, Vi,j,k)'l/Ji,/-l,k - min(O, Vi,j+l,k)'l/Ji,j+l,k]

+ ~: [max(O, Wi,j,k-d'I/Ji,j,k-l - min(O, Wi,j,k)'l/Ji,j,k+d

1)-1 (nl.*. _ nl.min ) = mt::.t [max(O U· . )nl.,:. - min(O U·· )nl.*. ] f.J... 'l-'z,J,k 'l-'z,J,k . , z+l,J,k 'l-'z,J,k , z,J,k 'l-'z,J,k

+ n t::.t [max(O, Vi,j+l,k)'l/Ji,j,k - min(O, Vi,j,k)'l/Ji,j,k]

+ ~: [max(O, Wi,j,k)'l/Ji,j,k - min(O, Wi,j,k-l)'l/Ji,j,k]

This picks up the value of'I/J to the left if u is positive, otherwise it picks up the value to the right

(Ff, is located at a u point on the grid).

Using these fluxes, we compute an upwind estimate of the new tracer:

'l/Ji:j,k = min('l/Ji,j,k, 'l/Ji-l,j,k, 'l/Ji+l,j,k, 'l/Ji,j-l,k1 'l/Ji,j+l,k, 'l/Ji,j,k-l, ~r/Ji,j,k+l,

'l/Ji,j,k' 'l/Ji-l,j,k' 'l/Ji+l,j,k' 'l/Ji,j-l,ki 'l/Ji,j+l,kl 'l/Ji,j,k-l' Ir/Ji,j,k+l)

'l/J0~k = max('l/Ji,j,k, 'l/Ji-l,j,k, 'l/Ji+l,j,k, 'l/Ji,j-l,k, 'l/Ji,j+l,k, 'l/Ji,j,k-l, Ir/Ji,j,k+l,

'l/Ji,j,kl 'l/Ji-l,j,k, 'l/Ji+l,j,kl 'l/Ji,j~ l,k' 'l/Ji,j+l,k' 'l/Ji,j,k-l, Ir/Ji,j,k+l)

This estimate is used to create an "antidiffusive" velocity:

In fact, this correction step can be repeated any number of times in order to improve the solution.

Thus far, we have described the original MPDATA scheme which maintains the positive. definite

nature of the tracer. However, it does not prevent under- and overshoots in a tracer such as salinity,

in which the initial range would be say 34 to 36 PSU. An optional nonoscillatory modification to

MPDATA is described in Smolarkiewicz [53]. It uses ideas from the FCT (flux-corrected transport)

scheme to limit the antidiffusive fluxes. This modification is based on multiplying the antidiffusive

fluxes by a coefficient C such that 0 :S C :S 1 and the resulting tracer value is bounded by

'l/Jmin :S 'l/Jn+ 1 :S 'l/Jmax.

First, we need to define 'l/Jmin and 'l/Jmax .

In SCRUM, these velocities are set to zero if the first term is smaller than the second term (there is

no such check in the ice MPDATA code, nor is there mention of such a thing in the Smolarkiewicz

papers). Using these velocities in another upwind iteration leads to a second estimatE~ of the new

tracer fields:
s Hzw [ 1 9 9 1 ] F = mn -16'l/Ji,j,k-l + 16'l/Ji,j,k + 16'l/Ji,j,k+l - 16 'l/Ji,j,k+2 (71)

One advantage of UTOPIA over MPDATA is that it can be used on variables having both

negative and positive values. Therefore, it can be used on velocity as well as scalars (is there a

reference for this?).' For the u-velocity, we have:

Finally, the flux-limited form of the antidiffusion velocities is:

u' = min(l, f3 4.i-l,j,k, f3ti,j,k) max(O, u) + min(l, f3ti-l,j,k, f3 4.i,j,k) min(O, u) (66)

v' = min(l, f3 4.i,j-l,k, f3ti,j,k) max(O, v) + min( 1, f3ti,j-l,k> f3 .!.i,j,k) min(O, v) (67)

w' = min(l, f3 4.i,j,k, f3ti,j,k+d max(O, w) + min(l, 13ti,j,k' f3 4.i,j,k+l) min(O, w) (68)

These velocities are used in equation (65) to produce the flux-limited estimate of'l/Jn+l.

The second derivative terms are centered on a p point in the grid, but are needed at a u or v point

in the flux. The upstream value is used [see equation (60)]. The value of, in the model is kwhile

that in Rasch [45] is i.

Because the third-order upwind scheme is designed to be two-dimensional, it is not used in the

vertical (though one might argue that we are simply performing one-dimensional operations here).

Instead, we use a centered fourth-order scheme in the vertical when the third-order upwind option

is turned on:

(77)

(76)

(75)

(74)

(73)

(72)

(69)

(70)

19

F~ = (u _ , 8

2

U) [Hzu _ ,~ (Hzu)] 8e n 8e n

P'/ = (u _ ,8

2

U) [Hzv _,~ (HzV)] 8~2 m 8~2 m

s Hzw [ 1 9 9 1 ] F = mn -16Ui ,j,k-l + 16Ui,j,k + 16Ui,j,k+l - 16Ui,j,k+2

while for the v-velocity we have:

F~ = (v _,8

2

V) [HzU _,~ (HtU)] 8~2 n 8~2 n

FT/ = (v _ ,8

2

V) [Hzv _,~ (Hzv)] 8~2 m 8~2 m

s Hzw [ 1 9 9 1 ] F = mn -16Vi,j,k-l + 16Vi,j,k + 16Vi,j,k+l - 16Vi,j,k+2

In all these terms, the second derivatives are evaluated at an upstream location.

3.7.2 Third-order Upwind

There is a class of third-order upwind advection schemes, both one-dimensional (Leonard [30])

and two-dimensional (Rasch [45]). This scheme is known as UTOPIA (Uniformly Third-Order

Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in Thuburn

[58], although it is not implemented in SCRUM., The two-dimensional formulation in Rasch contains

terms of order u2'l/J and u3'l/J, including cross terms (uv'l/J). The terms which are nonlinear in velocity

have been dropped in SCRUM, leaving one extra upwind term in the computation of the advective

fluxes:

•

•

•

•

•

•

•

•

•

•

e

,e

•

'.'.'.;.'.'.•

•

'.•

•

•

•

'.•

•

•

".'.•

•

,:.•

'.".•

•

e

•

•

•
20

3.10 Horizontal friction and diffusion

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(79)

(80)

HzO =! [6 (U!!~) + 6 (15YY') _6(u:!z() _6 (VH';7)] ds. mn n~ TJ ffi'1 n~ TJ mTJ

Solving for HzO/mn and using the semi-discrete notation of §3.3 we obtain:

This is the form traditionally used in sigma-coordinate models to account for the horizontal differences being taken along surfaces of constant s. This form can be shown to lead to significant errors

when IV'hl is large (Haney [20J; and Beckmann and Haidvogel [5]).

The pressure ¢ is computed by a vertical integration of the density field using equation (24).

Prior to the integration, a horizontal average ofthe density, p(z), is subtracted from p. As discussed

by Haney [20J and McCalpin [33], p does not contribute to the pressure gradient in the analytic

equations. However, when numerically computing its contribution to the pressure gradient, the

error from this term can be unacceptably large.

3.9 The pressure gradient terms

The pressure gradient terms in equations (19) and (20) are written in the form

3.8 Determination of the vertical velocity and density fields

Having obtained a complete specification of the u, v, T, and 5 fields at the next time level by the

methods outlined above, the vertical velocity and density fields can be calculated. The vertical

velocity is obtained by combining equations (26) and (53) to obtain:

In Chapter 2, the diffusive terms were written simply as 'Du , 'Dv , 'DT, and 'Ds. The vertical component of these terms was described in §3.4. Here we describe SCRUM's options for representing the

horizontal component of these terms.

The integral is actually computed as a sum from the bottom upwards and also as a sum from the

top downwards. The value used is a linear combination of the two, weighted so that the surface

down value is used near the surface while the other is used near the bottom.

The density is obtained from temperature and salinity via an equation of state. SCRUM

provides a choice of a nonlinear equation of state p = p(T, 5, z) or a linear equation of state

p = p(T). The nonlinear equation of state has been modified and now corresponds to the UNESCO

equation ofstate as derived by Jackett and McDougall [25J. It computes in situ density as a function

of potential temperature, salinity and pressure.

Warning: although we have used it quite extensively, McDougall (personal communication)

claims that the single-variable (p = p(T)) equation of state is not dynamically appropriate as is.

He has worked out the extra source and sink terms required, arising from vertical motions and the

compressibility of water. They are quite complicated and we have not impllemented them to see if

they alter the flow.
'1"'-'-

and is

(84)

(83)

(82)

(81)

where 'I/J is once again any of u, v, T, and S. Note that u and v are treated as independent scalar

quantities rather than as a vector. The complete Laplacian operator on a vector quantity ii contains

additional terms, including v terms in the u equation and vice versa. These extra terms were found

to be small in a test problem and have been left out of the model.

3.10.3 Rotated mixing tensors

Both the Laplacian and biharmonic terms above operate on surfaces of constant s and can contribute

substantially to the vertical mixing. However, the oceans are thought to mix along constant density

surfaces so this is not entirely satisfactory. Therefore, the option of using rotated mixing tensors

for the Laplacian and biharmonic operators has been added. Options exist to diffuse on constant

z surfaces (MIX_GP_UV, MIX_GP_TS), constant in situ density surfaces (MIX_EPLUV,

MIX-EPLTS), and constant potential density surfaces ((MIX.-lSO_UV, MIXJSO_TS)

The horizontal Laplacian diffusion operator is computed by finding the three components of

the flux of the quantity 'I/J. The and 7] components are locally horizontal, rather than along the

s surface. Here, MIX-EPI represents both the epineutral and isopycnal options. The diffusive

21

The biharmonic operator is \74 = \72\72; the corresponding term is computed using a temporary

variable Y:

3.10.2 Biharmonic

where 'I/J is any of u, v, T, and S. This form guarantees that the term does not contribute to the

volume-integrated equations, except when using no-slip boundaries in the momentum equations.

2: [8 (m 8'I/J) a(n 8'I/J)] \7 'I/J = \7 . \7'I/J = mn 8~ 8~ + 87] m 87]

This term in SCRUM is multiplied by 1I;n~z and becomes

3.10.1 Laplacian

The Laplacian of a scalar 'I/J in curvqinear coordinates is (see Appendix C): ,

•

•

•

•

•

•

•

•

'.,e:.'.I.•

'.:.•

•

•

•

•

•

•

•

:.•

•

•

'.•

•

•

•

•

'.;.•

•

•

".•

•

•

•
22

Finally, the flux divergence is calculated and is added to the right-hand-side term for the field

being computed:

•

•

•

•

e 

•

•

•

•

•

•

•

•

•

•

•

•

e

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(87)

(86)

(85)

(88)

at walls

at T] walls

at s = -1,0

V'

MILGP

(

az ) a'ljJ ma~' as

MIX...EPI

Ff. =0

P" == 0

F S = 0

a'ljJ m--

a~

a'ljJ n--

aT] [

az ) a'ljJ

naT] as

MIX...EPI , #

....

MILGP

FS= - -1(m-az +Sx ) Ff. - -1.(n-az +Sy ) P" Hz a~ "--" Hz aT] "--"

MIX...EPI MIX...EPI , ", , .#

MIX_GP Mli~GP

~(~Ff.)+~(~F'I)+~(~FS) a~ n aT] m as mn

The biharmonic rotated mixing tensors are computed much as the non-rotated biharmonic

mixing. We define a temporary variable Y based on equation (88):

No flux boundary conditions are easily imposed by setting

Y = mn .[~ (HzFf.) + (HzF'I) + (HzFS)J . (89) Hz a~ n aT] m as mn .

We then build up fluxes of Y as in equations (85)-(87). We then apply equation (88) to these Y

fluxes to obtain the biharmonic mixing tensors.

3.11 Vertical mixing schemes

SeRUM contains a variety of methods for setting the vertical viscous and diffusive coefficients.

The choices range from simply choosing fixed values to the KPP and Mellor-Yamada. turbulence

closure schemes. See Large [27] for a review of surface ocean mixing schemes. Many schemes have

a background molecular value which is used when the turbulent processes are assumed to be small

(such as in the interior).

where

fluxes are:
23

Pacanowski and Philander [40] developed a vertical mixing parameterization based on measurements in the equatorial oceans:

(97)

(98)

(96)

(95)

(93)

(94)

(91 )

(90)

(92)

= 10-7 1.0, C ,

if Rig < 0,

if Rig = 0,

if Rig> O.

Pb = KsN 2 ,

q3

~d = Bil

P, K m [(~)' + (:)'] ,

K - va K m - (1 + aRig)n + mbackground

K - K m K

s - (1 + aRig) + Sbackground

{

Co

K S = K Sbackground

, min(vmax , max(Vmin, C/ .j1'iI;))

K m -K- mbackground

where Ps is the shear production, Pb is the buoyant production and ~d is the dissipation of turbulent

kinetic energy. These terms are given by

3.11.3 Mellor-Yamada

One of the more popular closure schemes is that of Mellor and Yamada [36], [37]. They actually

present a hierarchy of closures of increasing complexity. SeRUM provides only the "Level 2.5"

closure with the Galperin et al. [14] modifications as described in Allen et al. [1]. This closure

scheme adds two prognostic equations, one for the turbulent kinetic energy (~q2) and one for the

turbulent kinetic energy times a length scale (q2l).

The turbulent kinetic energy eq\Jation is:

where B1 is a constant. One can also add a traditional horizontal Laplacian or biharmonic diffusion

(Dq ) to the turbulent kinetic energy equation. The form of this equation in the model coordinates

The constants here are va == .01, n == 2 and a = 5. The values of K m and Ks can get very large for

negative values of Rig, so we choose to limit these values to that obtained for Rig = O.

3.11.2 Pacanowski-Philander

Where Ks applies to any scalar. The constants in this expression are Co =

Vmin = 3 X 10-5 and Vmax = 4 X 10-4 , and the gradient Richardson number is

3.11.1 Brunt-Vaisala frequency scheme

One of the simplest schemes is to set the: vertical diffusion coefficients to be large when the water

column is vertically unstable. The vertical viscosity is uniform:

•

•

•

•

•

•

•

•

•

•

•

~.'.•

•

•

•

'.•

'.•

•

•

•

•

•

'.:,.•

•

•

••

•

•

•

•

•

•

•

•

•

•

•

•
The equation is timestepped much like the model tracer equations, including an implicit solve for

the vertical operations and an option for using the third-order upwind advection.

There is also an equation for the turbulent length scale l:

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(104)

(105)

(106)

(102)

(101)

(100)

H.fl = 0 mn

aq2 = !!..L.-' 2/3 '[()2T~ + (1"1/)2] mnH" as Po s s

H K (aU aV) = l (T~ T1/) z m az' az Po s, s

H K N 2 =-.!L z s PoCp

H;fl = 0 mn

-!iL~ = B~/3 :.[(T~)2 -I- (T1/)2] mnH; as Po' b b

H K ( au aV) _ 1 (~ 1/) z m az' az .- Po Tb , Tb

H zK sN 2 = 0

24

_ ( l )2 W = 1 +E2 kL

L-1 = _1_. + __1__

(-z H+z

K m = qlSm + KmbaCkgrOund

K s = qlSh + KSbaCkground

K q = qlSq + KQbaCkgrOund

top (z = ((x, y, t))

and bottom (z = -h(x, y))

where'Dql is the horizontal diffusion of the quantity q2l.

Given these solutions for q and l, the vertical viscosity and diffusivity coefficients are:

The form of this equation in the model coordinates becomes

where W is the wall proximity function:

becomes

The vertical boundary conditions are:
25

and the stability coefficients Sm, Sh and Sq are found by solving

(114)

(113)

(112)

(111)

(where /'1, = 0.4 is von Karman's contant and Bf is the surface buoyancy flux), and the shallowest

depth at which a critical bulk Richardson number is reached. The critical bulk Richardson number

(Ric) is typically in the range 0.25-0.5. The bulk Richardson number (Rib) is calculated as:

Surface Boundary layer depth The boundary layer depth hsbl is calculated as the minimum

of the Ekman depth, estimated as,'

Ss [1 (3A2B2 + 18AI A2)Gh] = A2 [1 - 6AI B1

I

] (107)

Sm [1 - 9AIA2Gh] - Ss [Gh(18Ai + 9AIA2)Gh] = Al [1 - 3CI - 6AI B1

I

] (108)

Z2N2

Gh = min(--2-,0.028). (109)

q

Sq = 0.41Sm (110)

where a is a non-dimensional coordinate ranging from°to 1 indicating depth within the surface

boundary layer. The x subscript stands for one of momentum, temperature and salinity.

where d is distance down from the surface, B is the buoyancy, Br is the buoyancy at a near surface

reference depth, V is the mean horizontal velocity, Vr the velocity at the near surface reference

depth and Vi is an estimate of the turbulent velocity contribution to velocity shear.

(where u*is the friction velocity u* = JT'1 + TJ/p ), the Monin-Obukhov depth:

Surface boundary layer The Large, McWilliams and Doney scheme (LMD) matches separate

parameterizations for vertical mixing of the surface boundary layer and the ocean interior. A

formulation based on boundary layer similarity theory is applied in the water column above a

calculated boundary layer depth hsbl. This parameterization is then matched at the interior with

schemes to account for local shear, internal wave and double diffusive mixing effects.

Viscosity and diffusivities at model levels above a calculated surface boundary layer depth

(hsbl ) are expressed as the product of the length scale hsbl, a turbulent velocity scale W x and a

non-dimensional shape function.

3.11.4 The Large, McWilliams and Doney parameterization

The vertical mixing parameterization introduced by Large, McWilliams and Doney [28] is a versatile

first order scheme which has been shown to perform well in open ocean settings. Its design facilitates

experimentation with additional or modified representations of specific turbulent processes.

The constants are set to (AI,A2,BI,B2,CI,EI,E2) == (0.92,0.74,16.6,10.1,0.08,1.8,1.33). The

quantities q2 and q2Z are both constrained to be no smaller than 10-8 while Z is set to be no larger

than 0.53q/N.

•

•

•

•

•

•

•

•

•

•

'.

Ie'.'.•

•

•

•

•

'.•

•

•

•

'.•

'.•

•

•

•

'.•

I,.•

•

••

•

'.•

•

•

;.•
The turbulent velocity shear term in this equation is given by LMD as,

In near-neutral unstable conditions common Businger-Dyer forms are used which mat.ch wit.h the

formulation for stable conditions at ( = O. Near neutral conditions are defined as

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(122)

(123)

(120)

(121)

(119)

(117)

(118)

(116)

(115)

cPx = 1 + 5(

-1.0 :::; ( < 0

-0.2:::; ( < 0

cPm == (1 - 16()1/4

cPs = (1 - 16()1/2

26

cPm = (1.26 - 8.38()1/3

cPs == ('-28.86 - 98.96()1/3

for scalars. The non dimensional flux profiles in this regime are,

for momentum and,

where ( is the surface layer stability parameter defined as z/L. cPx is a non-dimensional flux profile

which varies based on the stability of the boundary layer forcing. The stability parameter used in

this equation is assumed to vary over the entire depth of the boundary layer in stable and neutral

conditions. In unstable conditions it is assumed only to vary through the surface layer which is

defined as Ehsbl (where E is set at 0.10) . Beyond this depth ( is set equal to its value at Ehsbl ·

The flux profiles are expressed as analytical fits to atmospheric surface boundary layer data. In

stable conditions they vary linearly with the stability parameter ( as

turbulent velocity scale To estimate W x (where x is m - momentum or s - any scalar)

throughout the boundary layer, surface layer similarity theory is utilized. Following an argument

by Troen and Mahrt [59], Large et al. estimate the velocity scale as

In more unstable conditions cPx is chosen to match t.he Businger-Dyer forms and with the free

convective limit. Here the flux profiles are

where Cv is the ratio of interior N to N at the entrainment depth, f3T is ratio of entrainment flux

to surface buoyancy flux, Cs and E are constants, and W s is the turbulent velocity scale for scalars.

LMD derive (115) based on the expect~d behavior in the pure convective limit. The empirical rule

of convection states that the ratio of the surface buoyancy flux to that at the entrainment depth be

a constant. Thus the entrainment flux at the bottom of the boundary layer under such conditions

should be independent of the stratification at that depth. Without a turbulent shear term in the

denominator of the bulk Richardson number calculation, the estimated boundary layer depth is too

shallow and the diffusivity at the entrainment depth is too low to obtain the necessary entrainment

flux. Thus by adding a turbulent shear term proportional to the stratification in the denominator,

the calculated boundary layer depth will be deeper and will lead to a high enough diffusivity to

satisfy the empirical rule of convection.
Shear generated mixing The shear mixing term is calculated using a gradient Richardson

number formulation, with viscosity estimated as: .

with the coefficients specified to match surface boundary conditions and to smoothly blend with

the interior,

Countergradient flux term The second term of the LMD scheme's surface boundary layer

formulation is the non-local transport term, which can playa significant role in mixing during

surface cooling events. This is a redistribution term included in the tracer equation separate from

the diffusion term and is written as

(133)

(130)

(128)

(125)

(126)

(127)

(124)

Rig < 0,

0< Rig < Rio,

Rig> Rio.

27

-0 WTO+WTR

,T - s wT(a)h

where 110 is 5.0 X 10-3, Rio = 0.7.

for temperature and

8

-8zK,. (129)

LMD base their formulation for non-local scalar transport on a parameterization for pure free

convection from Mailhot and Benoit [32]. They extend this parameterization to cover any unstable

surface forcing conditions to give

wSo

,s = Os ws(a)h (131)

for salinity (other scalar quantities with surface fluxes can be treated similarly). LMD argue that

although there is evidence of non-local transport of momentum as well, the form the term would

take is unclear so they simply specify ,m = o.

The interior scheme The interior scheme of Large, McWilliams and Doney estimates the viscosity coefficient by adding the effects of several generating mechanisms: shear mixing, double-diffusive

mixing and internal wave generated mixing.

IIx(d) = II~ + II~ + II:: (132)

ao = 0

al = 1

a2 = -2 + 3I1x (hsbl) + 8xllx(h) + .....:lIx.:...;(_h~)8-=-u..,....w-"-x-'-(1-'-)

hwx(1) wx(1) hwi(1)

a3 =1 _ 2IJx(hsbl) _ 8xllx(h) _ IIx(h)8uwx(1)

'hwx(1) wx(1) hwi(1)

where IIx (h) is the viscosity calculated by the interior parameterization at the boundary layer depth.

The shape function The non-dimensional shape function G(a) is a third order polynomial

with coefficients chosen to match the interior viscosity at the bottom of the boundary layer and

Monin-Obukhov similarity theory approaching the surface. This function is defined as a 3rd order

polynomial.

•

•

•

•

•

•

•

'.•

:.'.I.'.'.•

•

'.".•

'.•

•

•

•

•

'.•

•

:.'.

•

•

'.•

•

•

•

•

•

•

•

•

I.

•
28

3.12.2 Radiation boundary condition

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(136)

(140)

(137)

(135)

(134)

(138)

(139)

for 1.0 < Rp < R~ = 1.9,

otherwise.

v~ = 1.0 x 1O-4m2s-1

v:: = 1.0 x 1O-5m2s~l

vg = (1.5- 6 )(0.90gexp(4.6exp[-0.54(R;1 -1)])

d_ {Vg(1.85 - 0.85Rp1)Rp for 0.5 <= R p < 1.0,

Vs - d

v1J0.15Rp otherwise.

Double diffusive processes The second component of the interior mixing parameterization represents double diffusive mixing. From limited sources of laboratory and field data LMD

parameterize the salt fingering case (Rp > 1.0)

r

In realistic domains, open boundary conditions can be extremely difficult to get right. There can

be situations where incoming flow and outgoing flow happen along the same boundary or even at

the same horizontal location. Orlanski [39] proposed a radiation scheme in which a local phase

velocity is computed and used to radiate things out (if it is indeed going out). This works well for

a wave propagating normal to the boundary, but has problems when waves approach the boundary

at an angle. Raymond and Kuo [46] have modified the scheme to account for propagation in all

three directions. In SCRUM, only the two horizontal directions are accounted for:

3.12.1 Gradient boundary condition

This boundary condition is extremely simple and consists of setting the gradient of a field to zero

at the edge. The outside value is set equal to the closest interior value. It is probably too simple

to be useful in realistic problems.

3.12 Open boundary conditions

Currently, SCRUM has two open boundary conditions options; a gradient condition and a radiation

condition. These options are available independently for each of the four sides and for each field

category. The four categories are 2-D momentum, 3-D momentum, free surface, and 3-D tracers.

Internal wave generated mixing Internal wave generated mixing serves as the background

mixing in the LMD scheme. It is specified as a constant for both scalars and momentum. Eddy

diffusivity is estimated based on the data of Ledwell et a1. [29]. While Peters et a1. [42] suggest

eddy viscosity should be 7 to 10 times larger than diffusivity for gradient Richardson numbers

below approximately 0.7. Therefore LMD use

For diffusive convection (0 < Rp < 1.0) LMD suggest several formulations from the literature

and choose the one with the most significant impact on mixing (Fedorov [11]).

for temperature. For other scalars,
These terms are evaluated at the closest interior point in a manner consistent with the timestepping

scheme used. The phase velocities are limited so that the local CFL condition is satisfied. They

are then applied to the boundary point using equation (140), again using a consistent timestepping

scheme. Raymond and Kuo give the form used for centered differencing and a leapfrog timestep

while SCRUM uses one~sided differences.

The radiation approach is appropriate for waves leaving the domain. A check is made to see

which way the phase velocity is headed. If it is entering the domain, a zero gradient condition is

applied.

29

•

•

•

•

'.•

•

•

•

•

•

'.•

;.'.:.

;'.•

i.

•

•

•

',.•

•

'.•

•

•

•

•

:.•

•

•

•

'.".•

'.•

'.••

Ie

where

p!!J!..

ex = 8~

(~r + (~r p81jJ

c 81)

y - (!!J!..)2 (!!J!..)2 8~ + 81)

P = _ a1/J

at

(141)

(142)

(143)
30

••.

:

•••••••••

••••••••••••••••••••••••••••••••••
..

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

..

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

.-

•

•

•

•

4 Details of the Code

4.1 Main subroutines

The main program is in serum.F. It calls the initialization routines and then calls either main3d

or main2d. A flow chart for main3d is shown in Fig. 6. The boxes refer to subroutines which are

described as follows:

depth2d Computes the evolving total depth of the water column that is associated with the 2-D

momentum equations. It also computes the coefficients used in the advection and viscosity

of 2-D momentum.

depth3d Computes the evolving depths of the model grid and its associated vertical transformation metric Hz. It also computes the coefficients which contain Hz and are used in the

horizontal advection of momentum and tracers and in the horizontal mixing.

fre2drhs Computes the forcing terms (RUslow' RVslow) for the 2-D momentum equations which are

held constant over the short time steps. These forcing terms contains the vertically

integrated terms from the 3-D momentum equations which are not considered in the 2~D

equations.

initial Does everything that needs to be done to start up the model run. It reads initial parameters and u, v, T, 5, and ( fields from disk or calls ana_initial. It then calculates the

remaining initial fields and opens the restart file. The flow chart for initial is shown in

Fig. 7.

inp_par Reads in input model parameters from standard input. It also writes out these parameters

to standard output and calls eheekdefs.

omega Calculates the scaled vertical velocity Hzn/mn according to equation (79).

prsgrd Calculates the horizontal pressure gradients according to equation (80).

rho_eos Calculates the density anomaly, p, using the equation of state (§3.8).

seLvbe Sets the vertical boundary conditions for momentum and tracers.

step2d Time steps free-surface and 2-D momentum equations. Also does the time-averaging

described in §3.5 and calls the 2-D boundary condition routines.

step3d Time steps the 3-D momentum and tracers (usually potential temperature and salinity),

using the tridiagonal solver described in §3.4. Couples 3-D and 2-D momentum fields. It

computes and removes the vertical means from the newly computed 3-D velocities and

replaces those means with the more accurate 2-D velocities. It also calls the boundary

condition routines.

trhs Calculates and stores contributions to the right-hand-side of the tracer equations (21)

and 22), where the advective terms have been moved to the right-hand-side.

u2drhs Calculates and stores contributions to the right-hand-side of equation (54), where all the

terms other than %t (~~) have been moved to the right-hand-side.

u3drhs Calculates and stores contributions to the right-hand-side of equation (19), where all the

terms other than %t (If.;::) have been moved to the right-hand-side.

v2drhs Calculates and stores contributions to the right-hand-side of equation (55), where all the

terms other than %t (~~) have been moved to the right-hand-side.

31
loop over short'steps

Figure 6: Flow chart of the model main program.

32

••••••••••••••••••••••••••••••••••••••••••••
'.••••••

~.••••

..•

'It••••••

'.•

'.••

4D

4D

4D

4D

4D

-D

'~D

4D

~D

4D

4D

4D

4D

4D

4D

.D

~.•

grid file?

Figure 7: Flow chart of the initial subroutine.

33
v3drhs Calculates and stores contributions to the right-hand-side of equation (20), where all the

terms other than gt (~::) have been moved to the right-hand-side.

verLmix Computes the vertical mixing coefficients for momentum (Akv) and tracers (Akt).

zetarhs Computes the right-hand~side of equation (53), where all the terms other than gt (-~n)

have been moved to the right-hand-side.

4.2 Other subroutines and functions

Initialization

ana_grid Sets up an analytic grid.

ana_initial Sets up analytic initial conditions.

ana_mask Sets up an analytic mask.

blkdat Initializes some variables and parameters stored in common blocks.

checkdefs Reports on which C preprocessor variables have been #defined and checks

their consistency.

geLgrid Reads in the curvilinear coordinate arrays as well as f and h from a grid netCDF

file.

geLinitial Reads initial fields from disk-either restart or initializing from a climatology.

geLmask Reads in the mask arrays from the grid netCDF file. It also adjusts pmask

as required for the free-slip/no-slip boundary conditions as described in §3.2.

metrics Computes the metric term combinations which do not depend on the surface

elevation and therefore remain constant in time.

seLscoord Sets and initializes relevant variables associated with the vertical transformation to nondimensional s-coordinate described in Appendix B.

zero_arrays Initializes (zeroes out) various arrays.

NetCDF I/O

deLavg Creates the SCRUM averages NetCDF file and defines its dimensions, attributes,

and variables.

deLhis Creates the SCRUM history NetCDF file and defines its dimensions, attributes,

and variables.

deLrst Creates the SCRUM restart NetCDF file and defines its dimensions, attributes,

and variables.

deLstation Creates the SCRUM station NetCDF file and defines its dimensions, attributes, and variables.

geLdate Gets today's date, day of the week and time called. It uses Sun's intrinsic date

routine by default.

lenstr Returns the character position of the last non-blank character in a "string" after

removing the leading blank characters, if any. Should not be called with a literal

string argument.

opencdf Opens an existing NetCDF file, inquires about its contents, and checks for consistency with model dimensions.

wrLavg Writes SCRUM time-averaged fields into the averages NetCDF file.

34

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

eo

e

•

•

.'•

•

•

•

•

•

•

•

•

•

•

•

•

•
wrLhis Writes requested SCRUM fields at requested levels into the history NetCDF file.

wrt...rst Writes SCRUM fields into the restart NetCDF file.

wrt....station Writes out data into the stations NetCDF file.

Forcing fields The file analytic.F contains analytical formulations for computing various forcings

and initializations. For more realistic problems these fields are read from NetCDF files.

ana_bmflux Computes analytic kinematic bottom momentum flux.

ana_btflux Computes analytic kinematic bottom flux of tracer type variables.

ana....smflux Computes analytic kinematic surface momentum flux (wind stress).

ana....srflux Computes analytic kinematic surface shortwave radiation.

ana....ssh Computes analytic sea surface height and dQdSST which are used in the surface

heat flux correction.

ana....sst Computes analytic sea surface temperature and dQdSST which are used in the

surface heat flux correction.

ana....stflux Computes analytic kinematic surface flux of tracer type variables.

ana_tclima Computes analytic tracer climatology fields.

geLbmflux Reads'bottom momentum flux (bottom stress) from the forcing NetCDF file,

and then linearly time~interpolates to current model time.

geLbtflux Reads bottom flux of tracer type variables from the forcing NetCDF file, and

then linearly time-interpolates to current model time.

geLcycle Determines relevant parameters for time cycling of data from the forcing NetCDF

file. For instance, you may wish to use monthly means for each year of a multiyear run.

get....smflux Reads surface momentum flux (wind stress) from the forcing NetCDF file,

and then linearly, time-interpolates to current model time.

get....srflux Reads shortwave radiation flux from the forcing NetCDF file, and then linearly

time-interpolates to current model time.

get....ssh Reads sea surface height from the forcing NetCDF file and then linearly timeinterpolates to current model time.

get....sst Reads sea surface temperature and surface net heat flux sensitivity to sea surface

temperature from the forcing NetCDF file and then linearly time-interpolates to

current model time.

get....stflux Reads surface flux of tracer type variables from the forcing NetCDF file, and

then linearly time-interpolates to current model time.

geLtclima Reads climatology of tracer type variables from the climatology NetCDF file,

and then linearly time-interpolates to current model time.

set_nudgcof Set nudging coefficients time-scales (l/s).

Horizontal mixing The horizontal mixing routines have options for doing Laplacian or biharmonic mixing, along surfaces of constant $, z, or density, as described in §3.1O. The

horizontal mixing of 2-D momentum is computed in u2drhs and v2drhs.

geLepislope Computes the epineutral slopes (nondimensional) used in the mixing tensor

rotation relative to geopotential surfaces.

geLisoslope Computes the isopycnal slopes (nondimensional) used in the mixing tensor

rotation relative to geopotential surfaces.

35
seLhmixing Sets horizontal mixing coefficients. If requested, it also scales horizontal

mixing by the grid size and enhances horizontal mixing in the sponge areas.

shapd Applies a 2-D Shapiro filter to an array.

smagorinsky This routine computes horizontal mixing coefficients for momentum using

Smagorinsky parameterization (Tag et a1. [57]).

t3dmix Computes horizontal mixing of tracer type variables.

u3dmix Computes horizontal mixing of the 3-D momentum component in the ~-direction.

v3dmix Computes horizontal mixing of the 3-D momentum component in the 71-direction.

Vertical mixing The model contains a variety of methods for computing the vertical mixing

coefficients Akt and Akv, including an analytic formula.

ana_vmix Computes analytic vertical mixing coefficients for momentum and tracers.

bv_freq Computes the squared Brunt-Viiisiilii frequency at w-points N 2 = - t ~.

Imd_vmix Computes vertical mixing coefficients for momentum and tracers at the ocean

interior using the Large, McWilliams and Doney [28] mixing scheme.

Imd_bldepth Determines the oceanic planetary boundary layer depth, hbl, as

the shallowest depth where the bulk Richardson number is equal to the

critical value, Ric.

Imd_blmix Sets the vertical mixing coefficients within the boundary layer.

Imd...swfrac Computes the fraction of solar shortwave flux penetrating to specified depth (times Zscale) due to exponential decay in Jerlov water type.

Imd_wscale Computes the turbulent velocity scale for momentum and tracers

using a 2-D lookup table as a function of ustar and zetahat.

my25_vmix Computes vertical mixing coefficients for momentum and tracers using the

Mellor and Yamada [36] mixing level 2.5 scheme with modifications described in

Galperin et a1. [14].

my25_q Solves the prognostic equation for turbulent energy variables used in the

Mellor-Yamada level 2.5 turbulent closure.

pp_vmix Computes vertical mixing coefficients for momentum and tracers using the

Pacanowski and Philander [40] mixing scheme which is based on the Richardson

number.

rLnumber Computes the gradient Richardson number for the vertical mixing schemes.

trisolver Solves the PDE A¢(k - 1) + B¢(k) + C¢(k + 1) = D for field ¢ using the

tridiagonal solver, also known as Thomas algorithm (Richtmeyer and Morton

[48]).

Bottom boundary-layer model The model has an optional bottom boundary layer based on

Styles and Glenn [56].

ana_bsedim Computes analytic bottom sediment grain size and density.

ana_wwave Computes analytic wind induced wave amplitude, direction and period.

geLbsedim Reads initial sediment grain size and density from the forcing NetCDF file.

geLwwave Reads wind induced wave amplitude, direction and period from the forcing

NetCDF file, and then linearly time-interpolates to current model time.

sg_bbl96 Computes kinematic bottom momentum stress using Styles and Glenn [56]

bottom boundary layer formulation.

36

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
Boundary conditions The files be~2d.F and bcs3d.F coritain horizontal boundary conditions

for the various 2-D and 3-D variables, respectively. The boundary routines are also called

to specify boundary conditions on "V2¢ for the horizontal biharmonic operator on the field

¢.

u2dbc

u3dbc

v2dbc

v3dbc

w3dbc

'.•

:.•

•

•

•

•

•

•

•

•

•

•

'.•

•

'.•r.,i.

Ie•

•

•

•

•

,e•

,.

,e

•

•

,.

•

•

•

•

•

•

•

•

••

•

Other

sg_ubab Computes maximum wave bottom velocity and excursion from wind induced

wave amplitude and period by solving the linear wave dispersion relation for a

given wave number.

inflow Processes prescribed inflow open boundary conditions from climatology data.

obc_volcons Computes integral mass flux across all the open boundaries. Then, it corrects 2D velocities across to enforce global mass conservation.

t3dbc Boundary conditions for 3-D tracer type variables.

Boundary conditions for 2-D u-type variables.

Boundary conditions for 3-D u-type variables.

Boundary conditions for 2-D v-type variables.

Boundary conditions for 3-D v-type variables.

Boundary conditions for 3-D w-type variables.

xtrbry Extracts and loads data into boundary field arrays. These boundary field arrays

are used in the treatment of the open boundaries via radiation conditions.

zetabc Boundary conditions for free-surface type variables.

ab..ratio Calculates the ratio of the thermodynamic expansion coefficients for potential

temperature and salinity, alpha/beta, at horizontal and vertical w-points from a

polynomial expression (Jackett and McDougall [25]).

alfabeta Computes thermal expansion and saline contraction coefficients as a function

of potential temperature, salinity, and pressure from a polynomial expression

(Jackett and McDougall [25]).

ana_diag Computes customized diagnostics.

ana_meanRHO Analytical mean density anomaly (rhobar).

crash Dies in the manner appropriate for your computer (stop or call exit). It also

closes any open NetCDF files.

day_code computes a code for the day of the week, given the date. This code is good for

dates after January 1, 1752 AD, the year the Gregorian calendar was adopted in

Britain and the American colonies.

diag Computes various diagnostic fields, such as the volume averaged kinetic and

potential energies.

smoLadv Evaluates horizontal and vertical advection terms for tracers using the Smolarkiewicz [50] advection scheme. It uses an upstream advection scheme with

a second corrective upstream step to reduce the implicit diffusion. An antidiffusion velocity is computed and used in the second pass through the advection

operator.

smoLadiff Computes the "anti-diffusion velocity" used to suppress the numerical

diffusion that is associated with the upstream differencing operator for

advection.

37
smoLups Computes a first-order upstream differencing operator for the 3-D advection of a tracer (scalar) field.

wvelocity Computes vertical velocity (w) from the model vertical velocity (nHz/mn).

4.3 C preprocessor variables

Before it can be compiled, the model must be run through the C preprocessor cpp, as described

in Appendix F. The C preprocessor has its own variables, which may be defined either with an

explicit #define command or with a command line option to cpp. We have chosen to define these

variables in an include file, cppdefs.h, except for some machine-dependent ones, which are defined

in the appropriate Makefiles. These variables allow you to conditionally compile sections of the

code. For instance, if MASKING is not defined then the masking code will not be seen by the

compiler, and the masking variables will not be declared. These cpp variables can be grouped into

several categories:

model test problems One of these can be defined to obtain an example test problem.

BASIN Define for the "Big Bad Basin" example.

CANYON-A Define for the Canyon A (homogeneous) example.

CANYON...B Define for the Canyon B (stratified) example.

GRAV_ADJ Define for the gravitational adjustment example.

GSYRONT Define for idealized Gulf Stream front example.

MUNK Define for Stommel/Munk wind driven ocean basin.

. OVERFLOW Define for the overflow example.

RIVERPLUME Define for the river plume example.

SEAMOUNT Define for the seamount example.

SOLITON Define for the equatorial Rossby soliton example.

UPWELLING Define for the upwelling/downwelling example described in §6.2.

momentum terms

BODYFORCE Define to apply the surface stresses as a body force.

CURVG~ID Define to compute the extra non~linear terms which arise when using curvilinear coordinates.

UV_ADV Define to compute the momentum advection terms.

UV_COR Define to compute the Coriolis term.

UV_GSCHEME Define for third-order upwind advection scheme.

UV-FRS Define to compute the horizontal pressure gradient term.

UV-FSOURCE Define for point sources/sinks.

UV_VIS2 Define to compute the horizontal Laplacian viscosity.

UV_VIS4 Define to compute the horizontal biharmonic viscositY.

WJ-FRS Define for weighted Jacobian pressure gradient.

tracers

DIAGNOSTIC Define for a diagnostic calculation in which the tracer fields do not

change in time.

38

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

ICE Define to use ice component of the model (see §8.

ICE_THERMO Define for ice thermodynamics. ,

NONLIN...EOS Define to use the nonlinear equation of state.

QCORRECTION Define to use the net heat flux correction.

SALINITY Define if salinity is used as one of the tracers.

SMOLARKIEWICZ Define to compute Smolarkiewicz advection.

TS.A.DV Define to compute the tracer advection terms.

TS-DIF2 Define to compute the horizontal Laplacian diffusion.

TS-DIF4 Define to compute the horizontal biharmonic diffusion.

TS_GSCHEME Define for third-order upwind advection scheme.

TS...PSOURCE Define for point sources/sinks.

general model configuration

AVERAGES Define to write out time-averaged model fields.

RMDOCINC Define to remove documentation in include files with the C preprocessor.

SOLVE3D Define to solve the 3-D primitive equations.

STATIONS Define to write out time-series information at specific points in the model.

TIME.A.VG Define to average over short timesteps as described in §3.5.

analytic fields

ANA-BMFLUX Define for an analytic bottom momentum stress.

ANA-BSEDIM Define for an analytic bottom sediment grain size and density.

ANA-BSFLUX Define for an analytic bottom salt flux.

ANA-BTFLUX Define for an analytic bottom heat flux.

ANA-DIAG Define for customized diagnostics.

ANA_GRID Define for an analytic model grid set-up.

ANAJNITIAL Define for analytic initial conditions.

ANA_MASK Define for an analytic mask.

ANA_MEANRHO Define for an analytic mean density anomaly.

ANA...PSOURCE Define for analytic point sources.

ANA_SMFLUX Define for an analytic kinematic surface momentum stress.

ANA_SRFLUX Define for an analytic kinematic surface shortwave radiation.

ANA_SSFLUX Define for an analytic kinematic surface freshwater flux.

ANA_SSH Define for an analytic sea surface height.

ANA_SST Define for an analytic SST and fJQ/8SST.

ANA_STFLUX Define for an analytic kinematidmtface heat flux.

ANA_TCLIMA Define for an analytic tracer climatology.

ANA_VMIX Define for analytic vertical mixing coefficients.

ANA_WWAVE Define for an analytic wind induced wave field.

horizontal mixing of momentum

39
MIX_GP_UV Define for viscosity along constant z (geopotential) surfaces.

MIX...EPLUV Define for viscosity along constant in situ density (epineutral)

surfaces.

MIX-lSO_UV Define for viscosity along constant potential density (isopycnal)

surfaces.

MIX_S_UV Define for viscosity along constant s surfaces.

SMAGORINSKY Define for Smagorinsky mixing parameterization.

VIS_GRID Define for horizontally variable viscosity coefficient.

horizontal mixing of tracers

DIF_GRID Define for horizontally variable diffusion coefficient.

MIX_GP_TS Define for diffusion along constant z (geopotential) surfaces.

MIX...EPLTS Define for diffusion along constant in situ density (epineutral)

surfaces.

MIX-lSO_TS Define for diffusion along constant potential density (epineutral)

surfaces.

MIX_S_TS Define for diffusion along constant s surfaces.

vertical mixing

BVF_MIXING Define to activate Brunt-ViiisiiHi frequency mixing.

LMD_MIXING Define to activate Large/McWilliams/Doney interior closure.

LMD_CONVEC Define to add convective mixing due to shear instabilities.

LMD.-DDMIX Define to add double-diffusive mixing.

LMD-KPP Define to add boundary layer mixing from a local K-Profile Parameterization (KPP).

LMD~ONLOCALDefine to add convective nonlocal transport.

LMD-RIMIX Define to add diffusivity due to shear instabilities.

MY25_MIXING Define to activate Mellor/Yamada Level-2.5 closure.

Q.-DIF2 Define for horizontal Laplacian diffusion of q.

Q.-DIF4 Define for horizontal biharmonic diffusion of q.

Q_GSCHEME Define for third-order upwind advection of q.

PP_MIXING Define to activate Pacanowski/Philander closure.

SG-BBL96 Define to activate Styles/Glenn bottom boundary layer formulation.

boundary conditions

EW--PERIODIC Define for periodic boundaries in the i direction.

INFLOW2D Define to use SSH climatology as 2D inflow data.

INFLOW3D Define to process 3D inflow from tracer climatology.

NS--PERIODIC Define for periodic boundaries in the j direction.

detailed eastern open boundary conditions-other sides have similar

EAST..FSGRADIENT Define for a gradient condition on the free surface.

EAST--M2GRADIENT Define for a gradient condition on the 2-D momentum.

40

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
EAST-M2RADIATION Define for a radiation condition on the 2-D momentum.

EAST-M3GRADIENT Define for a gradient condition on the 3-D momentum.

EAST-M3RADIATION Define for a radiation condition on the 3-D momentum.

EAST_TGRADIENT Define for a gradient condition on the tracers.

EAST_TRADIATION Define for a radiation condition on the tracers.

general

MASKING Define if there is land in the domain to be masked out.

M2NUDGING Define for nudging to 2-D momentum data.

M3NUDGING Define for nudging to 3-D momentum data.

TCLIMATOLOGY Define for processing the tracer climatology arrays.

TNUDGING Define for nudging to tracer climatology.

ZCLIMATOLOGY Define for processing the sea surface height climatology arrays.

ZNUDGING Define for nudging to sea surface height climatology.

precision These variables were introduced so that one code could be used for Crays and worksta~

tions, all using 64 bit precision.

BIGREAL This is the type of all floating point variables used in the model computations.

It must be defined to be something, such as real or real*8. If this is set to

double precision you should also use a compiler option for extending source

lines past 72 characters in width.

DBLEPREC For double precision arithmetic.

FLoaT This is used so that the correct intrinsic is called, either float, dfloat or real.

OUT-DOUBLE For double precision output.

NF-FOUT Either nLdouble or nLreal, depending on OUT-DOUBLE.

NF-F'TYPE Either nLdouble or nf...real, depending on model precision.

command line These are defined as a command line option in some of the Makefiles since they

are machine dependent.

NO-EXIT This will determine whether your program ends with a stop command or by

calling exit. I prefer exit on a Sun and stop on an IBM RS/6000. The RS/6000

will not properly close files when using call exit so it is possible to lose some of

your output unless you use stop.

Note that the exit subroutine on many computers does not require an argument.

The Sun exit subroutine uses the integer argument value as the return code from

SCRUM for use by the shell under which SCRUM is run.

AIX Most versions of cpp which are supplied by the vendor have some variables

automatically defined. For instance, on a SparcStation, sun, unix, and sparc

will all be defined. However, the RS/6000 cpp does not define anything useful

to check for so I have the RS/6000 Makefile define AIX. This is used because

both the SGI and the IBM RS/6000 will continue to compute if some variables

have become NaN. In order to stop the calculation, we check for NaN as the

error from diag,but the method of checking varies from one system to another.

Another system-dependent component of SCRUM is in the implementation of

get-date.

41
4.4 Important parameters

The following is a list of the important parameters in the model. The rest ofthe parameters'defined

in param.h are derived from L, M, and N.

isalt Index into tracer arrays for salinity.

itemp Index into tracer arrays for temperature.

L Number of grid points in the ~~direction.

M Number of grid points in the l1-direction.

N Number of grid points in the vertical.

NS Maximum number of output station points.

Nsrc Maximum number of point sources/sinks.

NT Number of tracer fields. Often NT = 2 for potential temperature and salinity.

There are a lot of parameters defined inpconst.h to represent literal constants of type BIGREAL.

It is much safer to use the parameters when these values are needed as subroutine arguments. The

names for the constants were chosen based on the following "rules":

• Use a prefix of c for whole real numbers (cO for zero and C1 for 1.0).

• Use a prefix of p for non repeating fractions (p5 for 0.5).

• Use a prefix of r for reciprocals (r3 for 1.0/3.0 and rIO for 0.1 which could also be pI).

• Combine use of the prefix and e for scientific notation (cle4 for 1.0e + 4 and cIem4 for

1.0e - 4).

• Use names when appropriate (pi for 1r = 3.14159265 ... ).

42

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

5 Support Programs for Initialization

5.1 Grid generation

On startup, SCRUM either reads a NetCDF file or calls ana_grid to find the location of the grid

points, the grid metrics, the bathymetry, the land/sea mask, and the Coriolis parameter f. If you

won't be using ana~grid, the grid file must be generated before SCRUM can be run, either with

ezgrid or with the programs in gridpak. The version of ezgrid which produces a NetCDF file is

available in

ftp://ahab.rutgers.edu/pub/gridpak/ezgrid.shar

5.1.1 ezgrid

ezgrid was written to generate a uniform rectangular grid with a simple bathymetry. It has two

modes, one for the upwelling example, and one for rectangular basins; the mode is determined by

the UPWELLING switch in cppdefs.h. If UPWELLING is not defined then the important

parameters are:

xl basin length in the ~-direction.

el basin width in the 1]-direction.

hO bottom depth.

£0, beta Coriolis parameter with the ,B-plane approximation, f = fo + ,By.

In either case you will have to also set the name of the gridfile, grdname, near the top of the

ezgrid.F file. Once these parameters are set to your chosen values, compile and run it:

make ezgrid

ezgrid

This should create a binary NetCDF file called grdname.

5.1.2 gridpak

SCRUM has been designed to be used with curvilinear orthogonal grids for boundary-following

domains, etc., so there are situations in which you want a more flexible grid-generation program

than ezgrid. We have been working on a suite of programs called gridpak, including xeoast,

an interactive boundary drawing program. See §1.1 for instructions on obtaining gridpak and its

documentation.

5.2 Masking

·5.2.1 The serum_mask program

SCRUM now supports the masking of land areas, for which it requires some new input arrays.

These arrays are read from the grid NetCDF file or computed in ana..mask. The mask is defined

on p-points; see Fig. 8 for an example of a small domain with an isolated island and a promontory

adjacent to the boundary. There are also arrays for the mask on u-points, v-points, and '¢I-points

which are derived from the p-point mask. The '¢I-point mask depends on the free-slip/no-slip option

chosen as described in §3.2.

The programs in gridpak find the p-point mask based on the bathymetry dataset. Elevations

at or above sea level are assumed to be in the land mask. You may choose to edit this mask, so

43
j = M ,---....---....---....---....---....---r-:--,.--,.-----,

j=1

i = 1

Figure 8: Small grid with masked regions

Hernan Arango has written a Matlab tool called serum_mask. It is an interactive tool which

requires Matlab as well as mexedf for reading and writing NetCDF files from Matlab. It is

available in

ftp://ahab.rutgers.edu/pub/scrum/rnatlab/rnask/

ftp://ahab.rutgers.edu/pub/scrum/tars/scrum3_rnatlab.tar.gz

and includes a README file. An example of its use is shown in Fig. 9. This represents the same

mask as in Fig. 8, with a circle for each p-point, including the boundary "image" points. The red

circles are land while the blue circles are ocean. Notice that I have made the "image" points have

the same mask value as the points they mirror.

5.3 Objective Analysis

[This section was contributed by Hernan Arango.]

The objective analysis (oa) package described here can be used to prepare initial, climatology,

update, and forcing fields for SCRUM. It maps oceanographic and atmospheric data to a specified

application grid. Currently, it processes the following fields: in situ temperature, potential temperature, in situ density anomaly, salinity, sigma-t, sound speed, dynamic height, surface net heat

flux (Q), surface freshwater flux, precipitation rate, evaporation rate, incoming solar shortwave radiation, surface momentum (wind) stress components, sea surface temperature (SST), and surface

net heat flux sensitivity to SST (8Q/8SST).

This oa package is derived from an earlier program which Hernan Arango and Carlos Lozano

wrote at Harvard University in 1993. The basic algorithm used by this package is described in

Carter and Robinson [8]. A comprehensive description of this methodology can also be found in

Gandin [15], Bretherton et al. [7], McWilliams et al. [34], Daley [9], Bennett [6], and oth~rs.

Given observations Si = S(Xi, ti) at location Xi, ti, i = 1, ... N an estimate <PE of a scalar <P is

derived for location X and time t. A linear unbiased estimate is given by:

44

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

-

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

for arbitrary Wi since ¢E = ¢. The associated variance of error is:

with W == (WI, . .. W N). The overbar denotes an expected or ensemble mean value. The minimizer

W*:

is

A -I

W* = P

with minimum error variance (Gauss-Markov):

Here, for convenience, matrix notation has been used. s = [SI,'" SN] is a column correlation vector,

p == (¢ - ¢)(s - 8), and A is the covariance matrix:

A == (s - 8)(S - :5)'

where the prime denotes a transpose.

Notice that A is symmetric. In what follows, excluding pathological cases, A is assumed to be

positive definite. The best linear estimate ¢* is then:

with error e;.

The essential information required is statistical; namely the spatial-temporal mean of the scalar

and observations, the covariance between observations, and the covariance between the scalar and

the observations.

The observations can be of different types, and different from the scalar which you are trying

to find. Their usefulness is measured by the fractional reduction of error:

In this package it is assumed that the covariance of the scalar is homogeneous in space and homogeneous and isotropic in time:

and errors at two different locations and times are uncorrelated:

Currently, an analytical, isotropic, Gaussian correlation function is assumed:

with

C(r,7) = exp [-GJ'] G(r)

45
where To is the time decorrelation scale, a is the zero crossing distance, and b is the spatial decorrelation scale.

This package uses a local solution to the oa equations. That is, only nnce influential observations are considered at each mapped grid point. This method is practical because it avoids inverting

large matrices when the number of observations is large. Observations that are too far apart in

space and time from the mapped point contribute very little to the estimate, as one might expect.

It is available from:

ftp://ahab.rutgers.edu/pub/scrum/tars/scrum3_oa.tar.gz

and includes a README file.

5.4 Forcing fields

There are options for calling either ana...smfiux or get...smfiux to get the surface momentum

forcing. If you do not have an analytic formulation for this field, you will have to create a NetCDF

forcing file which contains the surface momentum fluxes. It can either contain one point value or

a 2-D field of values. Likewise, the field can be constant in time or contain values for a series of

times. It is even possible to have a limited number of snapshots which get cycled over in time. For

instance, you can provide 12 monthly mean fields and tell it to cycle over these in a multi-year run.

The other forcing fields are treated in the same way and are also contained in the NetCDF

forcing file. These include surface and bottom heat and salt fluxes, the 8Q/ aT and Tref terms

from §2.2, the incoming shortwave radiation used by the Large et a1. mixing scheme, cLnd the wave

information used by the Styles and Glenn bottom boundary layer. The ice thermodynamics also

requires forcing fields such as air temperature and cloud fraction.

An example program which creates the forcing NetCDF file is provided by the files in

ftp://ahab.rutgers/pub/scrum/forcing

This program reads a file produced by the oa package.

5.5 Initial and climatology fields

The model will either read its initial fields from a NetCDF file or it will compute them in analytical.F.

If it is not computing them, the routine geLinitial will read a history file or a file produced by

the initial program. This program in turn is expecting to read the output of the oa program. The

initial program is in

ftp://ahab.rutgers/pub/scrum/initial

The model has the option of reading in 3-D climatology fields from a climate NetCDF file. This

file containsthe 3-D climatologies for the tracers, perhaps at a number of times. The subroutine

geLclima will read this file and do any necessary time interpolations. The climate file is also

produced by the initial program. The climatology could also be used for the boundary conditions,

both for the tracer values on inflow or for prescribed boundary conditions. In this ease' it would

make more sense to only store the 2-D arrays. We do not yet have the software for handling these

2-D arrays, but it would be a straightforward modification to the initial program.

46

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
;.,;.••••••••••••••••••••••••••••••••••••••••••

47
I

I

••••••••••••••

.,••

.'.'•••••••••••••••••••••••••
•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

6 Configuring SeRUM for a Specific Application

This chapter describes the parts of SCRUM for which the user is responsible when configuring it

for a given application. Section 6.1 describes the process in a generic fashion while §6.2 and §6.3

step through the application of SCRUM to upwelling/downwelling and wind-driven North Atlantic

problems, respectively. As distributed, SCRUM is ready to run quite a few examples, where the

C preprocessor flags determine which is to be executed. Some of these examples are described in

Haidvogel and Beckmann [19]. Some of them are listed here:

BASIN This is a rectangular, flat-bottomed basin with double-gyre wind forcing. When run, it

produces a western boundary current flowing into a central "Gulf Stream" which goes

unstable and generates eddies. The goal is to run adiabatically to study the homogenization of potential vorticity. It takes a long time and caused difficulties for SPEM 3 so we

call it the "Big Bad Basin" .

CANYON-A The canyon is a periodic channel with a steep shelf along one wall, where the

shelf contains a steep canyon. There is a periodic forcing which causes the water to

oscillate along the channel. The rotation and the shelf lead to non-zero mean flows,

especially near the canyon. Version A is homogeneous and can be executed with a 2-D

model. See Haidvogel and Beckmann [18] for a description of the canyon problems and

the gravitational adjustment problem.

CANYON.J3 This is like Canyon A, except that it is stratified.

GRAV_ADJ The gravitational adjustment problem takes place in a long narrow domain which

is initialized with dense water at one end and light water at the other. At time zero, the

water is released and it generates two propagating fronts as the light water rushes to fill

the top and the dense water rushes to fill the bottom. This configuration was used to test

various advection schemes.

OVERFLOW This configuration is similar to the GRAV_ADJ problem, but is initialized with

dense water in the shallow part of a domain with a sloping bottom.

SEAMOUNT The seamount test was used to test the pressure gradient errors. It has an idealized

seamount in a periodic channel. See Beckmann and Haidvogel [5] and McCalpin [33] for

more information.

UPWELLING The upwelling/downwelling example was contributed by Anthony Macks and Jason Middleton [31] and consists of a periodic channel with shelves on each side. There

is along-channel wind forcing and the Coriolis term leads to upwelling on one side and

downwelling on the other side. If you run it for several days, you end up with dense water

over light water.

The input files for the SCRUM examples are included in the file:

ftp://ahab.rutgers.edu/pub/scrum/tars/scrum3_examples.tar.gz

This file also contains sample output NetCDF files and plot files and is quite large.

6.1 Configuring SCRUM

The three main files you need to change in SCRUM are scrum.in, cppdefs.h, and analytical.F.

These provide the input, set the options you want, and provide analytic formulas for various fields,

respectively. If more realistic fields are desired, you will have to provide other input files as well,

for instance for the grid and the wind forcing.

49
6.1.1 cppdefs.h and checkdefs.F

For each of the cpp variables described in §4.3, decide whether or not you want it to be defined.

Each defined variable should have a line of the form:

#define SOME_VAR '

Note that any undefined variable need not be mentioned, but we leave pliaceholders for them in

cppdefs.h as a reminder that they are meaningful. These placeholders can be in any of the

following forms:

#undef SOME_VARl

c #define SOME_VAR2

! #define SOME_VAR3

We use the first of these.

When configuring SCRUM for your problem, it is recommended that you add a newcpp variable

for it. New cpp variables can be added to cppdefs.h and then used in the code with an #ifdef

statement. This is a simple way to keep track of pieces that you add for your application. For

instance, my simple ice test is called MMS-BOX:

#ifdef MMS_BOX

# define EW_PERIODIC

# define NS_PERIODIC

# define ICE

If it becomes necessary to update to a newer version of SCRpM, it is simple to find the parts of

the code which belong to the Arctic version and copy them to the new SCRUM.

For each new cpp variable, it is recommended that you also add the appropriate code to

checkdefs.F, such as:

#ifdef ICE

write(stdout,20) 'ICE',

& 'Coupled sea-ice model.'

is=lenstr(Coptibns)+l

Coptions(is:is+4)=' ICE,'

#endif 1* ICE *1

Note that the number "4" on the Coptions line must be set according to the length of the string

you are adding. In this case 4 is for "ICE,", including the comma.

6.1.2 Model domain

One of the first things the user must decide is how many grid points to use, and can be afforded.

There are three parameters in param.h which specify the grid size and one parameter for the

number of tracers:

L Number of finite-difference points in e.

M Number of finite-difference points in "I.

N Number of finite-difference points in the vertical.

NT Number of tracers.

There are no constraints on these except L 2: 2, M 2: 2, N 2: 2 and NT 2: 1. Land :M should be

at least 3 if the domain is periodic in that direction.

50

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
I.'.'.•

•

'.':.•

'.:.•

'.•

•

•

•

:.•

:.•

'.•

'.:.•

•

•

'.•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

6.1.3 x, y grid

The subroutine get-grid or ana_grid is called by initial to set the grid arrays, the bathymetry,

and the Coriolis parameter. Most of the simple test problems have their grid information specified

in ana_grid in the file analytical.F. More realistic problems require a NetCDF grid file, produced

by the grid generation programs described in Wilkin and Hedstrom [61]. The variables which are

read by get-grid are:

xl, el, spherical, f, h, pm, pn, x..rho, y..rho, lon_rho, lat..rho, angle.

If the grid is curved, get-grid will also read:

dndx, dmde.

6.1.4 ~,'f/ grid

Before providing initial conditions and boundary conditions, the user must understand the model

grid. The fields are laid out on an Arakawa C grid as in Fig. 2. The overall grid is shown in Fig. 10.

The thick outer line shows the position of the model boundary. The points inside this boundary

are those which are advanced in time using the model physics. The points on the boundary and

those on the outside must be supplied by the boundary conditions.

The three-dimensional model fields are carried in three-dimensional arrays, except the tracers

where the fourth array index tells which tracer is being referred to. For instance, itemp = 1

refers to potential temperature while isalt = 2 refers to salinity. The integers i, j, and k are used

throughout the model to index the three spatial dimensions:

z Index variable for the ~-direction.

J Index variable for the 'f/~direction.

k Index variable for the a-direction. k = 1 refers to the bottom

while k = N refers to the surface.

6.1.5 Initial conditions

The initial values for the model fields are provided by either anajnitial or get-initial. get-initial

is also used to read a restart file if the model is being restarted from a previous run.

Also in initial, rho_eos is called to initialize the density field. rho_eos in turn calls ana_meanRHO

to initialize the rhobar array. rhobar is a function of z only and should be more or less the hor~

izontal average of the density field. It is subtracted from p, before p is vertically integrated in

prsgrd, to reduce the errors in the pressure gradient terms.

The tracer climatology fields also require appropriate values if they are to be used, and are

provided by ana_tclima or get_tclima. Likewise, the surface height climatology is read by get-ssh

or provided by ana-ssh.

6.1.6 Equation of state

The equation of state is defined in the subroutine rho_eos. Two versions are provided in SCRUM:

a nonlinear p = p(T, S, z) from Jackett and McDougall [25] and a linear p(T, S). The linear form is

p = RO + Tcoef· T + Scoef· S or p = RO + Tcoef· T, depending on whether or not SALINITY

is defined. Specify which equation of state you would like to use by setting the NONLIN-EOS C

preprocessor flag in cppdefs.h. The linear coefficients RO, Tcoef and Scoef are set in scrum.in.

51
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

0 () 0 0 0 0

- -

0 0 0 0 0 0

- - -= ....

0 0 0 0 0 0

-

0 0 0 0 0 0

o 0

Figure 10: The whole grid.

x - u points

0- v points

0- p points

D D

52

i = 0 1 1 2 2 3 3 Lm Lm L L

2 D D

100

1 D D

200

D D

0 X 0 X 0 X 0 X 0 X 0 X 0 X 0

~mD D

~mo 0

j=O 0 x 0 x 0 x 0 x 0 x 0 x 0 x 0
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

••

•

•

•

•

•

•

'.'.•

6.1.7 Boundary conditions

The horizontal boundary conditions are provided by the subroutines in bcs3d' and bcs2d. They

are called every timestep and provide the boundary values for the fields u, v, u, v, T, Sand (. They

are currently configured for a closed basin, a periodic channel, a doubly periodic domain or an open

domain with radiation conditions. The open boundary conditions are evolving; the latest efforts

are in the ROMS model.

6.1.8 Model forcing

(a) Winds and thermal fluxes

There are two different ways to apply a wind forcing: as a surface momentum flux in the vertical

viscosity term, or as a body force over the upper water column. In the past, our vertical resolution

was relatively coarse and the vertical viscosity would have to have been unreasonably large for us

to resolve the surface Ekman layer. If that is your situation, define BODYFORCE in eppdefs.h

and provide a value for levsfre in serum.in. The forcing is applied over the levels from levsfre

to N. The above caution about vertical resolution also applies to the surface fluxes of T and S,

although BODYFORCE only refers to wind stress, not the surface tracer fluxes.

More recently, we have been setting the vertical s-coordinate parameters to retain some resolution near the surface and to apply the fluxes as boundary conditions to the vertical viscosity/ diffusivity. In either case, the surface and bottom fluxes are either defined analytically or

read from the forcing file. You must either edit the appropriate parts of analytieal.F or create

a NetCDF forcing file in the format expected by get..smflux, get..stflux and their friends. Note

that it is quite common to put the wind stress into the forcing file while having an analytic bottom

stress. '

(b) Climatology

One way to force the model is via a nudging to the tracer climatologies. This was is used in the

North Atlantic simulations in sponge layers along the northern and southern boundaries. Set the

climatologies in ana_tclima or in a file read by get-tclima, set NUDGING in eppdefs.h and

also set the array nudgeof in seLnudgeof.F.

6.1. 9 serum.in

SCRUM expects to read a number of variables on standard in. It is easiest to prepare an input file

and then run SCRUM as:

serum < serum.in > serum.out &

The input is organized as pairs of lines, the first with a number and then some text which is ignored,

the second with the values for the set of variables for which the first line provides the key. The

pairs of lines can be in any· order but are usually sorted numerically, The number 99 signals the

end of these pairs and the rest of the input file contains comments for the user. The input pairs

are as follows:

1 Time-stepping parameters.

ntimes Number of timesteps to evolve the 3-D equations in the current run. This is

actually the total number, including any previous segments of the same run. For

instance, if you already did a three-month run and wish to continue for another

three months, set ntimes to the number of steps needed for six months. If

you don't like this and would prefer to have the behavior of the SPEM variable

ntmes, modify main.F so that:

53
becomes

do iic=ntstart,ntimes

54

Isopycnal thicknesses.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

Constant mixing coefficient for the horizontal Laplacian diffusion of each tracer.

A value is expected for each of the NT tracers.

Constant mixing coefficient for the horizontal biharmonic diffusion of each tracer.

A value is expected for each of the NT tracers.

Record number of the restart file to read as the initial conditions.

Record number of the ice restart file to read as the iniltialconditions (coupled

ice model only).

Number of timesteps between writing of restart fields.

Number of timesteps between writing fields into the history file.

Starting timestep for the accumulation of output time-averaged data. For instance, you might want to average over the last day of a thirty-day run.

Number of timesteps between writing time-averaged data into the averages file.

Number of timesteps between' writing data into stations file.

Number of timesteps between printing a single line of diagnostic information to

the standard output.

Logical switch used to create the history file. If .true., a new history file is

created. If .false. and nlev > 0, data is appended to an existing history file.

Logical flag used to recycle time records in the restart file.' If .true., only the

latest two restart time records are retained. If .false., all restart fields are saved

every nrst timesteps without recycling.

tnu2

tnu4

nrst

nwrt

ldefhis

navg

nsta

ninfo

lcycle

ntsavg

Laplacian horizontal mixing of tracers.

Kdiff Isopycnal mixing thickness diffusivity for each tracer variable. A value is expected for each of the NT tracers. These isopycnal thicknesses are used when

the Gent/McWilliams isopycnic mixing is activated (not currently implimented).

Biharmonic horizontal mixing of tracers.

do iic=ntstart,ntimes-l+ntstirt

nrrec

nrreci

dt Timestep in seconds for the 3-D equations.

ndtfast Number of timesteps for the 2-D equations to be executed each dt.

Input/Output parameters. SCRUM has several possible output files. The output files

include a restart file, a history file, an averages file, and a statiOn file. The restart file

often contains only two records with the older record being overwritten during the next

write. The history file can contain a subset of the restart fields, for instance just the

surface elevation and the surface temperature. The averages file contains time-averages

of the model fields, for instance montly means, or yearly means, depending on navg. The

station file contains timeseries for specified points, possibly quite frequently since each

record is small.

4

3

5

2
12 Vertical s~coordinates parameters.

11 Various parameters.

9 Mellor-Yamada Level 2.5 parameters.

10 Bottom drag coefficients.

55

Mean density used in the Boussinesq approximation.

Linear bottom drag coefficient.

Quadratic bottom drag coefficient.

Bottom roughness.

rhoO

rdrg

rdrg2

Zo

akv_bak Background vertical mixing coefficient for momentum.

nmix_en Number of timesteps between computations of isopycnal slopes used in the

rotated mixing tensor.

adv_ord Order of advection scheme when using Smolarkiewicz advection. A value of

adv_ord = 2is recommended to suppress the diffusive nature of the "upwind"

scheme. A value of adv_ord = 1 will yield the standard "upwind" advection.

levsfrc Deepest level to apply surface momentum stresses as a body force. Used when

the C-preprocessor option BODYFORCE is defined.

levbfrc Shallowest level to apply bottom mbmemtum stresses as a body force. Used

when the C-preprocessor option BODYFORCE is defined.

akq_bak Background vertical mixing coefficient for turbulent kinetic energy.

q2nu2 Constant mixing coefficient for the horizontal Laplacian diffusion of turbulent

kinetic energy.

q2nu4 Constant mixing coefficient for the horizontal biharmonic diffusion of turbulent

kinetic energy.

theta....s s-coordinate surface control parameter, [0 < theta....s < 20].

theta_b s~coordinate bottom control parameter, [0 < theta_b < 1].

Tcline Width of the surface or bottom boundary layer in which higher vertical resolution

is required during stretching.

WARNING: Users need to experiment with these parameters. We have found

out that the model goes unstable with high values of theta....s. With steep and

v~ry tall topography, it is recommended that you use theta....s 3.0.

13 Mean Density and time stamp.

7 Vertical mixing coefficients for tracers.

akLbak Background vertical mixing coefficient for the tracers. A value is expected for

each of the NT tracers.

6 Horizontal viscosity coefficients.

uvnu2 Constant mixing c'oefficient for the horizontal Laplacian viscosity. ,.

uvnu4 Constant mixing coefficient for the horizontal biharmonic viscosity.

8 Vertical mixing coefficient for momentum.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

;.•

•

•

•

•

'.•

!.•

i.

•

•

'.'.•

•

•

'.'.'.•

i,.•

•

•
56

15 Linear equation of state parameters.

19 Logical switches to activate the writing of other fields into the NetCDF history file:

16 Slipperiness parameters.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

Write out tracer type variables: potential temperature, salinity, etc.

Background density value used in the linear equation of state.

Background potential temperature constant used in analytieal.F.

Background salinity constant used in analytical.F.

Thermal expansion coefficient in the linear equation of state.

Saline contraction coefficient in. the linear equation of state.

wrtT

RO

TO

SO

Teoef

Seoef

wrtU Write out 3-D u-velocity component.

wrtV Write out 3-D v-velocity component.

wrtW Write out 3~D w-velocity component.

wrtO Write out 3-D n vertical velocity.

wrtUBAR Write out 2-D u-velocity component.

wrtVBAR Write out 2-D v-velocity component.

wrtZ Write out free~surface.

dstart Time stamp assigned to model initialization (days). Usually a Calendar linear

coordinate, like modified Julian day. For example:

dstart = 10200 corresponds to May 1, 1996

It is called modified Julian day because an offset of 2440000 needs to be added.

rnudg Time scale (days) of nudging towards clima~ology at the interior and at the

boundaries.

gamma2 Slipperiness variable, either 1.0 (free slip) or -1.0 (no slip).

wall1 Logical switch for side 1 (i = 1), .true. if it is a wall, .false. if it is open.

wall2 Logical switch for side 2 (j = 1), .true. if it is a wall, .false. if it is open.

wa1l3 Logical switch for side 3 (i = L), .true. if it is a wall, .false. if it is open.

wa1l4 Logical switch for side 4 (j = M), .true. if it is a wall, .false. if it is open.

18 Logical switches to activate the writing of fields associated with the tracer equations into

the NetCDF history file. A value is expected for each of the NT tracers.

17 Logical switches to activate the writing of fields associated with the momentum equations

into the NetCDF history file:

14 Nudging/relaxation time scales.

Znudg Time scale (days) of nudging towards sea surface height climatology.

M2nudg Time scale (days) of nudging towards 2-D momenum climatology.

M3nudg Time scale (days) of nudging towards 3~D momenum climatology.

Tnudg Time scale (days) of nudging towards tracer climatology. A value is expected for

each of the NT tracers.
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

wrtRHO Write out density anomaly.

wrtAKV Write out vertical viscosity coefficient.

wrtAKT Write out vertical diffusion coefficient for temperature.

wrtAKS Write out vertical diffusion coefficient for salinity.

wrtHBLWrite out depth of the planetary boundary layer.

20 Number and Levels to output:

nlev Number of levels to write out to the history file for each activated 3~D field. If

nlev < 0, all model levels are written out. IF nlev = 0, the history file will not

be created.

lev If nlev > 0, levels to write out to the history file. nlev values are expected:

1 :::; lev(l : nlev) :::; N

Enter values in ascending numerical order.

21 String with a maximum of eighty characters.

title Title of the model run.

22 String with a maximum of eighty characters.

rstname Output restart file name (NetCDF).

23 String with a maximum of eighty characters.

hisname Output history file name (NetCDF).

24 String with a maximum of eighty characters.

avgname Name of the file for the averaged model fields (NetCDF).

25 String with a maximum of eighty characters.

staname Name of the file for the station output (NetCDF).

26 String with a maximum of eighty characters.

fltname Name of the file containing the float output (NetCDF). Not implemented yet.

27 String with a maximum of eighty characters.

grdname Name of the file containing the grid data (NetCDF).

28 String with a maximum of eighty characters.

ininame Name of the file containing the initial conditions. It can be a SCRUM restart

file (NetCDF).

29 String with a maximum of eighty characters.

frcname Name of the file containing the forcing fields (NetCDF).

30 String with a maximum of eighty characters.

clmname Name of the file containing the climatology fields (NetCDF).

57
31 String with a maximum of eighty characters.

assname Name of the file containing the assimilation fields (NetCDF).

32 String with a maximum of eighty characters.

aparnam Name of the file containing the assimilation parameters (ASCII).

33 String with a maximum of eighty characters.

sposnam Name of the file containing the stations positions (ASCII).

34 String with a maximum of eighty characters.

fposnam Name of the file containing the initial drifter positions {ASCII).

35 String with a maximum of eighty characters.

iininame Name of the ice input file (coupled ice model only).

36 String with a maximum of eighty characters.

irstnam Name of the ice restart file (coupled ice model only).

37 String with a maximum of eighty characters.

ihisnam Name of the ice history file (coupled ice model only).

38 String with a maximum of eighty characters.

iavgnam Name of the ice averages file (coupled ice model only).

39 String with a maximum of eighty characters.

usrname User's generic input file name.

An example input file without the trailing comments is:

1 NTIMES, DT (s), NDTFAST

1800 240.dO 20

2 NRREC, NRST, NWRT, NTSAVG, NAVG, NSTA, NINFO, lLDEFHIS, LCYCLE

0 360 360 1 360 1 1 T T

3 TNU2 [1 :NT] (m~2/s)

5.dO 5.dO

4 TNU4[1 :NT] (m~4/s)

1.0d+07 1.0d+07

5 Kdiff [1 :NT] (m2/s)

O.dO O.dO

6 UVNU2 (m~2/s), UVNU4 (m~4/s)

10.dO O.dO

7 AKT_BAK [1 :NT] (m~2/s)

1. Od-5 1.0d-5

8 AKV_BAK (m~2/s)

1.0d-4

9 AKQ_BAK (m~2/s) Q2NU2 (m~2/s), Q2NU4 (m~4/s)

1.0d-4 20.dO 1.0d+07

10 RDRG (m/s), RDRG2, 20 (m)

58

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

4. 5E-04 I 0 .dO 0 .dO

11 NMIX_EN, AOV_ORO, LEVSFRC, LEVBFRC

1 2 i 1

12 THETA_S, THETA_B, TCLINE (m)

3.dO O.dO 50.dO

13 RHOO (Kg/m-3), DSTART (days)

1025.dO O.dO

14 ZNUDG (days), M2NUDG (days), M3NUDG(days) , TNUOG[1 :NT] (days)

O.dO O.dO O.dO O.dO O.dO

15 RO (Kg/m-3) TO (deg C), SO (PSU) , TCOEF, SCOEF

1026.9524 10.dO 35.dO -1.67e-04 7.62e-04

16 GAMMA2, WALL1, WALL2, WALL3, WALL4

1.dO T F F F

17 wrtU, wrtV, wrtW, YrtO, wrtUBAR, wrtVBAR, wrtZ

T T T F F F T

lB wrtT(l:NT) (temperature, salinity, etc.)

T T

19 yrtRHO, wrtAKV, wrtAKT, wrtAKS, wrtHBL

F F F F F

20 NLEV, LEV(l:NLEV) in ascending order (if NLEV<O, all levels are saved)

-1 1 3 5

21 TITLE (aBO)

Scrum 4.0

22 RSTNAME (aBO): SCRUM output restart file name, if any.

scrum_rst.nc

23 HISNAME (aBO): SCRUM output history file name, if any.

scrum_his.nc

24 'AVGNAME (aBO): SCRUM output averages file name, if any.

scrum_avg.nc

25 STANAME (aBO): SCRUM output stations file name, if any.

scrum_sta.nc

26 FLTNAME (aBO): SCRUM output floats file name, if any.

scrum_flt.nc

27 GRDNAME (aBO): SCRUM input grid file name, if any.

scrum_grd.nc

2B ININAME (aBO): SCRUM input initial conditions file name, if any.

scrum_ini.nc

29 FRCNAME (aBO): SCRUM input forcing fields file name, if any.

scrum_frc.nc

30 CLMNAME (aBO): SCRUM input climatology fields file name, if any.

scrum_clm.nc

31 ASSNAME (aBO): SCRUM input assimilation fields file name, if any.

scrum_ass.nc

32 APARNAM (aBO): SCRUM input assimilation parameters file name, if any.

assimilation.dat

33 SPOSNAM (aBO): SCRUM input station positions file name, if any.

stations. dat

34 FPOSNAM (aBO): SCRUM input initial floats positions file name, if any.

floats .dat

35 USRNAME (aBO): USER's input/output generic file name, if any.

59
/dev/null

99 END of input data

6.1.10 User variables and subroutines

It is possible for the user to add new variables and common blocks appropriate to a given application.

It is also possible to add new subroutines, for instance to read in river inflow data. If you create

new source files they will have to be added to the Makefile or the Imakefile (see §H). Also, any

new #include statements will have to be listed in the Makefile dependencies. The simplest way

to add them is to run make depend.

6.2 Upwelling/Downwelling Example

One application for which SCRUM has been configured is a wind-driven upwelling and downwelling

example, described in Macks and Middleton [31J. There is a shelf on each wall of a periodic channel

and an along-channel wind forcing, which drives upwelling at one wall and downwelling at the other.

This problem depends on the Ekman layer, so a surface stress is used with vertical viscosity. The

Ekman depth is estimated to be 9 m if Av = 0.01m2 / s, so the vertical grid spacing must resolve

this. The maximum depth is 150 m and our choice of the vertical grid parameters leads to a surface

l::1z of 4.0 m.

6.2.1 cppdefs.h

The C preprocessor variable UPWELLING has been introduced to make sure that we can

#define UPWELLING and have a consistent upwelling configuration of the model. This is

done in part within cppdefs.h by

#ifdef UPWELLING

#define UV_ADV

#undef UV_VIS2

#define UV_PRS

#define UV_COR.

#define TS_ADV

#undef TS_DIF2

#undef NONLIN_EOS

#undef SALINITY

#undef CURVGRID

#define EW_PER.IODIC

#undef NS_PERIODIC

#define TIME_AVG

#undef BODYFORCE

#define ANA_GRID

#define ANA_INITIAL

#define ANA_MEANRHO

#define ANA_SMFLUX

#define ANA_STFLUX

#define ANA_SSFLUX

#define ANA_BTFLUX

#define ANA_BSFLUX

#define ANA_VMIX

#endif /* UPWELLING */

60

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

.1
•

•

•

•

•

•

•

•

•

:.•

,'.:.

•

:.'.

•

'.•

•

'.•

•

•

•

•

,.

•

•

'.•

•

•

•

•

•

•

'.•

•

I.•

•

•

Here we have declared that we want a periodic channel but no masking. There is neither salinity nor

climatology. The momentum equations have the Coriolis and pressure gradients, but no horizontal

viscosity. The only term in the tracJr equation, is the advection.

6.2.2 Model domain

The flow does not vary in x, so L can be small. Set the values for L, M, N and NT in param.h:

L = 42

M = 81

N = 16

NT=2.

6.2.3 ana_grid

For this geometry one has a choice of using the grid-generation programs described in Wilkin and

Hedstrom [61], or of using ana_grid to create the grid analytically. The ana_grid subroutine in

analytical.F was modified to produce a bathymetry with a shelf on both walls of the channel when

UPWELLING is defined. The fluid depth ranges from 27 m on the shelves to 150 m in the center

of the channel. The horizontal grid spacing is uniform at 1 km and the Coriolis parameter f is set

to a constant value suitable for Sydney, Australia.

6.2.4 Initial conditions and the equation of state

We would like the initial conditions to be a motionless fluid with an exponential stratification.

ana-initial is configured accordingly.

The stratification can be provided by either T or S, or by both T and S. For simplicity we

will only have an active temperature field and we will use the linear equation of state by setting

NONLIN.-EOS to #undef in rhsdefs.h. We want the density to be 26.35 at the bottom and

24.22 at the top with an e-folding scale of 50 meters. The initial temperature is set to 14 + 8ez/

50

in ana-initial. The linear equation of state parameters are set in serum.in.

Since density does not depend on salinity, we have a choice of how to handle the second tracer.

We can either use it as a passive tracer or not timestep on it at all by setting NT = 1. We will use

it as a passive tracer and initialize it to be a function of y.

We have set ana_meanRHO to the desired initial density field. The climatology fields are not

used and need not be initialized.

6.2.5 Boundary conditions

The periodic channel options have already been chosen in cppdefs.h. We do not have to do

anything else.

6.2.6 Model forcing

In this problem we want to resolve the surface Ekman layer and to use a surface wind stress rather

than a body force. We want the amplitude of the wind to ramp up with time so we modify

ana..smflux accordingly. The wind will build to an amplitude of 0.1 Pascals / Po, or 1O-4m2s-2.

We need to edit ana_vmix to make sure that the vertical viscosity Akv is set to the value we

want. This must be large at the surface (10- 2m2s-l) to create a thick Ekman layer, but has been

chosen to decrease with depth. We also need to check that ana..sbflux, ana..stflux, etc. are set

correctly.

61
6.2.7 scrum.in

6.2.8 Output

The model writes some information to standard out, after setting ninfo to 72:

The model has been set up to run for one day with an internal timestep of 120 s and an external

timestep of 12 s. We will write history and restart records every 1/4 day. The value of the linear

bottom friction coefficient rdrg is set to 4.5 x 10-4 and the channel walls are set to be free-slip.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

62

Number of timesteps to evolve 3-D equations.

Timestep size (s) for 3-D equations.

Number of timesteps for 2-D equations between ea.ch OT.

Number of restart records to read from disk.

Number of timesteps between storage of restart fields.

Number of timesteps between writing fields into

history file.

Number of timesteps between print of information

to standard output.

Switch to create a new history MetCOF file.

Switch to recycle time-records in restart NetCOF file.

Horizontal, Laplacian mixing coefficient (mA2/s)

for tracer 1.

Horizontal, Laplacian mixing coefficient (mA2/s)

for tracer 2.

Horizontal, Laplacian mixing coefficient (mA2/s)

for momentum.

Background vertical mixing coefficient (mA2/s)

for tracer 1.

Background vertical mixing coefficient (mA2/s)

for tracer 2.

Background vertical mixing coefficient (mA2/s)

for momentum.

Linear bottom drag coefficient (m/s).

Quadratic bottom drag coefficient.

S-coordinate surface control parameter.

S-coordinate bottom control parameter.

S-coordinate surface/bottom layer width (m) used

in vertical coordinate stretching.

Mean density (kg/mA3) used in Boussinesq approximation.

Time stamp assigned to model initialization (days).

Nudging/relaxation inverse time scale "c1/days)

for free-surface.

Nudging/relaxation inverse time scale (l/days)

for 20 momentum.

Nudging/relaxation inverse time scale O/days)

for 3D momentum.

Nudging/relaxation inverse time scale (l/days)

for tracer 1.

SCRUM input parameters:

720 ntimes

120.00 dt

10 ndffast

0 nrrec

180 nrst

180 nwrt

72 ninfo

T ldefhis

T lcyele

O.OOOE+OO tnu2(1)

O.OOOE+OO tnu2(2)

O.OOOE+OO uvnu2

O.OOOE+OO Akt_bak(1)

O.OOOE+OO Akt_bak(2)

1.000E-05 Akv_back

4.500E-04 rdrg

O.OOOE+OO rdrg2

3.000E+00 theta._s

O.OOOE+OO theta_b

50.0000 Teline

1000.0000 rhoO

0.0000 dstart

O.OOOE+OO Znudg

O.OOOE+OO M2nudg

O.OOOE+OO M3nudg

O.OOOE+OO Tnudg(l)
Activated C-preprocessing Options:

Upwelling/Downwelling Example on a Periodic Double Shelf Channel

Output/Input Files:

scrum_rst.nc

scrum_his.nc

/dev/null

Nudging/relaxation inverse time scale (l/days)

for tracer 2.

Background potential temperature (Celsius) constant.

Background salinity (PSU) constant.

Background density (kg/m- 3)' used in linear Equation

of State.

Thermal expansion coefficient (l/Celsius).

Saline contraction coefficient (l/PSU).

Slipperiness variable: free-slip (1.0) or

no-slip (-1. 0) .

Boundary for side 1 (i=l): wall/open (T!F).

Boundary for side 2 (j=l): wall/open (T/F).

Boundary for side 3 (i=L): wall/open (T/F).

Boundary for side 4 (j=M): wall/open (T/F).

Write out 3D U-momentum component (T/F) ,

Write out 3D V-momentum component (T/F).

Write out W-momentum component (T/F).

Write out omega vertical velocity (T/F).

Write out 2D U-momentum component (T/F).

Write out 2D V-momentum component (T/F).

Write out free-surface (T/F).

Write out tracer 1 (T/F).

Write out tracer 2 (T/F).

Write out density anomaly (T/F).

Write out vertical viscosity coefficient (T/F).

Write out vertical T-diffusion coefficient (T/F).

Write out vertical S-diffusion coefficient (T/F).

Number of levels to write out.

Levels to write out:

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16

Analytical kinematic bottom salt flux.

Analytical kinematic bottom heat flux.

Analytical grid set-up.

Analytical initial conditions.

Analytical mean density anomaly.

Analytical kinematic surface momentum flux.

Analytical kinematic freshwater (E-P) flux.

63

Tcoef

Scoef

garmna2

TO

SO

RO

Tnudg(2)

Output Restart File:

Output History File:

Input/Output USER File:

F wall1

T wall2

F wall3

T wall4

T wrtU

T wrtV

T wrtW

T wrtO

T wrtUBAR

T wrtVBAR

T wrtZ

T wrtT(1)

T wrtT(2)

T wrtRHO

F wrtAKV

F wrtAKT

F wrtAKS

16 Nlev

Lev

ANA_BSFLUX

ANA_BTFLUX

ANA_GRID

ANA_INITIAL

ANA_MEANRHO

ANA_SMFLUX

ANA_SSFLUX

0.0000

0.0000

30.3795

O.OOOE+OO

-2.800E-Ol

O.OOOE+OO

1.00
Vertical S-coordinate System:

level S-coord at hmin over slope at hmax

16 0.00 0.00 0.00 0.00

15 -0.06 -1.71 -2.87 -4.02

14 -0.12 -3.43 -5.78 -8.12

13 -0.19 -5.14 -8.77 -12.39

12 -0.25 -6.86 -11.89 -16.92

11 -0.31 -8.57 -15.18 -21.80

10 -0.38 -10.28 -18.71 -27.14 9 -0.44 -12.00 -22.54 -33.08 8 -0.50 -13.71 -26.74 -39.76 7 -0.56 -15.42 -31.40 -47.37 6 -0.62 -17 .14 -36.62 -56.09 5 -0.69 -18.85 -42.52 -66.20 4 -0.75 -20.57 -49.27 -77.97 3 -0.81 -22.28 -57.02 -91.76 2 -0.88 -23.99 -66.00 -108.01 1 -0.94 -25.71 -76.46 -127.21 0 -1.00 -27.42 -88.71 -150.00

MAIN- started time-stepping SCRUM:

Day= 0.100000 avgKE

= 7.090495E-19 avgPE

= 1.697521E-08

Day= 0.200000 avgKE

= 5.151371E-18 avgPE

= 1.697524E-08

Day= 0.300000 avgKE

= 1.904604E-17 avgPE

= 1.697527E-08

Day= 0.400000 avgKE

= 4.972789E-17 avgPE

= 1.697529E-08

Day= 0.500000 avgKE

= 1.058779E-16 avgPE

= 1.697532E-08

Day= 0.600000 avgKE

= 1.973340E-16 avgPE

= 1.697534£-08

Day= 0.700000 avgKE

= 3.345132E-16 avgPE

= 1.697535£-08

Day= 0.800000 avgKE

= 5.280177E-16 avgPE

= 1.697536£-08

Day= 0.900000 avgKE

= 7.890830E-16 avgPE

= 1.697537E-08

Day= 1.000000 avgKE

= 1. 130497E-15 avgPE

= 1.697536E-08

ANA_STFLUX

ANA_VMIX

DBLEPREC

EW_PERIDDIC

MIX_GP_TS

MIX_GP_UV

SOLVE2D

SOLVE3D

TIME_AVG

TS_ADV

UPWELLING

UV_ADV

UV_COR

UV_PRS

Analytical kinematic surface heat flux.

Analytical vertical mixing coefficients.

Double precision arithmetic.

East-West periodic boundaries.

Mixing of tracers along geopotential surfaces.

Mixing of momentum along geopotential surfaces.

Solving 2D Primitive Equations.

Solving 3D Primitive Equations.

Time averaging over two short timestep cycles.

Advection of tracers.

Upwelling/Downwelling Example.

Advection of momentum.

Coriolis term.

Hydrostatic pressure gradient term.

64

•.1•••••••••••••••••••••••

••••••••••••••••••••
Main - number of time records written in history file: 0005'

number of time records written in restart file: 0002

Main Done.

NetCDF history and restart files are also created, containing the model fields at the requested

times. We have asked it to save both history and restart records every 1/4 day. In this case, the

restart file has been told to "cycle", or to write over the second last record. The restart file at the

end of the run contains the fields at 3/4 day and 1 day. The history file contains records for 0, 1/4,

1/2, 3/4, and 1 day. Plots can be made from either file, using the plotting software described in

§7. Selected frames from such plots are shown in Fig. 11 to 14.

6.3 North Atlantic example

The upwelling/downwelling examples is one in which all the start-up fields are defined analytically.

The other extreme is one in which everything is read from files, as in our North Atlantic simulations.

6.3.1 cppdefs.h

The C preprocessor variable DAMEE-B has been introduced to make sure that we can #define DAMEE-B

and have a consistent configuration of the model. This is done in part within cppdefs.h by

# if defined DAMEE_B I I defined DAMEE_S

#define UV_ADV

#define UV_GSCHEME

#undef UV_VIS2

#undef UV_VIS4

#define UV_PRS

#define UV_COR

#undef MIX_GP_UV

#define TS_ADV

#define TS_GSCHEME

#undef TS_DIF2

#undef TS_DIF4

#undef MIX_GP_TS

#undef SMOLARKIEWICZ

#define NONLIN_EOS

#define SALINITY

#undef DIAGNOSTIC

#define QCORRECTION

#define CURVGRID

#define AVERAGES

#undef STATIONS

#undef OBC_EAST

#undef OBC_WEST

#define OBC_NORTH

#define OBC_SOUTH

#undef EW_PERIODIC

#undef NS_PERIODIC

#undef INFLOW

#define OBC_TPRESCRIBE

65
seRUM 3.0

Wind-Driven Upwelling/Downwelling Test overa Periodic 'Channel

Theta s=3, Theta b=O, Tcline=50; nu(v,t)=(O,O) m2/s

0.00 Day

Figure 11: The upwelling/downwelling bathymetry.

••••••••••••••••••••••••••••••••••••••••••••

Ho"day -,~p"'ll 28.1997- 5011001 PH

sCr'u..j"ls."C

119.4

143.9

88.7

f------j 33.5

·~~58.1

25

66

o km

Min= 2.7421E+Ol Max= 1.5000E+02

Bathymetry (m)

50
seRUM 3.0

Wind-Driven Upwelling/Downwelling Test over a Periodic Channel

Theta s=3, Theta b=O, Tcline=50; nu(v,t)=(O,O) rn2/s

1.00 Day

25

50

Min= 9.5674E+OO Max= 1.6064E+Ol

Total Velocity Vectors (crn/s) at Level 16

l1onda., • April 28. 1qq7 - 4019.'1 PM

lu~ru"J'I •. nc

Figure 12: Surface velocities after One day, showing the flow to the left of the wind (southern

hemisphere).

67
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

-268E

-175.

542.1

-4.2

1.4

3.0

I10nday - Apr-II 2B, 1'~q1 - 4,23.27 Pf1

.e,.um_~l •.nc

o km 50 Min=~2.B649E+03 Max= 7.2146E+02

Omega Vertical Velocity (em/day)

o km 50

Min=-4.6253E+OO Max= 3.3B06E+OO

Total V-velocity (cm/s)

100

-2670.3

-1035.5

-127.2

-15.4

-12.7

-9.3

-5.9

-3.2

seRUM 3.0

Wind-Driven Upwelling/Downwelling Test over a Periodic Channel

Theta s=3, Theta b=O, Tcline=50; nu(v,t)=(O,O) m2/s

1.00 Day

68

Figure 13: Constant slices of the u, v, wand n fields at day 1.

Max= 7.B099E+02

Vertical Velocity (em/day)

o km

Min: -1.6053E+O 1 Max= -2.5105E+OO

Total U-velocity (cm/s)

100

100
Figure 14: Constant slices of the T, S (tracer), kinetic energy and Ertel potential vorticity at

day 1.

1.4

0.5

0.5

Ho"da, - ,t,p,.ll 28, ,QQ7 - 4,23.27 Pr1

lIc,.ulll_nllll,nc

o km

Min= 2.B973E-IO Max= 4.7607E-09

Ertel Potential Vorticity (l/s/m)

50

Min= B.7490E-02 Max= 7.B562E+OO

Salinity (PSU)

100

100

16.3

14.9

0.10

x10- 5

seRUM 3.0

69

Wind-Driven Upwelling/Downwelling Test over a Periodic Channel

Theta s=3, Theta b=O, Tcline=50; nu(v,t)=(O,O) m2/s

1.00 Day

o km

Min= 1.4504E+OI Max= 2.IB46E+OI

Potential Temperature (C)

Kinetic Energy (Walls)

1=021 1=021

00 J=00~1~~~~~~~~~J~=079

100

100

.-

..

.:

•

..ie'.'.'.e;.'.•

I)'.I)

•

I.•

e

•

•

••

•

•

'.'.•

•

•

•

•

•

•

•

•

'.•

•

•

•

•i;_•
70

6.3.2 Model domain

A large number of horizontal grid points was chosen to resolve the domain at less than one degree.

Values for L, M, N, and NT are:

L = 129

M = 129

N = 20

NT=2.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

MASKING

TIME_AVG

BODYFORCE

BVF_MIXING

PP_MIXING

LMD_MIXING

LMD_RIMIX

LMD_CONVEC

LMD_DDMIX

LMD_KPP

CLIMATOLOGY

NUDGING

ANA_MEANRHO

ANA_SMFLUX

ANA_SSFLUX

ANA_STFLUX

ANA_SRFLUX

ANA_BSFLUX

ANA_BTFLUX

ANA_V2DBC

/* DAMEE_B I I DAMEE~S */

#define

#define

#define

#undef

#undef

#define

#define

#define

#undef

#define

#define

#define

#define

#undef

#undef

#undef

#undef

#define

#define

#undef

# endif

Here, we have declared that we want a closed basin (not periodic), masking, salinity, and the nonlinear equation of state. We want Laplacian viscosity and diffusion along constant z-surfaces and

the full non-linear, curvilinear momentum equations.

We also added the DAMEE flags to checkdefs.F:

#ifdef DAMEE_B

write(stdout,20) 'DAMEE_B',

& 'North Atlantic DAMEE Big Domain Application.'

is=lenstr(Coptions)+l

Coptions(is:is+9)=' DAMEE_B,'

iexample=iexample+l

#endif /* DAMEE_B */

#ifdef DAMEE S

write(stdout,20) 'DAMEE_S',

& 'North Atlantic DAMEE Small Domain Application.'

is=lenstr(Coptions)+l

Coptions(is:is+9)=' DAMEE_S,'

iexample=iexample+l

#endif /* DAMEE_S */
'.•

•

•

:.•

'.:.'.,i.'.'.'.

,.

•

•

'.•

•

:.•

•

•

•

•

•

•

•

•

•

•

,.

•

•

•

•

•

;.•

•

•

•

••

•

6.3.3 gridpak

The grid has uniform spacing on a Mercator projection so that both ~x and ~y get smaller as you

get farther from the equator. The grid was chosen to go from 30° S to 65° N and was generated with

sqgrid. We then found the latitude and longitude values with tolat and interpolated the etop05

bathymetry to the grid with bathtub. The grid is shown in Fig. 15 and the unsmoothed bathymetry

is shown in Fig. 16. It is clear that the unsmoothed bathymetry contains some incredibly steep

regions. We have not pushed SCRUM to see what its steepness limit is, but we also ran SPEM in

this configuration and its elliptic solver requires substantial smoothing at this resolution. We were

advised by Bernard Barnier to retain the shallow island arc in the Caribbean. We also had some

bad experiences with shelves that disappeared into the land mask, such as, at Cape Hatteras and

the Iberian peninsula. We filled in the Pacific and the Mediterranean and did some unspeakable

hacking to bathsuds to obtain the bathymetry shown in Fig. 17. We then ran sphere to obtain

the values of m and n suitable for a spherical Earth and ran Hernan Arango's mask editing tool

serum_mask.

6.3.4 Initial conditions

We would like the initial conditions to be a motionless fluid with temperature and salinity fields

from the Levitus 1994 February mean climatology. We prepared a NetCDF file with zero u, v and

( fields. The T and S fields were interpolated from the Levitus fields-we tried several different

interpolation/extrapolation techniques, including the oa program described in §5.3.

An analytic function for the mean density was added to anaJDeanRHO for this problem:

# elif defined DAMEE_B I I defined DAMEE_S

do k=l,N

do j=O,M

do i=O,L

rhobar(i,j,k)=30.5-0.004*z_r(i,j,k)-

& c4*exp(z_r(i,j,k)/2000.0)

enddo

enddo

enddo

6.3.5 Boundary conditions

The non-periodic option has already been chosen by not defining EW_PERIODIC or NS-PERIODIC

in cppdefs.h. After trying a number of options, we ended up with walls to the north and south

with nudging regions (see below).

6.3.6 Forcing

The forcing is provided by surface momentum, heat and salt fluxes from the COADS dataset. We

apply the heat flux correction (#define QCORRECTION), which is also provided in COADS.

We use the oa program to put the values onto the model grid for each of the twelve monthly means.

6.3.1 Climatology

We used the same Levitus temperature and salinity fields for the climatology as for the initial

conditions. The DAMEE problem was specified to have nudging to the climatology at the northern

and southern boundaries, as well as at the Straits of Gibraltar. We edited seLnudgcof.F to set

the nudgcof array accordingly.

71
72

••••••••••••••••••••••••••• ~-~- ••••

,,

•••••••••••••

I . I, , , , , , ,

, .1,

..

,,, , , , , , ,

Figure 15: The North Atlantic grid.

North Atlantic DAMEE #4
'.•••••••••••••

:.'.••

,i.••

'.••

'.••

'.••••

'.

,.•••

'.••••e

,e

Bottom Topography MIN DEPTH= 200,000

MAX DEPTH= 5500,0

5.00E+03

4.50E+03

4,OOE+03

.', 3.50E+03

3,OOE+03

2.50E+03

2,OOE+03 , ,50E+03

1.00E+03

5.00E+02

Figure 16: The raw bathymetry from etopo5.

73
74

Figure 17: The smoothed North Atlantic bathymetry.

Bottom Topography MIN DEPTH

- 200.000

MAX DEPTH- 5500.0

.1•••••••••••••••••••••••••••••••••••••••••••
6.3.8 serum.in

The model writes out information to standard out:

6.3.9 Output

75

Number of timesteps to evolve 3-D equations.

Timestep size (s) for 3-D equations.

Number of timesteps for 2-D equations between each OT.

Number of restart records to read from disk.

Number of timesteps between storage of restart fields.

Number of timesteps between writing fields into

history file.

Starting time step for the accumulation of output

time-averaged data.

Number of timesteps between writing of time-averaged

data into averages file.

Number of time steps between print of information

to standard output.

Switch to create a new history NetCOF file.

Switch to recycle time-records in restart NetCOF file.

Background vertical mixing coefficient (m-2/s)

for tracer 1.

Background vertical mixing coefficient (m-2/s)

for tracer 2.

Background vertical mixing coefficient (m-2/s)

for momentum.

Linear bottom drag coefficient (m/s).

Quadratic bottom drag coefficient.

Deepest level to apply surface stress as a body force.

Shallowest level to apply bottom stress as a body force.

S-coordinate surface control parameter.

S-coordinate bottom control parameter.

S-coordinate surface/bottom layer width (m) used

in vertical coordinate stretching.

Mean density (kg/m-3) used in Boussinesq approximation.

Time stamp assigned to model initialization (days).

Background potential temperature (Celsius) constant.

Background salinity (PSU) constant.

Slipperiness variable: free-slip (1.0) or

no-slip (-1. 0) .

Boundary for side 1 (i=l): wall/open (T/F).

Boundary for side 2 (j=l): wall/open (T/F).

Boundary for side 3 (i=L): wall/open (T/F).

Boundary for side 4 (j=M): wall/open (T/F).

SCRUM input parameters:

144000 ntimes

2160.00 dt

20 ndffast

0 nrrec

400 nrst

1200 nwrt

1 ntsavg

1200 navg

1 ninfo

T ldefhis

T lcycle

1.000E-05 Akt_bak(1)

1.000E-05 Akt_bak(2)

1.000E-04 Akv_back

3.000E-04 rdrg

0.000£+00 rdrg2

18 levsfrc

1 levbfrc

5.000E+00 theta_s

4.000E-Ol theta_b

200.0000 Teline

1000.0000 rhoO

30.0000 dstart

0.0000 TO

0.0000 SO

1.00 gamma2.

T walll

T wa1l2

T wa1l3

T wa1l4

We use an internal timestep of2160 s~?d an external timestep of 108 s. The horizontal viscosity and

diffusion is turned off (see above) sint~'-the third-order upwInd scheme is providing the smoothing.

The stretching parameters are(} = 5J'lJi_,=.4 and he = 200m.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

'.•

•

•

•

•

•

I.•

•

•

•

:.•

•

•

•

•

•

•

•

••

•
Scrum 3.0 - North Atlantic Damee 4: Annual Levitus 0.75 resolution

Activated C-preprocessing Options:

Output/Input Files:

Output Restart File:

Output Averages File:

Input Grid File:

Input Initial File:

Input Forcing File:

Input Climatology File:

Input/Output USER File:

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

scrum_rst.nc

scrum_avg.nc

damee_grid3.nc

damee_Iev_4feb.nc

frc_coads_4.nc

damee_clm_4L.nc

/dev/null

Write out 3D U-momentum component (T/F).

Write out 3D V-momentum component (T/F).

Write out W-momentum component (T/F).

Write out omega vertical velocity (T/F).

Write out 2D U-momentum component (T/F).

Write out 2D V-momentum component (T/F).

Write out free-surface (T/F).

Write out tracer 1 (T/F).

Write out tracer 2 (T/F).

Write out density anomaly (T/F).

Write out vertical viscosity coefficient (T/F).

Write out vertical T-diffusion coefficient (T/F).

Write out vertical S-diffusion coefficient (T/F).

Write out depth of mixed layer (T/F).

76

Analytical kinematic bottom salt flux.

Analytical kinematic bottom heat flux.

Analytical mean density anomaly.

Writing out time-averaged fields.

Momentum stresses as body-forces.

Orthogonal curvilinear grid.

North Atlantic DAMEE Big Domain Application.

Double precision arithmetic.

LMD convective mixing due to shear instability.

Large/McWilliams/Doney interior mixing.

Large/McWilliams/Doney boundary layer mixing.

LMD diffusivity due to shear instability.

Land/Sea masking.

Mixing of tracers along geopotentialsurfaces.

Mixing of momentum along geopotential surfaces.

Non-linear Equation of State for seawater.

Nudging toward climatology.

Open North boundary edge.

Open South boundary edge.

T wrtU

T wrtV

T wrtW

F wrtO

T wrtUBAR

T wrtVBAR

T wrtZ

T wrtT(1)

T wrtT(2)

F wrtRHO

T wrtAKV

T wrtAKT

T wrtAKS

T wrtHBL

ANA_BSFLUX

ANA_BTFLUX

ANA_MEANRHO

AVERAGES

BODYFORCE

CURVGRID

DAMEE_B

DBLEPREC

LMD_CONVEC

LMD_MIXING

LMD_KPP

LMD_RIMIX

MASKING

MIX_GP_TS

MIX_GP_UV

NONLIN_EOS

NUDGING

OBC_NORTH

OBC_SOUTH
77

MAIN - started time-stepping SCRUM:

GET INITIAL - Processing initial conditions for time = O.OOOOE+OO

GET SRFLUX - Read solar shortwave radiation for time = 345.0

GET STFLUX - Read surface flux of tracer 01 for time = 345.0

GET_STFLUX - Read surface flux of tracer 02 for time = 345.0

GET SMFLUX - Read surface momentum stresses for time = 345.0

GET_TCLIMA - Read climatology of tracer 01 for time = O.OOOOE+OO

GET_TCLIMA - Read climatology of tracer 02 for time = O.OOOOE+OO

Tracers, boundary reduced physics condition.

Surface net· heat flux correction.

Using salinity.

Solving 2D Primitive Equations.

Solving 3D Primitive Equations.

Processing tracer climatology data.

Time averaging over two short timestep cycles.

Nudging toward tracer climatology.

Advection of tracers.

Shchepetkin-McWilliams G-Scheme Advection of tracer.

Advection of momentum.

Coriolis term.

Shchepetkin-McWilliams G-Scheme Advection of momentum.

Hydrostatic pressure gradient term.

OBC_TREDUCED

QCORRECTION

SALINITY

SOLVE2D

SOLVE3D

TCLIMATOLOGY

TIME_AVG

TNUDGING

TS_ADV

TS_GSCHEME

UV_ADV

UV_COR

UV_GSCHEME

UV_PRS

Vertical S-coordinate System:

level S-coord at hmin over slope at hmax

20 0.00 0.00 0.00 0.00

19 -0.05 -10.00 -20.03 -30.05

18 -0.10 -20.00 -43.30 -66.60

17 -0.15 -30.00 -71.92 -113.84

16 -0.20 -40.00 -108.94 -177.89

15 -0.25 -50.00 -158.64 -267.27

14 -0.30 -60.00 -226.50 -393.01

13 -0.35 -70.00 -318.60 -567.19

12 -0.40 -80.00 -439.47 -798.94

11 -0.45 -90.00 -588.95 -1087.90

10 -0.50 -100.00 -759.64 -1419.28

9 -0.55 -110.00 -938.48 -1766.95

8 -0.60 -120.00 -1112.90 -2105.81

7 -0.65 -130.00 -1277.09 -2424.19

6 -0.70 -140.00 -1433.59 -2727.18

5 -0.75 -150.00 -1591.00 -3032.00

4 -0.80 -160.00 -1761.00 -3361.99

3 -0.85 -170.00 -1956.64 -3743.28

2 -0.90 -180.00 -2192.17 -4204.35

1 -0.95 -190.00 -2483.64 -4777.28

0 -1.00 -200.00 -2850.00 -5500.00
78

It also writes out NetCDF files for restart, history, and monthly averages. Plots can be made from

all three of these files; an example plot is shown in Fig. 18.

Day = Day = Day =

GET_SRFLUX - Read solar shortwave radiation for time = 15.00

GET_STFLUX - Read surface flux of tracer 01 for time 15.00

GET_STFLUX - Read surface flux of tracer 02 for time = 15.00

GET_SMFLUX - Read surface momentum stresses for time 15.00

0.025000 avgKE = 9.587033E-16 avgPE 7.947536E-07

0.050000 avgKE = 9.608381E-16 avgPE = 7. 947534E-07

0.075000 avgKE = 9.583708E-16 avgPE = 7.947537E-07

•

•

•

•

•

•

•

•

•

•

•

•

•f)

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
0.61

056

0.50

0,45

1 0 .

39

1--.......-1

1

o. 34

I!-----l 0.29

I!-----l 0.23

I!-----l 0 18

1f----1 O. 12

1f----1 0.07

jL----1002

10----1 -0,04

lL=-_-"J - 0 .09

-0.15

-020

-0.25

-0.31

1 -0.36

r.----ll- 0.42

-0.47

-052

-0.58

-063

-0.69

-074

-0,79

-085

-090

-0.96

79

Annual

Figure 18: The annual mean surface elevation for year 10.

North Atlantic Damee #4: Run 27

LMD mixing, G-Scheme, nu(v,t)=(O,O), Theta=(5,0.4), N=20

Sea surface elevation

ROMS 1.0

Min~-1.0090E+OO Max= 6.1071E-Ol

-100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10

Ie'.•

•

•

•

•

•

•

•

•

•

•

•

•

I.•

I.•

I.

IeIeI.•

I.

•

•

•

Ie•

•

•

'.•

•

•

•

•

•

•

•

•

•

•
o••••••8 1 ••••••••••

..•••••••••••••••••••••••••
81

cent creates color plots of the horizontal fields, including constant depth plots of the 3-D fields.

7 Plotting Programs for Postprocessing

1,2,3,4,5

o FRSTD

o LASTD

o DSKlP

o VlNTRP

0.0 PMlN

0.0 PMAX

1 lCNT

0.0 lSOVAL

1.2 VLWD

2.0 VLSCL

1 lVlNC

1 JVlNC

o lREF

25 !DOVER

1 LEVOVER:

0.0 RMlN

0.0 RMAX

10.0 LGR!D

4 lPROJ

-60.0 PLON

NLEVELS: number of levels and/or depths to plot (0 for all levels)

levels (>0) or depths «0) to plot: FLDLEV(l:NLEVELS)

NFlELDS: number of fields to plot. Line below, field(s) types:

field identification: FLD!D(l:NFlELDS)

NLEVELS: number of levels and/or depths to plot (0 for all levels)

levels (>0) or depths «0) to plot: FLDLEV(l:NLEVELS)

first day to plot

last day to plot

plot every other DSK!P days (0.0 plot at its own time frequency)

vertical interpolation scheme: O=linear, l:cubic splines

field minimum value for color palette (0.0 for default)

field maximum value for color palette (0.0 for default)

draw contours between color bands: O=no, l=yes

iso-surface value to process (see below)

vector line width (1.0 for default)

vector length scale (1.0 for default)

vector grid sampling in the X-direction (1 for default)

vector grid sampling in the Y-direction (1 for default)

secondary or reference field option (see below)

overlay field identification (for lREF=l,2 only)

level of the overlay field (set to 0 if same as current FLDLEV)

overlay field minimum value to consider (0.0 for default)

overlay field maximum value to consider (0.0 for default)

Desired longitude/latitude grid spacing (degrees)

map projection (see below).

projection Pole longitude (west values are negative).

1,2

1

sec creates black-and-white plots of vertical slices through the 3-D fields. It includes on option

of finding equal-spaced points along a straight line through the curvilinear grid.

esee creates color plots of vertical slices through the 3-D fields. It includes on option of finding

equal-spaced points along a straight line through the curvilinear grid.

All of these program come with example input files. For instance, the input file for cent is called

eent.in and is as follows:

8 NFlELDS: number of fields to plot. Line below, field(s) types:

42,45,46,47,48,49,50,121,122,123,124 field identification:

FLDlD(l:NFlELDS)

1

20

2

1996 -1 year and starting year-day (use yearday<O, for no time label)

SCRUM 3.0

Coarse Arctic ocean with Budgell ice dynamics

ice thermodynamics

Hernan Arango has provided SCRUM with some programs for creating plots from the NetCDF

history and restart files. Some example plots are shown in §6. There are four plotting programs:

ent creates black-and-white plots of the horizontal fields, including constant depth plots of

the 3-D fields.

.•

•

'.•

•

•

•

'.'.•

•

•

•

'.•

•••

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•
82

*** Above FILENAMES:

*** IPROJ: Map Projections option. Some of the values for th.e

projection Pole and rotation angle are suggested.

c

c===================~==================================================:=

C Copyright (c) 1996 Rutgers University ----

c=======================================================================

c

•

•

•

•

•

•

•

•

•

.i•

.1

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

projection Pole latitude (south ·values are negative).

projection rotation angle (clockyise; degrees).

flag to color mask land: [0] no, [1] yes

number of plots per page (currently 1, 2, or 4)

logical syitch to read in positions from grid NetCDF file.

logical 5Yitch dray Logo.

logical syitch to yrite out the plot header titles.

logical syitch to yrite out the plot bottom title.

logical syitch to yrite out data range values and CI.

PLAT

ROTA

LMSK

NPAGE

READGRD:

PLTLOGO:

WRTHDR :

WRTBLAB:

WRTRANG:

*** IREF: Secondary or reference field option:

-1 Overlay horizontal grid

o no secondary or reference field to plot

1 plot field overlay from primary file

2 plot field overlay from secondary file

3 primary - secondary file (field subtraction) '<

4 DayO - DayN (field subtraction)

1st line: input; variables ID file.

2nd line: input; color palette file.

3rd line: input; contour parameters.

4th line: input; primary NetCDF file.

5th line: input; secondary NetCDF file.

6th line: input; grid NetCDF file.

7th line: input; coastlines file.

T WRTFNAM: logical syitch to yrite out input primary filename.

T WRTDATE: logical syitch to yrite out current date.

T CST logical syitch to read and plot coastlines and islands.

50.0 50.0 bottom and top map latitudes (south values are negative).

-110.0 80.0 left and right map longitudes (yest values are negative).

/d2/kate/plot/varid.dat

/d1/arango/scrum3.0/plot/Palettes/gs1.pal

/d2/kate/plot/default.cnt

ice_rst.nc

scrum_rst_1.nc

/d2/kate/arctic/gridpak/arctic_grid_2.nc

/u1/coasts/coast.dat

90.0

0.0

o

-1

T

F

T

T

T
83

Check the NCAR manual for details.

Plotting Fields classification: (* derived fields)

total velocity component in the XI-direction (cm/s).

total velocity component in the ETA-direction (cm/s).

total velocity vectors (em/s).

total velocity vector magnitude (cm/s).

IDUTOT

IDVTOT

IDTVEC

IDTMAG

1]

2]

3]

4]

[1] Cylindrical equiaistant: PLON=O, PLAT~O, ROTA=O

[2] Mercator: PLON=?, PLAT=O, ROTA=O

[3] Lambert conformal conic: PLON=?, PLAT=?, ROTA=?

[4] Stereographic azimuthal: PLON=?, PLAT=90 or -90, ROTA=O

*** IVINC, JVINC: vector grid sampling. If either value is negative,

the velocity vectors at drawn at PSI-points. Otherwise,

if both values are positive, the vectors are drawn at

interior RHO-points.

*** NPAGE: Number of plots per page. Set this parameter to a negative

value (-1, -2, or -4) to activate preservation of the plot

aspect ratio.

*** LGRID: Longitude/latitude grid spacing. The grid is drawn at

LGRID spacing starting from specified map lower corner.

If LGRID is negative, the latitude labels in the map are

rotated 90 degrees to avoid label congestion, if any.

As you can see, there are comments describing what needs to be done. Please see the variable ID

file for the complete list of fields which can be plotted~this list changes as Hernan adds the ability

to plot new fields. Also, check your default.cnt file for other vector and contour parameters. The

palette file includes two number systems, one in the scale from 0 to 255 and the other from 0 to l.

The SeRUM plotting uses the first while the SEOM plotting uses the second.

•

•

•

'.'-•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
84

••••••••••••••••••••••••••••••••••••••••••••
8.1 Model structure

8 Ice Model Formulation

The nonlinear advection terms have been omitted, since they are usually much smaller than the

others. Nonlinear formulas are used for both the ocean-ice and air-ice surface stress:

(152)

(151)

(149)

(150)

(146)

(147)

(148)

(144)

(145)

fa = paCa\VlO\VlO

Ca = ~Cd [1 - cos(27l" min(hi + .1, .5)]

fw== PwCwlvw- vl(vw- v).

85

a [ ov au] a [ ( au ov ) ] :F = - , (77 + ()- + (( - 77)- - P/2 + - 77 -' + -

Y oy oy ax ax oy ax

a [ au ov ] a [ ( au ov ) ] :Fx = - (77 + ()- + (( - 77)-" - P/2 + - 77 - + - ax ax oy oy oy ax

with the result that

P

(1ij == 277Eij + (( - 77)Ekk 8ij - 28ij.

The viscous~plastic terms (151) and (152) are found by taking the divergence of the stress tensor:

The symbols used in these equations along with the values for the constants are listed in Table 3.

A key component of the momentum balance is the force due to the internal ice stress (:Fx and

Fy ). This force is based on a constitutive law which relates the ice stress to the strain rate and

ice strength (equation (149)). For this model, a viscous~plastic behavior is used. Rigid plastic

behavior is approximated in this law by allowing the ice to flow in a plastic manner for normal

strain rates and to creep in a linear viscous manner for small strain rates. The treatment of ice

as a viscous-plastic fluid was largely motivated by the desire to avoid the complexities associated

with elastic~plastic behavior under flow. The stress~strain relationship is given by

The overall structure consists of two principal components-the momentum equations and the ice

continuity equations. The momentum balance includes air and water stress, Coriolis force, internal

ice stress, inertial forces and ocean t~lt as shown in equations (144) and (145):

MOU o(w x x ,.,.. at = Mf v - M9 ax + T a + Tw+.rx

MOV = '-Mfu - Mgo(w + T Y + rY +:F . at oy a w Y

Hibler [22] has described a model fdr the simulation of sea ice circulation and thickness. The

model has been tested over a seasonal cycle and the results are also shown in that article. This

report was derived from the documentation for the sea ice model written by Hibler [23] since the

dynamical equations are much the same. The model itself has been rewritten by Paul Budgell and

is now implemented on an orthogonal, curvilinear Arakawa C-grid, has a new elliptic solver, and

the nonlinear advection of momentum has been omitted. The thermodynamics are derived from

Mellor and Kantha [35] (MK89 below). Sirpa Hakkinen allowed us to use her implementation of

MK89; we obtained it from the Norwegians, who call it "hakkis".

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•
Variable Value Description

A(x, y, t) ice concentration

a(A) ridging function

Ca nonlinear air drag coefficient

Cd 2.2 x 10-3 air drag coefficient

Cw 10 x 10-3 water drag coefficient

(1)h, 1)s, 1)A) diffusion terms

Eij(x, y, t) strain rate tensor

T/(x,y,t) nonlinear shear viscosity

e 2 eccentricity of the elliptical yield curve

(Fx,Fy ) internal ice stress

f(x, y) Coriolis parameter

g 9.8ms- 2 acceleration of gravity

H Heaviside function

hi(X, y, t) ice thickness of ice-covered fraction

ho 1m ice cutoff thickness

hs(x, y, t) snow thickness On ice-covered fraction

M(x, y, t) ice mass (density times thickness)

P(x, y, t) ice pressure or strength

(P*, C) (2.75 x 104 ,20) ice strength parameters

(Sh, Ss, SA) thermodynamic terms

aij(x, y, t) stress tensor

i a air stress

Tw water stress

(u,v) the (x, y) components of ice velocity v

("VlQ, vw) 10 meter air and surface water velocities

(Pa, Pw) (1.3 kg m-3 , 1025 kg m-3 ) air and water densities

((x, y, t) nonlinear bulk viscosity

(w(x, y, t) height of the ocean surface

Table 3: Variables used in the ice momentum equations

86

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
8.2 Horizontal curvilinear coordinates

where the nonlinear viscosities are given by

(161)

(160)

(159)

(156)

(158)

(157)

(155)

(154)

(153)

O~A~l.

O~A~l.

87

if·r,=v

(

1J = 2"'

e

p = P* Ahi exp[-C(1 - A)]H(-V7. if) .

and

8A = _ 8(uA) _ 8(vA) _ Aa(A) V7. ifH(-V7. if) + SA + 1)A

8t 8x 8y

where a(A) is an arbitrary function such that a(O) = 0, a(l) = 1, and 0 a(A) 1. The ridging

term leads to an increase in hi under convergent flow as would be produced by ridging. The function

a(A) should be chosen so that it is near zero until the ice concentration is large enough that ridging

is expected to occur, then should increase smoothly to one.

The first two equations represent the conservation of ice and snow. Equation 158 is discussed in

some detail in MK89, but represents the advection of ice blocks in which no ridging occurs as long

as there is any open water. An optional ridging term can be added (Gray and Killworth [17]):

Applying the curvilinear transformation used in §2.5 and described in Appendix C, we use a transformation to an orthogonal curvilinear coordinate system. Denoting the velocity components in

the new coordinate system by

The "pressure gradient" term is also modeled as a term in the internal ice stress. This term

represents the resistance which ice has to being compressed (ice strength) and is a function of ice

thickness and concentration:

The Heaviside function guarantees that the ice has no strength when the flow is divergent (Gray

and Killworth [16]).

The second major component of the model consists of continuity equations describing the evolution of the ice thickness characteristics. Three parameters are calculated: the ice thickness hi,

the snow thickness hs , and the compactness, A, which is defined as the fraction of area covered by

thick ice. Note that Hibler's hI variable is equivalent to our Ahi combination-his hI is the average

thickness over the whole gridbox while our hi is the average thickness over the ice-covered fraction

of the gridbox. The continuity equations describing the evolution of these parameters (equations

(156)-(158)) also include thermodynamic terms (Sh, Ss and SA), which will be described in §8.5:

8Ahi = _ 8(uAhi) _ 8(vAhd + Sh + 1)h

8t 8x 8y

8Ahs __ 8(uAhs) _ 8(vAhs) S 1)

8t - 8x 8y + s + s

8A = _8(uA) _ 8(vA) +SA +1)A

8t 8x 8y

and

'.:.'.

•

;.'.

•

•

,:.•

•

:.'.

•

•

•

'.•

,.

•

•

•

'.".•

'.•

•

•

•

•

•

•

•

•

I.•

•

•

•

•

i.'

•
88

A more general expression for derivatives of tensors is given by Aris [3] in terms of Christoffel

symbols.

the equations of motion (144), (145), (156)-(158) can be re-written (see, e.g., Arakawa and Lamb

[2]):

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(170)

(169)

(168)

(167)

(166)

(165)

(164)

(163)

(162)

E12 = e2i = ·[m B(nv) n 0 (mu)]

2 n B~ m Bry .

lOU = m Bu + vmn~(~) B~ ory m

Bv B (1) E22 == n Bry + umn B~ ;;:

BA = -mn [~(AU) + ~(AV)] + SA

Bt B~ n Bry m

Bu B(w 1 (c c ) - = fv - gm- + - T~ + T~ +:Fe Bt B~ Maw

Bv o(w 1 - = -fu - gn- + - (TTl + TTl + :F ) Bt Bry Maw TI

BAhs = -mn [~(Ahsu) + ~(AhsV)]. + Ss

Bt . B~ n Bry m .

--BAhi -

_ -mn .[-B(AhiU) -- + -0 (AhiV).] -- + Sh

Bt B~ n Bry m

A [BT12 BT22 B ( 1 ) +ry m-- + n--- Tumn- - o~ Bry Bry m

A [BTu BT12 B ( 1) V' . T = m-- + n-- + Tumn- - B~ Bry B~ n

In curvilinear coordinates the divergence of a symmetric tensor Tis:

Sh, Ss and SA remain unchanged.

The viscous-plastic terms can be derived from equation (149). In curvilinear coordinates the

strain rate tensor can be written as:
m8P ---

2 8~

89

(172)

(171) 28(mu) 8 1 , 28v 8 1 2 2[]28 1

+2'TJn ' -'(~) - 2'TJmn - -(-) - 2'TJum n -'(-) 8'TJ 8'TJ m 8'TJ 8~ n 8~ n

+m-8 [ m 8(nv) n 8(mu)] 8 [8V 8 (1)] 'TJ---'-+'TJ---- +2n-- 'TJn-+'TJumn- - , 8~ n 8~ m 8'TJ 8'TJ 8'TJ 8~ n

2 8U8 (1) 22[8(1)]2 28(nv) 8(1) - 2'TJm n 8~ 8'TJ m - 2'TJvm n 8'TJ m + 2'TJm ---ar 8~ ;;,

:F. = n-8[(( - 'TJ)mn-"8 (U)] - + n-8 [(( - 'TJ)mn-8 (-

V)] - --n 8P

'T/ 8'TJ 8~ n 8'TJ 8'TJ m 2 8'TJ

+2'TJm2n-- 8u 8 (1) 2 2 8 ( 1) 8 (1) 28(nv) 8 ( 1 ) - +2'TJvm n - - - - +2'TJm ---- - 8~ 8~ n 8'TJ m 8~ n 8~ 8'TJ m

8.3 Numerical Scheme

The viscous-plastic terms become:

+2m_8 ['TJm _

8u +'TJVmn -

8 (~)] +n_8 ['TJ m _

8(_nV_) +'TJ!!.._8(~m_u-,-)] 8~ 8~ 8'TJ m 8'TJ n 8~ m 8'TJ

The coupled nonlinear equations of the model are treated as an initial value problem using energyconserving finite difference techniques. Under this initial value approach, values of all parameters

at all grid points are required to start the integration, and boundary values of velocity are required

thereafter. The forcing fields for the model consist of wind stress fields as well as surface winds

for the thermodynamic quantities, air temperature, dew point temperature, cloud cover, and snow

precipitation. It also needs near-surface ocean currents and sea-surface elevation for computing

the ocean stress term and the surface tilt term. The inertial terms are retained in the dynamical

equations, while the nonlinear advection is assumed to be small enough to neglect. Because of

the very strong ice interaction term, explicit integration of the momentum equations would force

timesteps to be extremely small (of the order of a second or less). The ice thickness equations,

on the other hand, can be explicitly integrated over much longer timesteps. To avoid this severe

timestep limitation due to the ice interaction, the momentum equations are integrated implicitly.

This impl~cit integration does, however, require a relaxation solution of a set of simultaneous

equations at each timestep. It is currently implemented using calls to either the ESSL library or

the NSPCG library. We are using NSPCG because it is freely available.

Using this implicit scheme, the inertial terms will normally be significant for short timesteps

(of the order of 1 hour or less) but insignificant for long timesteps. Although not a formal stability

requirement, it is wise to choose timesteps that are small compared to the variability of the ice

forcing. The effect of different timestep magnitudes on the momentum balance is discussed in

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

,e

•

•
Appendix B of Hibler [22). This appendix also contains some examples which may help the user

decide on timestep magnitudes for particular applications. Because of this implicit treatment of the

momentum equation, the only formal stability requirement is a Courant-Friedrichs-Lewy condition

for the advection terms in the thickness equations:

A key feature of the numerical scheme is a staggered grid, known as the "Arakawa C grid",

where the velocity is defined on the sides of a grid box and variables such as thickness and viscosity

are defined at the center of the grid box, as shown in Fig. 2.

Because ofthe strong ice interaction (which in this model is dissipative in nature) thl~ momentum

equations are essentially parabolic in form and hence have few numerical instability problems over

long-term integrations. (While there are few numerical problems, it should be emphasized that

the dissipative ice interaction terms are highly nonlinear and can lead to unstable flow fields in

the absence of water drag. Such a feature is a physical characteristic of plastic flow and not a

numerical artifact). However, in principle it is possible for the ice interaction to be very small even

though there may be a finite ice mass. Under this situation, the momentum advection terms could

cause nonlinear instabilities. To ensure against such situations, Hibler put a lower limit on the bulk

viscosity parameter (and hence indirectly shear viscosity as well). It is never allowed to drop below

1.0 x 107 kg/s, a value which negligibly modifies the ice drift.

Nonlinear instabilities over long-term integrations can also arise from the nonlinear advection

terms in the thickness contInuity equations. To avoid this problem, the advection terms in equations

(156)-(158) can be computed with a choice of either Smolarkiewicz's MPDATA ([53]) or third-order

upwind (Leonard [30]). See §3.7 for a description of these schemes. It may be possible to omit the

diffusivity from these equations if the forcing fields are sufficiently smooth.

8.4 Horizontal boundary conditions

As mentioned above, initial conditions at all points and ice velocities at the boundaries are thereafter

required to initiate the integration of the system of equations forward in time. The most natural

boundary condition is to take the ice velocity to be zero on the boundaries. This can be done either

at a land boundary or at an ocean location where there is no ice. Note that the boundary condition

does not affect the ice motion in such circumstances since in the absence of ice the strength is

zero. More generally, as long as the ocean boundaries are removed from the ice edge, the coupled

nature of the model will cause a natural ice edge boundary condition to be created. However, it

is also possible to form an "open" boundary condition by setting the strength equal to zero near

a boundary. These gridboxes are called "outflow cells" in Hibler [22]. In the Arctic simulations,

these outflow cells are used at the open edge near Greenland.

The boundary conditions on the momentum equations are to set u and v to zero at all boundaries, including islands. This is accomplished by the elliptic solver during the implicit timestep.

The ice and snow thickness and ice concentration equations have no-flux boundary conditions

imposed along the mask boundaries. The outflow cells contain a radiation condition if the velocity

is outward and no change if the velocity is inward.

The primary characteristic of the outflow cells is that the ice strength goes to zero there. The

values of P, (, and 'f} are all set to zero in outflow cells.

8.5 Thermodynamics

The thermodynamics used is based on the algorithm described in Mellor and Kantha [35], who

have a useful description of the various melting and freezing processes, plus the coupling to a full

90

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
91

(175:, (176:

(177:

(178)

(174)

(173)

O~A~l.

Fr

Qao

Wao 0

Wao < 0

q> = 4.0

q> = 0.5

""".-----

7-/ ,

I A A A

I

Wai I

1

I

I Wao I

WTO

hi -------------TI

""'-_--:L.- -+ Qio T

o

--

Fr

Figure 20: Diagram of internal ice temperatures and fluxes. The hashed layer is the snow.

Figure 19: Diagram of the different locations where ice melting and freezing can occur.

Figure 20 shows the locations of the ice and snow temperatures and the heat fluxes. The

temperature profile-is assumed to be linear between adjacent temperature points. The interior of

the ice contains "brine pockets", leading to a prognostic equation for the temperature TI .

The surface flux to the air is:

three-dimensional ocean. Their form of equations (156) imd (158) is:

O;:i + o(uo~hd +_O(Vo~hd = :~[A(WiO - Wai ) +(1:'" A)'Wao + WfT]

oA _o(uA) o(vA) = PoA [q>(1 _ A)W (1 - A)W ] ot + ox + oy Pihi ao + fT

Here, the W variables are the freeze or melt rates as shown in Fig. 19 and Table 4. The frazil ice

growth WfT will be discussed further in §8.5.2-note that it contributes to changes in A as well

as to changes in hi. The other term that contributes to A is Wao' This term includes a factor q>

which Mellor and Kantha set to different values depending on whether ice is melting or freezing:

'.•

•

•

•

•

•

•

•

•

:.'.

,e

•

•

•

•

•

•

•

•

•

•

•

•

•

'.•

'.•

•

:.•

•

".•

'.•

•

•

•

•

'.'-------------------------

• ------
Variable Value Description

Q'w 0.10 shortwave albedo of water

Q'i 0.50 shortwave albedo of ice

Q's 0.75 shortwave albedo of snow

Ck snow correction factor

Cpi 2093 J kg-I K-I specific heat of ice

Cpo 3990 J kg-I K-I specific heat of water·

lOw 0.97 longwave emissivity of water

lOi 0.97 longwave emissivity of ice

lOs 0.99 longwave emissivity of snow

E(T,r) enthalpy of the ice/brine system

FTt heat flux from the ocean into the ice

Ht sensible heat

~w fraction of the solar heating transmitted

through a lead into the water below

ki 2.04 W m-I K-1 therma,l conductivity of ice

ks 0.31 W m-I K-1 thermal conductivity of snow

Li 302 MJ m-3 latent heat of fusion of ice

L s 110 MJ m-3 latent heat of fusion of snow

LEt latent heat

LWt incoming longwave radiation

m -0.054°C/PSU coefficient in linear Tf(S) = mS equation

<I> contribution to A equation from freezing water

Qai heat flux out of the snow/ice surface

Qao heat flux out of the ocean surface

Qi2 heat flux up out of the ice

Qio heat flux up into the ice

Qs heat flux up through the snow

r brine fraction in ice

Pi 910 m 3 /kg density of ice

Si 5 PSU salinity of the ice

SWt incoming shortwave radiation

a 5.67 x 10-8 W m-2 K- 4 Stefan-Boltzmann constant

To temperature of the bottom of the ice

TI temperature of the interior of the ice

T2 temperature at the upper surface of the ice

T3 temperature at the upper surface of the snow

Tf freezing temperature

TmelLi mSi melting temperature of ice

TmelLs 0° C melting temperature of snow

Wai melt rate on the upper ice/snow surface

Wao freeze rate at the air/water interface

Wfr rate of frazil ice growth

Wio freeze rate at the ice/water interface

W ro W ai rate of run-off of surface melt water

Table 4: Variables used in the ice thermodynamics

92

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
(185)

(186)

(184)

(183)

(182)

(181)

(179)

(180)

E(T, r) == r(Li + CpoT) + (1 - r)CpiT

where

Inside the snow, we have'

Substituting into (184), we get:

These can be set equal to each other to solve for T2

The heat conduction in the upper part of the ice layer is

Substituting in for r and differentiating gives:

93

Note that in the absence of snow, Ck becomes zero and we recover the formula for the no-snow

case in which T3 = T2 ·

The formulas for sensible heat, latent heat, and incoming longwave and shortwave radiations are

the same as in Parkinson and Washington [41] and are shown in Appendix E. The sensible heat is

a function of T3 , as is the heat flux~tlit6tigh the snow Qs< Setting Qai = Qs, we can solve for T3

by setting T;+l = Tr + b.T3 and linearizing in b.T3. The temperature T3 is found by an iterative

solution of the surface heat flux balance (using the previous value of T1 in equation 186).. As in

Parkinson and Washington, if T3 is found to be above the melting temperature, it is set to Tmelt

and the extra energy goes into melting the snow or ice:

Note that L3 = (1 - r)Li plus a small sensible heat correction. We are not storing water on

the surface in melt pools, so everything melted at the surface is assumed to flow into the ocean

(Wro = Wad·

Inside the ice there are brine pockets in which there is salt water at the in situ freezing temperature. It is assumed that the ice has a uniform overall salinity of Si and that the freezing

. temperature is a linear function of salinity. The brine fraction r is given by

Sim

r=--

T1

The enthalpy of the combined ice/brine system is given by

•

•

•

•

•

•

•

'.•

•

•

:.•

•

•

•

:.•

:.

•

•

•

'.•

•

•

•

•

•

•

•

•

•

',.•

I.'.i.

•

'.•

:.•

•
Table 5: Ocean surface variables

8.5.1 Ocean surface boundary conditions

The ocean receives surface stresses from both the atmosphere and the ice, according to the ice

concentration:

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(192)

(191)

(190)

(189)

(188)

(187)

z-tO

94

FT -- = -CTz (To - T)

PoCpo

Variable Value Definition

b 3.0 factor

E evaporation

k 0.4 von Karman's constant

1/ 1.8 x 1O-6m28 - 1 kinematic viscosity of seawater

P precipitation

Pro 13.0 molecular Prandtl number

Prt 0.85 turbulent Prandtl number

So surface salinity

Tio stress on the ocean from the ice

Tao stress on the ocean from the wind

T internal ocean temperature

Ur friction velocity ITio1 1/ 2p;;-1/2

Zo roughness parameter

auw Ax I-Ax K m = ~Tio + --Tao

uZ Po Po

avw _ A y 1 - A y K m - -Tio + --Tao

uZ Po Po

where the relevant variables are in table 5.

The surface ocean is assumed to be at the freezing temperature for the surface salinity (To = mS)

where we use the salinity from the uppermost model point at z= -!~z. From this, we can obtain

a vertical temperature gradient for the upper ocean to use in the heat flux formula:

where

At the bottom of the ice, we have

2ki

QI0 = --,;;(To - T1)

The difference between Q10 and Q12 goes into the enthalpy of the ice:

Pi hi [~~ +v· V'E] = Q10 - Q12

We can use the chain rule to obtain an equation for timestepping T1:
95

We also want the final temperature and salinity to be on the freezing line, which we approximate

as:

(208)

(206)

(204)

(205)

(199)

(200)

(201)

(202)

(203)

(196)

(197)

(198)

(195)

(193)

(194)

== AWio + (1 - A)Wao

Tf =mS+nz.

m W1 == m W2 + mi

m W1 (CpwTl + L) = mW2(CpwT2 + L) + miCpiT2

m w1 S1 = m w2 S2.

C UT

T = z Prtk- 1 In( -zjzo) + BT

", _ (ZOUT ) 1/2 2/3 BT-b -- Pr

v

Once we have a the value for FT' we can use it to find the ice growth rates:

1

Wio = -L(Qio - FT) Po 0

1

Wao = -L(Qao - FT) Po 0

The ocean model receives the following heat and salt fluxes:

FT = AQio + (1 - A)Qao - WoLo

Fs = (Wo - AWro)(Si - So) + (1 - A)SoCP - E)Wo

We can then solve for r:

-T1 +mSl +nz , = (207)

+T1 (1 - ~) - mSl . cpw cpw

The ocean is checked at each depth k and at each timestep for supercooling. If the water is below

freezing, the temperature and salinity are adjusted as in equations (204) and (205) and the ice

above is thickened by the amount:

The variables are defined in Table 6. Defining, = miJmw2 and dropping terms of order ,2 leads

to:

where

where

8.5.2 Frazil ice formation

Following Steele et al. [55], we check to see if any of the ocean temperatures are below freezing at

the end of each timestep. If so, frazil ice is assumed to form, changing the local temperature and

salinity. The ice that forms is assumed to instantly float up to the surface and add to the ice layer

there. We assume balances in the mass, heat, and salt before and after the ice is formed:

•

•

'.•

•

•

•

•

•

•

•

\~•

,.•

•

•'.•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
Variable Value Definition

Cpi 1994 J kg- 1 K- 1 specific heat of ice

C pw 3987 J kg- 1 K- 1 specific heat of water

, mi/mw2 fraction of water that froze

L 3.16e5 J kg- 1 latent heat of fusion

mi mass of ice formed

m W1 mass of water before freezing

m W2 mass of water after freezing

m -0.0543 constant in freezing equation

n 7.59 x 10-4 constant in freezing equation

Sl salinity before freezing

S2 salinity after freezing

T1 temperature before freezing

T2 temperature after freezing

Table 6: Frazil ice variables

8.5.3 Differences from Mellor and Kantha

We have tried to modify the hakkis model to more closely follow MK89. However, there are also

ways in which we have deviated from it.

• Add advection of snow.

• Add lateral melting of snow when ice is melting laterally.

• We iterate on the solution of T3·

• We took a shortcut in the solution of So, To for the surfac,e heat and salt fluxes. We also

apply them differently to the ocean model.

• We added various limiters:

- Ice concentration:A 2: A min , Amin = 0.02.

- Ice thickness: hi 2: hmin, hmin = 0.1.

- Brine fraction: r r max , r max = 0.2

96

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

.'•

•

•

•

•

•

•

•

•

•

•
9 Description of the Ice Model and the Coupling

9.1 Ice model structure

The flow of the main program for the ice model is shown in Fig. 21. The overall structure essentially

consists of two components-the momentum equations and the ice continuity equations. The

momentum balance includes air and water stresses, Coriolis force, internal ice stress, inertial forces

and ocean tilt (equations (144) and (145)). There is a choice of rheologies (form of internal ice

stress) including viscous-plastic (used here), free drift, Mohr-Coulomb, and cavitating fluid. A semiimplicit timestep is done on the momentum equations which are solved by calling a solver from

the NSPCG library. The ice viscosities are non-linear functiOlis of the velocity, so the solution is

iterated several times, with the viscosities being recomputed each iteration..

The ice continuity consists of advection (done in ice_mpdata or ice_advect), and thermodynamics (run_mk and ice_mk). We have also added the diffusion of the ice tracer variables in order

to obtain a smoother solution. This may be optional if the forcing fields are sufficiently smooth.

The main program calls ice-step, which in turn calls:

ice_adved Compute the ice advection according to a third-order upwind advection scheme. The

ice ridging and ice diffusion also happen here if they are enabled.

ice_bctrans Boundary conditions for A, hi and hs.

ice_bcuv Boundary conditions for ice u and v.

ice_cavitating Compute the ice viscosities according to a cavitating fluid rheology.

ice_freedrift Compute the ice parameters to model free drift (no internal ice stress contribution).

ice_gencoef Generate coefficients for the iterative solution to the ice momentum equations.

ice_mohrcoulomb Compute the ice viscosities according to a Mohr-Coulomb fluid rheology.

ice_mpdata Compute the advection of a scalar according to a Smolarkiewicz advection scheme.

iceJ"hs Gather the right-hand-side terms for the ice momentum equations prior to using the

solver.

ice-solver~SPCG Solve the implicit ice momentum equations using the NSPCG library.

run_mk Compute the changes in ice thickness and concentration due to the thermodynamics.

Also produce the surface heat and salt forcing for the ocean model.

ice_viscplast Compute the ice viscosities according to a viscous-plastic fluid rheology.

9.1.1 Thermodynamic subroutines

The thermodynamic subroutines used in the model are:

ice_mk Compute the net energy budget and change of thickness for ice and snow.

ice_frazil Compute the frazil ice growth due to supercooled water in the ocean (called from

main3d after the ocean timestep).

rads Compute the heat flux budgets over the water.

97
98

."•••••

.,•••••••••••••••••••••••••••••••••.1•

.;•

ice-freedrift

iterate3 times

MPDATA 3rd-order upwind,

Figure 21: Flow chart for the sea-ice model.

ice-lIlpdata

ice_viscplast
•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

9.1.2 Initialization

During initialization, the following routines are called:

deLice_avg Create the netCDF ice 'averages file.

defjce_his Create the netCDF ice history file.

defjce-rst Create the netCDF ice restart file.

freeze Make sure that no water temperatures are below freezing during initialization. It does

not form frazil ice.

getjce Read a record from the netCDF ice restart file.

hakkblkdat Block data initializing some parameters for the ice thermodynamics.

ice_blkdat Block data initializing some parameters for the Smolarkiewicz advection routine.

ice_init Initialize the ice model, either by reading a history file or by setting the initial values.

ice_user! Initialize some ice variables.

iniLhakkis Initialize some ice thermodynamic variables.

initjce Initialize some more ice variables.

9.1.3 Forcing fields

The ice model requires some extra forcing fields. These fields require their own routines for reading

them from the netCDF forcing files:

geLairp Reads surface air pressure from the forcing NetCDF file, and then linearly time-interpolates

to current model time.

geLairt Reads surface air temperature from the forcing NetCDF file, and then linearly timeinterpolates to current model time.

geLcloud Reads cloud fraction from the forcing NetCDF file, and then linearly time-interpolates

to current model time.

geLdewt Reads surface dewpoint temperature from the forcing NetCDF file, and then linearly

time-interpolates to current model time.

getJ.rflux Reads incoming longwave radiation from the forcing NetCDF file, and then linearly

time-interpolates to current model time.

geLprecip Reads precipitation from the forcing NetCDF file, and then linearly time-interpolates

to current model time.

geLwind Reads surface winds from the forcing NetCDF file, and then linearly time-interpolates

to current model time.

99
9.1.4 Other subroutines

The other subroutines in the ice model include:

ice_forcing Computes the forcing fields due to exchange of momentum by the ice and ocean.

smooth Smooths a field with a five-point Laplacian filter.

wrtjce_avg Writes a record to the netCDF ice averages file.

wrtjce_his Writes a record to the netCDF ice history file.

wrtjce-.rst Writes a record to the netCDF ice restart file.

9.2 Coupling strategy

The flow chart for the coupled ice-ocean model is shown in Fig. 22. The ice model is stepped

first since it provides the surface heat and salt fluxes to the ocean model. After the ocean model

is timestepped, the ocean temperatures are checked for supercooled water-if any is found it is

converted into frazil ice and the ocean temperature and salinity are adjusted to be at freezing. The

ocean equation of state is computed after this conversion.

9.3 C preprocessor variables

The ice model has two C preprocessor variables in cppdefs.h:

ICE Define to use ice component of the model.

ICE_THERMO Define for ice thermodynamics.

The ice model requires new forcing fields to be read or computed:

ANA-AIRT Define for an analytic air temperature.

ANA_CLOUD Define for an analytic cloud fraction.

ANA...DEWT Define for an analytic dew-point temperature.

ANA--LRFLUX Define for an analytic incoming longwave radiation.

ANA_SNOW Define for an analytic snow fall rate.

ANA_SLP Define for an analytic sea-level pressure.

ANA_WIND Define for analytic wind fields.

There are also some C preprocessor variables in icedefs.h:

NSPCG Use the nspcg library for the implicit solver. Either this or ESSL must be turned on.

ESSL Use the essllibrary for the implicit solver.

ANA...ICE.JNIT Initialize the ice fields analytically as opposed to reading them from a file.

ICE.-MOMENTUM Compute the ice momentum equations.

ICE-ADVECT Advect the ice tracers.

ice_GSCHEME Use a third-order upwind advection scheme for the tracers. This must be defined

for either of these to take effect:

ICE...DIFFUSION Add a Laplacian diffusion on the ice tracers.

ICE-RIDGING Add an ice ridging scheme from Gray and Killworth [16].

100

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
'.'.,e•••

:.'.••

I:.•

'.•

'.'.•

'.'.•

'.,:.•

:.'.•

'.;.••••••••••••••••

Figure 22: Flow chart for the coupled ice-ocean model.

101
102

•••••••••••••••••••••••••••••••••••••••••••

e:
103

A Model Timestep

where ¢ represents one of u, v, T, S or ( and F(t) represents all the right-hand-side terms. The

simplest approximation is the Euler timestep:

where the coefficients a, (3 and, are chosen to obtain a third-order estimate of ¢(t + b.t). We use

a Taylor series expansion:

(220)

(218)

(219)

(215)

(216)

(217)

(214)

(213)

(211)

(210)

(209)

(3 =

a =

, =

a¢(t) = F(t)

at

= ¢'

b.t2 ¢' - b.t¢" + __¢III + ... 2

= ¢' - 2b.t¢" + 2b.t2¢III + ...

¢(t + b.t) - ¢(t) = F(t)

b.t

¢(t + b.~~ - ¢(t) = aF(t) + (3F(t _ b.t) + ,F(t - 2b.t)

¢(t + b.t) - ¢(t) = ¢' + b.t ¢" + b.t2¢III + ... b.t 2 6

F(t)

F(t - b.t)

F(t - 2b.t)

We find that the coefficients are:

23

12

4

3

5

12

The model carries one time level for the physical fields and three time levels of the right-handside information. The initial fields are read in but the right-hand-sides are not stored; an Euler

timestep is used for the first two steps to get things going.

where

¢(t + b.~~ - ¢(t) = [F(t) + F*(t + b.t)]. (212)

This leapfrog-trapezoidal timestep is stable with respect to diffusion and it strongly damps the

computational mode. However, the right-hand-side terms are computed twice per timestep.

The timestep on SCRUM's full 3-D fields is dOhe with a third-order Adams-Bashforth step. It

uses three time-levels of the right-hand-side terms:

This timestep is more accurate, but it is unconditionally unstable with respect to diffusion. Also,

the even and odd timesteps tend to diverge in a computational mode. This computational mode

can be damped by taking correction. steps. SCRUM's timestep on the depth-integrated equations

is a leapfrog step with a trapezoidal correction on every step, which uses a leapfrog step to obtain

an initial guess of ¢(t + b.t). We will call the right-hand-side terms calculated from this initial

guess F*(t + b.t):

where you predict the next ¢ value based only on the current fields. This method is accurate to

first order in b.t; however, it is unconditionally unstable with respect to advection.

The leapfrog timestep is accurate to O(b.t2 ):

¢(t + b.t) - ¢(t - b.t) = 'T"( ) 2b.t .r t .

Numerical timestepping uses a discrete approximation to:

;.•

•

Ie•

•

,e•

•

•

•.'.'.'.

ie•

;.I.'.Ie

••

e

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
104

••••••••••••••••••••••••••••••••••••••••••••
105

• Some problems turn out to be sensitive to the value of 0 used.

The derivative of C(s) can be computed analytically:

(225)

(224)

(223)

z = ((1 + s) + hes + (h - he)C(s), ~l::;s::;O (221)

z = (( + h)(l + s) - h

which is the a-coordinate.

• For b = 0, the resolution all goes to the surface as 0 is increased.

• It has a linear dependence on ( and is infinitely differentiable in s.

• The larger the value of 0, the more resolution is kept above he·

• For b = 1, the resolution goes to both the surface and the bottom equally as 0 is increased.

C(s) = (1 _ b)sinh(Os) + btanh[O(s + !)] - tanh(!O) (222)

sinhO 2tanh(!0)

• It is a generalization of the a-coordinate system. Letting 0 go to zero and using L'Hopital's

rule, we get:

However, we choose to compute Hz discretely as b.z/b.s since this leads to the vertical sum of Hz

being exactly the total water depth D.

• For 0 =1= 0 there is a subtle mismatch in the discretization of the model equations, for instance

in the horizontal viscosity term. We recommend that you stick with "reasonable" values of

0, say 0 ::; 5.

Figure 23 shows the s-surfaces for several values of 0 and b for one of our domains. It was produced

by a Matlab tool written by Hernan'Arango which is available from our web site (see §l.1).

We find it convenient to define:

where 0 and b are surface and bottom control parameters. Their ranges are 0 < 0 < 20 and

o::; b ::; 1, respectively. Equation (221) leads to z = ( for s = 0 and z = h for s = -l.

Some features of this coordinate system:

Following Song and Haidvogel [54], the vertical coordinate has been chosen to be:

where he is either the minimum depth or a shallower depth above which we wish to have more

resolution. C(s) is defined as:

B The vertical s-coordinate

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

120

120

100

100

60 80

XI (god un~s)

60 80

,XI (god un~s)

Gr~ (W-po~ts) Section at ETA =66

Gr~ (W-points) Sectional ETA=66

40

40

20

20

-4500

-500

-5000

-3500

-4000

-1500

-2500

-2000

-1000

-3000

120

120

100

100

106

60 80

XI (god units)

60 80

XI (grid un~s)

Gr~ (W-points) Section at ETA =66

Gr~ (W-points) Section at ETA =66

40

40

20

20

Figure 23: The s-surfaces for the North Atlantic with (a) () = ·0.0001 and b = 0, (b) () = 8 and

b = 0, (c) () = 8 and b = 1. (d) The actual values used in this domain were () = 5 and b = 0.4.
(231)

(230)

(229)

(228)

(227)

(226)

h h k m n

V' x a= mn 0 0 0

of, °Tl OZ

2:.. 12. c m n

107

2 [ 0 (m o¢) 0 (n o¢)] \7 ¢ = V' . V'¢ = mn o( n o( + 07] m 07]

where ¢ is a scalar and ais a vector with components a, b, and c.

Here, m((,7]) and n((,7]) are the stale factors which relate the differential distances (.6.(, .6.7]) to

the actual (physical) arc lengths.

It is helpful to write the equations in vector notation and to use the formulas for div, grad, and

curl in curvilinear coordinates (see Batchelor, Appendix 2, [4]):

C Horizontal curvilinear coordinates

The requirement for a boundary-follQwJng, toordinate;systerp.-l1nd for a laterally variable grid resolution can both be met (for suitably smooth domains) by introducing an appropriate orthogonal

coordinate transformation in the hoHiontal. Let the new cdof'Clinates be ((x, y) and 7](x, y) where

the relationship of horizontal arc length to the differential distance is given by:

r•

.'

".'.'.'.'.•

•

•

'.'.'.•

•

•

•

•,'.

••'.".•

•

'.•

•

•

•

•

•

•

•

'.'.•

I:.•

•

•

•

•

•

•
108

-- .

••••••••••••

••••••••••••••••••••••••••••••••
109

D.2 Horizontal Diffusion

(233)

(232)

We have chosen anything from zero to the value of the horizontal viscosity for the horizontal

diffusion coefficient. One common choice is an order of magnitude smaller than the viscosity.

D.3 Vertical Viscosity and Diffusion

SeRUM stores the vertical mixing coefficients in arrays with three spatial dimensions called Akv

and Akt. Akt also has a fourth dimension specifying which tracer, so that temperature and salt

can have differing values. Both Akt and Akv are stored at w-points in the model; horizontal

averaging is done to obtain Akv at the horizontal u and v-points. The values for these coefficients

can be set in a number of ways, as described in §3.11.

for the Laplacian and biharmonic viscosity, respectively. We have found that the model

typically requires the boundary layer to be resolved with at least one grid cell. This leads

to coarse grids requiring large values of II.

This time should be short enough to damp out the numerical noise which is being generated but long enough On the larger scales to retain the features you are interested in.

This time should also be resolved by the model timestep.

boundary layer thickness The western boundary layer has a thickness proportional to

D Viscosity and Diffusion

D.l Horizontal viscosity

The horizontal viscosity and diffusion coefficients are scalars which are read in from serum.in.

Several factors to consider when choosing these values are:

spindown time The spindown time on wavenumber k is k'I'=1 for the Laplacian operator and ~l

for the biharmonic operator. The smallest wavenumber corresponds to the length 2~x

and is k = ;x' leading to

~·W ;.

••

•

'.'.;'.•

:.•

I.

•

;.••

,e:.

:.•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
110

.':;' ••••••••••••••••••••••••••••••.:.'

1

•.1•••••••••
E.! Shortwave radiation

where the variables are as in Table 7. The cosine of the zenith angle is computed using the formula:

(241)

(242)

(243)

(240)

(239)

(238)

(237)

(236)

(235)

(234)

111

611 x lOa(Tr273.16)/(Trb)

611 x lOa(Ts!c-273.16)/(Ts!c-b)

e

Q _' Scos2 Z

0- (cos Z + 2.7)e x 10-5 + 1.085 cos Z + 0.10

0= 23.44° x cos [(172 - day of year) x 271"/365]

LW.!-= (1 + 0.275c) F.!-.

The latent heat depends on the vapor pressure and the saturation vapor pressure given by:

E.4 Latent heat

cos Z = sin¢sino + cos ¢cos 0cos HA.

HA = (12 hours - solar time) x 71"/12.

The sensible heat is given by the standard aerodynamic formula:

E.3 Sensible heat

while the cloud correction is given by:

The clear sky formula for incoming longwave radiation is given by:

E.2 Longwave radiation

An estimate of the cloud fraction c will be provided by Jennifer Francis ([12]).

The correction for cloudiness is given by

and the hour angle is

The declination is

The Zillman equation for radiation under cloudless skies is:

As was seen in §8.5, the model thermodynamics requires fluxes of latent and sensible heat and long~

wave and shortwave radiation. We follow the lead of Parkinson and Washington [41] in computing

these terms.

E Radiant heat fluxes

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

'.•

!.••

:.•

'.•

'.•

•

•

'.'.•

•

•

e

e

.e

•

'.,e

•

e

e

e

e

'.
112

Table 7: Variables used in computing the incoming radiation and latent and sensible peat

Note that these need to be computed independently for the ice-covered and ice-free portions of each

gridbox since the empirical factors a and b and the factor L differ depending on the surface type.

•

•

•

•

•

•

.1•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

(245)

(244)

(246)

p-(l-E)e

Ees

p- (1- E)es

Variable Value Description

(a, b) (9.5, 7.66) vapor pressure constants over ice

(a, b) (7.5, 35.86) vapor pressure constants over water

c cloud cover fraction

CE 1.75 x 10-3 transfer coefficient for latent heat

Cn 1.75 x 10-3 transfer coefficient for sensible heat

cp 1004 J kg- 1 K- 1 specific heat of dry air

8 declination

e vapor pressure in pascals

es saturation vapor pressure

E 0.622 ratio of molecular weight of water to dry air

HA hour angle

L 2.5 x 106 J kg- 1 latent heat of vaporization

L 2.834 x 106 J kg- 1 latent heat of sublimation

</> latitude

Qo incoming radiation for cloudless skies

qs surface specific humidity

QlOm 10 meter specific humidity

Pa air density

S 1353 W m- 2 solar constant

a 5.67 x 10-8 W m- 2 K- 4 Stefan-Boltzmann constant

Ta air temperature

Td dew point temperature

TSfc surface temperature of the water/ice/snow

Vwg geostrophic wind speed

Z solar zenith angle

The latent heat is also given by a standard aerodynamic formula:

The vapor pressures are used to compute specific humidities according to:
113

F.1 File inclusion

F.2 Macro substitution

(.5 * «al) + (a2)))

name replacement text

#define av2(al,a2)

#define

A macro definition has the form

where name would be replaced with "replacement text" throughout the rest of the file. This is

used in SCRUM as a reasonably portable way to get 64-bit precision:

although this is riskier than the equivalent statement function

av2(al,a2) = .5 * (al + a2)

The statement function is preferable because it allows the compiler to do type checking and because

you don't have to be as careful about using enough parentheses.

The third form of macro has nO replacement text at all:

#define MASKING

In this case, MASKING will evaluate to true in the conditional tests described below.

We are using the C preprocessor style of inclu~e because many of the SCRUM include files are not

pure Fortran and must be processed by epp.

Placing common blocks in smaller files, which are then included in each subroutine, is the easiest

way to make sure that the common blocks are declared consistently. Many Fortran compilers

support an include statement where the compiler replaces the line

include 'file.h'

It is customary to use uppercase for epp macros-the C preprocessor is case sensitive.

It is also possible to define macros with arguments, as in

#define BIGREAL

with the contents of file.h; file.h is assumed to be in the current directory. The C preprocessor

has an equivalent include statement:

#include "file.h"

The C preprocessor, epp, is a standalone program which comes with most C compilers. On many

UNIX systems it is not in the default path, but in /lib or in / usr/lib. If you do not have a C

preprocessor then there are several v~rslons available via anonymous ftp. For instance, ftp.uu.net

has two in the /published/oreilly/nutshell/imake directory-I have built and used the one

from Der Mouse on a Cray. I have put this one in pubjutiljcpp.tar.gz on the ahab.rutgers.edu ftp

site since it supports the #elif construct. One also comes with gee, the gnu C compiler. If you

build this compiler, epp will have a path such as

/usr/local/lib/gcc-lib/sparc-sun-solaris2.5/2.7.2/cpp

where spare is the architecture, sun is the manufacturer; solaris2.5 is the operating system and

version, and 2.7.2 is gee's version number.

This Appendix describes the C preprocessor as used in SCRUM with the Fortran language. A

more complete description is given by Kernighan and Ritchie [26]. More practical advice on using

epp is given by Hazard [21].

FThe C preprocessor

'.,e

e

•

•

•

•

•

I-I.

•

-'.•

•

•

i.•

'.'.•

•

•

•

I'.!'.

I.

•

•

•

•

•

•

...

•

•

•

•

•

•

•

•

•

-
F.3 Conditional inclusion

It is possible to control which parts of the code are seen by the Fortran corripiler by the use of

cpp's conditional inclusion. For example, the statements

#ifdef MASKING

# include "mask.h Ol

#endif /* MASKING */

# ifdef MASKING

c

c Apply Land/Sea mask: slipperiness.

c

do j=1,M

do i=2,Lm

Uflux(i,j)=Uflux(i,j)*pmask(i,j)

enddo

enddo

# endif /* MASKING */

will not be in the Fortran source code if MASKING has not been defined. Likewise, #ifndef'

tests for a maCro being undefined:

#ifndef RMDOCINC

c rmask Mask at RHO-points (O=Land, 1=Sea).

c pmask Slipperiness mask at PSI-points (O=Land, 1=Sea,

c 1-gamma2=boundary).

c umask Mask at U-points (O=Land, 1=Sea).

c vmask Mask at V-points (O=Land, 1=Sea).

c

c=======================================================================

#endif

There are also #else and #elif (else if) statements, although #elif is newer and is not supported by all versions of cpp. An example using #else and #elif is shown:

#if defined BASIN

parameter (L=181, M=141, N=12, NT=1)

#elif defined CANYON_A

parameter (L=66, M=49 , N=10, NT=1)

#elif defined CANYON_B

parameter (L=66 , M=49 , N=15, NT=1)

#elif defined UPWELLING

parameter (L=42 , M=81, N=16, NT=2)

#else

parameter (L=???, M=???, N=??, NT=?)

#endif

Actually, #ifdef is a restricted version of the more general test

#if expression

where "expression" is a constant integer value. Zero evaluates to false and everything else is

considered true. Compound expressions may be built using the C logical operators:

114

.<J

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
i.

'.:.•

i.

•

•

•

'.•

•

•

•

•

,e

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

••

•

•

•

•

•

•

•

•

•

&& logical and

II logical or

logical not

These symbols would be used as in the following example:

#if defined CANYON_A I I defined CANYON_B

do j=O,M

do i=O,L

yc=c32000-c16000*(sin(pi*xr(i,j)/xl))**24

h(i,j)=c20+p5*(hmax-c20)*(cl+tanh«yr(i,j)-yc)/cl0000))

enddo

enddo

#endif

F.4 C comments

The C preprocessor will also delete C language comments starting with /* and ending with */ as

in:

#endif 1* MASKING *1

When mixed with Fortran code, it is safer to use a Fortran comment.

F.5 Potential problems

The use of the C preprocessor is not entirely free of problems, but many can be worked around or

avoided by using the Der Mouse version of epp.

1. Apostrophes in Fortran comments. epp does not know that it is in a comment and some

versions will complain about unmatched apostrophes in the following:

c Some useful comment about Green's functions.

The gnu version of epp (which comes with gee) has a -traditional option which makes it

more appropriate for use with Fortran.

2. C++ comments. Some of the newer versions of epp will remove C++ comments which go from

'//' to the end of the line. Some perfectly reasonable Fortran lines contain two consecutive

slashes, such as:

common II varl, var2

44 format U I)

and the new Fortran 90 string concatenation:

mystring = 'Hello, ' II 'World!'

3. Macro replacement. One feature of epp is that you can define macroS and have it perform

replacements. The code:

115
I started working on these ocean models before 1990, much less before Fortran 90 compilers were

generally available. Fortran 90's kind feature would be a better way to handle the BIGREAL type

declarations. On the other hand, Fortran 90 does not include conditional compilation. However,

it is deemed useful enough that the Fortran 2000 committee has a draft document describing how

Fortran might support conditional compilation. We might start using this in about ten years.

116

double precision really_long_variable, second_long_variable

#define REAL double prec~s~on

REAL really_long_variable, second_long_variable

Modern Fortran

and you run the risk of creating lines which are longer than 72 characters in length.

•

•

•

•

•

•

•

•

•

Also, make sure that your macros will not be found anywhere else in the code. I used to use •

#define DOUBLE for double precision until it was pointed out to me that DOUBLE PRECISION •

is perfectly valid Fortran. The macro processor would turn this into 1 PRECISION since •

something that is defined has a value of 1. •

•

•

•

•

•

•

•

•

•

•

.'•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

becomes

F.6
•

•

•

'.•

!.•

•

•

•

•

•

;

•

•

•

•

•

'.'.•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

G The patch program

We sometimes discover things in SCRUM which we would like to modify, either to fix bugs, or

to add new features. Hernan Arango keeps track of these changes and periodically sends patches

to the list of known SCRUM users so that they can update their versions. By sending out these

changes rather than the whole updated model, people can acquire bug fixes and still retain the

changes they have made to SCRUM for their own applications.

Larry Wall has written a program to take the output of cliff and automatically apply it to the

old version of a file to create the new version. This program is called patch and is available from

all the gnu archive sites. If the output of cliff has been saved in the file scrum.patch.20 then

patch would be used as follows:

patch < scrurn.patch.20

As patch updates the files, it leaves the original of file in file.orig. If it gets confused for some

reason (if you modified the lines of code patch wants to change) it will create a file.rej file. I often

check to see if any .rej files get created-these can be used to patch file by hand and can then be

deleted.

117
118

••••••.i•••••••.i•••••••••••••••••••••••••••••
119

H.I imake

H Makefiles

macros These are lines of the for:m eFT = £77, which in this case allows you to use one name

($(CFT» for the compiler even though each co:mpiler has a different name. The macros

in the Makefiles are defined in two separate sections:

ezgrid: ezgrid.o

<TAB>$(CFT) -0 ezgrid ezgrid.o

These are lines of the form:

where ezgrid is the target to be compiled, and ezgrid depends on ezgrid.o. make will

first check to ensure that ezgrid.o is up to date and then execute the commands on the

following lines (which must start with a <TAB> character).

machine dependent These macros give the name of the compiler and sensible flags for

that compiler, etc.

project dependent These macros depend on the project but not the computer, such as

the list of source files used to build the executable.

rules

dependencies These lines tell make which object files must be rebuilt when an include file is

modified. Also, if the C processor is being invoked specifically by the suffix rule for .F.0,

creating an intermediate .f file, then make must be told to recreate the .f file after its

include files are modified. The dependency lines are generated automatically by a perl

script, which searches the source files for #include directives (see Appendix I). Note

that if you add your own source files with #include statements, you will need to rerun

make depend to update the dependency list. If your files are not in the list of SCRUM

sources, they will have to ,be added to the depend: entry in the Makefilejlmakefile

first.

Since it is difficult to keep consistent Makefiles for several different computers, it was suggested

that we try imake, which is distributed with the X window system. It helps you to separate the

system dependent parts of a Makefile into configuration files (kept in a central location) and project

dependent parts called Imakefiles. Then, when you want to make SCRUM on a Cray computer,

One of the software development tools which comes with the UNIX operating system is called

make. make has many uses, but is most co:mmonly used to keep track of how a large program

should be compiled. You provide it with a list of your source files and instructions on how to

compile them. It will check the relative ages of the source and object files, only compiling those for

which the object file is out of date. It is assumed that make will be used to compile SCRUM and

its related programs. See Oram and Talbott [38] for a description of make that is easier to read

than the man page.

The file in which you provide make with its commands is usually called Makefile. Several

Makefiles are provided with SCRUM, one for each brand of computer to which I have easy access.

These Makefiles have become quite complex, but are organized into several sections:

suffix rules These are lines of the form .F.o:, followed by a rule telling make how to make a file

called foo.o from foo.F. This particular rule is used extensively and comes in two forms,

depending on whether or not the compiler will invoke the C preprocessor (cpp) for you.

If the compiler does not invoke cpp then make will do so, creating an intermediate foo.f

file.

'.

Ie'.

•

•

•ie

•

;.•

,.

•

•

•

•

'.•

•

•

•

•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
you combine the SCRUM Imakefile with the Cray configuration file to create the appropriate

Makefile. This is done by the shell script fmkmf, which takes the computer type as its argument:

fmkmf Cray

This will generate Makefile.Cray. If you provide an unrecognized computer type then the

generie.ef file will be used. The list of recognized computers is growing and currently includes:

Alpha DEC Alpha running OSF/1.

CM2 Connection Machine.

CF90 Cray with £90 and UNICOS.

Cray Cray with ef77 and UNICOS.

F90 The NAG Fortran 90 compiler in free format style.

HP Hewlett-Packard 9000/700 family.

Gnu Gnu Fortran.

NAG The NAG Fortran 90 compiler in fixed (old) format style.

RS6000 IBM RS/6000 with xlf and AIX.

RS60000ld IBM RS/6000 with an older xlf and AIX.

SGI SGI with IRIX and £17.

Solaris Sun Sparc with Solaris 2.x.

Sun Sun Spare with SunOS 4.x.

Titan Kubota Titan 3000 with fe.

It would be possible to have a different configuration file for each Cray you use, or for different

versions of the Sun compiler.

The fmkmf script executes imake followed by make depend, as shown in Fig. 24. It requires the $CONFIGDIR variable to be set to where the configuration files are kept. Also, the

make depend phase executes a perl script that requires the perl program. The configuration

files and the fmkmf script are distributed as described in §1.1 and perl is available from all the

gnu archive sites.

We could have chosen to have the fmkmf script try to determine which type of computer it

was being run on, and use the appropriate configuration file. However, it may be that we decide

to run SCRUM on a computer which does not have imake or perl. This way, we can generate

the Makefile on a computer which has the necessary support programs and then just transfer the

Makefile along with the SCRUMcode.

H.2 Your Makefile

If you are using one of the environments for which Makefiles are supplied, you are set, although

you may want to check the compile flags. Otherwise, you have the choice of copying and modifying

an existing Makefile or using imake and creating your own configuration file. In either case, you

will need to know certain things about your environment such as:

• The name of the Fortran compiler.

120

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
ie

•I.

Ie

•

I.•

'.•

•

•

•

•

'.•

•

•

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

$CONFIGDIR/

SCRUM.tmpl I--

imake make depend

$CONFIGDIR/ Makefile without Makefile with

machine.cf dependencies dependencies

Imakefile -

Figure 24: Creating Makefiles

• The compiler options you wish to use.

• How to link to the NCAR graphics libraries, ifthey exist and you wish to use them.

• Whether or not the Fortran compiler will invoke cpp and how to tell either the compiler or

make to do so.

• What file extension the compiler requires.

The biggest changes to the Makefile result from the way cpp is executed. If your compiler does

it for you, start from the Sun files, otherwise start from the old IBM RS/6000 files. Sometimes it

is better to tell make to execute cpp even if the compiler will do it. For instance, the old Cray

debugger became confused about line numbers if it did not have the intermediate .f files to work

with.

You will also have to edit the Makefile or Imakefile if you add source files to SCRUM. In this

case, you will have to add the new files to the OBJS and SReS macros, and make sure that the

dependencies are listed correctly for the new files.

121
122

•••••••••••••••-I.1

••••••••••••••••••••••-••••
~.

:.•

'.

Ie'.i.:.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

I Perl scripts for Fortran

Perl is a computer language, invented by Larry Wall, for manipulating text and other useful things.

It is fully described in Wall et al. [60], while a more tutorial approach is given by Schwartz [49].

Perl itself is available from your nearest CPAN archive site, for instance:

http://www.perl.com/CPAN/

I have several Perl scripts which I find useful when working with Fortran programs, and which

are available from:

http://marine.rutgers.edu/po/perl.html

It is not necessary to know Perl to use these scripts, but it must be installed on your system. To use

these scripts, simply place them somewhere in your path and make sure that they are executable:

chmod 755 relabel

(on a UNIX machine).

The following scripts modify your source code and usually work on the style of Fortran in

SeRUM, but have been known to do the wrong thing. Some of the scripts become confused when

part of an if or do statement is inside an #ifdef clause. The following will parse as two nested do .

loops, only one of which is terminated:

#ifdef EW_PERIODIC

do i=l,Lm

#else

do i=O,L

#endif

enddo

It is also extremely dangerous to run relabel on an arithmetic if.

Do not delete your original code before checking the new code.

1.1 redo

This program reformats do loops and was written to convert

do 10 i=1,20

10 sum = sum+i

to

do 10 i=1,20

sum = sum+i

10 continue

The -E option tells it to use enddo instead of continue as in

do i=1,20

sum =: sum+i

enddo

redo is used as follows:

redo < file.F > file.new

redo was written so that findent would work on SPEM.

123
1.2 findent

findent will indent your Fortran code two spaces for do loops and ifstatements. It will not correct

lines which extend beyond 72 characters, but will print out a warning for each one. It assumes that

each do loop ends with a continue or enddo. findent is used as follows:

findent < file.F > file.new

There is an option to change the number of spaces for each level ofindenting. To get an indent of

four instead of two, use:

findent -n 4 < file.F > file.new

See the comments at the top of the code for the more obscure options.

1.3 relabel

relabel was written by Sverre Froyen to replace the numbered Fortran labels with new sequentially

ordered labels. It was the first of these scripts and helped me to write the rest. relabel is used as

follows:

relabel < file.F > file.new

It does have some known bugs, however:

• No computed goto.

• No assigned goto or assign.

• No arithmetic if.

• No new-lines inside the parenthesis immediately following a read/write.

• Others not yet discovered.

All the source files in SPEM have been run through redo, findent, and relabel. In the C shell

(csh), a series of files can be processed at once:

ahab% foreach file (*.F)

foreach? redo < $file > $file.red

foreach? findent < $file.red > $file.fin

foreach? relabel < $file.fin >$file.rel

foreach? echo $file done

foreach? end

You can then rename $fi1e.rel to $file.F and get rid of the temporary files after checking to make

sure that all the files still look sensible.

1.4 unenddo

unenddo will turn all do-enddo loops into do-continue loops to comply with the Fortran 77

standard. It is used as follows:

unenddo < file.F > file.new

unenddo replaces enddo statements with labelled continue statements. It starts numbering

these statements at 2000 assuming that existing labels use only three digits. If desired, unenddo

can be told to start labelling with a different number by modifying the $labeLno-start variable.

124

.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
r.~.. . ,

•

".;.

•

'.•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

1.5 ifspace

When I am feeling particularly contentious I also run ifspace' on the code. It will convert

if(i.eq.O.or.j.eq.O) then

to

if (i .eq. 0 .or. j .eq. 0) then

Use as follows:

ifspace < file.F > file.new

1.6 sfmakedepend

The other Perl script I use with Fortran modifies the Makefile to include dependency information,

much like the Xll program makedepend. I originally wrote fmakedepend which was used with

traditional Fortran include statements. I later wrote a variant of it for use with the C preprocessor,

called sfmakedepend. The latest version of sfmakedepend does the job of both programs and

also searches for the dependencies introduced by Fortran 90 modules. It is used by the Makefiles

described in §H.

It recursively searches for Fortran style includes, for instance if file.f has the statement:

include 'commons.h'

the line

file.o: commons.h

will be added to the bottom of the Makefile. This tells make that file.o depends on commons.h

as well as file.f, and to recompile file.f whenever commons.h is modified. It.likewise searches

source files for C style includes such as

#include "commons.h"

and adds the corresponding dependencies to the Makefile. It has several options, including -s,

required for Fortran compilers which will not invoke the C preprocessor for you. In this case the

above dependency line would become

file.o: commons.h

file.f: commons.h

letting make know that the C preprocessor must be rerun on file.F whenever commons.h is

updated.

When using the C preprocessor, you can ask it to search directories other than the current

directory. Likewise, sfmakedepend can be instructed to search other directories with -I dir

options. Note that it is legal to have more than one -I dir option as in:

sfmakedepend -I /usr/local/include -I /home/me/include *.F

Fortran 90 introduces some interesting dependencies. Two compilers I have access to (NAG

f90 and IBM xlf) produce a private my...module.mod file if you define module My_Module in

file mod.f90. This file is used by the compiler when you use the module as a consistency check

(type-safe programming). If foo.f90 uses that module, you will need the following dependency

information:

125
126

foo.o: my_module.mod

my_module.mod: mod.o

This program warns about include files it can't find, but not if there is a "bogus" on the same line.

See the comments at the top of sfmakedepend for up-to-date information on the options. I

may someday get inspired to use a newer version of the getopt routine and rename the options to

have names like -SGI and -Cray.

This says that before compiling foo.f90 we need to have the file myJIlodule.mod. This file in

turn depends on mod.o, so that mod.f90 must be compiled before foo.f90. The sgi is similar

except that it uses the file MY_MODULE.kmo to store the private module information. Use

sfmakedepend -g on the SCI.

Rather than creating extra module files, the Cray and Parasoft compilers store the module

information in the object file and then files which use the modules need to be compiled with extra

flags pointing to the module object files. For instance, if foo.f90 uses My_Module which was

defined in mod.f90, then you will need to compile mod.f90 first and provide the Cray compiler with

the extra option -p mod.o when compiling foo.f90. When using the Cray, use sfmakedepend -c

to get the dependency information:

.~

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

/* bogus include */

#ifndef PLOTS

#include I must_define_PLOTS"

#endif

foo.o: mod.o

$(CFT) $(FFLAGS) -c -p mod.o foo.f90

$(CFT) and $(FFLAGS) are assumed to be previously defined as the name of the compiler and

the compiler options, respectively.

Note: These £90 module dependencies can confuse some versions of make, especially of the

System V variety. We use gnu make because it can follow these chained dependencies and do the

right thing.

sfmakedepend assumes that all the files using and defining modules are in the same directory

and are all in the list of files to besearched. It seems that the industry has not settled on a practical

way to deal with a separate modules directory, anyway.

I sometimes include non-existent files as a compile time consistency check:
,.

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

.~

•

..

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

References

[1] J. S. Allen, P. A. Newberger, and J. Federiuk. Upwelling circulation on the oregon continental

shelf. part i: Response to idealized forcing. J. Phys. Oceanogr., 25:1843-1866, 1995.

[2] A. Arakawa and V. R. Lamb. Methods of computational physics, volume 17, pages 174-265.

Academic Press, 1977.

[3] R. Aris. Vectors, tensors and the basic equations of fluid mechanics. Prentice-Hall, Englewood

Cliffs, NJ, 1962.

[4] G. K. Batchelor. An introduction to fluid dynamics. Cambridge University Press, 1967.

[5] A. Beckmann and D. B. Haidvogel. Numerical simulation of flow around a tall, isolated

seamount. part i: Problem formulation and model accuracy. J. Phys. Oceanogr., 23:1736-

1753, 1993.

[6] A. F. Bennett. Inverse methods in physical oceanography. Cambridge University Press, 1992.

[7] F. P. Bretherton, R. E. Davis, and C. B. Fandry. A technique for objective analysis and design

of oceanographic experiments applied to mode-73. Deep Sea Res., 23:559-582, 1976.

[8] E. F. Carter and A. R. Robinson. Analysis models for the estimation of oceanic fields. J.

Atmos. Ocean. Tech., 4:49-74, 1987.

[9] R. Daley. Atmospheric data analysis, chapter 5. Cambridge University Press, 1991.

[10] E. E. Ebert and J. A. Curry. An intermediate one-dimensional thermodynamic sea ice model

for investigating ice-atmosphere interactions. J. Geophys. Res., 98:10085-10109, 1993.

[11] K. N. Fedorov. Layer thicknesses and effective diffusivities in the diffusive thermocline convection in the ocean. In J. C. J. Nihoul and B. M. Jamart, editors, Small-scale turbulence and

mixing in the ocean, pages 471-479. Elsevier, New York, 1988.

[12] J. A. Francis and A. Schweiger. A new window opens on the arctic. Trans. Amer. Geophys.

Un., 81:77-83, 2000.

[13] N. G. Freeman, A. M. Hale, and M. B. Danard. A modified sigma equations' approach to the

numerical modeling of great lake hydrodynamics. J. Geophys. Res., 77(6):1050-1060, 1972.

[14] B. Galperin, L. H. Kantha, S. Hassid, and A. Rosati. A quasi-equilibrium turbulent energy

model for geophysical flows. J. Atmos. Sci., 45:55-62, 1988.

[15] L. S. Gandin. The objective analysis of meteorological fields. Hydrometeorological Publishing House, Leningrad, 1963. English translation: Israel Program for Scientific Translations,

Jerusalem, 1965.

[16] J. M. N. T. Gray and P. D. Killworth. Stability of the viscous-plastic sea ice rheology. J. Phys.

Oceanogr., 25:971-978, 1995.

[17] J. M. N. T. Gray and P. D. Killworth. Sea ice ridging schemes. J. Phys. Oceanogr., 26:2420-

2428, 1996.

[18] D. B. Haidvogel and A. Beckmann. Numerical models of the coastal ocean. The Sea, 10:457-

482, 1998.

127
[19] D. B. Haidvogel and A. Beckmann. Numerical Ocean Circulation Modeling. Imperial College

Press, 1999.

[20] R L. Haney. On the pressure gradient force over steep topography in sigma coord.inate ocean

models. J. Phys. Oceanogr., 21:610-619, 1991.

[21] W. P. Hazard. Using cpp to aid portability. Computer Language, 8(11):49-54, 1991.

[22] W. D. Hibler, III. A dynamic thermodynamic sea ice model. J. Phys. Oceanogr., 9:815-846,

1979.

[23] W. D. Hibler, III. Documentation for a two-level dynamic thermodynamic sea ice model.

Technical report, USACRREL, Hanover, NH, 1980. Special Report 80-8.

[24] M. Iskandarani, D.E. Haidvogel, and J.P. Boyd. A staggered spectral element model with

applications to the oceanic shallow water equations. Int. J. Num. Meth. Fl., 20:393-414, 1995.

[25] D. RJackett and T. J. McDougall. Stabilization of hydrographic data. J. Atmos. Ocean.

Tech., 12:381-389, 1995.

[26] B. W. Kernighan and D. M. Ritchie. The C Programming Language. Prentice Hall, Englewood

Cliffs, New Jersey 07632, second edition, 1988.

[27] W. G. Large. Modeling and parameterization ocean planetary boundary layers. In E. P. Chassignet and J. Verron, editors, Ocean Modeling and Parameterization, pages 81-120. Kluwer

Academic Publishers, 1998.

[28] W. G. Large, J. C. McWilliams, and S. C. Doney. Oceanic vertical mixing: a review and a

model with a nonlocal boundary layer parameterization. Rev. ,Geophys., 32:363-403, 1994.

[29] J. R Ledwell, A. J. Wilson, and C. S. Low. Evidence for slow mixing across the pycnocline

from an open-ocean tracer-release experiment. Nature, 364:701-703, 1993.

[30] B. P. Leonard.' A stable and accurate convective modelling procedure based on quadratic

upstream interpolation. Comput. Method Appl. Mech. Eng., 19:59-98, 1979.

[31] A. Macks and J. Middleton. Numerical modelling of wind-driven upwelling and downwelling.

University of New South Wales, 1993.

[32] J. Mailhot and R Benoit. A finite-element model of the atmospheric boundary layer suitable

for use with numerical weather prediction models. J. Atmos. Sci., 39:2249-2266, 1982.

[33] J. D. McCalpin. A comparison of second-order and fourth-order pressure gradient algorithms

in a a-coordinate ocean model. Int. J. Num. Meth. Fl., 18:361-383, 1994.

[34] J. C. McWilliams, W. B. Owens, and B. L. Hua. An objective analysis of the polymode

local dynamics experiment. part i: general formalism and statistical model selection. J. Phys.

Oceanogr., 16:483-504, 1986.

[35] G. L. Mellor and 1. Kantha. An ice-ocean coupled model. J. Geophys. Res., 94:10,937-10,954,

1989.

[36] G. L. Mellor and T. Yamada. A hierarchy of turbulence closure models for planetary boundary

layers. J. Atmos. Sci., 31:1791-1806, 1974. '

[37] G. L. Mellor and T. Yamada. Development of a turbulence closure model for geophysical fluid

problems. Rev. Geophys. Space Phys., 20:851-875, 1982.

128

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•
129

[49] R. 1. Schwartz. Learning perl. O'Reilly & Associates, Inc., Sebastopol, CA, 1993. the llama

book.

[38J A. Oram and S. Talbott. Managing Projects with make. O'Reilly & Associates, Inc., Sebastopol,

CA, 1991.

[39J 1. Orlanski. A simple boundary condition for unbounded hyperbolic flows. J. Compo Phys.,

21(3):251-269, July 1976.

[51] P. K. Smolarkiewicz. A fully multidimensional positive definite advection transport algorithm

with small implicit diffusion. J. Compo Phys., 54:325-362, 1984.

[52J P. K. Smolarkiewicz and T. L. Clark. The multidimensional positive definite advection transport algorithm: further development and applications. J. Compo Phys., 67:396-438, 1986.

[53] P. K. Smolarkiewicz arid W. W. Grabowski. The multidimensional positive definite advection

transport algorithm: non-oscillatory option. J. Compo Phys., 86:355-375, 1990.

[54] Y. Song and D. B. Haidvogel. A semi-implicit ocean circulation model using a generalized

topography-following coordinate system. J. Compo Phys., 115(1):228-244, 1994.

[55] M. Steele, G. L. Mellor, and M. G. McPhee. Role of the molecular sublayer in the melting or

freezing of sea ice. J. Phys. Oceanogr., 19:139-147, 1989.

[56] R. Styles and S. M. Glenn. Observation and modeling of sediment transport events in the

middle atlantic bight. In 8th International conference on Physics of estuaries and coastal seas,

1996. submitted.

[50] P. K. Smolarkiewicz. A simple positive definite advection scheme with small implicit diffusion.

Mon. Wea. Rev., 111:479-486, 1983.

[47] R. Rew, G. Davis, S. Emmerson, and H. Davies. NetCDF User's Guide. Unidata, University

Corporation for Atmospheric Research, Boulder, Colorado, 1996. Version 2.4.

[48J B-. D. Richtmeyer and K. W. Morton. Difference Methods for Initial- Value Problems. Interscience Publishers, J. Wiley and Sons, New York, New York, second edition, 1967.

[44] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes, The

Art of Scientific Computing. Cambridge University Press, 1986.

[45] P. J. Rasch. Conservative shape-preserving two-dimensional transport on a spherical reduced

grid. Mon. Wea. Rev., 122:1337-1350, 1994.

[46] W. H. Raymond and H. L. Kuo. A radiation boundary condition for multi-dimensional flows.

Quart. J. R. Met. Soc., 110:535-551, 1984.

[43] N. A. Phillips. A coordinate system having some special advantages for numerical forecasting.

J. Meteorology, 14(2):184-185, 1957.

[42] H. Peters, M. C. Gregg, and J. M. Toole. On the parameterization of equatorial turbulence.

J. Geophys. Res., 93:1199-1218, 1988.

[40J R. C. Pacanowski and G. H. Philander. Parameterization of vertical mixing in numerical

models of tropical oceans. J. Phys. Oceanogr., 11:1443-1451, 1981.

[41] C. L. Parkinson and W. M. Washington. A large~scale numerical model of sea ice. J. Geophys.

Res., 84:6565-6575, 1979.

~,~--------------~----------..-.III

•

•

•

•

•

..

•

•

'~•

•

•

•

•

•

•

•

•

•

•

•

•

•

..

•

•

'.•

•

•

•

r.•

'.•

•

•

•

•

•

•
[57] P. Tag. A comparison of several forms of eddy viscosity parameterization in a two-dimensional

cloud model. J. Appl. Meteor., 18:1429-1441, 1979.

[58] J. Thuburn. Multidimensional flux-limited advection schemes. J. Compo Phys., 123:74-83,

1996.

[59] 1. B. Troen and L. Mahrt. A simple model of the atmospheric boundary layer; sensitivity to

surface evaporation. Boundary-Layer Meteor., 37:129-148, 1986.

[60] L. Wall, T. Christiansen, and R. L. Schwartz. Programming perl. O'Reilly & Associates, Inc.,

Sebastopol, CA, second edition, 1996. the camel book.

[61] J. Wilkin and K. Hedstrom. User's manual for an orthogonal curvilinear grid-generation

package. Institute for Naval Oceanography, 1991.

130

.\111

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•

•.>

•

.,

•

•

•
'.

'

,.- .... ~.

".

," ~"\

......;

t .",

to ,"J~~ f"

'. 'f"

o',•

... ,".

. "':' ...... ' ... '1 " .;".\ /tt.: ~, ~.. .' .(: ...

". ,;.., '.'

I, '-> "> .• , " ,

t; I,

",

"'J •• ,!, ~..,

vd"

;. "w''''~

,

'I

.,~,.,. •• I. ',', , . l. t~)I •

.' ',.. ~,'

,t . ,

(,.' "

"j' ,·...i

"

".'"

• 'f'

." "'S' ,," ,

-.. .'

·r•

t:· , ~,

,~\.'" I'

',: ."

,,>,'. " ......

" ....,.

. \- i"t; .( ' . :;1 '....

'0'

" ~'I

,"

',;.'

';".

'''t

....

.

,'" "I' , '.

.' ~'.

" ••,"t

:.;.:' ';C t~" ,I

'",1

,, ~,~., ,\ ""

~', .'

""/.

,;;.

"

~\-,

.•..\:.; ~~,~'.

1\

.,

..

,t,'

',~ .',. ....

,..'.,;

Or. ;._ )';- '

'<',.-,-

~~ .;,'"

~;..

'1-•.••.,

. .. ·.t , •• !

/""

......'

. :f' ~"

',J.• ,'

',i,.

\.'<. <."

- <'

...1,"... ··(,

..... ,; ;/~" '\'..,.

> "'-. I, t" \. ., . ) .;~. .,,~ '. ':',

~;r';

,"

,'

~f

.

,~,.

. ;. '. I;' ".,

)I." '''; ,\' ••-"

":
~I

'"

!. .'

.'. ~, ...

-' .,.

..'. '-".'

,'.r ,,' J

"

.. '" ,',

", ,

't 1 J. . -.... ~~ .,i. ",I..

" ,,','

",

.\ '.

..,.' ,

t··..

' .., , '

" \,.,

.•...:.r."

.,

.

,i

1 '

. ..."!"