--- MITgcm_contrib/osse/codemod/ini_parms.F 2004/06/22 19:44:40 1.1 +++ MITgcm_contrib/osse/codemod/ini_parms.F 2004/06/22 20:37:41 1.2 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/osse/codemod/ini_parms.F,v 1.1 2004/06/22 19:44:40 afe Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/osse/codemod/ini_parms.F,v 1.2 2004/06/22 20:37:41 afe Exp $ C $Name: $ #include "CPP_OPTIONS.h" @@ -88,11 +88,36 @@ C :: used. coordsSet counts how many vertical coordinate systems have been C used to specify variables. coordsSet > 1 is an error. C + LOGICAL zCoordInputData LOGICAL pCoordInputData LOGICAL rCoordInputData INTEGER coordsSet +C Variables which have vertical coordinate system dependency. +C delZ :: Vertical grid spacing ( m ). +C delP :: Vertical grid spacing ( Pa ). +C viscAz :: Eddy viscosity coeff. for mixing of +C momentum vertically ( m^2/s ) +C viscAp :: Eddy viscosity coeff. for mixing of +C momentum vertically ( Pa^2/s ) +C diffKzT :: Laplacian diffusion coeff. for mixing of +C heat vertically ( m^2/s ) +C diffKpT :: Laplacian diffusion coeff. for mixing of +C heat vertically ( Pa^2/s ) +C diffKzS :: Laplacian diffusion coeff. for mixing of +C salt vertically ( m^2/s ) +C diffKpS :: Laplacian diffusion coeff. for mixing of +C salt vertically ( Pa^2/s ) + _RL delZ(Nr) + _RL delP(Nr) + _RL viscAz + _RL viscAp + _RL diffKzT + _RL diffKpT + _RL diffKzS + _RL diffKpS + C Retired main data file parameters. Kept here to trap use of old data files. C zonal_filt_lat - Moved to package "zonal_filt" C nRetired :: Counter used to trap gracefully namelists containing "retired" @@ -105,12 +130,15 @@ C-- Continuous equation parameters NAMELIST /PARM01/ & gravitySign, - & gravity, gBaro, rhonil, tAlpha, sBeta, f0, beta, omega, - & viscAh, viscAz, viscA4, cosPower, viscAstrain, viscAtension, + & gravity, gBaro, rhonil, tAlpha, sBeta, + & f0, beta, omega, rotationPeriod, + & viscAh, viscAhMax, viscAhGrid, viscC2leith, + & viscA4, viscA4Max, viscA4Grid, viscC4leith, + & viscAz, cosPower, viscAstrain, viscAtension, & diffKhT, diffKzT, diffK4T, & diffKhS, diffKzS, diffK4S, & tRef, sRef, eosType, integr_GeoPot, selectFindRoSurf, - & atm_Cp, atm_Rd, + & atm_Cp, atm_Rd, atm_Rq, & no_slip_sides,no_slip_bottom, & momViscosity, momAdvection, momForcing, useCoriolis, & momPressureForcing, metricTerms, vectorInvariantMomentum, @@ -119,21 +147,26 @@ & implicSurfPress, implicDiv2DFlow, & implicitFreeSurface, rigidLid, freeSurfFac, hFacMin, hFacMinDz, & exactConserv,uniformLin_PhiSurf,nonlinFreeSurf,hFacInf,hFacSup, + & select_rStar, & staggerTimeStep, & tempStepping, saltStepping, momStepping, tr1Stepping, & implicitDiffusion, implicitViscosity, + & tempImplVertAdv, saltImplVertAdv, momImplVertAdv, & viscAr, diffKrT, diffKrS, hFacMinDr, & viscAp, diffKpT, diffKpS, hFacMinDp, + & diffKrBL79surf, diffKrBL79deep, diffKrBL79scl, diffKrBL79Ho, & rhoConst, rhoConstFresh, buoyancyRelation, HeatCapacity_Cp, & writeBinaryPrec, readBinaryPrec, writeStatePrec, - & nonHydrostatic, quasiHydrostatic, globalFiles, - & allowFreezing, ivdc_kappa, + & nonHydrostatic, quasiHydrostatic, globalFiles, useSingleCpuIO, + & allowFreezing, useOldFreezing, ivdc_kappa, & bottomDragLinear,bottomDragQuadratic, - & usePickupBeforeC35, debugMode, + & usePickupBeforeC35, debugMode, debugLevel, & readPickupWithTracer, writePickupWithTracer, & tempAdvScheme, saltAdvScheme, tracerAdvScheme, & multiDimAdvection, useEnergyConservingCoriolis, - & useJamartWetPoints, useNHMTerms, + & useCDscheme, useJamartWetPoints, useJamartMomAdv, useNHMTerms, + & SadournyCoriolis, upwindVorticity, highOrderVorticity, + & useAbsVorticity, & useRealFreshWaterFlux, convertFW2Salt, & temp_EvPrRn, salt_EvPrRn, trac_EvPrRn, & zonal_filt_lat @@ -146,20 +179,21 @@ C-- Time stepping parammeters NAMELIST /PARM03/ - & nIter0, nTimeSteps, nEndIter, + & nIter0, nTimeSteps, nEndIter, pickupSuff, & deltaT, deltaTmom, deltaTtracer, deltaTfreesurf, & forcing_In_AB, abEps, tauCD, rCD, & startTime, endTime, chkPtFreq, - & dumpFreq, taveFreq, tave_lastIter, deltaTClock, diagFreq, - & monitorFreq, pChkPtFreq, cAdjFreq, - & tauThetaClimRelax, tauSaltClimRelax, tauTr1ClimRelax, + & dumpFreq, adjDumpFreq, taveFreq, tave_lastIter, deltaTClock, + & diagFreq, monitorFreq, pChkPtFreq, cAdjFreq, + & tauThetaClimRelax, tauSaltClimRelax, latBandClimRelax, + & tauTr1ClimRelax, & periodicExternalForcing, externForcingPeriod, externForcingCycle C-- Gridding parameters NAMELIST /PARM04/ & usingCartesianGrid, dxSpacing, dySpacing, delX, delY, delZ, & usingSphericalPolarGrid, phiMin, thetaMin, rSphere, - & usingCurvilinearGrid, bUseCylindricalGrid, + & usingCurvilinearGrid,bUseCylindricalGrid, & delP, delR, rkFac, Ro_SeaLevel, groundAtK1, delRc, & delXfile, delYfile @@ -170,7 +204,8 @@ & thetaClimFile, saltClimFile, & surfQfile, EmPmRfile, surfQswfile, & uVelInitFile, vVelInitFile, pSurfInitFile, - & dQdTFile, ploadFile, tCyl + & dQdTFile, ploadFile,tCyl + & mdsioLocalDir C _BEGIN_MASTER(myThid) @@ -192,6 +227,7 @@ C-- Iniialise retired parameters to unlikely value nRetired = 0 zonal_filt_lat = UNSET_RL + gravitySign = UNSET_RL C-- Open the parameter file OPEN(UNIT=scrUnit1,STATUS='SCRATCH') @@ -215,9 +251,11 @@ DO WHILE ( .TRUE. ) READ(modelDataUnit,FMT='(A)',END=1001) RECORD IL = MAX(ILNBLNK(RECORD),1) - IF ( RECORD(1:1) .NE. commentCharacter ) - & WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL) - WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL) + IF ( RECORD(1:1) .NE. commentCharacter ) THEN + CALL NML_SET_TERMINATOR( RECORD ) + WRITE(UNIT=scrUnit1,FMT='(A)') RECORD(:IL) + ENDIF + WRITE(UNIT=scrUnit2,FMT='(A)') RECORD(:IL) ENDDO 1001 CONTINUE CLOSE(modelDataUnit) @@ -266,9 +304,11 @@ diffKrS = UNSET_RL gBaro = UNSET_RL rhoConst = UNSET_RL + omega = UNSET_RL hFacMinDr = UNSET_RL hFacMinDz = UNSET_RL hFacMinDp = UNSET_RL + rhoConstFresh = UNSET_RL convertFW2Salt = UNSET_RL tAlpha = UNSET_RL sBeta = UNSET_RL @@ -298,11 +338,32 @@ IF ( rigidLid ) freeSurfFac = 0.D0 IF ( gBaro .EQ. UNSET_RL ) gBaro=gravity IF ( rhoConst .EQ. UNSET_RL ) rhoConst=rhoNil + IF ( rhoConstFresh .EQ. UNSET_RL ) rhoConstFresh=rhoConst + IF ( omega .EQ. UNSET_RL ) THEN + omega = 0. _d 0 + IF ( rotationPeriod .NE. 0. _d 0 ) + & omega = 2.D0 * PI / rotationPeriod + ELSEIF ( omega .EQ. 0. _d 0 ) THEN + rotationPeriod = 0. _d 0 + ELSE + rotationPeriod = 2.D0 * PI / omega + ENDIF IF (atm_Rd .EQ. UNSET_RL) THEN atm_Rd = atm_Cp * atm_kappa ELSE atm_kappa = atm_Rd / atm_Cp ENDIF +C-- On/Off flags for each terms of the momentum equation + nonHydrostatic = momStepping .AND. nonHydrostatic + quasiHydrostatic = momStepping .AND. quasiHydrostatic + momAdvection = momStepping .AND. momAdvection + momViscosity = momStepping .AND. momViscosity + momForcing = momStepping .AND. momForcing + useCoriolis = momStepping .AND. useCoriolis + useCDscheme = momStepping .AND. useCDscheme + momPressureForcing= momStepping .AND. momPressureForcing + momImplVertAdv = momAdvection .AND. momImplVertAdv + implicitViscosity= momViscosity .AND. implicitViscosity C-- Momentum viscosity on/off flag. IF ( momViscosity ) THEN vfFacMom = 1.D0 @@ -351,6 +412,8 @@ tempForcing = tempStepping .AND. tempForcing saltAdvection = saltStepping .AND. saltAdvection saltForcing = saltStepping .AND. saltForcing + tempImplVertAdv = tempAdvection .AND. tempImplVertAdv + saltImplVertAdv = saltAdvection .AND. saltImplVertAdv C-- z,p,r coord input switching. IF ( viscAz .NE. UNSET_RL ) zCoordInputData = .TRUE. IF ( viscAp .NE. UNSET_RL ) pCoordInputData = .TRUE. @@ -435,12 +498,13 @@ ELSE recip_gravity = 1.D0 / gravity ENDIF +C This flags are now passed to RW and MDSIO packages in ini_model_io.F C Set globalFiles flag for READ_WRITE_FLD package - CALL SET_WRITE_GLOBAL_FLD( globalFiles ) +c CALL SET_WRITE_GLOBAL_FLD( globalFiles ) C Set globalFiles flag for READ_WRITE_REC package - CALL SET_WRITE_GLOBAL_REC( globalFiles ) +c CALL SET_WRITE_GLOBAL_REC( globalFiles ) C Set globalFiles flag for READ_WRITE_REC package - CALL SET_WRITE_GLOBAL_PICKUP( globalFiles ) +c CALL SET_WRITE_GLOBAL_PICKUP( globalFiles ) C Check for retired parameters still being used nRetired = 0 @@ -455,6 +519,13 @@ & ' now read from file "data.zonfilt".' CALL PRINT_ERROR( msgBuf , myThid) ENDIF + IF ( gravitySign .NE. UNSET_RL ) THEN + nRetired = nRetired+1 + WRITE(msgBuf,'(A,A)') + & 'S/R INI_PARMS: "gravitySign" is set according to vertical ', + & ' coordinate and is no longer allowed in file "data".' + CALL PRINT_ERROR( msgBuf , myThid) + ENDIF C-- Elliptic solver parameters READ(UNIT=iUnit,NML=PARM02) !,IOSTAT=errIO) @@ -481,6 +552,7 @@ C-- Time stepping parameters rCD = -1.D0 + latBandClimRelax = UNSET_RL READ(UNIT=iUnit,NML=PARM03) !,IOSTAT=errIO) IF ( errIO .LT. 0 ) THEN WRITE(msgBuf,'(A)') @@ -553,12 +625,10 @@ CALL PRINT_ERROR( msgBuf , myThid) STOP 'ABNORMAL END: S/R INI_PARMS' ENDIF -C o CD coupling - IF ( tauCD .EQ. 0.D0 ) THEN - tauCD = deltaTmom - ENDIF - IF ( rCD .LT. 0. ) THEN - rCD = 1. - deltaTMom/tauCD + IF (useCDscheme) THEN +C o CD coupling (CD scheme): + IF ( tauCD .EQ. 0.D0 ) tauCD = deltaTmom + IF ( rCD .LT. 0. ) rCD = 1. _d 0 - deltaTMom/tauCD ENDIF C o Temperature climatology relaxation time scale IF ( tauThetaClimRelax .EQ. 0.D0 ) THEN @@ -782,10 +852,16 @@ STOP 'ABNORMAL END: S/R INI_PARMS' ENDIF IF ( goptCount .LT. 1 ) THEN +C- No horizontal grid is specified => use Cartesian grid as default: WRITE(msgBuf,'(A)') - & 'S/R INI_PARMS: No coordinate system requested' - CALL PRINT_ERROR( msgBuf , myThid) - STOP 'ABNORMAL END: S/R INI_PARMS' + & 'S/R INI_PARMS: No horizontal grid requested' + CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, + & SQUEEZE_RIGHT , myThid) + WRITE(msgBuf,'(A)') + & 'S/R INI_PARMS: => Use Cartesian Grid as default' + CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, + & SQUEEZE_RIGHT , myThid) + usingCartesianGrid = .TRUE. ENDIF C-- Make metric term settings consistent with underlying grid. IF ( usingCartesianGrid ) THEN @@ -804,9 +880,10 @@ useBetaPlaneF = .TRUE. WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; Cylinder OK' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, - & SQUEEZE_RIGHT , 1) - + & SQUEEZE_RIGHT , 1) + ENDIF + IF ( usingSphericalPolarGrid ) THEN useConstantF = .FALSE. useBetaPlaneF = .FALSE. @@ -818,6 +895,12 @@ metricTerms = .FALSE. useNHMTerms = .FALSE. ENDIF +C-- Set default for latitude-band where relaxation to climatology applies + IF ( latBandClimRelax .EQ. UNSET_RL) THEN + IF ( usingCartesianGrid ) latBandClimRelax = delY(1)*Ny*Ny + IF ( usingSphericalPolarGrid ) latBandClimRelax= 180. _d 0 + IF ( usingCurvilinearGrid ) latBandClimRelax= 180. _d 0 + ENDIF C-- set cell Center depth and put Interface at the middle between 2 C setCenterDr = .FALSE. IF (delRc(1).NE.UNSET_RL) setCenterDr=.TRUE. @@ -833,9 +916,8 @@ IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delP(K) IF ( delR(K) .EQ. UNSET_RL ) delR(K) = delRDefault(K) IF (.NOT.setCenterDr .AND. delR(K).EQ.delRDefault(K) ) THEN - WRITE(msgBuf,'(A,I4,I4.2)') - & 'S/R INI_PARMS: No value for delZ/delP/delR at K = ',K, - & dXspacing + WRITE(msgBuf,'(A,I4)') + & 'S/R INI_PARMS: No value for delZ/delP/delR at K = ',K CALL PRINT_ERROR( msgBuf , 1) STOP 'ABNORMAL END: S/R INI_PARMS' ELSEIF ( setCenterDr .AND. delR(K).NE.delRDefault(K) ) THEN @@ -857,6 +939,14 @@ STOP 'ABNORMAL END: S/R INI_PARMS' ENDIF +C-- When using the dynamical pressure in EOS (with Z-coord.), +C needs to activate specific part of the code (restart & exchange) +c useDynP_inEos_Zc = .FALSE. + useDynP_inEos_Zc = ( buoyancyRelation .EQ. 'OCEANIC' + & .AND. ( eosType .EQ. 'JMD95P' .OR. + & eosType .EQ. 'UNESCO' .OR. + & eosType .EQ. 'MDJWF' ) ) + C-- Input files READ(UNIT=iUnit,NML=PARM05) !,IOSTAT=errIO) IF ( errIO .LT. 0 ) THEN @@ -898,19 +988,24 @@ rkFac = 1.D0 horiVertRatio = 1.D0 ENDIF - convertEmP2rUnit = 1. _d 0 - IF (buoyancyRelation.EQ.'ATMOSPHERIC') - & horiVertRatio = Gravity * rhoConst + gravitySign = -1. _d 0 + IF (buoyancyRelation.EQ.'ATMOSPHERIC') THEN + gravitySign = 1. _d 0 + horiVertRatio = Gravity * rhoConst + ENDIF IF (buoyancyRelation.EQ.'OCEANICP') THEN + gravitySign = 1. _d 0 horiVertRatio = Gravity * rhoConst - convertEmP2rUnit = Gravity * rhoConstFresh ENDIF + convertEmP2rUnit = rhoConstFresh*recip_rhoConst*horiVertRatio IF ( rkFac .EQ. UNSET_RS ) rkFac=rkFacDefault recip_rkFac = 1.D0 / rkFac recip_horiVertRatio = 1./horiVertRatio IF ( zCoordInputData ) usingZCoords = .TRUE. IF ( pCoordInputData ) usingPCoords = .TRUE. +c-- gradually replacing debugMode by debugLevel + IF ( debugMode ) debugLevel = debLevB C CLOSE(iUnit)