Home

Technical Report No. 19 PIPE ISSN 0940-9327

1. 103 5 1 3 Jnstallati n Step NOx ae ae E A O ete Mee nde Sy ne 104 5 1 3 1 Installing the Plot System a s u5 A as 104 5 1 3 2 Installing the Postprocessing 0 0 104 da SOI dip PIPE cao be gee SA a oe RS Se ke eee es 105 524e Th Code Preprocessor sica ry a ay Ba Sea Be MA AS 105 522 o th Gidig io A Goad ee alae ee Bla pe ety ee ee Ge a 106 5 2 3 Specification of the Dimensions 00 8 106 5 2 4 Defining a Focus e ae Gece A ee gow kk ee es 106 5 2 5 Tuning Topography and Coastline 2 2 020 000 2020 107 5 2 6 Selecting the Forcing Data 24 003 toa a a ee ee 108 5 2 7 Defining the Output rc A A ARA 108 5 3 Examples for Specific Model Layouts o o e 109 53 1 Predefined Grids pr a baa aaa A A 109 5 3 2 Self defined Grid 020 0002000 0000000000000 00804 110 5 4 9 0 5 3 3 OR Modela e ua ES aa a A Be E eee ee 110 Rubin PIPE ie sio A ls A A 110 5 4 1 How to Generate an Executable Model Version 110 54 2 Start from Scratch ss 4 40e 2 A E A eth eee AA ee ek 110 543 Data PEEPTOCESSIN E lt o i bal A A Bee ee Ee Eck a 111 Hoel Data Bank Lit bo kes IES AA eee eS 111 5 4 3 1 1 Ocean Initialization 2 0 2 2 a 111 BA sul GOADS Data e awe E AA A 113 5 4 3 1 3 ECMWF Data ocio Sed a 115 5 4 3 1 4 Precipitation Data i 6 40 4 ek aaa 117 5 4 3 1 5 Land Sea Masks ecb oe ak ak eB e 117 5 4 3 1 6 Open Boundary Forcing o
2. phS ph6 ALF Che 2 31 P phS y AtE s 2 32 2 where puh represents the mass flux ph the mass content ph0 the heat content and phS the salt content All equations are discretized in time using an Euler scheme G represents all those terms that are treated implicitly and F those terms that are explicitly computed with forward steps In the predictor step the pressure gradient the Coriolis term and the vorticity part of the advection in the momentum equation the flux divergence in the continuity equation and the advection and diffusion terms in the equation for potential temperature and salinity are treated implicitly All remaining terms are collected in F In order to demonstrate how the technique works it is explained here in detail how to derive the wave equation A simplified derivation of the semi implicit scheme can be found in Appendix A In the first step the momentum equation is discretized in time by using a centered Euler scheme However the code allows for an arbitrary choice for the forward and backward weights The x and y components of the momentum Eqn 1 1 then read de Al a At de puh t puh PAP LAR O lo Foun 2 33 Sn A a At ia a pvh 1 pvh h TS f puh a Ph 2 34 2 Oy 2 where Fr and Fen represent all the remaining terms taken at time level n In the next step the equations are solved for puh or pvh by eliminating the
3. e lt TRIONE gt computes the solution of NV independent pairs of coupled linear equations of block tridiagonal type with NX number of equations It is used be lt SEAICE gt to find the solutions for the rheology terms in the momentum equation Similar to lt TRITWO gt the pair of equations couple the x and y components of the ice flux in order to obtain also a fast convergence for the bulk and shear viscosity terms e lt TRITWOX gt computes the solution of NV NZ independent pairs of coupled lin ear equations of block tridiagonal type with NX number of equations It is used by lt MOLVERX gt and is required to couple x and y components of the advection and cur vature terms since flux components are cross referenced in the momentum equations e lt TRITWOY gt computes the solution of NV NZ independent pairs of coupled lin ear equations of block tridiagonal type with NY number of equations It is used by lt MOLVERY gt and is required to couple x and y components of the advection and cur vature terms since flux components are cross referenced in the momentum equations e lt TRIBLCKX gt computes the solution of NV independent sets of linear equations of block tridiagonal type with NX NZ number of equations Each linear equation system solves the matrix on a xz section NV NY 1 2 is the result of the used line relaxation method e lt TRIBLCKY gt computes the solution of NV independent sets of linear equations of block tridiagon
4. lt jpi_hold gt performs a barrier Part III User s Manual Chapter 5 How to Use PIPE 5 1 Installing PIPE The PIPE model comes along with a tar file containing codes scripts and data Assuming that the file is untarred in a directory PIPE then the directory structure becomes as follows PIPE Source Source Manual Source PVP PVP Jobs PVP Model Source RISC RISC Jobs RISC Model Source MPP MPP Jobs MPP Model Model bin Model lib Model code Source DOS PIPE PostPro PostPro Codes PostPro Scripts PostPro Help PostPro Objects PIPE Graphic Graphic Manual Graphic GrADS PIPE Tools PIPE Data PIPE bin full source codes PostScript file for Techn Rep No 7 preprocessed source code for global T42 scripts for compiling and running the code example directory for specific model version preprocessed source code for global T42 scripts for compiling and running the code example directory for specific model version unpreprocessed source code scripts for compiling and running the code example directory for specific model version directory containing makefiles and executables directory containing libraries directory containing fsplit d code preprocessed code for DOS home directory of postprocessing software codes for modules scripts for modules documentation files for modules object files for modules graphic library manual for graphic library GrADS alternat
5. 5 2 7 Defining the Output The model delivers the following files 1 File OCDAT contains all data required to restart the model from an earlier state It uses binary format although an ASCII file can alternatively be created see switch ISW31 During the start the model is able to read the binary or the formatted version of OCDAT The data format is recognized automatically 2 File TRDAT holds the tracer distribution for restarting the tracer model like OCDAT holds the distribution of the dynamical quantities TRDAT is generated at the end of the initialization and is used for restarting the tracer model 3 File FORCES contains all monthly mean atmospheric forcing data as well as the monthly mean SST and SSS on the model grid It is automatically created from the chosen global data sets of atmospheric forcing sea surface temperature and sea surface salinity when OCDAT is found to be empty However it must be ensured that all re quired global data are available see model sun or model cri Depending on SW 42 FORCES also can be created if an OCDAT exists or it need not be created despite OCDAT being empty This is useful if FORCES has already been created but if one wants to start the ocean state from initial conditions 4 File OBNDSEC holds the boundary data i e the temperatures salinities and barotropic transports along the xz and yz sections on the model boundaries If the tide model is used also
6. If the tide equations for simplicity written in non dimensional form are discretized in time by t ut FP a ar a id 2 47 E At a t F gt 3 y or apa 2 48 ie At yn QO pS Fe a a as Hao 2 49 where Fa F and Fp are defined by At Fa u p fu 2 50 ut Sp fo 2 50 At 0 w 2 51 e EEr fu 251 ie o o Fy p u v 2 52 and if F 18 defined as At o ee Fa Fo 2 53 2 a Se el Oy then the equations for the pressure and the velocities at the new time level become AP o a At n l A m EL n 1 F A pnl 2 54 p Par aes ape NS a 2 54 3 pynt a n en Fi SF 2p mn i 2 55 ees At At n41 _ At f pnl jar y PEREA A 2 56 e At i a The latter three equations form a system of equations for the three unknowns at the new time level 1 e the pressure and the two velocity components Because the factor aif lt 1 for tidal applications the vorticity term in Eqn 2 54 can be updated during the iteration for the pressure by using Eqns 2 55 and 2 56 for the velocity components at the new time step Convergence is guaranteed if the time step is shorter than the inertial period The time step of the tide model DTIDE is hardcoded in routine lt POTTIDE gt and currently is set to 900 seconds 2 5 Open Boundary Formulation The basic philosophy to install open boundary conditions into PIPE is to feed the model
7. e lt ECCLD gt reads one month of the ECMWF cloudiness from CLD and interpolates it onto the model grid e lt ECSLP gt reads one month of the ECMWF surface pressure from SLP and interpolates it onto the model grid e lt ECHUM gt reads one month of the ECMWF surface air humidity from HUM and interpolates it onto the model grid e lt LEGATES gt reads one month of the Legates amp Willmott precipitation data set from RAIN1DEG and interpolates it onto the model grid e lt SEATEMP gt reads the Levitus monthly mean surface temperature from TEMP MON from the local file SSTEMP and interpolates it onto the model grid lt SEASALT gt reads the Levitus annual mean sea surface salinity from SSSALT or from the monthly mean sea surface salinity from SALT MON from the loacl file SSSALT and interpolates it onto the model grid lt WIND gt reads one month of the Hellermann and Rosenstein surface wind stress from UVGLOB and interpolates it onto the model grid lt UPDATE gt computes the new current forcing of all the interpolated forcing data from the current forcing and the next stored monthly mean lt WRITEF gt saves the interpolated atmospheric forcing on FORCES Entry lt READF gt reads the forcing from FORCES 4 4 4 Initialization of Ocean State lt LAYOUTI1 gt initializes the x and y coordinates reads the topography from TOP1DEG interpolates it onto the model grid and arranges the land sea mask
8. 102 be u T AO 1 7 x 10 T Aq 1 62 where Ag is the difference between the specific humidity of air and the sea surface and AO is the potential temperature difference The only change from Large and Pond s work is that cmv is not fitted against data using ad hoc chosen curves but by tuning the Charnock constant The equations 1 51 1 54 and 1 55 describe the dependence of the neutral drag coefficient cyy on the friction velocity w the Karman constant the height of the measurement Z and Charnock s constant Char In order to obtain a drag coefficient of about 1 15 x 1073 for neutral conditions at 10 MfS Cehar should be set to 0 0064 However as outlined by Oberhuber 1988 Cenar in fact is set to 0 032 to compensate for the underestimation of the transfer coefficient resulting from the application of monthly mean values instead of instantaneous values Note that the standard forcing of PIPE is based on monthly mean winds 1 3 3 Evaluation of the Net Fresh Water Flux A Newtonian relaxation is used to force the salinity equations with fresh water fluxes see also JMO The simplest way is to use relaxation towards the observed salinity Optionally data for precipitation evaporation which are computed from latent heat and river runoff may be used too The relation is obs 5 Qr Fi pie AR ith E 1 63 U A E a Pucon ae where is the time constant with which the actual salinity relaxes towards the obse
9. lt DECAY gt lt DECAY gt lt BUBBLES gt time scale for detrainment tuning parameter penetration depth of solar radiation transmission coefficient parameter for salt ejection lt lt lt MIXEXP gt lt lt Table 4 10 Parameters for Mixed Layer Variable Value Default Unit Reference XGL XGR YGU YGO IDENS JDENS ICENTER JCENTER IBAND JBAND XCOM YCOM IORIGIN JORIGIN any any any any any any any any any any any any any any LEAVOUT any gt 0 Variable ALPHA BETA GAMMA ISWITCH 0 1 Meaning left margin of model domain right margin of model domain lower margin of model domain upper margin of model domain number of extra grid points in x dir number of extra grid points in y dir x index location of focus centre y index location of focus centre index width used for focus in x dir index width used for focus in y dir refinement factor in x direction refinement factor in y direction x index origin of focus y index origin of focus number of ignored latitudes Table 4 11 Parameters for Horizontal Grid Definition Value any any any Unit deg deg deg Meaning first longitudinal rotation angle second latitudinal rotation angle third longitudinal rotation angle switch for activating rotation Table 4 12 Threshold Variables lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY
10. o 118 543 2 Model Forcing caras dador PR APES BA ES 119 5 4 4 The Multi Grid Method o e e e 119 D45 Diagnostic Output ses es anra naa aa a AAA 119 5 4 6 How to Apply the Coupling Interface o o aaa aaa 120 5 4 7 How to Use the Bias Correction aooaa aa a a a 121 5 4 8 About Numerical Stability of PIPE aaa 122 D amp I Numerical Filters se a A oa A Be a 123 D410 Error C nditions s s 2 45 ach Some 49 de a Ae OO BGs AS a 124 5 4 11 Conceptual Problems oi odas an 24463 244 63 124 5 4 12 About Code Consistencia a a a eee A AS 125 POS processo e a A A AAA 125 Dol Code WennitiOns s tae Mesa ita ia Be AAA AA ya 125 5 5 2 PPSF Binary Formats adan ale ee sa ee RS 125 5 5 3 Portable Compressed Format o 0208 129 DOA Tist of PPSESEMES errada aa A aaa 129 5 5 4 1 The Quick Look System o a 129 5 5 4 1 1_ File PPSF_SUR lt 2er ias a 129 Sota File PPSFICE erone ea shee he kee 4 129 5541 3 File PESE BAS adds oe PRA eS BS 130 so dido File PS E SEO at A AAA a 130 5 5 4 1 5 File PPSF_VER xicos isa ee Sat es 130 SoLo File PPSPELX 726 2024 A 130 5 5 4 1 7 File PPSF_FOR 5 2 24 44 aga Oe ace a amp 130 5 5418 File PPSF TID sera 2e 2 Dae eis ae was 130 55419 File PPSE TOP 2 ain e a e bea ea ho 130 5 04 10 File PESEINT os 4 4 4 eee ace a Bit Bae 4 131 pov Fhe History Miles e a eal eo a Be ace A Ae ee Re 131 5042
11. Az and Ay values are needed both on vector and scalar points Since X and Y are defined on vector points Az on a scalar point is obtained by Ars coste yin Xy Xy 1 1 2 3 while Az on a vector point is Xy r41 Xv 1 Azry r J c0s si 5 2 4 Ay is obtained on scalar points by Aysu Yvan Yvu 1 2 5 and for Ay on vector points AYv 1 J ES en 2 6 2 For better vectorization values of Az are precalculated and stored into the arrays HDXZX and HDXX for differences at scalar points and into VDXZX and VDXX for differences at vector points 2 1 1 1 Boundary Conditions Basically all horizontal boundary conditions are derived from the array IFLG which is defined on vector points This array contains the information about the coastline geometry However it also carries the locations where a physically existing layer with non zero thickness continues as a zero layer At these locations a boundary condition has to be computed see JMO Mathematically these time variable boundary conditions are introduced in the same way as the time constant coastline geometry The only generalization is that JFLG contains a 3rd index for the layers While JFLG is constant for the 1st layer IFLG is computed for all deeper layers through the algorithm that sets correct boundary conditions near vanishing or reappearing layers IFLG contains 0 for land and 1 for ocean Boundary conditions are set at vector points by multiplying each
12. Ur Up h hg 1 h h Pu 9 94 where g _ and g are the stabilities across the upper and lower interfaces and h is a tuning coefficient Compared with JMO the weight W has been introduced to make the vertical diffusivity less dependent on the stability The expression for W enhances vertical mixing in the deep ocean while mixing is unchanged in the upper ocean The available turbulent energy 2mou3 used for internal mixing is defined as 30 4 1 q a a 1 71 218 Uy Uy 1 where UY is a constant background energy u is the turbulent energy from Eqn 1 64 and qi the sea ice concentration This parameterization for TKE available for internal mixing is made dependent on the ice concentration in order to express that TKE cannot be produced if ice decouples the ocean from the energy source which is the wind The factor 3 in Eqn 1 71 parameterizes the effect that an internal energy source also might exist The 2nd term in Eqn 1 71 parameterizes the effect that wind available from data and used by the mixed layer model also causes TKE below the mixed layer through braking internal waves This energy available below the mixed layer partly reduces the effective entrainment rate into the mixed layer This parameterization was added in order to avoid unrealistic deep mixed layers in the subtropics i e under conditions of strong wind and strong stratification The problem of vertical mixing is underdeter
13. e lt ROLLX gt shifts 1 dimensional real data in x direction in order to improve convergence e lt ROLLZX gt shifts 2 dimensional real data on a xz section in x direction in order to improve convergence 4 5 11 Routines to Bookkeep Budgets e lt QBDGT_INIT gt initializes arrays that boodkeep the budget fields Entry lt QBDGT_ADD gt adds current budget terms onto already existing integrals Entry lt QBDGT_PPSF gt writes out budget terms at a specified frequency 4 6 Passive Tracer Model ocetrac f e lt TRAINIT gt reads initialization data from TRDAT If it is empty the tracer fields are initialized e lt TRAFORC gt determines the kind of forcing Please adjust to your needs e lt TRASTEP gt organizes initialization time stepping and postprocessing e lt TRAPOST gt performs the postprocessing e lt TRASTOP gt cleans up and saves the data in TRDAT e lt TRACPAS gt performs the tracer forecast based on a scheme similar to lt TRACTIV gt for the active tracers e lt TRAVECT gt performs the vertical convection when stratification becomes instable e lt TREAD gt prepares the ocean data like flow layer thickness density vertical exchange rates and others to provide data for the time stepping Dependent on whether the tracer model is running in online or offline mode its forcing data either are taken from the ocean model s common blocks or are read from TRAFORCE e lt CSOURCE gt is an empty routine
14. filter q KODEQ s KODES 1 KODEL Input_File Output_file Purpose Extract data set from Input_file and write onto Output_File Definition s Integer KODEQ KODES KODEL File Input_File Output_File Optional Parameter s q KODEQ where KODEQ is the Quantity Code s KODES where KODES is the Section Code 1 KODEL where KODEL is the Location Code Remark s If one of the KODEs matches codes in the data set these are written to Output_File The following commands consist of scripts codes and documentation files The scripts are tested on a SUN station and on a CRAY using CShell Dependent on the operating system the script decides whether the assign or the link command is used to link the data file with the local file in order to avoid extensive copying of data across directories or disks They also use different loader commands 5 5 6 1 Environment Variables In order to find the subsequent scripts for execution the search path for in the example HOME OPYC PostPro Scripts needs to be defined in profile cshrc or login under UNIX Furthermore all scripts expect an environment variable OPYC_PPS_OBJ to be defined as e g SHOME OPYC PostPro Objects in order to find the object files 5 5 6 2 Utilities for Manipulating Data Sequence e Script append merges a sequence of input files e Script delete deletes specified codes on the input file e Script extract extracts blocks of data from an input file by de
15. perature and salinity for an xz section lt INTEGRAY lt performs the vertical integration of potential density potential tem perature and salinity for an yz section lt ROTVEC lt rotates a vector from the geographical coordinates onto a rotated grid 4 7 3 Routines for Pressure Calculation lt STATETRX gt computes the in situ density for an array with dimension NX Entry lt STATETPX gt computes the potential density for an array with dimension NX lt STATETRY gt computes the in situ density for an array with dimension NY Entry lt STATETPY gt computes the potential density for an array with dimension NY lt STATEVRX gt computes the in situ density for an array with dimension NZ NX Entry lt STATEVPX gt computes the potential density for an array with dimension NZ NX lt STATEVRY gt computes the in situ density for an array with dimension NZ N Y Entry lt STATEVPY gt computes the potential density for an array with dimension NZ NY lt TTOTETX gt computes the potential temperature for an array with dimension NX lt TTOTETY gt computes the potential temperature for an array with dimension NY lt REPRESSX gt computes the pressure for an array with dimension NZ NX lt REPRESSY gt computes the pressure for an array with dimension NZ NY lt NEWDENSX gt computes iterativly in situ density from potential temperature and salinity for an array with dimension NZ NX lt NEWDENSY gt computes iterat
16. 0 0000 0 139 5 6 1 11 Executing Program OCEPLOT 139 5 6 2 The Underlying Library JMOPLANE 140 50 3 Example lots poa a ble pw ao wa BAL ee a E A E A 140 Doe Ser DOGS EEE A a e Aa Ae 140 6 Appendices 145 6 1 Appendix A Prognostic Pressure Equation o o 145 6 2 Appendix B Pressure Boundary Condition 0 o 145 6 3 Appendix C Constants for Equation of State 146 IV References 147 7 How to Find References 149 fal Code emer OnGes a St e ip ia here Se Re GO oe e ey OR 149 hake Acknowledgements gt ccd ance A AS E AAA 154 hese Iteratute Referentes m isane A tas Peat pena ta a a a R me ae 154 TSE Text RETGTER CES ft tk Se A em are ck ga ge Ak he lee ees 154 7 3 2 Key Publications with PIPE sorted in time 156 List of Tables 1 1 Constants of Major Tidal Modes 0 0 0 0 38 4 1 Time Step Control Parameters 0 0200000004 64 4 2 Coupling Interface Parameters o 65 4 3 Parameters for Annual Mean Flux Adjustmedt 65 4 4 Heat Flux Parameters e 66 4 5 Parameters for Horizontal Diffusion o o 66 4 6 Snow and Sea Ice Parameters o 67 4 7 Parameters for Coordinate Generation 67 4 8 Some Physical Parameters 0 0 0 20000002 eee 68 4 9 Parameters for Diapycnal
17. In some cases 2 d arrays are extracted out of a 3 d array and scattered back if required e Equations that must be solved either iteratively or directly are treated such that direct solutions are obtained by an alternating direction line relaxation ADLR scheme It con sists of two steps first a direct solution on the zx plane and an iteration in y direction and a second direct solution on the yz plane and an iteration in x direction 2 6 1 PVP Architecture If the code is setup for a PVP architecture iterations are performed such that a parallelization is performed over each second latitude or longitude respectively for the ADLR scheme For instance a coarse resolution T42 global ocean version of the model achieves 450 Mflops on a CRAY C90 and requires only 6MWords of memory Parallelization works best when granularity is coarse i e a lot of work is carried out in parallel regions in order to avoid latency time of frequent initializations of them Therefore the code contains micro tasking directives which tell the compiler that a loop should be executed in parallel Routines which cause more overhead than speedup due to parallelization as well as those which are called inside a parallel region are excluded from parallelization Due to fine granularity of the code the speedup is more moderate than excellent The code achieves for instance a speedup of 3 2 on 4 processors or 5 4 on 8 processors for a coarse resolution T42 version 2 6
18. TE Toqi TA qi 1 91 This yields the snow surface temperature T or the ice surface temperature Tp by simplifying the total heat flux Q from 1 33 through a Taylor expansion around Ta dQ Q T Qa Gp T Ta 1 92 This equation is used to formulate the two equations for T and Th as two linear equations that can be solved easily Since the skin temperatures influence the transfer coefficients via the stability dependent Large and Pond coefficients the final solution for the skin temperature is obtained by iteration The resulting T and T are used to compute the heat flux Q that is converted into melting snow and the heat flux Q that is converted into melting or freezing sea ice ks Ps Q QTn SEEN Tn 1 93 ki 4 Qi Q T a Th T1 1 94 Thus the heat flux Q used for melting snow is the difference between the atmospheric heat flux and the conductive heat flux through the snow while the heat flux Q used for melting or freezing sea ice is the total heat flux Q T into the ice free ocean plus the conducted heat flux through the ice floe Note that if T lt Tm then Q 0 This is guaranteed by 1 88 however this equation is not used when 1 89 and 1 90 yield T gt Tm In this case the error in 1 88 is interpreted as heat flux used for melting snow The linearity of 1 87 and 1 88 ensures that this flux is always downward The resulting change in the local ice thickness Fr due to thermodyn
19. switch for black amp white or color plots switch for rotation of picture on sheet number of grid points for x shifting NQUAL MUE NUE IPART JPART NFILTER NSIGN NFILM NOFFSET NTRAJ eT Toa 1 switch for quality of contour plots any gt 1 any gt 1 any gt 1 any gt 1 0 1 NPRETTY 0 1 0 1 0 1 any any x refinement factor for contours y refinement factor for contours subdividing figure in x direction subdividing figure in y direction switch for spatial filter switch from model mask to real mask switch from solid to dashed contours switch from sheet to movy mode time offset to model time counter switch for trajectories instead of vectorrs Oooo CO FR e e me re Table 5 10 Definition of Main Switches JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 5 6 1 3 Block Data STRINGS lt STRINGS gt contains the relation of the first code number to the physical quantity Data are transferred into lt OCEPLOT gt via JMOSTRI and are required to define the header for each plot This program unit contains a complete list of code definitions Variable Value Unit Default Meaning Reference XLEFT any deg 0 location of left margin JMOPLT2 XRIGHT any gt XLEFT deg location of right margin JMOPLT2 YLOW any deg 90 location of lower margin JMOPLT2 YUPP any gt YLOW deg 90 location of upper margin JMOPLT2 NXS any gt 2 13 num
20. vh of the ice thickness h of the ice concentration q and of the snow depth s are with the latter taken as equivalent ice thickness 0 gt a rg T P Ta To 3 ay 0h V AV wh Fx wh ghi VOZ he D n A diri k 1 a 1 gt gt gt ail V A Vh V vh Fn Fa 1 78 Zi gt ae AN Ge Va vq Fy 1 79 gt gt gt a V A Vs V us F Fa 1 80 where 7 and 7 are the surface wind stress and the stress at the bottom of the ice and Shs hy D is the sea surface elevation In these equations the transport of momentum has been neglected A is a constant diffusion coefficient fis the Coriolis parameter E Fr Fa Fs and F are the forcing functions for the ice momentum ice mass ice concentration snow mass and the aging process respectively In more detail F represents the internal ice stress determined by a viscous plastic ice rheology F the ice thickness change due to freezing or melting F the change of the ice concentration due to external heat fluxes F the change of the snow mass due to snowfall or melting and F the conversion rate from snow to ice that describes the aging process of snow The sea ice rheology reads Ou Ov Ps Ou Ov Po ee ee a ae 1 81 lA Sn Sy 181 Ov Ou Py Ou Ov Fy o ae A 1 82 A 1 82 Iui Avian 1 e Ou via 1 u Ov e 1 2 1 A ay a e E
21. 1 Block Data SWITCH defines switches that control how the forcing is com puted and whether observed or artificial data should be used see 4 18 The switch IS W48 is noteworthy With SW48 1 the model computes the running mean fluxes of heat fresh water Variable Value Unit Meaning Reference CP 4180 Wskg K specific heat capacity of water JMOPAR GRAV 9 81 ms earth accelleration JMOPAR DRAGT 5 E 6 s7 time constant for salinity forcing lt SURFLUX gt DRAGS 1 5E 6 drag coefficient at surface lt STRESS gt DRAGB 1 E 4 drag coefficient at bottom lt STRESS gt Table 4 8 Some Physical Parameters Variable Value Unit Meaning Reference EPSILON 0 10 m s vertical mixing damping coef lt CROSMIX gt RICRIT 25 critical Richardson Number lt CROSMIX gt HVMIX 500 m tuning coefficient lt CROSMIX gt TURBLEV 4 E 7 m3s 3 tuning coefficient lt CROSMIX gt Table 4 9 Parameters for Diapycnal Diffusion and stress if NINIT 1 If the length of the run is a multiple of a year the arrays HEATQ will contain the mean heat flux PMEQ the mean fresh water flux and TAUXQ and TAUYQ the mean stresses If the user switches to SW 48 2 the model will continue to compute the running mean however it will also use this mean value for forcing the surface temperature and salinity After every finished year these running means are stored onto HEATM PMEM TAUXM and TAUYM If the run with SW 8 2 is long enough and again a multiple of a ye
22. 1 File PP SE ALD asias wis aa e 131 90422 le PRBSE UIS amnesia ura boa ea hots 131 dao geti and Put se ama Be ek ee at Ae See BS QA Se BOS 131 5 5 6 Utilities for Data Analysis 2 4 q 4 ae gale De ace e eee 4 131 5 5 6 1 Environment Variables saaa p aaa e ra k ek 2 020048 132 5 5 6 2 Utilities for Manipulating Data Sequence aaou aa 132 5 5 6 3 Utilities for Manipulating Data Contents 133 5 5 6 4 Time Series Analysis lr a ee aon ee See 134 5 5 6 5 Frequency Analysis o o 4 134 00 00 Data Piltene a wae Gee a a ek 135 5 5 6 7 Data Reformatting 44 4 te a e a A 135 A Viewing A RE Ad Be aod Bia 135 5 5 6 9 Data Preprocessor cine de qa Ga eed Be ee we Bes 136 5 5 6 10 Data Postprocessor 2 64402206 68 a4 b a es 136 DOr Ds Flot SystemC Scion Gag tcp teeta Aa A Mae te Re oy Nee ee Bene ae at emi ng 136 als Si be Plog Prosa ride te E A Re hn ri 136 5 6 1 1 Main Program PICTURE ses Gane ak a aioe es ta 137 OO Block Data PLEOTPAR ax 4 a2 8 aea Pace ee tds ida 137 5 6 1 3 Block Data STRINGS exec ye a a aa 137 5 6 1 4 Block Data MARKS 12 2 dio bebe ee Re ee 138 5 6 9 Block Data UNITS 40 das a ee A HS 4 138 5 6 1 6 Block Data ALLIND e oa te aoe aa e a 138 5 6 1 7 Block Data MINMAXS 4002 eb 44S a bw es 139 5 6 1 8 Block Data GVA RIA aooaa aaa o a ee a ae 139 5 6 1 9 Block Data ADDFACS x ston ok eee ee Ok eS 139 5 6 1 10 Block Data QZONMER 0
23. 12 DO MONTH 1 12 DO J 1 160 READ 2014 C ISTD I I 1 360 DO I 1 360 STD I J MONTH ISTD 1 1000 ENDDO ENDDO ENDDO 5 4 3 1 4 Precipitation Data e File RAINIDEG contains the precipitation data set of Legates amp Willmott 1990 However it has been averaged onto a 1 x 1 grid from the original 0 5 x 0 5 grid It is read by lt LEGATES gt from the local file MONRAIN with the following simplified statements DIMENSION IRAIN 360 RAIN 360 180 12 DO MONTH 1 12 DO J 1 180 READ 2014 IRAIN 1 I 1 360 DO I 1 360 RAIN I J MONTH IRAIN I 1000 2592000 ENDDO ENDDO ENDDO 5 4 3 1 5 Land Sea Masks In order to make the setup of various global models simpler a number of readily tuned masks are offered This is done because certain combinations of land and sea points lead to B grid splitting The following masks don t have this problem Thus a tuning of the land sea mask is not needed when these masks are taken These files are read by lt FLAGMOD gt from the local file LSMASK However one has to ensure that the script setup is consistent with the script model sun or model cri which has to copy the corresponding grid file onto the local file LSMASK The available pre defined grids are File T21GRID has dimension 64 32 File T32GRID has dimension 96 48 File T63GRID has dimension 192 96 l File T42GRID has dimension 128 64 l File T84GRID has
24. 2 RISC Architecture The RISC benefits from the optimization for the MPP architecture because PIPE has been developed for a CRAY T3D which is based on RISC chips from DEC There exist many code sections that are alternatively formulated for cache based processors The key word used by the cpp or fpp preprocessor is RISC Basically RISC processors have a memory bandwidth problem i e memory is slow Because PIPE is a memory intensive code the slow memory does not allow a high speedup when the code is compiled as parallel code on a workstation Therefore it does not make sense to use more than 4 processors 2 6 3 MPP Architecture The model is parallalized for MPP computers Oberhuber and Ketelsen 1995 In order to execute the model in parallel on many processors like possible on a CRAY T3D the code is split into four executables see title page of this report These executables have the task to 1 initialize the data structure on the server to be able to efficiently supply the client with data which are read from the files OCDAT FORCES and OBNDSEC 2 perform the model integration by running two executables e one running on the server is supplying the client with data if required e the other one is running on the client performs the time integration and returns data to the server for e g later analysis 3 clean up the intermediate results and return a data structure identical to a model run on either a PVP or RISC machine i
25. 6 and the relation for the in situ pressure 1 7 and 1 8 Eqn 2 40 determines h The density at the new time level is updated during the iteration by using the last guess for the layer thickness and the previously determined potential temperature and salinity Finally after having found the solution for layer thickness and density the mass fluxes are obtained from Eqns 2 35 through 2 36 The problem is formulated as a linear equation with frozen matrix coefficients during one iteration step Between iteration steps the matrix coefficients are updated with a method that avoids oscillations around the solution The layer equation is formulated as a system of linear equations in the x and z direction and iterated in y which is an application of the line relazation method Optionally the equations are formulated directly in y and z direction in addition in order to improve convergence The layer equation consists of quadratic terms which represent the inertia gravity waves and mixed derivatives which correspond to Rossby waves Because of the large time steps used in this model At f 2 gt 1 the mixed derivative terms dominate the quadratic terms so that simple iteration methods used for Laplacian type equations cannot be used Based on the fact that an iteration converges if the mean diagonal elements dominate over the neighbouring diagonals a method has been developed which calculates an optimal relaxation coefficient for ever
26. An e Ox y e 1 58 P PM exp ei 1 a 1 84 t Fo TANA 1 85 j 1 86 P is a proportionality coefficient e is the excentricity of the viscous plastic rheology ey is a threshold for the computed bulk viscosity e is a decay parameter and p is a constant ice density Parts of the operators on the right hand side can be identified as a diffusion operator with two diffusion coefficients the bulk viscosity and the shear viscosity 7 which are highly dependent on the flow field However additional terms occur that cross couple the velocity components 1 5 2 The Thermodynamic Equations The model is forced by atmospheric heat fluxes as well as fresh water fluxes given by precip itation and evaporation Two prognostic variables for the heat content of snow and ice are introduced the temperature between snow and ice T and the temperature between snow and the atmosphere T The internal temperature profile within the ice and the snow is idealized to be linear The underlying approach is to convert the difference of the heat flux through the ice and snow layer into an accumulation of the heat content of both the snow and ice layers Furthermore the heat flux at the snow surface is set equal to the heat flux through the snow layer in the case that no heat is used to melt snow The system of equations that results is as follows pih O q y Pisi o kiqi ksqips T T Tr T h 1 Sipi aon s Cp
27. B grid L represents land points S represents sea points Circles denote scalar points crosses denote vector points The thick line marks the boundary at a certain time level The dashed line shows an example of how the boundary changes if a vector point is switched from a sea to a land point at the edge of a zero layer marked with the subscript V K I J for a 3 d array 2 d arrays are similarly marked by S T J or V T J Vector quantities are defined on vector points while scalar quantities are defined on scalar points If a scalar quantity is marked with the subscript V K J it means that the value is obtained by averaging the surrounding scalar values onto the vector point Thus the location of the vector point is located in the center of the cell that is surrounded by the scalar Al points S 1 J S T 1 J SU J 1 and S I 1 J 1 This means that if a scalar quantity P is needed on a vector point it is computed by Ps Psqrrr Psa J41 Pscr4s s41 Pirn 7 2 1 The analog for a vector quantity needed on a scalar point is Pyan Pya Pra Praia Poun 2 7 I 1 J 7 1 J 1 1 1 J 1 2 2 The model coordinates X VX for the x direction and Y NY for the y direction are defined on vector points Thus the array X contains the x values on the equator The local grid distance on the sphere is obtained by multiplying Az on the equator by cos yp see array SCALEX NY which also is defined on vector points
28. Command Pipe ses Figure 2 2 Concept of Server Client Interaction In order to improve convergence of various iteratively solved equations the starting latitude is interchanged on subsequent processors in longitude A further technique to speedup conver gence is to shift the data relative to the processor address space This data shifting is setup such that the amount of sent data is minimized The latter technique allows to resign on the solution of linear equation systems across processors Therefore there exist no sequential operations Because the data transfer is minimized the code scales well up to many processors Even in high resolution experiments there has not been observed a significant reduction of convergence speed A further trick is to perform the iterations on alternating latitudes with increasing pe_z as shown in Fig 2 4 In order to speedup the data postprocessing i e the transfer of data to disk the data are gathered in x direction on 1 pe y and then written onto files A signal is sent to the server which contains the amount of data already written to disk Then the server reads these distributed data and merges them onto one file in PPSF format Data that finally are saved on file OCDAT are returned to the servier via the parallel pipes and written to disk by the server Finally there are no data files left on disk which are needed later Due to the small network time and the negligible amount of code sections that hav
29. DIMENSION ISST 180 SST 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 ISST I I 1 180 DO I 1 180 SST I J MONTH ISST I 100 273 16 ENDDO ENDDO ENDDO e File AIRGLOB contains the monthly mean surface air temperature on a 2 x 2 grid The data are merged from the COADS whenever possible from the data of Shea 1986 for the Arctic and the Taljaard data for the southern hemisphere It is read by lt ATTEMP gt from the local file AIRTEMP with the following simplified statements DIMENSION TAIRT 180 AIRTEMP 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 IAIRT 1 I 1 180 DO I 1 180 AIRTEMP 1 J MONTH IAIRT 1 100 60 273 16 ENDDO ENDDO ENDDO e File COVER contains the monthly mean fractional cloud cover from COADS on a 2 x 2 erid It is read by lt CUMULUS gt with the following simplified statements DIMENSION ICLOUDS 180 CLOUDS 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 ICLOUDS I I 1 180 DO I 1 180 CLOUDS I J MONTH ICLOUDS 1 80 ENDDO ENDDO ENDDO e File WETNESS contains the monthly mean relative humidity from COADS on a 2 x 2 erid It is read by lt VAPOR gt with the following simplified statements DIMENSION IVAPOR 180 VAPOR 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 C IVAPOR I I 1 180 DO I 1 180 VAPOR I J MONTH IVAPOR I 1000 ENDDO ENDDO ENDDO e File PRESURF contains the monthly mean sea level pressure fro
30. Diffusion o 68 4 10 Parameters for Mixed Layer o e 69 4 11 Parameters for Horizontal Grid Definition 69 4 12 Threshold Variables 0 e 69 4 13 Parameters for Vertical Grid Definition 70 4 14 Definition of Eliminating Points 70 4 15 Grid Tunime Witches cua sedan ee bp ae ee be ed e 70 4 16 Parameters for Quick Look System o val 4 17 Switches for Listing 2 ote ek SA EES Bee BEES EM Se HR 71 4 18 Switches for Type of Forcing aaa aa aaa a a 72 4 19 Switches for Model Physics aaa aaa aaa a 72 4 20 Switches for Numerical Schemes aaa aaa aaa a 72 21 Switches for Output rec a BA a Be OG OES ORE eS 73 4 22 Switches for Data Sources 0 0 0 02000000000 eee 73 4 23 Time Weights for Implicit Scheme o o 73 4 24 Numerical Tuning Variables 0 13 4 25 Threshold Variables 4 4 2 4 4 Re a ee a ada 74 5 1 Definition of Quantity Code Prognostic Variables 126 5 2 Definition of Quantity Code Level Quantities 126 5 3 Definition of Quantity Code Diagnostic Variables 127 5 4 Definition of Quantity Code Forcing Variables 127 5 5 Definition of Quantity Code Flux Variables 128 5 6 Definition of Quantity Code Transport Variables 128 5 7 Definition of Qu
31. J M Oberhuber 1994a Interdecadal variability of the Pacific Ocean model response to observed heat fluxes and wind stress anomalies Clim Dyn Vol 9 287 302 Miller A J D R Cayan T P Barnett N E Graham and J M Oberhuber 1994b The 1976 77 Climate Shift of the Pacific Ocean Oceanography Vol 7 21 26 Oberhuber J M 1995 Parallelization of an ocean general circulation model on the CRAY T3D system CRAY CHANNELS Vol 17 No 2 22 25 Holland D M R G Ingram L A Mysak and J M Oberhuber 1995 A numerical simulation of the sea ice cover in the northern Greenland Sea J Geophys Res Vol 100 No C3 4751 4760 Oberhuber J M and K Ketelsen 1995 Parallelization of an OGCM on the CRAY T3D Proceedings of the Sixth ECMWF Workshop on the Use of Parallel Processors in Meteorology Title Coming of Age Editor G R Hoffmann and N Kreitz World Scientific Publishing Aukrust T and J M Oberhuber 1995 Modeling of the Greenland Iceland and Nor wegian Seas with a coupled sea ice mixed layer isopycnal ocean model J Geophys Res Vol 100 No C3 4771 4789 Duffy P B D E Eliason A J Bourgeois and C C Covey 1995 Simulation of bomb radiacarbon in two global ocean general circulation models J Geophys Res Vol 100 No C11 22 545 22 563 Oberhuber J M 1996 Surface flux fields for ocean models WCRP Report WMO TD No 762 53 58 Holland D M L A Mysak and J M Oberhuber 1996a An
32. Range Weather Forecasts Reading UK 72pp Hellermann S and M Rosenstein 1983 Normal monthly wind stress over the world ocean with error estimates J Phys Oceanogr 13 1093 1104 Hibler W D III 1979 A Dynamic Thermodynamic Sea Ice Model J Phys Oceanogr 9 815 846 Kleeman R and S B Power 1995 A simple atmospheric model of the surface heat flux for use in ocean modeling studies J Phys Oceanogr 25 92 105 Kraus E B and J S Turner 1967 A one dimensional model of the seasonal thermocline Tellus 1 88 97 Large W G and S Pond 1981 Open Ocean Momentum Measurements in Moderate to Strong Winds J Phys Oceanogr 11 324 336 Large W G and S Pond 1982 Sensible and Latent Heat Flux Measurements over the Sea J Phys Oceanogr 12 464 482 Legates D R and C J Willmott 1990 Mean seasonal and spatial variability in gauge corrected global precipitation Internat J Climatol 9 111 127 Martin P J 1985 Simulation of the Mixed Layer at OWS November and Papa With Sev eral Models J Geoph Res 90 903 916 Mikolajewicz U and E Maier Reimer 1994 Mixed boundary conditions in ocean general circulation models and their influence on the stability of the model s conveyor belt J Geophys Res 99 22633 22644 Millero F J 1978 Freezing point of seawater Eighth Report of the Joint Panel on Oceano graphic Tables and Standards Unesco Tech Pap in Mar Sci No 28 Annex 6
33. Since the mixed layer depth MLD is not only influenced by local mixing but also by horizon tal convergence of mass or heat the ML model invokes the full dynamics of equations 1 1 to 1 4 combined with a parameterization for the vertical transfer of mass and related quantities across the ML base The ML model always has nonzero thickness and carries an arbitrary instantaneous potential density c 1 2 2 The Mixed Layer Model Equations Entrainment and detrainment are treated differently While entrainment enters the continuity equation i e the prognostic equation for the layer thickness as a transfer rate between two adjacent layers detrainment is treated diagnostically The equation for the entrainment rate w is whimaz g g wRierit uy 014 01 014 m 2m mu 1 18 1 D hmams B yB Tmom B 1 T ymomsB h 1 exp h hg 2hg 1 exp h hg 12 Ma MoTafly mah where g yA TI gT 1 19 O61 g CpP B 8B P E R 1 20 ag AEP B R 1 20 ag B 0 770 do 1 21 g is the reduced gravity between the mixed layer and the next layer below g a threshold for Y Riesi the critical Richardson number B the total buoyancy flux through the surface Q the total heat flux into the mixed layer P E R denotes precipitation minus evaporation plus river runoff B the buoyancy flux due to the solar radiative heat flux Qs bot amp and are the analytically determine
34. Tide simulated with a T106 equivalent version of the barotropic tide model coupled to the 3 dimensional PIPE model help to improve regional simulations In order to simulate mean flows that are significantly modified as result of nonlinear interaction with tides a barotropic tide model was developed If Gide 18 the anomalous sea surface elevation Utide and Uride the respective zonal and meridional flow anomalies H gt gt hy the mean ocean depth as simulated by the 3 dimensional model and Frae the gravitational forcing then the equations are written in flux form OU e H ide 0 OttidelH Cria gp Guide Bride Gride f ride H ride Ot Cptttide Utide Un Utide UN 1 116 On e H ade tid F G d gt g H ride 5 Bride Cride f Utiae H Cride y CpUrndey Utide Un vride UN 1 117 0 0 By Stide i fee ride E Cride By Via A Cride 1 118 where wy and vy are the respective velocities in the lowermost present layer of the isopycnal model Because the tide model typically uses a much smaller time step than the 3 dimensional model the tidal flow anomalies are averaged in time and added as additional flow to the bottom friction parameterization This residual flow is also used in all transport equations in the 3 dimensional ocean model as for temperature salinity and momentum Because typical tidal frequencies are in the range of a time step of the 3 dimensio
35. UNESCO Paris Niiler P P 1975 Deepening of the wind mixed layer J Mar Res 33 405 422 Niiler P P and E B Kraus 1977 One dimensional model of the seasonal thermocline The Sea Vol VI Wiley Interscience 97 115 North G R J G Mengel and D A Short 1983 Simple Energy Balance Model Re solving the Seas and the Continents Application to the Astronomical Theory of Ice Ages J Geophys Res 88 6576 6586 Oberhuber J M 1986 About Some Numerical Methods Used in an Ocean General Circu lation Model with Isopycnic Coordinates Advanced Physical Oceanographic Numerical Modelling NATO ASI Series Series C Mathematical and Physical Sciences Vol 186 511 522 Oberhuber J M 1988 An atlas based on the COADS data set The budgets of heat buoyancy and turbulent kinetic energy at the surface of the global ocean Max Planck Institute for Meteorology Hamburg Report No 15 199pp Oberhuber J M 1990 Simulation of the Atlantic Circulation with a coupled Sea Ice Mixed Layer Isopycnal General Circulation Model Max Planck Institute for Meteorol ogy Hamburg Report No 59 86pp Oberhuber J M 1993a Simulation of the Atlantic Circulation with a Coupled Sea Ice Mixed Layer Isopycnal General Circulation Model Part I Model Description J Phys Oceanogr Vol 23 No 5 808 829 Oberhuber J M 1993b Simulation of the Atlantic Circulation with a Coupled Sea Ice Mixed Layer Isopycnal General Circu
36. and for JBAND and YCOM depend on each other 5 2 5 Tuning Topography and Coastline If the model initializes itself during the first model run see below there will possibly be sea points that should be declared as land points and vice versa With the switches KSW1 KSW20 one is able to treat the standard problems like closing the Panama channel etc These switches activate pre defined land bridges or channels Undesired lakes or decoupled parts of the world ocean in a regional model can be deleted by the eliminating points that define array indices The procedure is such that all sea points around an eliminating point are changed to a land point until all sea points within a connected area are redefined If the eliminating points are improperly chosen it can happen that all sea points vanish That will result in an overflow during divides by the number of sea points Therefore the first run to initialize the model should be performed with KSW1 KSW20 and KILL1I KILL6J set to zero Then these switches have no effect An explicit way to change the land sea mask is to put statements into lt FLAGMOD gt that redefine certain array elements of IFLG which carries the land sea mask on vector points An alternative option to control the assignment of land or sea to a grid point is to save the land sea mask during the first initialization For this use SW45 1 Then the file LS5MASK which contains the land sea distribution on scalar points can be edit
37. are plotted for better orientation The picture is created interactively using plot b p e q 69 PPSF_SUR Chapter 6 Appendices 6 1 Appendix A Prognostic Pressure Equation In order to present the strategy to solve the model equations with an implicit time integration scheme the momentum and pressure equations are written in non dimensional form If the pressure gradient and the flux divergence are equally weighted in time the finite difference representation in space is At At utt y TD v a LPP 6 3 At At PA ey 6 4 where Dj and D denote the algebraic representations of the gradient on vector points and divergence on scalar points respectively At is the time step and the superscripts l and 1 1 mark the old and new time levels The boundary conditions are introduced via a flag T that is PT 0 for land and T 1 for the ocean This implies no slip coastal boundary conditions Because Eqn 6 1 and 6 2 are interconnected by the pressure gradient and the flux divergence one of the two options is to remove u t in the pressure equation This results in an equation for the pressure p only at the new time level At At At pt Es p Peel Dap t p pS gt Da tl p Pee T Dap 6 5 After having found the solution for p t it is used to calculate the flux u at the new time level This simple strategy to obtain a stable solution for any chosen time step is used to derive a wave equation for
38. be done to create a model forcing by hand except for exotic demands 5 4 3 1 Data Bank Next it is briefly outlined how to read all the data files The FORTRAN statements neither contain any interpolation to the model grid nor do they show how undefined data points are treated 5 4 3 1 1 Ocean Initialization e File TOP1DEG contains the so called Scripps topography on a 1 x 1 grid Extra latitudes for the South and North Pole have been included It is read by lt READTOP gt from the local file TOPOG with the following simplified statements DIMENSION IDEPTH 360 DEPTH 360 180 DO J 1 180 READ 1X 12 15 1X CIDEPTH I I 1 360 DO I 1 360 DEPTH I J HOTAL IDEPTH T ENDDO ENDDO e File TOP5MIN contains the NOAA topography which is based on a 5 minutes res olution It is read by lt READTOF gt from the local file TOPOG with the following simplified statements DIMENSION IDEPTH 4320 DEPTH 4320 2160 DO J 1 2160 READ 1216 CIDEPTH I I 1 4320 DO I 1 4320 DEPTH I J MONTH HOTAL IDEPTH TI ENDDO ENDDO e File SALTEMP contains the monthly mean temperature and salinity of Levitus for the global ocean on a 1 x 1 grid It is read by lt INTERPO gt with the following simplified statements DIMENSION ITEMP 32 ISALT 32 TEMP 360 180 32 SALT 360 180 32 DO MONTH 1 12 DO J 1 180 DO I 1 360 READ 3214 ISALT K K 1 32 READ 3214 ITEMP K K 1 32 DO K 1
39. coastline geometry of some rotated model domain This code is usefull to optimize the Eulerian angles The resulting object file must be loaded together with jmoplane o and utility o as shown in bin paint sea Surface Velocity January Mean l Z IR ANNA Fe Ta a i Sree aN A AA ACC x ub PENS f Y y zZ o Ne Z o on SEs EOL gt gt y Ze he ES ELTZL A AA vt N r ef y AOS lt MG EOS bh ELL y Ne KNSS SSS z See ery Py oT Nx AA ere Ls N aa Ke SZ ZS mee nro Do raro h RP GRA Y R PT 7 A PAS Fe fogs Y ES A A LAG LI ETS SLL S y a IAN OO O NA Cada Y MN IMA S T T T T T T T T T T T pe P 30E 60E 90E 120E 150E 180 150W 120W 90W 60W 30W 0 30E COMMENTS 3 8 cm s t 3 8 Josef M Oberhuber DKRZ Wed Mar 10 14 42 15 1999 C JMOPlane Figure 5 1 Example showing a vector plot for the surface flow of PIPE T42 The picture is created interactively using plot b e q 69 a Annual Mean PPSF_SUR Sea Surface Elevation nnual Mean 60N 90 30N 305 605 905 T T T T 30E 60E 90E 120E 150E 180 150W 120W 90W 60W 30W 0 30E T T T T T T T T T T l 210 190 170 150 130 110 90 70 50 30 10 10 30 50 70 90 110 130 150 Unit cm Josef M Oberhuber DKRZ Wed Mar 10 16 07 53 1999 C JMOPlane Figure 5 2 Example showing a greyscale colour plot for the sea level of the PIPE T42 The picture is create
40. defined on a xz section lt DRUCKHY gt prints out a scalar quantity defined on a yz section lt DRUDATX gt prints out observed ocean data defined on a xz section lt DRUDATY gt prints out observed ocean data defined on a yz section 4 8 Parallel Code oceutil f 4 8 1 Pre and Postprocessing lt MPPSEND gt is the main program that preprocesses the input data and organizes the preparation of the subsequent execution of the parallel code lt MPPSTEP gt is the main program that serves to provide the parallel job with data while 1t is running lt MPPSTOP gt is the main program that postprocesses the data left by the parallel job on the disk lt t3d_main gt is the main program of the parallel job lt t3d_setup gt sets up the topology of the parallel job lt RECEIVE gt as part of the parallel job has the task to receive the data from the server and to distribute the data to the respective processors lt SEND_C90 gt as part of the parallel job has the task to collect the data and return them to the server 4 8 2 Data Transfer 4 8 2 1 Routines on the Server Side lt SPIPE_OPN gt opens the command pipe The command pipe is used to synchronize the tasks on the server and the client Entry lt SPIPE_CLS gt closes the command pipe e lt SERVER_OPN gt opens the data pipes Entry lt SERVER_CLS gt closes the data pipes e lt SERVER_GET gt receives real data from the client Entry lt SERVER_PUT gt
41. for UNICOS CRAY or SUNOS operating systems Before using them adjust the directory path names to your account The script files are 1 setup is required to generate the code for a specific computer with the help of the precompiler precomp f and changes the initial model dimensions to the desired values for NX NY NZ NXOLD and NYOLD For the latter two parameters see also Multi Grid Method 2 make sun or make cri generates the object files for ocestep f oceproc f ocepost f ocemods f and oceplot f and saves them in the prescribed directory The script con tains switches for how the FORTRAN codes should be compiled e g for autotasking debugging or performance tests 3 model sun or model cri runs the model A number of switches control e g which data have to be used for initialization whether the model is started from scratch which output files are saved etc 4 plot is an example for how to create plots It plots data written onto one of the postprocessing files whose names start by default with PPSF_ 5 4 1 How to Generate an Executable Model Version The grid dimensions NX NY NZ NXOLD and NYOLD have to be determined before the final FORTRAN source code can be compiled In the next step setup creates the FORTRAN code The batch job model sun or model cri requires at least object files of ocestep f oceproc f ocepost f and optionally ocemods f These files are created by make sun
42. grid points yields a non elliptic equation which can be solved nevertheless with an appropriate underrelaxation scheme 6 3 Appendix C Constants for Equation of State 0 36504 aro 0 017439 0 083198 ar11 0 00029778 0 89306 a210 0 0041057 0 031626 0 00054065 a292 0 00021987 0 0000040274 0 16056 0 0050484 Table 6 1 Constants for Bryden s Formula The subscripts denote the exponent for pres sure first index salinity second index and temperature third index See also Eqn 1 16 ao 999 842594 8 24493 101 19652 21 54 6746 a 0 06793952 4 0899 1073 148 4206 0 603459 a gt 9 09529 1073 7 6438 1075 2 327105 0 0109987 az 1 001685 1074 8 2467 1077 0 01360477 6 167 1075 aa 1 120083 1074 5 3875 107 5 155288 107 as 6 536332 107 5 72466 107 7 944 107 3 239908 1 0227 1074 1 6483 107 1 43713 107 1 6546 1076 5 3009 1074 1 16092 1074 5 77905 1077 d ie 107 7 190 10 m 9 9348 107 i 2 2838 107 ky 8 50935 107 mo 2 0816 107 41 1 0981 107 ky 6 12293 1076 m3 9 1697 107 ig 1 6078 1076 k 5 2787 1078 Table 6 2 Constants for UNESCO Formula The subscripts denote the exponent of re spective variables See also Eqn 1 17 Part IV References 147 Chapter 7 How to Find References 7 1 Code References Name in Code Name in Formulae Meaning Unit Equation UH zonal mass flux VH meridional mass flux HEIGTH layer thickness TEMP potential temperature SALT salini
43. gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt lt MAKEXY gt Reference lt ROMATIN gt lt ROMATIN gt lt ROMATIN gt lt ROMATIN gt Variable TEMPTOP any TEMPBOT any SALTTOP SALTBOT Value any any Default 299 16 273 16 35 00 34 90 Unit Kelvin temperature for top layer Kelvin temperature for bottom layer psu salinity for top layer Meaning Reference lt MAKEVER gt lt MAKEVER gt lt MAKEVER gt lt MAKEVER gt Variable KILL1I KILL1J KILL2I KILL2J KILL3I KILL3J KILL4I KILL4J KILLSI KILL5J KILL6I psu salinity for bottom layer Table 4 13 Parameters for Vertical Grid Definition Value any index any index any index any index any index any index any index any index any index any index any index Default Meaning x index for eliminating point 1 y index for eliminating point 1 x index for eliminating point 2 y index for eliminating point 2 x index for eliminating point 3 y index for eliminating point 3 x index for eliminating point 4 y index for eliminating point 4 x index for eliminating point 5 y index for eliminating point 5 x index for eliminating point 6 Reference lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN lt ELIMIN OOOO OS OV OS OOS VVVVVVVVVV VV KILL6J any index y index for eliminating point 6 lt ELIMIN T
44. have been rewritten The net longwave radiation at the surface is parameterized as Qi bot e0 T To 1 39 0 4 A E es E 5s ye 1 40 ff E 100 0021 ren 1 40 where e is the emissivity of water o the Stefan Boltzmann constant and n the relative cloud cover The insolation is calculated from the daily averaged heat flux at the top of the atmosphere and is then corrected after Zillmann 1972 for relative humidity and inclination The insola tion is reduced by a combined cloudiness and humidity correction The resulting relations to compute the daily mean downward shortwave radiative flux Qs bot are cosy sindsiny cos cos y cost 1 41 SiN noon sin sin Y cos cos p 1 42 1 0 0 001 99100 1 A A 1 43 1 0 0 00191egua 1 0 0 79 1 0 0 0004re 2 n2 ky l w fte S cos d PATE a E a 1 44 277 t cos y 2 7 r5 1 085 cos y 0 1 d where S is the solar constant 7 the solar elevation Nequa 90 deg and w the albedo d denotes the distance between the sun and earth and d its annual average Following Paltridge and Platt 1976 the ratio d d is estimated in terms of the Julian Day for the present day orbit d Pi 1 00011 0 00128 sin 5 0 034221 cos 0 000077 sin 28 0 000719 cos 28 1 45 The declination in radians which is needed to compute 7 is given by 0 006918 0 070257 sin 8 0 399912 cos 8 0 000907 sin 2 0 006758 cos 28 0 00148 sin 38 0 0
45. in order to interrupt the execution In this situation the conclusion is that the model had problems in finding clean solutions In other words the algorithms are overcharged with the combination of physical situation and layout of model grid choice of time step external forcing etc Model Time Threshold Reached If the model time counter has reached the model time as defined by NHOUR NDAY NMONTH and NYEAR then the model stops independent on the setting of MAXREP If one tries to continue the run then the routine lt SUICIDE gt is called which creates a floating point exception If you want to continue the run increase the model time thresholds Other Floating Point Exceptions So far the only known possibilities are that either the model time threshold is reached however lt SUICIDE gt was not called or that a numerical instability appeared In the first case the time step DT becomes zero which results in a divide by zero in subsequent routines The first routine then called is lt CROSMIX gt In the second case the routine lt SURFLUX gt is the first one called during a time step which seems to respond sensitively to e g erroneous temperatures Stop statement If an unexpected stop statement is executed follow the instructions This is a way to educate people that it is worth to read manuals These stop statements appear in lt OCEINIT gt in order to issue nonsense switch or parameter settings and in lt RALEIGH gt and lt BAROT
46. investigation of the gen eral circulation of the Arctic Ocean using an isopycnal model Tellus Vol 48A 138 157 Holland D M L A Mysak and J M Oberhuber 1996b Simulation of the mixed layer circulation in the Arctic Ocean J Geophys Res Vol 101 No C1 1111 1128 Cherniawsky J Y and J M Oberhuber 1996 The seasonal cycle of mixed layer tem peratures in a global ocean general circulation model Clim Dyn Vol 12 171 183 Lunkeit F R Sausen and J M Oberhuber 1996 Climate simulations with the global coupled atmosphere ocean model ECHAM2 OPYC Clim Dyn Vol 12 195 212 Roeckner E J M Oberhuber A Bacher M Christoph and I Kirchner 1996 ENSO variability and atmospheric response in a global coupled atmosphere ocean GCM Clim Dyn Vol 12 737 754 Stossel A J M Oberhuber and E Maier Reimer 1996 On the representation of sea ice in global general circulation models J Geophys Res Vol 101 No C8 18 193 18 212 Kauker F and J M Oberhuber 1997 An Isopycnal Ocean Circulation Model of the North Sea for Dynamical Downscaling GKSS 97 E 47 GKSS Forschungszentrum GmbH Geesthacht Tyler R H L A Mysak and J M Oberhuber 1997 Electromagnetic fields generated by a three dimensional global ocean circulation J Geophys Res Vol 102 No C3 5531 5551 Tyler R H T B Sanford and J M Oberhuber 1997 Geophysical Challenges in Us ing Large Scale Ocean Generated EM Fields to Deter
47. lt NEWDENS gt computes the in situ density the potential density and the in situ pressure from potential temperature and salinity lt REPRESS gt computes the in situ pressure lt STATETR gt computes the in situ density from temperature salinity and in situ pres sure for one layer only Entry lt STATETP gt performs the same calculation for potential density lt STATEVR gt computes the in situ density from temperature salinity and in situ pres sure for all layers Entry lt STATEVP gt performs the same calculation for potential density e lt THERMO gt computes the layer thickness changes at a given mass content tempera ture and salinity This routine is required to update the expansion effects while solving the wave equation lt TTOTET gt computes the potential temperature from temperature for one layer lt TETATOT gt computes the temperature from potential temperature by inverting the formula of Bryden for all layers 4 3 8 Model Forcing lt SURFLUX gt computes the heat fluxes fresh water fluxes and buoyancy fluxes from the atmospheric data the model SST and model sea surface salinity Depending on the switches either simple Newtonian type parameterizations or more detailed one are taken It also contains the calculation of the fluxes through ice and snow and computes by iteration the snow skin temperature as well as the snow ice interface temperature as prognostic variable Note that this routine is responsible f
48. lt SYSEC gt prints out scalar field e lt VZSEC gt writes out a one dimensional vector in z direction for certain level Entry lt SZSEC gt prints out scalar field e lt VXYSEC gt writes out an internally computed xy section at a given depth Entry lt SXYSEC gt prints out scalar field e lt VXZSEC gt writes out an internally computed xz section at a given latitude Entry lt SXZSEC gt prints out scalar field e lt VYZSEC gt writes out an internally computed yz section at a given longitude Entry lt SYZSEC gt prints out scalar field 4 5 5 Tridiagonal Solvers e lt TRIDIAX gt computes the solution of a single linear equation of tridiagonal type for NX number of equations It is used by lt TIDES gt e lt TRIDIAY gt computes the solution of a single linear equation of tridiagonal type for NY number of equations It is used by lt TIDES gt e lt TRIDIA2 gt computes the solution of two independent linear equations of tridiagonal type It is used by lt EBM gt to compute the solution of the two tridiagonal systems for temperature and moisture e lt TRIDIA3 gt computes the solution of three independent linear equations of tridiagonal type It is used by lt SEAICE gt to compute the solution of the three tridiagonal systems for ice thickness ice concentration and snow cover e lt TRIDIAV gt computes the solution of NV NZ linear equations of tridiagonal type for NX number of equations It is used by lt TRACTIV gt
49. mean fresh water PMEM Frias fresh water flux bias correction TAUXQ Ta ti running mean zonal wind stress TAUXM Ta bias zonal wind stress bias correction TAUXQ Totun running mean meridional wind stress Nm PMEM Ty bids meridional wind stress bias correction Nm FACDEC run relaxation factor Table 7 4 Names for Bias Correction Terms Name in Code Name in Formulae Meaning Unit Equation UICE zonal flux of ice VICE meridional flux of ice HICE ice thickness COMPACT ice compactness HSNOW snow depth TICE ice temperature TSNOW snow temperature DMELT melting of ice amp snow per At Table 7 5 Names for Prognostic Snow and Sea Ice Variables Name in Code Name in Formulae Meaning Unit Equation VERTIC entrainment rate HEQU Monin Obukhov length ENERGY turbulent kinetic energy TAUX Garwood term Table 7 6 Names for Diagnostic Mixed Layer Variables Name in Code Name in Formulae Meaning Unit Equation DIFFUSV Ay factor for momentum diffusion DIFMINV threshold for momentum diffusion DIFFUSS A factor for scalar diffusion DIFMINS threshold value for scalar diffusion DIADIFO A 0 diffusion for layer thickness DIADIF1 AM factor for layer thickness diffusion TURBLEV mous TKE for vertical mixing HMIXV hy reference thickness DRAGT time constant for salinity forcing DRAGS Ca drag coefficient at surface DRAGB Ca drag coefficient at bottom Table 7 7 Names for Horizontal and Vertical Diffusion Paramete
50. numbers are documented in OUTLIST as in the example that shows five time steps The numbers appearing in rows 1 6 denote used iterations within the major iterative problems The smallest number that the model has to take is defined by NITMIN Each of the algorithms decides by itself whether more iterations are necessary The required accuracy can be adjusted with the variables ACCUR If the maximal number of iterations defined by NITMAX is exceeded the model diagnoses a numerical problem Depending on the switch SW35 PIPE either creates an error exit or stops the code without saving any state Note that the shown numbers always are half of what is really used as number of iterations Row 7 shows the time step that the model has chosen If DT MAX DTMIN then the model does not try to vary the time step however if DTMIN lt DTMAX then the model searches that maximal time step that still ensures a sufficient convergence of all iterative algorithms In the case that the model suddenly blows up the model reduces the time step and retries the same time step once more in the hope that the convergence conditions are satisfied by the reduced time step Row 8 is the average of the absolute value of the acceleration over each grid point It is a sensitive measure for numerical noise Any instability will immediately appear as big number These are slightly problem depdendent but must never exceed values of about 5000 The following two rows 9 and 10 show t
51. ocestep f contains all routines to carry out the time steps oceproc f those routines required to initialize the model and the final operations of the model run and ocepost f contains all routines for the postprocessing ocemods f optionally contains routines of ocestep f oceproc f or ocepost f that require some specific changes and overwrite the original version of those routines while the code is loaded In general ocestep f oceproc f and ocepost f contain only the dimensions as specification for the model setup All parameters are defined in one of the BLOCK DATA statements This has the advantage that as long as the model dimensions are unchanged all executable routines have to be compiled only once Any change to one of the executable routines can be put into ocemods f Any change of switches or physical parameters are implemented in ocemain f as it contains all BLOCK DATA statements Therefore only ocemain f and optionally ocemods f have to be compiled before each model run 4 2 Contents of ocemain f 4 2 1 Main Program MY_OGCM The main program consists of four calls to the routines lt OCEINIT gt lt OCESTEP gt lt OCEPOST gt and lt OCESTOP gt When the model is coupled to an AGCM the main pro gram also contains calls to the atmosphere and to routines that carry out the communication between both systems From the view point of the main program the ocean appears as a black box Any data exchange with t
52. perform a model initialization 7 Run Source RISC Model model mpp 8 Adjust Source RISC Model model mpp to perform a continuation run 9 Run Source RISC Model model mpp In order to setup an arbitrary version for a MPP computer there are rules to setup the dimen sions During the installation the file Source parammpp h is inserted into the code through an include statement The dimensions used on each node are the result of the parameters 1 nz_tot which is the full dimension in x direction of the model domain 2 npe_x which is the number of processors used in x direction 3 ny tot which is the full dimension in y direction of the model domain 4 npe_y which is the number of processors used in y direction The dimensions nz and ny used for each node are then computed through nz nx_tot 2 npe_x 2 and ny ny_tot 2 npe_y 2 inside parammpp h It must be guaranteed that nx_tot 2 npe_x and ny ny tot 2 npe_y has no remainder In addition it must be guaranteed that the variable NPIPE in the script pipe is equal npe_y That s because each processor in y direction must be connected via UNIX pipes with the server In general npe_z and npe_y with the product npe_x npe_y equal the number of used processors should be chosen in such a way that npe_r npe_y in order to minimize data traffic If this is not possible choose npe x lt npe_y which results in better convergence in x direction because of a smalle
53. relative changes become too big at one grid point The procedure to define a focus is the following 1 Define the model margins and compute the required dimensions to obtain a certain grid distance see 5 2 3 This yields a grid distance that will be used outside the focus area 2 Choose the number of grid points that will additionally be added to the grid in the region of higher resolution the focus area The final dimensions are XGR XGL NX LA 2424 DENS 5 3 NY KORYO 42428 JDENS 5 4 y where now Az and Ay are the grid distances outside the focus and IDENS and JDENS are the number of grid points used to increase the resolution on both sides of the focus center In the next step the index for the focus center has to be determined If XGC and YGC are the x and y locations of the focus center ICENTER and JCENTER are given by XGC XGL ICENTER 414 IDENS 5 5 T YGC YGU JCENTER EZ 1 JDENS 5 6 y If AX and AY are the total widths of the focus in the x and y direction respectively then the widths in terms of index range are given by 2x AX PAND GR KGL ae 2x AY BAND id YGO Y GU 2 8 3 Choose the refinement factors XCOM and YCOM for the x and y direction respectively These represent the ratio between the grid distance outside the focus and in the focus center The focus is switched off if IDENS JDENS 0 IBAND JBAND 0 and XCOM YCOM 1 Note that the choices for IBAND and XCOM
54. scale quantity JMOADFA FACTOR any multiplicative constant to scale quantity JMOADFA NDEC any switch for selecting the label type JMOADFA Table 5 15 Definition of Scaling Parameters 5 6 1 10 Block Data QZONMER lt QZONMER gt contains definitions for additional curves for means in x direction or and y direction and for ignoring a contour with some specified index see also table 5 16 Variable Value Meaning Reference NMERID any number of contour colour intervals shown as x curve ZONMER NZONAL any number of contour colour intervals shown as y curve ZONMER IGNORE any index of contour which is not plotted ZONMER Table 5 16 Definition of Zonal and or Meridional Means 5 6 1 11 Executing Program OCEPLOT lt OCEPLOT gt uses the plot library jmoplane f Basically it analyses the header of each data block and decides how a data set has to be treated The access to the key parameters is provided by the BLOCK DATA structures explained above 5 6 2 The Underlying Library JMOPLANE moplane f is a plot library developed by the author It is optimized for the purpose of plotting ocean quantities with land sea masks being easily treated For any detailed information there is an online manual available located on Graphic Manual which describes each routine It either expects some GKSOA level GKS library or writes the output in PostScript format if the option DNCAR is not used at compile time 5 6 3 Example Plot
55. sea ice SEAICE solve pressure equation DYNAMIC correct for mixed layer MIXIMP correct for momentum MODIMP TRACER CONVECT main program MY_OGCM Figure 3 3 Main structure of lt OCESTEP gt which is the driving routine for carrying out time steps update atmospheric forcing aa create quick look file PPSF_SUR create quick look file PPSF_ICE create quick look file PPSF_BAS create quick look file PPSF_SEC create quick look file PPSF_VER create quick look file PPSF_FLX create quick look file PPSF_FOR create quick look file PPSF_TID main program E MY_OGCM Figure 3 4 Main structure of lt OCEPOST gt which is the driving routine for data postprocessing and generation of quick look files final I O OCESTOP save restart file OCDAT WRITEM create listing LISTING create model statistics do final data postprocessing main program y MY_OGCM lt create CPU time statistics CPU time statistics Figure 3 5 Main structure of lt OCESTOP gt which is the driving routine for final calculations Chapter 4 Model Code Description 4 1 General Remarks The ocean code is divided up into 5 files which are ocemain f ocestep f oceproc f ocepost f and ocemods f ocemain f consists of a simple main program and all parameter definitions via BLOCK DATA statements
56. switch for output of monthly data lt OCESTEP gt Table 4 21 Switches for Output Variable Value Default Meaning Reference restart from other OCDAT lt OCEINIT gt H amp R or ECMWF stresses lt OCEINIT gt COADS ECMWF or Levitus SST data lt OCEINIT gt Scripps or NOAA topography lt OCEINIT gt initialization mode for FORCES and OBNDSEC lt OCEINIT gt COADS or ECMWE air temperature lt OCEINIT gt import export of land sea mask lt OCEINIT gt COADS or ECMWF winds lt OCEINIT gt Table 4 22 Switches for Data Sources Variable Value Meaning Reference ALPHA 0 75 time weight for flux divergence JMOBACK BETA 0 75 time weight for pressure gradient JMOBACK GAMMA 0 75 time weight for Coriolis force JMOBACK DELTA 0 75 time weight for momentum advection JMOBACK ETA 0 75 time weight for momentum diffusion JMOBACK Table 4 23 Time Weights for Implicit Scheme Variable Value Unit Meaning Reference OVER 1 25 overrelaxation coefficient JMOTUNE UNDERI 1 00 underrelaxation coefficient except for ice JMOTUNE UNDER2 1 00 underrelaxation coefficient for ice JMOTUNE NITMIN 3 minimal number of iterations JMOITER NITMAX 10 maximal number of iterations JMOITER ACCURU 0 1 kgs m required accuracy for momentum lt MODIMP gt ACCURH 0 01 m required accuracy for thickness lt SOLVER gt ACCURT 0 00002 K required accuracy for temperature lt TRACER gt ACCURS 0 00001 gkg required accuracy for sa
57. temperature and the 0Q OT out of the common data area e lt OADICE gt reads the ice thickness snow thickness and ice compactness out of the common data area e g for later heat flux computations in an atmospheric model 4 5 6 3 Internal Data Transfer The following routines supply the model with data from the common data area which is JMOCOUP for the fluxes and JMOBASE for the underlying observational data however they should not be used explicitly to feed the model with data Instead use the following routines e lt HEATMOD gt extracts the total surface heat flux from the common data area See lt SURFLUX gt e lt SOLMOD gt extracts the solar radiation from the common data area See lt SURFLUX gt e lt PMEMOD gt extracts the fresh water flux from the common data area See lt SURFLUX gt e lt TAUMOD gt extracts the wind stress components from the common data area See lt STRESS gt e lt MIXMOD gt extracts the turbulent kinetic energy for the mixed layer from the common data area See lt MIXEXP gt and lt MIXIMP gt e lt ICEMOD gt extracts the snow surface and snow ice interface temperature from the common data area See lt SURFLUX gt e lt SNOMOD gt extracts the heat fluxes that concern ice and snow out of the common block area See lt DIADICE gt e lt TKEMOD gt extracts the turbulent kinetic energy out of the common data area See lt MECHAN gt e lt AIRTMOD gt extracts t
58. term with FLG while for scalar points four surrounding values of IFLG have to be considered thus requiring detailed knowledge about the origin of each expression The interfaces in the isopycnic part of the model physically disappear either into the mixed layer where the isopycnals are vertical or into the topography e g at the continental shelf or at the sea floor Formally a layer that has physically disappeared is defined in the model as a layer with zero thickness Grid points in this layer do not contain mass but still hold dummy values for temperature salinity and other quantities An appropriate boundary condition is the key to decouple a zero layer from the ocean Furthermore physical processes that change the location where an isopycnal disappears reappears or which create water masses with a not yet existing potential density must have their counterpart in a technique that allows the shifting of boundaries or the flooding of zero layers This is explained in JMO Simply explained a finite mass cell remains part of the ocean as long as it contains mass and is switched off from the ocean only if it has already lost its entire mass and is still losing mass thus avoiding a negative mass content Massless cells become part of the ocean if convergence is predicted and remain massless if the flow is divergent 2 1 1 2 Operators The equations for momentum mass heat salinity and tracer concentration are formulated in spherical coor
59. that allows to add an arbitrary forcing to the model without the need to understand the time stepping 4 7 Open Boundaries oceopen f 4 7 1 Routines for Preprocessing Data lt INTCUT gt is the main driving routine to create the boundary data It is called by lt OCEINIT gt lt STREAM_T106 gt reads the barotropic transport from data and interpolates it onto the model grid lt SELECOBND gt reads monthly 3 dimensional ocean data for temperature and salinity from Levitus reads the barotropic transports that have been created by a T106 version of PIPE interpolates them onto the model grid and saves them on file OBNDSEC If the tide model is switched on data are read interpolated and saved onto OBNDSEC Therefore a file OBNDSEC is not interchangable between models versions that don t use and those which use a tide model although the model geometry and grid is identical lt WRITEOBND gt writes temperature salinity and transports along the southern northern western and eastern boundary of the model domain onto OBNDSEC lt INTIDES gt reads tidal data and interpolates them onto the model grid lt WRITETIDE gt writes amplitude and phases of various tidal modes onto OBNSEC lt READOBND gt reads current boundary data lt READOBNDOLD gt reads boundary data of last month lt STRDATI1X gt writes southern boundary of observed ocean data on a 3 dimensional array onto a 2 dimensional array containing an xz section
60. the color table for the type of projection chosen for all vector contour or color plots and a parameter which defines whether contours surround all equally colored areas see also tabel 5 13 Variable Value Meaning Reference IROT 0 1 switch for rotation of the color table JMOIND KEYCOL 1 2 3 4 switch to control the margins of color table JMOIND IPOL 0 1 2 3 switch to select a projection JMOIND ICON any gt 0 switch for contours in color plots JMOIND Table 5 13 Definition of Layout Parameters 5 6 1 7 Block Data MINMAXS lt MINMAXS gt contains the definitions for the lowest and highest level that is plotted by con tours or colored areas see also table 5 14 Variable Value Meaning Reference SMIN any lowest level used for contouring SMAX any gt SMIN highest level used for contouring DMIN any as SMIN but for negative section code DMAX any gt SMAX as SMAX but for negative code numbers IMON 2 1 0 1 2 switch for shading in contour plots IKOL any number of colours used for shading Table 5 14 Definition of Interval Parameters 5 6 1 8 Block Data CVARIAB lt CVARIAB gt allows to specify non equidistant contour or color intervals 5 6 1 9 Block Data ADDFACS lt ADDFACS gt contains definitions that rescale quantities e g to avoid long labels in a contour plot and a parameter that determines the layout of a label see also table 5 15 Variable Value Meaning Reference ADD any additive constant to
61. the hydrostatic approximation that prohibits its application down to eddy resolving scales Ignoring for the moment dissipative processes the basic quantities which should be con served are momentum energy mass potential vorticity heat and salt content Momentum mass heat and salt content are easily conserved if the flux form of the equations is chosen Use of the flux form is crucial otherwise conservation is hard to maintain in a moving coordinate system Potential vorticity and energy can also be conserved after some manipulation of the momentum transport and Coriolis terms See further discussion in JMO In brief the philosophy of PIPE is to 1 not assume that the model a priori uses isopycnals Therefore the model is allowed to deviate from isopycnal coordinates if it is necessary to avoid an otherwise occuring conflict in the concept 2 use two densities which is the in situ density allowing for mass conservation in con junction with a proper continuity equation and the potential density which is needed to formulate pressure gradients in the curved isopycnal coordinates and to control the coordinate system through mixing 1000 750 500 250 0 905 605 305 0 30N 60N 90N Figure 1 1 Meridional cross section of interface height for the example of a T42 global ocean simulation The section is located at 120W The interfaces are plotted as steps in order to show where grid cells are located Each step denotes the location o
62. the model s equations for the momentum and pressure 6 2 Appendix B Pressure Boundary Condition A subject of controversy is often the use of an elliptic equation for predicting the pressure It must be stressed that the analytical derivation of Eqn 2 40 after eliminating the velocities according to section 2 1 2 does not automatically provide the correct pressure boundary con dition In order to obtain a correct discretization for pressure points near the boundary first 145 a no slip boundary condition is applied for the momentum Eqn 1 1 and a no flux condition for the continuity Eqn 1 2 So far there is no difference to explicit time stepping models The discrete prognostic equations are then modified such that specific terms are labelled with the new time level index 1 rather than the old time level Then the prognostic equation for the pressure p t at the new time level is obtained using algebraic manipulations which do not change the solution at all and thus do not introduce inconsistencies in boundary con ditions This way of deriving a discrete equation for the new pressure provides the algebraic representation of the prognostic pressure equation near the boundary This equation replaces otherwise necessary and mostly intuitive diagnostic pressure boundary conditions that might lead to erroneous physical behaviour In analogy to the derivation of Eqn 6 5 the derivation of the model s pressure equation on boundary
63. the open boundary conditions I also appreciate the valuable comments by Arthur Miller at the initial state of the model and Andreas Bacher for his comments on numerous coupling aspects and related improvements 7 3 Literature References 7 3 1 Text References Bleck R and D B Boudra 1981 Initial testing of a numerical ocean circulation model using a hybrid quasi isopycnic vertical coordinate J Phys Oceanogr 1 755 770 Bryden H 1973 New polynomials for thermal expansion adiabatic temperature gradient and potential temperature of sea water Deep Sea Res 20 401 408 Denman K L and M Miyake 1973 Upper Layer Modification at Ocean Station Papa Observations and Simulations J Phys Oceanogr 3 185 196 Garwood R W 1977 An oceanic mixed layer model capable of simulating cyclic states J Phys Oceanogr 7 455 468 Garwood R W Jr P C Gallacher and P M ller 1985a Wind Direction and Equi librium Mixed layer Depth General Theory J Phys Oceanogr 15 1525 1531 Garwood R W Jr P C Gallacher and P Miller 1985b Wind Direction and Equi librium Mixed Layer Depth in the Tropical Pacific Ocean J Phys Oceanogr 15 1532 1538 Gaspar P 1988 Modeling the Seasonal Cycle of the Upper Ocean J Phys Oceanogr 18 161 180 Gibson J K O Kallberg S Uppala A Hernandez A Nomura and E Serrano 1997 ECMWF reanalysis project report series 1 ERA description European Centre for Medium
64. to lt TOUCHLW gt but uses the flag field IF LG to detect a physical layer Entry lt CONTUP gt works similarly to lt TOUCHUP gt but uses the flag field JF LG to detect a physical layer 4 4 12 Routines for Data Manipulation lt SELECT gt extracts a 2 dimensional real array out of a 3 dimensional one lt KELECT gt extracts a 2 dimensional integer array out of a 3 dimensional one lt SAMMELN gt extracts a 2 dimensional real array out of a 3 dimensional one with the vertical dimension of the number of data levels i e from Levitus Entry lt VERTEIL gt performs the inverse operation lt COLLECT gt extracts a 2 dimensional real array out of a 3 dimensional one for a variable vertical index lt VECMAX gt sets the 2 dimensional input field to the maximum of itself and a specified constant Entry lt VECMIN gt uses the minimum instead lt FELDMAX gt sets the 3 dimensional input field to the maximum of itself and a spec ified constant Entry lt FELDMIN gt uses the minimum instead lt STOREXZ gt stores xz section onto another array required for open boundaries lt STOREYZ gt stores yz section onto another array required for open boundaries e lt STORECI1X gt stores xz section onto southern boundary of 3 dimensional field re quired for open boundaries e lt STORECNX gt stores xz section onto northern boundary of 3 dimensional field re quired for open boundaries lt STOREC1Y gt stores yz section onto
65. western boundary of 3 dimensional field required for open boundaries lt STORECNY gt stores yz section onto eastern boundary of 3 dimensional field required for open boundaries 4 4 13 Routines for Filtering e gt RUBBS1 gt filters 2 dimensional scalar quantities Entry lt RUBBV1 gt filters 2 dimensional vetcor quantities e gt RUBBSNZ gt filters 3 dimensional scalar quantities Entry lt RUBBVNZ gt filters 3 dimensional vector quantities 4 4 14 Routines for Filling Up e lt FILLUP gt fills up arrays containing undefined data e lt FILLUP1 gt fills up arrays containing undefined data 4 4 15 Routines for Error Detection e lt INFORM gt outputs an information that certain data exceeded a specified range e lt CAUTION gt outputs an information that certain data exceeded a specified range e lt CHECK gt outputs an information that certain data exceeded a specified range e lt SUICIDE gt causes an error exit through floating point exception upon request 4 5 Postprocessing ocepost f The vertical coordinates in the model are defined by HTOTAL This is the depth on which the origin of the vertical coordinate is located Therefore the topography depth is zero at a depth of HTOTAL This determines how the depth values are stored in the 2d array DEPTH The sea level will vary around HTOTAL which is important to know if one converts the layer thicknesses and topography depth into interface depth values via lt HTOZ
66. with sections of temperature and salinity i e with the pressure field while leaving it to the physics of the model to generate any flow at and near the open boundaries This strategy avoids an overdetermined system Therefore one may resign on a swamp near the open boundaries as required in many other ocean and atmosphere models when the pressure and the flow fields are prescribed on and near the open boundaries Because 3 dimensional ocean data do not provide an information on the barotropic mode the computation of the full pressure is performed in the following order 1 transform a temperature and salinity profile into layer thicknesses with integrated tem perature and salinity as property 2 compute the baroclinic pressure i e the internal isopycnal interfaces and the according sea level elevation in order to obtain a zero pressure gradient at the bottom 3 take the barotropic mode from a global higher resolution ocean simulation and add its barotropic sea level gradient onto the sea level induced by the baroclinic modes The sea level gradient induced by the barotropic mode is assumed to be in geostrophic balance with the barotropic flow which is the data source 2 6 Computer Performance The model concept is to use isopycnal coordinates for the interior ocean The decision to use the flux form of the equations in order to be able to conserve mass which is the spatial integral of in situ density and to conserve momentum results in
67. y z t lt PPSFOUT gt model run means on xy section lt PPSFOUT gt model run means on xz section lt PPSFOUT gt model run means on yz section lt PPSFOUT gt section on xt plane lt PPSFOUT gt section on yt plane lt PPSFOUT gt section on zt plane lt PPSFOUT gt integral on xy section over t lt PPSFOUT gt integral on xz section over t lt PPSFOUT gt integral on yz section over t lt PPSFOUT gt point section along t lt PPSFOUT gt annual means on xy section lt PPSFOUT gt annual means on xz section lt PPSFOUT gt annual means on yz section lt PPSFOUT gt Table 5 9 Definition of Section Code Part II 5 5 4 1 10 File PPSF_INI contains all atmospheric quantities for forcing the ocean as well as the sea surface temperature and sea surface salinity Fields are written out for all months except the salinity which appears only once as an annual mean This file should be used to verify the correct initialization of the model forcing or to detect possible data problems Note that this file is created only when the atmospheric forcing is created and saved in FORCES Caution The file contains a large amount of data 5 5 4 2 The History Files 5 5 4 2 1 File PPSF_ALL is the history file which contains accumulated model quan tities that are written out with some frequency but in unsorted order Which quantities are written out by lt OCEPST1 gt lt OCEPST2 gt lt OCEPST3 gt lt OCEPST4 gt and lt OCEPST5 gt depends on
68. 02697 cos 38 1 46 Variations in the distance between sun and earth account for slightly more than 3 of the variations in the net global solar radiation Eqn 1 43 for is used in connection with the COADS cloudiness However because the statistics of the ECMWF reanalysis cloudiness differs significantly from that of the COADS data an alternative parameterization for x is needed if cloudiness is taken from the ECMWF reanalysis project The relation was found by tuning and is 1 0 0 0019 no0n 1 roy _ A 1 47 1 0 0 0019egua 1 0 2 8 1 0 0 0001re 2 n4 au Kb The net longwave radiation at the surface stays independent on the use of COADS or ECMWF reanalysis cloudiness 1 3 2 Transfer Coefficients The bulk coefficients were taken from Large and Pond 1981 1982 CHN q kyen E CLN O E 1 49 KJeMN Y CM 1 __ 1 50 CMN 1 AL y 2 K 1 51 K In cry 0 0346 1 53 q hn oe Can 1 54 a 1 54 ul cyu 1 55 T T 1 17x10 Tg 1 56 Here cy Cy and cr are the transfer coefficients for momentum sensible and latent heat respectively The subscript N denotes the transfer coefficient for neutral conditions For stable conditions it is used Z Ym Y F 1 57 2 Gp 0 2 5 x 10 F Ag 1 58 while for unstable conditions it is used vm 2n 1 X 2 In 1 X 2 2arctanX 7 2 1 59 da Y1 2n 1 X 2 1 60 X a 166 1 61 Z
69. 1 1 puh y 1 3 1 cos p y J 1 cos y yn Ays y Jg PL sa By san 2 15 2 1 1 5 Horizontal Diffusion of Momentum Aph 141 7 Aph sr Bv 141 3 Uvas 2Az 5141 AT y 1 Aph sr AP rr van Z Uv a 1 1 2A75 1 Aty r Aph sr AP 9 141 741 08 741 Uv 1 141 Uva Ays s41 Ayy 1 cos sy c08 5 741 CAph s 1 1 Aph sir41 7 eos si Ty r y Uy 1 J 1 2 17 Aysi Avy sy costo sy cos p s 5 1 If the older version for diffusing momentum instead of velocity is used then Aph is replaced by A and 7 by uph g A gga 2 16 2 1 1 6 Horizontal Advection of Momentum Depending on the value of the switch SW32 in JMOFLAG one of the following advection schemes is chosen 2 1 1 6 1 The Crowley Scheme The scheme by Crowley 1968 has the following form At gge puh y z uy uyay uy WUD Saga 2 18 E pub r 1 7 puh y r 7 2Az g 741 At Yuan uva n uy aa WD Azs puh y r z puh y r 1 7 2Az 1 At puhia wan tevin vas VEAD Rasa 2 19 cos P y 141 Pub yr 741 g cos p y y puh y r y 2Ays 741 cos e v s41 coste y 5 dy At 2A9Ys cos p y y puh y 1 5 cos y y 1 Puh yr 7 1 2Ays y cos e v s41 cos Y va vv g 1 vva a Crag tin 2 1 1 6 2 The Potential Vorticity and Energy Conserving Scheme The basic idea used to obtain a p
70. 1 The Dynamic Equations sir eee eee ee eae RAE ERS 1 5 2 The Thermodynamic Equations aoaaa aa 1 5 3 Coupling Sea Ice with Ocean Salinity aaa a 1 6 Atmospheric Energy Balance Model aaau aaa 2 000040 1 6 1 Transport Equations for Temperature and Humidity 1 6 2 Diagnostic Relations 1 ia a OS A AES 1 6 3 Radiation aa kaiera Ot A aia ot EA A 1 6 4 River Runoff Model Equations 2 0 2 2 aaa a 1 6 5 Glacier Model Equations a a 1 6 6 Soil Model Equations a pe e e A AA BS eS 1 6 7 Tuning of the TO Me ne cae 25 te a en a hh He A gm k ide Modelo ico ss oo a Wee ee oe et ere eB ee td ed LS racer Model ts AS oF dee ahs Ad Set ug ee Ae ae eS 2 Model Numerics 2 1 The Interior Ocean and the Mixed Layer Model 2 0 02 2 1 1 Discretization in Space on the Arakawa B Grid 1 13 15 15 16 19 20 20 20 23 23 23 25 26 26 26 27 28 28 29 30 32 32 33 33 34 30 30 30 30 36 38 2 1 1 1 Boundary Conditions srl rear se ee Ge aa ka aes 42 o 2 4 amp 2 A ee Bie ee T S 43 2 1 1 3 The Pressure Gradient 4044 a O Rae a 43 2 1 1 4 The Flux Divergence 2 data Bb Ad ee a A AA 44 2 1 1 5 Horizontal Diffusion of Momentum 44 2 1 1 6 Horizontal Advection of Momentum 44 2 1 1 6 1 The Crowley Scheme cocer 49s 284 44 2 1 1 6 2 The Potential Vorticity and Energy Conserving Scheme 45 2 1 1 7 Horizontal Diffu
71. 32 TEMP 1 J ITEMP K 100 273 16 SALT I J ISALT K 100 ENDDO ENDDO ENDDO ENDDO e File SSSALT contains the monthly mean sea surface salinity from Levitus on a 1 x 1 grid It is read by lt SEASALT gt with the following simplified statements DIMENSION ISSS 360 SSS 360 180 12 DO MONTH 1 12 DO J 1 180 READ 2014 ISSS I I 1 360 DO I 1 360 SSS 1 J MONTH ISSS 1 100 35 ENDDO ENDDO ENDDO e File TEMP MONTH contains the monthly mean sea surface temperature from Lev itus on a 1 x 1 grid It is read by lt SEATEMP gt from the local file SSTEMP with the following simplified statements DIMENSION ISSS 360 SSS 360 180 12 DO MONTH 1 12 DO J 1 180 READ 2014 ISSS I I 1 360 DO I 1 360 SSS 1 J MONTH ISSS 1 100 35 ENDDO ENDDO ENDDO e File SALT MONTH contains the monthly mean sea surface salinity from Levitus on a 1 x 1 grid It is read by lt SEASALT gt with the following simplified statements DIMENSION ISSS 360 SSS 360 180 12 DO MONTH 1 12 DO J 1 180 READ 2014 ISSS I I 1 360 DO I 1 360 SSS 1I J MONTH ISSS 1 100 35 ENDDO ENDDO ENDDO 5 4 3 1 2 COADS Data e File SSTGLOB contains a merged data set from the COADS and the Reynolds data set on a 2 x 2 grid The Reynolds data are used only if the COADS contain undefined values It is read by lt OCTEMP gt from the local file SSTEMP with the following simplified statements
72. 52 net de RA beh A ee ee 91 4 5 8 Routines for Vertical Interpolation o 91 4 5 9 Routines for Data Manipulation 91 4 5 10 Routines to Shift Data on edo a a a le ee Bie eee o 92 4 5 11 Routines to Bookkeep Budgets 2 0 0200 0 0 0 0 92 4 6 Passive Tracer Model ocetrae f our eee eh dae ae hae as 92 4 7 Open Boundaries oceopen f Ps Aa 93 4 7 1 Routines for Preprocessing Data aaa 93 4 7 2 Routines for Coordinate Conversion a oaoa a a 20000004 93 4 7 3 Routines for Pressure Calculation o 2 000004 94 4 7 4 Routines for Preparing Data puras kee ee Se a a de a 94 dio Ro tines for Listing ey tcs as ia ae A A A 96 4 8 Parallel Code oceutil f lo aaa 96 4 8 1 Pre and Postprocessing Lt di be hoe Bee A AAA 96 4 8 2 Data Transfer oi tl 2 aaoi a a ioa i a aiaa AIR a A th 96 4 8 2 1 Routines on the Server Side o o 96 4 8 2 2 Routines on the Client Side o o 97 4 8 3 Josef s Private Interface ccoo a 97 III User s Manual 99 5 How to Use PIPE 101 ac Wnstalling PIPE 4 3 ec lot rada ee eS AO BOS RE 101 5 1 1 Installation Step No 1 o 102 5 2 Installation Step Nos 2 ns a sen trat AK BS a 102 5 1 2 1 Installation of a RISC versi0dM o 102 5 1 2 2 Installation of a PVP version o 103 5 1 2 3 Installation of a MPP version o
73. 8 and NY 50 e T42 Grid NX 130 and NY 66 e T63 Grid NX 194 and NY 98 e T84 Grid NX 258 and NY 130 e T106 Grid NX 322 and NY 162 3 Use ISW45 2 in order to use the according land sea mask 4 Open cyclic boundary conditions must be enabled by SW15 1 5 It is recommended to set ISW1 1 ISW5 1 ISW8 1 ISW13 1 ISW18 0 ISW19 0 ISW41 0 6 Switch full physics on 7 Ensure SW4S 0 8 Ensure that Eulerian angles are switched off 9 Choose optionally an open North Pole with JSW28 1 or ISW28 2 Note that the variables for defining the model margins XGL XGR YGU and YGO are not used when one of the T21 T32 T42 T63 T84 or T106 grids is activated However a focus can be included But then the dimensions for NX or and NY must be increased according to IDENS or and JDENS 5 3 2 Self defined Grid The only difference from the pre defined grid is that now XGL XGR YGU and YGO can be chosen arbitrarily The switches SW15 for the cyclic boundary condition and JSW 28 must be chosen according to the application 5 3 3 Box Model A model that is driven by artificial data can be realized by setting ISW1 0 ISW5 0 ISW8 0 and SW41 1 Then the model expects the topography definition from lt DEFTOPO gt the atmospheric forcing from lt DEFORCE gt and uses an initial state defined by TEMMEAN and SALMEAN 5 4 Running PIPE In order to run the model a number of batch jobs and scripts are available All these scripts are written
74. 843 0 72523 ho 90 Q Elliptical Lunar 0 019256 0 64959 ho 389 90 Mf Fortnightly Lunar 0 041742 0 053234 259 Mm Monthly Lunar 0 022026 0 026392 s 9 po Ssa Semiannual Solar 0 019446 0 0038921 2hg Table 1 1 Constants of Major Tidal Modes As an example of the quality of PIPE s tide submodel Fig 1 3 shows the simulated am plitudes of the M2 tide and Fig 1 4 the M2 phases The model used the 5 minutes NOAA topography without any further tuning The tide submodel used a time step of 10 minutes while the 3 dimensional ocean model used a time step of 6 hours The residual current required by the bottom stress formulation of the 3 dimensional ocean model and by the transport equation for acyive and passive tracers is computed by time averaging the tidal flow through pes At fa Uriaedt 1 8 Tracer Model In order to be able to in addition simulate passive tracers PIPE has a separate module that can be run in an online or offline mode This means that the tracer model can be run synchroneously with the ocean model by using those model variables that are used in parallel for the transport of temperature and salinity In the offline mode however it is assumed that the forcing data are saved on a file This is created while running the ocean model once Then these data are used to run the tracer model while the rest of the ocean model is disabled Thus there is a considerable saving in CPU time Because the ocean mod
75. ALL CLOCK LTIME CALL TIMECON KTOTAL MSEC MINUTE MHOUR MDAY MONTH MYEAR WRITE 40 KODEQ KODES N1 N2 N3 KODEL KTOTAL MSEC MINUTE MHOUR MDAY MONTH MYEAR LDATE LTIME EXPERIM WRITE 40 CHARCO CHAR 1 CHARCIFLG I J I 1 N1 J 1 N2 WRITE 40 X1 D I 1 MNM1 WRITE 40 X2 J J 1 N2 WRITE 40 FIELD I J K I 1 N1 J 1 N2 K 1 N3 Code Number Variable Unit Reference velocity stream function m s absolute stress m s entrainment detrainment rate ms VERTIC cross isopycnal mixing ms VERTUP VERTDN downward bubble salt flux ms SST error K maximal depth of surface convection m lt QCONV gt Table 5 3 Definition of Quantity Code Diagnostic Variables Code Number Variable Unit Reference observed sea level pressure Nm lt OCESTEP gt topography m DEPTH observed sea surface temperature SST observed air temperature AIRTEMP observed salinity SSS observed wind stress TAUX TAUY observed wind speed UVABSOL observed standard deviation UVDEV observed humidity HUMID observed cloudiness CLOUDS observed precipitation RAIN Table 5 4 Definition of Quantity Code Forcing Variables where the variables have the following meaning 1 10 KODEQ is the first code number and defines which quantity is written onto e g PPSF ALL See also table 5 1 through table 5 6 KODES describes the type of section on which the quantity is written out See also table 5 7 and 5 8 N1 is the first
76. AQFLX gt exports the model computed heat fluxes lt AOTAU gt writes the wind stress components into the common data area Entry lt OATAU gt exports the model computed wind stress lt AOUSTAR gt writes the turbulent kinetic energy u3 into the common data area Entry lt OAUSTAR gt exports the model computed turbulent kinetic energy lt AOTICE gt writes the snow surface and snow ice interface temperature into the com mon data area Entry lt OATICE gt exports the model computed snow surface and snow ice interface temperature lt AOAIRT gt writes the surface air temperature into the common data area Entry lt OAAIRT gt exports the model computed surface air temperature lt AOOCET gt writes the ocean surface temperature into the common data area Entry lt OAOCET gt exports the model computed ocean surface temperature lt AOOCES gt writes the ocean surface salinity into the common data area Entry lt OAOCES gt exports the model computed ocean surface salinity lt AOAHUM gt writes the surface air humidity into the common data area Entry lt OAHUM gt exports the model computed surface air humidity lt AOACLD gt writes the atmospheric cloudiness into the common data area Entry lt OACLD gt exports the model computed atmospheric cloudiness lt AOAUVA gt writes the atmospheric scalar surface wind into the common data area Entry lt OAAUVA gt exports the model computed atmospheric scalar surface wind lt AOAU
77. DDO ENDDO e File ECSLP contains the surface pressure from the ECMWF reanalysis Project on a T106 grid It is read by lt ECSLP gt from the local file MSLP with the following simplified statements DIMENSION ISLP 320 SLP 320 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 CISLP I 1 1 320 DO I 1 320 SLP I J MONTH ISLP 1 95000 ENDDO ENDDO ENDDO e File ECTAU contains the surface wind stress from the ECMWF reanalysis Project on a T106 grid It is read by lt ECTAU gt from the local file UVGLOB with the following simplified statements DIMENSION ITAUX 360 ITAUY 360 TAUX 360 160 12 TAUY 3600 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 CITAUX I 1 1 360 READ 2014 ITAUY I I 1 360 DO I 1 144 TAUX 1 J MONTH ITAUX I 4000 5000 TAUY 1I J MONTH CITAUY I 4000 5000 ENDDO ENDDO ENDDO e File ECUSC contains the surface scalar wind from the ECMWF reanalysis Project on a T106 grid It is read by lt ECUSC gt from the local file UVABS with the following simplified statements DIMENSION IUSC 360 USC 360 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 1USC I I 1 360 DO I 1 360 USC I J MONTH IUSC 1 500 ENDDO ENDDO ENDDO e File ECSTD contains the surface scalar wind from the ECMWF reanalysis Project on a T106 grid It is read by lt ECSTD gt from the local file UVDEV with the following simplified statements DIMENSION ISTD 360 STD 360 160
78. EP gt lt OCEPOST gt and lt OCESTOP gt The computational work is pushed into lower level routines Thus the driving routines widely consist only of calls to these lower level routines Nevertheless lt OCEINIT gt and lt OCESTEP gt are pretty long 57 code data main program initialization global data sets OCEINIT TRAINIT MY_OGCM files see text parameters time stepping model forcing Block Data OCESTEP FORCES TRASTEP OBNDSEC restart files postprocessing OCEPOST OCDAT TRAPOST TRDAT final VO history files OCESTOP TRASTOP PPSF_ALL PPSF_HIS quick look files see text protocol OUTLIST Figure 3 1 Main structure of main program lt MY OGCM gt which describes the interaction between the code and the data base OCEINIT read restart files s READM initialize ocean state initialize topography LAYOUT INTLAY initialize atmosphere initialize atmospheric forcing interpolate to new grid grid interpolation OLDINIT o MY_OGCM update atmospheric forcing atmospheric forcing Figure 3 2 Main structure of lt OCEINIT gt which is the driving routine for initialization sine SpE SURFLUX compute explicit part compute tides compute tides POTTIDE gt solve for
79. ETA gt and performs the reverse transformation via lt ZETATOH gt However for the output onto the postprocessing files the vertical coordinate starts at the mean sea level and is positive downward see various definitions of default levels by COORDIN in ocepost f 4 5 1 Postprocessing e lt OCEPST1 gt writes time averaged data onto the history file This routine is called at every time step Entry lt OCEPST2 gt writes data onto the history file It is called at each 15th of a month Entry lt OCEPST3 gt writes data onto the history file at each n th time step where n is defined by SW17 Entry lt OCEPST5 gt writes data onto the history file for SW17 0 This routine is used to create e g movies by writing data at each time step 4 5 2 General Output Routine The output of data is not straightforward in a model that uses Lagrangian coordinates Since the data postprocessing is made easier if the quantities are available on levels instead as layers a routine is offered for each of the standard problems that transforms quantities from Lagrangian coordinates to levels of constant depth The available routines are listed below e lt PPSFOUT gt These routines write data in well defined format see 5 5 2 onto the history file The data format henceforth is called PPSF Post Processing System Format 4 5 3 Output of Averaged Quantities For many purposes it is sufficient to provide data on a time averaged basis e g when quant
80. How to Apply the Coupling Interface The following example shows how to introduce data from outside the ocean model without making any changes to the bulk of the code The listed changes apply only for lt MY_OGCM gt Assume that data with daily resolution are prepared on FORCING with the sequence as defined in the example below and that the data are already interpolated to the model grid Furthermore assume that the model run should cover one month of prediction and uses a time step of a day In this case set MAXREP 2 and DTMAX 43200 in order to sychronize the data flow from FORCING with the model time stepping After having set NINIT 1 the parameters NHEAT NFRESH NSTRESS and NMIX must be chosen dependent on whether the forcing is stored as absolute value or as anomaly on FORCING Note that the ice coupling can only be switched on if the snow and ice surface temperature are consistently computed from the model snow and ice cover This can only be ensured in a coupled AGCM OGCM where the AGCM uses the same algorithm to compute the snow and ice temperature as lt SURFLUX gt The example is as follows PARAMETER NX 130 NY 63 DIMENSION USTERNC NX NY TAUX NX NY TAUY NX NY DIMENSION SOLAR NX NY FRESH NX NY HEAT NX NY CALL OCEINIT MSEC MINUTE MHOUR MDAY MONTH MYEAR OPEN 90 FILE FORCING FORM UNFORMATTED DO LREP 1 30 READ 90 USTERN I J I 1 NX J 1 NY READ 90 TAUX 1 J 1 1 NX J 1 NY READ 90 TAUY 1 J
81. I 1 NX J 1 NY READ 90 SOLAR I J 1 1 NX J 1 NY READ 90 FRESH I J 1 1 NX J 1 NY READ 90 HEAT 1 J 1 1 NX J 1 NY CALL AOUSTARCUSTERN CALL AOTAU CTAUX TAUY CALL AOQFLX HEAT SOLAR CALL AOPME FRESH CALL OCESTEP MSEC MINUTE MHOUR MDAY MONTH MYEAR ENDDO CLOSE 90 CALL OCEPOST CALL OCESTOP In this style all routines can be called which start with AO stands for atmosphere to ocean for importing data into the ocean model as well as all those routines which start with OA stands for ocean to atmosphere for exporting data out of the ocean model The latter must be called after lt OCESTEP gt 5 4 7 How to Use the Bias Correction Bias correction is a technique that has been developed for recent coupled experiments with ECHAM4 and OPYC3 In opposite to flux correction based on monthly means bias correction means that an annual mean correction for heat fresh water and stress is applied In order to use this technique the model has to pass three steps These are for 1 spinning up the model with full feedbacks for temperature and salinity During the first period the switches NHEAT NFRESH and NSTRS must be set to unity as well as KHEAT 1 KFRESH 1 and KSTRESS 1 This means that the model expects fluxes from e g data or an atmospheric model that are fed through the coupling interface If ISW48 1 then the bias correction is computed as temporal average however because model physics is nonlinear the result
82. It also creates the EBM s orography lt READTOP gt reads the file TOP1DEG lt LAYOUT2 gt initializes the x and y coordinates reads the topography from TOP5MIN interpolates it onto the model by using an envelope algorithm to take care of the variance of the bottom topography and arranges the land sea mask It also creates the EBM s orography lt READTOF gt reads the file TOP5MIN lt DEFTOPO gt allows the user to construct his own topography and land sea distribu tion It is called by choosing SW8 0 lt REMEDY gt fills up isolated depressions in the orography which is essential to avoid lakes in the soil model of the EBM lt INTLAY gt is the underlying routine to interpolate the observed ocean temperature and salinity onto the model grid and to generate the layer thicknesses and temperature and salinity for each layer lt INTERPO gt reads the Levitus ocean temperature and salinity from SALTEMP and interpolates it onto the model grid lt LAYDEPT gt optimizes the layer thicknesses and their respective temperature and salinity in order to yield the prescribed potential density within each layer lt INTEGRA gt computes the vertically averaged potential density between two given depth levels from the Levitus data set lt CREATIV gt initializes the vertical coordinates that are stored in POTSOLL and also initializes the whole density field from the initial temperature and salinity profile 4 4 5 Ge
83. J W 1972 A study of some aspects of the radiation and the heat budgets of the southern hemisphere oceans Meteorol Stud 26 562pp Bur of Meteorol Dept of the Interior Canberra Australia 7 32 Key Publications with PIPE sorted in time Oberhuber J M 1988 An atlas based on the COADS data set The budgets of heat buoyancy and turbulent kinetic energy at the surface of the gloabl ocean Max Planck Institut fur Meteorologie Hamburg Report No 15 199pp Miller A J J M Oberhuber N E Graham and T P Barnett 1992 Tropical Pacific Ocean Response to observed Winds in a Layered General Circulation Model J Geophys Res Vol 97 No C5 7117 7340 Holland D M L A Mysak D K Manak and J M Oberhuber 1993 Sensitivity Study of a Dynamic Thermodynamic Sea Ice Model J Geophys Res Vol 98 No C2 2561 2586 Oberhuber J M 1993 The OPYC Ocean General Circulation Model Technical Report No 7 Deutsches Klimarechenzentrum GmbH Hamburg Oberhuber J M 1993a Simulation of the Atlantic Circulation with a Coupled Sea Ice Mixed Layer Isopycnal General Circulation Model Part I Model Description J Phys Oceanogr Vol 23 No 5 808 829 Oberhuber J M 1993b Simulation of the Atlantic Circulation with a Coupled Sea Ice Mixed Layer Isopycnal General Circulation Model Part II Model Experiment J Phys Oceanogr Vol 23 No 5 830 845 Miller A J D R Cayan T P Barnett N E Graham and
84. JSW11 0 1 potential vorticity lt LISTING gt JSW12 0 1 velocity lt LISTING gt JSW13 0 1 potential density lt LISTING gt JSW14 0 1 SST error lt LISTING gt JSW15 0 1 sea ice and snow lt LISTING gt JSW16 0 1 topography lt LISTING gt JSW17 0 1 index field lt LISTING gt SP Orr Or OOGrrF O00BeRRI KF KF D Table 4 17 Switches for Listing Default 0 Variable ISW1 ISW5 ISW8 ISW13 ISW15 ISW18 ISW19 ISW28 ISW41 ISW48 Value Meaning 0 1 0 1 0 1 0 1 0 1 switch from artificial to observed forcing switch from artificial to observed initial state switch from artificial to observed topography switch from linear to parameterized heat fluxes switch from closed cyclic or open boundaries switch from linear to parameterized salt fluxes switch from diurnal cycle switch for closed open North Pole switch from artificial to observed atmosphere 0 0 0 0 0 1 0 0 1 0 0 1 2 0 0 1 0 0 1 2 3 0 Table 4 18 Switches for Type of Forcing Variable Value Default Meaning switch for heat fluxes switch for layer diffusion switch for mixed layer physics switch for river runoff when EBM absent switch for Coriolis force switch for spherical coordinates switch for diurnal cycle switch for updating isopycnal coordinates switch for convection switch for forcing through sea level pressure switch for tracer model switch for EBM switch for salt calculation switch for tides switch for Greg Holloway s Ne
85. Model Run 5 Adjust Source PVP Model model cri to perform a model initialization 6 Run Source PVP Model model cri 7 Adjust Source PVP Model model cri to perform a continuation run 8 Run Source PVP Model model cri In order to setup an arbitrary version for a PVP computer the file Source parampvp h is inserted into the code through an include statement The dimensions nz_tot and nz for the x direction and the dimensions ny_tot and ny for the y direction each have identical meaning These dimensions are set through the script Source PVP Model setup There is no need to modify parampvp h during the installation 5 1 2 3 Installation of a MPP version 1 Adjust Source MPP Model setup in order to obtain FORTRAN codes that use the proper dimensions and other options as documented in setup 2 Run Source MPP Model setup 3 Adjust code essentially ocemain f to obtain a model domain and proper choice of the model parameters 4 Run Source MPP Model setup again The files for compilation now are located on Source MPP Model code Transfer these files if the setup was not executed on the main frame of the MPP 5 Run Source MPP Model make_all job in order to create libraries which by default are copied to the subdirectory Source RISC Model lib and to create an executable copied to the subdirectory Source RISC Model bin 6 Adjust Source RISC Model model mpp to
86. ORIMIX 2 64 ADD 20 Block Data IGE piati 2 6249 4 lee a ae de 65 4 2 2 2 4 Block Data PHYSICS 5344 aaa ae etd 65 4 2 2 3 Block Data which concern Initialization 65 4 2 2 3 1 Block Data LAYOUT 2501000304 BR Rae RAY 65 4 2 2 3 2 Block Data ROTATE 66 4 3 4 4 4 5 4 2 2 4 Block Data which concern Data Postprocessing 66 4 2 2 4 1 Block Data POSTPRO 66 4 2 2 4 2 Block Data PRINTER 67 4 2 2 4 3 Block Data QBDGT 67 4 2 2 4 4 Block Data DEFINE 67 4 2 2 5 Block Data which concern Model Layout 67 4 2 2 5 1 Block Data SWITCH 67 4 2 2 6 Block Data which concern Numerics 68 4 2 2 6 1 Block Data TUNING 68 The Ocean Model ocestep f 4 rr 2 4 E AA A AA A 74 4 3 1 Interior Ocean Model Subroutines o 74 4 3 2 Mixed Layer Model Subroutines aooaa aa a a 000004 76 4 3 3 Active Tracer Subroutines f o 2 a foe k Roe YARD EY SAR YS 76 dal The SCA leo Model 0203 aii e a peti Se EE AES E TT 4 3 5 Atmospheric Energy Balance Model aoaaa aaa a TT A30 Wide Model 25 49 50 ce Dhar a i of Sak BaF AE BOG BLES 77 4 3 7 The Equation of State us cs ise eee eee ee es TT 4 3 8 Model Morcing ic a a pee Oe AMER ANE A 78 Preprocessing oceproct f2 2 4 4 Ana AS AA ESR ee 78 AAL Driving Routines y AS ink O
87. Ox 4 Azs 1 1 phu ly 1 1 1 phu v r 3 UY 1 3 UV 11 3 4 Az 57 ph u _ phu yir z413 Phu yir n coste vr Uv 1 341 Cos P y Uy 1 3 2 23 2 Oy 2 cos p y y costo 141 AYs s 1 phu yira phu y r s 1 cos P y J Uy 1 1 cos P y J 1 UV 1 31 2 cos p y 7 1 coste 1 Ays The local vorticity is discretized by Gig US 1 J 1 AL SR Utd 2 24 Tv I _ Us 1 1 3 US 1 J US 1 1 3 1 US 1 J 1 2AYy 7 This involves a horizontal filter that removes efficiently B grid computational modes from the local vorticity 2 1 1 7 Horizontal Diffusion of Scalar Variables 2 Fg Len Asai san Osun 2 25 Ox Ox 2Ary nAz 1 Asg ir1 3 Asa n Osa Os 1 1 3 2Ary 1 Ars I DED ie E Asa J41 Asian Os r 1 1 Os 1 3 esley 2 26 Oy Oy Aysi Ayy y coslp y J 1 coslp y 5 Asin Agir s 1 OSa res eer Aysin Ayy s 1 costp y J 1 cos p y y The same discretization is used for the salinity S 2 1 1 8 Horizontal Advection of Scalar Variables The Crowley Scheme The advection scheme is similar to that discussed by Smolarkiewicz 1982 It is an implicit version of the scheme presented by Crowley 1968 The scheme is less diffusive than the upstream scheme and reduces the tendency of centered difference schemes to overshoot It was modified so that it converges towards the upstream scheme for CFL gt 1 This has to be done as otherwise the artificial diffus
88. RISC version worked well on a DOS machine if compiled and loaded by hand 5 1 2 Installation Step No 2 In each of the directories Source PVP Source RISC and Source MPP there exists a subdirectory Model where a further script setup is located These scripts install by default a T42 model with 130x 66 horizontal grid points and 11 layers In the directories for the PVP RISC and MPP version files ocemain f ocestep f oceproc f ocepost f ocemods f and ocetrac f appear These files need to be changed according the rules as described in the next section in order to setup a different model version than the T42 For the MPP version a further step is to split these FORTRAN files into files containing one subroutine each only These are copied to Source MPP Model code and are used by the makefiles to create libraries and finally executables for a CRAY T3D The compile jobs can be found in Source PVP Jobs Source RISC Jobs and Source MPP Jobs During a first installation of the T42 model it is recommended to leave the script Source RISC Model setup unchanged also leave unchanged the subsequently installed file Source RISC Model ocemain f and test whether you succeed in running the T42 model If this step is successfull start to setup your own model application 5 1 2 1 Installation of a RISC version 1 Adjust Source RISC Model setup in order to obtain FORTRAN codes that use the proper di
89. ROP gt for adjustments that are required for the open boundaries 5 4 11 Conceptual Problems PIPE has been developed to simulate the global ocean circulation with isopycnal coordinates Therefore a line relaxation scheme has been chosen in order to deal with converging grids without the need for any spatial strong filters The disadvantage is that the matrices for tridi agonal or block tridiagonal linear equation systems have to be setup and solved in an iterative process at each time step This computationally demanding concept however combines the excellent wave propagation properties of a layer model see Oberhuber et al 1998 e g needed for regional applications with an uncritical treatment of polar problems in a global model see Roeckner et al 1996 Fancy grid rotations are not needed in PIPE However due to the use of an implicit time stepping scheme one has to be aware of specific properties that are important to know An analytic solution of a linear model solved with an implicit scheme has the following properties e The wave speed always is limited by the ratio between wave length and time step Thus long waves might travel nearly at their theoretical speed while short but theoretically fast waves are slowed down significantly In consequence slow modes like baroclinic ones unlikely are slowed down while fast external modes i e the barotropic modes are slowed down significantly if they are short Long barotropic waves are le
90. S BWSR ASG 78 4 4 2 Grid Initialization 3 Bae te ee A AA 78 4 4 3 Initialization of Atmospheric Forcing o 80 4 4 4 Initialization of Ocean State 2 4 44 65h a bee ee ee ee es 81 4 4 5 Generation of Isopycnal Coordinates o 82 4 4 6 Conservation Properties ue xao ea a oak ee 82 4 4 7 Boundary Conditions 5 24 a e AAA 82 4 4 8 Routines for Listing ade AAA AA 83 4 4 9 Routines for Model I O ahaaa fatty ae Fae a a AA 83 4 4 10 Routines for Model and CPU Time aaa a 84 4 4 11 Routines for Vertical Flag Calculations a oaoa aaa aaa 84 4 4 12 Routines for Data Manipulation o 84 4 4 13 Routines for Filtering 0 pia bw we Se ee Bie Dee ee 85 4 4 14 Routines for Filling Up eaaa aa aaa e 85 4 4 15 Routines for Error Detection 44202000 4 a aa a 85 Postprocessing ocepost f oaoa aa A A 85 Lol JPOSEPFOCESSINe 204 e ir a a a ele Bid Dae es 86 4 5 2 General Output Routine evocar eae a a we 86 4 5 3 Output of Averaged Quantities 2 a e 86 4 5 4 Output of Instantaneous Fields 0 0 87 4 5 5 Tridiagonal Solvers o 03 amp es dat Seis a a a a 87 0 0 Coupling Interface ca 5 0 2 Samoa res ea bea ea da 88 4 5 6 1 Data Transfer into the Ocean Model 88 4 5 6 2 Data Transfer out of the Ocean Model 90 4 5 6 3 Internal Data Transfer 0 0 020202 020008 4 90 4 5 7 Coupling with Tracers 2
91. SNOW 0 3 Wk m heat conductivity for snow lt SURFLUX gt EDDYS 2000 m s scalar diffusion coefficient lt SEAICE gt EDDYUV 2000 m s vector diffusion coefficient lt SEAICE gt RHOICE 900 kgm ice density lt SURFLUX gt RHOSNOW 300 kgm snow density lt SURFLUX gt COMPMAX 0 99 maximal ice concentration lt SEAICE gt Table 4 6 Snow and Sea Ice Parameters Variable Value Unit Meaning Reference KTOTAL any 8 initial model time lt OCEINIT gt KTSHIFT 1296000 s length of half a month lt OCEINIT gt SALMEAN 35 0 g kg mean salinity lt DEFORCE gt TEMMEAN 298 15 K mean temperature lt DEFORCE gt HTOTAL 6000 m reference depth of ocean lt OCEINIT gt UMFANG 40030174 m earth s circumference lt MAKEXY gt Table 4 7 Parameters for Coordinate Generation 4 2 2 4 2 Block Data PRINTER defines switches that control which quantities are writ ten to OUTLIST A focus allows to print every grid point within the some area see table 4 17 4 2 2 4 3 Block Data QBDGT defines the vertical layer index KBDGT used to bookkeep the budgets of various terms and the frequency at which the budgets are written out The latter is done each DTBDGT seconds model time 4 2 2 4 4 Block Data DEFINE defines all coefficients required for the equation of state UNESCO 1981 and for the equation of Bryden 1973 that converts temperature into potential temperature see Appendix C 4 2 2 5 Block Data which concern Model Layout 4 2 2 5
92. ST gt 1 generate file PPSF_BAS lt OCEPOST gt 0 1 1 generate file PPSF_SEC lt OCEPOST gt 0 1 1 generate file PPSF_VER lt OCEPOST gt 0 1 1 generate file PPSF_FLX lt OCEPOST gt 1 1 1 0 ISWPOST ISWPOST ISWPOST 8 0 1 0 1 generate file PPSF_FOR lt OCEPOST gt generate file PPSF_INI lt OCEPOST gt ISWPOST 0 1 generate file PPSF_TID lt OCEPOST gt ISEC 10 1 lt ISEC lt NX x indices for cross sections lt OCEPOST gt JSEC 10 1 lt JSEC lt NY 0 y indices for cross sections lt OCEPOST gt VSEC 10 any in m depth of horizontal sections lt OCEPOST gt KSEC 1 lt KSEC lt NZ 0 z index for layer sections lt OCEPOST gt AVERTIM any gt 1 2592000 averaging period lt OCEPOST gt 4 5 6 7 8 9 Table 4 16 Parameters for Quick Look System Variable Value Default Meaning Reference KSW 0 1 focus used for listing lt DRUCK gt MAXIM 1 lt MAXIM lt NZ number of printed layer lt LISTING gt ILEFT any lt IRIGHT left index used for focus lt LISTING gt IRIGHT any gt ILEFT right index used for focus lt LISTING gt JLOW any lt JUPP lower index used for focus lt LISTING gt JUPP any gt JLOW upper index used for focus lt LISTING gt JSW1 0 1 horizontal fluxes lt LISTING gt JSW2 0 1 surface and interfaces lt LISTING gt JSW3 0 1 layer thickness lt LISTING gt JSW4 0 1 temperature lt LISTING gt JSW5 0 1 salinity lt LISTING gt JSW9 0 1 diffusion coefficient lt LISTING gt
93. Sea Ice with Ocean Salinity During freezing salt is ejected from the ice but is not confined to the ML alone It is assumed that a fraction of this ejected salt penetrates into the subsurface ocean with a length scale that depends on the stratification in the surface layers Therefore a higher stability should give a shorter length scale The following ad hoc parameterization for salt transfer out of the mixed layer into all deeper layers k 2 N is chosen as O61 k gt 091 aS p exp 2 WZS i 1 100 P hy Ri max Fy 0 p1 S1 S exp where h 20kgm is a free parameter tuned to obtain a reasonable model response The salinity budget in the ML is the result of the salinity gain due to freezing ice and the salinity loss due to the downward transfer of a fraction of the salinity gain The salinity change in the ML is given by Fr En erp 2492 pp 1 101 Ri pil Si S Fy h P Thus water formed from melting ice F lt 0 is completely mixed within the ML since Ri 0 for k 2 N but freezing ice injects a fraction of the salt into deeper layers This allows the model to build up a salinity stratification in the Arctic basin although in the annual mean there is a net transport of salt from the sea ice through the ML into the deeper ocean The balance between downward transport and upward diffusion yields an equilibrium state for the salinity stratification in areas covered by sea ice 1 6 A
94. Technical Report No 19 Description of the Parallel Isopycnal Primitive Equation OGCM PIPE with models for Sea Ice and Snow Mixed Layer and Tides coupled to an wet atmospheric energy balance model with river runoff soil and glaciers Josef Maximilian Oberhuber Deutsches Klimarechenzentrum GmbH Bundestrafe 55 D 20146 Hamburg Edited by Deutsches Klimarechenzentrum GmbH Modellbetreuungsgruppe Hamburg April 1999 Revision No 1 ISSN 0940 9327 Contents I Model Description 1 Model Physics and Dynamics 1 1 The Interior Ocean Model ida a a 1 1 1 Ocean Model Equations 2 3 aaa Ee eee ee 1 1 2 Equation of State A as a a ia diria a y 1 2 The Mixed Layer Model sos erer mais da a a a A A 1 2 1 Physical Background cairo A ee a RES 1 2 2 The Mixed Layer Model Equations 00 0 1 2 3 Coupling of the Mixed Layer to the Interior Ocean 1 3 Parameterization of the Surface Fluxes 2 2 ee 1 3 1 Parameterization of the Surface Heat Fluxes 1 3 2 Transter Coefficients yla pta Soest See bd ee Deas ae A a 13 3 Evaluation of the Net Fresh Water Flux 1 3 4 Estimate of Turbulent Kinetic Energy Input 1 4 Parameterizations of Internal Diffusion o o 1 4 1 Vertical Mixing Coordinate Maintenance o o o LADA ACONVSCUION 2 9m cose ayes da A IA A ASA ATA EA 1 5 The Snow Sea ce Model gt La ERE ADA AA A 1 5
95. V its standard deviation This equation can be compared with Kolmogorov s idea to relate a diffusion coefficient with eddy length scale and turbulent energy Here we may identify the mean speed V with the eddy length scale and the wind s standard deviation o V with the turbulent energy The EBM s internal forcing is the result of latent heat release due to condensation which is assumed to fall out immediately as rain Like in reality the heat release is compensated by cooling due to upward motion which is an unresolved dynamical effect of heating The internal forcing is parameterized as Qaim 4Ly P 5V 7 1 106 In order to parameterize the precipitation P the idea is to consider it as sum of convective and large scale precipitation While the first one is typical for low latitudes where near surface convergence of moisture is driving cloud formation the second one is asssumed to be created by large scale variability due to low pressure systems Therefore it is the concept to parameterize rain with the help of wind convergence for low latitudes and standard deviation of the wind for higher latitudes The relation is P Veum ego V erminfo V maz 0 r 0 6 1 107 Assuming a typical thickness of the humid boundary layer H the relative humidity r is com puted from wV E 1 108 pad where q is the specific humidity see Eqn 1 38 The evaporation Fis computed from the latent heat flux see Eqn 1 36 An exam
96. VS gt writes the atmospheric standard deviation of the scalar surface wind into the common data area Entry lt OAAUVS gt exports the model computed atmospheric standard deviation of the scalar wind lt AOSURP gt writes the atmospheric surface pressure into the common data area Entry lt OASURP gt exports the model computed surface pressure lt SFOR2MOD gt interpolates scalar data from an arbitrary grid onto the model grid It is required if data defined on another grid have to be prepared before one of the routines lt AOPME gt to lt AOTICE gt can be called lt SFOR2MOD gt also is able to fill gaps in the data source using a Laplacian filter lt VFOR2MOD gt interpolates vector data from an arbitrary grid onto the model grid It is required if data defined on another grid have to be prepared before one of the routines lt AOPME gt to lt AOICE gt can be called lt VFOR2MOD gt also is able to fill gaps in the data source using a Laplacian filter 4 5 6 2 Data Transfer out of the Ocean Model The output data of the following routines are defined on the ocean model grid This routines typically are required to supply an atmospheric model with data that are required for either grid interpolation or flux calculations e lt OAGRID gt provides the land sea mask and the x and y coordinates of the ocean model to be used e g for interpolation between two different model or data grids e lt OASST gt reads the model sea surface
97. able 4 14 Definition of Eliminating Points Variable Value Default Reference lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt lt GRITUNE gt Meaning o switch for Panama bridge switch for Florida Peninsula switch for Thule land bridge or channels switch for Hudson Bay closed or channel switch for cut off of Ochotskic Sea switch for bridge Spitzbergen Norway switch for bridge Alaska Russia or chan switch for bridge Gibraltar or channel switch for Weddell Sea Peninsula switch for Cuba Thaiti or wipeout switch for bridge Malaysia and Indonesia switch for channel Met Sea and Black Sea switch for bridge Australia Indonesia switch for bridge Greenland Canada switch for island Novaja Semlja switch for island Severnaja Semlja switch for island Japan switch for island New Zealand switch for bridge Baltic North Sea switch for bridge Red Sea and Indic 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Table 4 15 Grid Tuning Switches Variable Value Default Meaning Reference ISWPOST 1 ISWPOST 2 ISWPOST 3 ISWPOST ISWPOST 0 1 0 1 0 1 1 generate file PPSF_SUR lt OCEPOST gt 1 generate file PPSF_ICE lt OCEPO
98. al 2 at 2 dt cet T Q T 0 1 88 SAP APSO LB Here T is the mixed layer temperature k the conductivity of snow k the conductivity of ice and Q T s the heat flux through the snow or at the ice surface c the specific heat capacity of ice and p is the density of snow and h q is the effective thickness of an ice floe and p s p q is the effective snow depth Note that s is a cell averaged snow depth at the ice density The equation for the heat accumulation 1 87 contains two time derivatives The first involves the ice surface temperature T assuming that the time derivative of the sea surface temperature Ti is negligible The second time derivative term measures the vertically averaged heat content change of the snow layer where T and 7 are used as surface and bottom temperature of the snow layer In order to treat the special cases of no snow and no ice the equations are rewritten to give Ta T for no snow and to give Ta T T for no snow and ice For that 1 88 is used to modify 1 87 and the resulting equations are multiplied by the ice or snow mass This finally yields pihi Dn pisi 3 n A k 4 8 O higu Da Tr hiQ Ls hie T T 0 189 o 1 90 4 In the case that T exceeds the melting temperature Tm 273 16 of ice and snow T Tm is used in 1 87 to compute the ice surface temperature Tp The cell averaged skin temperature Ta is defined by
99. al file is STMONO9 File SALTEMP OCT is the 3 dimensional observed temperature and salinity for Oc tober Format is identical with SALTEMP The local file is STMON10 File SALTEMP NOV is the 3 dimensional observed temperature and salinity for November Format is identical with SALTEMP The local file is STMON11 File SALTEMP DEC is the 3 dimensional observed temperature and salinity for De cember Format is identical with SALTEMP The local file is STMON12 File TIDES ASS contains the observed tide data The file is read by lt INTIDES gt The local file is TIDES File BAROFLOW contains the barotropic transport computed with a T106 version of the PIPE model The file is read by lt SELECOBND gt and lt STREAM_T106 gt The local file is BAROFLOW 5 4 3 2 Model Forcing The model forcing is derived from the global data set by bilinear interpolation onto the model grid and is saved in FORCES The global data sets are required only once A continuation run needs only a restart file OCDAT and a forcing file FORCES If other data sets should be used there are two choices 1 If the data are monthly means then the responsible routine for interpolating that quantity onto the model grid can be adjusted for the other data set 2 If the data are not monthly means then the simplest way is to use the coupling interface routines If open boundary conditions are used the file OBNDSEC is needed If the trac
100. al type with NY x NZ number of equations Each linear equation system solves the matrix on a yz section NV NX 1 2 is the result of the used line relaxation method 4 5 6 Coupling Interface The following routines allow easy use of the model either in coupled mode or with simple forcing from atmospheric data that are not provided by the model s data bank The idea is to define a common data area for all fluxes required to force the model and for those quantities another model needs to analyse the ocean state By calling these routines the user does not need to know where the data are stored Thus the model can be used as black box The further advantage is that the code ensures that the forcing is correctly read from this common data area whenever needed Note also that the fluxes can be interpreted either as absolute values or as anomalies In the latter case the climatology is taken from the model 4 5 6 1 Data Transfer into the Ocean Model The data passed into the subsequent routines must already be defined on the ocean model grid and must have the same dimensions as the model grid Each routine also contains an entry that performs the opposite operation 1 e to export data from the ocean model common blocks lt AOPME gt writes the fresh water flux into the common data area Entry lt OAPME gt exports the model computed fresh water flux lt AOQFLX gt writes the heat flux and the solar radiation into the common data area Entry lt O
101. amics is then given by the heat flux Q and the heat flux induced by the entrainment rate w pa e Sen Cp m PiCpm where Cpm is the latent heat of fusion 1 95 The thermodynamic part of the ice concentration equation differs slightly from Hibler s formulation The change of the ice concentration F is evaluated through mes Fyqi 2h if Fp lt 0 ee where hy is interpreted as ice thickness that immediately builds up within leads during freezing conditions The aging process of snow is expressed as a rate F at which snow is converted into ice The first term parameterizes the metamorphosis of snow cristalls to ice by assuming a simple snow depth and time scale dependent y The second term adjusts the snow depth when snow suppresses the ice surface below the sea level Thus this mechanism represents an upper threshold for the snow depth at a given ice thickness The change in snow thickness F is 1 8 hs g h Mi At p Pi where y reflects mean values for the conversion rate The forcing of the snow depth equation F depends on precipitation minus evaporation melting and loss of snow mass due to reduction of the ice concentration Q ares Si A F qR ey min Fy 0 1 98 pm Gi 1 97 Following Millero 1978 the equation used for the salinity dependent freezing point of sea water is Ty Tm 0 0575 S 0 001710523 S 0 0002154996 S 1 99 which is used to set a threshold to the SST 1 5 3 Coupling
102. and saves only the object files while the script and the domumentation files are already installed after the command tar xuf pipe tar 5 2 Setting up PIPE In this section it will be described how to adjust the code for a certain purpose PIPE is not developed only for one application It is programmed to allow fundamental changes in the application with only few adjustments of parameters Standard FORTRAN 7 7 is used so the model will run on many different types of computers The model was tested on workstations like SUN IBM RISC DEC or HP on mini computers like MicroVax and on main frame computers like IBM3090 Cyber205 NEC SX3 4 Fujitsu or CRAY type computers like CRAY C90 or CRAY T3D 5 2 1 The Code Preprocessor The code is written for scalar vector or parallel computers with shared or distributed memory The optimization was primarily done for CRAY type computers Since the same code cannot be used for all of type of computers control statements have been included that are used by precomp f to either activate or disable certain FORTRAN statements These control words are written into columns 73 80 The precompiler precomp f distinguishes 4 types of statements e Single or Double Precision Since the code has to work with 64 bit words to ensure an accurate solution of the matrix for the wave equation a double precision version has to be set up for a 32 bit workstation The code optionally declares all REAL variables with 64 bi
103. antity Code EBM Variables 129 5 8 Definition of Section Code Partl o 130 5 9 Definition of Section Code Part II 131 5 10 5 11 5 12 5 13 5 14 5 15 5 16 6 1 6 2 7 1 7 2 7 3 7 4 7 5 7 6 T T 7 8 7 9 7 10 7 11 7 12 Definition of Main Switches e 137 Definition of Margins 3 2 aah Ga ae sa As a Aa Oe 138 Definition of Eulerian Angles o 138 Definition of Layout Parameters 138 Definition of Interval Parameters 139 Definition of Scaling Parameters o 139 Definition of Zonal and or Meridional Means 139 Constants for Bryden s Formula The subscripts denote the exponent for pressure first index salinity second index and temperature third index See A Te LO spot tae NA 146 Constants for UNESCO Formula The subscripts denote the exponent of respective variables See also Eqn Lira ko Ye a ed ee Ao 146 Names for Prognostic Ocean Quantities o 149 Names for Diagnostic Ocean Quantities 149 Names for Surface Flux Parameters 0 0 0 0 0000000584 150 Names for Bias Correction Terms 0 00004 150 Names for Prognostic Snow and Sea Ice Variables 150 Names for Diagnostic Mixed Layer Variables 151 Names for Horizontal and Vertical Diff
104. antity Code EBM Variables 5 5 3 Portable Compressed Format If PPSF files are created as binary files these usually are not portable to another UNIX platform unless the files are created as 32 bit IEEE format which is the default on many workstations However if the PPSF files are created e g on a CRAY one needs to convert them in order to be able to postprocess data on a workstation There are two possibilities The script ppsf2i3e can be used to convert the binary files to 32 bit IEEE which gives a compression effect of nearly 50 Another option is to use the script implode It converts the binary format into a character string with 3 byte resolution If the script explode is installed onto some system then it converts the result of implode back into the local binary format on that particular system Because explode and implode use standard ANSI language only the combination implode and explode allow the transfer of data that originally are in PPSF binary format onto all computer platforms 5 5 4 List of PPSF Files 5 5 4 1 The Quick Look System The model delivers a number of quick look files in PPSF Post Processing System Format They contain model data from the end of a model run The data are sorted in the quick look files according to their meaning The contents of these files can be easily adjusted to any requirements See lt OCEPOST gt for details 5 5 4 1 1 File PPSF_SUR contains surface quantities for veloci
105. ar one can start to use JSW48 38 In this mode the HEATM PMEM TAUXM and TAUYM is used as ocean forcing however the annual mean is not updated any more This procedure allows to the model to find that fresh water flux which is required for a stable integration but avoids the negative feedback due to the Newton relaxation Other switches control the physi cal content see table 4 19 select numerical schemes see 4 20 control the output see table 4 21 and select the data source used to create FORCES see table 4 22 The usual choice is ISW82 2 which means that the potential vorticity conserving scheme is taken instead of the Crowley scheme for momentum JSW32 1 Note that it is possible to choose an open pole for the sea ice model only with JSW28 1 while SW28 2 means that an open pole is used for sea ice and ocean This choice is meaningful since currently the sea ice model is numerically more stable than the ocean model with open pole boundary conditions SW 42 determines the ini tialization mode If SW42 0 then the model tries to initialize OCDAT and FORCES when an empty OCDAT is detected If SW42 1 then the model creates FORCES independent of whether an OCDAT exists In case of JS W42 2 the model initializes the ocean state by delivering a new OCDAT but it assumes that a FORCES is already available These choices are useful during the model setup to save CPU time 4 2 2 6 Block Data which concern Numerics 4 2 2 6 1 Blo
106. be used to activate the relevant parts of the code e lt ELIMIN gt deletes sea points in the land sea mask JF LG which are connected to the defined eliminating point indices KILLI KILLJ by sea points This routine is useful for deleting undesired lakes or sea points that are not connected to the ocean region of interest However if an index is defined which by mistake sits inside the desired model domain all active ocean points may be erased Even a floating point exception may appear when dividing through a zero number of active grid points lt GRITUNE gt calls lt ELIMIN gt lt CHANNEL gt and lt SHUTOFF gt respectively to delete undesired parts of the ocean after they have been cut off by a land bridge to introduce channels to ensure that the ocean is deep enough a narrow strait is left in the model grid and to build desired land bridges lt FLAGMOD gt allows to change the land sea mask definition after it has been derived from the topography data set lt YCORT21 gt computes the y coordinates of the Gaussian grid for the T21 resolution lt YCORT32 gt computes the y coordinates of the Gaussian grid for the T32 resolution lt YCORT42 gt computes the y coordinates of the Gaussian grid for the T42 resolution lt YCORT63 gt computes the y coordinates of the Gaussian grid for the T63 resolution lt YCORT84 gt computes the y coordinates of the Gaussian grid for the T84 resolution lt YCORT16 gt computes the y coordina
107. ber of labels for x direction JMOPLT2 LXO any 0 first label for x direction JMOPLT2 LX any label difference for x direction JMOPLT2 NYS any gt 2 7 number of labels for y direction JMOPLT2 LYO 0 1 90 first label for y direction JMOPLT2 LY 0 1 label difference for y direction JMOPLT2 BMINA any gt 0 length of plot excluding frame JMOPLT2 DMINC any gt 0 height of plot excluding frame JMOPLT2 Table 5 11 Definition of Margins Variable Value Unit Default Meaning Reference LANGLE 0 1 switch for Eulerian angles JMOPLT6 ALPHA any deg lst Eulerian rotation angle in longitude JMOPLTG BETA any deg 2nd Eulerian rotation angle in latitude JMOPLT6 GAMMA any deg 3rd Eulerian rotation angle in longitude JMOPLT6 Table 5 12 Definition of Eulerian Angles 5 6 1 4 Block Data MARKS lt MARKS gt contains the relation of the section code to the type of section used to write the data Data are transferred into lt OCEPLOT gt via JMOMARK and are required to switch between different layouts of the plots as between xy xz or yz sections 5 6 1 5 Block Data UNITS lt UNITS gt contains the definition of the unit for each quantity and a different definition in case the section code is negative Data are transferred into lt OCEPLOT gt via JMOUNIT and are required for the plot 5 6 1 6 Block Data ALLIND lt ALLIND gt contains the definition for reordering the sequence of colors in the color table for the layout of
108. cal average over the ML and the underlying layer UL The procedure how the coordinate system is reorganized when an isopycnal layer becomes too heavy is described in JMO A well known problem is the determination of the stability when potential densities are defined with respect to some reference level If the reference pressure is too low the model might analyse a weakly stable stratification while a potential density defined with respect to the considered depth yields an unstable stratification To be consistent in the model only potential densities with respect to one common pressure level are used e g for horizontal pressure gradients as well as for convection and vertical mixing In the isopycnal part of the model layers are initially stably stratified If in rare cases unstable stratification arises through a combination of vertical mixing and a combined effect of the nonlinearity of the equation of state and temperature salinity transport then the concern ing water mass is assigned to one of the neighbouring layers such that the isopycnal coordinate system becomes valid again 1 5 The Snow Sea Ice Model Sea ice is an important boundary condition for the ocean at high latitudes The seasonal cycle of ice thickness and ice extent influences the heat budget at the ocean surface and the internal stratification During cold periods freezing ice ejects salt into the mixed layer and thereby contributes to the production of heavy water Duri
109. ch a term inconsistently so that its explicit does not fit with its implicit formulation the result is unpredictable This con cerns mainly the routines lt MODEXP gt and lt MODIMP gt and lt MIXEXP gt and lt MIXIMP gt The same is true for the ADLR scheme If changes in lt SOLVERX gt are not consistent with those in lt SOLVERY gt or lt MOLVERX gt becomes inconsistent with lt MOLVERY gt then the model might behave unphysical Basically these pairs of routines must solve the identical mathematical problem 5 5 Postprocessing The data postprocessing is based on the self developed PPSF which means Post Processing System Format The goal is to write all key quantities onto a block of data that allows an unique identification of each quantity For the purpose of quick postprocessing see quick look files data are written out in simple binary format To allow transfers to other computers via a network like ethernet a conversion to 32 bit IEEE format is possible see script ppsf2i3e Another possibility is to compress data to a 3 byte accurate portabe file format see script implode This allows the user to reduce the file size by a factor of about 3 The major underlying idea is to use several types of codes to identify a quantity to include the land sea mask for data analysis and plotting purposes and to include the coordinates that are required to determine the location of each grid point Postprocessing routines such as the plot p
110. ck Data TUNING contains the weighting coefficients that define how the old and new time level is weighted within the implicit algorithm Note that these parameters must be between 0 5 and 1 0 to ensure stability see table 4 23 The over and underrelaxation coefficients are of empirical origin NITMIN defines the number of iterations that are used under all conditions In fact twice the number of iterations is taken everywhere NITMAX defines the iteration count that enables the model to detect insufficient convergence of one of the equations and dependent on SW35 causes an error exit see table 4 24 Threshold parameters finally allow one to limit certain variables mainly for stability purposes see table 4 25 Variable CMIXO CMIX1 CMIX2 CMIX3 CMIX4 CMIX5 CMIX6 CMIX7 CMIX8 CMIX9 Value 0 1 E6 0 25 2 5 1 25 0 1 E6 0 5 0 01 1 0 Unit m s m 3 Meaning linear damping term of ML constant TKE decay length scale critical Richardson number Ekman dissipation for TKE wave energy conversion factor Garwood term tuning coefficient constant buoyancy decay length scale Ekman dissipation for buoyancy flux threshold for hg tuning coefficient for buoyancy flux Reference lt MIXEXP gt lt MIXEXP gt lt MIXEXP gt lt MIXEXP gt lt MIXEXP gt lt MIXEXP gt lt MIXEXP gt MIXEXP gt MIXEXP gt DETIME TURBEN TURBID TURBREL PENET 0 0 15 0 42 20 e MIXEXP gt MIXEXP gt
111. cross referenced flux components This yields 4 p At pn 1 2p 1_ AEF OR 12m 1 aoe puh poh 1 2 i 2 4 where Fin and Fion are the abbreviations of sh puh Foun Sf 0 ovh Foun 2 37 puh At A a ae vh uh p p EF J af Oo Foun 2 38 a f 02 If the continuity equation is discretized in the same manner as the momentum equation this yields At o h t A N Ayr Fr oh oh 5 ou T where Fin represents all the remaining terms taken at time level n If the flux divergence now is eliminated by using 2 35 and 2 36 an equation for the layer thickness h only is obtained nytt 2 39 E areta Pee hg hg nytt h e ae FR ppt oe pe 2 40 At a ae OF sts At ale 44 A ATREA FAO Ada AV ART ALF 2 0 4 01 AR f4 O2 40y n 1 A Byl AP F O2 4 ay APO REAL EOI 4 dy1 AC f 02 400 n The following identity is used to separate p and h from the product ph Artl h oe ans p A h het Peo ee ae 2 dol 2 This relation is also used to split the heat and salt content up into their basic quantities which are potential temperature and salinity The quantity Fi denotes the explicit part of the continuity equation in the semi implicit technique defined by At an Pr 2 42 ph ph 2 a uh T y pah Pht ph P t 2 41 xn n Together with the equation of state 1
112. ct f and a collection of routines in oceplot f The object code is created by the command modgen sun plot or by modgen cri plot Currently the whole plot system requires the plot library jmoplane f a small C program utility c and a GKS library like provided by NCAR Thereby the simplest version of a GKS ststem the GKS0A level is sufficient If a GKS system is not available jmoplane f can be compiled without the directive DNCAR Then the resulting executable creates a PostScript file which might be viewed with e g ghostscript 5 6 1 1 Main Program PICTURE The main program lt PICTURE gt consists only of one call to lt OCEPLOT gt which is the driving subroutine It reads data block by data block and creates a plot that is dependent on the code numbers required to identify a quantity 5 6 1 2 Block Data PLOTPAR The parameters of table 5 10 control the overall behaviour of the plot program Features as number of to be plotted pictures plot quality type of plotter laser or color plotter etc can be chosen The parameters of table 5 11 define in more detail the size content and lables used for plotting Finally table 5 12 allows to define the Eulerian angles essentially to add lines of geographical longitude and latitude on each frame Reference JMOPLT1 JMOPLT1 JMOPLT1 JMOPLT1 Variable Value Default NPIC any gt 1 999 NCOL 1 0 1 2 1 NTURN 0 1 0 NSHIFT any 0 Meaning number of pictures to be plotted
113. ction term In the case of open boundaries this routine might be used to define aswamp for momentum A stop statement remembers the user to adjust this entry to the user s application e lt CROSMIX gt determines the upward and downward mass flux due to diapycnal mix ing The result is added to VERTUP and VERTDN These arrays are later used in lt TRACTIV gt to compute the related changes in temperature and salinity concentration The contribution to the momentum budget is added to RESTU and RESTV e lt ISODIFF gt computes the diffusion coefficient from the shear flow with a specified threshold minimum value e lt DIADIFF gt computes the diffusion coefficient for the layer diffusion e lt DIVADD gt adds the flux divergence of the continuity equation to RESTH e lt DIVFOR gt computes parts of the right hand side of the wave Eqn 2 40 and adds the contribution to RESTH Compare with Eqn 2 37 and 2 38 e lt FLUXES gt computes the new fluxes from the pressure gradient at the new time level and the already determined right hand side of the momentum equations that are not treated implicitly Compare with Eqn 2 35 and 2 36 e lt VELOC gt computes the velocities from fluxes lt LIMVEL gt limits velocities to UVMAX Entry lt LIMFLX gt limits the according mass flux lt SETGRADX gt computes as differences on two latitudes the part of the pressure gradient that results from layer thickness gradients The differenc
114. d point 4 5 9 Routines for Data Manipulation lt STORE1 gt copies a 2 dimensional real array onto another one lt STOREN gt copies a 3 dimensional real array onto another one lt KTORE1 gt copies a 2 dimensional integer array onto another one lt KTOREN gt copies a 3 dimensional integer array onto another one lt VALADD1 gt adds a 2 dimensional real array onto another one lt VALADDN gt adds a 3 dimensional real array onto another one lt VALSHFT gt adds a constant onto a 2 dimensional real array lt VALMUL1 gt multiplies a 2 dimensional real array with some factor lt VALMULN gt multiplies a 3 dimensional real array with some factor lt SETVAL1 gt sets a 2 dimensional real array to a constant value lt SETVALN gt sets a 3 dimensional real array to a constant value lt KETVAL1 gt sets a 2 dimensional integer array to a constant value lt KETVALN gt sets a 3 dimensional integer array to a constant value lt SELECTIDE gt extracts a specific tidal mode from an array containing all tides 4 5 10 Routines to Shift Data e lt ROLLZXY gt shifts 3 dimensional real data in x direction in order to improve conver gence e lt KOLLZXY gt shifts 3 dimensional integer data in x direction in order to improve con vergence e lt ROLLXY gt shifts 2 dimensional real data in x direction in order to improve conver gence e lt KOLLXY gt shifts 2 dimensional integer data in x direction in order to improve con vergence
115. d expansion coefficients with respect to temperature and salinity 7 the zonal wind stress Q the y component of the earth s angular velocity and I 1 in the presence of ice otherwise I 0 The variables m1 ma m3 M4 M5 and mg are tuning coefficients which in the code are called CMIX with an extension The favorite parameter choice of JMO is m3 M4 Me 0 and m 1 The parameter g limits the entrainment rate to a finite value when g 0 The entrainment Eqn 1 18 contains numerous parameters These allow the use of many different types of mixed layer models that have been used in literature As an ex ample the simplest possible one from Kraus and Turner 1967 requires the following pa rameters CMIX0 CMIX2 CMIX3 CMIX5 CMIX7 0 CMIX9 1 CMIX1 CMIX6 1 E6 TURBEN TURBREL 0 This parameter set reduces the entrainment equation to wg h 2m u2 hB The entrainment detrainment rate w is related to the transfer rates wi in the Eqns 1 1 to 1 4 by Kias wa u 122 The sea ice model described subsequently has some influence on the relation for entrainment If sea ice is present B is taken as the buoyancy flux and g as reduced gravity The buoyancy flux B from the atmosphere into the mixed layer is expressed as g S1 Si Da B qQ 1 A _ 1 2 O Fen 1 28 where S is the ice salinity Q is the heat flux from the ice into the mixed layer see Eqn 1 94 Q is the heat flux from the atmosphere into t
116. d interactively using plot b p e q 68 a Annual Mean PPSF_SUR Potential Temperature NZ 0 Id Ro OOnNDOFP WHR L 10 11 o 12 13 14 L 15 L 16 L 17 L 18 L 19 L 20 21 L 22 L 23 L 24 L 25 L 26 ao L 28 L 29 L 30 31 L 32 T 905 605 305 0 30N 60N 90N Unit Celsius SS 1000 750 500 250 0 3500 6000 Figure 5 3 Example showing a yz section for the potential temperature at 120W from the PIPE T42 The picture is created interactively using plot b p e q 13 1 240 PPSF_VER Sea Ice Thickness Surface O1FEB 007 T T T T T T I T T 00 05 10 15 20 25 30 35 40 45 50 55 Unit m Josef M Oberhuber DKRZ Fri Mar 12 13 38 14 1999 C JMOPlane Figure 5 4 Example showing a north polar projection for the sea ice thickness of PIPE T42 The picture is created interactively using plot b p e r q 18 N PPSF_ICE Sea Surface Velocity Surface O1JAN 003 57 I I I 0 0 10 20 30 40 50 60 70 80 90 10 0 11 0 12 0 13 0 14 0 15 0 Unit cm s Josef M Oberhuber DKRZ Fri Mar 12 13 22 40 1999 C JMOPlane Figure 5 5 Example showing the surface flow in a North Atlantic Arctic model using Eule rian angles ALPHA 40 BETA 50 GAMMA 0 Geographical longitude and latitude
117. dimension 256 128 File T106GRID has dimension 320 160 5 4 3 1 6 Open Boundary Forcing File SALTEMP JAN is the 3 dimensional observed temperature and salinity for Jan uary Format is identical with SALTEMP The local file is STMONO1 File SALTEMP FEB is the 3 dimensional observed temperature and salinity for February Format is identical with SALTEMP The local file is STMONO02 File SALTEMP MAR is the 3 dimensional observed temperature and salinity for March Format is identical with SALTEMP The local file is STMON03 File SALTEMP APR is the 3 dimensional observed temperature and salinity for April Format is identical with SALTEMP The local file is STMONO0A4 File SALTEMP MAY is the 3 dimensional observed temperature and salinity for May Format is identical with SALTEMP The local file is STMONO5 File SALTEMP JUN is the 3 dimensional observed temperature and salinity for June Format is identical with SALTEMP The local file is STMONO6 File SALTEMP JUL is the 3 dimensional observed temperature and salinity for July Format is identical with SALTEMP The local file is STMONO7 File SALTEMP AUG is the 3 dimensional observed temperature and salinity for Au gust Format is identical with SALTEMP The local file is STMONO8 File SALTEMP SEP is the 3 dimensional observed temperature and salinity for September Format is identical with SALTEMP The loc
118. dimension of the array N2 is the second dimension of the array N3 is the third dimension of the array Note that this allows to write out several 2d arrays that require the same land sea mask and have the same coordinates KODEL is a key number that defines the location of a section e g as depth in meter or longitude latitude in degree IFLG is the land sea mask with dimension N1 N2 In the ocean model IFLG is defined as INTEGER array while on an PPSF file it is written out as CHARACTER 1 KTOTAL is the model time in seconds MSEC is the model time in seconds after the full minute MINUTE is the model time in minutes after the full hour 11 12 13 14 15 16 17 18 19 Code Number Variable total surface heat flux solar radiation longwave radiation sensible heat flux latent heat flux net fresh water flux as salt flux surface buoyancy flux solar buoyancy flucx friction velocity drag coefficient Reference QFLUX QSOLAR QLONG QSENSIB QLATENT SFLUX BFLUX BSOLAR USTERN DRAG lt VARDRAG gt lt VARDRAG gt DQDT makeall f makeall f makeall f makeall f lt SURFLUX gt lt SURFLUX gt lt STRESS gt transfer coefficient for sensible heat transfer coefficient for latent heat feedback coefficient vertical velocity layer relative vorticity layer absolute vorticity layer potential vorticity heat flux bias correction fresh water flux bias correcti
119. dinates The differential operators are taken along the interfaces They run essentially in the horizontal The divergence operator V is given by v sp 2 7 r cos pad r cos pay dl where A is the longitude y the latitude and r the earth s radius Similarly the gradient V is 3 0 a 2 8 r cos yO rOp Due to the use of Lagrangian coordinates vertical derivatives do not appear The Laplacian operator for scalar quantities is 3 1 ge 1 ne Ne r2cos p IN E r2cos p Be ap oe The diffusion operator of a vector quantity u v on the sphere is vn 1 o E 1 sn v E r2cos p gr r2cos p Op A ap 2 1 _ iatan y u 2sin p ie 2 10 qa r2cos p OX 3 1 o 1 V Vvu 0 r2cos p ar ij r2cos p Fas T 1 tan p u i 2sin p 0 2 11 r2 r cos g T 2 1 1 3 The Pressure Gradient hs r41 3 Rs BSU J PS0 J41 e O A a PO o MR A ee 2 12 Ox VLJ 2ATV 7 J hs hs r 1 Agen hscrsa re I J 1 141 341 LJ I 1 J 2 13 Oy 2Ayy 1 3 The same discretization is used for og Note that the pressure gradient within each layer is the sum of layer thickness gradients topography gradients and potential density gradients 2 1 1 4 The Flux Divergence puh y r 1 1 ouh y r y 2Azr5 7 E puh y 11 11 puh y r 1 7 2 14 2Az 57 cos y yi evh y r1 7 evh a cos y s 1 costp y 1 Ays y cos p y 3 1 pvh y r1
120. e the file OCDAT and all PPSF files appear in the same style The data transfer between server and client is carried out using parallel pipes see Fig 2 2 If the domain is decomposed into pe_z pe_y subdomains each executed on one processor the horizontal grid is separated such that each data piece contains only those latitudes which 1 pe_y is responsible for The data piece still has full dimension in x direction If these data pieces are sent to the according processor 1 pe_y this processor only needs to further split the data in x direction such that these data pieces exactly for with the grid on the processors pe_z pe_y There does not exist severe bottlenecks on the network Therefore using the concept of parallel pipes allows an efficient data transfer between server and client and does not cause pronounced latency times till communication is finished On the MPP each processor holds one extra longitude to the west and to the east and an extra latitude to the south and to the north Whenever there is a call to one of the boundary routines in the PVP or RISC version of the code these are the locations where data have to exchange in order to provide updated boundary values for computing gradients or divergences The southernmost latitude on the southernmost processor row and the northernmost latitude on the northernmost processor row are treated differently Fig 2 3 shows a sketch for how data are exchanged Data Pipes CFLs
121. e File ECSST contains the observed SST data from the ECMWF reanalysis Project on a T106 grid It is read by lt ECSST gt from the local file SSTEMP with the following simplified statements DIMENSION ISST 360 SST 360 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 CISST I I 1 360 DO I 1 360 SST 1I J MONTH ISST I 200 273 16 ENDDO ENDDO ENDDO e File ECAIR contains surface air temperature data from the ECMWF reanalysis Project on a T106 grid It is read by lt ECTSC gt from the local file AIRTEMP with the following simplified statements DIMENSION ITAIRT 360 AIRT 360 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 IAIRT 1 I 1 360 DO I 1 360 AIRT I J MONTH IAIRT I 80 203 15 ENDDO ENDDO ENDDO e File ECCLD contains the cloudiness from the ECMWF reanalysis Project on a T106 grid It is read by lt ECCLD gt from the local file COVER with the following simplified statements DIMENSION ICLD 360 CLD 360 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 ICLD I I 1 360 DO I 1 360 CLD I J MONTH ICLD 1 10000 ENDDO ENDDO ENDDO e File ECHUM contains the surface humidity from the ECMWF reanalysis Project on a T106 grid It is read by lt ECHUM gt from the local file WETNESS with the following simplified statements DIMENSION IHUM 320 HUM 320 160 12 DO MONTH 1 12 DO J 1 160 READ 2014 IHUM I I 1 320 DO I 1 320 HUM I J MONTH IHUM 1I 10000 ENDDO EN
122. e is created Note however that all keywords vanish as well Further control words occur in the code that are not touched by precomp f These are OPTION to mark optional changes in the code and DANGER for statements that are crucial for the model numerics precomp f is called by the script setup These UNIX scripts have to be edited once to define the model setup see comments in setup In addition they change the dimensions of the arrays by the batch editor sed see further comments on how to specify the array dimensions Finally these scripts generate the code for a certain computer and with a specified dimension 5 2 2 Defining the Grid If one of the T21 T32 T42 T63 T84 or T106 grids is chosen the work is already done by precomp f Only the proper dimension have to be chosen see section 5 3 1 If the grid is self defined then the corresponding switch in precomp f has to be set accordingly The model margins are then defined by XGL XGR YGU and YGO Note that the arrays X and Y that contain the x and y coordinates are defined on velocity points and that X 2 XGL X NX 1 XGR Y 2 YGU and Y NY 1 YGO 5 2 3 Specification of the Dimensions For the case of a constant grid spacing Az and Ay the following procedure is used to determine the dimensions NX and NY XGR XGL NS AAA i Na T 5 1 YGO Y Na ALLE TOO ys 5 2 Ay For the purpose of specifying the boundary conditions there is always an e
123. e to run sequentially on an MPP PIPE scales excellently If overall 64 processors are used even a low resolution version of PIPE like the T42 with 128 x 64 grid points on a CRAY T3D achieves a nearly theoretical speedup when comparing with a one processor version A higher resolution version with 512 x 256 grid points e g achieves a speedup of 436 on 512 processors with about a CRAY C90 equivalent Mflop rate of 24 The reasons for this excellent behaviour are the use of an implicit scheme which performs computational intensive work and requires little communication Currently there exists the restriction that the MPP version of PIPE does not allow for open boundaries and does not allow to use the tracer module These pieces of the code will be ported to the MPP architecture upon request Coupling of Sub Domains Figure 2 3 Concept for Exchanging Data between Processors parallel even odd line relaxation on neighbouring PEs l parallel solution of tridiagonal systems step gt Figure 2 4 Concept for Iterating on Neighbouring Processors Part II The PIPE System Components Chapter 3 Flow Diagrams The PIPE model has a pretty simple and easy to understand program flow It to some extend is controlled by switches consisting of JS Wand an appended number Either these switches have to do with initialization postprossing the physics or the numerical scheme The driving routines are lt OCEINIT gt lt OCEST
124. ecause parameters depend on each other a hand based multi variational approach was used to find an optimal set of parameters Only annual mean fields have been used to optimize the parameters Well known drawbacks of the EBM are that 1 due to a missing surface boundary layer heat and fresh water fluxes are too strong in the subtropical gyres 2 due to missing upper atmosphere northward heat transport due to the Icelandic Low the heat fluxes into the Greenland Iceland Norwegian Sea are heavily underestimated 3 due to missing vegetation soil becomes dry and hot during sommer 4 due to missing winds above the atmospheric boundary layer the heat and water transport over larger distances might be underestimated 1 7 Tide Model 90 N T 30 N 30 S T 60 S Y Lo O T T T T T T 30 E 60E 90E 120 E150 E180 150 W120 W90W 60W 30W 0 30 E Figure 1 3 Amplitudes of the M2 Tide simulated with a T106 equivalent version of the barotropic tide model coupled to the 3 dimensional PIPE model Tides represent high frequency variability that might have important impacts on regional scales This is because tides become more important with decreasing water depth due to increased interaction between the large scale mean flow tidal flow and topography Thus they might 90 N T 60 N T 30 N I T T T T JOE 60E 90E 120 E150 E180 150 W120 W90W 60W 30W 0 30 E Figure 1 4 Phases of the M2
125. ecision but the postprocessing should be done with single precision Script dpr2spr converts a binary PPSF format that has been created with double pre cision into a binary PPSF format with single precision Script implode converts a binary PPSF format to a portable file with reduced accuracy i e 3 byte resolution This file can be ported to any computer system and unpacked with explode Script explode converts the portable file with reduced accuracy i e 3 byte resolution to a binary PPSF format Dependent on the operating system it either will create a 32 bit IEEE format or like on PVP computers a 64 bit binary format Script ppsf2i3e converts a binary PPSF format to a 32 bit IEEE format in order to make PPSF files created on a PVP computer portable to a 32 bit operating system This script makes sense only to execute it a PVP computer Script i3e2ppsf converts a 32 bit IEEE format to a binary PPSF format This script makes sense only to execute it a PVP computer Script modify modifies one of the quantity section or location codes in the file headers of the PPSF format Script reduce extracts a fraction of the model domain and creates new data in PPSF format with adjusted dimensions in the header and the data records Script fill fills in extra latitudes south of the first latitude of a data set This makes sense only for horizontal fields where for the sake of saving CPU time the southernmost latitudes on An
126. ed After having set ISW45 2 the file LSMASK is read and used to overwrite the internally initialized land sea mask This switch setting also must be used if one of the T21 through T106 land sea masks are used in connection with the according pre defined grid see section 5 3 1 In order to prepare the orography on land to allow the river runoff model to produce meaningfull results lakes i e depressions that are not connected with the ocean are eliminated Furthermore the standard deviation of the fine resolution raw data is used to properly define rivers 5 2 6 Selecting the Forcing Data As explained later in detail the model offers the choice to force the model with either the COADS data or ECMWF analysis Two topography data sets are offered one with 1 degree and the other with 5 minutes resolution The choice of the topography data set depends on whether the roughest possible topography should be used in an experiment The COADS data set can be used for regional models as for the North Atlantic The drawback of the COADS is that they are not spatially complete Even if the model fills up holes within a data set by artificial but physical meaningfull data a global model should use the ECMWF data The chosen data must be defined as switches in the code see table 4 21 as well as switches in the batch job model sun or model cri that runs the model The data format in general is 2014 For details see references given in 5 4 3 1
127. el is disabled in that mode the time step can be chosen significantly higher because of less restrictive CFL numbers for the tracer transport equation This contributes to a further significant saving effect The equations are a r y E OL V 7 C ph x Oe Dres V gt A ph V Cx RE we wpC k wpC i wpC k 1 124 where Cph is the tracer mass R is the tracer forcing and Upes is the time averaged residual flow of the barotropic tide model All other variables are already explained in section 1 1 1 By default an early version of the tracer model is part of the collection of codes This version is not conservative A newer conservative version is available however as it is delicate to use it will be distributed only upon a request and guided with advice Chapter 2 Model Numerics 2 1 The Interior Ocean and the Mixed Layer Model 2 1 1 Discretization in Space on the Arakawa B Grid The notation is that K is the 1st index J the 2nd and J the 3rd A grid point with the indices 1 1 1 is located in the top layer the mixed layer and in the south west corner of the model grid The B grid distinguishes between scalar and vector points as shown in Fig 2 1 In the subsequent formulae scalar quantities have a subscript S K I J while a vector point is O O scalar grid points yL vector grid points land x vector grid points sea O oO O O Figure 2 1 Schematic layout of the
128. emperature is related with in situ temperature T in situ pressure P and salinity S through the formula by Bryden 1973 0 T P aw awo T ao T da T aiio S a S T Pe a200 4210 S a201 T a202 T P a300 a301 T 1 16 where P x S S 35 and 7 T 273 16 All coefficients a with subscripts are defined in Appendix C The temperature T is computed from the prognostic potential temperature 0 by inverting the formula by Bryden using a Newtonian technique The in situ density is computed using the UNESCO formula for sea water which is defined through the following relations Pw atT a 4 T a2 T a3 T a4 T as kw T e1 T eo T e3 T ea dy ho T hy T ha T hs by ko T ki Tka b by S mo T m T ma a Aw S io Ta T i2 joS ky Lo Tf TL T f3 5 90 T g1 T g2 z ro P a P b Po Put S bo T b T bo T bg T b4 doS S coT c1 T c2 p LL 1 17 o oll II Tp where now P mm S Vy S and T T 273 16 All unexplained coefficients are defined in Appendix C In order to compute the potential density o the in situ density p is computed using the UNESCO formula however on the reference pressure level REF PRES In order to use the correct temperature 7 also the formula by Bryden is inverted using the same reference pressure level Whenever a heat and fresh
129. er model is used the file TRDAT is needed in addition If the tracer model is running in offline mode then also file TRAFORCE is needed 5 4 4 The Multi Grid Method The model allows the user to jump between different grids within a model run For instance this can be useful if a spin up run is carried out with a coarse resolution but the production experiments require higher resolution Since the coarse resolution model can be run to a stationary state much faster than a fine resolution model this provides a convenient way to save CPU time Before the model is compiled the variables NXOLD and NY OLD must be set to the dimensions of NX and NY that have been used for the coarse grid This can be defined in setup while creating the model source code Note that currently there is the restriction that the vertical dimension NZ of the old and new model grid must be identical However widely untested modifications already exist for allowing to decrease or increase NZ Simply try and if necessary debug lt OLDINIT gt In order to use this method initialize the model for the new model grid as explained above However ensure that the OCDAT of the old run is saved as OLDFILE on the data directory that the parameter JU MGRID 1 in model sun or model cri and that J SW12 1 in ocemain f This causes the code to interpolate the old grid onto the new grid and onto the new topography If there are no data available from OLDFILE for interpola
130. es are used to compute the right hand side of the linear equation system of the wave equation It also determines the weighting coefficients that are required to compute the matrix for the wave equation This routine is used by lt SOLVERX gt only lt SETGRADY gt computes as differences on two latitudes the part of the pressure gradient that results from layer thickness gradients The differences are used to compute the right hand side of the linear equation system of the wave equation It also determines the weighting coefficients that are required to compute the matrix for the wave equation This routine is used by lt SOLVERY gt only lt GRDADD gt computes explicitly that part of the pressure gradient that is due to layer thickness and density gradients The result is added to RESTU and RESTV lt LAPRADD gt computes the contribution to the right hand side of the wave equation that contains the gradients of density as a Laplacian operator The routine updates its contribution for the new layer thicknesses during the iteration lt VORABS gt updates the potential vorticity lt LAYDIFF gt determines explicitly the layer thickness changes that result from the layer thickness diffusion and adds the result to RESTH The same term is treated implicitly as part of the matrix within lt SOLVERX gt and lt SOLVERY gt lt LAYFLUX gt corrects the horizontal transports by those contributions that result from layer thickness diffusion Thi
131. esidual error that can lead to instability Based on the experience of the author a few suggestions on how to set the numerical parameters properly should be given To visualize the models behaviour the following lines appear in the protocol OUTLIST a gt Wave Equation gt gt Momentum Advection Diffusion gt Temperature Equation t gt gt Salinity Equation gt Iterations gt Ice Thickness amp Concentration gt Ice Momentum gt gt Used Time Step gt Mean Explicit Acceleration gt Computed Change of Flux gt Computed Change of H gt Total Momentum 5 w www w LLMM w www w o _ w www w w www w 0 500 865 166 281 1759 174 5 9 11 071 006 47 85 2976 0 500 904 182 282 1760 169 5 9 11 071 006 75 55 2934 0 500 935 179 280 1761 170 5 9 11 071 006 76 95 3009 0 500 924 175 279 1763 170 5 9 11 071 006 79 71 2984 0 500 935 182 278 1765 168 5 9 11 071 006 65 91 2918 Found negative H Points lt Maximal found negative H lt Z X Y Index of Point lt Number of Cells being opened lt Number of Cells being closed lt Number of total Convection Events lt Each time step writes a line of critical numbers into OUTLIST These
132. exp h hie 127 exp h hue if B lt 0 i exp h Rsuo if B gt 0 1 28 where hike and hy are the dissipation length scales for negative and positive buoyancy fluxes In the retreat phase of the mixed layer the depth is determined by setting w 0 in Eqn 1 18 and solving for h The resulting equation for the Monin Obukhov length hy for the case ma M4 Mg 0 and ms 1 is 2m u7 hu B yBs y1Bs huU ezxp hy hp 2hp 1 exp hy hp 0 1 29 The solution for hy is determined iteratively with Newton s method A second diagnostic calculation is carried out as soon as the flow becomes unstable due to excessive vertical shear In this case a minimum depth hp is defined by Uy wi4 01 114 hri Riera 7 g 1 30 As a result there are two constraints that limit the MLD First the MLD cannot be deeper than the Monin Obukhov length and second it cannot be shallower that hp If the two constraints contradict each other hp is taken as criterion In addition the mixed layer is not allowed to be shallower than some threshold value hmin so that finally hy max hmin hri min hy hi 1 31 where hj is the mixed layer depth computed prognostically through Eqn 1 18 This means that as long as the MLD is thinner than hy and thicker than hp the MLD is not altered by these diagnostic calculations but as soon as the MLD becomes thinner than hp or thicker than hy the corresponding diagnostic quanti
133. f a vector point The brightest shades denote sections through Antarctica and North America 1 1 1 Ocean Model Equations The basic equations are formulated in flux form as conservation equations for the vertical mean of the mass flux pvh the mass content ph the heat content ph and the salt content Sph y in a column of the k th layer gt EE 5 gt gt gt A Sesha a gh V puh hy Vpn f x puh V A ph V 0 Ta EE cos Saik gt k EIE F wpv wpv a WPU py wpv _ Te Ty FT 1 1 oh F poh A paa Rh phe edi pda wo 1 2 A nw S F wo wp0 k wp Wwp0 k 1 3 2 sph DS ph a Br F APM VS RE wpS f gt wp5 k wpS k_ we 1 4 The terms wpv wp wp0 and wpS describe exchange processes of momentum fluxes mass itself heat and salt content between neighbouring layers g are barotropic velocities due to Greg Holloway s Neptun effect and v es is the time averaged residual flow computed by the barotropic tide model The different kinds of exchange processes are entrainment detrainment vertical exchange as cross isopycnal mixing and convection The terminology indicates a transfer from the l th layer into the k th layer where l k represents the next upper and I k the next lower physically present layer N is the number of layers the index k starting with k 1 in the uppermost layer which always
134. fining the first block number the last one and an increment e Script filter extracts all data blocks with a prescribed quantity code KODEQ section code KODES and location code KODEL from an input file e Script insert inserts data at a specified location of an input file and writes the result to an output file e Script cycle performs a cyclic move of the data from the end to the start of the input file or vice versa e Script annual computes the climatological seasonal cycle from an input file containing many years of data e Script seasons computes seasonal means like DJF MAM JJA and SON from an input file containing many years of data Script area computes an area average for non zero data Script curve extracts a time series of means Script blank extracts either the Atlantic the Pacific or the Indic Ocean basin data from an input file containing e g global data Script clean sets the coordinates and flag fields of all subsequent data to those of the first data record All data must be of the same type 5 5 6 3 Utilities for Manipulating Data Contents Script add adds the contents of two input files and writes the result to an output file Script mult multiplies the contents of two input files and writes the result to an output file Script divide divides the contents of the first input file by the contents of the second one and writes the result to an output file Script differ computes t
135. fline mode in order to compute further tracers Adjustment to specific computer architectures e The code now runs efficiently on vector and RISC processors e The code now runs on parallel shared and distributed memory architectures Changes in Pre and Postprocessing e The installation is highly simplified e The data postprocessing now is supported by a big number of usefull tools Changes in Supplied Data e The ECMWF reanalysis data are available now Part I Model Description Chapter 1 Model Physics and Dynamics Remark This manual is based on the model initially published by Oberhuber 1990 1993a 1993b hence forth called JMO An updated version was documented by Oberhuber 1993c Updated ver sions have of OPYC have been used for a couple of key studies like those published by Miller et al 1992 Miller et al 1994a b Aukrust et al 1995 Holland et al 1996a b Lunkeit et al 1996 Cherniawsky and Oberhuber 1996 Roeckner et al 1996 Bacher et al 1998 Zhang et al 1998 Timmermann et al 1998 1999 or Oberhuber et al 1998 or have produced significant results during a BMBF project Nies et al 1998 for the problem of ra dioactive tracers in the Arctic oceans This manual describes the full source code of OPYC Cycle 3 10 which is called the Parallel Isopycnal Primitive Equation PIPE model In the first chapter the model is described in the second chapter the model numerics is discussed in the
136. gnose a variety of model quantities e Script makeall computes the horizonal transport stream function and the vertical over turning stream function Optionally the overturning stream functions for the Atlantic the Pacific and the Indian Ocean are written onto an output file It also computes hori zontal and vertical sections of temperature salinity and many other quantities Use the command makeall h in order to obtain an overview As then listed makeall also computes quantities like electro streamfunction or 3 dimensional magnetic potentials as part of a project with APL Seattle through Robert Tyler Warning The dimension statements must be adjusted to the particular settings of the current model version e Script ppsf2grd converts the PPSF format to the data structure expected by the GrADS System This possibility is offered but not supported Warning The dimension statements must be adjusted to the particular settings of the current model version 5 6 The Plot System The plot software reads files in binary PPSF format The header of each data block specifies how a quantity will be plotted The intention for the plot system is to work without additional changes to produce any chosen picture unless the subsequently listed control parameters have been specified by a new model layout 5 6 1 The Plot Program Similar to the ocean model the plot program is split up into a simple main PROGRAM and several BLOCK DATA BLOCKS in ocepi
137. he black box is accomplished with the help of the coupling interface routines If tracers are computed then the routines lt TRAINIT gt lt TRASTEP gt lt TRAPOST gt and lt TRASTOP gt are called as well 4 2 2 Definition of Control Parameters The meaning of the control parameters is listed in each of the BLOCK DATAs Henceforth it is distinguished between default values and constants Default values are marked with a and should be understood as recommended and well tuned values for most purposes If any appears in the column for Value then there is no limitation except when further restrictions are listed Constants such as the specific heat capacity of water are listed without a and 63 never should be changed Values such as 0 1 2 represent lists of allowed values for switches Other choices lead to unpredictable results In most cases routines are referenced however if a variable is used at many locations of the code the COMMON block that encloses the variable is referenced In general the code strictly follows the FORTRAN convention that variables starting with 1 J K L M or N are INTEGER variables while all other variables are REAL except for a few cases of CHARACTER definitions The subsequent tables are not complete For a complete overview see the descriptions in ocemain f and don t be too lazy to simply look where some variable switch or whatever appears in the code Use an ASCII editor to understa
138. he difference of the contents of a first input file with that of a second one and writes the result to an output file Script triade computes the sum of the first two input files multiplies the result with the contents of a third input file and writes the result to an output file Script addc adds a constant value to all quantities of an input file and writes the result to an output file Script multc multiplies a constant to all quantities of an input file and writes the result to an output file Script curl computes the rotation of vector input data Script div computes the divergence of vector input data Script grad computes the gradient of scalar input data Script norm computes the norm of vector input data i e the result of yu v if u v is the input vector Script abs computes the absolute value Script min computes the minimum of subsequent data of the same type Script max computes the maximum of subsequent data of the same type Script mean computes the mean of subsequent data of the same type Script int computes the integral of subsequent data of the same type Script average computes the mean of subsequent data In opposite to mean the input file now may contain different quantities The output file contains the averages for each quantity Script merge computes the mean of subsequent ocean layer raw data only separately for each quantity Other quantities are forwarded to the output fi
139. he mean of the absolute values of the computed change of the flux and layer thickness The next row shows the total content of momentum The rows 12 14 are important Since the model works with time variable boundary conditions it formally may happen that the layer thickness becomes negative due to a too strong divergence at a grid cell that has already lost its entire mass The numbers are extracted from the predictor step In the corrector steps negative values are cancelled However these first guesses are an indication of a too large time step It is important that even residual negative layer thicknesses are listed in the sum A rule should be that the negative layer thicknesses must never appear as since this indicates fundamental stability problems The indices in the order z x y indicate where the problems are located Finally the next two rows list the number of cells that have been opened or closed per time step This is also some measure for numerical noise The last row shows how many horizontal grid points experienced a convection event 5 4 9 Numerical Filters In order to suppress instabilities the model has a number of mechanisms to reduce numerical noise These are e The parameter SMOOTH enables filtering of all interface heights This parameter must be small otherwise one violates the first law of ocean modelling namely to never touch the pressure field One can create arbitrary nonsense if SMOOTH is too big Howeve
140. he mixed layer see Eqn 1 34 q is the sea ice concentration see Eqn 1 78 and 007 05 is the expansion coefficient with respect to salinity This expression implies that the heat flux through the ice is associated with a fresh water flux because ice keeps the mixed layer temperature at the freezing point Furthermore with B 0 it is assumed that no solar radiation penetrates through the ice In the presence of ice i e T 1 the definition of the reduced gravity g through Eqn 1 19 used in the entrainment Eqn 1 18 is corrected by ce 091 4 ee C Sy Si 004 O06 p KAPON y AA 1 24 oh A og Se B a ll 0 1 24 Compared with the entrainment equation without sea ice additional mechanisms have to be considered Entrainment provides some heat flux into the mixed layer Under the presence of sea ice this flux is exactly balanced by a heat flux due to melting or freezing of ice which keeps the mixed layer temperature at the freezing point The connected fresh water flux induces a buoyancy flux that changes the amount of turbulence available for mixing The first term on the left hand side of Eqn 1 18 describes the production of mean po tential energy MPE and the second describes the production of mean kinetic energy MKE by entrainment On the right hand side the first term represents the TKE production due to wind stirring u denotes the friction velocity and the following treats the influence of the surface buoyancy fluxe
141. he surface air temperature from the common data area e lt OCETMOD gt extracts the sea surface temperature from the common data area e lt OCESMOD gt extracts the sea surface salinity from the common data area e lt AHUMMOD gt extracts the surface air humidity from the common data area e lt ACLDMOD gt extracts the atmospheric cloudiness from the common data area e lt AUVAMOD gt extracts the scalar surface wind from the common data area lt AUVSMOD gt extracts the standard deviation of the scalar surface wind from the common data area lt ASURMOD gt extracts the atmospheric surface pressure from the common data area 4 5 7 Coupling with Tracers lt TWRIT gt generates time averaged data required by the offline version of the tracer model 4 5 8 Routines for Vertical Interpolation lt TRASXY gt computes the horizontal section of some scalar quantity at some depth lt TRASZX gt computes the zx section of some scalar quantity at some depth lt TRASZY gt computes the zy section of some scalar quantity at some depth lt TRAVXY gt computes the horizontal section of some vector quantity at some depth lt TRAVZX gt computes the zx section of some vector quantity at some depth lt TRAVZY gt computes the zy section of some vector quantity at some depth lt TRAVONE gt computes a vertical profile of a vector quantity at one xy grid point Entry lt TRASONE gt computes a vertical profile of a scalar quantity at one xy gri
142. her from e g an atmospheric model or from data fed through the coupling interface Unfortunately the method is not as easy to apply than it reads The difficulty is that during the second step the updating process for the bias correction may produce oscillations such that the bias correction does not converge monotonously to some final field In order to save the situation it has been found that an EOF analysis of the bias correction fields yields PCs for the mean trend which decay continuously and further PCs that contain the oscillations If the trend EOF is used to reconstruct the time series for the bias correction then the oscillations are cancelled In a final step the arrays HEATM PMEM and TAUXM TAUYM are overwritten in the file OCDAT before the production experiment starts Alternatively one may use an exponential fit of the time series for the bias corrections in order to determine the stationary corrections in the limit case of an infinite long spinup run This method was used for the successfull experiments by Bacher et al 1998 5 4 8 About Numerical Stability of PIPE The CFL number does not provide a necessary condition for numerical stability as PIPE is a highly nonlinear model Furthermore e g stability properties of implicit schemes are based on the assumption that a solution of any either iteratively or directly solved equation is found exactly In practice this is never the case Any iterative algorithm leaves a r
143. ies e lt UVPOLEN gt Similar to above but used for full model arrays Entry lt SPOLEN gt is used for scalar quantities e lt BNDFLGS gt performs the boundary conditions for a land sea mask on scalar points e lt BNDFLGV gt performs the boundary conditions for a land sea mask on vector points e lt OBH_MOD gt overwrites boundary data of layer thickness in case open boundaries are used e lt OBT_MOD gt overwrites boundary data of potential temperature in case open bound aries are used e lt OBS_MOD gt overwrites boundary data of salinity in case open boundaries are used e lt OBZH_MOD gt overwrites boundary data of tidal sea level anomalies in case open boundaries are used e lt SAVOPEN gt saves open boundary data for later use e lt SEALAND gt sets vector quantities to zero on land 4 4 8 Routines for Listing e lt LISTING gt generates parts of the protocol that appears on OUTLIST Depending on the switches various fields are printed out mainly to provide a simple preview system e lt DRUCKV gt These routines print vector fields on the protocol as integer arrays with land set to blank fields If the model dimensions are too large for the width of the output only grid points with some increment are printed Entry lt DRUCKH gt is used for scalar quantities e lt EXTRAV gt Same as before but able to print all grid points within a portion of the model domain Entry lt EXTRAH gt is used for scalar qua
144. ine an artificial forcing when SW1 0 In this case no atmospheric data are used by PIPE e lt ABSOL gt reads one month of the COADS surface wind speed from UVABS and interpolates 1t onto the model grid e lt VARIAN gt reads one month of the COADS standard deviation of the absolute wind speed from UVDEV and interpolates it onto the model grid e lt OCTEMP gt reads one month of the COADS sea surface temperature from SSTEMP and interpolates it onto the model grid e lt ATTEMP gt reads one month of the COADS air temperature from AIRGLOB and interpolates it onto the model grid e lt CUMULUS gt reads one month of the COADS cloud cover from COVER and inter polates it onto the model grid e lt VAPOR gt reads one month of the COADS surface air relative humidity from WET NESS and interpolates it onto the model grid e lt ECUSC gt reads one month of the ECMWF absolute wind speed from USC and iter polates it onto the model grid e lt ECSTD gt reads one month of the ECMWF standard deviation of the wind from STDV_USC and interpolates it onto the model grid e lt ECTAU gt reads one month of the ECMWF surface wind stress from TAUXY and interpolates it onto the model grid e lt ECTSC gt reads one month of the ECMWF surface air temperature from TSC and interpolates it onto the model grid e lt ECSST gt reads one month of the ECMWF surface temperature from SST and inter polates it onto the model grid
145. ing bias correction is not exactly that which is expected when Newtonian feedbacks are finally switched off 2 learning what the bias correction is for zero feedbacks Setting SW48 2 reduces the Newtonian feedback continuously until it is negligible 3 freezing the bias correction and starting the production experiment During the first and second updating process at each time step the following recursion is used to improve the bias correction Qin Gran 2 Qi t 5 9 FH Fla D F 5 10 gt lA ll 5 11 fu aa ea 5 12 After each completed year the variables Qpias Frias Tx bias AN Ty pias are set to their respective current annual means Qrun Fran Tz run and Tyrun From the correction terms the final fluxes are computed by Qi Q Qbias 1 run T Tots Or 5 run 5 13 Fi FY Fhias 1 run S Sobs 85 run 5 14 Ta ET Bias 5 15 Ty Ty Ty bias 5 16 where dun is a variable see lt SURFLUX gt that decays from unity to zero during the updating process T and S are the model s surface temperature and salinity with Toss and Soss their respective observed fields r and s are feedback coefficients that induce a relaxation towards the observed fields at least in the annual mean Qj is the preliminary heat flux computed from radiative and turbulent fluxes F is the preliminary fresh water flux as computed through relaxation and P E R and 7 and 7 are wind stresses eit
146. ion introduced by the Crowley scheme dominates the physics for CFL gt 1 It can be understood as a quasi Lagrangian scheme in which higher order terms are neglected in the Taylor expansion of the transport equation At zz Peh sa pub yiray pub ay 1 uvr 7 1 WUD a 2 27 R Os 141 3 s 1 3 4Azy 1 At ouh yr1 1 1 pu yaa A Uy r1 1 1 UNED ea TV 1 1 A OSI Os 1 1 3 4Azy 1 1 a Opuh puh puh ha 2 2 28 U U Ti U SC es u gt y ay PUN T 7 P V I 1 J pun y 1 J V I 1 J V 1 J 2AY YN 4 OS 1 3 1 Zz OSJ cos P y y 2AYy 3 cos p y 3 1 cos p y y At puh yir 1 7 1 lead y A eya 1 Ur 1 gt _ 2AYy 1 1 OSI Os 1 1 1 cos y 7 1 2Ayy s 1 cos p y J 1 cos p y y 2 1 2 Discretization in Time 2 1 2 1 Predictor Step In order to allow larger time steps a predictor corrector technique is adopted Each of these steps is based on a semi implicit method For details see Robert et al 1972 The method yields an unconditionally stable time integration scheme with respect to all external and internal waves advection diffusion and mixed layer physics Coupling of the semi implicit steps has been described in Oberhuber 1986 The predictor step can be written as A ER At SN 7 At Aan Han RDU Say a GY ae Gyo 2 29 n At n n At XN KN ue 2 ze on ph gt ph At Fh 2 30 ph8 G PhS gt Fons pho
147. is the mixed layer All terms in which 1 or N occur are set to zero except for the term 7 7 which represents the surface wind stress and the term Fj which represents the bottom stress The forcing function Q represents the heat flux R the fresh water forcing and R the forcing due to the sea ice ocean coupling cp is the specific heat capacity of water The parameter f in Eqn 1 1 is the Coriolis parameter defined as 27 20 si it 1 f 20sin p with Q 36164 1 5 where 2 is the angular velocity of the rotating earth The stress 7 Tz Ty between neighbouring layers is set to zero Note that vertical diffusion causes stress between layers due to vertical mixing of momentum Since observed surface stresses are used the surface drag coefficient does not need to be specified however because vector wind is not provided through observed data however required for a few operation stress is converted to wind and vice versa assuming some surface drag coefficient Bottom stress is parameterized also using some drag coefficient In order to complete the equations the in situ density p and the potential density dy are related with temperature 7 salinity S and pressure P through the equation of state for sea water UNESCO 1981 written symbolically as Pk P Sp Px and O6k ple Sk Pref 1 6 where P is the reference pressure After discretization the in situ pressure in the first layer is given as a vertical mean b
148. ities do not vary much in time Therefore routines are offered that average a quantity and write out the average at a chosen rate For the subsequent routines it is essential that they are called at each time step as they use some internal bookkeeping The variable DTMONTH is used to determine the averaging period The default setting is 1 month which means that the subsequently described routines compute and write monthly means onto the history file e lt QENTXY gt writes out the time averaged entrainment and detrainment e lt QFLUXY gt writes out the time averaged heat flux fresh water flux and wind stress components e lt QICEXY gt writes out time averaged ice flow components ice thickness and ice com pactness e lt QMLDXY gt writes out the time averaged mixed layer depth onto the history file e lt QSSTXY gt writes out the time averaged sea surface temperature onto the history file e lt QSALTXY gt writes out the time averaged sea surface salinity onto the history file e lt QSLDXY gt writes out the time averaged sea level onto the history file e lt QUVXY gt writes out the time averaged sea surface velocity components onto the his tory file e lt QCONVXY gt bookkeeps and writes out maximal convection depth e lt QEBM gt writes out a number of quantities characterizing the EBM s physics e lt QTIDES gt writes out the anomalous sea level due to tides e lt QFTTID gt performs a modal analysis and writes o
149. ive for graphics a few tiny helpfull codes data for initializing PIPE usefull scripts There exist 3 tested installations for a 1 Parallel Vector Processor machine like a CRAY C90 2 Massively Parallel Processor machine like CRAY T3D 3 RISC Processor based workstations like from DEC HP IBM or SUN The code also works on a PC under Linux and has been successfully tested excluding scripts and postprocess ing on a PC using DOS6 22 and the GCC G77 compilers from the Gnu Project 101 5 1 1 Installation Step No 1 In order to arrive at a ready to compile model there exist the e Script Source setup pvp to create a PVP version on the directory Source PVP e Script Source setup rsc to create a RISC version on the directory Source RISC e Script Source setup mpp to create a MPP version on the directory Source MPP e Script Source setup bat to create a DOS version on the directory Source DOS The model source code contains preprocessing directives that are used by the scripts setup pvp setup rsc and setup bat These extract a PVP or a RISC version of the code The result is copied to Source PVP and Source RISC respectively The MPP version however ist not preprocessed by setup mpp It is only copied to the directory Source MPP The MPP codes are preprocessed at compile time later The installation of a DOS version is not complete All relevant scripts are missing However the
150. ivly in situ density from potential temperature and salinity for an array with dimension NZ N Y lt GRDX gt computes the pressure gradient in x direction lt GRDY gt computes the pressure gradient in y direction 4 7 4 Routines for Preparing Data lt UPDATEOBND gt organizes the updating of the boundary data in time lt UPDATEX gt updates an array with dimension NX lt UPDATEY gt updates an array with dimension NY lt UPDATEXZ gt updates an array with dimension NZ NX lt UPDATEYZ gt updates an array with dimension NZ NY lt STRLAYCUT gt organizes the computation of layer thicknesses and vertically inte grated temperatures and salinities from the profiles of temperature and salinity in z coordinates and finally computes sea level from the baroclinic and barotropic contribu tions lt LAYX gt computes the layer thickness and vertically integrated temperatures and salin ities for arrays with dimension NZ NX lt LAYY gt computes the layer thickness and vertically integrated temperatures and salin ities for arrays with dimension NZ NY lt WRITEC1X gt writes southern boundary data of a 3 dimensional array onto a 2 dimensional array containing an xz section lt WRITECNX gt writes northern boundary data of a 3 dimensional array onto a 2 dimensional array containing an xz section lt WRITEC1Y gt writes western boundary data of a 3 dimensional array onto a 2 dimensional array containing an yz
151. l using Eu lerian angles ALPHA 40 BETA 50 GAMMA 0 Geographical longitude and latitude are plotted for better orientation The picture is created interac tively using plot b p e q 69 PPSF_SUR 24 eee i oe ee we es 144 Summary Page Short description of the model The name PIPE is derived from Parallel Isopycnal Primitive Equation model Its origin is the OPYC model It has been renamed because it now is fully portable between parallel shared and distributed memory computers as well as between vector and scalar ones In addition it experienced significant changes in numerics and physical complexity The concept of using isopycnals as the vertical coordinate system for an OGCM is based on the observation that the interior ocean behaves rather like a conservative fluid Even over long distances the origin of water masses can be traced back by considering the distribution of active or passive tracers Isopycnal coordinates have the advantage that diapycnal mixing by numerical discretization errors is explicitly excluded which may cause a too diffusive thermocline in conventional z coordinate models that don t use advanced correction techniques Treating the ocean as a conservative fluid fails in areas of significant turbulence activity such as the surface boundary layer Therefore a surface mixed layer is coupled to the isopycnic interior ocean in order to represent near surface vertical mixing and to improve the response time scales to at
152. lated on land as well In the present version however calving of glaciers or glacier flow is not yet considered This could be included easily based on the same idea as applied for the river runoff model however with some diagnostic equation for the ice flow which replaces that for the water flow 1 6 6 Soil Model Equations The equation for the soil s temperature Toii 18 To ebm Teil Qem 1 115 Ot Chis fells Cop Pw Webm where cp is the soil s heat capacity p the soil s density and H the soil thickness Thus the soil includes the amount of water heated or cooled however the heat transport of water is ignored 1 6 7 Tuning of the EBM All parameters g are tuned in such a way that annual mean temperatures fit well with observations The tuning strategy was to 1 adjust the parameters for short and longwave radiation to ERBE data and ECHAM4 AMIP experiments 2 tune the annual zonal mean surface temperature to the EC MWE analysis by choosing a proper eddy diffusivity that exports enough heat from low into high latitudes 3 tune the effective height used for moisture in order to obtain a realistic vertically inte grated water content 4 tune convective precipitation so that enough water is left to be exported into higher latitudes 5 tune adiabatic heating in such a way that the sum of latent heat release due to conden sation and adiabatic heating due to vertical motion yields a spatially smooth field B
153. lation Model Part I Model Experiment J Phys Oceanogr Vol 23 No 5 830 845 Oberhuber J M 1993c The OPYC Ocean General Circulation Model Technical Report No 7 Deutsches Klimarechenzentrum GmbH Hamburg Paulson C A J J Simpson 1977 Irradiance measurements in the upper ocean J Phys Oceanogr 7 952 956 Robert A J Henderson and C Turnbell 1972 An implicit time step scheme for baroclinic models of the atmosphere Mon Wea Rev 100 329 335 Schwiderski E W 1979 Ocean Tides Part I Global Tidal Equations Marine Geodesy 3 161 179 Shea D J 1986 Climatological Atlas 1950 1979 National Center for Atmospheric Re search Boulder Colorado Smagorinsky J S 1963 General circulation experiments with the primitive equations I The basic experiment Mon Wea Rev 91 99 164 Smolarkiewicz P K 1982 The multi dimensional Crowley Advection Scheme Mon Wea Rev 110 1968 1983 UNESCO 1981 The Practical Salinity Scale 1978 and the International Equation of State of Seawater 1980 Unesco Techn Pap in Mar Sci No 36 13 21 Woodruff S D R J Slutz R L Jenne and P M Steurer 1987 A Comprehensive Ocean Atmosphere Data Set Bull Amer Met Soc 68 1239 1250 Wright P 1988 An Atlas based on the COADS data set Fields of mean wind cloudiness and humidity at the surface of the global ocean Max Planck Institute for Meteorology Hamburg Report No 14 70pp Zillmann
154. le e Script anomaly computes the anomaly of subsequent data of the same type It is a macro that calls other scripts e Script stdv computes the standard deviation of subsequent data of the same type 5 5 6 4 Time Series Analysis e Script xtime extracts data for a Hovmueller diagram x direction and time from an input file that contains a sequence of some quantity ordered in time The location of the section is defined by one parameter e Script ytime extracts data for a Hovmueller diagram y direction and time from an input file that contains a sequence of some quantity ordered in time The location of the section is defined by one parameter e Script regress computes the linear regression coefficient from subsequent data of the same type e Script degress reconstructs a time series of subsequent data of the same type by using the linear regression coefficients that have been generated by regress e Script ppsf2eof computes the empirical orthogonal functions EOFs and the principal components PCs from subsequent data of the same type e Script eof2ppsf reconstructs the time series from empirical orthogonal functions EOFs and according principal components PCs e Script eof2cca performs a canonical correlation analysis CCA using EOFs and PCs as generated by ppsf2eof Warning The dimension statements must be adjusted to the particular settings of the current model version e Script ppsf2cca performs a cano
155. linity lt TRACER gt ACCURF 0 00001 required accuracy for tracer lt TRACER gt Table 4 24 Numerical Tuning Variables Variable Value Unit Meaning Reference STABMIN 0 01 kgm threshold for density difference lt STABIL gt STABDIF 0 005 ms threshold for g in mixed layer lt MIXEXP gt HMIXMIN 5 m threshold for mixed layer depth lt MIXEXP gt UVMAX 5 ms 7 threshold for velocity lt LIMVEL gt TOPOMIN 20 m threshold for ocean depth lt OCEINIT gt Table 4 25 Threshold Variables 4 3 The Ocean Model ocestep f 4 3 1 Interior Ocean Model Subroutines e lt OCESTEP gt is the driving routine to carry out the time step At each call of it MAX REP time steps are carried out unless other thresholds as defined by the maximal CPU time CALCMAX or the final model time NYEAR etc are reached It can be called several times from the main program e lt CONVECT gt computes the exchange of momentum heat and salt concentration be tween two adjacent layers in case of unstable stratification In the first part the mixed layer exchanges information with the next adjacent layer and in the second part two lay ers in the deep ocean exchange their properties In the latter case the coordinates are rearranged if the coordinate error exceeds a certain threshold e lt CORFORC gt computes the Coriolis force and adds its contribution to the momentum change to the arrays RESTU and RESTV Entry lt RALEIGH gt adds on a Raleigh fri
156. lt STRDATNX gt writes northern boundary of observed ocean data on a 3 dimensional array onto a 2 dimensional array containing an xz section lt STRDATIY gt writes western boundary of observed ocean data on a 3 dimensional array onto a 2 dimensional array containing an yz section lt STRDATNY gt writes eastern boundary of observed ocean data on a 3 dimensional array onto a 2 dimensional array containing an yz section lt BAROTROP gt merges the baroclinic and the barotropic sea level anomaly and con structs a sea level which is connected around corners of the model domain This routine contains a stop statement in order to remember a user to adjust this routine to the user s application 4 7 2 Routines for Coordinate Conversion lt CONVDX gt converts the definition of the vertical coordinate s origin for an array with dimension NX lt CONVDY gt converts the definition of the vertical coordinate s origin for an array with dimension NY lt LAYDEPTX lt performs the data conversion from level to layer coordinates i e com putes layer thicknesses from ocean data and the vertical integrals for temperature and salinity for an xz section lt LAYDEPTY lt performs the data conversion from level to layer coordinates i e com putes layer thicknesses from ocean data and the vertical integrals for temperature and salinity for an yz section lt INTEGRAX lt performs the vertical integration of potential density potential tem
157. m COADS on a 2 x 2 grid This file is not used e File UVGLOB It contains the monthly mean surface wind stress and its annual mean from Hellermann and Rosenstein as pseudo wind stresses on a 2 x 2 grid The first block of data is the annual mean the second one starts with January It is read by lt WIND gt with the following simplified statements DIMENSION ITAUX 180 ITAUY 180 TAUX 180 72 13 TAUY 180 72 13 DO MONTH 1 13 DO J 1 90 READ 2014 CITAUX I 1 1 180 READ 2014 CITAUY I 1 1 180 DO I 1 180 TAUXM ITAUX 1 5000 5000 TAUYM ITAUY 1 5000 5000 TAUX 1 J MONTH TAUXM SQRT TAUXM 2 TAUYMx x 2 TAUY 1 J MONTH TAUYM SQRT TAUXM 2 TAUYM 2 ENDDO ENDDO ENDDO e File UVABS contains the monthly mean wind speed from the COADS on a 2 x 2 grid It is read by lt ABSOL gt with the following simplified statements DIMENSION IABSOL 180 UVABSOL 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 IABSOL I I 1 180 DO I 1 180 UVABSOL 1 J MONTH IABSOL 1 10 ENDDO ENDDO ENDDO e File UVDEV contains the monthly mean wind speed standard deviation from the COADS on a 2 x 2 grid It is read by lt VARIAN gt with the following simplified state ments DIMENSION IVARIAN 180 UVARIAN 180 90 12 DO MONTH 1 12 DO J 1 90 READ 2014 CIVARIAN I I 1 180 DO I 1 180 UVARIAN 1 J MONTH IVARIAN I 100 ENDDO ENDDO ENDDO 5 4 3 1 3 ECMWF Data
158. mass transport stream function the vertical overturning stream function and many other diagnostic variables such as vertical velocity potential vorticity etc 5 5 4 1 4 File PPSF_SEC contains the flow temperature salinity and vertical velocity on those depth levels which are defined by VSEC see table 4 13 5 5 4 1 5 File PPSF_VER contains xz and yz sections of temperature salinity zonal and meridional velocity components and the interface distribution on longitudes and latitudes that are defined by ISEC and JSEC as E and Jindices of the model arrays 5 5 4 1 6 File PPSF_FLX contains the fluxes into the ocean These are the surface wind stress the friction velocity and the fluxes of heat fresh water and buoyancy Dependent on ISW4S8 it also can contain the long term mean heat and fresh water flux or the bias correction terms if enabled It also contains the EBM fluxes 5 5 4 1 7 File PPSF_FOR contains the model topography the absolute wind speed and its standard deviation the air temperature the sea surface temperature and the sea surface salinity at the end of the model run and some EBM output 5 5 4 1 8 File PPSF_TID contains quick look data of the tide model 5 5 4 1 9 File PPSF_TOP contains the model topography if the ocean state has been initialized Note that this file should be used to verify a correct initialization of the grid and coastline geometry Code Number Type of Section Reference intergral over x
159. mensions and other options as documented in setup 2 Run Source RISC Model setup 3 Adjust code essentially ocemain f to obtain a model domain and proper choice of the model parameters 4 Run Source RISC Model make sun in order to create object files which by default are copied to the subdirectory Source RISC Model Run 5 Adjust Source RISC Model model sun to perform a model initialization 6 Run Source RISC Model model sun 7 Adjust Source RISC Model model sun to perform a continuation run 8 Run Source RISC Model model sun In order to setup an arbitrary version for a RISC computer the file Source parampvp h is inserted into the code through an include statement The dimensions nz_tot and nz for the x direction and the dimensions ny tot and ny for the y direction each have identical meaning These dimensions are set through the script Source RISC Model setup There is no need to modify parampvp h during the installation 5 1 2 2 Installation of a PVP version 1 Adjust Source PVP Model setup in order to obtain FORTRAN codes that use the proper dimensions and other options as documented in setup 2 Run Source PVP Model setup 3 Adjust code essentially ocemain f to obtain a model domain and proper choice of the model parameters 4 Run Source PVP Model make cri in order to create object files which by default are copied to the subdirectory Source PVP
160. mine the Ocean Flow J Geomag Geoelectr Vol 49 1351 1372 Nies H ILH Harms M J Karcher D Dethleff C Bahe G Kuhlmann J M Oberhuber J O Backhaus E Kleine P L we D Matishov A Stephanov and O F Vasiliev 1998 Anthropogenic radiaoactivity in the Nordic Seas and the Arctic Ocean Results of a joint project Project Report BMBF Christoph M T B Barnett and E Roeckner 1998 The Antarctic Circumpolar Wave in a Coupled Ocean Atmosphere GCM J Climate Vol 11 1659 1672 Cabos Narvaez W M J OrtizBevi and J M Oberhuber 1998 The variability of the tropical Atlantic J Geophys Res Vol 103 No C4 7475 7489 Bacher A J M Oberhuber and E Roeckner 1998 ENSO dynamics and seasonal cy cle in the tropical Pacific as simulated by the ECHAM4 OPYC3 coupled general circu lation model Clim Dyn Vol 14 432 450 Zhang X H J M Oberhuber A Bacher and E Roeckner 1998 Interpretation of interbasin exchange in an isopycnal ocean model Clim Dyn Vol 14 725 740 Timmermann A J M Oberhuber A Bacher M Esch M Latif and E Roeck ner 1998 ENSO Response to Greenhouse Warming Max Planck Institut f r Mete orologie Hamburg Report No 251 Nature in press Oberhuber J M E Roeckner M Christoph M Esch and M Latif 1998 Predicting the 97 El Nino event with a global climate model Geophys Res Let Vol 25 No 13 2273 2276
161. mined since one is free to choose the amount of water to be entrained from above or from below If a free parameter a for this unknown ratio is introduced the equation for the total entrainment rate w is given by wy 1 a wh awe 1 72 with individual entrainment rates we 1 ajw 1 73 wet awit 1 74 After substituting these entrainment rates into the momentum mass heat and salt equations the free parameter a can now be chosen to maintain the potential density within some layer at or nearby a prescribed value oj while at the same time compensating some potential density drift due to discretization errors in the advection and diffusion formulation In order to balance these errors must be made space and time dependent The relevant balance equation is given by the approximate equation for the potential density at the new time level hy At wy wi hpoop Ati oon AWET Oop 1 75 where the potential density at the new time level should be equal the prescribed value oj This equation was derived under the assumption that the potential density of a mixed water mass is the layer thickness weighted average of the potential densities of the unmixed water masses This yields the final equation for a _ Paloi oon Atow 05 on 1 76 Atwj 05 Cor At 05 Font 1 4 2 Convection If the stratification becomes unstable it is removed by vertical mixing All quantities are set to their verti
162. mospheric forcing which is controlled by the mixed layer thickness In addi tion parameterizations for vertical mixing in the interior ocean and a convection mechanism completes the model Since the model is designed for climate studies on large scales a snow and sea ice model with rheology is included and serves the purpose of decoupling the ocean from extreme high latitude winter conditions and promoting a realistic treatment of the salinity forcing due to melting or freezing of sea ice An atmospheric energy balance model EBM is forced by atmospheric winds and cloudiness and predicts some vertical mean temperature and humidity Horizontal eddy diffusivity is parameterized Further contributions to the heat budgets as radiative fluxes diabatic and adiabatic heating and turbulent surface fluxes are diagnostic variables In order to close the global heat budget evaporation and precipitation are linked with the sensible heat budget Simple models of river runoff and glaciers on land close the global heat and water budget A barotropic tide model coupled to the 3 dimensional ocean via bottom drag and residual currents improves the circulation in regional application The open boundary formulation allows to feed the model with temperature and salinity from 3 dimensional ocean data By assuming some barotropic mode along the open boundaries multiple disconnected boundaries can be used to compute the ocean state in a fraction of an ocean basin Auth
163. ms to add channels or land bridges across topography and or coastline geometry They are disabled if set to zero see table 4 15 4 2 2 3 2 Block Data ROTATE defines the three Eulerian angles which are responsible for a rotation in longitude then in latitude and finally in longitude The rotation becomes active if ISWITCH 1 see table 4 12 4 2 2 4 Block Data which concern Data Postprocessing 4 2 2 4 1 Block Data POSTPRO defines parameters that allow one to control which vertical or horizontal sections appear on the quick look files Up to 10 different horizontal or vertical sections can be selected see table 4 16 Variable Value Unit Meaning Reference TMELT 2713 16 K freezing temperature of fresh water lt TFREEZ gt CPMELT 334000 Wskg heat of freezing lt SEAICE gt CPICE 2090 Wskg K specific heat capacity of ice lt SEAICE gt SALTICE 5 0 gkg salinity of sea ice lt SEAICE gt CMELT 1 0 tuning coefficient for melting lt SEAICE gt CFREEZ 2 0 tuning coefficient for freezing lt SEAICE gt EXCENT 2 0 excentricity of ice rheology lt SEAICE gt PRESICE 25 0 Wskg proportionality for ice pressure lt SEAICE gt PDECAY 2 0 decay coefficient for ice pressure lt SEAICE gt HNULL 0 5 m tuning coefficient for freezing lt SEAICE gt EPSNULL 1 E 6 ms coefficient that limits ice pressure lt SEAICE gt AGING m1 E 7 s7 snow to ice conversion rate lt SEAICE gt CICE 2 0 WK im heat conductivity for ice lt SURFLUX gt C
164. nal ocean sampling problems are avoided through further averaging in time in order to obtain smooth residual currents Tides are forced through a gravitational anomaly E following Schwiderski 1979 If ho sy and po are the mean longitudes of sun moon and lunar perigee at Greenwich midnight respectively then these are defined through T 27392 500528 1 0000000356 D 1 119 ho 279 696680 36000 768930485 T 0 000303 T 1 120 So 270 434358 481267 88314137 T 0 001133 T 0 0000019 T 1 121 Po 334 329653 4069 0340329575 T 0 010325 T 0 0000120 T 1 122 where D is the number of days starting with D 0 at 01 01 1975 This results in an equation for the gravitational anomaly induced by the sun and moon as equivalent sea level anomaly 4 Pyae 0 69 cos y Y Kn cos ont 2 Xn n 1 8 sin 2y Y Kn cos ont Xn n 5 3cos y 2 y Kn cos ont Xn 1 123 n 9 where p is the latitude A the longitude and t the time in 10 s The partial tide s amplitudes K frequencies o and astronomical arguments x are defined in table 1 1 Tidal Mode K o X M Principal Lunar 0 242334 1 40519 2 ho so S Principal Solar 0 112841 1 45444 0 No Elliptical Lunar 0 046398 1 37880 2ho 389 Po K Declination Luni Solar 0 030704 1 45842 2ho K Declination Luni Solar 0 141565 0 72921 ho 90 O Principal Lunar 0 100574 0 67598 ho 289 90 P Principal Solar 0 046
165. nd PIPE 4 2 2 1 Block Data which concern Flow Control 4 2 2 1 1 Block Data SCREW defines the time step and its limits The maximal number of time steps carried out during one call of lt OCESTEP gt is MAXREP unless CALCMAX initiates an earlier exit see also table 4 1 in which CPUs denotes seconds CPU time and MODs denotes seconds of model time Switches control the coupling interface NHEAT NFRESH NSTRS NMIX and NICE are active only if NINIT 1 Depending on how the switches NHEAT through NICE are set the fluxes are either interpreted as total fields or as anomalies see table 4 2 Also there exists the possibility to feed the model with atmospheric raw data which then are used to calculate fluxes in the same style as performed with the data available on the file FORCES Finally the variables KHEAT through LENAVER determine whether a bias correction is applied see table 4 3 Variable Value Unit Meaning Reference MAXREP any gt 0 maximal number of time steps lt OCESTEP gt CALCMAX any gt 0 CPUs CPU time after which model stops lt OCESTEP gt DTMIN any lt DTMAX MODs smallest allowed time step lt OCESTEP gt DTMAX any gt 0 MODs largest allowed time step lt OCESTEP gt DTSTOP any lt DTMAX MODs time step at which model exits lt OCESTEP gt NHOUR any gt 0 MODh _ completed hours for model stop lt OCESTEP gt NDAY any gt 0 MODd completed days for model stop lt OCESTEP gt NMONTH any gt 0 MODm completed months for model stop l
166. neration of Isopycnal Coordinates MAKENEW gt organizes the generation of the model s isopycnal coordinates MAKEVER gt computes the vertical coordinate system MAKETET gt computes potential density for a single profile lt lt lt lt MAKETET gt computes in situ density for a single profile lt OLDINIT gt is the driver for generating a 3 dimensional ocean state for initialization from data imported from some OCDAT however whose data use a different horizontal and vertical grid 1 e number of layers NZ lt REFORMH gt interpolates scalar model quantities between two different grids Entry lt REFORMU gt interpolates vector quantities between two different grids 4 4 6 Conservation Properties lt MASSCON gt computes the mass content of the ocean lt HEATCON gt computes the heat content of the ocean lt SALTCON gt computes the salt content of the ocean lt HGTCOR gt removes residual negative layer thicknesses lt MEAN2D gt computes the 2 dimensional area mean lt MEAN3D gt computes the 3 dimensional area mean lt ADJUST gt optionally corrects for the mass heat and salt content when one of the budgets is not closed lt ZEROLAY gt fills values in cells that have physically closed and have zero thickness This routine is only needed for cosmetics to provide some numbers for a zero layer 4 4 7 Boundary Conditions The basic idea behind the model s structure is that boundary conditions never appear e
167. ng warm periods melting ice decreases the salinity in the mixed layer and therefore contributes to a stabilization of the upper ocean Ba sically sea ice has the task to convert an atmospheric heat flux into an oceanic fresh water flux While some heat flux results in a vanashing buoyancy flux due to the low thermal expansion of sea water around the freezing point a fresh water flux induces a high buoyancy flux as result of brine rejection Snow cover modifies the heat fluxes that occur in the presence of an ice cover This is because snow has both a lower conductivity and a higher albedo than ice Furthermore the albedo depends on the snow type and for a thin snow cover also on its thickness 1 5 1 The Dynamic Equations Hibler 1979 proposed a rheology for a dynamical sea ice model that can be used for a wide range of space and time scales For a number of technical reasons such as the introduction of spherical coordinates and of the momentum and mass conserving flux form of the equations which permits an easier treatment of the ice edge the model has been completely rewritten from scratch but based on the same physics as in Hibler s model A snow model has been added based on a continuity equation for snow and a parameterization for the aging process of snow The heat capacity of snow and ice is included via two prognostic temperatures for the skin temperature over snow and ice The basic equations for the cell averages of the ice flux
168. nical correlation analysis CCA from subsequent data of the same type It is a macro that calls other scripts Warning The dimension statements must be adjusted to the particular settings of the current model version e Script detrend detrend subsequent data of the same type using an EOF analysis It is a macro that calls other scripts Warning The dimension statements must be adjusted to the particular settings of the current model version 5 5 6 5 Frequency Analysis e Script spectrum performs a spectral analysis and outputs data either as curves or as Hovmueller diagram in time and either in x or in y direction e Script fourier performs a fourier analysis and outputs amplitude and phases at a spec ified period e Script wavelet performs a wavelet analysis e Script fttid performs a fourier analysis for tidal frequencies Warning The dimension statements must be adjusted to the particular settings of the current model version 5 5 6 6 Data Filtering Script smooth filters the data using a Helmholtz equation Thus the result is widely independent on the actual grid resolution Script rmean computes the running mean over a specified period 5 5 6 7 Data Reformatting Script spr2dpr converts a binary PPSF format that has been created with single pre cision into a binary PPSF format that should contain double precision words This pro cedure is useful if the model runs on a computer with 32bit CPU with double pr
169. nonlinear equations The disadvantageous consequence is that the matrix coefficients need to be updated during each iteration step in order to find the solution of the nonlinear equations However ocean models are not made to spend CPU time only but should produce good oceanographic results PIPE has demonstrated it by a number of key publications see reference list The full source code of PIPE contains code sections that are either commonly or alterna tively used by PVP MPP or RISC versions of the model These versions are created using preprocessing directives System commands like cpp or fpp extract the model version from the full source code The strategy to achieve a high Mflop rates on all these architectures is as follows e 3 d arrays are always dimensioned as NZ NX NY where NZ is the vertical index NX is the zonal index and NY the meridional index Since there are no vertical derivatives to compute the 1st index always runs over its full dimension The 2nd index typically runs between 2 and NX 1 which allows to write many statements with vector length NZ NX 2 If the 2nd index runs over its full dimension then the 3rd index allows even longer vectors for instance of length NZ NX NY e If a 2 d operation is carried out with a 3 d array these operations always run over NX and NY This means that vector elements are not contiguous but have an increment of NZ Therefore on a CRAY NZ must be odd otherwise bank conflicts occur
170. ns by lt SURFLUX gt The routine subtracts from QFLUX that part which is used for melting or freezing The same applies for the salt flux SFLUX lt DIADICE gt computes the thermodynamic forcing of the sea ice model lt ENSURE gt makes ice extent consistent with SST distribution lt PFREEZ gt computes the dependence of the freezing point of sea water on salinity 4 3 5 Atmospheric Energy Balance Model lt EBM gt computes the horizontal transport of temperature and moisture It also predicts the soil temperature over land and below a glacier if existing The radiation calculation and that of the turbulent fluxes is performed in lt SURFLUX gt lt RIVERS gt performs the river runoff on land 4 3 6 Tide Model lt POTTIDE gt is the driver for the tide model It organizes the timestepping 1 e the subcycling of the tide model It also bookkeeps the residual currents and performs the data postprocessing lt TIDES gt predicts the anomalous tidal flow and sea surface elevation lt EQUIPOT gt computes the gravitational forcing of the tide model 4 3 7 The Equation of State lt EXPANDT gt computes the expansion coefficients with respect to temperature The UNESCO formula is differentiated analytically Entry lt EXPANDS gt performs the same calculation for salinity Entry lt EXPANDP gt performs the same calculation for in situ pressure Entry lt EXPANDT gt performs the same calculation for potential temperature
171. ntire model It performs one complete iteration of the wave equation and yields the layer thickness as the result It uses a line relaxation method that solves a linear equation system directly on yz slices and iterates in the x direction lt STRESS gt computes the contribution of the surface stress the stresses between the layers and the bottom stress to the right hand side of the momentum equation which is added to RESTU and RESTV lt VTOTAU gt converts the vector surface wind to stress Entry lt TAUTOV gt performs the inverse operation lt TRHTOU gt computes mean layer thickness on velocity grid points with additional weights The output variable HQLIN is a key variable needed for the flux form of the dynamical equations lt STABIL gt limits the stability to a marginally positive value lt HTOZETA gt converts layer thickness to interface height above reference level Entry lt ZETATOH gt performs the inverse operation 4 3 2 Mixed Layer Model Subroutines lt MIXEXP gt computes the mixed layer physics explicitly forward It computes the entrainment and detrainment rates diagnostically and delivers changes of layer thickness and momentum to lt OCESTEP gt Entry lt MIXIMP gt performs the same calculations It treats the mixed layer implicitly by an iterative procedure to find the new layer thickness for the corrector step lt EQUILIB gt computes the Monin Obukhov length by Newton s method to find zeros Values a
172. ntities e lt DRUCKF gt prints out the land sea mask 4 4 9 Routines for Model I O e lt WRITEM gt writes all the model fields onto the restart file OCDAT that are required to restart the model for continuing a model run Optionally the output can be saved as an ASCII file Entry lt READM gt reads the first part of OCDAT In case the file is an ASCII version of OCDAT lt READM gt uses formatted I O automatically Entry lt READMF gt reads the second part of the data from OCDAT that have been created by lt WRITEM gt If the file is an ASCII version of OCDAT lt READMF gt also uses formatted I O automatically lt FILOPEN gt opens various files Entry lt FILCLOS gt closes various files lt PPSFAST gt writes data onto quick look files The data format is called PPSF Post Processing System Format 4 4 10 Routines for Model and CPU Time lt COUNTER gt converts the variable KTOT AL into seconds hours days months and years lt TIMECON gt converts the variable KTOTAL into seconds minutes hours days months and years lt PRINCPU gt prints out the bookkeeped CPU times 4 4 11 Routines for Vertical Flag Calculations lt TOUCHLW gt computes the indices of the next deeper physical layer for each horizon tal grid box within a given layer It uses a critical number EPSHGT to detect a physical layer Entry lt TOUCHUP gt computes indices for next upper physical layer lt CONTLW gt works similarly
173. olic temperature distribution where TEMMEAN is used as surface temperature and 273 16 K is used as bottom temperature see table 4 17 The potential density coordinates are stored onto POTSOLL Note that DRAGS is used to convert the stress data into wind vectors ATU and ATV The same drag coefficient is used to convert the wind vectors back to stress values see table 4 18 The definition of TURBLEV allows to tune the upper ocean thermocline while HMIX is more responsible for the deep ocean see table 4 13 The mixed layer parameters are well tuned see table 4 9 4 2 2 3 Block Data which concern Initialization 4 2 2 3 1 Block Data LAYOUT defines the model boundaries and the parameters used to create a focus 1 e an area with enhanced resolution Furthermore the Eulerian angles allow one to rotate the model longitudes and latitudes on the sphere This is useful for rotating a Variable Value Unit Meaning Reference HEATLAT 2500800 Wskg heat of fusion lt SURFLUX gt SIGMA 5 67E 8 WsK Stefan Boltzmann constant lt SURFLUX gt SOLAR 1367 Wm solar constant lt SURFLUX gt GAS 287 m K157 gas constant of dry air lt SURFLUX gt CPAIR 1005 Wskg K 1 specific heat capacity of air lt SURFLUX gt SLPRESS 100000 Nm mean air surface pressure lt SURFLUX gt CONSTS 0 0327 constant for sensible transfer coefficient lt VARDRAG gt CONSTL 0 0346 constant for latent tranfer coefficient lt VARDRAG gt CHARNCK 0 016 Charnock constan
174. ompactness diffusion according to 2 1 1 7 and the mass and compactness convergence ac cording to 2 1 1 4 The discretization of the quadratic and mixed derivatives of the rheology terms is based on 9 point formulae Thus the derivatives of the velocity components are car ried out on scalar points according to 2 1 1 4 and the result is differentiated on vector points according to 2 1 1 3 Note that the 9 point formulae do not diffuse B grid computational waves which appear as checker board pattern However these spurious waves do not appear in the solutions anyway because the formulae are consistent with the divergence terms and thus permit an explicit Eulerforward step for the mass and compactness convergence 2 2 2 Discretization in Time A predictor corrector method in connection with the semi implicit technique is applied First diffusion of ice thickness and ice concentration are determined implicitly In this step the flux divergence of ice thickness and ice concentration are still taken explicitly forward In the next step the stress and Coriolis term are treated implicitly In the final correction step the sea ice rheology in flux form is determined implicitly The coupling of the predictor and two corrector steps has already been outlined in the contest of the ocean model All iterations are carried out in the y direction only since the system of linear equations is solved directly in x The matrix coefficients in the rheology part a
175. on wind stress bias correction Table 5 5 Definition of Quantity Code Flux Variables Code Number Variable Unit barotropic stream function meridional overturning overturning in Atlantic Ocean Reference makeall f makeall f makeall f makeall f makeall f overturning in Indic Ocea overturning in Pacific Ocean Table 5 6 Definition of Quantity Code Transport Variables MHOUR is the model time in hours after the full day MONTH is the model time in months after the full year MYEAR is the model time in years after the start LDATE contains the date at which the data set was created LTIME contains the time at which the data set was created EXPERIM It is a character string declared as CHARACTER 40 and contains the identification of an experiment X1 contains the coordinates of the first index of the array FIELD X2 contains the coordinates of the second index of the array FIELD FIELD contains the quantity to be saved Code Number Variable Unit Reference surface air relative humidity surface water on land atmospheric water content diagnostic precipitation surface air temperature top shortwave radiation bottom shortwave radiation top longwave radiation bottom longwave radiation sea surface and soil temperature sea surface temperature ice and snow cover horizontal heat transport vertical heat release orography surface turbulent fluxes Table 5 7 Definition of Qu
176. or make cri After having adjusted all control parameters in ocemain f the entire source code is ready for compilation In case any modification of some routine is necessary save it onto ocemods f 5 4 2 Start from Scratch There are several steps to develop a model version that can be restarted and used for production 1 Grid Check The purpose of the first run is to verify that the grid definition was carried out properly For that run set JSW11 1 in ocemain f when starting model sun or model cri for the first time Since the x and y coordinates as well as the topography are initialized first SW11 allows the model to stop after short a CPU time Before that it saves PPSF_TOP if SAVE POST 1 in model sun or model cri Ensure that TOPFORC in model sun or model cri is consistently defined with 5W38 The OUTLIST contains the x and y coordinates as well as the relative changes of the grid spacing With the help of this listing the model margins and the focus parameters can be optimized For this run also set AGEDATA 0 UPDATED 0 INITIAL 1 and SAV FORC 0 2 Grid Manipulation After having obtain a first guess for the land sea mask it is com monly necessary to modify the land sea mask The earlier still available mechanism is to use the tools which are accessible through switches defined in ocemaon f Meanwhile the switch JSW45 allows a better control of the land sea mask which is recommended 3 Model Initiali
177. or of the model Josef Maximilian Oberhuber Deutsches Klimarechenzentrum GmbH Person responsible for model support at the DKRZ Josef Maximilian Oberhuber DKRZ Model Support Group Ivalid till 31 July 1999 11 Changes relative to OP YC2 17 see DKKZ Rep No 7 Changes in Physics e The parameterizations for the mixed layer and internal mixing have been altered to avoid too deep mixed layers in the subtropics e The assignment of water masses to specific isopycnal layers has been reformulated for all vertical exchange processes such as for mixed layer internal mixing and convection e The radiation parameterizations are improved in order to obtain a more realistic climate sensitivity Additional Physic Modules e An atmospheric energy balance model EBM has been included which considers a closed budget of sensible and latent heat and contains models for soil river runoff and glaciers e A tide model is coupled via bottom drag and residual currents to the 3 dimension ocean model Changes in Numerics e Alternating direction line relaxation scheme for momentum and mass instead of line relaxation in x direction only for better convergence e Exact mass conservation provided that the wave equation has converged Additional Technical Tools e Open boundaries now allow to specify multiple disconnected vertical sections on which ocean data like temperature and salinity can be prescribed e A tracer model can be run in an online or of
178. or predicting the thermody namics of the snow amp sea ice model while lt SEAICE gt is responsible for the dynamical part lt SURFLUX gt also computes radiative and turbulent fluxes for the EBM e lt VARDRAG gt computes the drag coefficient and the transfer coefficients for sensible and latent heat from wind speed boundary layer stability and humidity 4 4 Preprocessing oceproc f 4 4 1 Driving Routines e lt OCEINIT gt interpolates global data sets onto the model grid and stores them in FORCES via lt WRITEF gt The file names used as data source from the following routines are the same names as those of the permanent files and not those used as local file names Furthermore it reads file OCDAT and initializes the 3 dimensional ocean state and it allows to initialize the ocean by interpolating data of an existing OCDAT onto data on a different grid e lt OCEPOST gt is the driving routine for postprocessing performed after the model run Entry lt OCEPST4 gt writes out data at the end of the model run This routine is called from lt OCESTOP gt e lt OCESTOP gt cleans up the model run finally 4 4 2 Grid Initialization e lt MAKEXY gt initializes the x and y coordinates of the model using either the individual definitions for the model boundary XGL XGR YGU and YGO or the Gaussian grid for the T21 T32 T42 T63 T84 or T106 resolutions To obtain one of the latter grids the preprocessor precomp f has to
179. ospheric sea level pressure HUMID relative humidity DEPTH topography depth Table 7 10 Names of Observed Quantities Name in Code Name in Formulae Meaning Unit Equation TEBM 1 welght for temperature transport HEBM E2 weight for moisture transport DIFWGHT 3 weight for eddy diffusion 4 welght for precipitation 5 weight for latent heat conversion 6 weight for large scale precipitation 7 weight for convective precipitation FLOWFAC g weight for river runoff velocity H thickness of the soil m HHUMY H effective thicknes of the moist layer m QINSOL GO actin top downward solar heat flux Wm QINLON O tsp top longwave radiative heat flux Wm Table 7 11 Names of EBM Parameters Name in Code Name in Formulae Meaning Unit Equation UHTIDE uh tide zonal barotropic tidal flow mes VHTIDE uh tide meridional barotropic tidal flow m s ZHTIDE Cride tidal sea level anomaly UTIDE Upes zonal residual flow due to tides VTIDE Urea meridional residual flow due to tides U NZ 1 J UN zonal flow in the lowest layer V NZ 1 J UN meridional flow in the lowest layer EQUTIDE tide gravitational tidal forcing Table 7 12 Names of Tide Parameters 7 2 Acknowledgements I would like to thank Trond Aukrust Bergen Scientific Centre Bergen for his contributions to the grid rotation by Eulerian angles David Holland McGill University Montreal for the bookkeeping of the budget terms and Frank Kauker GKSS Geesthacht for
180. otential vorticity and energy conserving scheme for momentum transport Bleck and Boudra 1981 is to rearrange the momentum advection and Coriolis terms so that the terms on the right hand side of the momentum Eqn 1 1 can be rewritten as for the x component o utan p gg 22h By e fuph vph o ph u av Ouph Ovph utan y a TA E ETA ul e Dy uph gt f uph 2 20 and for the y component o utang zz Ph ae fuph uph e ph du dv Ouph Ovph utang where is the relative vorticity The right hand sides consist of terms representing the energy gradient momentum convergence and curvature and of an altered Coriolis term that now contains the absolute vorticity instead of only the Coriolis parameter This term is included in the implicit formulation of the layer Eqn 2 40 which is the predictor step The residual terms are treated implicitly in the corrector step which also contains the momentum diffusion terms Since the momentum components are cross coupled in the Eqns 2 33 and 2 34 a linear equation in x is formulated simultaneously for both velocity components Alternatively a linear equation for the y direction may improve performance The discretizations for the flux divergence are those in 2 1 1 4 the energy gradient depending on u are approximated by the same formulae are valid for v ph u E ohu y r yy phu y 141 3 UY 141 3 UV 1 J 2 22 2
181. ple for a simulated precipitation is shown in Fig 1 2 EBM Precipitation Annual Mean T T T 30E 60E 90E 120E 150E 180 150W 120W 90W 60W 30W 0 30E 0 5 10 15 20 25 30 35 40 45 50 55 60 65 Unit cm month Josef M Oberhuber DKRZ Fri Mar 26 10 51 05 1999 C JMOPlane Figure 1 2 Annual mean precipitation as simulated by the EBM in a T42 resolution 1 6 3 Radiation Radiation is a key process to control the EBM s temperature distribution The difference of the short and longwave radiation at the top and the bottom of the atmosphere is the most important heat source Bulk parameterizations for the radiation at the surface are known to give reasonable spatial distributions However sensitivity experiments showed that their climate sensitivity is unrealistic Therefore new formulae have been developed that give similar realistic spatial distributions than earlier bulk formulae however behave more realistic in climates of higher water vapour contents see also Eqns 1 39 and 1 44 The formulae for the top longwave radiation is Qitop EoT ess 1 109 0 82 1 0 0 03 Vren where Tepp is considered to be some effective temperature at which the atmosphere radiates into space The idea is to let the atmosphere radiate at some weighted average between the rather constant temperature at tropopause level which here i
182. processes are connected with the problem of maintaining the coordinates in an isopycnal OGCM was given by JMO however the strategy has changed slightly in that the potential density in each layer is allowed to vary between the lower threshold 0 _ and the upper threshold o Tri Ton a O 1 66 Tora TGR Oy Hh 1 67 As long as the potential density cg varies between these threshold values only vertical mixing tries to attract the coordinate system towards o If the potential density Gop exceeds one of these thresholds then the convection procedure assigns water of that particular layer to a neigbouring layer such that the newly organized layer fullfills the limits Basically a layer is allowed to have a potential density error in a range between og Dori and og Con 4h However vertical diffusion is efficient enough to keep the error nearby zero unless a layer is very thin so that the assumption of linearity of the mixing process leads to some potential density error 1 4 1 Vertical Mixing Coordinate Maintenance Following the mixed layer parameterization it is assumed that turbulent kinetic energy is converted into mean potential energy via some u term If water is entrained only from above or only from below the equations for the entrainment rates wi and w are 2m u Co _ _ 14W 1 68 e Ge B Rieria Uk tix T 0 TED l 2m u A 1 69 gy4 hk Rieri Ug upt
183. ptun effect switch for vertical diffusion switch for mass budget gt sea level switch for mass conservation switch for sea ice calculation Table 4 19 Switches for Model Physics Variable ISW21 ISW27 ISW32 ISW37 Value 0 1 0 1 0 1 2 0 1 2 Default 1 Meaning switch for ADLR scheme switch for diapycnal pressure gradient switch for choice of transport scheme switch for momentum diffusion 1 2 0 Table 4 20 Switches for Numerical Schemes switch for mode of heat and fresh water forcing Reference lt OCEINIT gt lt OCEINIT gt lt OCEINIT gt lt SURFLUX gt lt OCEINIT gt lt SURFLUX gt lt SURFLUX gt lt BND gt lt OCEINIT gt lt OCESTEP gt Reference lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCEINIT gt lt OCEINIT gt lt OCEINIT gt lt OCEINIT gt lt OCESTEP gt lt OCESTEP gt lt EBM gt lt EBM gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt lt OCESTEP gt Reference lt OCESTEP gt lt OCESTEP gt lt MODEXP gt lt OCESTEP gt Variable Value Default Meaning Reference ISW11 switch for exit after topography generation lt OCEINIT gt ISW17 switch for generation of history file lt OCESTEP gt ISW21 switch for generation of vertical velocity lt OCESTOP gt ISW31 switch from binary to formatted OCDAT lt WRITEM gt ISW35 switch for emergency exit lt OCESTEP gt ISW46
184. qs 1 36 where Ta is the air temperature T the sea surface temperature qa the air specific humidity qs the specific humidity close to the surface which is assumed to be the saturated value pa is the surface air density and V is the scalar wind see also Eqn 1 64 Cp aiy is the specific heat of air and L the latent heat of evaporation The bulk coefficients cz and cz are calculated as proposed by Large and Pond 1982 The specific humidity q is given by the water vapor pressure e and the atmospheric surface pressure p In the following equations 1 37 1 44 e and p are given in Pascal and T is in K e 6ll x 1075 T 273 16 T 35 86 1 37 0 622 0 622 ds ces and q AS 1 38 p 0 378e p 0 378re where r is the relative humidity The respective formulae are given in the next section Radiative fluxes have been rewritten in order to achieve a realistic climate sensitivity This work has been done as response to greenhouse experiments with the atmospheric energy balance model EBM see later Basically it was found that factors of the type 1 an where a is some coefficient and n is a function of cloudiness and humidity yields to unrealistic fluxes if humidity is higher than in the present climate However if the former factor is considered as 1 Taylor expansion of 75 then a reversal of the factor s sign is avoided Then a super greenhouse effect is not simulated In the following all radiative fluxes
185. r it should be non zero in order to suppress B grid waves in the pressure field e The parameters DIFFUSV and DIFMINV specify the strength of momentum diffusion Note that momentum diffusion may be strong without severe consequences It might even cause a physical improvement e The parameters DIADIFO and DIADIF1 control the strength of the layer thickness diffusion A small value is physically justified and like SMOOTH is helpfull to remove B Grid waves from the layer thickness field e The parameters DIFFUSS and DIFMINS are required to smooth the temperature and salinity fields However physically meaningfull small values are sufficient to suppress numerical noise e The parameter NITMIN controls the quality of the iterative solutions If NITMIN is too small and or the time step too long then residual errors might accumulate over many time steps such that the model blows up finally The parameter NITMAX limits the number of iterations If the iteration count exceeds NITMAX for the wave equation then an error i e floating point exception is created This error is generated by the routine lt SUICIDE gt 5 4 10 Error Conditions The model contains a number of stop statements All of them give an indication on the kind of error It also may happen that the model exits with a floating point exception If the latter happens the possible reasons are No Convergence in Wave Equation As a consequence the routine lt SUICIDE gt is called
186. r partitioning and because the data traffic between server and clients is splitted over a higher number of npe_y UNIX pipes For more details please also read README 5 1 3 Installation Step No 3 5 1 3 1 Installing the Plot System In order to create pictures from model output the graphic library jmoplane f has to be installed For this execute the scripts libgen sun or libgen cri as examples for a SUN workstation or a CRAY on the directory Graphic Basically the script plot expects a file moplane o and utility o on the Graphic subdirectory Furthermore the files COAST LIN and COAST COL which contain the coastline geometry must be located on the same direc tory To be able to access these data adjust the open statement in jmoplane f Finally don t forget to adjust the file plot f in such a way that model geometry appears properly on the pictures 5 1 3 2 Installing the Postprocessing In order to install the postprocessing the scripts modgen sun or modgen cri must be executed separately for each utility The task of these scripts is to create object files which are saved Later when some postprocessing command is executed these object files are loaded and an executable is generated If the postprocessing is to be installed on a remote computer the script modgen cri gives an example on how scripts object files and documentation files are installed On the local machine the script modgen sun generates
187. re saved on HEQU and used either in lt MIXEXP gt or lt MIXIMP gt lt MECHAN gt computes the turbulent kinetic energy ENERGY the energy source due to the Garwood term TAUX and the inverse Ekman length scale EK MAN lt SHEAR gt computes the shear of the flow between mixed layer and the adjacent layer The result is used to calculate the shear instability 4 3 3 Active Tracer Subroutines lt TRACER gt organizes the time step for temperature and salinity lt TRACTIV gt computes the new distribution of temperature salinity or tracer con centration Since the underlying equations for all these quantities are identical and only differ in their forcing the routine is used for all different type of tracers The only dif ference is how the forcing is computed for the various tracers This is organized within lt TRACER gt lt BUBBLES gt computes how the ejected salt from freezing ice is distributed among the near surface layers The result is used to predict the new salinity via lt TRACER gt lt DECAY gt determines the heat flux into each layer that is the result of the penetrating solar radiation The result is used when lt TRACTIV gt is called to determine the new temperature distribution 4 34 The Sea Ice Model lt SEAICE gt computes the ice flow as well as the snow and ice cover from forcings due to stress heat flux and fresh water flux Note that the heat flux stored in QFLUX is already adjusted to the ice conditio
188. re updated during the iteration to account for the extreme nonlinearities in the viscous plastic rheology This means that bulk and shear viscosities are taken partly at the new time step For consistency with the continuity equation 9 point formulae are taken for the rheology The formulation of the model ensures that the small diffusion coefficient A in equations 1 77 to 1 80 is sufficient to guarantee computational stability 2 3 The Atmospheric Energy Balance Model The energy balance model works on the same grid and with the same time step as the ocean model As transport scheme the Crowley scheme is used for temperature transport Humidity is transported through flux form equations The according divergence operator is implicitly formulated as part of a matrix system 2 4 Tide Model The tide model is formulated on an Arakawa C grid however scalar points coincide with those of the ocean model In order to allow tides to be resolved in time an internally defined small time step is used The tide model is subcycled to synchronize it with the ocean model The iterative scheme is an alternating direction line relaxation ADLR scheme which gives good results even at surprising long time steps of e g 10 minutes at 1 degree resolution The derivation of the wave equation differs to that of the layer thickness equation such that the Coriolis terms are not eliminated from the pressure equation Instead the final equations are derived as follows
189. rogram are then independent of any other data source required to analyse the model fields 5 5 1 Code Definitions For a complete list of code definitions see block data statements in ocepict f or execute show 5 5 2 PPSF Binary Format The following write statements demonstrate how a data block for one quantity is generated Code Number Variable Unit Reference layer velocity components layer thickness HEIGTH layer potential temperature TEMP layer temperature layer salinity al SALT layer in situ density e DENSITY layer potential density z POTDENS layer in situ pressure 3 PRESS sea ice velocity components T UICE VICE sea ice thickness HICE sea ice compactness COMPACT snow surface temperature TSNOW snow ice interface temperature TICE snow cover ice equivalent HSNOW Table 5 1 Definition of Quantity Code Prognostic Variables Code Number Variable Unit Reference level potential density level velocity components U V level interface height m HEIGTH level potential temperature level in situ temperature level salinity sea surface temperature K TEMP I mixed layer depth HEIGTH sea surface elevation HEIGTH 1 sea surface velocity components z U 1 1 J V 1 1 J sea surface salinity SALT had if Lads J J J Table 5 2 Definition of Quantity Code Level Quantities DIMENSION IFLG N1 N2 X1 N1 X2 N2 FIELD N1 N2 N3 CHARACTER 8 LDATE LTIME CHARACTER 40 EXPERIM CALL DATE LDATE C
190. rs Name in Code Name in Formulae Meaning Unit Equation CMIXO linear damping term of ML CMIX1 constant TKE decay length scale CMIX2 critical Richardson number CMIX3 Ekman dissipation for TKE CMIX4 wave enery conversion factor CMIX5 Garwood term tuning coefficient CMIX6 constant buoyancy decay length scale CMIX7 Ekman dissipation for buoyancy CMIX8 threshold for hg CMIX9 tuning coefficient for buoyancy flux TURBEN tuning parameter HMIXMIN minimal mixed layer thickness Table 7 8 Names for Mixed Layer Parameters Name in Code Name in Formulae Meaning Unit Equation EDDYUV A diffusion coefficient EDDYS A diffusion coefficient RHOICE Pi ice density RHOSNOW Ps snow density CICE ki heat conductivity for ice CSNOW ks heat conductivity for snow CPICE Cpi specific heat capacity of ice CPMELT i heat of fusion EXCENT excentricity PRESICE A proportionality PDECAY proportionality HNULL tuning coefficient EPSNULL tuning coefficient AGING snow aging parameter SALTICE i salinity of sea ice PENET parameter for salt ejection Table 7 9 Names of Snow and Sea Ice Parameters Name in Code Name in Formulae Meaning Unit Equation ATU u zonal surface wind component ATV v meridional surface wind component UVABSOL V wind speed UVARIAN a V standard deviation of wind speed AIRT Ta surface air temperature SST sea surface temperature SSS S obs sea surface salinity CLOUDS cloudiness HUMID relative humidity SURPRES atm
191. rved salinity Sobs P is the precipitation E is the evaporation deduced from the latent heat flux Qz Cp the heat of fusion p the density of fresh water and R the river runoff The precipitation data are taken from Legates and Willmott 1990 1 3 4 Estimate of Turbulent Kinetic Energy Input The monthly mean absolute wind V and its standard deviation o V are available from the COADS by Wright 1988 or from the ECMWF analyses The standard deviation is required to accurately evaluate the time averaged third power of the friction velocity u which occurs in the relation for the entrainment rate Since 7 determined from the monthly mean absolute wind V only is much smaller than the required u3 an effective u3 must be determined by the additional use of the monthly mean standard deviation of the absolute wind By assuming that the amplitude of the fluctuations is not too large compared with the mean wind the effective u3 can be approximated by assuming a sinusoidal fluctuation around the mean value see JMO w be y 307 V 1 64 A similar relation is required to compute u for the neutral drag coefficient w She yey 20 V 1 65 1 4 Parameterizations of Internal Diffusion A parameterization for vertical diffusion in the interior ocean is introduced by allowing mass to be transferred between neighbouring layers An explicit procedure for mixing by convection is also included A first discussion of how the vertical mixing
192. s Fig 5 1 shows a typical example for a vector plot Each second arrow is taken out in order to improve the layout Here the land sea mask of the model is used to mark the continents Fig 5 2 is a standard example for a contour plot on a xy section Fig 5 3 is an example for the layout of a vertical section Vertical sections are plotted twice for the upper 500 m and for the whole depth However they are split up into two linear representations one for the upper 1000 mand one between 1000 m and 6000 m Fig 5 4 shows an example for a polar projection This is achieved by setting JPOL 3 in ALLIND to obtain a geometrical projection Otherwise a scalar quantity would be plotted as shown by Fig 5 2 If a global field is plotted the south and north polar projection are plotted on two separate sheets The south polar projection is not shown here The same projection is available also for vector plots Fig 5 5 shows an example for a fine resolution North Atlantic Arctic model in rotated coordinates 5 6 4 Usefull Tools The directory Tools contains 3 simple codes which help to define the model s coordinates in ocemain f These are e File makexy f which creates the horizontal grid using the same block data statements than used in ocemain f e File makever f which generates the potential density coordinates using the same block data statements than used in ocemain f e File plotgeo f which generates a plot that shows the
193. s assumed to be 50 C and some near surface temperature Thus with increasing humidity and or cloudiness the effective radiative temperature becomes lower The top shortwave radiation is parameterized by a Pon 22315 1 110 Ke ft S cos 7 d 2 ate Ss dt 1 111 Qatop 27 1 cos y 2 7 r5 1 085 cos 7 0 1 D 1 0 0 0019 noon 1 a AA a ar a 1 112 1 0 0 0019egua 1 0 0 59 1 0 0 03e n where Nequa 90 deg 1 6 4 River Runoff Model Equations The river runoff model follows the simple concept of a downhill flow whose speed is dependent on the topography gradient and the amount of water The prognostic equation for the water height Webm on land is Qi OW etm V Webm P E 1 113 Y GWeom rra 1 113 where Pis the precipitation E the evaporation and the last term being the melting freezing rate compare with Eqn 1 95 which converts ice or snow to water or vice versa The evaporation on land assumes that the surface is saturated with water for some Wem with a linear reduction of evaporation with decreasing Wim The water s flow speed is defined through the following approximation VO B Wotm VO VO 1 114 where O is the orography height 1 6 5 Glacier Model Equations So far the glacier model is realized by extending the functionality of the sea ice model While its dynamics is disabled on land the ice and snow mass as well as their thermodynamics are calcu
194. s for the ice covered and ice free conditions respectively The last term which is only non zero for ice free cases represents the influence of penetrating solar radiation on the total buoyancy flux Denman and Miyake 1973 The free parameter m represents the effi ciency of how turbulence available for mixing is produced by the mean wind stress Following Paulson and Simpson 1977 y is the fraction of solar radiation that penetrates through the ocean surface and hpg is the depth at which the penetrating radiation has decayed to 1 e In the relation for the entrainment rate the weighting coefficients a b are defined as exponential decay functions m exp hf 5u 1 25 exp hf Kus if B lt 0 i P Qu if B gt 0 1 26 Two different length scales are chosen a larger scale for those terms that create deepening and a smaller one for those that reduce the MLD This follows from the heuristic argument that wind stirring generates only ww terms at the surface whereas w T terms are also produced as a result of positive buoyancy fluxes They act as a source of TKE due to the generation of unstable stratification This kind of turbulence penetrates into the ocean and becomes a source of u w turbulence Thus positive buoyancy fluxes are assumed to be more efficient than wind induced turbulence or negative buoyancy fluxes An alternative is to specify a dissipation length scale which only depends on the buoyancy flux mi
195. s found that enlarges the diffusion in areas where the Kelvin or Rossby deformation radius is smaller than the grid resolution According to Munk s theory the diffusion coefficient is responsible for the width of the western boundary currents If the width is enlarged to a value comparable to the grid resolution a discrete model then is able to resolve western boundary currents Therefore if A once is tuned it is an universal constant of the model and thus there is no need to retune it if the model uses another resolution A wellcome side effect of momentum diffusion being bigger than scalar diffusion is that the two separate solutions on the B grid are coupled to each other A similar null mode does not exist for the variables 9ph and Sph In order to use the correct diffusion operators these have been modified such that e for momentum V AR pub has been exchanged by V A ph Vi while for heat and salt the respective operators V 4 ph VO and V A ph VS have already been used all the time however not been correctly documented in JMO The diffusion operator for momentum is still available as option to the newer momentum conserving velocity diffusion see Eqn 1 1 The new one is recommended because it gives a significant improvement of eddy physics In case the barotropic tide model is activated the bottom stress is computed from 7 PwC U Utide V Vide T Viide 1 1 2 Equation of State The potential t
196. s routine is needed to provide horizontal transport fields U H and VH containing diffusive mass fluxes for later postprocessing lt MODEXP gt determines the momentum advection and diffusion explicitly and adds their contribution to RESTU and RESTV Entry lt MODIMP gt performs the same calculation however in implicit formulation It uses the line relaxation method to treat these terms directly in x and z direction and iterates in y direction only or optionally solves these terms also in y and z direction and iterates in x direction lt MOLVERX gt performs an direct solution in x direction and an iterative one in y direction It is called by lt MODIMP gt lt MOLVERY gt performs an direct solution in y direction and an iterative one in x direction It is called by lt MODIMP gt lt DYNAMIC gt organizes the solution of the wave equation computes the contribution of expansion updates potential vorticity and calls lt SOLVERX gt and lt SOLVERY gt de pendent on whether the alternating direction line relaxation scheme is used or the solution is performed directly only in xz direction lt SOLVERX gt is one of the major routines of the entire model It performs one complete iteration of the wave equation and yields the layer thickness as the result It uses a line relaxation method that solves a linear equation system directly on xz slices and iterates in the y direction lt SOLVERY gt is one of the major routines of the e
197. salinity is parameterized by being related to the deformation of the flow field Our Our Our Ov vo s 1 2 2 v 0 A A om Dy Pay a A 1 10 Ou Ov Ou Ov AR ee AR E ae E e A 1 11 Ge ag ae 1 11 ur Ov ur v so _ s 1 en 2 2 8 0 A A Cor ay Cay ae A 1 12 The values for A t A and A are constants and have to be chosen according to the grid resolution see also table 4 5 The diffusion coefficients A A 0 and A depend on the grid spacing which are variable in longitude and latitude For momentum the underlying idea is to adjust the momentum diffusion such that the resulting frictional boundary layer is thick enough to be resolved by the grid This is achieved by relating the actual grid distance with Kelvin and Rossby deformation radii In the limit case of infinite resolution diffusion for vector and scalar quantities becomes comparable For scalar quantities such an automatic adjustment is not necessary The approach is 2 2 Ao a 28 44 Gedy 1 13 2 2 an Aia SP 4 Oud 1 14 A A 1 15 where A 470 A Here Az and Ay are the grid distances in the spherical coordinate system c is some phase velocity and is the planetary vorticity This formulation ensures a reasonable and automatic choice of diffusion coefficients in low and high resolution areas as well as in low and high latitudes By relating the grid distance with the deformation radii for Kelvin and Rossby waves a factor i
198. section lt WRITECNY gt writes eastern boundary data of a 3 dimensional array onto a 2 dimensional array containing an yz section lt STROBND gt organizes the copy of boundary data lt STRDATX gt copies an array with dimension NX lt STRDATY gt copies an array with dimension NY lt STRDATY gt copies an array with dimension NY lt STRDATXZ gt copies an array with dimension NZ NX lt STRDATYZ gt copies an array with dimension NZ NY lt SAMMELNX gt extracts an array with dimension NX from a 2 dimensional zx section Entry lt VERTEILX gt performs the reverse operation lt SAMMELNY gt extracts an array with dimension NY from a 2 dimensional zy section Entry lt VERTEILY gt performs the reverse operation lt SETVAL1X gt sets an array with dimension NX to a specified constant lt SETVAL1Y gt sets an array with dimension NY to a specified constant lt SETVALNX gt sets an array with dimension NZ NX to a specified constant lt SETVALNY gt sets an array with dimension NZ NY to a specified constant lt STORENX gt copies an array with dimension NZ NX lt STORENY gt copies an array with dimension NZ NY lt SELHMLX gt extract the mixed layer thickness on an x section from the 3 dimensional layer thickness array lt SELHMLY gt extract the mixed layer thickness on an y section from the 3 dimensional layer thickness array 4 7 5 Routines for Listing lt DRUCKHX gt prints out a scalar quantity
199. sends real data to the client e lt KERVER_GET gt receives integer data from the client Entry lt KERVER_PUT gt sends integer data to the client 4 8 2 2 Routines on the Client Side e lt cpipe_opn gt opens the command pipe The command pipe is used to synchronize the tasks on the server and the client Entry lt cpipe_cls gt closes the command pipe e lt client_opn gt opens the data pipes Entry lt client_cls gt closes the data pipes e lt client_get gt receives real data from the server Entry lt client_put gt sends real data to the server e lt klient_get gt receives integer data from the server Entry lt client_put gt sends integer data to the server 4 8 3 Josef s Private Interface Josef s Private Interface is a collection of routines that perform a data transfer either via the SHMEM or the MPI library The routines do the hard work for communication between the processors of an MPP Subsequently lx is some arbitrary vertical dimension nz_tot is the model s full dimension in x direction ny_tot is the model s full dimension in y direction and nx and ny are the x and y dimensions on each node e lt jpi_col_y gt collects data with dimension ny onto processor 1 1 and constructs data with dimension ny_tot Entry lt jpi_col_y_1 gt collects data with dimension ny onto processor 1 pe_y e lt jpi_col_x gt collects data with dimension nz onto processor 1 1 and constructs data with dimension nx_
200. sion of Scalar Variables 2 2 45 2 1 1 8 Horizontal Advection of Scalar Variables The Crowley Scheme 46 2A 2 Wiseretigation in TIMES Dad teal as AE RPO ARTE 46 2ZL21 Predictor Step aor eeo D A A Oe SS 46 2k22 Corrector Step ais a ar AAA oh h 48 22 Th Sea lce Mod l ia tas hot an he 4 de eo ai adas BAC BAD 49 2 2 1 Discretization in Space 4 gk A A A a Re eRe ee 49 2 2 2 Discretization in Time a as ata cad a ia a a 49 2 3 The Atmospheric Energy Balance Model 0 o o o 49 2A Dile Motel lt a bia e dl po dle OOK ARA As AR A AY 49 2 5 Open Boundary Formulation 22 44 5624 44540 48 4d 4b ee de sd 50 2 0 Computer Permermance waa ce a bie 4 0 4 Seek bee ek bw eee Se ele oe ee 51 ZO PVP ArcHitect r a 2 A BAe aoe we A ee d a Bie eee h 51 20 2 RISC AATEMLECTUES 2 a lada a ale tk ARA ita a A wey 52 20 3 MPP Architecture SATA EE RAE ARTE 4 52 II The PIPE System Components 55 3 Flow Diagrams 57 4 Model Code Description 63 All General Remarks 3 436 90 Se oP oe Se PAE ee AECE OO OS 63 4 2 Contents of ocemain f 256 2 1 a ai 63 4 2 1 Main Program MY 2OGO Moy IS a eae ee or AA AAA 63 4 2 2 Definition of Control Parameters 00202004 63 4 2 2 1 Block Data which concern Flow Control 64 4 2 2 1 1 Block Data SCREW 64 4 2 2 2 Block Data which concern Physical Parameters 64 4 2 2 2 1 Block DataHEATFLX 64 4 2 2 2 2 Block Data H
201. ss affected e The ratio between the theoretical speed versus the speed allowed by the implicit scheme indirectly appears as contribution to the matrices relative to the main diagonal element which is created by the time derivative The consequence is that the equations converge fast for slow modes while they converge much slower for fast waves This is because for fast waves the equations have the nature of a Laplace equation which is known to converge much slower than a Helmholtz equation e Due to the need to iteratively solve the equations for e g pressure the solutions are accurate for short waves but might contain non negligible errors on large scales Even the recently introduced ADLR scheme cannot remedy this problem fully As aresult of the reduced convergence for large scale fast modes the transport potentials which usually is a tiny residual of the flow might contain non negligible values Because errors are large scale the local velocity error should always be negligible Anyway the transport potential can be used to test the quality of the solution besides the stationarity of the ocean state To test stationarity analyse annual means only 5 4 12 About Code Consistency If one of the terms should be modified which is part of the predictor corrector scheme then one has to perform the modifications with care That s because those terms are formulated twice namely in an explicit and then in an implicit manner If one changes su
202. t lt VARDRAG gt REFHGT 10 reference height for transfer coefficient lt VARDRAG gt RKARMAN 0 4 Karman constant lt VARDRAG gt ALBEDO 0 07 albedo over water lt SURFLUX gt ABSORB 0 97 emissivity for longwave radiation lt SURFLUX gt ALBICE 0 50 albedo over snow free dry sea ice lt SURFLUX gt ALBSNOW 0 60 albedo over dry snow lt SURFLUX gt ALBJUMP 0 30 albedo jump from dry to wet ice lt SURFLUX gt Table 4 4 Heat Flux Parameters Variable Value Unit Meaning Reference DIFFUSS 1 E 9 m proportionality coefficient lt ISODIFF gt DIFFUSV 1000 s scaling coefficient for momentum lt MODEXP IMP gt DIFMINS 2000 s smallest scalar diffusion lt TRACER gt DIFMINV 5000 m s smallest vector diffusion lt MODEXP IMP gt DIADIFO 200 m s layer diffusion coefficient lt DIADIFF gt DIADIF1 1 E 10 m proportionality coefficient lt DIADIFF gt SLBGRID 0 2 m sea level smoothing parameter lt SOLVERX Y gt BAROINS Q m barotropic instability parameter lt ISODIFF gt Table 4 5 Parameters for Horizontal Diffusion coordinate pole out of the model area to avoid the coordinate convergence and the resulting deteriorating performance near the pole see table 4 11 Surface and bottom temperatures and salinities see table 4 13 define the vertical i e isopycnal coordinates Furthermore eliminating points allow one to remove lakes from the model area They are disabled if set to zero see table 4 14 Switches allow for standard proble
203. t OCESTEP gt NYEAR any gt 0 MODy completed years for model stop lt OCESTEP gt Table 4 1 Time Step Control Parameters 4 2 2 2 Block Data which concern Physical Parameters 4 2 2 2 1 Block Data HEATFLX defines all variables required to compute heat and fresh water fluxes see table 4 4 4 2 2 2 2 Block Data HORIMIX defines parameters required for computing the hori zontal diffusion coefficients for scalar and vector quantities see table 4 5 These parameters MUST be adjusted to the requirements of numerical stability which depends on the grid resolu tion Use Eqns 1 13 and 1 14 to determine which diffusion coefficients the model computes internally See also lt ISODIFF gt lt MODEXP gt lt MODIMP gt and lt TRACTIV gt Note that too large values for DIADIFO and DIADIF1 will result in physical rather than numerical in stability 4 lt 500m2s7 Variable Value Default Meaning Reference NINIT 0 1 2 switch for coupling interface JMOGLOB NDRUCK 0 1 switch for print JMOGLOB NHEAT 0 1 2 switch for heat coupling lt SURFLUX gt NSOL 0 1 2 switch for solar radiation coupling lt SURFLUX gt NFRESH 0 1 2 switch for salt coupling lt SURFLUX gt NSTRS 0 1 2 switch for stress coupling lt STRESS gt NMIX 0 1 2 switch for mixing energy coupling lt SURFLUX gt NICE 0 1 2 switch for snow amp sea ice coupling lt SURFLUX gt NAUVA 0 1 switch for scalar wind lt SURFLUX gt NAUVS 0 1 switch for standard deviation l
204. t SURFLUX gt NAIRT 0 1 switch for surface air temperature lt SURFLUX gt NAHUM 0 1 switch for surface air humidity lt SURFLUX gt NACLD 0 1 switch for cloudiness lt SURFLUX gt NOCET 0 1 switch for sea surface temperature lt SURFLUX gt NOCES 0 1 switch for sea surface salinity lt SURFLUX gt NPRES 0 1 switch for surface pressure lt SURFLUX gt o oooocoocooooooo 0 Table 4 2 Coupling Interface Parameters Variable Value Default Meaning Reference KHEAT 0 1 2 0 switch for heat fluxes lt SURFLUX gt KFRESH 0 1 2 0 switch for fresh water fluxes lt SURFLUX gt KSTRESS 0 1 2 0 switch for wind stress lt STRESS gt LENAVER any gt 0 100 switch for averaging period lt SURFLUX gt Table 4 3 Parameters for Annual Mean Flux Adjustment 4 2 2 2 3 Block Data ICE defines parameters for the snow and sea ice model see table 4 6 Note that CMELT and CFREEZ are tuning coefficients that allow changes in the rate at which heating causes the opening of leads or heat loss causes ice to form within leads Note that it is recommended to guarantee EDDYS EDDYUV see table 4 6 4 2 2 2 4 Block Data PHYSICS defines numerous parameters KTOTAL defines the model time when initialized TEMMEAN and SALMEAN determine the prescribed potential density that is used as potential density coordinate The layer properties if not initialized using Levitus data are chosen in lt CREATIV gt by assuming water with salinity of SALMEAN and a parab
205. t words via an IMPLICIT DOUBLE PRECISION statement The keywords for the precompiler are DOUBLE and _SINGLEW e Scalar or Vector Dependent on whether a computer uses a RISC or Vector processor DO loops are written differently to help compilers to obtain the best performance The keywords are SCALAR and VECTOR e CRAY or Not CRAY Computer Depending on the type of computer it is sometimes necessary to specify statements Therefore PIPE distinguishes between directives used by CRAY and non CRAY computers If a directive is used for a CRAY then the keyword is _CRADIR_ Currently other directives don t appear However if it is necessary to add directives for another PVP machine then precomp f already has the proper structure to allow further kind of directives e Choice of the Grid There exists the possibility to use pre defined grids as for T21 T32 T42 T63 T84 and T106 Gaussian grid Alternatively one may define a self defined grid see Table 4 11 The precompiler either activates or disables statements that are required for one of the choices The keywords are MYGRID for self defined grid _T21MOD _ for the T21 grid T32MOD for the T32 grid T42MOD_ for the T42 grid T63MOD for the T63 grid T84MOD for the T84 grid and _TI6MOD_ for the T106 grid e Cleaning from Disabled Codes To make the code more readable it is possible to delete all disabled statements while the new cod
206. tarctica are not part of the model domain however are needed to adjust the data format with that of other model data Script fit merges data on two input files containing data of the same type such that e g two scalar data sets are merged to one vector data set 5 5 6 8 Viewing Data Script contents prints the contents of the headers Script brief prints the statistics of each data set i e the means minima maxima and the standard deviations Script show prints the code definitions e Script list prints coordinates flag fields and data e Script plot creates plots of all quantities from an input file plot f needs to be adjusted to the model version by choosing parameters that define model geometry as shown by the plot Warning The dimension statements must be adjusted to the particular settings of the current model version 5 5 6 9 Data Preprocessor e Script lev2mod computes horizontal sections of temperature and salinity from the Levitus data file SALTEMP and delivers these quantities on the model grid The output file can be used to compute e g differences between the Levitus data and model data Warning The dimension statements must be adjusted to the particular settings of the current model version A script does not exist 5 5 6 10 Data Postprocessor The postprocessor consisting of the script makeall the code makeall f and the documentation file makeall doc is a general purpose tool to dia
207. tes of the Gaussian grid for the T106 resolution lt SPHERIC gt computes the weights of the Gaussian latitudes This routine is copy of the ECHAM model lt BESSEL gt computes the Bessel coefficients for lt SPHERIC gt This routine is a copy of the ECHAM model lt SHUTOFF gt switches off sea points between two points in geographical coordinates and thus creates a land brige between these two points lt CHANNEL gt opens a channel between two points in geographical coordinates and ensures a minimal depth along the line that connects these two points lt DEFMASK gt corrects the land sea mask for various boundary conditions lt GEOTOMO gt transforms a pair of x y coordinates from geographical into model coordinates on the sphere These routines are activated if the grid rotation is used They are needed to transform the forcing data onto the rotated model grid to compute the local Coriolis parameter and the daily mean insolation for each grid point Entry lt MOTOGEO gt performs the inverse operation lt ROMATIN gt computes the rotation matrix lt MO2GEOV gt performs a rotation from model to geographical coordinates for vector quantities Entry lt MO2GEOS gt performs a rotation from model to geographical coordinates for scalar quantities 4 4 3 Initialization of Atmospheric Forcing Each atmospheric quantity used to compute fluxes into the ocean is initialized by its own routine e lt DEFORCE gt is used to def
208. the particular model experiment With the help of the postprocessing output rou tines see section 4 5 it can be decided which quantities are saved on the history file and the frequency of storing them 5 5 4 2 2 File PPSF_HIS is the history file of the tracer model 5 5 5 getf and putf The scripts getf and putf are used in numerous scripts Both scripts perform a ftp to another machine where getf is trying to receive a file while putf is trying to send a file In both scripts the length of the local file and the remote file are compared with each other If both have the same size the transfer is finished Possibly the scripts have to be adjusted to the operating system so that the file size in bytes is correctly cut with the awk command from the line created by the ftp s dir command 5 5 6 Utilities for Data Analysis The following routines allow to manipulate and convert the contents of a PPSF file and to perform mathematical operations See also the file README in PostPro for further infor mation Unless noticed otherwise all codes below work as long as the dimensions used in the PPSF format does not exceed the dimensions used in the underlying routines A few routines expect that the dimensions are identical with that used in the input data see warnings All utilities have a simple online documentation Use extension h to view the help file instead of performing any operation For instance filter h results in
209. third chapter flow diagrams show how the model is organized in the fourth chapter the model code is described and finally the fifth chapter explains how PIPE is to be installed and applied to self defined applications Part of this last chapter also is the usage of the postprocessing software Notation The subsequent descriptions reference the code in the following way e a variable is written in italic type e g TEMP for potential temperature e a SUBROUTINE is referenced by enclosing the name within lt and gt eg lt SOLVERX gt for the solution of the block tridiagonal system in x direction e a BLOCK DATA program unit is abbreviated by enclosing the name by e g SCREW for defining the time step a COMMON block name is enclosed by e g JMOBAS which contains the prog nostic variables of the model e a file or directory name is enclosed by and e g ocestep f for the kernel of PIPE 1 1 The Interior Ocean Model PIPE is a general circulation model based on the following ideas Isopycnals are used as La grangian vertical coordinates a realistic equation of state is included the primitive equations together with the hydrostatic approximation are applied and a surface mixed layer and a snow and sea ice model are coupled to the interior ocean The equation of state together with the 15 sea ice model allows the model to be used on global scales No approximation is made except
210. thus are constructed in such a way that they automatically conserve heat and water independent on a proper parameter tuning 1 6 1 Transport Equations for Temperature and Humidity The EBM s vertically integrated temperature Tem and water vapour contents Vs are predicted through the following equations OT in Qebm at VTam V Actin VT com 1 102 9 1U bm Actm bm bs peli Ot 2V Vebm V Act V Vebm E P 1 103 where Qem 1s the sum of all the heat fluxes at the top and the bottom of the atmosphere and the internal heat fluxes Qebm Jaden Qi top Qs bot Qi bot Qu F Qatm 1 104 where Q top is the top solar radiation see Eqn 1 111 Qitop the top longwave radiation see Eqn 1 109 Qs sor the solar radiation at the surface see Eqn 1 44 Qi bot the longwave radiation at the surface see Eqn 1 39 Qg the sensible heat flux into the ocean or soil see Eqn 1 35 Qatm is the sum of latent heat release due to condensation precipitation and adiabatic cooling due to vertical motion under the presence of some assumed atmospheric vertical stratification E and P are the evaporation and precipitation at the earth s surface respectively Asp 1s a turbulent horizontal eddy diffusion coefficient Hy is the height of the atmosphere 1 6 2 Diagnostic Relations The atmospheric eddy diffusion coefficient A is parameterized as Aebm Vo V 1 105 with V the observed scalar wind and o
211. tidal sea level anomalies from data are saved 5 File TRAFORCE contains data for an offline running tracer model which is basically built around the routine lt TRACPAS gt The tracer model obtains flow and isopycnal data from the file TRAFORCE as monthly forcing data 6 File PPSF_ALL is the history file which is created depending on the switch SW17 Other files with an initial PPSF_ are created if one of the switches ISWPOST is set for details see Quick Look System 5 3 Examples for Specific Model Layouts In this section it is demonstrated how simple it is to setup PIPE for any geometry and forcing Starting with any model version of PIPE the code preprocessor precomp f has to be called by setup to ensure that the correct grid type pre or self defined is enabled Also dimensions have to be specified The final changes only concern ocemain f The definition of numerical parameters such as the time step or physical parameters such as the horizontal diffusion coefficient have to be chosen according to the application and numerical limitations High resolution grids require some adjustments 5 3 1 Pre defined Grids In addition to the rules to set up a model version the following parameter setting is required or recommended 1 One of the T21 T32 T42 T63 T84 or T106 grids has to be enabled by setup 2 Ensure that the following dimensions have been chosen e T21 Grid NX 66 and NY 34 e T32 Grid NX 9
212. tion the initial data on these points are taken from Levitus Also ensure that the initial model time defined by KTOT AL agrees with the model time of the old restart file OLDFILE Otherwise a discontinuity of the forcing between the model time of OLDFILE and the model time during the initialization may result Note that a proper definition of KTOTAL allows the definition of the model time during initialization in general 5 4 5 Diagnostic Output During the model run a number of formatted files appear that allow a quick check of the model physics These are 1 File OUTLIST contains the protocol of the model After each time step crucial num bers mainly for monitoring stability properties are written out See description in OUTLIST 2 File LSMASK contains the land sea mask written out during the initialization of the topography if S W45 1 This file can be read again and is used to overwrite the land sea mask during a second topography initialization if SW45 2 3 File MEANS contains layer averages for temperature salinity in situ density and potential density These data are written out by lt STATIST gt at each time step It is called by lt OCESTEP gt 4 File PERFORM contains an overview of CPU times from the most CPU time con suming routines 5 File BUDGET contains the heat and fresh water budget as well as conservation errors for each time step The data are written out from lt OCESTEP gt 5 4 6
213. tmospheric Energy Balance Model The atmospheric energy balance model represents a cheap replacement of todays full blown atmospheric GCMs The neglect of dynamical feedbacks between pressure and flow fields as well as the resignment on detailed vertical resolution yields a one layer atmospheric model that predicts the horizontal distribution of temperature and humidity only In order to keep such a model sufficiently realistic the EBM is data driven by prescribing wind and cloudiness All other fields required to close the global energy budget as due to radiation precipitation sensible and latent heat exchange between ocean and atmosphere and diabatic and adiabatic heating are diagnosed In order to conserve energy soil and river runoff models predict temperature water content in the soil and water flow towards the oceans At high latitudes the variables of the sea ice model are used to simulate glaciers with snow on top So far various approaches have been realized to obtain a simple atmosphere model without dynamics but with more realistic feedbacks than possible with a Newtonian relaxation North et al 1983 Kleeman and Power 1995 which also don t contain the surprising sensitivity of mixed boundary conditions to arbitrary parameter choices as discussed by Mikolajewicz and Maier Reimer 1994 The goal is to develop an EBM that has a closed hydrological cycle and conserves the sum of sensible and latent heat Some of the parameterizations
214. tor point is switched from a sea to a land point at the edge of a zero layer dara daa a A Ee Be a 41 Concept of Server Client Interaction oaoa 0 0 0 53 Concept for Exchanging Data between Processors oouo a 54 Concept for Iterating on Neighbouring Processors oaoa aaa a 54 General Overview of Flow 6 Sok lia bot eal Pa a 58 Flow of Routine OCEINIT 2 Litas Bode Se eek rt Race wok 4 59 Flow of Routine OCESTEP a asses o wat eae A eas 60 Flow of Routine OCGEPOST g4ep e a eee kee ek ee AA 61 Flow of routine OCHS TOP 2 447 08 bed Gt bee Ee Eee EE te 61 Example showing a vector plot for the surface flow of PIPE T42 The picture is created interactively using plot b e q 69 a Annual Mean PPSF_SUR 141 Example showing a greyscale colour plot for the sea level of the PIPE T42 The picture is created interactively using plot b p e q 68 a Annual Mean TAPES CTs aue O a we deo Be Enel AS Ge ae BOS htm RIG dae ein Hah eh 141 Example showing a yz section for the potential temperature at 120W from the PIPE T42 The picture is created interactively using plot b p e q 13 1 240 FEES VTE TE S pai te 18 a odos ihe non Be BE aks LS Pie Roe Ge de Pe Was ele hose 142 Example showing a north polar projection for the sea ice thickness of PIPE T42 The picture is created interactively using plot b p e r q 18 N PPSF_ICE 143 9 5 5 Example showing the surface flow in a North Atlantic Arctic mode
215. tot e lt jpi_col_xy gt collects data with dimension lz ng ny onto processor 1 1 and constructs data with dimension lz nx_tot ny_tot Entry lt jpi_col_xy_1 gt collects data with dimension lx nx ny onto processor 1 pe_y and constructs data with dimension lz nx_tot ny_tot e lt jpi_dis_y gt distributes data with dimension ny tot from processor 1 1 and constructs data with dimension ny Entry lt jpi dis_y_1 gt distributes data with dimension ny_tot from processor 1 pe y e lt jpi_dis_x gt distributes data with dimension nx_tot from processor 1 1 and constructs data with dimension nz e lt jpi_dis_xy gt distribute data with dimension lx nx_tot ny tot from processor 1 1 and constructs data with dimension lx nx ny Entry lt jpi_dis_xy_1 gt distributes data with dimension lz ng tot ny_tot from processor 1 pe_y and constructs data with dimension lz nz ny e lt jpi_zon gt performs a zonal data exchange e lt jpi_mer gt performs a meridional data exchange e lt jpi_mer_s gt performs a meridional data exchange for scalar quantities e lt jpi_mer_v gt performs a meridional data exchange for vector quantities lt jpi_shift gt performs a zonal data shift lt jpi sum gt performs a global sum based on a binary tree algorithm lt jpiz max gt performs a global minimum based on a binary tree algorithm lt jpi_flood gt performs a broadcast from processor 1 1 to all others lt jpi_merge gt collects data
216. ty Table 7 1 Names for Prognostic Ocean Quantities Name in Code Name in Formulae Meaning Unit Equation U zonal velocity V meridional velocity DENSITY in situ density POTDENS potential density PRESS in situ pressure DIFFUSE s scalar diffusion coefficient VERTUP upward diffusion VERTDN downward diffusion LEVUPP index of next upper layer LEVLOW index of next lower layer RESTU Fah right side of zonal flux RESTV Boh right side of meridional flux RESTH Fon right side of wave equation POTSOLL o potential density coordinate Table 7 2 Names for Diagnostic Ocean Quantities 149 Name in Code Name in Formulae Meaning Unit Equation TAUX zonal wind stress TAUY meridional wind stress QFLUX net surface heat flux QSOLAR downward solar heat flux QSENSIB surface sensible heat flux QLATENT surface latent heat flux QLONG longwave radiation budget SFLUX equivalent salt flux BFLUX net surface buoyancy flux BSOLAR solar buoyancy flux USTERN friction velocity HEATLAT heat of fusion SIGMA Stefan Boltzmann constant SOLAR solar constant CPAIR specific heat capacity CHARNCK Charnock constant RKARMAN von Karman constant ALBEDO albedo ABSORB emmisivity of water TURBID penetration depth of Qs tor m TURBREL transmission coefficient Table 7 3 Names for Surface Flux Parameters Name in Code Name in Formulae Meaning Unit Equation HEATQ Qrun running mean heat flux HEATM ias heat flux bias correction PMEQ Fes running
217. ty sea level mixed layer depth temperature and salinity 5 5 4 1 2 File PPSF_ICE contains the ice drift the ice thickness the ice compactness the snow cover the snow surface temperature and the snow ice interface temperature Code Number Type of Section Reference section for xy plane of instantaneous variable lt PPSFOUT gt section for xz plane of instantaneous variable lt PPSFOUT gt section for yz plane of instantaneous variable lt PPSFOUT gt integral of yz plane over x lt PPSFOUT gt integral of xz plane over y lt PPSFOUT gt integral of xy plane over z lt PPSFOUT gt monthly means on xy section lt PPSFOUT gt monthly means on xz section lt PPSFOUT gt monthly means on yz section lt PPSFOUT gt line section in x direction lt PPSFOUT gt line section in y direction lt PPSFOUT gt line section in z direction lt PPSFOUT gt integral along x over yz section lt PPSFOUT gt integral along y over xz section lt PPSFOUT gt integral along z over xy section lt PPSFOUT gt monthly means along line section in x direction lt PPSFOUT gt monthly means along line section in y direction lt PPSFOUT gt monthly means along line section in z direction lt PPSFOUT gt CON MDM OB WM eR Table 5 8 Definition of Section Code Part I 5 5 4 1 3 File PPSF_BAS contains the momentum fluxes layer thicknesses temperature and salinity in original model coordinates This file can be used by makeall f to compute the horizontal
218. ty is taken The detrainment procedure is implemented so that during the predictor step a detrainment rate is computed from hig h 1 32 w At and is used to force the continuity equation while in the corrector step it is ensured that hit hi7 from which an effective entrainment detrainment rate is determined 1 2 3 Coupling of the Mixed Layer to the Interior Ocean Entrainment is a process that does not conflict with the maintenance of the isopycnal coordinate system The conceptual problem however is to select one of the layers below the ML into which the detrained water is pumped The strategy is described in JMO and can be found in lt MIXEXP gt lt MIXIMP gt and lt CONVECT gt In detail the strategy has changed in that mixed layer water with some potential density dj is detrained into that layer such that Trt lt 04 lt Oeti see Eqn 1 66 and 1 67 However first it is tried to reduce the potential density error dy 0 in the first physically existing layer below the mixed layer 1 3 Parameterization of the Surface Fluxes The surface fluxes required by the model are the fluxes of momentum heat fresh water tur bulent kinetic energy and buoyancy The required data sets are surface air temperature sea surface temperature relative humidity cloudiness the time averaged absolute wind speed and its standard deviation wind stress precipitation air surface pressure and surface salinity The flux of moment
219. um is taken from Hellermann and Rosenstein 1983 or from ECMWF as global stresses Based on the COADS Comprehensive Ocean Atmosphere Data Set Woodruff et al 1987 Wright 1988 prepared data on a 2 x 2 grid which is sufficiently fine for forcing a coarse resolution OGCM These fields represent a 30 year climatology from 1950 to 1979 and have been extended to 1986 Alternatively the same quantities are also made available from the ECMWF Reanalysis Project Gibson et al 1997 which are available on a T106 Gaussian grid 1 3 1 Parameterization of the Surface Heat Fluxes The total heat flux into the mixed layer Q is given by Q1 Qu Qr Quito Qs bor 1 yezp h1 hg 1 33 and consists of the sensible heat flux Q y the latent heat flux Qz the net heat flux by longwave radiation Qi bo and the heat flux Q so due to insolation y defines the fraction of the insolation that is not immediately absorbed at the surface but penetrates into the ocean Following Paulson and Simpson 1977 hg is the decay length scale for the absorption of solar radiation This means that when all upper layers are shallow enough deeper layers gain heat due to insolation The heating rates of all deeper layers are defined by Qt Quialezp Y hu hs exp Sha hs 7 1 34 The turbulent surface heat fluxes namely the sensible and latent heat fluxes are estimated by the bulk formulae Qu Patri ca V Ta A ES 1 35 QL PalwerV Ga
220. usion Parameters 151 Names for Mixed Layer Parameters o 151 Names of Snow and Sea Ice Parameters 152 Names of Observed Quantities 0 aa a a 152 Names of EBM Parameters o 153 Names of Tide Parameters aoaaa aa a a a a 153 List of Figures 1 1 1 2 1 3 1 4 2 1 2 2 2 3 2 4 3 1 3 2 3 3 3 4 3 5 5 1 5 2 5 3 5 4 Meridional cross section of interface height for the example of a T42 global ocean simulation The section is located at 120W The interfaces are plotted as steps in order to show where grid cells are located Each step denotes the location of a vector point The brightest shades denote sections through Antarctica and NOU AMELIA nas e el a Beep aa ek A Bea ee Be a 16 Annual mean precipitation as simulated by the EBM in a T42 resolution 34 Amplitudes of the M2 Tide simulated with a T106 equivalent version of the barotropic tide model coupled to the 3 dimensional PIPE model 36 Phases of the M2 Tide simulated with a T106 equivalent version of the barotropic tide model coupled to the 3 dimensional PIPE model 2 37 Schematic layout of the B grid L represents land points S represents sea points Circles denote scalar points crosses denote vector points The thick line marks the boundary at a certain time level The dashed line shows an example of how the boundary changes if a vec
221. ut the results e lt QTIDST gt bookkeeps and writes out residual currents 4 5 4 Output of Instantaneous Fields The following routines can be called to write out model fields after any transformation Some of them already compute spatial averages and thus write out fields with reduced dimension The first character of the name denotes whether a vector or a scalar is treated The quantity code number must be specified however the routines specify and store the correct section code internally e lt VALL gt writes out all layers of a model array without any transformation and should be used if one wants to reconstruct model quantities during the data postprocessing e lt SALL gt writes out scalar quantities e lt VXYALL gt writes out all levels of vertically interpolated vector quantities Entry lt SXYALL gt writes out all levels of vertically interpolated scalar quantities e lt VXYONE gt writes out a precalculated xy section at a given depth Entry lt SXYONE gt prints out scalar field e lt VXONE gt writes out a precalculated one dimensional vector in x direction Entry lt SXONE gt prints out scalar field e lt VYONE gt writes out a precalculated one dimensional vector in y direction Entry lt SYONE gt prints out scalar field e lt VXSEC gt vrite outs a one dimensional vector in x direction Entry lt SXSEC gt prints out scalar field e lt VYSEC gt writes out a one dimensional vector in y direction Entry
222. water flux needs to be converted to a buoyancy flux derivatives of in situ density relative to temperature or salinity are required These derivatives are com puted analytically as both equations for potential temperature and in situ density are simple polynomials 1 2 The Mixed Layer Model An important aspect of the model implementation is that a large variety of mixed layer models can be invoked JMO described only one possible choice Various models such as Kraus Turner 1967 Niiler 1975 Garwoord 1977 Niiler and Kraus 1977 Garwood et al 1985a b or Gaspar 1988 can be approximately obtained by an appropriate choice of numerous parame ters 1 2 1 Physical Background A mixed layer ML is the result of turbulence produced by wind stirring and surface buoyancy fluxes Temperature salinity velocities and tracers are then uniformly distributed in the vertical Kraus and Turner 1967 The turbulent kinetic energy TKE is partly converted into mean potential energy MPE and partly dissipated One of the standard problems is that typical mixed layer models are tuned to particular local conditions e g ocean station Papa see Martin 1985 There is at present no general formulation with a satisfying theoretical justification that is able to treat simultaneously the significantly different mixed layer regimes in different parts of the ocean A summary of some difficulties concerning earlier ML models is given in Gaspar 1988
223. xplicitly in the model All procedures assume that 1 J NX J 1 and J NY contain the correct values for any differentiation The subsequent routines are called whenever required to store the proper values on the grid points along the model margin The land sea mask F LG allows the model to deduce the proper boundary conditions inside the model area lt BNDH1 gt There are three different versions of the routines that provide the boundary conditions along the model margin lt BNDH1 gt is used for scalar quantities on scalar grid points Note that the extension 1 is taken to express that only arrays containing one layer are treated by these routines Entry lt BNDU1 gt is used for scalar quantities on vector grid points Entry lt BNDV1 gt is used for vector quantities on vector grid points lt BNDHNZ gt The following routines are similar to the above ones however full model arrays are treated with these routines Entry lt BNDUNZ gt works like lt BNDU1 gt but for all NZ layers Entry lt BNDHNZ gt works like lt BNDH1 gt but for all NZ layers e lt UVPOLE1 gt These routines are called from lt BND gt and compute the northern boundary values on vector points in case that the open North Pole is switched on The UV marks the routine used for vector points and S the one for scalar points The 1 at the end is used since these routines treat one layer only Entry lt SPOLE1 gt is used for scalar quantit
224. xtra grid point to the east west south and north of the model domain It also must be ensured that the focus is switched off first To avoid bank conflicts on the CRAY there are two limitations 1 The number of layers NZ must be an odd integer 2 The number of latitudes NY must be chosen in such a way that the FORTRAN expression NV NY 1 2 yields an odd value for NV Other computers may have other restrictions Performance on Cache type computers suffers mainly through the large address increments of non contiguous vectors In most cases N Z is the address increment If the model is implemented on a cache based i e RISC processor then alternative formulations ensure that the address increments are contiguous This is ensured during the installation of PIPE as a RISC version 5 2 4 Defining a Focus A focus is an area of enhanced resolution The code allows for any definition of the x and y coordinates as long as the values increase monotonically For numerical stability however it is necessary to change the grid distance as smooth as possible This means that the ratio of the left hand and right hand side grid distance should be nearly unity This is ensured to some extent by a spline technique to compute the coordinates The protocol OUTLIST provides the x and y coordinates and the relative changes of the grid distances when the model initializes itself The focus parameters and the dimensions for NX and NY should be changed if the
225. y Pe ph PP ice a where g is the gravitational constant and Paty is the atmospheric pressure which optionally can be added to the oceanic pressure field In all deeper layers the in situ pressure P is recursively defined by Py Prat oh 1 ph e 1 8 The pressure gradient in a layer k is computed from the sum of gradients due to the sea level gradient the gradients of the interface heights above the layer k and the potential density gradients integrated from the surface down into the center of layer k The diagnostic equation for the pressure gradient is N mw k 1 S h 3 Ve 90 Whi VD gE S vinos Von 1 9 l k 6k j 4 where D is the height of the topography above some reference level The concept is to let the pressure gradient be valid for an arbitrary distribution of density and potential density in each layer As explained in JMO this results in keeping the model dynamics valid in situations when isopycnal coordinates cannot be maintained Therefore gradients of potential density dy are kept which is not necessary if one a priori assumes that the vertical coordinates always are isopycnal However the latter a priori assumption leads to conceptual difficulties in the formulation of vertical mixing processes Al and Aj represent the spatial and time dependent diffusion coefficients for vector layer and tracer quantities Following Smagorinsky 1963 the horizontal diffusion coefficient for temperature and
226. y grid point by changing the mean diagonal element artificially without changing the final solution An issue worth noting is that there exists a concept to derive the correct boundary conditions used for solving the pressure Eqn 2 40 The concept is explained in Appendix B with the help of Appendix A The same concept is used as part of the tide model which is an application that demonstrates that external gravity waves are realistically reflected along coasts see also Fig 1 1 and 1 2 2 1 2 2 Corrector Step The entrainment and detrainment rate and the resulting changes in the mass fluxes are treated implicitly in the corrector step The solution for the mixed layer is obtained with Newton s method to find zeros of the resulting nonlinear equation Because no spatial derivatives occur the following equations can be solved pointwise gt n l At gt n l At puh Gent puh 5G 2 43 ony Eos oh Eo 2 44 prayer Samt phy Ear 2 45 hs Ea hs Lar 2 46 gt x wn 1 puh ph RL pho and phS are the corrected guesses for the new time level n 1 2 2 The Sea Ice Model 2 2 1 Discretization in Space The sea ice model works on the same grid and with the same time step as the ocean model It is written in spherical coordinates This means that the momentum diffusion is discretized according to 2 1 1 5 the sea level pressure gradient according to 2 1 1 3 the ice mass and c
227. zation After setting SW11 0 SAVFORC 1 UPDATED 1 and hav ing chosen VELFORC TEMFORC and SSTFORC consistent with SW30 ISW43 and ISW47 run model sun or model cri After completion of this short run it delivers FORCES and OCDAT as well as all quick look files 4 Continuation Run After having set INITIAL 0 SAVFORC 0 AGEDATA 1 and UPDATED 2 the run can be continued from the previously computed state saved on OCDAT Note that if OCDAT is empty the model tries to initialize the topography Also note that the initial time step DT MIN should be small The model will increase the time step according to an algorithm that uses the convergence of the wave equation as measure The model can be stopped at a certain date by specifying NYEAR NMONTH etc 5 Production Run The easiest way to synchronize the data output with the time stepping is to choose DTMIN DTMAX This ensures that the model does not vary the time step according to some numerical measure However it must be ensured that the model works properly for the entire production run with these chosen values With MAXREP the user can specify the integration time per job run 5 4 3 Data Preprocessing In order to avoid extra work for computing the forcing for individual model layouts all forcing data are held only as global fields Depending on the chosen grid the model creates its forcing during initialization and saves its individual forcing on FORCES Thus no work has to

Technical Report No. 19 PIPE ISSN 0940-9327

Contents

Download Pdf Manuals

Related Search

Related Contents