/[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.157 - (show annotations) (download)
Thu Jun 9 15:53:55 2005 UTC (18 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57i_post
Changes since 1.156: +2 -2 lines
add flag for momentum vertical advection (upwindShear)

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_parms.F,v 1.156 2005/05/31 14:49:38 adcroft Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
8 CBOP
9 C !ROUTINE: INI_PARMS
10 C !INTERFACE:
11 SUBROUTINE INI_PARMS( myThid )
12
13 C !DESCRIPTION:
14 C Routine to set model "parameters".
15
16 C !USES:
17 IMPLICIT NONE
18 #include "SIZE.h"
19 #include "EEPARAMS.h"
20 #include "PARAMS.h"
21 #include "GRID.h"
22 #include "EOS.h"
23
24 C !INPUT/OUTPUT PARAMETERS:
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 deltaTtracer :: Timestep for tracer equations ( s )
34 C tmp4delX,tmp8delX - temporary arrays to read in delX
35 C tmp4delY,tmp8delY - temporary arrays to read in delY
36 C goptCount - Used to count the nuber of grid options
37 C (only one is allowed! )
38 C msgBuf - Informational/error meesage buffer
39 C errIO - IO error flag
40 C iUnit - Work variable for IO unit number
41 C record - Work variable for IO buffer
42 C K, I, J - Loop counters
43 C xxxDefault - Default value for variable xxx
44 _RL dxSpacing
45 _RL dySpacing
46 _RL deltaTtracer
47 REAL*4 tmp4delX(Nx), tmp4delY(Ny)
48 REAL*8 tmp8delX(Nx), tmp8delY(Ny)
49 CHARACTER*(MAX_LEN_FNAM) delXfile
50 CHARACTER*(MAX_LEN_FNAM) delYfile
51 CHARACTER*(MAX_LEN_MBUF) msgBuf
52 CHARACTER*(MAX_LEN_PREC) record
53 INTEGER goptCount
54 INTEGER K, i, j, IL, iUnit
55 INTEGER errIO
56 INTEGER IFNBLNK
57 EXTERNAL IFNBLNK
58 INTEGER ILNBLNK
59 EXTERNAL ILNBLNK
60 C Default values for variables which have vertical coordinate system
61 C dependency.
62 _RL viscArDefault
63 _RL diffKrTDefault
64 _RL diffKrSDefault
65 _RL hFacMinDrDefault
66 _RL delRDefault(Nr)
67 _RS rkFacDefault
68 C zCoordInputData :: Variables used to select between different coordinate systems.
69 C pCoordInputData :: The vertical coordinate system in the rest of the model is
70 C rCoordInputData :: written in terms of r. In the model "data" file input data can
71 C coordsSet :: be interms of z, p or r.
72 C :: e.g. delZ or delP or delR for the vertical grid spacing.
73 C :: The following rules apply:
74 C :: All parameters must use the same vertical coordinate system.
75 C :: e.g. delZ and viscAz is legal but
76 C :: delZ and viscAr is an error.
77 C :: Similarly specifyinh delZ and delP is an error.
78 C :: zCoord..., pCoord..., rCoord... are used to flag when z, p or r are
79 C :: used. coordsSet counts how many vertical coordinate systems have been
80 C used to specify variables. coordsSet > 1 is an error.
81 C
82
83 LOGICAL zCoordInputData
84 LOGICAL pCoordInputData
85 LOGICAL rCoordInputData
86 INTEGER coordsSet
87 LOGICAL diffKrSet
88
89 C Variables which have vertical coordinate system dependency.
90 C delZ :: Vertical grid spacing ( m ).
91 C delP :: Vertical grid spacing ( Pa ).
92 C viscAz :: Eddy viscosity coeff. for mixing of
93 C momentum vertically ( m^2/s )
94 C viscAp :: Eddy viscosity coeff. for mixing of
95 C momentum vertically ( Pa^2/s )
96 C diffKzT :: Laplacian diffusion coeff. for mixing of
97 C heat vertically ( m^2/s )
98 C diffKpT :: Laplacian diffusion coeff. for mixing of
99 C heat vertically ( Pa^2/s )
100 C diffKzS :: Laplacian diffusion coeff. for mixing of
101 C salt vertically ( m^2/s )
102 C diffKpS :: Laplacian diffusion coeff. for mixing of
103 C salt vertically ( Pa^2/s )
104 _RL delZ(Nr)
105 _RL delP(Nr)
106 _RL viscAz
107 _RL viscAp
108 _RL diffKzT
109 _RL diffKpT
110 _RL diffKrT
111 _RL diffKzS
112 _RL diffKpS
113 _RL diffKrS
114
115 C Retired main data file parameters. Kept here to trap use of old data files.
116 C tracerAdvScheme :: tracer advection scheme (old passive tracer code)
117 C trac_EvPrRn :: tracer conc. in Rain & Evap (old passive tracer code)
118 C saltDiffusion :: diffusion of salinity on/off (flag not used)
119 C tempDiffusion :: diffusion of temperature on/off (flag not used)
120 C zonal_filt_lat :: Moved to package "zonal_filt"
121 C nRetired :: Counter used to trap gracefully namelists containing "retired"
122 C :: parameters. These are parameters that are either no-longer used
123 C or that have moved to a different input file and/or namelist.
124 LOGICAL tempDiffusion, saltDiffusion
125 INTEGER tracerAdvScheme
126 _RL trac_EvPrRn
127 _RL zonal_filt_lat
128 INTEGER nRetired
129
130 C-- Continuous equation parameters
131 NAMELIST /PARM01/
132 & gravitySign, nh_Am2,
133 & gravity, gBaro, rhonil, tAlpha, sBeta,
134 & f0, beta, omega, rotationPeriod,
135 & viscAh, viscAhW, viscAhMax,
136 & viscAhGrid, viscAhGridMax, viscAhGridMin,
137 & viscC2leith, viscC4leith,
138 & useFullLeith, useAnisotropicViscAGridMax,
139 & viscC2leithD, viscC4leithD, viscC2Smag,
140 & viscAhD, viscAhZ, viscA4D, viscA4Z,
141 & viscA4, viscA4W,
142 & viscA4Max, viscA4Grid, viscA4GridMax, viscA4GridMin,
143 & viscAz, cosPower, viscAstrain, viscAtension,
144 & diffKhT, diffKzT, diffK4T,
145 & diffKhS, diffKzS, diffK4S,
146 & tRef, sRef, eosType, integr_GeoPot, selectFindRoSurf,
147 & atm_Cp, atm_Rd, atm_Rq,
148 & no_slip_sides,no_slip_bottom,
149 & momViscosity, momAdvection, momForcing, useCoriolis,
150 & useConstantF,
151 & momPressureForcing, metricTerms, vectorInvariantMomentum,
152 & tempDiffusion, tempAdvection, tempForcing,
153 & saltDiffusion, saltAdvection, saltForcing,
154 & implicSurfPress, implicDiv2DFlow,
155 & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz,
156 & exactConserv,uniformLin_PhiSurf,nonlinFreeSurf,hFacInf,hFacSup,
157 & select_rStar,
158 & staggerTimeStep,
159 & tempStepping, saltStepping, momStepping,
160 & implicitDiffusion, implicitViscosity,
161 & tempImplVertAdv, saltImplVertAdv, momImplVertAdv,
162 & viscAr, diffKrT, diffKrS, diffKrNrT, diffKrNrS, hFacMinDr,
163 & viscAp, diffKpT, diffKpS, hFacMinDp,
164 & diffKrBL79surf, diffKrBL79deep, diffKrBL79scl, diffKrBL79Ho,
165 & rhoConst, rhoConstFresh, buoyancyRelation, HeatCapacity_Cp,
166 & writeBinaryPrec, readBinaryPrec, writeStatePrec,
167 & nonHydrostatic, quasiHydrostatic, globalFiles, useSingleCpuIO,
168 & allowFreezing, useOldFreezing, ivdc_kappa,
169 & bottomDragLinear,bottomDragQuadratic,
170 & usePickupBeforeC35, usePickupBeforeC54, debugMode, debugLevel,
171 & tempAdvScheme, tempVertAdvScheme,
172 & saltAdvScheme, saltVertAdvScheme, tracerAdvScheme,
173 & multiDimAdvection, useEnergyConservingCoriolis,
174 & useCDscheme, useJamartWetPoints, useJamartMomAdv, useNHMTerms,
175 & SadournyCoriolis, upwindVorticity, highOrderVorticity,
176 & useAbsVorticity, upwindShear,
177 & useRealFreshWaterFlux, convertFW2Salt,
178 & temp_EvPrRn, salt_EvPrRn, trac_EvPrRn,
179 & zonal_filt_lat,
180 & inAdExact
181
182 C-- Elliptic solver parameters
183 NAMELIST /PARM02/
184 & cg2dMaxIters, cg2dChkResFreq, cg2dTargetResidual,
185 & cg2dTargetResWunit, cg2dpcOffDFac, cg2dPreCondFreq,
186 & cg3dMaxIters, cg3dChkResFreq, cg3dTargetResidual
187
188 C-- Time stepping parammeters
189 NAMELIST /PARM03/
190 & nIter0, nTimeSteps, nEndIter, pickupSuff,
191 & deltaT, deltaTmom, deltaTtracer, dTtracerLev, deltaTfreesurf,
192 & forcing_In_AB, abEps, alph_AB, beta_AB, startFromPickupAB2,
193 & tauCD, rCD,
194 & baseTime, startTime, endTime, chkPtFreq,
195 & dumpFreq, adjDumpFreq, taveFreq, tave_lastIter, deltaTClock,
196 & diagFreq, monitorFreq, adjMonitorFreq, pChkPtFreq, cAdjFreq,
197 & outputTypesInclusive,
198 & tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax,
199 & tauTr1ClimRelax,
200 & periodicExternalForcing, externForcingPeriod, externForcingCycle,
201 & calendarDumps
202
203 C-- Gridding parameters
204 NAMELIST /PARM04/
205 & usingCartesianGrid, dxSpacing, dySpacing, delX, delY, delZ,
206 & usingSphericalPolarGrid, phiMin, thetaMin, rSphere,
207 & usingCurvilinearGrid,usingCylindricalGrid,
208 & delP, delR, rkFac, Ro_SeaLevel, groundAtK1, delRc,
209 & delXfile, delYfile
210
211 C-- Input files
212 NAMELIST /PARM05/
213 & bathyFile, topoFile, hydrogThetaFile, hydrogSaltFile,
214 & zonalWindFile, meridWindFile,
215 & thetaClimFile, saltClimFile,
216 & surfQfile, surfQnetFile, surfQswFile, EmPmRfile, saltFluxFile,
217 & lambdaThetaFile, lambdaSaltFile,
218 & uVelInitFile, vVelInitFile, pSurfInitFile,
219 & dQdTFile, ploadFile,tCylIn,tCylOut,
220 & eddyTauxFile, eddyTauyFile,
221 & mdsioLocalDir,
222 & the_run_name
223 CEOP
224
225 C
226 _BEGIN_MASTER(myThid)
227
228 C Defaults values for input parameters
229 CALL SET_DEFAULTS(
230 O viscArDefault, diffKrTDefault, diffKrSDefault,
231 O hFacMinDrDefault, delRdefault, rkFacDefault,
232 I myThid )
233
234 C-- Initialise "which vertical coordinate system used" flags.
235 zCoordInputData = .FALSE.
236 pCoordInputData = .FALSE.
237 rCoordInputData = .FALSE.
238 coordsSet = 0
239
240 C-- Initialise retired parameters to unlikely value
241 nRetired = 0
242 tempDiffusion = .TRUE.
243 saltDiffusion = .TRUE.
244 tracerAdvScheme = UNSET_I
245 trac_EvPrRn = UNSET_RL
246 zonal_filt_lat = UNSET_RL
247 gravitySign = UNSET_RL
248
249 C-- Open the parameter file
250 OPEN(UNIT=scrUnit1,STATUS='SCRATCH')
251 OPEN(UNIT=scrUnit2,STATUS='SCRATCH')
252 OPEN(UNIT=modelDataUnit,FILE='data',STATUS='OLD',
253 & IOSTAT=errIO)
254 IF ( errIO .LT. 0 ) THEN
255 WRITE(msgBuf,'(A)')
256 & 'S/R INI_PARMS'
257 CALL PRINT_ERROR( msgBuf , 1)
258 WRITE(msgBuf,'(A)')
259 & 'Unable to open model parameter'
260 CALL PRINT_ERROR( msgBuf , 1)
261 WRITE(msgBuf,'(A)')
262 & 'file "data"'
263 CALL PRINT_ERROR( msgBuf , 1)
264 CALL MODELDATA_EXAMPLE( myThid )
265 STOP 'ABNORMAL END: S/R INI_PARMS'
266 ENDIF
267
268 DO WHILE ( .TRUE. )
269 READ(modelDataUnit,FMT='(A)',END=1001) RECORD
270 IL = MAX(ILNBLNK(RECORD),1)
271 IF ( RECORD(1:1) .NE. commentCharacter ) THEN
272 CALL NML_SET_TERMINATOR( RECORD )
273 WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL)
274 ENDIF
275 WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL)
276 ENDDO
277 1001 CONTINUE
278 CLOSE(modelDataUnit)
279
280 C-- Report contents of model parameter file
281 WRITE(msgBuf,'(A)')
282 &'// ======================================================='
283 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
284 & SQUEEZE_RIGHT , 1)
285 WRITE(msgBuf,'(A)') '// Model parameter file "data"'
286 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
287 & SQUEEZE_RIGHT , 1)
288 WRITE(msgBuf,'(A)')
289 &'// ======================================================='
290 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
291 & SQUEEZE_RIGHT , 1)
292 iUnit = scrUnit2
293 REWIND(iUnit)
294 DO WHILE ( .TRUE. )
295 READ(UNIT=iUnit,FMT='(A)',END=2001) RECORD
296 IL = MAX(ILNBLNK(RECORD),1)
297 WRITE(msgBuf,'(A,A)') '>',RECORD(:IL)
298 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
299 & SQUEEZE_RIGHT , 1)
300 ENDDO
301 2001 CONTINUE
302 CLOSE(iUnit)
303 WRITE(msgBuf,'(A)') ' '
304 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
305 & SQUEEZE_RIGHT , 1)
306
307
308 C-- Read settings from model parameter file "data".
309 iUnit = scrUnit1
310 REWIND(iUnit)
311
312 C-- Set default "physical" parameters
313 viscAhW = UNSET_RL
314 viscA4W = UNSET_RL
315 viscAhD = UNSET_RL
316 viscAhZ = UNSET_RL
317 viscA4D = UNSET_RL
318 viscA4Z = UNSET_RL
319 viscAz = UNSET_RL
320 viscAr = UNSET_RL
321 viscAp = UNSET_RL
322 diffKzT = UNSET_RL
323 diffKpT = UNSET_RL
324 diffKrT = UNSET_RL
325 diffKzS = UNSET_RL
326 diffKpS = UNSET_RL
327 diffKrS = UNSET_RL
328 DO k=1,Nr
329 diffKrNrT(k) = UNSET_RL
330 diffKrNrS(k) = UNSET_RL
331 ENDDO
332 gBaro = UNSET_RL
333 rhoConst = UNSET_RL
334 omega = UNSET_RL
335 hFacMinDr = UNSET_RL
336 hFacMinDz = UNSET_RL
337 hFacMinDp = UNSET_RL
338 rhoConstFresh = UNSET_RL
339 convertFW2Salt = UNSET_RL
340 tAlpha = UNSET_RL
341 sBeta = UNSET_RL
342 tempVertAdvScheme = 0
343 saltVertAdvScheme = 0
344 C-- z,p,r coord input switching.
345 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; starts to read PARM01'
346 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
347 & SQUEEZE_RIGHT , 1)
348 READ(UNIT=iUnit,NML=PARM01) !,IOSTAT=errIO)
349 IF ( errIO .LT. 0 ) THEN
350 WRITE(msgBuf,'(A)')
351 & 'S/R INI_PARMS'
352 CALL PRINT_ERROR( msgBuf , 1)
353 WRITE(msgBuf,'(A)')
354 & 'Error reading numerical model '
355 CALL PRINT_ERROR( msgBuf , 1)
356 WRITE(msgBuf,'(A)')
357 & 'parameter file "data"'
358 CALL PRINT_ERROR( msgBuf , 1)
359 WRITE(msgBuf,'(A)')
360 & 'Problem in namelist PARM01'
361 CALL PRINT_ERROR( msgBuf , 1)
362 CALL MODELDATA_EXAMPLE( myThid )
363 STOP 'ABNORMAL END: S/R INI_PARMS'
364 ELSE
365 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; read PARM01 : OK'
366 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
367 & SQUEEZE_RIGHT , 1)
368 ENDIF
369
370 C- set the type of vertical coordinate and type of fluid
371 C according to buoyancyRelation
372 usingPCoords = .FALSE.
373 usingZCoords = .FALSE.
374 fluidIsAir = .FALSE.
375 fluidIsWater = .FALSE.
376 IF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
377 usingPCoords = .TRUE.
378 fluidIsAir = .TRUE.
379 ELSEIF ( buoyancyRelation.EQ.'OCEANICP') THEN
380 usingPCoords = .TRUE.
381 fluidIsWater = .TRUE.
382 ELSEIF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
383 usingZCoords = .TRUE.
384 fluidIsWater = .TRUE.
385 ELSE
386 WRITE(msgBuf,'(2A)') 'S/R INI_PARMS:',
387 & ' Bad value of buoyancyRelation '
388 CALL PRINT_ERROR( msgBuf , myThid)
389 STOP 'ABNORMAL END: S/R INI_PARMS'
390 ENDIF
391
392 IF ( implicitFreeSurface ) freeSurfFac = 1.D0
393 IF ( rigidLid ) freeSurfFac = 0.D0
394 IF ( gBaro .EQ. UNSET_RL ) gBaro=gravity
395 IF ( rhoConst .EQ. UNSET_RL ) rhoConst=rhoNil
396 IF ( rhoConstFresh .EQ. UNSET_RL ) rhoConstFresh=rhoConst
397 IF ( omega .EQ. UNSET_RL ) THEN
398 omega = 0. _d 0
399 IF ( rotationPeriod .NE. 0. _d 0 )
400 & omega = 2.D0 * PI / rotationPeriod
401 ELSEIF ( omega .EQ. 0. _d 0 ) THEN
402 rotationPeriod = 0. _d 0
403 ELSE
404 rotationPeriod = 2.D0 * PI / omega
405 ENDIF
406 IF (atm_Rd .EQ. UNSET_RL) THEN
407 atm_Rd = atm_Cp * atm_kappa
408 ELSE
409 atm_kappa = atm_Rd / atm_Cp
410 ENDIF
411 C-- On/Off flags for each terms of the momentum equation
412 nonHydrostatic = momStepping .AND. nonHydrostatic
413 quasiHydrostatic = momStepping .AND. quasiHydrostatic
414 momAdvection = momStepping .AND. momAdvection
415 momViscosity = momStepping .AND. momViscosity
416 momForcing = momStepping .AND. momForcing
417 useCoriolis = momStepping .AND. useCoriolis
418 useCDscheme = momStepping .AND. useCDscheme
419 momPressureForcing= momStepping .AND. momPressureForcing
420 momImplVertAdv = momAdvection .AND. momImplVertAdv
421 implicitViscosity= momViscosity .AND. implicitViscosity
422 C-- Momentum viscosity on/off flag.
423 IF ( momViscosity ) THEN
424 vfFacMom = 1.D0
425 ELSE
426 vfFacMom = 0.D0
427 ENDIF
428 C-- Momentum advection on/off flag.
429 IF ( momAdvection ) THEN
430 afFacMom = 1.D0
431 ELSE
432 afFacMom = 0.D0
433 ENDIF
434 C-- Momentum forcing on/off flag.
435 IF ( momForcing ) THEN
436 foFacMom = 1.D0
437 ELSE
438 foFacMom = 0.D0
439 ENDIF
440 C-- Coriolis term on/off flag.
441 IF ( useCoriolis ) THEN
442 cfFacMom = 1.D0
443 ELSE
444 cfFacMom = 0.D0
445 ENDIF
446 C-- Pressure term on/off flag.
447 IF ( momPressureForcing ) THEN
448 pfFacMom = 1.D0
449 ELSE
450 pfFacMom = 0.D0
451 ENDIF
452 C-- Metric terms on/off flag.
453 IF ( metricTerms ) THEN
454 mTFacMom = 1.D0
455 ELSE
456 mTFacMom = 0.D0
457 ENDIF
458 C-- Non-hydrostatic/quasi-hydrostatic
459 IF (nonHydrostatic.AND.quasiHydrostatic) THEN
460 WRITE(msgBuf,'(A)')
461 & 'Illegal: both nonHydrostatic = quasiHydrostatic = TRUE'
462 CALL PRINT_ERROR( msgBuf , myThid)
463 STOP 'ABNORMAL END: S/R INI_PARMS'
464 ENDIF
465 C-- Advection and Forcing for Temp and salt on/off flags
466 tempAdvection = tempStepping .AND. tempAdvection
467 tempForcing = tempStepping .AND. tempForcing
468 saltAdvection = saltStepping .AND. saltAdvection
469 saltForcing = saltStepping .AND. saltForcing
470 tempImplVertAdv = tempAdvection .AND. tempImplVertAdv
471 saltImplVertAdv = saltAdvection .AND. saltImplVertAdv
472 IF (tempVertAdvScheme.EQ.0) tempVertAdvScheme = tempAdvScheme
473 IF (saltVertAdvScheme.EQ.0) saltVertAdvScheme = saltAdvScheme
474 C-- horizontal viscosity for vertical moments
475 IF ( viscAhW .EQ. UNSET_RL ) viscAhW = viscAh
476 IF ( viscA4W .EQ. UNSET_RL ) viscA4W = viscA4
477 C-- horizontal viscosity (acting on Divergence or Vorticity)
478 IF ( viscAhD .EQ. UNSET_RL ) viscAhD = viscAh
479 IF ( viscAhZ .EQ. UNSET_RL ) viscAhZ = viscAh
480 IF ( viscA4D .EQ. UNSET_RL ) viscA4D = viscA4
481 IF ( viscA4Z .EQ. UNSET_RL ) viscA4Z = viscA4
482 C-- useAnisotropicViscAGridMax requires that viscA*GridMax be set
483 IF ( useAnisotropicViscAGridMax ) THEN
484 IF ( viscAhGridMax .EQ. 1.D21 ) viscAhGridMax = 0.25
485 IF ( viscA4GridMax .EQ. 1.D21 ) viscA4GridMax = 0.25
486 ENDIF
487 C-- z,p,r coord input switching.
488 IF ( viscAz .NE. UNSET_RL ) zCoordInputData = .TRUE.
489 IF ( viscAp .NE. UNSET_RL ) pCoordInputData = .TRUE.
490 IF ( viscAr .NE. UNSET_RL ) rCoordInputData = .TRUE.
491 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAz
492 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscAp
493 IF ( viscAr .EQ. UNSET_RL ) viscAr = viscArDefault
494
495 IF ( diffKzT .NE. UNSET_RL ) zCoordInputData = .TRUE.
496 IF ( diffKpT .NE. UNSET_RL ) pCoordInputData = .TRUE.
497 IF ( diffKrT .NE. UNSET_RL ) rCoordInputData = .TRUE.
498 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKzT
499 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKpT
500 IF ( diffKrT .EQ. UNSET_RL ) diffKrT = diffKrTDefault
501 diffKrSet = .TRUE.
502 DO k=1,Nr
503 IF ( diffKrNrT(k).EQ. UNSET_RL ) diffKrSet = .FALSE.
504 ENDDO
505 IF ( .NOT.diffKrSet ) THEN
506 DO k=1,Nr
507 diffKrNrT(k) = diffKrT
508 ENDDO
509 ELSEIF ( diffKrT.NE.diffKrTDefault ) THEN
510 WRITE(msgBuf,'(2A)') 'S/R INI_PARMS: Cannot set both ',
511 & 'diffKrNrT and diffKrT (or Kp,Kz) in input file data'
512 CALL PRINT_ERROR( msgBuf , myThid)
513 STOP 'ABNORMAL END: S/R INI_PARMS'
514 ENDIF
515
516 IF ( diffKzS .NE. UNSET_RL ) zCoordInputData = .TRUE.
517 IF ( diffKpS .NE. UNSET_RL ) pCoordInputData = .TRUE.
518 IF ( diffKrS .NE. UNSET_RL ) rCoordInputData = .TRUE.
519 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKzS
520 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKpS
521 IF ( diffKrS .EQ. UNSET_RL ) diffKrS = diffKrSDefault
522 diffKrSet = .TRUE.
523 DO k=1,Nr
524 IF ( diffKrNrS(k).EQ. UNSET_RL ) diffKrSet = .FALSE.
525 ENDDO
526 IF ( .NOT.diffKrSet ) THEN
527 DO k=1,Nr
528 diffKrNrS(k) = diffKrS
529 ENDDO
530 ELSEIF ( diffKrS.NE.diffKrSDefault ) THEN
531 WRITE(msgBuf,'(2A)') 'S/R INI_PARMS: Cannot set both ',
532 & 'diffKrNrS and diffKrS (or Kp,Kz) in input file data'
533 CALL PRINT_ERROR( msgBuf , myThid)
534 STOP 'ABNORMAL END: S/R INI_PARMS'
535 ENDIF
536
537 IF ( hFacMinDz .NE. UNSET_RL ) zCoordInputData = .TRUE.
538 IF ( hFacMinDp .NE. UNSET_RL ) pCoordInputData = .TRUE.
539 IF ( hFacMinDr .NE. UNSET_RL ) rCoordInputData = .TRUE.
540 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDz
541 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDp
542 IF ( hFacMinDr .EQ. UNSET_RL ) hFacMinDr = hFacMinDrDefault
543
544 IF (convertFW2Salt.EQ.UNSET_RL) THEN
545 convertFW2Salt = 35.
546 IF (useRealFreshWaterFlux) convertFW2Salt=-1
547 ENDIF
548
549 IF ( ivdc_kappa .NE. 0. .AND. .NOT. implicitDiffusion ) THEN
550 WRITE(msgBuf,'(A,A)')
551 & 'S/R INI_PARMS: To use ivdc_kappa you must enable implicit',
552 & ' vertical diffusion.'
553 CALL PRINT_ERROR( msgBuf , myThid)
554 STOP 'ABNORMAL END: S/R INI_PARMS'
555 ENDIF
556
557 coordsSet = 0
558 IF ( zCoordInputData ) coordsSet = coordsSet + 1
559 IF ( pCoordInputData ) coordsSet = coordsSet + 1
560 IF ( rCoordInputData ) coordsSet = coordsSet + 1
561 IF ( coordsSet .GT. 1 ) THEN
562 WRITE(msgBuf,'(A)')
563 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
564 CALL PRINT_ERROR( msgBuf , myThid)
565 STOP 'ABNORMAL END: S/R INI_PARMS'
566 ENDIF
567 IF ( rhoConst .LE. 0. ) THEN
568 WRITE(msgBuf,'(A)')
569 & 'S/R INI_PARMS: rhoConst must be greater than 0.'
570 CALL PRINT_ERROR( msgBuf , myThid)
571 STOP 'ABNORMAL END: S/R INI_PARMS'
572 ELSE
573 recip_rhoConst = 1.D0 / rhoConst
574 ENDIF
575 IF ( rhoNil .LE. 0. ) THEN
576 WRITE(msgBuf,'(A)')
577 & 'S/R INI_PARMS: rhoNil must be greater than 0.'
578 CALL PRINT_ERROR( msgBuf , myThid)
579 STOP 'ABNORMAL END: S/R INI_PARMS'
580 ELSE
581 recip_rhoNil = 1.D0 / rhoNil
582 ENDIF
583 IF ( HeatCapacity_Cp .LE. 0. ) THEN
584 WRITE(msgBuf,'(A)')
585 & 'S/R INI_PARMS: HeatCapacity_Cp must be greater than 0.'
586 CALL PRINT_ERROR( msgBuf , myThid)
587 STOP 'ABNORMAL END: S/R INI_PARMS'
588 ELSE
589 recip_Cp = 1.D0 / HeatCapacity_Cp
590 ENDIF
591 IF ( gravity .LE. 0. ) THEN
592 WRITE(msgBuf,'(A)')
593 & 'S/R INI_PARMS: gravity must be greater than 0.'
594 CALL PRINT_ERROR( msgBuf , myThid)
595 STOP 'ABNORMAL END: S/R INI_PARMS'
596 ELSE
597 recip_gravity = 1.D0 / gravity
598 ENDIF
599 C This flags are now passed to RW and MDSIO packages in ini_model_io.F
600 C Set globalFiles flag for READ_WRITE_FLD package
601 c CALL SET_WRITE_GLOBAL_FLD( globalFiles )
602 C Set globalFiles flag for READ_WRITE_REC package
603 c CALL SET_WRITE_GLOBAL_REC( globalFiles )
604 C Set globalFiles flag for READ_WRITE_REC package
605 c CALL SET_WRITE_GLOBAL_PICKUP( globalFiles )
606
607 C Check for retired parameters still being used
608 nRetired = 0
609 IF ( zonal_filt_lat .NE. UNSET_RL ) THEN
610 nRetired = nRetired+1
611 WRITE(msgBuf,'(A,A)')
612 & 'S/R INI_PARMS: Paramater "zonal_filt_lat" is',
613 & ' no longer allowed in file "data".'
614 CALL PRINT_ERROR( msgBuf , myThid)
615 WRITE(msgBuf,'(A,A)')
616 & 'S/R INI_PARMS: Paramater "zonal_filt_lat" is',
617 & ' now read from file "data.zonfilt".'
618 CALL PRINT_ERROR( msgBuf , myThid)
619 ENDIF
620 IF ( gravitySign .NE. UNSET_RL ) THEN
621 nRetired = nRetired+1
622 WRITE(msgBuf,'(A,A)')
623 & 'S/R INI_PARMS: "gravitySign" is set according to vertical ',
624 & ' coordinate and is no longer allowed in file "data".'
625 CALL PRINT_ERROR( msgBuf , myThid)
626 ENDIF
627 IF ( tracerAdvScheme .NE. UNSET_I ) THEN
628 nRetired = nRetired+1
629 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: "tracerAdvScheme" ',
630 & '(old passive tracer code) is no longer allowed in file "data"'
631 CALL PRINT_ERROR( msgBuf , myThid)
632 ENDIF
633 IF ( trac_EvPrRn .NE. UNSET_RL ) THEN
634 nRetired = nRetired+1
635 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: "trac_EvPrRn" ',
636 & '(old passive tracer code) is no longer allowed in file "data"'
637 CALL PRINT_ERROR( msgBuf , myThid)
638 ENDIF
639 IF ( .NOT. tempDiffusion ) THEN
640 nRetired = nRetired+1
641 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: "tempDiffusion" ',
642 & 'is no longer allowed in file "data"'
643 CALL PRINT_ERROR( msgBuf , myThid)
644 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: to turn off diffusion',
645 & ' => set diffusivity to zero'
646 CALL PRINT_ERROR( msgBuf , myThid)
647 ENDIF
648 IF ( .NOT. saltDiffusion ) THEN
649 nRetired = nRetired+1
650 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: "saltDiffusion" ',
651 & 'is no longer allowed in file "data"'
652 CALL PRINT_ERROR( msgBuf , myThid)
653 WRITE(msgBuf,'(A,A)') 'S/R INI_PARMS: to turn off diffusion',
654 & ' => set diffusivity to zero'
655 CALL PRINT_ERROR( msgBuf , myThid)
656 ENDIF
657
658 C-- Elliptic solver parameters
659 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; starts to read PARM02'
660 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
661 & SQUEEZE_RIGHT , 1)
662 READ(UNIT=iUnit,NML=PARM02) !,IOSTAT=errIO)
663 IF ( errIO .LT. 0 ) THEN
664 WRITE(msgBuf,'(A)')
665 & 'S/R INI_PARMS'
666 CALL PRINT_ERROR( msgBuf , 1)
667 WRITE(msgBuf,'(A)')
668 & 'Error reading numerical model '
669 CALL PRINT_ERROR( msgBuf , 1)
670 WRITE(msgBuf,'(A)')
671 & 'parameter file "data".'
672 CALL PRINT_ERROR( msgBuf , 1)
673 WRITE(msgBuf,'(A)')
674 & 'Problem in namelist PARM02'
675 CALL PRINT_ERROR( msgBuf , 1)
676 CALL MODELDATA_EXAMPLE( myThid )
677 STOP 'ABNORMAL END: S/R INI_PARMS'
678 ELSE
679 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; read PARM02 : OK'
680 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
681 & SQUEEZE_RIGHT , 1)
682 ENDIF
683
684 C-- Time stepping parameters
685 rCD = -1.D0
686 latBandClimRelax = UNSET_RL
687 deltaTtracer = 0. _d 0
688 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; starts to read PARM03'
689 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
690 & SQUEEZE_RIGHT , 1)
691 READ(UNIT=iUnit,NML=PARM03) !,IOSTAT=errIO)
692 IF ( errIO .LT. 0 ) THEN
693 WRITE(msgBuf,'(A)')
694 & 'S/R INI_PARMS'
695 CALL PRINT_ERROR( msgBuf , 1)
696 WRITE(msgBuf,'(A)')
697 & 'Error reading numerical model '
698 CALL PRINT_ERROR( msgBuf , 1)
699 WRITE(msgBuf,'(A)')
700 & 'parameter file "data"'
701 CALL PRINT_ERROR( msgBuf , 1)
702 WRITE(msgBuf,'(A)')
703 & 'Problem in namelist PARM03'
704 CALL PRINT_ERROR( msgBuf , 1)
705 CALL MODELDATA_EXAMPLE( myThid )
706 STOP 'ABNORMAL END: S/R INI_PARMS'
707 ELSE
708 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; read PARM03 : OK'
709 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
710 & SQUEEZE_RIGHT , 1)
711 ENDIF
712 C Process "timestepping" params
713 C o Time step size
714 IF ( deltaTtracer .NE. dTtracerLev(1) .AND.
715 & deltaTtracer .NE. 0. .AND. dTtracerLev(1) .NE. 0. ) THEN
716 WRITE(msgBuf,'(A)')
717 & 'S/R INI_PARMS: deltaTtracer & dTtracerLev(1) not equal'
718 CALL PRINT_ERROR( msgBuf , myThid)
719 STOP 'ABNORMAL END: S/R INI_PARMS'
720 ELSEIF ( dTtracerLev(1) .NE. 0. ) THEN
721 deltaTtracer = dTtracerLev(1)
722 ENDIF
723 IF ( deltaT .EQ. 0. ) deltaT = deltaTmom
724 IF ( deltaT .EQ. 0. ) deltaT = deltaTtracer
725 IF ( deltaTmom .EQ. 0. ) deltaTmom = deltaT
726 IF ( deltaTtracer .EQ. 0. ) deltaTtracer = deltaT
727 IF ( deltaTClock .EQ. 0. ) deltaTClock = deltaT
728 DO k=1,Nr
729 IF (dTtracerLev(k).EQ.0.) dTtracerLev(k)= deltaTtracer
730 ENDDO
731 C Note that this line should set deltaFreesurf=deltaTmom
732 C but this would change a lot of existing set-ups so we are
733 C obliged to set the default inappropriately.
734 C Be advised that when using asynchronous time stepping
735 C it is better to set deltaTreesurf=deltaTtracer
736 IF ( deltaTfreesurf .EQ. 0. ) deltaTfreesurf = deltaTmom
737 IF ( periodicExternalForcing ) THEN
738 IF ( externForcingCycle*externForcingPeriod .EQ. 0. ) THEN
739 WRITE(msgBuf,'(A)')
740 & 'S/R INI_PARMS: externForcingCycle,externForcingPeriod =0'
741 CALL PRINT_ERROR( msgBuf , 1)
742 STOP 'ABNORMAL END: S/R INI_PARMS'
743 ENDIF
744 IF ( INT(externForcingCycle/externForcingPeriod) .NE.
745 & externForcingCycle/externForcingPeriod ) THEN
746 WRITE(msgBuf,'(A)')
747 & 'S/R INI_PARMS: externForcingCycle <> N*externForcingPeriod'
748 CALL PRINT_ERROR( msgBuf , 1)
749 STOP 'ABNORMAL END: S/R INI_PARMS'
750 ENDIF
751 IF ( externForcingCycle.lt.externForcingPeriod ) THEN
752 WRITE(msgBuf,'(A)')
753 & 'S/R INI_PARMS: externForcingCycle < externForcingPeriod'
754 CALL PRINT_ERROR( msgBuf , 1)
755 STOP 'ABNORMAL END: S/R INI_PARMS'
756 ENDIF
757 IF ( externForcingPeriod.lt.deltaTclock ) THEN
758 WRITE(msgBuf,'(A)')
759 & 'S/R INI_PARMS: externForcingPeriod < deltaTclock'
760 CALL PRINT_ERROR( msgBuf , 1)
761 STOP 'ABNORMAL END: S/R INI_PARMS'
762 ENDIF
763 ENDIF
764 C o Convection frequency
765 IF ( cAdjFreq .LT. 0. ) THEN
766 cAdjFreq = deltaTClock
767 ENDIF
768 IF ( ivdc_kappa .NE. 0. .AND. cAdjFreq .NE. 0. ) THEN
769 WRITE(msgBuf,'(A,A)')
770 & 'S/R INI_PARMS: You have enabled both ivdc_kappa and',
771 & ' convective adjustment.'
772 CALL PRINT_ERROR( msgBuf , myThid)
773 STOP 'ABNORMAL END: S/R INI_PARMS'
774 ENDIF
775 IF (useCDscheme) THEN
776 C o CD coupling (CD scheme):
777 IF ( tauCD .EQ. 0.D0 ) tauCD = deltaTmom
778 IF ( rCD .LT. 0. ) rCD = 1. _d 0 - deltaTMom/tauCD
779 ENDIF
780 C o Temperature climatology relaxation time scale
781 IF ( tauThetaClimRelax .EQ. 0.D0 ) THEN
782 doThetaClimRelax = .FALSE.
783 ELSE
784 doThetaClimRelax = .TRUE.
785 ENDIF
786 C o Salinity climatology relaxation time scale
787 IF ( tauSaltClimRelax .EQ. 0.D0 ) THEN
788 doSaltClimRelax = .FALSE.
789 ELSE
790 doSaltClimRelax = .TRUE.
791 ENDIF
792 C o Tracer 1 climatology relaxation time scale
793 IF ( tauTr1ClimRelax .EQ. 0.D0 ) THEN
794 doTr1ClimRelax = .FALSE.
795 lambdaTr1ClimRelax = 0.D0
796 ELSE
797 doTr1ClimRelax = .TRUE.
798 lambdaTr1ClimRelax = 1./tauTr1ClimRelax
799 ENDIF
800
801 C o Base time
802 IF ( nIter0.NE.0 .AND. startTime.NE.0. .AND. baseTime.EQ.0. )
803 & baseTime = startTime - deltaTClock*float(nIter0)
804 C o Start time
805 IF ( nIter0 .NE. 0 .AND. startTime .EQ. 0. )
806 & startTime = baseTime + deltaTClock*float(nIter0)
807 C o nIter0
808 IF ( nIter0 .EQ. 0 .AND. startTime .NE. baseTime )
809 & nIter0 = NINT( (startTime-baseTime)/deltaTClock )
810
811 C o nTimeSteps 1
812 IF ( nTimeSteps .EQ. 0 .AND. nEndIter .NE. 0 )
813 & nTimeSteps = nEndIter-nIter0
814 C o nTimeSteps 2
815 IF ( nTimeSteps .EQ. 0 .AND. endTime .NE. 0. )
816 & nTimeSteps = NINT((endTime-startTime)/deltaTclock)
817 C o nEndIter 1
818 IF ( nEndIter .EQ. 0 .AND. nTimeSteps .NE. 0 )
819 & nEndIter = nIter0+nTimeSteps
820 C o nEndIter 2
821 IF ( nEndIter .EQ. 0 .AND. endTime .NE. 0. )
822 & nEndIter = NINT((endTime-baseTime)/deltaTclock)
823 C o End Time 1
824 IF ( endTime .EQ. 0. .AND. nTimeSteps .NE. 0 )
825 & endTime = startTime + deltaTClock*float(nTimeSteps)
826 C o End Time 2
827 IF ( endTime .EQ. 0. .AND. nEndIter .NE. 0 )
828 & endTime = baseTime + deltaTClock*float(nEndIter)
829
830 C o Consistent?
831 IF ( startTime .NE. baseTime+deltaTClock*float(nIter0) ) THEN
832 WRITE(msgBuf,'(A)')
833 & 'S/R INI_PARMS: startTime, baseTime and nIter0 are inconsistent'
834 CALL PRINT_ERROR( msgBuf , 1)
835 WRITE(msgBuf,'(A)')
836 & 'S/R INI_PARMS: Perhaps more than two were set at once'
837 CALL PRINT_ERROR( msgBuf , 1)
838 STOP 'ABNORMAL END: S/R INI_PARMS'
839 ENDIF
840 IF ( nEndIter .NE. nIter0+nTimeSteps ) THEN
841 WRITE(msgBuf,'(A)')
842 & 'S/R INI_PARMS: nIter0, nTimeSteps and nEndIter are inconsistent'
843 CALL PRINT_ERROR( msgBuf , 1)
844 WRITE(msgBuf,'(A)')
845 & 'S/R INI_PARMS: Perhaps more than two were set at once'
846 CALL PRINT_ERROR( msgBuf , 1)
847 STOP 'ABNORMAL END: S/R INI_PARMS'
848 ENDIF
849 IF ( nTimeSteps .NE. NINT((endTime-startTime)/deltaTClock) )
850 & THEN
851 WRITE(msgBuf,'(A)')
852 & 'S/R INI_PARMS: both endTime and nTimeSteps have been set'
853 CALL PRINT_ERROR( msgBuf , 1)
854 WRITE(msgBuf,'(A)')
855 & 'S/R INI_PARMS: but are inconsistent'
856 CALL PRINT_ERROR( msgBuf , 1)
857 STOP 'ABNORMAL END: S/R INI_PARMS'
858 ENDIF
859
860 C o Monitor (should also add CPP flag for monitor?)
861 IF (monitorFreq.LT.0.) THEN
862 monitorFreq=0.
863 IF (dumpFreq.NE.0.) monitorFreq=dumpFreq
864 IF (diagFreq.NE.0..AND.diagFreq.LT.monitorFreq)
865 & monitorFreq=diagFreq
866 IF (taveFreq.NE.0..AND.taveFreq.LT.monitorFreq)
867 & monitorFreq=taveFreq
868 IF (chkPtFreq.NE.0..AND.chkPtFreq.LT.monitorFreq)
869 & monitorFreq=chkPtFreq
870 IF (pChkPtFreq.NE.0..AND.pChkPtFreq.LT.monitorFreq)
871 & monitorFreq=pChkPtFreq
872 IF (monitorFreq.EQ.0.) monitorFreq=deltaTclock
873 ENDIF
874
875 C-- Grid parameters
876 C In cartesian coords distances are in metres
877 rkFac = UNSET_RS
878 DO K =1,Nr
879 delZ(K) = UNSET_RL
880 delP(K) = UNSET_RL
881 delR(K) = UNSET_RL
882 ENDDO
883 C In spherical polar distances are in degrees
884 recip_rSphere = 1.D0/rSphere
885 dxSpacing = UNSET_RL
886 dySpacing = UNSET_RL
887 delXfile = ' '
888 delYfile = ' '
889 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; starts to read PARM04'
890 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
891 & SQUEEZE_RIGHT , 1)
892 READ(UNIT=iUnit,NML=PARM04,IOSTAT=errIO)
893 IF ( errIO .LT. 0 ) THEN
894 WRITE(msgBuf,'(A)')
895 & 'S/R INI_PARMS'
896 CALL PRINT_ERROR( msgBuf , 1)
897 WRITE(msgBuf,'(A)')
898 & 'Error reading numerical model '
899 CALL PRINT_ERROR( msgBuf , 1)
900 WRITE(msgBuf,'(A)')
901 & 'parameter file "data"'
902 CALL PRINT_ERROR( msgBuf , 1)
903 WRITE(msgBuf,'(A)')
904 & 'Problem in namelist PARM04'
905 CALL PRINT_ERROR( msgBuf , 1)
906 CALL MODELDATA_EXAMPLE( myThid )
907 STOP 'ABNORMAL END: S/R INI_PARMS'
908 ELSE
909 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; read PARM04 : OK'
910 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
911 & SQUEEZE_RIGHT , 1)
912 ENDIF
913
914 C X coordinate
915 IF ( delXfile .NE. ' ' ) THEN
916 IF ( delX(1) .NE. UNSET_RL .OR. dxSpacing .NE. UNSET_RL ) THEN
917 WRITE(msgBuf,'(A,A)') 'Too many specifications for delX:',
918 & 'Specify only one of delX, dxSpacing or delXfile'
919 CALL PRINT_ERROR( msgBuf , 1)
920 STOP 'ABNORMAL END: S/R INI_PARMS'
921 ELSE
922 _BEGIN_MASTER( myThid )
923 IF (readBinaryPrec.EQ.precFloat32) THEN
924 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
925 & ACCESS='DIRECT',RECL=WORDLENGTH*Nx)
926 READ(37,rec=1) tmp4delX
927 #ifdef _BYTESWAPIO
928 call MDS_BYTESWAPR4( Nx, tmp4delX )
929 #endif
930 CLOSE(37)
931 DO i=1,Nx
932 delX(i) = tmp4delX(i)
933 ENDDO
934 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
935 OPEN(37,FILE=delXfile,STATUS='OLD',FORM='UNFORMATTED',
936 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Nx)
937 READ(37,rec=1) tmp8delX
938 #ifdef _BYTESWAPIO
939 call MDS_BYTESWAPR8( Nx, tmp8delX )
940 #endif
941 CLOSE(37)
942 DO i=1,Nx
943 delX(i) = tmp8delX(i)
944 ENDDO
945 ENDIF
946 _END_MASTER(myThid)
947 ENDIF
948 ENDIF
949 IF ( dxSpacing .NE. UNSET_RL ) THEN
950 DO i=1,Nx
951 delX(i) = dxSpacing
952 ENDDO
953 ENDIF
954 C Y coordinate
955 IF ( delYfile .NE. ' ' ) THEN
956 IF ( delY(1) .NE. UNSET_RL .OR. dySpacing .NE. UNSET_RL ) THEN
957 WRITE(msgBuf,'(A,A)') 'Too many specifications for delY:',
958 & 'Specify only one of delY, dySpacing or delYfile'
959 CALL PRINT_ERROR( msgBuf , 1)
960 STOP 'ABNORMAL END: S/R INI_PARMS'
961 ELSE
962 _BEGIN_MASTER( myThid )
963 IF (readBinaryPrec.EQ.precFloat32) THEN
964 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
965 & ACCESS='DIRECT',RECL=WORDLENGTH*Ny)
966 READ(37,rec=1) tmp4delY
967 #ifdef _BYTESWAPIO
968 call MDS_BYTESWAPR4( Ny, tmp4delY )
969 #endif
970 CLOSE(37)
971 DO j=1,Ny
972 delY(j) = tmp4delY(j)
973 ENDDO
974 ELSEIF (readBinaryPrec.EQ.precFloat64) THEN
975 OPEN(37,FILE=delYfile,STATUS='OLD',FORM='UNFORMATTED',
976 & ACCESS='DIRECT',RECL=WORDLENGTH*2*Ny)
977 READ(37,rec=1) tmp8delY
978 #ifdef _BYTESWAPIO
979 call MDS_BYTESWAPR8( Ny, tmp8delY )
980 #endif
981 CLOSE(37)
982 DO j=1,Ny
983 delY(j) = tmp8delY(j)
984 ENDDO
985 ENDIF
986 _END_MASTER(myThid)
987 ENDIF
988 ENDIF
989 IF ( dySpacing .NE. UNSET_RL ) THEN
990 DO i=1,Ny
991 delY(i) = dySpacing
992 ENDDO
993 ENDIF
994 C
995 IF ( rSphere .NE. 0 ) THEN
996 recip_rSphere = 1.D0/rSphere
997 ELSE
998 recip_rSphere = 0.
999 ENDIF
1000 C-- Check for conflicting grid definitions.
1001 goptCount = 0
1002 IF ( usingCartesianGrid ) goptCount = goptCount+1
1003 IF ( usingSphericalPolarGrid ) goptCount = goptCount+1
1004 IF ( usingCurvilinearGrid ) goptCount = goptCount+1
1005 IF ( usingCylindricalGrid ) goptCount = goptCount+1
1006 IF ( goptCount .GT. 1 ) THEN
1007 WRITE(msgBuf,'(A)')
1008 & 'S/R INI_PARMS: More than one coordinate system requested'
1009 CALL PRINT_ERROR( msgBuf , myThid)
1010 STOP 'ABNORMAL END: S/R INI_PARMS'
1011 ENDIF
1012 IF ( goptCount .LT. 1 ) THEN
1013 C- No horizontal grid is specified => use Cartesian grid as default:
1014 WRITE(msgBuf,'(A)')
1015 & 'S/R INI_PARMS: No horizontal grid requested'
1016 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
1017 & SQUEEZE_RIGHT , myThid)
1018 WRITE(msgBuf,'(A)')
1019 & 'S/R INI_PARMS: => Use Cartesian Grid as default'
1020 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
1021 & SQUEEZE_RIGHT , myThid)
1022 usingCartesianGrid = .TRUE.
1023 ENDIF
1024 C-- Make metric term settings consistent with underlying grid.
1025 IF ( usingCartesianGrid ) THEN
1026 usingSphericalPolarMterms = .FALSE.
1027 metricTerms = .FALSE.
1028 useNHMTerms = .FALSE.
1029 mTFacMom = 0.
1030 useBetaPlaneF = .TRUE.
1031 ENDIF
1032 C-- Make metric term settings consistent with underlying grid.
1033 IF ( usingCylindricalGrid) THEN
1034 usingSphericalPolarMterms = .FALSE.
1035 metricTerms = .FALSE.
1036 useNHMTerms = .FALSE.
1037 mTFacMom = 1.
1038 useBetaPlaneF = .TRUE.
1039 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; Cylinder OK'
1040 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1041 & SQUEEZE_RIGHT , 1)
1042
1043 ENDIF
1044
1045 IF ( usingSphericalPolarGrid ) THEN
1046 useSphereF = .TRUE.
1047 usingSphericalPolarMterms = metricTerms
1048 ENDIF
1049 IF ( usingCurvilinearGrid ) THEN
1050 useSphereF = .TRUE.
1051 metricTerms = .FALSE.
1052 useNHMTerms = .FALSE.
1053 ENDIF
1054 C-- Set default for latitude-band where relaxation to climatology applies
1055 IF ( latBandClimRelax .EQ. UNSET_RL) THEN
1056 IF ( usingCartesianGrid ) latBandClimRelax = delY(1)*Ny*Ny
1057 IF ( usingSphericalPolarGrid ) latBandClimRelax= 180. _d 0
1058 IF ( usingCurvilinearGrid ) latBandClimRelax= 180. _d 0
1059 ENDIF
1060 C-- set cell Center depth and put Interface at the middle between 2 C
1061 setCenterDr = .FALSE.
1062 IF (delRc(1).NE.UNSET_RL) setCenterDr=.TRUE.
1063 DO K=1,Nr+1
1064 IF (delRc(K).EQ.UNSET_RL) setCenterDr = .FALSE.
1065 ENDDO
1066 C-- p, z, r coord parameters
1067 DO K = 1, Nr
1068 IF ( delZ(K) .NE. UNSET_RL ) zCoordInputData = .TRUE.
1069 IF ( delP(K) .NE. UNSET_RL ) pCoordInputData = .TRUE.
1070 IF ( delR(K) .NE. UNSET_RL ) rCoordInputData = .TRUE.
1071 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delZ(K)
1072 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delP(K)
1073 IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delRDefault(K)
1074 IF (.NOT.setCenterDr .AND. delR(K).EQ.delRDefault(K) ) THEN
1075 WRITE(msgBuf,'(A,I4)')
1076 & 'S/R INI_PARMS: No value for delZ/delP/delR at K = ',K
1077 CALL PRINT_ERROR( msgBuf , 1)
1078 STOP 'ABNORMAL END: S/R INI_PARMS'
1079 ELSEIF ( setCenterDr .AND. delR(K).NE.delRDefault(K) ) THEN
1080 WRITE(msgBuf,'(2A,I4)') 'S/R INI_PARMS:',
1081 & ' Cannot specify both delRc and delZ/delP/delR at K=', K
1082 CALL PRINT_ERROR( msgBuf , 1)
1083 STOP 'ABNORMAL END: S/R INI_PARMS'
1084 ENDIF
1085 ENDDO
1086 C Check for multiple coordinate systems
1087 CoordsSet = 0
1088 IF ( zCoordInputData ) coordsSet = coordsSet + 1
1089 IF ( pCoordInputData ) coordsSet = coordsSet + 1
1090 IF ( rCoordInputData ) coordsSet = coordsSet + 1
1091 IF ( coordsSet .GT. 1 ) THEN
1092 WRITE(msgBuf,'(A)')
1093 & 'S/R INI_PARMS: Cannot mix z, p and r in the input data.'
1094 CALL PRINT_ERROR( msgBuf , myThid)
1095 STOP 'ABNORMAL END: S/R INI_PARMS'
1096 ENDIF
1097
1098 C-- When using the dynamical pressure in EOS (with Z-coord.),
1099 C needs to activate specific part of the code (restart & exchange)
1100 c useDynP_inEos_Zc = .FALSE.
1101 useDynP_inEos_Zc = ( fluidIsWater .AND. usingZCoords
1102 & .AND. ( eosType .EQ. 'JMD95P' .OR.
1103 & eosType .EQ. 'UNESCO' .OR.
1104 & eosType .EQ. 'MDJWF' ) )
1105
1106 C-- Input files
1107 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; starts to read PARM05'
1108 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1109 & SQUEEZE_RIGHT , 1)
1110 READ(UNIT=iUnit,NML=PARM05) !,IOSTAT=errIO)
1111 IF ( errIO .LT. 0 ) THEN
1112 WRITE(msgBuf,'(A)')
1113 & 'Error reading numerical model '
1114 CALL PRINT_ERROR( msgBuf , 1)
1115 WRITE(msgBuf,'(A)')
1116 & 'parameter file "data"'
1117 CALL PRINT_ERROR( msgBuf , 1)
1118 WRITE(msgBuf,'(A)')
1119 & 'Problem in namelist PARM05'
1120 CALL PRINT_ERROR( msgBuf , 1)
1121 CALL MODELDATA_EXAMPLE( myThid )
1122 STOP 'ABNORMAL END: S/R INI_PARMS'
1123 ELSE
1124 WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; read PARM05 : OK'
1125 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1126 & SQUEEZE_RIGHT , 1)
1127 ENDIF
1128
1129 C-- Set factors required for mixing pressure and meters as vertical coordinate.
1130 C rkFac is a "sign" parameter which is used where the orientation of the vertical
1131 C coordinate (pressure or meters) relative to the vertical index (K) is important.
1132 C rkFac = 1 applies when K and the coordinate are in the opposite sense.
1133 C rkFac = -1 applies when K and the coordinate are in the same sense.
1134 C horiVertRatio is a parameter that maps horizontal units to vertical units.
1135 C It is used in certain special cases where lateral and vertical terms are
1136 C being combined and a single frame of reference is needed.
1137 IF ( zCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
1138 rkFac = 1.D0
1139 horiVertRatio = 1.D0
1140 ENDIF
1141 IF ( pCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
1142 C- jmc: any time P-coordinate is used (ocean,atmos), it requires rkFac=1
1143 c rkFac = -1.D0
1144 horiVertRatio = Gravity * rhoConst
1145 ENDIF
1146 IF ( rCoordInputData .AND. rkFac .EQ. UNSET_RS ) THEN
1147 rkFac = 1.D0
1148 horiVertRatio = 1.D0
1149 ENDIF
1150 gravitySign = -1. _d 0
1151 IF ( usingPCoords ) THEN
1152 gravitySign = 1. _d 0
1153 horiVertRatio = Gravity * rhoConst
1154 ENDIF
1155 convertEmP2rUnit = rhoConstFresh*recip_rhoConst*horiVertRatio
1156 IF ( rkFac .EQ. UNSET_RS ) rkFac=rkFacDefault
1157 recip_rkFac = 1.D0 / rkFac
1158 recip_horiVertRatio = 1./horiVertRatio
1159
1160 c-- gradually replacing debugMode by debugLevel
1161 IF ( debugMode ) debugLevel = debLevB
1162
1163 c-- flag for approximate adjoint
1164 IF ( inAdExact ) THEN
1165 inAdTrue = .FALSE.
1166 inAdFALSE = .FALSE.
1167 ELSE
1168 inAdTrue = .TRUE.
1169 inAdFALSE = .FALSE.
1170 ENDIF
1171 C
1172 CLOSE(iUnit)
1173
1174 C-- Check whether any retired parameters were found.
1175 C-- Stop if they were
1176 IF ( nRetired .GT. 0 ) THEN
1177 WRITE(msgBuf,'(A)')
1178 & 'Error reading '
1179 CALL PRINT_ERROR( msgBuf , 1)
1180 WRITE(msgBuf,'(A)')
1181 & 'parameter file "data"'
1182 CALL PRINT_ERROR( msgBuf , 1)
1183 WRITE(msgBuf,'(A)')
1184 & 'some out of date parameters were found in the namelist'
1185 CALL PRINT_ERROR( msgBuf , 1)
1186 STOP 'ABNORMAL END: S/R INI_PARMS'
1187 ENDIF
1188
1189 _END_MASTER(myThid)
1190
1191 C-- Everyone else must wait for the parameters to be loaded
1192 _BARRIER
1193 C
1194 RETURN
1195 END

  ViewVC Help
Powered by ViewVC 1.1.22