/[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.61 - (show annotations) (download)
Wed Jun 6 14:55:45 2001 UTC (22 years, 11 months ago) by adcroft
Branch: MAIN
Changes since 1.60: +2 -2 lines
Added a debugMode that uses same statistics stuff as monitor.F
Can be disabled with -DEXCLUDE_DEBUGMODE. Turn on at run-time
with debugMode=.true.  Default is enabled but off.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_parms.F,v 1.60 2001/05/30 19:33:18 adcroft Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 SUBROUTINE INI_PARMS( myThid )
7 C /==========================================================\
8 C | SUBROUTINE INI_PARMS |
9 C | o Routine to set model "parameters" |
10 C |==========================================================|
11 C | Notes: |
12 C | ====== |
13 C | The present version of this routine is a place-holder. |
14 C | A production version needs to handle parameters from an |
15 C | external file and possibly reading in some initial field |
16 C | values. |
17 C \==========================================================/
18 IMPLICIT NONE
19
20 C === Global variables ===
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GRID.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_FNAM) delXfile
46 CHARACTER*(MAX_LEN_FNAM) delYfile
47 CHARACTER*(MAX_LEN_MBUF) msgBuf
48 CHARACTER*(MAX_LEN_PREC) record
49 INTEGER goptCount
50 INTEGER K, I, J, IL, iUnit
51 INTEGER errIO
52 INTEGER IFNBLNK
53 EXTERNAL IFNBLNK
54 INTEGER ILNBLNK
55 EXTERNAL ILNBLNK
56 C Default values for variables which have vertical coordinate system
57 C dependency.
58 _RL viscArDefault
59 _RL diffKrTDefault
60 _RL diffKrSDefault
61 _RL hFacMinDrDefault
62 _RL delRDefault(Nr)
63 _RS rkFacDefault
64 C zCoordInputData - These are used to select between different coordinate systems.
65 C pCoordInputData The vertical coordinate system in the rest of the model is
66 C rCoordInputData written in terms of r. In the model "data" file input data can
67 C coordsSet be interms of z, p or r.
68 C e.g. delZ or delP or delR for the vertical grid spacing.
69 C The following rules apply:
70 C All parameters must use the same vertical coordinate system.
71 C e.g. delZ and viscAz is legal but
72 C delZ and viscAr is an error.
73 C Similarly specifyinh delZ and delP is an error.
74 C zCoord..., pCoord..., rCoord... are used to flag when z, p or r are
75 C used. coordsSet counts how many vertical coordinate systems have been
76 C used to specify variables. coordsSet > 1 is an error.
77 C
78 LOGICAL zCoordInputData
79 LOGICAL pCoordInputData
80 LOGICAL rCoordInputData
81 INTEGER coordsSet
82
83 C-- Continuous equation parameters
84 NAMELIST /PARM01/
85 & gravity, gBaro, rhonil, tAlpha, sBeta, f0, beta, omega,
86 & viscAh, viscAz, viscA4, cosPower,
87 & diffKhT, diffKzT, diffK4T,
88 & diffKhS, diffKzS, diffK4S,
89 & tRef, sRef, eosType,
90 & no_slip_sides,no_slip_bottom,
91 & momViscosity, momAdvection, momForcing, useCoriolis,
92 & momPressureForcing, metricTerms,
93 & tempDiffusion, tempAdvection, tempForcing,
94 & saltDiffusion, saltAdvection, saltForcing,
95 & implicSurfPress, implicDiv2DFlow,
96 & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz,
97 & exactConserv, nonlinFreeSurf, hFacInf, hFacSup,
98 & staggerTimeStep,
99 & tempStepping, saltStepping, momStepping,
100 & implicitDiffusion, implicitViscosity,
101 & viscAr, diffKrT, diffKrS, hFacMinDr,
102 & viscAp, diffKpT, diffKpS, hFacMinDp,
103 & rhoConst, buoyancyRelation, HeatCapacity_Cp,
104 & writeBinaryPrec, readBinaryPrec, writeStatePrec,
105 & nonHydrostatic, globalFiles,
106 & allowFreezing, ivdc_kappa,
107 & zonal_filt_lat, zonal_filt_sinpow, zonal_filt_cospow,
108 & bottomDragLinear,bottomDragQuadratic,
109 & usePickupBeforeC35, debugMode
110
111 C-- Elliptic solver parameters
112 NAMELIST /PARM02/
113 & cg2dMaxIters, cg2dChkResFreq, cg2dTargetResidual,
114 & cg2dTargetResWunit, cg2dpcOffDFac,
115 & cg3dMaxIters, cg3dChkResFreq, cg3dTargetResidual
116
117 C-- Time stepping parammeters
118 NAMELIST /PARM03/
119 & nIter0, nTimeSteps, nEndIter, deltaT, deltaTmom, deltaTtracer,
120 & abEps, tauCD, rCD,
121 & startTime, endTime, chkPtFreq, dumpFreq, taveFreq, deltaTClock,
122 & monitorFreq,
123 & pChkPtFreq, cAdjFreq, tauThetaClimRelax, tauSaltClimRelax,
124 & periodicExternalForcing, externForcingPeriod, externForcingCycle
125
126 C-- Gridding parameters
127 NAMELIST /PARM04/
128 & usingCartesianGrid, dxSpacing, dySpacing, delX, delY, delZ,
129 & usingSphericalPolarGrid, phiMin, thetaMin, rSphere,
130 & usingCurvilinearGrid,
131 & delP, delR, rkFac, Ro_SeaLevel, groundAtK1,
132 & delXfile, delYfile
133
134 C-- Input files
135 NAMELIST /PARM05/
136 & bathyFile, hydrogThetaFile, hydrogSaltFile,
137 & zonalWindFile, meridWindFile,
138 & thetaClimFile, saltClimFile,
139 & surfQfile, EmPmRfile, surfQswfile,
140 & uVelInitFile, vVelInitFile, pSurfInitFile,
141 & dQdTFile
142
143 C
144 _BEGIN_MASTER(myThid)
145
146 C Defaults values for input parameters
147 CALL SET_DEFAULTS(
148 O viscArDefault, diffKrTDefault, diffKrSDefault,
149 O hFacMinDrDefault, delRdefault, rkFacDefault,
150 I myThid )
151
152 C-- Initialise "which vertical coordinate system used" flags.
153 zCoordInputData = .FALSE.
154 pCoordInputData = .FALSE.
155 rCoordInputData = .FALSE.
156 usingPCoords = .FALSE.
157 usingZCoords = .FALSE.
158 coordsSet = 0
159
160 C-- Open the parameter file
161 OPEN(UNIT=scrUnit1,STATUS='SCRATCH')
162 OPEN(UNIT=scrUnit2,STATUS='SCRATCH')
163 OPEN(UNIT=modelDataUnit,FILE='data',STATUS='OLD',
164 & IOSTAT=errIO)
165 IF ( errIO .LT. 0 ) THEN
166 WRITE(msgBuf,'(A)')
167 & 'S/R INI_PARMS'
168 CALL PRINT_ERROR( msgBuf , 1)
169 WRITE(msgBuf,'(A)')
170 & 'Unable to open model parameter'
171 CALL PRINT_ERROR( msgBuf , 1)
172 WRITE(msgBuf,'(A)')
173 & 'file "data"'
174 CALL PRINT_ERROR( msgBuf , 1)
175 CALL MODELDATA_EXAMPLE( myThid )
176 STOP 'ABNORMAL END: S/R INI_PARMS'
177 ENDIF
178
179 DO WHILE ( .TRUE. )
180 READ(modelDataUnit,FMT='(A)',END=1001) RECORD
181 IL = MAX(ILNBLNK(RECORD),1)
182 IF ( RECORD(1:1) .NE. commentCharacter )
183 & WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL)
184 WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL)
185 ENDDO
186 1001 CONTINUE
187 CLOSE(modelDataUnit)
188
189 C-- Report contents of model parameter file
190 WRITE(msgBuf,'(A)')
191 &'// ======================================================='
192 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
193 & SQUEEZE_RIGHT , 1)
194 WRITE(msgBuf,'(A)') '// Model parameter file "data"'
195 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
196 & SQUEEZE_RIGHT , 1)
197 WRITE(msgBuf,'(A)')
198 &'// ======================================================='
199 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
200 & SQUEEZE_RIGHT , 1)
201 iUnit = scrUnit2
202 REWIND(iUnit)
203 DO WHILE ( .TRUE. )
204 READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD
205 IL = MAX(ILNBLNK(RECORD),1)
206 WRITE(msgBuf,'(A,A)') '>',RECORD(:IL)
207 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
208 & SQUEEZE_RIGHT , 1)
209 ENDDO
210 2001 CONTINUE
211 CLOSE(iUnit)
212 WRITE(msgBuf,'(A)') ' '
213 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
214 & SQUEEZE_RIGHT , 1)
215
216
217 C-- Read settings from model parameter file "data".
218 iUnit = scrUnit1
219 REWIND(iUnit)
220
221 C-- Set default "physical" parameters
222 viscAz = UNSET_RL
223 viscAr = UNSET_RL
224 viscAp = UNSET_RL
225 diffKzT = UNSET_RL
226 diffKpT = UNSET_RL
227 diffKrT = UNSET_RL
228 diffKzS = UNSET_RL
229 diffKpS = UNSET_RL
230 diffKrS = UNSET_RL
231 gBaro = UNSET_RL
232 rhoConst = UNSET_RL
233 hFacMinDr = UNSET_RL
234 hFacMinDz = UNSET_RL
235 hFacMinDp = UNSET_RL
236 READ(UNIT=iUnit,NML=PARM01) !,IOSTAT=errIO)
237 IF ( errIO .LT. 0 ) THEN
238 WRITE(msgBuf,'(A)')
239 & 'S/R INI_PARMS'
240 CALL PRINT_ERROR( msgBuf , 1)
241 WRITE(msgBuf,'(A)')
242 & 'Error reading numerical model '
243 CALL PRINT_ERROR( msgBuf , 1)
244 WRITE(msgBuf,'(A)')
245 & 'parameter file "data"'
246 CALL PRINT_ERROR( msgBuf , 1)
247 WRITE(msgBuf,'(A)')
248 & 'Problem in namelist PARM01'
249 CALL PRINT_ERROR( msgBuf , 1)
250 CALL MODELDATA_EXAMPLE( myThid )
251 STOP 'ABNORMAL END: S/R INI_PARMS'
252 ENDIF
253 IF ( implicitFreeSurface ) freeSurfFac = 1.D0
254 IF ( rigidLid ) freeSurfFac = 0.D0
255 IF ( gBaro .EQ. UNSET_RL ) gBaro=gravity
256 IF ( rhoConst .EQ. UNSET_RL ) rhoConst=rhoNil
257 C-- Momentum viscosity on/off flag.
258 IF ( momViscosity ) THEN
259 vfFacMom = 1.D0
260 ELSE
261 vfFacMom = 0.D0
262 ENDIF
263 C-- Momentum advection on/off flag.
264 IF ( momAdvection ) THEN
265 afFacMom = 1.D0
266 ELSE
267 afFacMom = 0.D0
268 ENDIF
269 C-- Momentum forcing on/off flag.
270 IF ( momForcing ) THEN
271 foFacMom = 1.D0
272 ELSE
273 foFacMom = 0.D0
274 ENDIF
275 C-- Coriolis term on/off flag.
276 IF ( useCoriolis ) THEN
277 cfFacMom = 1.D0
278 ELSE
279 cfFacMom = 0.D0
280 ENDIF
281 C-- Pressure term on/off flag.
282 IF ( momPressureForcing ) THEN
283 pfFacMom = 1.D0
284 ELSE
285 pfFacMom = 0.D0
286 ENDIF
287 C-- Metric terms on/off flag.
288 IF ( metricTerms ) THEN
289 mTFacMom = 1.D0
290 ELSE
291 mTFacMom = 0.D0
292 ENDIF
293 C-- z,p,r coord input switching.
294 IF ( viscAz .NE. UNSET_RL ) zCoordInputData = .TRUE.
295 IF ( viscAp .NE. UNSET_RL ) pCoordInputData = .TRUE.
296 IF ( viscAr .NE. UNSET_RL ) rCoordInputData = .TRUE.
297 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAz
298 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAp
299 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscArDefault
300
301 IF ( diffKzT .NE. UNSET_RL ) zCoordInputData = .TRUE.
302 IF ( diffKpT .NE. UNSET_RL ) pCoordInputData = .TRUE.
303 IF ( diffKrT .NE. UNSET_RL ) rCoordInputData = .TRUE.
304 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKzT
305 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKpT
306 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKrTDefault
307
308 IF ( diffKzS .NE. UNSET_RL ) zCoordInputData = .TRUE.
309 IF ( diffKpS .NE. UNSET_RL ) pCoordInputData = .TRUE.
310 IF ( diffKrS .NE. UNSET_RL ) rCoordInputData = .TRUE.
311 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKzS
312 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKpS
313 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKrSDefault
314
315 IF ( hFacMinDz .NE. UNSET_RL ) zCoordInputData = .TRUE.
316 IF ( hFacMinDp .NE. UNSET_RL ) pCoordInputData = .TRUE.
317 IF ( hFacMinDr .NE. UNSET_RL ) rCoordInputData = .TRUE.
318 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDz
319 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDp
320 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDrDefault
321
322 IF ( ivdc_kappa .NE. 0. .AND. .NOT. implicitDiffusion ) THEN
323 WRITE(msgBuf,'(A,A)')
324 & 'S/R INI_PARMS: To use ivdc_kappa you must enable implicit',
325 & ' vertical diffusion.'
326 CALL PRINT_ERROR( msgBuf , myThid)
327 STOP 'ABNORMAL END: S/R INI_PARMS'
328 ENDIF
329
330 IF ( implicitFreeSurface .AND. rigidLid ) THEN
331 WRITE(msgBuf,'(A,A)')
332 & 'S/R INI_PARMS: Cannot select both implicitFreeSurface',
333 & ' and rigidLid.'
334 CALL PRINT_ERROR( msgBuf , myThid)
335 STOP 'ABNORMAL END: S/R INI_PARMS'
336 ENDIF
337 IF ( (implicSurfPress.NE.1. .OR. implicDiv2DFlow.NE.1.)
338 & .AND. nonHydrostatic ) THEN
339 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: nonHydrostatic',
340 & ' NOT SAFE with non-fully implicit Barotropic solver'
341 CALL PRINT_ERROR( msgBuf , myThid)
342 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: To by-pass this',
343 & 'STOP, comment this test and re-compile ini_params'
344 CALL PRINT_ERROR( msgBuf , myThid)
345 STOP 'ABNORMAL END: S/R INI_PARMS'
346 ENDIF
347
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 ( HeatCapacity_Cp .LE. 0. ) THEN
375 WRITE(msgBuf,'(A)')
376 & 'S/R INI_PARMS: HeatCapacity_Cp must be greater than 0.'
377 CALL PRINT_ERROR( msgBuf , myThid)
378 STOP 'ABNORMAL END: S/R INI_PARMS'
379 ELSE
380 recip_Cp = 1.D0 / HeatCapacity_Cp
381 ENDIF
382 IF ( gravity .LE. 0. ) THEN
383 WRITE(msgBuf,'(A)')
384 & 'S/R INI_PARMS: gravity must be greater than 0.'
385 CALL PRINT_ERROR( msgBuf , myThid)
386 STOP 'ABNORMAL END: S/R INI_PARMS'
387 ELSE
388 recip_gravity = 1.D0 / gravity
389 ENDIF
390 C Set globalFiles flag for READ_WRITE_FLD package
391 CALL SET_WRITE_GLOBAL_FLD( globalFiles )
392 C Set globalFiles flag for READ_WRITE_REC package
393 CALL SET_WRITE_GLOBAL_REC( globalFiles )
394 C Set globalFiles flag for READ_WRITE_REC package
395 CALL SET_WRITE_GLOBAL_PICKUP( globalFiles )
396
397 C-- Elliptic solver parameters
398 READ(UNIT=iUnit,NML=PARM02) !,IOSTAT=errIO)
399 IF ( errIO .LT. 0 ) THEN
400 WRITE(msgBuf,'(A)')
401 & 'S/R INI_PARMS'
402 CALL PRINT_ERROR( msgBuf , 1)
403 WRITE(msgBuf,'(A)')
404 & 'Error reading numerical model '
405 CALL PRINT_ERROR( msgBuf , 1)
406 WRITE(msgBuf,'(A)')
407 & 'parameter file "data".'
408 CALL PRINT_ERROR( msgBuf , 1)
409 WRITE(msgBuf,'(A)')
410 & 'Problem in namelist PARM02'
411 CALL PRINT_ERROR( msgBuf , 1)
412 CALL MODELDATA_EXAMPLE( myThid )
413 STOP 'ABNORMAL END: S/R INI_PARMS'
414 ENDIF
415
416 C-- Time stepping parameters
417 rCD = -1.D0
418 READ(UNIT=iUnit,NML=PARM03) !,IOSTAT=errIO)
419 IF ( errIO .LT. 0 ) THEN
420 WRITE(msgBuf,'(A)')
421 & 'S/R INI_PARMS'
422 CALL PRINT_ERROR( msgBuf , 1)
423 WRITE(msgBuf,'(A)')
424 & 'Error reading numerical model '
425 CALL PRINT_ERROR( msgBuf , 1)
426 WRITE(msgBuf,'(A)')
427 & 'parameter file "data"'
428 CALL PRINT_ERROR( msgBuf , 1)
429 WRITE(msgBuf,'(A)')
430 & 'Problem in namelist PARM03'
431 CALL PRINT_ERROR( msgBuf , 1)
432 CALL MODELDATA_EXAMPLE( myThid )
433 STOP 'ABNORMAL END: S/R INI_PARMS'
434 ENDIF
435 C Process "timestepping" params
436 C o Time step size
437 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
438 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
439 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
440 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
441 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
442 IF ( periodicExternalForcing ) THEN
443 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
444 WRITE(msgBuf,'(A)')
445 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
446 CALL PRINT_ERROR( msgBuf , 1)
447 STOP 'ABNORMAL END: S/R INI_PARMS'
448 ENDIF
449 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
450 & externForcingCycle/externForcingPeriod ) THEN
451 WRITE(msgBuf,'(A)')
452 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
453 CALL PRINT_ERROR( msgBuf , 1)
454 STOP 'ABNORMAL END: S/R INI_PARMS'
455 ENDIF
456 IF ( externForcingCycle.le.externForcingPeriod ) THEN
457 WRITE(msgBuf,'(A)')
458 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
459 CALL PRINT_ERROR( msgBuf , 1)
460 STOP 'ABNORMAL END: S/R INI_PARMS'
461 ENDIF
462 IF ( externForcingPeriod.lt.deltaTclock ) THEN
463 WRITE(msgBuf,'(A)')
464 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
465 CALL PRINT_ERROR( msgBuf , 1)
466 STOP 'ABNORMAL END: S/R INI_PARMS'
467 ENDIF
468 ENDIF
469 C o Convection frequency
470 IF ( cAdjFreq .LT. 0. ) THEN
471 cAdjFreq = deltaTClock
472 ENDIF
473 IF ( ivdc_kappa .NE. 0. .AND. cAdjFreq .NE. 0. ) THEN
474 WRITE(msgBuf,'(A,A)')
475 & 'S/R INI_PARMS: You have enabled both ivdc_kappa and',
476 & ' convective adjustment.'
477 CALL PRINT_ERROR( msgBuf , myThid)
478 STOP 'ABNORMAL END: S/R INI_PARMS'
479 ENDIF
480 C o CD coupling
481 IF ( tauCD .EQ. 0.D0 ) THEN
482 tauCD = deltaTmom
483 ENDIF
484 IF ( rCD .LT. 0. ) THEN
485 rCD = 1. - deltaTMom/tauCD
486 ENDIF
487 C o Temperature climatology relaxation time scale
488 IF ( tauThetaClimRelax .EQ. 0.D0 ) THEN
489 doThetaClimRelax = .FALSE.
490 lambdaThetaClimRelax = 0.D0
491 ELSE
492 doThetaClimRelax = .TRUE.
493 lambdaThetaClimRelax = 1./tauThetaClimRelax
494 ENDIF
495 C o Salinity climatology relaxation time scale
496 IF ( tauSaltClimRelax .EQ. 0.D0 ) THEN
497 doSaltClimRelax = .FALSE.
498 lambdaSaltClimRelax = 0.D0
499 ELSE
500 doSaltClimRelax = .TRUE.
501 lambdaSaltClimRelax = 1./tauSaltClimRelax
502 ENDIF
503
504 C o Start time
505 IF ( nIter0 .NE. 0 .AND. startTime .EQ. 0. )
506 & startTime = deltaTClock*float(nIter0)
507 C o nIter0
508 IF ( nIter0 .EQ. 0 .AND. startTime .NE. 0. )
509 & nIter0 = INT( startTime/deltaTClock )
510
511 C o nTimeSteps 1
512 IF ( nTimeSteps .EQ. 0 .AND. nEndIter .NE. 0 )
513 & nTimeSteps = nEndIter-nIter0
514 C o nTimeSteps 2
515 IF ( nTimeSteps .EQ. 0 .AND. endTime .NE. 0. )
516 & nTimeSteps = int(0.5+(endTime-startTime)/deltaTclock)
517 C o nEndIter 1
518 IF ( nEndIter .EQ. 0 .AND. nTimeSteps .NE. 0 )
519 & nEndIter = nIter0+nTimeSteps
520 C o nEndIter 2
521 IF ( nEndIter .EQ. 0 .AND. endTime .NE. 0. )
522 & nEndIter = int(0.5+endTime/deltaTclock)
523 C o End Time 1
524 IF ( endTime .EQ. 0. .AND. nTimeSteps .NE. 0 )
525 & endTime = startTime + deltaTClock*float(nTimeSteps)
526 C o End Time 2
527 IF ( endTime .EQ. 0. .AND. nEndIter .NE. 0 )
528 & endTime = deltaTClock*float(nEndIter)
529
530 C o Consistent?
531 IF ( nEndIter .NE. nIter0+nTimeSteps ) THEN
532 WRITE(msgBuf,'(A)')
533 & 'S/R INI_PARMS: nIter0, nTimeSteps and nEndIter are inconsistent'
534 CALL PRINT_ERROR( msgBuf , 1)
535 WRITE(msgBuf,'(A)')
536 & 'S/R INI_PARMS: Perhaps more than two were set at once'
537 CALL PRINT_ERROR( msgBuf , 1)
538 STOP 'ABNORMAL END: S/R INI_PARMS'
539 ENDIF
540 IF ( nTimeSteps .NE. int(0.5+(endTime-startTime)/deltaTClock) )
541 & THEN
542 WRITE(msgBuf,'(A)')
543 & 'S/R INI_PARMS: both endTime and nTimeSteps have been set'
544 CALL PRINT_ERROR( msgBuf , 1)
545 WRITE(msgBuf,'(A)')
546 & 'S/R INI_PARMS: but are inconsistent'
547 CALL PRINT_ERROR( msgBuf , 1)
548 STOP 'ABNORMAL END: S/R INI_PARMS'
549 ENDIF
550
551 C o Monitor (should also add CPP flag for monitor?)
552 IF (monitorFreq.LT.0.) THEN
553 monitorFreq=0.
554 c IF (dumpFreq.NE.0.) monitorFreq=dumpFreq
555 c IF (taveFreq.NE.0..AND.taveFreq.LT.monitorFreq)
556 c & monitorFreq=taveFreq
557 c IF (chkPtFreq.NE.0..AND.chkPtFreq.LT.monitorFreq)
558 c & monitorFreq=chkPtFreq
559 c IF (pChkPtFreq.NE.0..AND.pChkPtFreq.LT.monitorFreq)
560 c & monitorFreq=pChkPtFreq
561 IF (monitorFreq.EQ.0.) monitorFreq=deltaTclock
562 ENDIF
563
564 C o If taveFreq is finite, then we must make sure the diagnostics
565 C code is being compiled
566 #ifndef ALLOW_TIMEAVE
567 IF (taveFreq.NE.0.) THEN
568 WRITE(msgBuf,'(A)')
569 & 'S/R INI_PARMS: taveFreq <> 0 but you have'
570 CALL PRINT_ERROR( msgBuf , 1)
571 WRITE(msgBuf,'(A)')
572 & 'not compiled the model with the diagnostics routines.'
573 CALL PRINT_ERROR( msgBuf , 1)
574 WRITE(msgBuf,'(A,A)')
575 & 'Re-compile with: #define ALLOW_TIMEAVE',
576 & ' or -DALLOW_TIMEAVE'
577 CALL PRINT_ERROR( msgBuf , 1)
578 STOP 'ABNORMAL END: S/R INI_PARMS'
579 ENDIF
580 #endif
581
582 C-- Grid parameters
583 C In cartesian coords distances are in metres
584 rkFac = UNSET_RS
585 DO K =1,Nr
586 delZ(K) = UNSET_RL
587 delP(K) = UNSET_RL
588 delR(K) = UNSET_RL
589 ENDDO
590 C In spherical polar distances are in degrees
591 recip_rSphere = 1.D0/rSphere
592 dxSpacing = UNSET_RL
593 dySpacing = UNSET_RL
594 delXfile = ' '
595 delYfile = ' '
596 READ(UNIT=iUnit,NML=PARM04) !,IOSTAT=errIO)
597 IF ( errIO .LT. 0 ) THEN
598 WRITE(msgBuf,'(A)')
599 & 'S/R INI_PARMS'
600 CALL PRINT_ERROR( msgBuf , 1)
601 WRITE(msgBuf,'(A)')
602 & 'Error reading numerical model '
603 CALL PRINT_ERROR( msgBuf , 1)
604 WRITE(msgBuf,'(A)')
605 & 'parameter file "data"'
606 CALL PRINT_ERROR( msgBuf , 1)
607 WRITE(msgBuf,'(A)')
608 & 'Problem in namelist PARM04'
609 CALL PRINT_ERROR( msgBuf , 1)
610 CALL MODELDATA_EXAMPLE( myThid )
611 STOP 'ABNORMAL END: S/R INI_PARMS'
612 ENDIF
613
614 C X coordinate
615 IF ( delXfile .NE. ' ' ) THEN
616 IF ( delX(1) .NE. UNSET_RL .OR. dxSpacing .NE. UNSET_RL ) THEN
617 WRITE(msgBuf,'(A,A)') 'Too many specifications for delX:',
618 & 'Specify only one of delX, dxSpacing or delXfile'
619 CALL PRINT_ERROR( msgBuf , 1)
620 STOP 'ABNORMAL END: S/R INI_PARMS'
621 ELSE
622 _BEGIN_MASTER( myThid )
623 IF (readBinaryPrec.EQ.precFloat32) THEN
624 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
625 & ACCESS='DIRECT',RECL=WORDLENGTH*Nx)
626 READ(37,rec=1) delX
627 #ifdef _BYTESWAPIO
628 call MDS_BYTESWAPR4( Nx, delX )
629 #endif
630 CLOSE(37)
631 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
632 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
633 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Nx)
634 READ(37,rec=1) delX
635 #ifdef _BYTESWAPIO
636 call MDS_BYTESWAPR8( Nx, delX )
637 #endif
638 CLOSE(37)
639 ENDIF
640 _END_MASTER(myThid)
641 ENDIF
642 ENDIF
643 IF ( dxSpacing .NE. UNSET_RL ) THEN
644 DO i=1,Nx
645 delX(i) = dxSpacing
646 ENDDO
647 ENDIF
648 C Y coordinate
649 IF ( delYfile .NE. ' ' ) THEN
650 IF ( delY(1) .NE. UNSET_RL .OR. dySpacing .NE. UNSET_RL ) THEN
651 WRITE(msgBuf,'(A,A)') 'Too many specifications for delY:',
652 & 'Specify only one of delY, dySpacing or delYfile'
653 CALL PRINT_ERROR( msgBuf , 1)
654 STOP 'ABNORMAL END: S/R INI_PARMS'
655 ELSE
656 _BEGIN_MASTER( myThid )
657 IF (readBinaryPrec.EQ.precFloat32) THEN
658 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
659 & ACCESS='DIRECT',RECL=WORDLENGTH*Ny)
660 READ(37,rec=1) delY
661 #ifdef _BYTESWAPIO
662 call MDS_BYTESWAPR4( Ny, delY )
663 #endif
664 CLOSE(37)
665 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
666 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
667 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Ny)
668 READ(37,rec=1) delY
669 #ifdef _BYTESWAPIO
670 call MDS_BYTESWAPR8( Ny, delY )
671 #endif
672 CLOSE(37)
673 ENDIF
674 _END_MASTER(myThid)
675 ENDIF
676 ENDIF
677 IF ( dySpacing .NE. UNSET_RL ) THEN
678 DO i=1,Ny
679 delY(i) = dySpacing
680 ENDDO
681 ENDIF
682 C
683 IF ( rSphere .NE. 0 ) THEN
684 recip_rSphere = 1.D0/rSphere
685 ELSE
686 recip_rSphere = 0.
687 ENDIF
688 C-- Initialize EOS coefficients (3rd order polynomial)
689 IF (eostype.eq.'POLY3') THEN
690 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
691 READ(37,*) I
692 IF (I.NE.Nr) THEN
693 WRITE(msgBuf,'(A)')
694 & 'ini_parms: attempt to read POLY3.COEFFS failed'
695 CALL PRINT_ERROR( msgBuf , 1)
696 WRITE(msgBuf,'(A)')
697 & ' because bad # of levels in data'
698 CALL PRINT_ERROR( msgBuf , 1)
699 STOP 'Bad data in POLY3.COEFFS'
700 ENDIF
701 READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nr)
702 DO K=1,Nr
703 READ(37,*) (eosC(I,K),I=1,9)
704 ENDDO
705 CLOSE(37)
706 ENDIF
707 C-- Check for conflicting grid definitions.
708 goptCount = 0
709 IF ( usingCartesianGrid ) goptCount = goptCount+1
710 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
711 IF ( usingCurvilinearGrid ) goptCount = goptCount+1
712 IF ( goptCount .GT. 1 ) THEN
713 WRITE(msgBuf,'(A)')
714 & 'S/R INI_PARMS: More than one coordinate system requested'
715 CALL PRINT_ERROR( msgBuf , myThid)
716 STOP 'ABNORMAL END: S/R INI_PARMS'
717 ENDIF
718 IF ( goptCount .LT. 1 ) THEN
719 WRITE(msgBuf,'(A)')
720 & 'S/R INI_PARMS: No coordinate system requested'
721 CALL PRINT_ERROR( msgBuf , myThid)
722 STOP 'ABNORMAL END: S/R INI_PARMS'
723 ENDIF
724 C-- Make metric term settings consistent with underlying grid.
725 IF ( usingCartesianGrid ) THEN
726 usingSphericalPolarMterms = .FALSE.
727 metricTerms = .FALSE.
728 mTFacMom = 0.
729 useBetaPlaneF = .TRUE.
730 ENDIF
731 IF ( usingSphericalPolarGrid ) THEN
732 useConstantF = .FALSE.
733 useBetaPlaneF = .FALSE.
734 useSphereF = .TRUE.
735 usingSphericalPolarMterms = metricTerms
736 ENDIF
737 IF ( usingCurvilinearGrid ) THEN
738 useSphereF = .TRUE.
739 ENDIF
740 C-- p, z, r coord parameters
741 DO K = 1, Nr
742 IF ( delZ(K) .NE. UNSET_RL ) zCoordInputData = .TRUE.
743 IF ( delP(K) .NE. UNSET_RL ) pCoordInputData = .TRUE.
744 IF ( delR(K) .NE. UNSET_RL ) rCoordInputData = .TRUE.
745 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delZ(K)
746 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delP(K)
747 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delRDefault(K)
748 IF ( delR(K) .EQ. 0. ) THEN
749 WRITE(msgBuf,'(A,I4)')
750 & 'S/R INI_PARMS: No value for delZ/delP/delR at K = ',K
751 CALL PRINT_ERROR( msgBuf , 1)
752 STOP 'ABNORMAL END: S/R INI_PARMS'
753 ENDIF
754 ENDDO
755 C Check for multiple coordinate systems
756 CoordsSet = 0
757 IF ( zCoordInputData ) coordsSet = coordsSet + 1
758 IF ( pCoordInputData ) coordsSet = coordsSet + 1
759 IF ( rCoordInputData ) coordsSet = coordsSet + 1
760 IF ( coordsSet .GT. 1 ) THEN
761 WRITE(msgBuf,'(A)')
762 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
763 CALL PRINT_ERROR( msgBuf , myThid)
764 STOP 'ABNORMAL END: S/R INI_PARMS'
765 ENDIF
766
767 C-- Input files
768 READ(UNIT=iUnit,NML=PARM05) !,IOSTAT=errIO)
769 IF ( errIO .LT. 0 ) THEN
770 WRITE(msgBuf,'(A)')
771 & 'Error reading numerical model '
772 CALL PRINT_ERROR( msgBuf , 1)
773 WRITE(msgBuf,'(A)')
774 & 'parameter file "data"'
775 CALL PRINT_ERROR( msgBuf , 1)
776 WRITE(msgBuf,'(A)')
777 & 'Problem in namelist PARM05'
778 CALL PRINT_ERROR( msgBuf , 1)
779 CALL MODELDATA_EXAMPLE( myThid )
780 STOP 'ABNORMAL END: S/R INI_PARMS'
781 ENDIF
782
783 C
784 C-- Set factors required for mixing pressure and meters as vertical coordinate.
785 C rkFac is a "sign" parameter which is used where the orientation of the vertical
786 C coordinate (pressure or meters) relative to the vertical index (K) is important.
787 C rkFac = 1 applies when K and the coordinate are in the opposite sense.
788 C rkFac = -1 applies when K and the coordinate are in the same sense.
789 C horiVertRatio is a parameter that maps horizontal units to vertical units.
790 C It is used in certain special cases where lateral and vertical terms are
791 C being combined and a single frame of reference is needed.
792 IF ( zCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
793 rkFac = 1.D0
794 horiVertRatio = 1.D0
795 ENDIF
796 IF ( pCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
797 rkFac = -1.D0
798 horiVertRatio = Gravity * rhoConst
799 ENDIF
800 IF ( rCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
801 rkFac = 1.D0
802 horiVertRatio = 1.D0
803 ENDIF
804 IF (buoyancyRelation.EQ.'ATMOSPHERIC')
805 & horiVertRatio = Gravity * rhoConst
806 IF ( rkFac .EQ. UNSET_RS ) rkFac=rkFacDefault
807 recip_rkFac = 1.D0 / rkFac
808 recip_horiVertRatio = 1./horiVertRatio
809 IF ( zCoordInputData ) usingZCoords = .TRUE.
810 IF ( pCoordInputData ) usingPCoords = .TRUE.
811
812 C
813 CLOSE(iUnit)
814
815 _END_MASTER(myThid)
816
817 C-- Everyone else must wait for the parameters to be loaded
818 _BARRIER
819 C
820
821 RETURN
822 END
823

  ViewVC Help
Powered by ViewVC 1.1.22