/[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.26 - (show annotations) (download)
Sat Aug 22 17:51:08 1998 UTC (25 years, 9 months ago) by cnh
Branch: MAIN
Changes since 1.25: +9 -9 lines
Isomorphism consistency changes

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_parms.F,v 1.25 1998/08/19 13:37:52 cnh 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,GMmaxval,
59 & tRef, sRef, eosType,
60 & momViscosity, momAdvection, momForcing, useCoriolis,
61 & momPressureForcing, metricTerms,
62 & tempDiffusion, tempAdvection, tempForcing,
63 & saltDiffusion, saltAdvection, saltForcing,
64 & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz,
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,Nr
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 GMmaxval=2500.d0
175 tAlpha=2.d-4
176 sBeta=7.4d-4
177 eosType='LINEAR'
178 implicitFreeSurface = .TRUE.
179 rigidLid = .FALSE.
180 freeSurfFac = 1. _d 0
181 hFacMin = 0. _d 0
182 hFacMinDz = 0. _d 0
183 momViscosity = .TRUE.
184 momAdvection = .TRUE.
185 momForcing = .TRUE.
186 useCoriolis = .TRUE.
187 momPressureForcing = .TRUE.
188 momStepping = .TRUE.
189 tempStepping = .TRUE.
190 saltStepping = .TRUE.
191 metricTerms = .TRUE.
192 implicitDiffusion = .FALSE.
193 READ(UNIT=iUnit,NML=PARM01,IOSTAT=errIO,err=3)
194 IF ( errIO .GE. 0 ) GOTO 4
195 3 CONTINUE
196 WRITE(msgBuf,'(A)')
197 & 'S/R INI_PARMS'
198 CALL PRINT_ERROR( msgBuf , 1)
199 WRITE(msgBuf,'(A)')
200 & 'Error reading numerical model '
201 CALL PRINT_ERROR( msgBuf , 1)
202 WRITE(msgBuf,'(A)')
203 & 'parameter file "data"'
204 CALL PRINT_ERROR( msgBuf , 1)
205 WRITE(msgBuf,'(A)')
206 & 'Problem in namelist PARM01'
207 CALL PRINT_ERROR( msgBuf , 1)
208 CALL MODELDATA_EXAMPLE( myThid )
209 STOP 'ABNORMAL END: S/R INI_PARMS'
210 4 CONTINUE
211 IF ( implicitFreeSurface ) freeSurfFac = 1. _d 0
212 IF ( rigidLid ) freeSurfFac = 0. _d 0
213 IF ( momViscosity ) THEN
214 vfFacMom = 1. _d 0
215 ELSE
216 vfFacMom = 0. _d 0
217 ENDIF
218 IF ( momAdvection ) THEN
219 afFacMom = 1. _d 0
220 ELSE
221 afFacMom = 0. _d 0
222 ENDIF
223 IF ( momForcing ) THEN
224 foFacMom = 1. _d 0
225 ELSE
226 foFacMom = 0. _d 0
227 ENDIF
228 IF ( useCoriolis ) THEN
229 cfFacMom = 1. _d 0
230 ELSE
231 cfFacMom = 0. _d 0
232 ENDIF
233 IF ( momPressureForcing ) THEN
234 pfFacMom = 1. _d 0
235 ELSE
236 pfFacMom = 0. _d 0
237 ENDIF
238 IF ( metricTerms ) THEN
239 mTFacMom = 1. _d 0
240 ELSE
241 mTFacMom = 1. _d 0
242 ENDIF
243
244 IF ( implicitFreeSurface .AND. rigidLid ) THEN
245 WRITE(msgBuf,'(A)')
246 & 'S/R INI_PARMS: Cannot select implicitFreeSurface and rigidLid.'
247 CALL PRINT_ERROR( msgBuf , myThid)
248 STOP 'ABNORMAL END: S/R INI_PARMS'
249 ENDIF
250
251 C-- Elliptic solver parameters
252 cg2dMaxIters = 150
253 cg2dTargetResidual = 1. _d -7
254 cg2dChkResFreq = 1
255 cg2dpcOffDFac = 0.51 _d 0
256 READ(UNIT=iUnit,NML=PARM02,IOSTAT=errIO,err=5)
257 IF ( errIO .GE. 0 ) GOTO 6
258 5 CONTINUE
259 WRITE(msgBuf,'(A)')
260 & 'S/R INI_PARMS'
261 CALL PRINT_ERROR( msgBuf , 1)
262 WRITE(msgBuf,'(A)')
263 & 'Error reading numerical model '
264 CALL PRINT_ERROR( msgBuf , 1)
265 WRITE(msgBuf,'(A)')
266 & 'parameter file "data".'
267 CALL PRINT_ERROR( msgBuf , 1)
268 WRITE(msgBuf,'(A)')
269 & 'Problem in namelist PARM02'
270 CALL PRINT_ERROR( msgBuf , 1)
271 CALL MODELDATA_EXAMPLE( myThid )
272 STOP 'ABNORMAL END: S/R INI_PARMS'
273 6 CONTINUE
274
275 C-- Time stepping parameters
276 startTime = 0.
277 nTimeSteps = 0
278 endTime = 0.
279 nIter0 = 0
280 deltaT = 0.
281 deltaTClock = 0.
282 deltaTtracer = 0.
283 deltaTMom = 0.
284 abEps = 0.01
285 pchkPtFreq = 0.
286 chkPtFreq = 3600.*25
287 dumpFreq = 3600.*100
288 taveFreq = 0.
289 writeStatePrec = precFloat32
290 nCheckLev = 1
291 checkPtSuff(1) = 'ckptA'
292 checkPtSuff(2) = 'ckptB'
293 cAdjFreq = -1. _d 0
294 rCD = -1. _d 0
295 tauCD = 0. _d 0
296 tauThetaClimRelax = 0. _d 0
297 doThetaClimRelax = .FALSE.
298 tauSaltClimRelax = 0. _d 0
299 doSaltClimRelax = .FALSE.
300 periodicExternalForcing = .FALSE.
301 externForcingPeriod = 0.
302 externForcingCycle = 0.
303 READ(UNIT=iUnit,NML=PARM03,IOSTAT=errIO,err=7)
304 IF ( errIO .GE. 0 ) GOTO 8
305 7 CONTINUE
306 WRITE(msgBuf,'(A)')
307 & 'S/R INI_PARMS'
308 CALL PRINT_ERROR( msgBuf , 1)
309 WRITE(msgBuf,'(A)')
310 & 'Error reading numerical model '
311 CALL PRINT_ERROR( msgBuf , 1)
312 WRITE(msgBuf,'(A)')
313 & 'parameter file "data"'
314 CALL PRINT_ERROR( msgBuf , 1)
315 WRITE(msgBuf,'(A)')
316 & 'Problem in namelist PARM03'
317 CALL PRINT_ERROR( msgBuf , 1)
318 CALL MODELDATA_EXAMPLE( myThid )
319 STOP 'ABNORMAL END: S/R INI_PARMS'
320 8 CONTINUE
321 C Process "timestepping" params
322 C o Time step size
323 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
324 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
325 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
326 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
327 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
328 IF ( periodicExternalForcing ) THEN
329 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
330 WRITE(msgBuf,'(A)')
331 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
332 CALL PRINT_ERROR( msgBuf , 1)
333 STOP 'ABNORMAL END: S/R INI_PARMS'
334 ENDIF
335 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
336 & externForcingCycle/externForcingPeriod ) THEN
337 WRITE(msgBuf,'(A)')
338 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
339 CALL PRINT_ERROR( msgBuf , 1)
340 STOP 'ABNORMAL END: S/R INI_PARMS'
341 ENDIF
342 IF ( externForcingCycle.le.externForcingPeriod ) THEN
343 WRITE(msgBuf,'(A)')
344 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
345 CALL PRINT_ERROR( msgBuf , 1)
346 STOP 'ABNORMAL END: S/R INI_PARMS'
347 ENDIF
348 IF ( externForcingPeriod.lt.deltaTclock ) THEN
349 WRITE(msgBuf,'(A)')
350 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
351 CALL PRINT_ERROR( msgBuf , 1)
352 STOP 'ABNORMAL END: S/R INI_PARMS'
353 ENDIF
354 ENDIF
355 C o Convection frequency
356 IF ( cAdjFreq .LT. 0. ) THEN
357 cAdjFreq = deltaTClock
358 ENDIF
359 C o CD coupling
360 IF ( tauCD .EQ. 0. _d 0 ) THEN
361 tauCD = deltaTmom
362 ENDIF
363 IF ( rCD .LT. 0. ) THEN
364 rCD = 1. - deltaTMom/tauCD
365 ENDIF
366 C o Temperature climatology relaxation time scale
367 IF ( tauThetaClimRelax .EQ. 0. _d 0 ) THEN
368 doThetaClimRelax = .FALSE.
369 lambdaThetaClimRelax = 0. _d 0
370 ELSE
371 doThetaClimRelax = .TRUE.
372 lambdaThetaClimRelax = 1./tauThetaClimRelax
373 ENDIF
374 C o Salinity climatology relaxation time scale
375 IF ( tauSaltClimRelax .EQ. 0. _d 0 ) THEN
376 doSaltClimRelax = .FALSE.
377 lambdaSaltClimRelax = 0. _d 0
378 ELSE
379 doSaltClimRelax = .TRUE.
380 lambdaSaltClimRelax = 1./tauSaltClimRelax
381 ENDIF
382 C o Time step count
383 IF ( endTime .NE. 0 ) THEN
384 IF ( deltaTClock .NE. 0 ) nTimeSteps =
385 & INT((endTime-startTime)/deltaTClock)
386 ENDIF
387 C o Finish time
388 IF ( endTime .EQ. 0. ) endTime = FLOAT(nTimeSteps)*deltaTClock
389
390 C o If taveFreq is finite, then we must make sure the diagnostics
391 C code is being compiled
392 #ifndef ALLOW_DIAGNOSTICS
393 IF (taveFreq.NE.0.) THEN
394 WRITE(msgBuf,'(A)')
395 & 'S/R INI_PARMS: taveFreq <> 0 but you have'
396 CALL PRINT_ERROR( msgBuf , 1)
397 WRITE(msgBuf,'(A)')
398 & 'not compiled the model with the diagnostics routines.'
399 CALL PRINT_ERROR( msgBuf , 1)
400 WRITE(msgBuf,'(A)')
401 & 'Re-compile with: #define ALLOW_DIAGNOSTICS or -DALLOW_DIAGNOSTICS'
402 CALL PRINT_ERROR( msgBuf , 1)
403 STOP 'ABNORMAL END: S/R INI_PARMS'
404 ENDIF
405 #endif
406
407 C-- Grid parameters
408 C In cartesian coords distances are in metres
409 usingCartesianGrid = .TRUE.
410 DO K =1,Nr
411 delZ(K) = 100. _d 0
412 ENDDO
413 dxSpacing = 20. _d 0 * 1000. _d 0
414 dySpacing = 20. _d 0 * 1000. _d 0
415 DO i=1,Nx
416 delX(i) = dxSpacing
417 ENDDO
418 DO j=1,Ny
419 delY(j) = dySpacing
420 ENDDO
421 C In spherical polar distances are in degrees
422 usingSphericalPolarGrid = .FALSE.
423 phiMin = -5.0
424 thetaMin = 0.
425 rSphere = 6370. * 1. _d 3
426 recip_rSphere = 1. _d 0/rSphere
427 IF ( usingSphericalPolarGrid ) THEN
428 dxSpacing = 1.
429 dySpacing = 1.
430 DO I=1,Nx
431 delX(I) = dxSpacing
432 ENDDO
433 DO J=1,Ny
434 delY(J) = dySpacing
435 ENDDO
436 ENDIF
437
438 READ(UNIT=iUnit,NML=PARM04,IOSTAT=errIO,err=9)
439 IF ( errIO .GE. 0 ) GOTO 10
440 9 CONTINUE
441 WRITE(msgBuf,'(A)')
442 & 'S/R INI_PARMS'
443 CALL PRINT_ERROR( msgBuf , 1)
444 WRITE(msgBuf,'(A)')
445 & 'Error reading numerical model '
446 CALL PRINT_ERROR( msgBuf , 1)
447 WRITE(msgBuf,'(A)')
448 & 'parameter file "data"'
449 CALL PRINT_ERROR( msgBuf , 1)
450 WRITE(msgBuf,'(A)')
451 & 'Problem in namelist PARM04'
452 CALL PRINT_ERROR( msgBuf , 1)
453 CALL MODELDATA_EXAMPLE( myThid )
454 STOP 'ABNORMAL END: S/R INI_PARMS'
455 10 CONTINUE
456
457 IF ( rSphere .NE. 0 ) THEN
458 recip_rSphere = 1. _d 0/rSphere
459 ELSE
460 recip_rSphere = 0.
461 ENDIF
462
463 C Initialize EOS coefficients (3rd order polynomial)
464 IF (eostype.eq.'POLY3') THEN
465 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
466 READ(37,*) I
467 IF (I.NE.Nr) THEN
468 WRITE(0,*) 'ini_parms: attempt to read POLY3.COEFFS failed'
469 WRITE(0,*) ' because bad # of levels in data'
470 STOP 'Bad data in POLY3.COEFFS'
471 ENDIF
472 READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nr)
473 DO K=1,Nr
474 READ(37,*) (eosC(I,K),I=1,9)
475 write(0,'(i3,13f8.3)') K,eosRefT(K),eosRefS(K),eosSig0(K),
476 & (eosC(I,K),I=1,9)
477 ENDDO
478 CLOSE(37)
479 ENDIF
480
481 goptCount = 0
482 IF ( usingCartesianGrid ) goptCount = goptCount+1
483 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
484 IF ( goptCount .NE. 1 ) THEN
485 WRITE(msgBuf,'(A)')
486 & 'S/R INI_PARMS: More than one coordinate system requested'
487 CALL PRINT_ERROR( msgBuf , myThid)
488 STOP 'ABNORMAL END: S/R INI_PARMS'
489 ENDIF
490
491 IF ( usingCartesianGrid ) THEN
492 usingSphericalPolarMterms = .FALSE.
493 metricTerms = .FALSE.
494 mTFacMom = 0
495 useBetaPlaneF = .TRUE.
496 ENDIF
497 IF ( usingSphericalPolarGrid ) THEN
498 useConstantF = .FALSE.
499 useBetaPlaneF = .FALSE.
500 useSphereF = .TRUE.
501 omega = 2. _d 0 * PI / ( 3600. _d 0 * 24. _d 0 )
502 usingSphericalPolarMterms = .TRUE.
503 metricTerms = .TRUE.
504 mTFacMom = 1
505 ENDIF
506
507 C-- Input files
508 bathyFile = ' '
509 hydrogSaltFile = ' '
510 hydrogThetaFile = ' '
511 zonalWindFile = ' '
512 meridWindFile = ' '
513 thetaClimFile = ' '
514 saltClimFile = ' '
515 READ(UNIT=iUnit,NML=PARM05,IOSTAT=errIO,err=11)
516 IF ( errIO .GE. 0 ) GOTO 12
517 11 CONTINUE
518 WRITE(msgBuf,'(A)')
519 & 'S/R INI_PARMS'
520 CALL PRINT_ERROR( msgBuf , 1)
521 WRITE(msgBuf,'(A)')
522 & 'Error reading numerical model '
523 CALL PRINT_ERROR( msgBuf , 1)
524 WRITE(msgBuf,'(A)')
525 & 'parameter file "data"'
526 CALL PRINT_ERROR( msgBuf , 1)
527 WRITE(msgBuf,'(A)')
528 & 'Problem in namelist PARM05'
529 CALL PRINT_ERROR( msgBuf , 1)
530 CALL MODELDATA_EXAMPLE( myThid )
531 STOP 'ABNORMAL END: S/R INI_PARMS'
532 12 CONTINUE
533
534 C
535 CLOSE(iUnit)
536
537 _END_MASTER(myThid)
538
539 C-- Everyone else must wait for the parameters to be loaded
540 _BARRIER
541 C
542
543 RETURN
544 END
545

  ViewVC Help
Powered by ViewVC 1.1.22