/[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.80 - (show annotations) (download)
Fri Jun 21 18:36:05 2002 UTC (21 years, 10 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint45d_post, checkpoint46, checkpoint46a_pre
Changes since 1.79: +9 -2 lines
Added new parameter: deltaTfreesurf

Previously, the free-surface equation was intergrated forward
synchronously with the momentum equations. It is more consistent
to use the tracer time-step. This increases the number of
iterations required but strengthens the damping.

We *SHOULD* make the default time-step equal to the tracer time-step.
However, we don't for backward compatibility. At some point in the
future we need to change the default behaviour.

It turns out that the reason for the "reduced stability" encountered
in large-scale runs seems to be related to excess variability in
the free surface which in turn happens when the waves aren't damped.
Using a longer time-step fixes this.

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

  ViewVC Help
Powered by ViewVC 1.1.22