model/src/set_parms.F

4e66ab0b67 Oliv*0001 #include "PACKAGES_CONFIG.h"
809bdccbfc Jean*0002 #include "CPP_OPTIONS.h"
0301aae365 Jean*0003 #ifdef ALLOW_MOM_COMMON
                0004 # include "MOM_COMMON_OPTIONS.h"
                0005 #endif
809bdccbfc Jean*0006
                0007 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0008 CBOP
                0009 C     !ROUTINE: SET_PARMS
                0010 C     !INTERFACE:
                0011       SUBROUTINE SET_PARMS( myThid )
                0012
                0013 C     !DESCRIPTION:
                0014 C     Set model "parameters" that might depend on the use of some pkgs;
                0015 C     called from INITIALISE_FIXED, after INI_PARMS & PACKAGES_READPARAMS
                0016 C     NOTES: After leaving this S/R, parameters will not change anymore.
                0017
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
c7ad17745a Jean*0023 #include "EOS.h"
0301aae365 Jean*0024 #ifdef ALLOW_MOM_COMMON
                0025 # include "MOM_VISC.h"
                0026 #endif
809bdccbfc Jean*0027
                0028 C     !INPUT/OUTPUT PARAMETERS:
3fcd8a21e5 Jean*0029 C     myThid :: My Thread Id number
809bdccbfc Jean*0030       INTEGER myThid
                0031
1cfd78e5a7 Jean*0032 C     !FUNCTIONS:
809bdccbfc Jean*0033 c     INTEGER  ILNBLNK
                0034 c     EXTERNAL ILNBLNK
1cfd78e5a7 Jean*0035
                0036 C     !LOCAL VARIABLES:
                0037       CHARACTER*(MAX_LEN_MBUF) msgBuf
aa6b2555c8 Jean*0038       INTEGER errCount, k
1cfd78e5a7 Jean*0039       _RL tmpVar
809bdccbfc Jean*0040 CEOP
                0041
ac2fc522e8 Jean*0042 C--   Set (or reset) On/Off flags :
809bdccbfc Jean*0043
ac2fc522e8 Jean*0044 C-    For off-line calculation, switch off Momentum and Active-tracers (=T,S):
                0045 #ifdef ALLOW_OFFLINE
809bdccbfc Jean*0046       IF ( useOffLine ) THEN
ac2fc522e8 Jean*0047         CALL OFFLINE_RESET_PARMS( myThid )
809bdccbfc Jean*0048       ENDIF
ac2fc522e8 Jean*0049 #endif /* ALLOW_OFFLINE */
809bdccbfc Jean*0050
ac2fc522e8 Jean*0051       _BEGIN_MASTER(myThid)
aa6b2555c8 Jean*0052       errCount = 0
809bdccbfc Jean*0053
427e24e121 Jean*0054 C--   Make metric term & Coriolis settings consistent with underlying grid.
                0055       IF ( usingCartesianGrid ) THEN
                0056         selectMetricTerms = 0
                0057         useNHMTerms   = .FALSE.
                0058       ENDIF
                0059       IF ( usingCylindricalGrid ) THEN
                0060         useNHMTerms   = .FALSE.
                0061         WRITE(msgBuf,'(A)') 'S/R INI_PARMS ; Cylinder OK'
                0062         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0063      &                     SQUEEZE_RIGHT, myThid )
                0064       ENDIF
                0065       IF ( usingCurvilinearGrid ) THEN
31fb0e0e6d Jean*0066 C-    Horizontal metric terms not implemented for Curvilinear-Grid
427e24e121 Jean*0067         selectMetricTerms = 0
                0068       ENDIF
31fb0e0e6d Jean*0069       IF ( usingCylindricalGrid ) THEN
                0070 C-    No current alternative metric terms formulation for Cylindrical-Grid
                0071         selectMetricTerms = MIN( selectMetricTerms, 1 )
                0072       ENDIF
427e24e121 Jean*0073       IF ( selectCoriMap.EQ.-1 ) THEN
                0074         IF ( usingCartesianGrid.OR.usingCylindricalGrid ) THEN
                0075 C       default is to use Beta-Plane Coriolis
                0076           selectCoriMap = 1
                0077         ELSE
                0078 C       default for other grids is to use Spherical Coriolis
                0079           selectCoriMap = 2
                0080         ENDIF
                0081       ENDIF
                0082       IF ( .NOT.(nonHydrostatic.OR.quasiHydrostatic) )
                0083      &                          select3dCoriScheme = 0
                0084       IF ( (selectCoriMap.EQ.0 .OR.selectCoriMap.EQ.1)
                0085      &     .AND. fPrime.EQ.0. ) select3dCoriScheme = 0
                0086
809bdccbfc Jean*0087 C--   On/Off flags for each terms of the momentum equation
                0088       nonHydrostatic   = momStepping .AND. nonHydrostatic
                0089       quasiHydrostatic = momStepping .AND. quasiHydrostatic
59b35dd864 Jean*0090       momAdvection     = momStepping .AND. momAdvection
                0091       momViscosity     = momStepping .AND. momViscosity
                0092       momForcing       = momStepping .AND. momForcing
                0093       momTidalForcing  = momForcing  .AND. momTidalForcing
                0094       useCoriolis      = momStepping .AND. useCoriolis
427e24e121 Jean*0095       IF ( .NOT.useCoriolis ) select3dCoriScheme = 0
59b35dd864 Jean*0096       useCDscheme      = momStepping .AND. useCDscheme
809bdccbfc Jean*0097       momPressureForcing= momStepping .AND. momPressureForcing
                0098       implicitIntGravWave=momPressureForcing .AND. implicitIntGravWave
                0099       momImplVertAdv   = momAdvection .AND. momImplVertAdv
31fb0e0e6d Jean*0100       useNHMTerms      = momAdvection .AND. useNHMTerms
                0101       IF ( .NOT.momAdvection ) selectMetricTerms = 0
809bdccbfc Jean*0102       implicitViscosity= momViscosity .AND. implicitViscosity
dc3adfb09b Jean*0103       useSmag3D        = momViscosity .AND. useSmag3D
809bdccbfc Jean*0104       use3Dsolver      = nonHydrostatic.OR. implicitIntGravWave
4947dcf697 Jean*0105       calc_wVelocity   = momStepping .OR. exactConserv
3fcd8a21e5 Jean*0106
0301aae365 Jean*0107 #ifndef ALLOW_3D_VISCAH
aa6b2555c8 Jean*0108       IF ( viscAhDfile.NE.' ' .OR. viscAhZfile.NE.' ' ) THEN
1a6ecc2dd5 Jean*0109         WRITE(msgBuf,'(2A)') 'SET_PARMS: ',
0301aae365 Jean*0110      &      'viscAhDfile and viscAhZfile cannot be used with'
                0111         CALL PRINT_ERROR( msgBuf, myThid )
1a6ecc2dd5 Jean*0112         WRITE(msgBuf,'(2A)') 'SET_PARMS: ',
0301aae365 Jean*0113      &      '"#undef ALLOW_3D_VISCAH" in MOM_COMMON_OPTIONS.h'
                0114         CALL PRINT_ERROR( msgBuf, myThid )
aa6b2555c8 Jean*0115         errCount = errCount + 1
                0116       ENDIF
0301aae365 Jean*0117 #endif
                0118 #ifndef ALLOW_3D_VISCA4
aa6b2555c8 Jean*0119       IF ( viscA4Dfile.NE.' ' .OR. viscA4Zfile.NE.' ' ) THEN
1a6ecc2dd5 Jean*0120         WRITE(msgBuf,'(2A)') 'SET_PARMS: ',
0301aae365 Jean*0121      &      'viscA4Dfile and viscA4Zfile cannot be used with'
                0122         CALL PRINT_ERROR( msgBuf, myThid )
1a6ecc2dd5 Jean*0123         WRITE(msgBuf,'(2A)') 'SET_PARMS: ',
0301aae365 Jean*0124      &      '"#undef ALLOW_3D_VISCA4" in MOM_COMMON_OPTIONS.h'
                0125         CALL PRINT_ERROR( msgBuf, myThid )
aa6b2555c8 Jean*0126         errCount = errCount + 1
                0127       ENDIF
0301aae365 Jean*0128 #endif
                0129
                0130 #ifdef ALLOW_MOM_COMMON
                0131 C-    On/Off flags for viscosity coefficients
                0132       useVariableVisc   =
                0133      &      viscAhGrid  .NE.zeroRL .OR. viscA4Grid  .NE.zeroRL
                0134      &  .OR. viscC2smag .NE.zeroRL .OR. viscC4smag  .NE.zeroRL
                0135      &  .OR. viscC2leith.NE.zeroRL .OR. viscC2leithD.NE.zeroRL
f59d76b0dd Ed D*0136      &  .OR. viscC2LeithQG.NE.zeroRL
0301aae365 Jean*0137      &  .OR. viscC4leith.NE.zeroRL .OR. viscC4leithD.NE.zeroRL
                0138      &  .OR. viscAhDfile.NE.' '    .OR.  viscAhZfile.NE.' '
                0139      &  .OR. viscA4Dfile.NE.' '    .OR.  viscA4Zfile.NE.' '
                0140
                0141       useHarmonicVisc   = viscAh .NE.zeroRL
                0142      &  .OR. viscAhD    .NE.zeroRL .OR. viscAhZ     .NE.zeroRL
                0143      &  .OR. viscAhGrid .NE.zeroRL .OR. viscC2smag  .NE.zeroRL
                0144      &  .OR. viscC2leith.NE.zeroRL .OR. viscC2leithD.NE.zeroRL
f59d76b0dd Ed D*0145      &  .OR. viscC2LeithQG.NE.zeroRL
0301aae365 Jean*0146      &  .OR. viscAhDfile.NE. ' '   .OR. viscAhZfile .NE. ' '
                0147
                0148       useBiharmonicVisc = viscA4.NE.zeroRL
                0149      &  .OR. viscA4D    .NE.zeroRL .OR. viscA4Z     .NE.zeroRL
                0150      &  .OR. viscA4Grid .NE.zeroRL .OR. viscC4smag  .NE.zeroRL
                0151      &  .OR. viscC4leith.NE.zeroRL .OR. viscC4leithD.NE.zeroRL
                0152      &  .OR. viscA4Dfile.NE. ' '   .OR. viscA4Zfile .NE. ' '
                0153
                0154       useVariableVisc   = momViscosity .AND. useVariableVisc
                0155       useHarmonicVisc   = momViscosity .AND. useHarmonicVisc
                0156       useBiharmonicVisc = momViscosity .AND. useBiharmonicVisc
                0157 #endif /* ALLOW_MOM_COMMON */
ab47de63dc Mart*0158       IF ( ( bottomDragQuadratic.EQ.0. .AND. zRoughBot.EQ.0. )
                0159      &     .OR. .NOT.momViscosity ) selectBotDragQuadr = -1
0301aae365 Jean*0160
3fcd8a21e5 Jean*0161 C--   Free-surface & pressure method
2ee74bea59 Jean*0162       uniformFreeSurfLev = usingZCoords
283763fbbb Jean*0163 C- Note: comment line below to revert to full-cell hydrostatic-pressure
                0164 C        calculation in surface grid-cell below ice-shelf
                0165       uniformFreeSurfLev = usingZCoords .AND. .NOT.useShelfIce
8233d0ceb9 Jean*0166      &                                  .AND. topoFile.EQ.' '
3fcd8a21e5 Jean*0167       IF ( selectNHfreeSurf.NE.0 .AND.
                0168      &      ( .NOT.nonHydrostatic .OR. usingPCoords
                0169      &        .OR. .NOT.exactConserv
                0170      &      ) ) THEN
                0171         WRITE(msgBuf,'(2A)') '** WARNING ** SET_PARMS: ',
                0172      &                       'reset selectNHfreeSurf to zero'
                0173         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0174      &                      SQUEEZE_RIGHT, myThid )
                0175         selectNHfreeSurf = 0
                0176       ENDIF
02d90fb24c Jean*0177 #ifdef ALLOW_AUTODIFF
627a21a418 Jean*0178       doResetHFactors = .TRUE.
                0179 #endif
                0180 #ifndef NONLIN_FRSURF
                0181       doResetHFactors = .FALSE.
                0182 #endif
aa3dd96e14 Jean*0183
7fe6343684 Jean*0184 C--   Set default Vorticity-Term Scheme:
                0185       IF ( vectorInvariantMomentum ) THEN
                0186         IF ( selectVortScheme.EQ.UNSET_I ) THEN
                0187           selectVortScheme = 1
                0188           IF ( upwindVorticity )    selectVortScheme = 0
                0189           IF ( highOrderVorticity ) selectVortScheme = 0
                0190         ENDIF
                0191       ELSEIF ( selectVortScheme.NE.UNSET_I ) THEN
                0192         WRITE(msgBuf,'(A,A)') '** WARNING ** SET_PARMS: ',
                0193      &   'Vector-Invariant Momentum unused => ignore selectVortScheme'
3fcd8a21e5 Jean*0194         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0195      &                      SQUEEZE_RIGHT, myThid )
7fe6343684 Jean*0196       ENDIF
9293d3c672 Hajo*0197       useAbsVorticity = vectorInvariantMomentum .AND. useAbsVorticity
809bdccbfc Jean*0198 C--   Momentum viscosity on/off flag.
                0199       IF ( momViscosity        ) THEN
3fcd8a21e5 Jean*0200        vfFacMom = 1. _d 0
809bdccbfc Jean*0201       ELSE
3fcd8a21e5 Jean*0202        vfFacMom = 0. _d 0
809bdccbfc Jean*0203       ENDIF
                0204 C--   Momentum advection on/off flag.
                0205       IF ( momAdvection        ) THEN
3fcd8a21e5 Jean*0206        afFacMom = 1. _d 0
809bdccbfc Jean*0207       ELSE
3fcd8a21e5 Jean*0208        afFacMom = 0. _d 0
809bdccbfc Jean*0209       ENDIF
                0210 C--   Momentum forcing on/off flag.
                0211       IF ( momForcing ) THEN
3fcd8a21e5 Jean*0212        foFacMom = 1. _d 0
809bdccbfc Jean*0213       ELSE
3fcd8a21e5 Jean*0214        foFacMom = 0. _d 0
809bdccbfc Jean*0215       ENDIF
                0216 C--   Coriolis term on/off flag.
                0217       IF ( useCoriolis ) THEN
3fcd8a21e5 Jean*0218        cfFacMom = 1. _d 0
809bdccbfc Jean*0219       ELSE
3fcd8a21e5 Jean*0220        cfFacMom = 0. _d 0
809bdccbfc Jean*0221       ENDIF
                0222 C--   Pressure term on/off flag.
                0223       IF ( momPressureForcing ) THEN
3fcd8a21e5 Jean*0224        pfFacMom = 1. _d 0
809bdccbfc Jean*0225       ELSE
3fcd8a21e5 Jean*0226        pfFacMom = 0. _d 0
809bdccbfc Jean*0227       ENDIF
                0228 C--   Metric terms on/off flag.
427e24e121 Jean*0229       IF ( selectMetricTerms.GE.1 ) THEN
3fcd8a21e5 Jean*0230        mTFacMom = 1. _d 0
809bdccbfc Jean*0231       ELSE
3fcd8a21e5 Jean*0232        mTFacMom = 0. _d 0
809bdccbfc Jean*0233       ENDIF
                0234
                0235 C--   Advection and Forcing for Temp and salt  on/off flags
bf6138bedc Jean*0236       tempVertDiff4 = .FALSE.
                0237       saltVertDiff4 = .FALSE.
                0238       DO k=1,Nr
                0239         tempVertDiff4 = tempVertDiff4 .OR. ( diffKr4T(k).GT.0. _d 0 )
                0240         saltVertDiff4 = saltVertDiff4 .OR. ( diffKr4S(k).GT.0. _d 0 )
                0241       ENDDO
809bdccbfc Jean*0242       tempAdvection = tempStepping .AND. tempAdvection
bf6138bedc Jean*0243       tempVertDiff4 = tempStepping .AND. tempVertDiff4
809bdccbfc Jean*0244       tempForcing   = tempStepping .AND. tempForcing
                0245       saltAdvection = saltStepping .AND. saltAdvection
bf6138bedc Jean*0246       saltVertDiff4 = saltStepping .AND. saltVertDiff4
809bdccbfc Jean*0247       saltForcing   = saltStepping .AND. saltForcing
                0248       tempImplVertAdv = tempAdvection .AND. tempImplVertAdv
                0249       saltImplVertAdv = saltAdvection .AND. saltImplVertAdv
7f31549cd8 Jean*0250       doThetaClimRelax = ( tempForcing .OR.( useOffLine.AND.useKPP ) )
                0251      &             .AND. ( tauThetaClimRelax.GT.0. _d 0 )
                0252       doSaltClimRelax  = ( saltForcing .OR.( useOffLine.AND.useKPP ) )
                0253      &             .AND. ( tauSaltClimRelax .GT.0. _d 0 )
809bdccbfc Jean*0254
c7ad17745a Jean*0255 C--   Dynamically Active Tracers : set flags
                0256       tempIsActiveTr = momPressureForcing .AND. tempAdvection
                0257       saltIsActiveTr = momPressureForcing .AND. saltAdvection
aa3dd96e14 Jean*0258       IF ( eosType.EQ.'IDEALG' .AND. atm_Rq.EQ.0. ) THEN
c7ad17745a Jean*0259         saltIsActiveTr = .FALSE.
                0260       ELSEIF ( eosType.EQ.'LINEAR' ) THEN
                0261         IF ( tAlpha.EQ.0. ) tempIsActiveTr = .FALSE.
                0262         IF ( sBeta .EQ.0. ) saltIsActiveTr = .FALSE.
                0263       ENDIF
                0264
427e24e121 Jean*0265 C--   Set default for latitude-band where relaxation to climatology applies
                0266 C     note: done later (once domain size is known) if using CartesianGrid
                0267       IF ( latBandClimRelax .EQ. UNSET_RL) THEN
                0268         IF ( usingSphericalPolarGrid ) latBandClimRelax= 180. _d 0
                0269         IF ( usingCurvilinearGrid )    latBandClimRelax= 180. _d 0
                0270       ENDIF
                0271
6ef71429db Jean*0272       IF ( usingZCoords ) THEN
                0273 C--   Select which pressure to use in EOS:
                0274 C     set default according to EOS type (as it was until chkpt65t)
                0275         IF ( selectP_inEOS_Zc.EQ.UNSET_I ) THEN
                0276           IF ( eosType .EQ. 'JMD95P' .OR.  eosType .EQ. 'UNESCO'
                0277      &    .OR. eosType .EQ. 'MDJWF'  .OR.  eosType .EQ. 'TEOS10'
                0278      &       ) THEN
                0279            selectP_inEOS_Zc = 2
                0280           ELSE
                0281            selectP_inEOS_Zc = 0
                0282           ENDIF
                0283         ELSEIF ( selectP_inEOS_Zc.LT.0
                0284      &      .OR. selectP_inEOS_Zc.GT.3 ) THEN
1a6ecc2dd5 Jean*0285           WRITE(msgBuf,'(A,I9,A)') 'SET_PARMS: selectP_inEOS_Zc=',
6ef71429db Jean*0286      &                  selectP_inEOS_Zc, ' : invalid selection'
                0287           CALL PRINT_ERROR( msgBuf, myThid )
aa6b2555c8 Jean*0288           errCount = errCount + 1
6ef71429db Jean*0289         ELSEIF ( .NOT.nonHydrostatic ) THEN
                0290           selectP_inEOS_Zc = MIN( selectP_inEOS_Zc, 2 )
                0291         ENDIF
                0292         IF ( ( eosType .EQ. 'LINEAR' .OR.  eosType .EQ. 'POLY3 ' )
                0293      &      .AND. selectP_inEOS_Zc.NE.0  ) THEN
1a6ecc2dd5 Jean*0294           WRITE(msgBuf,'(A,I9,2A)') 'SET_PARMS: selectP_inEOS_Zc=',
4382c26089 Jean*0295      &     selectP_inEOS_Zc, ' : invalid with eosType=', eosType
                0296           CALL PRINT_ERROR( msgBuf, myThid )
aa6b2555c8 Jean*0297           errCount = errCount + 1
6ef71429db Jean*0298         ENDIF
                0299       ELSE
                0300         selectP_inEOS_Zc = -1
                0301       ENDIF
809bdccbfc Jean*0302 C--   When using the dynamical pressure in EOS (with Z-coord.),
                0303 C     needs to activate specific part of the code (restart & exchange)
6ef71429db Jean*0304       storePhiHyd4Phys = selectP_inEOS_Zc.GE.2
1861cae8e7 Jean*0305 C-    pkg/atm_phys uses main-model geopotential:
6ef71429db Jean*0306       storePhiHyd4Phys = storePhiHyd4Phys .OR. useAtm_Phys
809bdccbfc Jean*0307
1cfd78e5a7 Jean*0308 C--   Adjust parameters related to length of the simulation
                0309
                0310 C-    Need to adjust endTime for sub-timestep mismatch , since in
                0311 C     several places, test for last iteration with time==endTime :
d987e9fc35 Jean*0312       tmpVar = startTime + deltaTClock*DFLOAT(nTimeSteps)
1cfd78e5a7 Jean*0313       IF ( endTime.NE.tmpVar ) THEN
                0314        IF ( ABS(endTime-tmpVar).GT.deltaTClock*1. _d -6 ) THEN
                0315         WRITE(msgBuf,'(A,A)') '** WARNING ** SET_PARMS: ',
                0316      &   '(endTime-baseTime) not multiple of time-step'
                0317         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0318      &                      SQUEEZE_RIGHT, myThid )
                0319         WRITE(msgBuf,'(2A,1PE20.13)') '** WARNING ** SET_PARMS: ',
                0320      &   'Previous endTime=', endTime
                0321         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0322      &                      SQUEEZE_RIGHT, myThid )
                0323         WRITE(msgBuf,'(2A,1PE20.13)') '** WARNING ** SET_PARMS: ',
                0324      &   'Adjusted endTime=', tmpVar
                0325         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0326      &                      SQUEEZE_RIGHT, myThid )
                0327        ENDIF
                0328        endTime = tmpVar
                0329       ENDIF
                0330
4e66ab0b67 Oliv*0331 #ifdef ALLOW_LONGSTEP
                0332       IF ( usePTRACERS ) THEN
                0333         CALL LONGSTEP_CHECK_ITERS(myThid)
                0334       ENDIF
                0335 #endif /* ALLOW_LONGSTEP */
                0336
aa6b2555c8 Jean*0337       IF ( OLx.LE.0 .OR. OLy.LE.0 ) THEN
                0338 C--   Overlap-size will be checked later in CONFIG_CHECK. This minimal
                0339 C     check here is just to set-up the model grid (INI_GRID) safely.
                0340         WRITE(msgBuf,'(2A)') 'SET_PARMS: ',
                0341      &      'model unusable with no overlap (OLx,OLy = 0)'
                0342         CALL PRINT_ERROR( msgBuf, myThid )
                0343         errCount = errCount + 1
                0344       ENDIF
                0345
                0346       IF ( errCount.GE.1 ) THEN
                0347         WRITE(msgBuf,'(A,I3,A)')
                0348      &       'SET_PARMS: detected', errCount,' fatal error(s)'
                0349         CALL PRINT_ERROR( msgBuf, myThid )
                0350         CALL ALL_PROC_DIE( 0 )
                0351         STOP 'ABNORMAL END: S/R SET_PARMS'
                0352       ENDIF
                0353
                0354 C--   After this point, main model parameters are not supposed to be modified.
809bdccbfc Jean*0355        WRITE(msgBuf,'(A,A)') 'SET_PARMS: done'
                0356        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0357      &                      SQUEEZE_RIGHT , 1)
                0358
                0359       _END_MASTER(myThid)
                0360
                0361 C--   Everyone else must wait for the parameters to be set
                0362       _BARRIER
                0363
                0364       RETURN
                0365       END