/[MITgcm]/MITgcm/model/src/ini_parms.F
ViewVC logotype

Contents of /MITgcm/model/src/ini_parms.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.36 - (show annotations) (download)
Fri Nov 6 22:44:47 1998 UTC (25 years, 6 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint18
Changes since 1.35: +5 -4 lines
Changes to allow for atmospheric integration builds of the code

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_parms.F,v 1.35 1998/11/02 03:34:12 cnh Exp $
2
3 #include "CPP_OPTIONS.h"
4
5 SUBROUTINE INI_PARMS( myThid )
6 C /==========================================================\
7 C | SUBROUTINE INI_PARMS |
8 C | o Routine to set model "parameters" |
9 C |==========================================================|
10 C | Notes: |
11 C | ====== |
12 C | The present version of this routine is a place-holder. |
13 C | A production version needs to handle parameters from an |
14 C | external file and possibly reading in some initial field |
15 C | values. |
16 C \==========================================================/
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "EEPARAMS.h"
21 #include "PARAMS.h"
22 #include "GRID.h"
23 #include "CG2D.h"
24
25 C === Routine arguments ===
26 C myThid - Number of this instance of INI_PARMS
27 INTEGER myThid
28
29 C === Local variables ===
30 C dxSpacing, dySpacing - Default spacing in X and Y.
31 C Units are that of coordinate system
32 C i.e. cartesian => metres
33 C s. polar => degrees
34 C goptCount - Used to count the nuber of grid options
35 C (only one is allowed! )
36 C msgBuf - Informational/error meesage buffer
37 C errIO - IO error flag
38 C iUnit - Work variable for IO unit number
39 C record - Work variable for IO buffer
40 C K, I, J - Loop counters
41 C xxxDefault - Default value for variable xxx
42 _RL dxSpacing
43 _RL dySpacing
44 CHARACTER*(MAX_LEN_MBUF) msgBuf
45 CHARACTER*(MAX_LEN_PREC) record
46 INTEGER goptCount
47 INTEGER K, I, J, IL, iUnit
48 INTEGER errIO
49 INTEGER IFNBLNK
50 EXTERNAL IFNBLNK
51 INTEGER ILNBLNK
52 EXTERNAL ILNBLNK
53 C Default values for variables which have vertical coordinate system
54 C dependency.
55 _RL viscArDefault
56 _RL diffKrTDefault
57 _RL diffKrSDefault
58 _RL hFacMinDrDefault
59 _RL delRDefault
60 C zCoordInputData - These are used to select between different coordinate systems.
61 C pCoordInputData The vertical coordinate system in the rest of the model is
62 C rCoordInputData written in terms of r. In the model "data" file input data can
63 C coordsSet be interms of z, p or r.
64 C e.g. delZ or delP or delR for the vertical grid spacing.
65 C The following rules apply:
66 C All parameters must use the same vertical coordinate system.
67 C e.g. delZ and viscAz is legal but
68 C delZ and viscAr is an error.
69 C Similarly specifyinh delZ and delP is an error.
70 C zCoord..., pCoord..., rCoord... are used to flag when z, p or r are
71 C used. coordsSet counts how many vertical coordinate systems have been
72 C used to specify variables. coordsSet > 1 is an error.
73 C
74 LOGICAL zCoordInputData
75 LOGICAL pCoordInputData
76 LOGICAL rCoordInputData
77 INTEGER coordsSet
78
79 C-- Continuous equation parameters
80 NAMELIST /PARM01/
81 & gravity, gBaro, rhonil, tAlpha, sBeta, f0, beta,
82 & viscAh, viscAz, viscA4,
83 & diffKhT, diffKzT, diffK4T,
84 & diffKhS, diffKzS, diffK4S,
85 & GMmaxslope,GMlength,GMalpha,GMdepth,GMkbackground,GMmaxval,
86 & tRef, sRef, eosType,
87 & momViscosity, momAdvection, momForcing, useCoriolis,
88 & momPressureForcing, metricTerms,
89 & tempDiffusion, tempAdvection, tempForcing,
90 & saltDiffusion, saltAdvection, saltForcing,
91 & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz,
92 & tempStepping, saltStepping, momStepping, implicitDiffusion,
93 & viscAr, diffKrT, diffKrS, hFacMinDr,
94 & viscAp, diffKpT, diffKpS, hFacMinDp,
95 & rhoConst, buoyancyRelation
96
97 C-- Elliptic solver parameters
98 NAMELIST /PARM02/
99 & cg2dMaxIters, cg2dChkResFreq, cg2dTargetResidual, cg2dpcOffDFac,
100 & cg3dMaxIters, cg3dChkResFreq, cg3dTargetResidual
101
102 C-- Time stepping parammeters
103 NAMELIST /PARM03/
104 & nIter0, nTimeSteps, deltaT, deltaTmom, deltaTtracer, abEps,
105 & tauCD, rCD,
106 & startTime, endTime, chkPtFreq, dumpFreq, taveFreq, deltaTClock,
107 & pChkPtFreq, cAdjFreq, tauThetaClimRelax, tauSaltClimRelax,
108 & periodicExternalForcing, externForcingPeriod, externForcingCycle
109
110 C-- Gridding parameters
111 NAMELIST /PARM04/
112 & usingCartesianGrid, delZ, dxSpacing, dySpacing, delX, delY,
113 & usingSphericalPolarGrid, phiMin, thetaMin, rSphere,
114 & l, m, n, delP, delR, rkFac
115
116 C-- Input files
117 NAMELIST /PARM05/
118 & bathyFile, hydrogThetaFile, hydrogSaltFile,
119 & zonalWindFile, meridWindFile, thetaClimFile,
120 & saltClimFile
121
122
123 C
124 _BEGIN_MASTER(myThid)
125
126 C-- Initialise "which vertical coordinate system used" flags.
127 zCoordInputData = .FALSE.
128 pCoordInputData = .FALSE.
129 rCoordInputData = .FALSE.
130 usingPCoords = .FALSE.
131 usingZCoords = .FALSE.
132 coordsSet = 0
133
134 C-- Open the parameter file
135 OPEN(UNIT=scrUnit1,STATUS='SCRATCH')
136 OPEN(UNIT=scrUnit2,STATUS='SCRATCH')
137 OPEN(UNIT=modelDataUnit,FILE='data',STATUS='OLD',
138 & IOSTAT=errIO)
139 IF ( errIO .LT. 0 ) THEN
140 WRITE(msgBuf,'(A)')
141 & 'S/R INI_PARMS'
142 CALL PRINT_ERROR( msgBuf , 1)
143 WRITE(msgBuf,'(A)')
144 & 'Unable to open model parameter'
145 CALL PRINT_ERROR( msgBuf , 1)
146 WRITE(msgBuf,'(A)')
147 & 'file "data"'
148 CALL PRINT_ERROR( msgBuf , 1)
149 CALL MODELDATA_EXAMPLE( myThid )
150 STOP 'ABNORMAL END: S/R INI_PARMS'
151 ENDIF
152
153 DO WHILE ( .TRUE. )
154 READ(modelDataUnit,FMT='(A)',END=1001) RECORD
155 IL = MAX(ILNBLNK(RECORD),1)
156 IF ( RECORD(1:1) .NE. commentCharacter )
157 & WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL)
158 WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL)
159 ENDDO
160 1001 CONTINUE
161 CLOSE(modelDataUnit)
162
163 C-- Report contents of model parameter file
164 WRITE(msgBuf,'(A)')
165 &'// ======================================================='
166 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
167 & SQUEEZE_RIGHT , 1)
168 WRITE(msgBuf,'(A)') '// Model parameter file "data"'
169 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
170 & SQUEEZE_RIGHT , 1)
171 WRITE(msgBuf,'(A)')
172 &'// ======================================================='
173 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
174 & SQUEEZE_RIGHT , 1)
175 iUnit = scrUnit2
176 REWIND(iUnit)
177 DO WHILE ( .TRUE. )
178 READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD
179 IL = MAX(ILNBLNK(RECORD),1)
180 WRITE(msgBuf,'(A,A)') '>',RECORD(:IL)
181 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
182 & SQUEEZE_RIGHT , 1)
183 ENDDO
184 2001 CONTINUE
185 CLOSE(iUnit)
186 WRITE(msgBuf,'(A)') ' '
187 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
188 & SQUEEZE_RIGHT , 1)
189
190
191 C-- Read settings from model parameter file "data".
192 iUnit = scrUnit1
193 REWIND(iUnit)
194
195 C-- Set default "physical" parameters
196 DO K =1,Nr
197 tRef(K) = 30.D0 - FLOAT( K )
198 ENDDO
199 gravity = 9.81D0
200 gBaro = gravity
201 rhoNil = 999.8D0
202 rhoConst = 999.8D0
203 f0 = 1.D-4
204 beta = 1.D-11
205 viscAh = 1.D3
206 diffKhT = 1.D3
207 diffKhS = 1.D3
208 viscArDefault = 1.D-3
209 viscAz = UNSET_RL
210 viscAr = UNSET_RL
211 viscAp = UNSET_RL
212 diffKrTDefault = 1.D-5
213 diffKzT = UNSET_RL
214 diffKpT = UNSET_RL
215 diffKrT = UNSET_RL
216 diffKrSDefault = 1.D-5
217 diffKzS = UNSET_RL
218 diffKpS = UNSET_RL
219 diffKrS = UNSET_RL
220 viscA4 = 0.
221 diffK4T = 0.
222 diffK4S = 0.
223 GMmaxslope = 1.D-2
224 GMlength = 200.D3
225 GMalpha = 0.D0
226 GMdepth = 1000.D0
227 GMkbackground= 0.D0
228 GMmaxval = 2500.D0
229 tAlpha = 2.D-4
230 sBeta = 7.4D-4
231 eosType = 'LINEAR'
232 buoyancyRelation = 'OCEANIC'
233 implicitFreeSurface = .TRUE.
234 rigidLid = .FALSE.
235 freeSurfFac = 1.D0
236 hFacMin = 0.D0
237 hFacMinDrDefault = 0.D0
238 hFacMinDr = UNSET_RL
239 hFacMinDz = UNSET_RL
240 hFacMinDp = UNSET_RL
241 momViscosity = .TRUE.
242 momAdvection = .TRUE.
243 momForcing = .TRUE.
244 useCoriolis = .TRUE.
245 momPressureForcing = .TRUE.
246 momStepping = .TRUE.
247 tempStepping = .TRUE.
248 saltStepping = .TRUE.
249 metricTerms = .TRUE.
250 implicitDiffusion = .FALSE.
251 READ(UNIT=iUnit,NML=PARM01,IOSTAT=errIO)
252 IF ( errIO .LT. 0 ) THEN
253 WRITE(msgBuf,'(A)')
254 & 'S/R INI_PARMS'
255 CALL PRINT_ERROR( msgBuf , 1)
256 WRITE(msgBuf,'(A)')
257 & 'Error reading numerical model '
258 CALL PRINT_ERROR( msgBuf , 1)
259 WRITE(msgBuf,'(A)')
260 & 'parameter file "data"'
261 CALL PRINT_ERROR( msgBuf , 1)
262 WRITE(msgBuf,'(A)')
263 & 'Problem in namelist PARM01'
264 CALL PRINT_ERROR( msgBuf , 1)
265 CALL MODELDATA_EXAMPLE( myThid )
266 STOP 'ABNORMAL END: S/R INI_PARMS'
267 ENDIF
268 IF ( implicitFreeSurface ) freeSurfFac = 1.D0
269 IF ( rigidLid ) freeSurfFac = 0.D0
270 C-- Momentum viscosity on/off flag.
271 IF ( momViscosity ) THEN
272 vfFacMom = 1.D0
273 ELSE
274 vfFacMom = 0.D0
275 ENDIF
276 C-- Momentum advection on/off flag.
277 IF ( momAdvection ) THEN
278 afFacMom = 1.D0
279 ELSE
280 afFacMom = 0.D0
281 ENDIF
282 C-- Momentum forcing on/off flag.
283 IF ( momForcing ) THEN
284 foFacMom = 1.D0
285 ELSE
286 foFacMom = 0.D0
287 ENDIF
288 C-- Coriolis term on/off flag.
289 IF ( useCoriolis ) THEN
290 cfFacMom = 1.D0
291 ELSE
292 cfFacMom = 0.D0
293 ENDIF
294 C-- Pressure term on/off flag.
295 IF ( momPressureForcing ) THEN
296 pfFacMom = 1.D0
297 ELSE
298 pfFacMom = 0.D0
299 ENDIF
300 C-- Metric terms on/off flag.
301 IF ( metricTerms ) THEN
302 mTFacMom = 1.D0
303 ELSE
304 mTFacMom = 1.D0
305 ENDIF
306 C-- z,p,r coord input switching.
307 IF ( viscAz .NE. UNSET_RL ) zCoordInputData = .TRUE.
308 IF ( viscAp .NE. UNSET_RL ) pCoordInputData = .TRUE.
309 IF ( viscAr .NE. UNSET_RL ) rCoordInputData = .TRUE.
310 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAz
311 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAp
312 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscArDefault
313
314 IF ( diffKzT .NE. UNSET_RL ) zCoordInputData = .TRUE.
315 IF ( diffKpT .NE. UNSET_RL ) pCoordInputData = .TRUE.
316 IF ( diffKrT .NE. UNSET_RL ) rCoordInputData = .TRUE.
317 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKzT
318 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKpT
319 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKrTDefault
320
321 IF ( diffKzS .NE. UNSET_RL ) zCoordInputData = .TRUE.
322 IF ( diffKpS .NE. UNSET_RL ) pCoordInputData = .TRUE.
323 IF ( diffKrS .NE. UNSET_RL ) rCoordInputData = .TRUE.
324 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKzS
325 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKpS
326 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKrSDefault
327
328 IF ( hFacMinDz .NE. UNSET_RL ) zCoordInputData = .TRUE.
329 IF ( hFacMinDp .NE. UNSET_RL ) pCoordInputData = .TRUE.
330 IF ( hFacMinDr .NE. UNSET_RL ) rCoordInputData = .TRUE.
331 IF ( hFacMinDz .EQ. UNSET_RL ) hFacMinDr = hFacMinDz
332 IF ( hFacMinDp .EQ. UNSET_RL ) hFacMinDr = hFacMinDp
333 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDrDefault
334
335 IF ( implicitFreeSurface .AND. rigidLid ) THEN
336 WRITE(msgBuf,'(A,A)')
337 & 'S/R INI_PARMS: Cannot select both implicitFreeSurface',
338 & ' and rigidLid.'
339 CALL PRINT_ERROR( msgBuf , myThid)
340 STOP 'ABNORMAL END: S/R INI_PARMS'
341 ENDIF
342 coordsSet = 0
343 IF ( zCoordInputData ) coordsSet = coordsSet + 1
344 IF ( pCoordInputData ) coordsSet = coordsSet + 1
345 IF ( rCoordInputData ) coordsSet = coordsSet + 1
346 IF ( coordsSet .GT. 1 ) THEN
347 WRITE(msgBuf,'(A)')
348 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
349 CALL PRINT_ERROR( msgBuf , myThid)
350 STOP 'ABNORMAL END: S/R INI_PARMS'
351 ENDIF
352 IF ( rhoConst .LE. 0. ) THEN
353 WRITE(msgBuf,'(A)')
354 & 'S/R INI_PARMS: rhoConst must be greater than 0.'
355 CALL PRINT_ERROR( msgBuf , myThid)
356 STOP 'ABNORMAL END: S/R INI_PARMS'
357 ELSE
358 recip_rhoConst = 1.D0 / rhoConst
359 ENDIF
360 IF ( gravity .LE. 0. ) THEN
361 WRITE(msgBuf,'(A)')
362 & 'S/R INI_PARMS: gravity must be greater than 0.'
363 CALL PRINT_ERROR( msgBuf , myThid)
364 STOP 'ABNORMAL END: S/R INI_PARMS'
365 ELSE
366 recip_gravity = 1.D0 / gravity
367 ENDIF
368
369 C-- Elliptic solver parameters
370 cg2dMaxIters = 150
371 cg2dTargetResidual = 1.D-7
372 cg2dChkResFreq = 1
373 cg3dMaxIters = 150
374 cg3dTargetResidual = 1.D-7
375 cg3dChkResFreq = 1
376 cg2dpcOffDFac = 0.51D0
377 READ(UNIT=iUnit,NML=PARM02,IOSTAT=errIO)
378 IF ( errIO .LT. 0 ) THEN
379 WRITE(msgBuf,'(A)')
380 & 'S/R INI_PARMS'
381 CALL PRINT_ERROR( msgBuf , 1)
382 WRITE(msgBuf,'(A)')
383 & 'Error reading numerical model '
384 CALL PRINT_ERROR( msgBuf , 1)
385 WRITE(msgBuf,'(A)')
386 & 'parameter file "data".'
387 CALL PRINT_ERROR( msgBuf , 1)
388 WRITE(msgBuf,'(A)')
389 & 'Problem in namelist PARM02'
390 CALL PRINT_ERROR( msgBuf , 1)
391 CALL MODELDATA_EXAMPLE( myThid )
392 STOP 'ABNORMAL END: S/R INI_PARMS'
393 ENDIF
394
395 C-- Time stepping parameters
396 startTime = 0.
397 nTimeSteps = 0
398 endTime = 0.
399 nIter0 = 0
400 deltaT = 0.
401 deltaTClock = 0.
402 deltaTtracer = 0.
403 deltaTMom = 0.
404 abEps = 0.01
405 pchkPtFreq = 0.
406 chkPtFreq = 3600.*25
407 dumpFreq = 3600.*100
408 taveFreq = 0.
409 writeStatePrec = precFloat64
410 nCheckLev = 1
411 checkPtSuff(1) = 'ckptA'
412 checkPtSuff(2) = 'ckptB'
413 cAdjFreq = -1.D0
414 rCD = -1.D0
415 tauCD = 0.D0
416 tauThetaClimRelax = 0.D0
417 doThetaClimRelax = .FALSE.
418 tauSaltClimRelax = 0.D0
419 doSaltClimRelax = .FALSE.
420 periodicExternalForcing = .FALSE.
421 externForcingPeriod = 0.
422 externForcingCycle = 0.
423 READ(UNIT=iUnit,NML=PARM03,IOSTAT=errIO)
424 IF ( errIO .LT. 0 ) THEN
425 WRITE(msgBuf,'(A)')
426 & 'S/R INI_PARMS'
427 CALL PRINT_ERROR( msgBuf , 1)
428 WRITE(msgBuf,'(A)')
429 & 'Error reading numerical model '
430 CALL PRINT_ERROR( msgBuf , 1)
431 WRITE(msgBuf,'(A)')
432 & 'parameter file "data"'
433 CALL PRINT_ERROR( msgBuf , 1)
434 WRITE(msgBuf,'(A)')
435 & 'Problem in namelist PARM03'
436 CALL PRINT_ERROR( msgBuf , 1)
437 CALL MODELDATA_EXAMPLE( myThid )
438 STOP 'ABNORMAL END: S/R INI_PARMS'
439 ENDIF
440 C Process "timestepping" params
441 C o Time step size
442 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
443 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
444 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
445 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
446 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
447 IF ( periodicExternalForcing ) THEN
448 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
449 WRITE(msgBuf,'(A)')
450 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
451 CALL PRINT_ERROR( msgBuf , 1)
452 STOP 'ABNORMAL END: S/R INI_PARMS'
453 ENDIF
454 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
455 & externForcingCycle/externForcingPeriod ) THEN
456 WRITE(msgBuf,'(A)')
457 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
458 CALL PRINT_ERROR( msgBuf , 1)
459 STOP 'ABNORMAL END: S/R INI_PARMS'
460 ENDIF
461 IF ( externForcingCycle.le.externForcingPeriod ) THEN
462 WRITE(msgBuf,'(A)')
463 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
464 CALL PRINT_ERROR( msgBuf , 1)
465 STOP 'ABNORMAL END: S/R INI_PARMS'
466 ENDIF
467 IF ( externForcingPeriod.lt.deltaTclock ) THEN
468 WRITE(msgBuf,'(A)')
469 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
470 CALL PRINT_ERROR( msgBuf , 1)
471 STOP 'ABNORMAL END: S/R INI_PARMS'
472 ENDIF
473 ENDIF
474 C o Convection frequency
475 IF ( cAdjFreq .LT. 0. ) THEN
476 cAdjFreq = deltaTClock
477 ENDIF
478 C o CD coupling
479 IF ( tauCD .EQ. 0.D0 ) THEN
480 tauCD = deltaTmom
481 ENDIF
482 IF ( rCD .LT. 0. ) THEN
483 rCD = 1. - deltaTMom/tauCD
484 ENDIF
485 C o Temperature climatology relaxation time scale
486 IF ( tauThetaClimRelax .EQ. 0.D0 ) THEN
487 doThetaClimRelax = .FALSE.
488 lambdaThetaClimRelax = 0.D0
489 ELSE
490 doThetaClimRelax = .TRUE.
491 lambdaThetaClimRelax = 1./tauThetaClimRelax
492 ENDIF
493 C o Salinity climatology relaxation time scale
494 IF ( tauSaltClimRelax .EQ. 0.D0 ) THEN
495 doSaltClimRelax = .FALSE.
496 lambdaSaltClimRelax = 0.D0
497 ELSE
498 doSaltClimRelax = .TRUE.
499 lambdaSaltClimRelax = 1./tauSaltClimRelax
500 ENDIF
501 C o Time step count
502 IF ( endTime .NE. 0 ) THEN
503 IF ( deltaTClock .NE. 0 ) nTimeSteps =
504 & INT((endTime-startTime)/deltaTClock)
505 ENDIF
506 C o Finish time
507 IF ( endTime .EQ. 0. ) endTime = FLOAT(nTimeSteps)*deltaTClock
508
509 C o If taveFreq is finite, then we must make sure the diagnostics
510 C code is being compiled
511 #ifndef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
512 IF (taveFreq.NE.0.) THEN
513 WRITE(msgBuf,'(A)')
514 & 'S/R INI_PARMS: taveFreq <> 0 but you have'
515 CALL PRINT_ERROR( msgBuf , 1)
516 WRITE(msgBuf,'(A)')
517 & 'not compiled the model with the diagnostics routines.'
518 CALL PRINT_ERROR( msgBuf , 1)
519 WRITE(msgBuf,'(A,A)')
520 & 'Re-compile with: #define INCLUDE_DIAGNOSTICS_INTERFACE_CODE',
521 & ' or -DINCLUDE_DIAGNOSTICS_INTERFACE_CODE'
522 CALL PRINT_ERROR( msgBuf , 1)
523 STOP 'ABNORMAL END: S/R INI_PARMS'
524 ENDIF
525 #endif
526
527 C-- Grid parameters
528 C In cartesian coords distances are in metres
529 rkFac = UNSET_I
530 usingCartesianGrid = .TRUE.
531 delRDefault = 1.D2
532 DO K =1,Nr
533 delZ(K) = UNSET_RL
534 delP(K) = UNSET_RL
535 delR(K) = UNSET_RL
536 ENDDO
537 dxSpacing = 20.D0 * 1000.D0
538 dySpacing = 20.D0 * 1000.D0
539 DO i=1,Nx
540 delX(i) = dxSpacing
541 ENDDO
542 DO j=1,Ny
543 delY(j) = dySpacing
544 ENDDO
545 C In spherical polar distances are in degrees
546 usingSphericalPolarGrid = .FALSE.
547 phiMin = -5.0
548 thetaMin = 0.
549 rSphere = 6370. * 1.D3
550 recip_rSphere = 1.D0/rSphere
551 IF ( usingSphericalPolarGrid ) THEN
552 dxSpacing = 1.
553 dySpacing = 1.
554 DO I=1,Nx
555 delX(I) = dxSpacing
556 ENDDO
557 DO J=1,Ny
558 delY(J) = dySpacing
559 ENDDO
560 ENDIF
561
562 READ(UNIT=iUnit,NML=PARM04,IOSTAT=errIO)
563 IF ( errIO .LT. 0 ) THEN
564 WRITE(msgBuf,'(A)')
565 & 'S/R INI_PARMS'
566 CALL PRINT_ERROR( msgBuf , 1)
567 WRITE(msgBuf,'(A)')
568 & 'Error reading numerical model '
569 CALL PRINT_ERROR( msgBuf , 1)
570 WRITE(msgBuf,'(A)')
571 & 'parameter file "data"'
572 CALL PRINT_ERROR( msgBuf , 1)
573 WRITE(msgBuf,'(A)')
574 & 'Problem in namelist PARM04'
575 CALL PRINT_ERROR( msgBuf , 1)
576 CALL MODELDATA_EXAMPLE( myThid )
577 STOP 'ABNORMAL END: S/R INI_PARMS'
578 ENDIF
579 C
580 IF ( rSphere .NE. 0 ) THEN
581 recip_rSphere = 1.D0/rSphere
582 ELSE
583 recip_rSphere = 0.
584 ENDIF
585 C-- Initialize EOS coefficients (3rd order polynomial)
586 IF (eostype.eq.'POLY3') THEN
587 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
588 READ(37,*) I
589 IF (I.NE.Nr) THEN
590 WRITE(msgBuf,'(A)')
591 & 'ini_parms: attempt to read POLY3.COEFFS failed'
592 CALL PRINT_ERROR( msgBuf , 1)
593 WRITE(msgBuf,'(A)')
594 & ' because bad # of levels in data'
595 CALL PRINT_ERROR( msgBuf , 1)
596 STOP 'Bad data in POLY3.COEFFS'
597 ENDIF
598 READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nr)
599 DO K=1,Nr
600 READ(37,*) (eosC(I,K),I=1,9)
601 ENDDO
602 CLOSE(37)
603 ENDIF
604 C-- Check for conflicting grid definitions.
605 goptCount = 0
606 IF ( usingCartesianGrid ) goptCount = goptCount+1
607 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
608 IF ( goptCount .NE. 1 ) THEN
609 WRITE(msgBuf,'(A)')
610 & 'S/R INI_PARMS: More than one coordinate system requested'
611 CALL PRINT_ERROR( msgBuf , myThid)
612 STOP 'ABNORMAL END: S/R INI_PARMS'
613 ENDIF
614 C-- Make metric term settings consistent with underlying grid.
615 IF ( usingCartesianGrid ) THEN
616 usingSphericalPolarMterms = .FALSE.
617 metricTerms = .FALSE.
618 mTFacMom = 0
619 useBetaPlaneF = .TRUE.
620 ENDIF
621 IF ( usingSphericalPolarGrid ) THEN
622 useConstantF = .FALSE.
623 useBetaPlaneF = .FALSE.
624 useSphereF = .TRUE.
625 omega = 2.D0 * PI / ( 3600.D0 * 24.D0 )
626 usingSphericalPolarMterms = .TRUE.
627 metricTerms = .TRUE.
628 mTFacMom = 1
629 ENDIF
630 C-- p, z, r coord parameters
631 DO K = 1, Nr
632 IF ( delZ(K) .NE. UNSET_RL ) zCoordInputData = .TRUE.
633 IF ( delP(K) .NE. UNSET_RL ) pCoordInputData = .TRUE.
634 IF ( delR(K) .NE. UNSET_RL ) rCoordInputData = .TRUE.
635 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delZ(K)
636 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delP(K)
637 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delRDefault
638 ENDDO
639 C Check for multiple coordinate systems
640 coordsSet = 0
641 IF ( zCoordInputData ) coordsSet = coordsSet + 1
642 IF ( pCoordInputData ) coordsSet = coordsSet + 1
643 IF ( rCoordInputData ) coordsSet = coordsSet + 1
644 IF ( coordsSet .GT. 1 ) THEN
645 WRITE(msgBuf,'(A)')
646 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
647 CALL PRINT_ERROR( msgBuf , myThid)
648 STOP 'ABNORMAL END: S/R INI_PARMS'
649 ENDIF
650
651 C-- Input files
652 bathyFile = ' '
653 hydrogSaltFile = ' '
654 hydrogThetaFile = ' '
655 zonalWindFile = ' '
656 meridWindFile = ' '
657 thetaClimFile = ' '
658 saltClimFile = ' '
659 READ(UNIT=iUnit,NML=PARM05,IOSTAT=errIO)
660 IF ( errIO .LT. 0 ) THEN
661 WRITE(msgBuf,'(A)')
662 & 'S/R INI_PARMS'
663 CALL PRINT_ERROR( msgBuf , 1)
664 WRITE(msgBuf,'(A)')
665 & 'Error reading numerical model '
666 CALL PRINT_ERROR( msgBuf , 1)
667 WRITE(msgBuf,'(A)')
668 & 'parameter file "data"'
669 CALL PRINT_ERROR( msgBuf , 1)
670 WRITE(msgBuf,'(A)')
671 & 'Problem in namelist PARM05'
672 CALL PRINT_ERROR( msgBuf , 1)
673 CALL MODELDATA_EXAMPLE( myThid )
674 STOP 'ABNORMAL END: S/R INI_PARMS'
675 ENDIF
676
677 C
678 C-- Set factors required for mixing pressure and meters as vertical coordinate.
679 C rkFac is a "sign" parameter which is used where the orientation of the vertical
680 C coordinate (pressure or meters) relative to the vertical index (K) is important.
681 C rkFac = 1 applies when K and the coordinate are in the opposite sense.
682 C rkFac = -1 applies when K and the coordinate are in the same sense.
683 C horiVertRatio is a parameter that maps horizontal units to vertical units.
684 C It is used in certain special cases where lateral and vertical terms are
685 C being combined and a single frame of reference is needed.
686 IF ( zCoordInputData .AND. rkFac .EQ. UNSET_I ) THEN
687 rkFac = 1.D0
688 horiVertRatio = 1.D0
689 ENDIF
690 IF ( pCoordInputData .AND. rkFac .EQ. UNSET_I ) THEN
691 rkFac = -1.D0
692 horiVertRatio = Gravity * rhoConst
693 ENDIF
694 IF ( rCoordInputData .AND. rkFac .EQ. UNSET_I ) THEN
695 rkFac = 1.D0
696 horiVertRatio = 1.D0
697 ENDIF
698 recip_rkFac = 1.D0 / rkFac
699 recip_horiVertRatio = 1./horiVertRatio
700 IF ( zCoordInputData ) usingZCoords = .TRUE.
701 IF ( pCoordInputData ) usingPCoords = .TRUE.
702 C
703 CLOSE(iUnit)
704
705 _END_MASTER(myThid)
706
707 C-- Everyone else must wait for the parameters to be loaded
708 _BARRIER
709 C
710
711 RETURN
712 END
713

  ViewVC Help
Powered by ViewVC 1.1.22