/[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.57 - (show annotations) (download)
Sun Mar 25 22:33:52 2001 UTC (23 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: c37_adj
Changes since 1.56: +4 -3 lines
Modifications and additions to enable automatic differentiation.
Detailed info's in doc/notes_c37_adj.txt

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

  ViewVC Help
Powered by ViewVC 1.1.22