/[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.38 - (show annotations) (download)
Wed Dec 9 16:11:52 1998 UTC (25 years, 5 months ago) by adcroft
Branch: MAIN
Changes since 1.37: +10 -1 lines
Added IMPLICIT NONE in a lot of subroutines.
Also corrected the recip_Rhonil bug: we didn't set it in ini_parms.F

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

  ViewVC Help
Powered by ViewVC 1.1.22