/[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.22 - (show annotations) (download)
Wed Jul 1 19:57:22 1998 UTC (25 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint11, checkpoint10, checkpoint12, branch-point-rdot
Branch point for: branch-rdot
Changes since 1.21: +2 -2 lines
Changed the #include at the top from CPP_EEOPTIONS.h to CPP_OPTIONS.h
This ought to be done through-out (?) but was necessary in at least
the_model_main.F, dynamics.F, ini_parms.F because of the new
macro ALLOW_DIAGNOSTICS.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_parms.F,v 1.21 1998/07/01 19:49:36 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
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "EEPARAMS.h"
21 #include "PARAMS.h"
22 #include "CG2D.h"
23
24 C === Routine arguments ===
25 C myThid - Number of this instance of INI_PARMS
26 INTEGER myThid
27
28 C === Local variables ===
29 C dxSpacing, dySpacing - Default spacing in X and Y.
30 C Units are that of coordinate system
31 C i.e. cartesian => metres
32 C s. polar => degrees
33 C goptCount - Used to count the nuber of grid options
34 C (only one is allowed! )
35 C msgBuf - Informational/error meesage buffer
36 C errIO - IO error flag
37 C iUnit - Work variable for IO unit number
38 C record - Work variable for IO buffer
39 C K, I, J - Loop counters
40 REAL dxSpacing
41 REAL dySpacing
42 CHARACTER*(MAX_LEN_MBUF) msgBuf
43 CHARACTER*(MAX_LEN_PREC) record
44 INTEGER goptCount
45 INTEGER K, I, J, IL, iUnit
46 INTEGER errIO
47 INTEGER IFNBLNK
48 EXTERNAL IFNBLNK
49 INTEGER ILNBLNK
50 EXTERNAL ILNBLNK
51
52 C-- Continuous equation parameters
53 NAMELIST /PARM01/
54 & gravity, gBaro, rhonil, tAlpha, sBeta, f0, beta,
55 & viscAh, viscAz, viscA4,
56 & diffKhT, diffKzT, diffK4T,
57 & diffKhS, diffKzS, diffK4S,
58 & GMmaxslope,GMlength,GMalpha,GMdepth,GMkbackground,
59 & tRef, sRef, eosType,
60 & momViscosity, momAdvection, momForcing, useCoriolis,
61 & momPressureForcing, metricTerms,
62 & tempDiffusion, tempAdvection, tempForcing,
63 & saltDiffusion, saltAdvection, saltForcing,
64 & implicitFreeSurface, rigidLid, freeSurfFac,
65 & tempStepping, saltStepping, momStepping, implicitDiffusion
66
67 C-- Elliptic solver parameters
68 NAMELIST /PARM02/
69 & cg2dMaxIters, cg2dChkResFreq, cg2dTargetResidual, cg2dpcOffDFac
70
71 C-- Time stepping parammeters
72 NAMELIST /PARM03/
73 & nIter0, nTimeSteps, deltaT, deltaTmom, deltaTtracer, abEps, tauCD, rCD,
74 & startTime, endTime, chkPtFreq, dumpFreq, taveFreq, deltaTClock,
75 & pChkPtFreq, cAdjFreq, tauThetaClimRelax, tauSaltClimRelax,
76 & periodicExternalForcing, externForcingPeriod, externForcingCycle
77
78 C-- Gridding parameters
79 NAMELIST /PARM04/
80 & usingCartesianGrid, delZ, dxSpacing, dySpacing, delX, delY,
81 & usingSphericalPolarGrid, phiMin, thetaMin, rSphere,
82 & l, m, n
83
84 C-- Input files
85 NAMELIST /PARM05/
86 & bathyFile, hydrogThetaFile, hydrogSaltFile,
87 & zonalWindFile, meridWindFile, thetaClimFile,
88 & saltClimFile
89
90 C
91 _BEGIN_MASTER(myThid)
92
93 C-- Open the parameter file
94 OPEN(UNIT=scrUnit1,STATUS='SCRATCH')
95 OPEN(UNIT=scrUnit2,STATUS='SCRATCH')
96 OPEN(UNIT=modelDataUnit,FILE='data',STATUS='OLD',err=1,IOSTAT=errIO)
97 IF ( errIO .GE. 0 ) GOTO 2
98 1 CONTINUE
99 WRITE(msgBuf,'(A)')
100 & 'S/R INI_PARMS'
101 CALL PRINT_ERROR( msgBuf , 1)
102 WRITE(msgBuf,'(A)')
103 & 'Unable to open model parameter'
104 CALL PRINT_ERROR( msgBuf , 1)
105 WRITE(msgBuf,'(A)')
106 & 'file "data"'
107 CALL PRINT_ERROR( msgBuf , 1)
108 CALL MODELDATA_EXAMPLE( myThid )
109 STOP 'ABNORMAL END: S/R INI_PARMS'
110 2 CONTINUE
111
112 1000 CONTINUE
113 READ(modelDataUnit,FMT='(A)',END=1001) RECORD
114 IL = MAX(ILNBLNK(RECORD),1)
115 IF ( RECORD(1:1) .NE. commentCharacter )
116 & WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL)
117 WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL)
118 GOTO 1000
119 1001 CONTINUE
120 CLOSE(modelDataUnit)
121
122 C-- Report contents of model parameter file
123 WRITE(msgBuf,'(A)')
124 &'// ======================================================='
125 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, SQUEEZE_RIGHT , 1)
126 WRITE(msgBuf,'(A)') '// Model parameter file "data"'
127 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, SQUEEZE_RIGHT , 1)
128 WRITE(msgBuf,'(A)')
129 &'// ======================================================='
130 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
131 & SQUEEZE_RIGHT , 1)
132 iUnit = scrUnit2
133 REWIND(iUnit)
134 2000 CONTINUE
135 READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD
136 IL = MAX(ILNBLNK(RECORD),1)
137 WRITE(msgBuf,'(A,A)') '>',RECORD(:IL)
138 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, SQUEEZE_RIGHT , 1)
139 GOTO 2000
140 2001 CONTINUE
141 CLOSE(iUnit)
142 WRITE(msgBuf,'(A)') ' '
143 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
144 & SQUEEZE_RIGHT , 1)
145
146
147 C-- Read settings from model parameter file "data".
148 iUnit = scrUnit1
149 REWIND(iUnit)
150
151 C-- Set default "physical" parameters
152 DO K =1,Nz
153 tRef(K) = 30.D0 - FLOAT(K)
154 ENDDO
155 gravity = 9.81 D0
156 gBaro = gravity
157 rhoNil = 999.8 D0
158 f0=1.D-4
159 beta = 1. _d -11
160 viscAh=1.d3
161 diffKhT=1.0d3
162 diffKhS=1.0d3
163 viscAz=1.d-3
164 diffKzT=1.d-5
165 diffKzS=1.d-5
166 viscA4=0.
167 diffK4T=0.
168 diffK4S=0.
169 GMmaxslope=1.d-2
170 GMlength=200.d3
171 GMalpha=0.
172 GMdepth=1000.
173 GMkbackground=0.
174 tAlpha=2.d-4
175 sBeta=7.4d-4
176 eosType='LINEAR'
177 implicitFreeSurface = .TRUE.
178 rigidLid = .FALSE.
179 freeSurfFac = 1. _d 0
180 momViscosity = .TRUE.
181 momAdvection = .TRUE.
182 momForcing = .TRUE.
183 useCoriolis = .TRUE.
184 momPressureForcing = .TRUE.
185 momStepping = .TRUE.
186 tempStepping = .TRUE.
187 saltStepping = .TRUE.
188 metricTerms = .TRUE.
189 implicitDiffusion = .FALSE.
190 READ(UNIT=iUnit,NML=PARM01,IOSTAT=errIO,err=3)
191 IF ( errIO .GE. 0 ) GOTO 4
192 3 CONTINUE
193 WRITE(msgBuf,'(A)')
194 & 'S/R INI_PARMS'
195 CALL PRINT_ERROR( msgBuf , 1)
196 WRITE(msgBuf,'(A)')
197 & 'Error reading numerical model '
198 CALL PRINT_ERROR( msgBuf , 1)
199 WRITE(msgBuf,'(A)')
200 & 'parameter file "data"'
201 CALL PRINT_ERROR( msgBuf , 1)
202 WRITE(msgBuf,'(A)')
203 & 'Problem in namelist PARM01'
204 CALL PRINT_ERROR( msgBuf , 1)
205 CALL MODELDATA_EXAMPLE( myThid )
206 STOP 'ABNORMAL END: S/R INI_PARMS'
207 4 CONTINUE
208 IF ( implicitFreeSurface ) freeSurfFac = 1. _d 0
209 IF ( rigidLid ) freeSurfFac = 0. _d 0
210 IF ( momViscosity ) THEN
211 vfFacMom = 1. _d 0
212 ELSE
213 vfFacMom = 0. _d 0
214 ENDIF
215 IF ( momAdvection ) THEN
216 afFacMom = 1. _d 0
217 ELSE
218 afFacMom = 0. _d 0
219 ENDIF
220 IF ( momForcing ) THEN
221 foFacMom = 1. _d 0
222 ELSE
223 foFacMom = 0. _d 0
224 ENDIF
225 IF ( useCoriolis ) THEN
226 cfFacMom = 1. _d 0
227 ELSE
228 cfFacMom = 0. _d 0
229 ENDIF
230 IF ( momPressureForcing ) THEN
231 pfFacMom = 1. _d 0
232 ELSE
233 pfFacMom = 0. _d 0
234 ENDIF
235 IF ( metricTerms ) THEN
236 mTFacMom = 1. _d 0
237 ELSE
238 mTFacMom = 1. _d 0
239 ENDIF
240
241 IF ( implicitFreeSurface .AND. rigidLid ) THEN
242 WRITE(msgBuf,'(A)')
243 & 'S/R INI_PARMS: Cannot select implicitFreeSurface and rigidLid.'
244 CALL PRINT_ERROR( msgBuf , myThid)
245 STOP 'ABNORMAL END: S/R INI_PARMS'
246 ENDIF
247
248 C-- Elliptic solver parameters
249 cg2dMaxIters = 150
250 cg2dTargetResidual = 1. _d -7
251 cg2dChkResFreq = 1
252 cg2dpcOffDFac = 0.51 _d 0
253 READ(UNIT=iUnit,NML=PARM02,IOSTAT=errIO,err=5)
254 IF ( errIO .GE. 0 ) GOTO 6
255 5 CONTINUE
256 WRITE(msgBuf,'(A)')
257 & 'S/R INI_PARMS'
258 CALL PRINT_ERROR( msgBuf , 1)
259 WRITE(msgBuf,'(A)')
260 & 'Error reading numerical model '
261 CALL PRINT_ERROR( msgBuf , 1)
262 WRITE(msgBuf,'(A)')
263 & 'parameter file "data".'
264 CALL PRINT_ERROR( msgBuf , 1)
265 WRITE(msgBuf,'(A)')
266 & 'Problem in namelist PARM02'
267 CALL PRINT_ERROR( msgBuf , 1)
268 CALL MODELDATA_EXAMPLE( myThid )
269 STOP 'ABNORMAL END: S/R INI_PARMS'
270 6 CONTINUE
271
272 C-- Time stepping parameters
273 startTime = 0.
274 nTimeSteps = 0
275 endTime = 0.
276 nIter0 = 0
277 deltaT = 0.
278 deltaTClock = 0.
279 deltaTtracer = 0.
280 deltaTMom = 0.
281 abEps = 0.01
282 pchkPtFreq = 0.
283 chkPtFreq = 3600.*25
284 dumpFreq = 3600.*100
285 taveFreq = 0.
286 writeStatePrec = precFloat32
287 nCheckLev = 1
288 checkPtSuff(1) = 'ckptA'
289 checkPtSuff(2) = 'ckptB'
290 cAdjFreq = -1. _d 0
291 rCD = -1. _d 0
292 tauCD = 0. _d 0
293 tauThetaClimRelax = 0. _d 0
294 doThetaClimRelax = .FALSE.
295 tauSaltClimRelax = 0. _d 0
296 doSaltClimRelax = .FALSE.
297 periodicExternalForcing = .FALSE.
298 externForcingPeriod = 0.
299 externForcingCycle = 0.
300 READ(UNIT=iUnit,NML=PARM03,IOSTAT=errIO,err=7)
301 IF ( errIO .GE. 0 ) GOTO 8
302 7 CONTINUE
303 WRITE(msgBuf,'(A)')
304 & 'S/R INI_PARMS'
305 CALL PRINT_ERROR( msgBuf , 1)
306 WRITE(msgBuf,'(A)')
307 & 'Error reading numerical model '
308 CALL PRINT_ERROR( msgBuf , 1)
309 WRITE(msgBuf,'(A)')
310 & 'parameter file "data"'
311 CALL PRINT_ERROR( msgBuf , 1)
312 WRITE(msgBuf,'(A)')
313 & 'Problem in namelist PARM03'
314 CALL PRINT_ERROR( msgBuf , 1)
315 CALL MODELDATA_EXAMPLE( myThid )
316 STOP 'ABNORMAL END: S/R INI_PARMS'
317 8 CONTINUE
318 C Process "timestepping" params
319 C o Time step size
320 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
321 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
322 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
323 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
324 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
325 IF ( periodicExternalForcing ) THEN
326 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
327 WRITE(msgBuf,'(A)')
328 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
329 CALL PRINT_ERROR( msgBuf , 1)
330 STOP 'ABNORMAL END: S/R INI_PARMS'
331 ENDIF
332 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
333 & externForcingCycle/externForcingPeriod ) THEN
334 WRITE(msgBuf,'(A)')
335 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
336 CALL PRINT_ERROR( msgBuf , 1)
337 STOP 'ABNORMAL END: S/R INI_PARMS'
338 ENDIF
339 IF ( externForcingCycle.le.externForcingPeriod ) THEN
340 WRITE(msgBuf,'(A)')
341 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
342 CALL PRINT_ERROR( msgBuf , 1)
343 STOP 'ABNORMAL END: S/R INI_PARMS'
344 ENDIF
345 IF ( externForcingPeriod.lt.deltaTclock ) THEN
346 WRITE(msgBuf,'(A)')
347 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
348 CALL PRINT_ERROR( msgBuf , 1)
349 STOP 'ABNORMAL END: S/R INI_PARMS'
350 ENDIF
351 ENDIF
352 C o Convection frequency
353 IF ( cAdjFreq .LT. 0. ) THEN
354 cAdjFreq = deltaTClock
355 ENDIF
356 C o CD coupling
357 IF ( tauCD .EQ. 0. _d 0 ) THEN
358 tauCD = deltaTmom
359 ENDIF
360 IF ( rCD .LT. 0. ) THEN
361 rCD = 1. - deltaTMom/tauCD
362 ENDIF
363 C o Temperature climatology relaxation time scale
364 IF ( tauThetaClimRelax .EQ. 0. _d 0 ) THEN
365 doThetaClimRelax = .FALSE.
366 lambdaThetaClimRelax = 0. _d 0
367 ELSE
368 doThetaClimRelax = .TRUE.
369 lambdaThetaClimRelax = 1./tauThetaClimRelax
370 ENDIF
371 C o Salinity climatology relaxation time scale
372 IF ( tauSaltClimRelax .EQ. 0. _d 0 ) THEN
373 doSaltClimRelax = .FALSE.
374 lambdaSaltClimRelax = 0. _d 0
375 ELSE
376 doSaltClimRelax = .TRUE.
377 lambdaSaltClimRelax = 1./tauSaltClimRelax
378 ENDIF
379 C o Time step count
380 IF ( endTime .NE. 0 ) THEN
381 IF ( deltaTClock .NE. 0 ) nTimeSteps =
382 & INT((endTime-startTime)/deltaTClock)
383 ENDIF
384 C o Finish time
385 IF ( endTime .EQ. 0. ) endTime = FLOAT(nTimeSteps)*deltaTClock
386
387 C o If taveFreq is finite, then we must make sure the diagnostics
388 C code is being compiled
389 #ifndef ALLOW_DIAGNOSTICS
390 IF (taveFreq.NE.0.) THEN
391 WRITE(msgBuf,'(A)')
392 & 'S/R INI_PARMS: taveFreq <> 0 but you have'
393 CALL PRINT_ERROR( msgBuf , 1)
394 WRITE(msgBuf,'(A)')
395 & 'not compiled the model with the diagnostics routines.'
396 CALL PRINT_ERROR( msgBuf , 1)
397 WRITE(msgBuf,'(A)')
398 & 'Re-compile with: #define ALLOW_DIAGNOSTICS or -DALLOW_DIAGNOSTICS'
399 CALL PRINT_ERROR( msgBuf , 1)
400 STOP 'ABNORMAL END: S/R INI_PARMS'
401 ENDIF
402 #endif
403
404 C-- Grid parameters
405 C In cartesian coords distances are in metres
406 usingCartesianGrid = .TRUE.
407 DO K =1,Nz
408 delZ(K) = 100. _d 0
409 ENDDO
410 dxSpacing = 20. _d 0 * 1000. _d 0
411 dySpacing = 20. _d 0 * 1000. _d 0
412 DO i=1,Nx
413 delX(i) = dxSpacing
414 ENDDO
415 DO j=1,Ny
416 delY(j) = dySpacing
417 ENDDO
418 C In spherical polar distances are in degrees
419 usingSphericalPolarGrid = .FALSE.
420 phiMin = -5.0
421 thetaMin = 0.
422 rSphere = 6370. * 1. _d 3
423 rRsphere = 1. _d 0/rSphere
424 IF ( usingSphericalPolarGrid ) THEN
425 dxSpacing = 1.
426 dySpacing = 1.
427 DO I=1,Nx
428 delX(I) = dxSpacing
429 ENDDO
430 DO J=1,Ny
431 delY(J) = dySpacing
432 ENDDO
433 ENDIF
434
435 READ(UNIT=iUnit,NML=PARM04,IOSTAT=errIO,err=9)
436 IF ( errIO .GE. 0 ) GOTO 10
437 9 CONTINUE
438 WRITE(msgBuf,'(A)')
439 & 'S/R INI_PARMS'
440 CALL PRINT_ERROR( msgBuf , 1)
441 WRITE(msgBuf,'(A)')
442 & 'Error reading numerical model '
443 CALL PRINT_ERROR( msgBuf , 1)
444 WRITE(msgBuf,'(A)')
445 & 'parameter file "data"'
446 CALL PRINT_ERROR( msgBuf , 1)
447 WRITE(msgBuf,'(A)')
448 & 'Problem in namelist PARM04'
449 CALL PRINT_ERROR( msgBuf , 1)
450 CALL MODELDATA_EXAMPLE( myThid )
451 STOP 'ABNORMAL END: S/R INI_PARMS'
452 10 CONTINUE
453
454 IF ( rSphere .NE. 0 ) THEN
455 rRSphere = 1. _d 0/rSphere
456 ELSE
457 rRSphere = 0.
458 ENDIF
459
460 C Initialize EOS coefficients (3rd order polynomial)
461 IF (eostype.eq.'POLY3') THEN
462 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
463 READ(37,*) I
464 IF (I.NE.Nz) THEN
465 WRITE(0,*) 'ini_parms: attempt to read POLY3.COEFFS failed'
466 WRITE(0,*) ' because bad # of levels in data'
467 STOP 'Bad data in POLY3.COEFFS'
468 ENDIF
469 READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nz)
470 DO K=1,Nz
471 READ(37,*) (eosC(I,K),I=1,9)
472 write(0,'(i3,13f8.3)') K,eosRefT(K),eosRefS(K),eosSig0(K),
473 & (eosC(I,K),I=1,9)
474 ENDDO
475 CLOSE(37)
476 ENDIF
477
478 goptCount = 0
479 IF ( usingCartesianGrid ) goptCount = goptCount+1
480 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
481 IF ( goptCount .NE. 1 ) THEN
482 WRITE(msgBuf,'(A)')
483 & 'S/R INI_PARMS: More than one coordinate system requested'
484 CALL PRINT_ERROR( msgBuf , myThid)
485 STOP 'ABNORMAL END: S/R INI_PARMS'
486 ENDIF
487
488 IF ( usingCartesianGrid ) THEN
489 usingSphericalPolarMterms = .FALSE.
490 metricTerms = .FALSE.
491 mTFacMom = 0
492 useBetaPlaneF = .TRUE.
493 ENDIF
494 IF ( usingSphericalPolarGrid ) THEN
495 useConstantF = .FALSE.
496 useBetaPlaneF = .FALSE.
497 useSphereF = .TRUE.
498 omega = 2. _d 0 * PI / ( 3600. _d 0 * 24. _d 0 )
499 usingSphericalPolarMterms = .TRUE.
500 metricTerms = .TRUE.
501 mTFacMom = 1
502 ENDIF
503
504 C-- Input files
505 bathyFile = ' '
506 hydrogSaltFile = ' '
507 hydrogThetaFile = ' '
508 zonalWindFile = ' '
509 meridWindFile = ' '
510 thetaClimFile = ' '
511 saltClimFile = ' '
512 READ(UNIT=iUnit,NML=PARM05,IOSTAT=errIO,err=11)
513 IF ( errIO .GE. 0 ) GOTO 12
514 11 CONTINUE
515 WRITE(msgBuf,'(A)')
516 & 'S/R INI_PARMS'
517 CALL PRINT_ERROR( msgBuf , 1)
518 WRITE(msgBuf,'(A)')
519 & 'Error reading numerical model '
520 CALL PRINT_ERROR( msgBuf , 1)
521 WRITE(msgBuf,'(A)')
522 & 'parameter file "data"'
523 CALL PRINT_ERROR( msgBuf , 1)
524 WRITE(msgBuf,'(A)')
525 & 'Problem in namelist PARM05'
526 CALL PRINT_ERROR( msgBuf , 1)
527 CALL MODELDATA_EXAMPLE( myThid )
528 STOP 'ABNORMAL END: S/R INI_PARMS'
529 12 CONTINUE
530
531 _END_MASTER(myThid)
532
533 C-- Everyone else must wait for the parameters to be loaded
534 _BARRIER
535 C
536
537 RETURN
538 END
539

  ViewVC Help
Powered by ViewVC 1.1.22