/[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.50 - (show annotations) (download)
Tue Apr 11 13:51:10 2000 UTC (24 years, 1 month ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint28
Changes since 1.49: +3 -3 lines
Corrected specification of hFacMinDr. If hFacMinDp or hFacMinDz were
used, it used to do something strange...

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_parms.F,v 1.49 2000/04/05 17:52:37 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 #ifdef ALLOW_KPP
26 #include "KPPMIX.h"
27 #endif
28
29 C === Routine arguments ===
30 C myThid - Number of this instance of INI_PARMS
31 INTEGER myThid
32
33 C === Local variables ===
34 C dxSpacing, dySpacing - Default spacing in X and Y.
35 C Units are that of coordinate system
36 C i.e. cartesian => metres
37 C s. polar => degrees
38 C goptCount - Used to count the nuber of grid options
39 C (only one is allowed! )
40 C msgBuf - Informational/error meesage buffer
41 C errIO - IO error flag
42 C iUnit - Work variable for IO unit number
43 C record - Work variable for IO buffer
44 C K, I, J - Loop counters
45 C xxxDefault - Default value for variable xxx
46 _RL dxSpacing
47 _RL dySpacing
48 CHARACTER*(MAX_LEN_FNAM) delXfile
49 CHARACTER*(MAX_LEN_FNAM) delYfile
50 CHARACTER*(MAX_LEN_MBUF) msgBuf
51 CHARACTER*(MAX_LEN_PREC) record
52 INTEGER goptCount
53 INTEGER K, I, J, IL, iUnit
54 INTEGER errIO
55 INTEGER IFNBLNK
56 EXTERNAL IFNBLNK
57 INTEGER ILNBLNK
58 EXTERNAL ILNBLNK
59 C Default values for variables which have vertical coordinate system
60 C dependency.
61 _RL viscArDefault
62 _RL diffKrTDefault
63 _RL diffKrSDefault
64 _RL hFacMinDrDefault
65 _RL delRDefault(Nr)
66 _RS rkFacDefault
67 C zCoordInputData - These are used to select between different coordinate systems.
68 C pCoordInputData The vertical coordinate system in the rest of the model is
69 C rCoordInputData written in terms of r. In the model "data" file input data can
70 C coordsSet be interms of z, p or r.
71 C e.g. delZ or delP or delR for the vertical grid spacing.
72 C The following rules apply:
73 C All parameters must use the same vertical coordinate system.
74 C e.g. delZ and viscAz is legal but
75 C delZ and viscAr is an error.
76 C Similarly specifyinh delZ and delP is an error.
77 C zCoord..., pCoord..., rCoord... are used to flag when z, p or r are
78 C used. coordsSet counts how many vertical coordinate systems have been
79 C used to specify variables. coordsSet > 1 is an error.
80 C
81 LOGICAL zCoordInputData
82 LOGICAL pCoordInputData
83 LOGICAL rCoordInputData
84 INTEGER coordsSet
85
86 C-- Continuous equation parameters
87 NAMELIST /PARM01/
88 & gravity, gBaro, rhonil, tAlpha, sBeta, f0, beta,
89 & viscAh, viscAz, viscA4, cosPower,
90 & diffKhT, diffKzT, diffK4T,
91 & diffKhS, diffKzS, diffK4S,
92 & GMmaxslope,GMlength,GMalpha,GMdepth,GMkbackground,GMmaxval,
93 & tRef, sRef, eosType,
94 & no_slip_sides,no_slip_bottom,
95 & momViscosity, momAdvection, momForcing, useCoriolis,
96 & momPressureForcing, metricTerms,
97 & tempDiffusion, tempAdvection, tempForcing,
98 & saltDiffusion, saltAdvection, saltForcing,
99 & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz,
100 & tempStepping, saltStepping, momStepping,
101 & implicitDiffusion, implicitViscosity,
102 & viscAr, diffKrT, diffKrS, hFacMinDr,
103 & viscAp, diffKpT, diffKpS, hFacMinDp,
104 & rhoConst, buoyancyRelation, HeatCapacity_Cp,
105 & writeBinaryPrec, readBinaryPrec, writeStatePrec,
106 & openBoundaries, nonHydrostatic, globalFiles,
107 & allowFreezing, ivdc_kappa
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 & l, m, n, delP, delR, rkFac,
127 & delXfile, delYfile
128
129 C-- Input files
130 NAMELIST /PARM05/
131 & bathyFile, hydrogThetaFile, hydrogSaltFile,
132 & zonalWindFile, meridWindFile,
133 & thetaClimFile, saltClimFile,
134 & surfQfile, EmPmRfile,
135 & uVelInitFile, vVelInitFile, pSurfInitFile
136
137 C-- Open Boundaries
138 NAMELIST /PARM06/
139 & OB_Jnorth, OB_Jsouth, OB_Ieast, OB_Iwest
140
141 #ifdef ALLOW_KPP
142 C-- KPP vertical mixing parameters
143 NAMELIST /PARM07/
144 & usingKPPmixing, KPPmixingMaps, KPPwriteState,
145 & minKPPviscAz, maxKPPviscAz, minKPPghat, maxKPPghat,
146 & minKPPdiffKzT, maxKPPdiffKzT, minKPPdiffKzS, maxKPPdiffKzS,
147 & epsln, epsilon, vonk, conc1, conam, concm, conc2, zetam,
148 & conas, concs, conc3, zetas,
149 & Ricr, cekman, cmonob, concv, hbf, Vtc,
150 & zmin, zmax, umin, umax,
151 & num_v_smooth_Ri, num_v_smooth_BV,
152 & num_z_smooth_sh, num_m_smooth_sh,
153 & Riinfty, BVSQcon, difm0, difs0, dift0,
154 & difmiw, difsiw, diftiw, difmcon, difscon, diftcon,
155 & cstar, cg
156 #endif
157
158 C
159 _BEGIN_MASTER(myThid)
160
161 C Defaults values for input parameters
162 CALL SET_DEFAULTS(
163 O viscArDefault, diffKrTDefault, diffKrSDefault,
164 O hFacMinDrDefault, delRdefault, rkFacDefault,
165 I myThid )
166
167 C-- Initialise "which vertical coordinate system used" flags.
168 zCoordInputData = .FALSE.
169 pCoordInputData = .FALSE.
170 rCoordInputData = .FALSE.
171 usingPCoords = .FALSE.
172 usingZCoords = .FALSE.
173 coordsSet = 0
174
175 C-- Open the parameter file
176 OPEN(UNIT=scrUnit1,STATUS='SCRATCH')
177 OPEN(UNIT=scrUnit2,STATUS='SCRATCH')
178 OPEN(UNIT=modelDataUnit,FILE='data',STATUS='OLD',
179 & IOSTAT=errIO)
180 IF ( errIO .LT. 0 ) THEN
181 WRITE(msgBuf,'(A)')
182 & 'S/R INI_PARMS'
183 CALL PRINT_ERROR( msgBuf , 1)
184 WRITE(msgBuf,'(A)')
185 & 'Unable to open model parameter'
186 CALL PRINT_ERROR( msgBuf , 1)
187 WRITE(msgBuf,'(A)')
188 & 'file "data"'
189 CALL PRINT_ERROR( msgBuf , 1)
190 CALL MODELDATA_EXAMPLE( myThid )
191 STOP 'ABNORMAL END: S/R INI_PARMS'
192 ENDIF
193
194 DO WHILE ( .TRUE. )
195 READ(modelDataUnit,FMT='(A)',END=1001) RECORD
196 IL = MAX(ILNBLNK(RECORD),1)
197 IF ( RECORD(1:1) .NE. commentCharacter )
198 & WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL)
199 WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL)
200 ENDDO
201 1001 CONTINUE
202 CLOSE(modelDataUnit)
203
204 C-- Report contents of model parameter file
205 WRITE(msgBuf,'(A)')
206 &'// ======================================================='
207 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
208 & SQUEEZE_RIGHT , 1)
209 WRITE(msgBuf,'(A)') '// Model parameter file "data"'
210 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
211 & SQUEEZE_RIGHT , 1)
212 WRITE(msgBuf,'(A)')
213 &'// ======================================================='
214 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
215 & SQUEEZE_RIGHT , 1)
216 iUnit = scrUnit2
217 REWIND(iUnit)
218 DO WHILE ( .TRUE. )
219 READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD
220 IL = MAX(ILNBLNK(RECORD),1)
221 WRITE(msgBuf,'(A,A)') '>',RECORD(:IL)
222 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
223 & SQUEEZE_RIGHT , 1)
224 ENDDO
225 2001 CONTINUE
226 CLOSE(iUnit)
227 WRITE(msgBuf,'(A)') ' '
228 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
229 & SQUEEZE_RIGHT , 1)
230
231
232 C-- Read settings from model parameter file "data".
233 iUnit = scrUnit1
234 REWIND(iUnit)
235
236 C-- Set default "physical" parameters
237 viscAz = UNSET_RL
238 viscAr = UNSET_RL
239 viscAp = UNSET_RL
240 diffKzT = UNSET_RL
241 diffKpT = UNSET_RL
242 diffKrT = UNSET_RL
243 diffKzS = UNSET_RL
244 diffKpS = UNSET_RL
245 diffKrS = UNSET_RL
246 gBaro = UNSET_RL
247 rhoConst = UNSET_RL
248 hFacMinDr = UNSET_RL
249 hFacMinDz = UNSET_RL
250 hFacMinDp = UNSET_RL
251 READ(UNIT=iUnit,NML=PARM01) !,IOSTAT=errIO)
252 IF ( errIO .LT. 0 ) THEN
253 WRITE(msgBuf,'(A)')
254 & 'S/R INI_PARMS'
255 CALL PRINT_ERROR( msgBuf , 1)
256 WRITE(msgBuf,'(A)')
257 & 'Error reading numerical model '
258 CALL PRINT_ERROR( msgBuf , 1)
259 WRITE(msgBuf,'(A)')
260 & 'parameter file "data"'
261 CALL PRINT_ERROR( msgBuf , 1)
262 WRITE(msgBuf,'(A)')
263 & 'Problem in namelist PARM01'
264 CALL PRINT_ERROR( msgBuf , 1)
265 CALL MODELDATA_EXAMPLE( myThid )
266 STOP 'ABNORMAL END: S/R INI_PARMS'
267 ENDIF
268 IF ( implicitFreeSurface ) freeSurfFac = 1.D0
269 IF ( rigidLid ) freeSurfFac = 0.D0
270 IF ( gBaro .EQ. UNSET_RL ) gBaro=gravity
271 IF ( rhoConst .EQ. UNSET_RL ) rhoConst=rhoNil
272 C-- Momentum viscosity on/off flag.
273 IF ( momViscosity ) THEN
274 vfFacMom = 1.D0
275 ELSE
276 vfFacMom = 0.D0
277 ENDIF
278 C-- Momentum advection on/off flag.
279 IF ( momAdvection ) THEN
280 afFacMom = 1.D0
281 ELSE
282 afFacMom = 0.D0
283 ENDIF
284 C-- Momentum forcing on/off flag.
285 IF ( momForcing ) THEN
286 foFacMom = 1.D0
287 ELSE
288 foFacMom = 0.D0
289 ENDIF
290 C-- Coriolis term on/off flag.
291 IF ( useCoriolis ) THEN
292 cfFacMom = 1.D0
293 ELSE
294 cfFacMom = 0.D0
295 ENDIF
296 C-- Pressure term on/off flag.
297 IF ( momPressureForcing ) THEN
298 pfFacMom = 1.D0
299 ELSE
300 pfFacMom = 0.D0
301 ENDIF
302 C-- Metric terms on/off flag.
303 IF ( metricTerms ) THEN
304 mTFacMom = 1.D0
305 ELSE
306 mTFacMom = 1.D0
307 ENDIF
308 C-- z,p,r coord input switching.
309 IF ( viscAz .NE. UNSET_RL ) zCoordInputData = .TRUE.
310 IF ( viscAp .NE. UNSET_RL ) pCoordInputData = .TRUE.
311 IF ( viscAr .NE. UNSET_RL ) rCoordInputData = .TRUE.
312 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAz
313 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAp
314 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscArDefault
315
316 IF ( diffKzT .NE. UNSET_RL ) zCoordInputData = .TRUE.
317 IF ( diffKpT .NE. UNSET_RL ) pCoordInputData = .TRUE.
318 IF ( diffKrT .NE. UNSET_RL ) rCoordInputData = .TRUE.
319 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKzT
320 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKpT
321 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKrTDefault
322
323 IF ( diffKzS .NE. UNSET_RL ) zCoordInputData = .TRUE.
324 IF ( diffKpS .NE. UNSET_RL ) pCoordInputData = .TRUE.
325 IF ( diffKrS .NE. UNSET_RL ) rCoordInputData = .TRUE.
326 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKzS
327 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKpS
328 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKrSDefault
329
330 IF ( hFacMinDz .NE. UNSET_RL ) zCoordInputData = .TRUE.
331 IF ( hFacMinDp .NE. UNSET_RL ) pCoordInputData = .TRUE.
332 IF ( hFacMinDr .NE. UNSET_RL ) rCoordInputData = .TRUE.
333 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDz
334 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDp
335 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDrDefault
336
337 IF ( ivdc_kappa .NE. 0. .AND. .NOT. implicitDiffusion ) THEN
338 WRITE(msgBuf,'(A,A)')
339 & 'S/R INI_PARMS: To use ivdc_kappa you must enable implicit',
340 & ' vertical diffusion.'
341 CALL PRINT_ERROR( msgBuf , myThid)
342 STOP 'ABNORMAL END: S/R INI_PARMS'
343 ENDIF
344
345 IF ( implicitFreeSurface .AND. rigidLid ) THEN
346 WRITE(msgBuf,'(A,A)')
347 & 'S/R INI_PARMS: Cannot select both implicitFreeSurface',
348 & ' and rigidLid.'
349 CALL PRINT_ERROR( msgBuf , myThid)
350 STOP 'ABNORMAL END: S/R INI_PARMS'
351 ENDIF
352 coordsSet = 0
353 IF ( zCoordInputData ) coordsSet = coordsSet + 1
354 IF ( pCoordInputData ) coordsSet = coordsSet + 1
355 IF ( rCoordInputData ) coordsSet = coordsSet + 1
356 IF ( coordsSet .GT. 1 ) THEN
357 WRITE(msgBuf,'(A)')
358 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
359 CALL PRINT_ERROR( msgBuf , myThid)
360 STOP 'ABNORMAL END: S/R INI_PARMS'
361 ENDIF
362 IF ( rhoConst .LE. 0. ) THEN
363 WRITE(msgBuf,'(A)')
364 & 'S/R INI_PARMS: rhoConst must be greater than 0.'
365 CALL PRINT_ERROR( msgBuf , myThid)
366 STOP 'ABNORMAL END: S/R INI_PARMS'
367 ELSE
368 recip_rhoConst = 1.D0 / rhoConst
369 ENDIF
370 IF ( rhoNil .LE. 0. ) THEN
371 WRITE(msgBuf,'(A)')
372 & 'S/R INI_PARMS: rhoNil must be greater than 0.'
373 CALL PRINT_ERROR( msgBuf , myThid)
374 STOP 'ABNORMAL END: S/R INI_PARMS'
375 ELSE
376 recip_rhoNil = 1.D0 / rhoNil
377 ENDIF
378 IF ( HeatCapacity_Cp .LE. 0. ) THEN
379 WRITE(msgBuf,'(A)')
380 & 'S/R INI_PARMS: HeatCapacity_Cp must be greater than 0.'
381 CALL PRINT_ERROR( msgBuf , myThid)
382 STOP 'ABNORMAL END: S/R INI_PARMS'
383 ELSE
384 recip_Cp = 1.D0 / HeatCapacity_Cp
385 ENDIF
386 IF ( gravity .LE. 0. ) THEN
387 WRITE(msgBuf,'(A)')
388 & 'S/R INI_PARMS: gravity must be greater than 0.'
389 CALL PRINT_ERROR( msgBuf , myThid)
390 STOP 'ABNORMAL END: S/R INI_PARMS'
391 ELSE
392 recip_gravity = 1.D0 / gravity
393 ENDIF
394 C Set globalFiles flag for READ_WRITE_FLD package
395 CALL SET_WRITE_GLOBAL_FLD( globalFiles )
396 C Set globalFiles flag for READ_WRITE_REC package
397 CALL SET_WRITE_GLOBAL_REC( globalFiles )
398 C Set globalFiles flag for READ_WRITE_REC package
399 CALL SET_WRITE_GLOBAL_PICKUP( globalFiles )
400
401 C-- Elliptic solver parameters
402 READ(UNIT=iUnit,NML=PARM02) !,IOSTAT=errIO)
403 IF ( errIO .LT. 0 ) THEN
404 WRITE(msgBuf,'(A)')
405 & 'S/R INI_PARMS'
406 CALL PRINT_ERROR( msgBuf , 1)
407 WRITE(msgBuf,'(A)')
408 & 'Error reading numerical model '
409 CALL PRINT_ERROR( msgBuf , 1)
410 WRITE(msgBuf,'(A)')
411 & 'parameter file "data".'
412 CALL PRINT_ERROR( msgBuf , 1)
413 WRITE(msgBuf,'(A)')
414 & 'Problem in namelist PARM02'
415 CALL PRINT_ERROR( msgBuf , 1)
416 CALL MODELDATA_EXAMPLE( myThid )
417 STOP 'ABNORMAL END: S/R INI_PARMS'
418 ENDIF
419
420 C-- Time stepping parameters
421 rCD = -1.D0
422 READ(UNIT=iUnit,NML=PARM03) !,IOSTAT=errIO)
423 IF ( errIO .LT. 0 ) THEN
424 WRITE(msgBuf,'(A)')
425 & 'S/R INI_PARMS'
426 CALL PRINT_ERROR( msgBuf , 1)
427 WRITE(msgBuf,'(A)')
428 & 'Error reading numerical model '
429 CALL PRINT_ERROR( msgBuf , 1)
430 WRITE(msgBuf,'(A)')
431 & 'parameter file "data"'
432 CALL PRINT_ERROR( msgBuf , 1)
433 WRITE(msgBuf,'(A)')
434 & 'Problem in namelist PARM03'
435 CALL PRINT_ERROR( msgBuf , 1)
436 CALL MODELDATA_EXAMPLE( myThid )
437 STOP 'ABNORMAL END: S/R INI_PARMS'
438 ENDIF
439 C Process "timestepping" params
440 C o Time step size
441 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
442 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
443 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
444 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
445 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
446 IF ( periodicExternalForcing ) THEN
447 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
448 WRITE(msgBuf,'(A)')
449 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
450 CALL PRINT_ERROR( msgBuf , 1)
451 STOP 'ABNORMAL END: S/R INI_PARMS'
452 ENDIF
453 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
454 & externForcingCycle/externForcingPeriod ) THEN
455 WRITE(msgBuf,'(A)')
456 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
457 CALL PRINT_ERROR( msgBuf , 1)
458 STOP 'ABNORMAL END: S/R INI_PARMS'
459 ENDIF
460 IF ( externForcingCycle.le.externForcingPeriod ) THEN
461 WRITE(msgBuf,'(A)')
462 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
463 CALL PRINT_ERROR( msgBuf , 1)
464 STOP 'ABNORMAL END: S/R INI_PARMS'
465 ENDIF
466 IF ( externForcingPeriod.lt.deltaTclock ) THEN
467 WRITE(msgBuf,'(A)')
468 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
469 CALL PRINT_ERROR( msgBuf , 1)
470 STOP 'ABNORMAL END: S/R INI_PARMS'
471 ENDIF
472 ENDIF
473 C o Convection frequency
474 IF ( cAdjFreq .LT. 0. ) THEN
475 cAdjFreq = deltaTClock
476 ENDIF
477 IF ( ivdc_kappa .NE. 0. .AND. cAdjFreq .NE. 0. ) THEN
478 WRITE(msgBuf,'(A,A)')
479 & 'S/R INI_PARMS: You have enabled both ivdc_kappa and',
480 & ' convective adjustment.'
481 CALL PRINT_ERROR( msgBuf , myThid)
482 STOP 'ABNORMAL END: S/R INI_PARMS'
483 ENDIF
484 C o CD coupling
485 IF ( tauCD .EQ. 0.D0 ) THEN
486 tauCD = deltaTmom
487 ENDIF
488 IF ( rCD .LT. 0. ) THEN
489 rCD = 1. - deltaTMom/tauCD
490 ENDIF
491 C o Temperature climatology relaxation time scale
492 IF ( tauThetaClimRelax .EQ. 0.D0 ) THEN
493 doThetaClimRelax = .FALSE.
494 lambdaThetaClimRelax = 0.D0
495 ELSE
496 doThetaClimRelax = .TRUE.
497 lambdaThetaClimRelax = 1./tauThetaClimRelax
498 ENDIF
499 C o Salinity climatology relaxation time scale
500 IF ( tauSaltClimRelax .EQ. 0.D0 ) THEN
501 doSaltClimRelax = .FALSE.
502 lambdaSaltClimRelax = 0.D0
503 ELSE
504 doSaltClimRelax = .TRUE.
505 lambdaSaltClimRelax = 1./tauSaltClimRelax
506 ENDIF
507
508 C o Start time
509 IF ( nIter0 .NE. 0 .AND. startTime .EQ. 0. )
510 & startTime = deltaTClock*float(nIter0)
511 C o nIter0
512 IF ( nIter0 .EQ. 0 .AND. startTime .NE. 0. )
513 & nIter0 = INT( startTime/deltaTClock )
514
515 C o nTimeSteps 1
516 IF ( nTimeSteps .EQ. 0 .AND. nEndIter .NE. 0 )
517 & nTimeSteps = nEndIter-nIter0
518 C o nTimeSteps 2
519 IF ( nTimeSteps .EQ. 0 .AND. endTime .NE. 0. )
520 & nTimeSteps = int(0.5+(endTime-startTime)/deltaTclock)
521 C o nEndIter 1
522 IF ( nEndIter .EQ. 0 .AND. nTimeSteps .NE. 0 )
523 & nEndIter = nIter0+nTimeSteps
524 C o nEndIter 2
525 IF ( nEndIter .EQ. 0 .AND. endTime .NE. 0. )
526 & nEndIter = int(0.5+endTime/deltaTclock)
527 C o End Time 1
528 IF ( endTime .EQ. 0. .AND. nTimeSteps .NE. 0 )
529 & endTime = startTime + deltaTClock*float(nTimeSteps)
530 C o End Time 2
531 IF ( endTime .EQ. 0. .AND. nEndIter .NE. 0 )
532 & endTime = deltaTClock*float(nEndIter)
533
534 C o Consistent?
535 IF ( nEndIter .NE. nIter0+nTimeSteps ) THEN
536 WRITE(msgBuf,'(A)')
537 & 'S/R INI_PARMS: nIter0, nTimeSteps and nEndIter are inconsistent'
538 CALL PRINT_ERROR( msgBuf , 1)
539 WRITE(msgBuf,'(A)')
540 & 'S/R INI_PARMS: Perhaps more than two were set at once'
541 CALL PRINT_ERROR( msgBuf , 1)
542 STOP 'ABNORMAL END: S/R INI_PARMS'
543 ENDIF
544 IF ( nTimeSteps .NE. int(0.5+(endTime-startTime)/deltaTClock) )
545 & THEN
546 WRITE(msgBuf,'(A)')
547 & 'S/R INI_PARMS: both endTime and nTimeSteps have been set'
548 CALL PRINT_ERROR( msgBuf , 1)
549 WRITE(msgBuf,'(A)')
550 & 'S/R INI_PARMS: but are inconsistent'
551 CALL PRINT_ERROR( msgBuf , 1)
552 STOP 'ABNORMAL END: S/R INI_PARMS'
553 ENDIF
554
555 C o If taveFreq is finite, then we must make sure the diagnostics
556 C code is being compiled
557 #ifndef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
558 IF (taveFreq.NE.0.) THEN
559 WRITE(msgBuf,'(A)')
560 & 'S/R INI_PARMS: taveFreq <> 0 but you have'
561 CALL PRINT_ERROR( msgBuf , 1)
562 WRITE(msgBuf,'(A)')
563 & 'not compiled the model with the diagnostics routines.'
564 CALL PRINT_ERROR( msgBuf , 1)
565 WRITE(msgBuf,'(A,A)')
566 & 'Re-compile with: #define INCLUDE_DIAGNOSTICS_INTERFACE_CODE',
567 & ' or -DINCLUDE_DIAGNOSTICS_INTERFACE_CODE'
568 CALL PRINT_ERROR( msgBuf , 1)
569 STOP 'ABNORMAL END: S/R INI_PARMS'
570 ENDIF
571 #endif
572
573 C-- Grid parameters
574 C In cartesian coords distances are in metres
575 rkFac = UNSET_RS
576 DO K =1,Nr
577 delZ(K) = UNSET_RL
578 delP(K) = UNSET_RL
579 delR(K) = UNSET_RL
580 ENDDO
581 C In spherical polar distances are in degrees
582 recip_rSphere = 1.D0/rSphere
583 dxSpacing = UNSET_RL
584 dySpacing = UNSET_RL
585 delXfile = ' '
586 delYfile = ' '
587 READ(UNIT=iUnit,NML=PARM04) !,IOSTAT=errIO)
588 IF ( errIO .LT. 0 ) THEN
589 WRITE(msgBuf,'(A)')
590 & 'S/R INI_PARMS'
591 CALL PRINT_ERROR( msgBuf , 1)
592 WRITE(msgBuf,'(A)')
593 & 'Error reading numerical model '
594 CALL PRINT_ERROR( msgBuf , 1)
595 WRITE(msgBuf,'(A)')
596 & 'parameter file "data"'
597 CALL PRINT_ERROR( msgBuf , 1)
598 WRITE(msgBuf,'(A)')
599 & 'Problem in namelist PARM04'
600 CALL PRINT_ERROR( msgBuf , 1)
601 CALL MODELDATA_EXAMPLE( myThid )
602 STOP 'ABNORMAL END: S/R INI_PARMS'
603 ENDIF
604
605 C X coordinate
606 IF ( delXfile .NE. ' ' ) THEN
607 IF ( delX(1) .NE. UNSET_RL .OR. dxSpacing .NE. UNSET_RL ) THEN
608 WRITE(msgBuf,'(A,A)') 'Too many specifications for delX:',
609 & 'Specify only one of delX, dxSpacing or delXfile'
610 CALL PRINT_ERROR( msgBuf , 1)
611 STOP 'ABNORMAL END: S/R INI_PARMS'
612 ELSE
613 _BEGIN_MASTER( myThid )
614 IF (readBinaryPrec.EQ.precFloat32) THEN
615 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
616 & ACCESS='DIRECT',RECL=WORDLENGTH*Nx)
617 READ(37,rec=1) delX
618 #ifdef _BYTESWAPIO
619 call MDS_BYTESWAPR4( Nx, delX )
620 #endif
621 CLOSE(37)
622 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
623 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
624 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Nx)
625 READ(37,rec=1) delX
626 #ifdef _BYTESWAPIO
627 call MDS_BYTESWAPR8( Nx, delX )
628 #endif
629 CLOSE(37)
630 ENDIF
631 _END_MASTER(myThid)
632 ENDIF
633 ENDIF
634 IF ( dxSpacing .NE. UNSET_RL ) THEN
635 DO i=1,Nx
636 delX(i) = dxSpacing
637 ENDDO
638 ENDIF
639 C Y coordinate
640 IF ( delYfile .NE. ' ' ) THEN
641 IF ( delY(1) .NE. UNSET_RL .OR. dySpacing .NE. UNSET_RL ) THEN
642 WRITE(msgBuf,'(A,A)') 'Too many specifications for delY:',
643 & 'Specify only one of delY, dySpacing or delYfile'
644 CALL PRINT_ERROR( msgBuf , 1)
645 STOP 'ABNORMAL END: S/R INI_PARMS'
646 ELSE
647 _BEGIN_MASTER( myThid )
648 IF (readBinaryPrec.EQ.precFloat32) THEN
649 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
650 & ACCESS='DIRECT',RECL=WORDLENGTH*Ny)
651 READ(37,rec=1) delY
652 #ifdef _BYTESWAPIO
653 call MDS_BYTESWAPR4( Ny, delY )
654 #endif
655 CLOSE(37)
656 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
657 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
658 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Ny)
659 READ(37,rec=1) delY
660 #ifdef _BYTESWAPIO
661 call MDS_BYTESWAPR8( Ny, delY )
662 #endif
663 CLOSE(37)
664 ENDIF
665 _END_MASTER(myThid)
666 ENDIF
667 ENDIF
668 IF ( dySpacing .NE. UNSET_RL ) THEN
669 DO i=1,Ny
670 delY(i) = dySpacing
671 ENDDO
672 ENDIF
673 C
674 IF ( rSphere .NE. 0 ) THEN
675 recip_rSphere = 1.D0/rSphere
676 ELSE
677 recip_rSphere = 0.
678 ENDIF
679 C-- Initialize EOS coefficients (3rd order polynomial)
680 IF (eostype.eq.'POLY3') THEN
681 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
682 READ(37,*) I
683 IF (I.NE.Nr) THEN
684 WRITE(msgBuf,'(A)')
685 & 'ini_parms: attempt to read POLY3.COEFFS failed'
686 CALL PRINT_ERROR( msgBuf , 1)
687 WRITE(msgBuf,'(A)')
688 & ' because bad # of levels in data'
689 CALL PRINT_ERROR( msgBuf , 1)
690 STOP 'Bad data in POLY3.COEFFS'
691 ENDIF
692 READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nr)
693 DO K=1,Nr
694 READ(37,*) (eosC(I,K),I=1,9)
695 ENDDO
696 CLOSE(37)
697 ENDIF
698 C-- Check for conflicting grid definitions.
699 goptCount = 0
700 IF ( usingCartesianGrid ) goptCount = goptCount+1
701 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
702 IF ( goptCount .NE. 1 ) THEN
703 WRITE(msgBuf,'(A)')
704 & 'S/R INI_PARMS: More than one coordinate system requested'
705 CALL PRINT_ERROR( msgBuf , myThid)
706 STOP 'ABNORMAL END: S/R INI_PARMS'
707 ENDIF
708 C-- Make metric term settings consistent with underlying grid.
709 IF ( usingCartesianGrid ) THEN
710 usingSphericalPolarMterms = .FALSE.
711 metricTerms = .FALSE.
712 mTFacMom = 0
713 useBetaPlaneF = .TRUE.
714 ENDIF
715 IF ( usingSphericalPolarGrid ) THEN
716 useConstantF = .FALSE.
717 useBetaPlaneF = .FALSE.
718 useSphereF = .TRUE.
719 omega = 2.D0 * PI / ( 3600.D0 * 24.D0 )
720 usingSphericalPolarMterms = .TRUE.
721 metricTerms = .TRUE.
722 mTFacMom = 1
723 ENDIF
724 C-- p, z, r coord parameters
725 DO K = 1, Nr
726 IF ( delZ(K) .NE. UNSET_RL ) zCoordInputData = .TRUE.
727 IF ( delP(K) .NE. UNSET_RL ) pCoordInputData = .TRUE.
728 IF ( delR(K) .NE. UNSET_RL ) rCoordInputData = .TRUE.
729 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delZ(K)
730 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delP(K)
731 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delRDefault(K)
732 IF ( delR(K) .EQ. 0. ) THEN
733 WRITE(msgBuf,'(A,I4)')
734 & 'S/R INI_PARMS: No value for delZ/delP/delR at K = ',K
735 CALL PRINT_ERROR( msgBuf , 1)
736 STOP 'ABNORMAL END: S/R INI_PARMS'
737 ENDIF
738 ENDDO
739 C Check for multiple coordinate systems
740 CoordsSet = 0
741 IF ( zCoordInputData ) coordsSet = coordsSet + 1
742 IF ( pCoordInputData ) coordsSet = coordsSet + 1
743 IF ( rCoordInputData ) coordsSet = coordsSet + 1
744 IF ( coordsSet .GT. 1 ) THEN
745 WRITE(msgBuf,'(A)')
746 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
747 CALL PRINT_ERROR( msgBuf , myThid)
748 STOP 'ABNORMAL END: S/R INI_PARMS'
749 ENDIF
750
751 C-- Input files
752 READ(UNIT=iUnit,NML=PARM05) !,IOSTAT=errIO)
753 IF ( errIO .LT. 0 ) THEN
754 WRITE(msgBuf,'(A)')
755 & 'S/R INI_PARMS'
756 CALL PRINT_ERROR( msgBuf , 1)
757 WRITE(msgBuf,'(A)')
758 & 'Error reading numerical model '
759 CALL PRINT_ERROR( msgBuf , 1)
760 WRITE(msgBuf,'(A)')
761 & 'parameter file "data"'
762 CALL PRINT_ERROR( msgBuf , 1)
763 WRITE(msgBuf,'(A)')
764 & 'Problem in namelist PARM05'
765 CALL PRINT_ERROR( msgBuf , 1)
766 CALL MODELDATA_EXAMPLE( myThid )
767 STOP 'ABNORMAL END: S/R INI_PARMS'
768 ENDIF
769
770 C
771 C-- Set factors required for mixing pressure and meters as vertical coordinate.
772 C rkFac is a "sign" parameter which is used where the orientation of the vertical
773 C coordinate (pressure or meters) relative to the vertical index (K) is important.
774 C rkFac = 1 applies when K and the coordinate are in the opposite sense.
775 C rkFac = -1 applies when K and the coordinate are in the same sense.
776 C horiVertRatio is a parameter that maps horizontal units to vertical units.
777 C It is used in certain special cases where lateral and vertical terms are
778 C being combined and a single frame of reference is needed.
779 IF ( zCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
780 rkFac = 1.D0
781 horiVertRatio = 1.D0
782 ENDIF
783 IF ( pCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
784 rkFac = -1.D0
785 horiVertRatio = Gravity * rhoConst
786 ENDIF
787 IF ( rCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
788 rkFac = 1.D0
789 horiVertRatio = 1.D0
790 ENDIF
791 IF ( rkFac .EQ. UNSET_RS ) rkFac=rkFacDefault
792 recip_rkFac = 1.D0 / rkFac
793 recip_horiVertRatio = 1./horiVertRatio
794 IF ( zCoordInputData ) usingZCoords = .TRUE.
795 IF ( pCoordInputData ) usingPCoords = .TRUE.
796
797 C-- OBCS
798 IF (openBoundaries) THEN
799 READ(UNIT=iUnit,NML=PARM06)
800 DO J=1,Ny
801 if (OB_Ieast(J).lt.0) OB_Ieast(J)=OB_Ieast(J)+Nx+1
802 ENDDO
803 DO I=1,Nx
804 if (OB_Jnorth(I).lt.0) OB_Jnorth(I)=OB_Jnorth(I)+Ny+1
805 ENDDO
806 write(0,*) 'OB Jn =',OB_Jnorth
807 write(0,*) 'OB Js =',OB_Jsouth
808 write(0,*) 'OB Ie =',OB_Ieast
809 write(0,*) 'OB Iw =',OB_Iwest
810 ENDIF
811
812 #ifdef ALLOW_KPP
813 C-- KPP vertical mixing parameters
814 usingKPPmixing = .FALSE.
815 KPPmixingMaps = .FALSE.
816 KPPwriteState = .FALSE.
817 minKPPghat = 0.
818 maxKPPghat = 1.e10
819 minKPPviscAz = 0.
820 DO K =1,Nr
821 maxKPPviscAz(K) = 100.
822 ENDDO
823 minKPPdiffKzT = 0.
824 maxKPPdiffKzT = 100.
825 minKPPdiffKzS = 0.
826 maxKPPdiffKzS = 100.
827
828 c-----------------------------------------------------------------------
829 c define some non-dimensional constants and
830 c the vertical mixing coefficients in m-k-s units
831 c-----------------------------------------------------------------------
832
833 epsln = 1.e-20 !
834 epsilon = 0.1 ! nondimensional extent of surface layer
835 vonk = 0.40 ! von Karmans constant
836 conc1 = 5.0 ! scalar coefficients
837 conam = 1.257 ! .
838 concm = 8.380 ! .
839 conc2 = 16.0 ! .
840 zetam = -0.2 ! .
841 conas = -28.86 ! .
842 concs = 98.96 ! .
843 conc3 = 16.0 ! .
844 zetas = -1.0 ! .
845
846 c parameters for subroutine "bldepth"
847
848 Ricr = 0.30 ! critical bulk Richardson Number
849 cekman = 0.7 ! coefficient for ekman depth
850 cmonob = 1.0 ! coefficient for Monin-Obukhov depth
851 concv = 1.8 ! ratio of interior buoyancy frequency to
852 ! buoyancy frequency at entrainment depth
853 hbf = 1.0 ! fraction of bounadry layer depth to
854 ! which absorbed solar radiation
855 ! contributes to surface buoyancy forcing
856 Vtc = concv * sqrt(0.2/concs/epsilon) / vonk**2 / Ricr
857 ! non-dimensional coefficient for velocity
858 ! scale of turbulant velocity shear
859
860 c parameters and common arrays for subroutines "kmixinit" and "wscale"
861
862 zmin = -4.e-7 ! minimum limit for zehat in table (m3/s3)
863 zmax = 0.0 ! maximum limit for zehat in table (m3/s3)
864 umin = 0.0 ! minimum limit for ustar in table (m/s)
865 umax = .04 ! maximum limit for ustar in table (m/s)
866
867 c parameters for subroutine "Ri_iwmix"
868
869 num_v_smooth_Ri = 1 ! number of times Ri is vertically smoothed
870 num_v_smooth_BV = 1 ! number of times BV is vertically smoothed
871 num_z_smooth_sh = 0 ! number of times shear is zonally smoothed
872 num_m_smooth_sh = 0 ! number of times shear is meridionally smoothed
873 Riinfty = 0.7 ! local Richardson Number limit for shear instability
874 BVSQcon = -0.2e-4 ! Brunt-Vaisala squared (1/s^2)
875 difm0 = 0.005 ! max visc due to shear instability (m^2/s)
876 difs0 = 0.005 ! max tracer diffusivity .. (m^2/s)
877 dift0 = 0.005 ! max heat diffusivity .. (m^2/s)
878 difmiw = 0.001 ! viscosity from background internal wave (m^2/s)
879 difsiw = 0.00003 ! tracer diffusivity .. (m^2/s)
880 diftiw = 0.00003 ! heat diffusivity .. (m^2/s)
881 difmcon = 0.1 ! viscosity due to convective instability (m^2/s)
882 difscon = 0.1 ! tracer diffusivity .. (m^2/s)
883 diftcon = 0.1 ! heat diffusivity .. (m^2/s)
884
885 c parameters for subroutine "blmix"
886
887 cstar = 10. ! proportionality coefficient for nonlocal transport
888
889 ! non-dimensional coefficient for counter-gradient term
890 cg = cstar * vonk * (concs * vonk * epsilon)**(1./3.)
891
892 READ(UNIT=iUnit,NML=PARM07,IOSTAT=errIO,err=13)
893 IF ( errIO .GE. 0 ) GOTO 14
894 13 CONTINUE
895 WRITE(msgBuf,'(A)')
896 & 'S/R INI_PARMS'
897 CALL PRINT_ERROR( msgBuf , 1)
898 WRITE(msgBuf,'(A)')
899 & 'Error reading numerical model '
900 CALL PRINT_ERROR( msgBuf , 1)
901 WRITE(msgBuf,'(A)')
902 & 'parameter file "data"'
903 CALL PRINT_ERROR( msgBuf , 1)
904 WRITE(msgBuf,'(A)')
905 & 'Problem in namelist PARM07'
906 CALL PRINT_ERROR( msgBuf , 1)
907 CALL MODELDATA_EXAMPLE( myThid )
908 STOP 'ABNORMAL END: S/R INI_PARMS'
909 14 CONTINUE
910 #endif
911
912 C
913 CLOSE(iUnit)
914
915 _END_MASTER(myThid)
916
917 C-- Everyone else must wait for the parameters to be loaded
918 _BARRIER
919 C
920
921 RETURN
922 END
923

  ViewVC Help
Powered by ViewVC 1.1.22