C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_readparms.F,v 1.23 2017/08/09 15:23:36 mlosch Exp $
C $Name: $
#include "THSICE_OPTIONS.h"
CBOP
C !ROUTINE: THSICE_READPARMS
C !INTERFACE:
SUBROUTINE THSICE_READPARMS( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R THSICE_READPARMS
C | o Routine to initialize THSICE parameters and constants
C *==========================================================*
C | Initialize Th-Sea-ICE parameters, read in data.ice
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "THSICE_PARAMS.h"
#ifdef ALLOW_MNC
#include "MNC_PARAMS.h"
#endif
#ifdef ALLOW_COST
#include "THSICE_COST.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C myThid :: My Thread Id. number
INTEGER myThid
CEOP
#ifdef ALLOW_THSICE
C === Local variables ===
C msgBuf :: Informational/error message buffer
C iUnit :: Work variable for IO unit number
CHARACTER*(MAX_LEN_MBUF) msgBuf
CHARACTER*(20) fmtR, fmtI, fmtL, fmtC
INTEGER iUnit
C-- Th-Sea-ICE parameter
NAMELIST //THSICE_CONST
& rhos, rhoi, rhosw, rhofw,
& cpIce, cpWater,
& kIce, kSnow,
& bMeltCoef, Lfresh, qsnow,
& albColdSnow, albWarmSnow, tempSnowAlb,
& albOldSnow, hNewSnowAge, snowAgTime,
& albIceMax, albIceMin, hAlbIce, hAlbSnow,
& i0swFrac, ksolar, dhSnowLin,
& saltIce, S_winton, mu_Tf,
& Tf0kel, Terrmax, nitMaxTsf,
& hIceMin, hiMax, hsMax, iceMaskMax, iceMaskMin,
& fracEnMelt, fracEnFreez, hThinIce, hThickIce, hNewIceMax
NAMELIST //THSICE_PARM01
& startIceModel, stepFwd_oceMxL, thSIce_calc_albNIR,
& thSIce_skipThermo, thSIce_deltaT, thSIce_dtTemp,
& ocean_deltaT, tauRelax_MxL, tauRelax_MxL_salt,
& hMxL_default, sMxL_default, vMxL_default,
& thSIce_diffK, thSIceAdvScheme, stressReduction,
& thSIceBalanceAtmFW,
& thSIce_taveFreq, thSIce_diagFreq, thSIce_monFreq,
& thSIce_tave_mnc, thSIce_snapshot_mnc, thSIce_mon_mnc,
& thSIce_pickup_read_mnc, thSIce_pickup_write_mnc,
& thSIceFract_InitFile, thSIceThick_InitFile,
& thSIceSnowH_InitFile, thSIceSnowA_InitFile,
& thSIceEnthp_InitFile, thSIceTsurf_InitFile
#ifdef ALLOW_COST
NAMELIST //THSICE_COST
& mult_thsice, thsice_cost_ice_flag
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF ( .NOT.useThSIce ) THEN
C- pkg THSICE is not used
_BEGIN_MASTER(myThid)
C- Track pkg activation status:
C print a (weak) warning if data.ice is found
CALL PACKAGES_UNUSED_MSG( 'useThSIce', ' ', 'ice' )
_END_MASTER(myThid)
RETURN
ENDIF
_BEGIN_MASTER(myThid)
WRITE(msgBuf,'(A)') ' THSICE_READPARMS: opening data.ice'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
CALL OPEN_COPY_DATA_FILE(
I 'data.ice', 'THSICE_READPARMS',
O iUnit,
I myThid )
C-- Default values (constants)
rhos = 330. _d 0
rhoi = 900. _d 0
rhosw = rhoConst
rhofw = rhoConstFresh
cpIce = 2106. _d 0
cpWater = HeatCapacity_Cp
kIce = 2.03 _d 0
kSnow = 0.30 _d 0
bMeltCoef=0.006 _d 0
Lfresh = 3.34 _d 5
qsnow = Lfresh
albColdSnow= 0.85 _d 0
albWarmSnow= 0.70 _d 0
tempSnowAlb= -10. _d 0
albOldSnow = 0.55 _d 0
albIceMax = 0.65 _d 0
albIceMin = 0.20 _d 0
hAlbIce = 0.50 _d 0
hAlbSnow = 0.30 _d 0
hNewSnowAge= 2. _d -3
snowAgTime = 50. _d 0 * 86400. _d 0
i0swFrac = 0.3 _d 0
ksolar = 1.5 _d 0
dhSnowLin= 0. _d 0
saltIce = 4. _d 0
S_winton = 1. _d 0
mu_Tf = 0.054 _d 0
Tf0kel = celsius2K
Terrmax = 5.0 _d -1
nitMaxTsf= 20
hIceMin = 1. _d -2
hiMax = 10. _d 0
hsMax = 10. _d 0
iceMaskMax = 1. _d 0
iceMaskMin = 0.1 _d 0
fracEnMelt = 0.4 _d 0
fracEnFreez= 0. _d 0
hThinIce = 0.2 _d 0
hThickIce = 2.5 _d 0
hNewIceMax = UNSET_RL
C-- Default values (parameters)
stepFwd_oceMxL = .FALSE.
thSIce_skipThermo = .FALSE.
thSIce_calc_albNIR = .FALSE.
startIceModel = 0
thSIce_deltaT = dTtracerLev(1)
thSIce_dtTemp = UNSET_RL
ocean_deltaT = dTtracerLev(1)
tauRelax_MxL = 0. _d 0
tauRelax_MxL_salt = UNSET_RL
hMxL_default = 50. _d 0
sMxL_default = 35. _d 0
vMxL_default = 5. _d -2
thSIce_diffK = 0. _d 0
thSIceAdvScheme = 0
stressReduction = 1. _d 0
IF ( useSEAICE ) stressReduction = 0. _d 0
thSIceBalanceAtmFW = 0
thSIce_taveFreq = taveFreq
thSIce_diagFreq = dumpFreq
thSIce_monFreq = monitorFreq
#ifdef ALLOW_MNC
thSIce_tave_mnc = timeave_mnc
thSIce_snapshot_mnc = snapshot_mnc
thSIce_mon_mnc = monitor_mnc
thSIce_pickup_read_mnc = pickup_read_mnc
thSIce_pickup_write_mnc = pickup_write_mnc
#else
thSIce_tave_mnc = .FALSE.
thSIce_snapshot_mnc = .FALSE.
thSIce_mon_mnc = .FALSE.
thSIce_pickup_read_mnc = .FALSE.
thSIce_pickup_write_mnc = .FALSE.
#endif
thSIceFract_InitFile = ' '
thSIceThick_InitFile = ' '
thSIceSnowH_InitFile = ' '
thSIceSnowA_InitFile = ' '
thSIceEnthp_InitFile = ' '
thSIceTsurf_InitFile = ' '
#ifdef ALLOW_COST
thsice_cost_ice_flag = 1
mult_thsice = 0. _d 0
#endif
C-- Read parameters from open data file
READ(UNIT=iUnit,NML=THSICE_CONST)
WRITE(msgBuf,'(A)') ' THSICE_READPARMS: read THSICE_CONST'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
READ(UNIT=iUnit,NML=THSICE_PARM01)
WRITE(msgBuf,'(A)') ' THSICE_READPARMS: read THSICE_PARM01'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
#ifdef ALLOW_COST
READ(UNIT=iUnit,NML=THSICE_COST)
WRITE(msgBuf,'(A)') ' THSICE_READPARMS: read THSICE_COST'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , 1)
#endif
C-- Close the open data file
#ifdef SINGLE_DISK_IO
CLOSE(iUnit)
#else
CLOSE(iUnit,STATUS='DELETE')
#endif /* SINGLE_DISK_IO */
C- neutral default:
IF ( hNewIceMax .EQ. UNSET_RL ) hNewIceMax = hiMax
C If using the same time step for both icetop temp solver
C and ice thickness/growth, use thSIce_deltaT value
IF ( thSIce_dtTemp .EQ. UNSET_RL ) thSIce_dtTemp=thSIce_deltaT
C- If undef, set salt relax to temperature relax
IF ( tauRelax_MxL_salt .EQ. UNSET_RL ) THEN
tauRelax_MxL_salt = tauRelax_MxL
ENDIF
C- Define other constants (from previous ones):
Tmlt1=-mu_Tf*S_winton
floodFac = (rhosw - rhoi)/rhos
C Set I/O parameters
thSIce_tave_mdsio = .TRUE.
thSIce_snapshot_mdsio = .TRUE.
thSIce_mon_stdio = .TRUE.
thSIce_pickup_write_mdsio = .TRUE.
#ifdef ALLOW_MNC
IF (useMNC) THEN
IF ( .NOT.outputTypesInclusive
& .AND. thSIce_tave_mnc ) thSIce_tave_mdsio = .FALSE.
IF ( .NOT.outputTypesInclusive
& .AND. thSIce_snapshot_mnc )
& thSIce_snapshot_mdsio = .FALSE.
IF ( .NOT.outputTypesInclusive
& .AND. thSIce_mon_mnc ) thSIce_mon_stdio = .FALSE.
IF ( .NOT.outputTypesInclusive
& .AND. thSIce_pickup_write_mnc )
& thSIce_pickup_write_mdsio = .FALSE.
ENDIF
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Check/reset parameters:
IF ( useSEAICE .AND. stressReduction.NE.0. _d 0 ) THEN
C-- If useSEAICE=.true., the stress is computed in seaice_model,
C-- so that it does not need any further reduction
WRITE(msgBuf,'(2A)') '** WARNING ** THSICE_READPARMS:',
& ' reset stressReduction to zero'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(2A)') 'THSICE_READPARMS: useSEAICE=T =>',
& ' stress is be computed by SEAICE pkg => no reduction'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
stressReduction = 0. _d 0
ENDIF
IF ( fluidIsAir .AND. thSIceBalanceAtmFW.NE.0 ) THEN
WRITE(msgBuf,'(2A)') '** WARNING ** THSICE_READPARMS:',
& ' reset thSIceBalanceAtmFW to zero'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(2A)') 'THSICE_READPARMS:',
& ' since it is not available in Atmospheric set-up'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
thSIceBalanceAtmFW = 0
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
iUnit = standardMessageUnit
fmtR = '(A6,A20,1PE21.13)'
fmtI = '(A6,A20,I10)'
fmtL = '(A6,A20,L10)'
fmtC = '(A26,1X,A)'
c CALL MDSFINDUNIT( iUnit, myThid )
c OPEN(iUnit,file='thsice_check_params',status='unknown')
WRITE(iUnit,fmtR) 'ThSI:', 'rhos =', rhos
WRITE(iUnit,fmtR) 'ThSI:', 'rhoi =', rhoi
WRITE(iUnit,fmtR) 'ThSI:', 'rhosw =', rhosw
WRITE(iUnit,fmtR) 'ThSI:', 'rhofw =', rhofw
WRITE(iUnit,fmtR) 'ThSI:', 'floodFac =', floodFac
WRITE(iUnit,fmtR) 'ThSI:', 'cpIce =', cpIce
WRITE(iUnit,fmtR) 'ThSI:', 'cpWater =', cpWater
WRITE(iUnit,fmtR) 'ThSI:', 'kIce =', kIce
WRITE(iUnit,fmtR) 'ThSI:', 'kSnow =', kSnow
WRITE(iUnit,fmtR) 'ThSI:', 'bMeltCoef =', bMeltCoef
WRITE(iUnit,fmtR) 'ThSI:', 'Lfresh =', Lfresh
WRITE(iUnit,fmtR) 'ThSI:', 'qsnow =', qsnow
WRITE(iUnit,fmtR) 'ThSI:', 'albColdSnow =', albColdSnow
WRITE(iUnit,fmtR) 'ThSI:', 'albWarmSnow =', albWarmSnow
WRITE(iUnit,fmtR) 'ThSI:', 'tempSnowAlb =', tempSnowAlb
WRITE(iUnit,fmtR) 'ThSI:', 'albOldSnow =', albOldSnow
WRITE(iUnit,fmtR) 'ThSI:', 'hNewSnowAge =', hNewSnowAge
WRITE(iUnit,fmtR) 'ThSI:', 'snowAgTime =', snowAgTime
WRITE(iUnit,fmtR) 'ThSI:', 'albIceMax =', albIceMax
WRITE(iUnit,fmtR) 'ThSI:', 'albIceMin =', albIceMin
WRITE(iUnit,fmtR) 'ThSI:', 'hAlbIce =', hAlbIce
WRITE(iUnit,fmtR) 'ThSI:', 'hAlbSnow =', hAlbSnow
WRITE(iUnit,fmtR) 'ThSI:', 'i0swFrac =', i0swFrac
WRITE(iUnit,fmtR) 'ThSI:', 'ksolar =', ksolar
WRITE(iUnit,fmtR) 'ThSI:', 'dhSnowLin =', dhSnowLin
WRITE(iUnit,fmtR) 'ThSI:', 'saltIce =', saltIce
WRITE(iUnit,fmtR) 'ThSI:', 'S_winton =', S_winton
WRITE(iUnit,fmtR) 'ThSI:', 'mu_Tf =', mu_Tf
WRITE(iUnit,fmtR) 'ThSI:', 'Tf0kel =', Tf0kel
WRITE(iUnit,fmtR) 'ThSI:', 'Tmlt1 =', Tmlt1
WRITE(iUnit,fmtR) 'ThSI:', 'Terrmax =', Terrmax
WRITE(iUnit,fmtI) 'ThSI:', 'nitMaxTsf =', nitMaxTsf
WRITE(iUnit,fmtR) 'ThSI:', 'hIceMin =', hIceMin
WRITE(iUnit,fmtR) 'ThSI:', 'hiMax =', hiMax
WRITE(iUnit,fmtR) 'ThSI:', 'hsMax =', hsMax
WRITE(iUnit,fmtR) 'ThSI:', 'iceMaskMax =', iceMaskMax
WRITE(iUnit,fmtR) 'ThSI:', 'iceMaskMin =', iceMaskMin
WRITE(iUnit,fmtR) 'ThSI:', 'fracEnMelt =', fracEnMelt
WRITE(iUnit,fmtR) 'ThSI:', 'fracEnFreez =', fracEnFreez
WRITE(iUnit,fmtR) 'ThSI:', 'hThinIce =', hThinIce
WRITE(iUnit,fmtR) 'ThSI:', 'hThickIce =', hThickIce
WRITE(iUnit,fmtR) 'ThSI:', 'hNewIceMax =', hNewIceMax
WRITE(iUnit,fmtR) 'ThSI:','stressReduction =',stressReduction
WRITE(iUnit,fmtL) 'ThSI:','thSIce_skipThermo =',thSIce_skipThermo
WRITE(iUnit,fmtI) 'ThSI:','thSIceAdvScheme =',thSIceAdvScheme
WRITE(iUnit,fmtI) 'ThSI:','thSIceBalanceAtmFW=',thSIceBalanceAtmFW
WRITE(iUnit,fmtR) 'ThSI:','thSIce_diffK =',thSIce_diffK
WRITE(iUnit,fmtR) 'ThSI:','thSIce_deltaT =',thSIce_deltaT
WRITE(iUnit,fmtR) 'ThSI:','ocean_deltaT =',ocean_deltaT
WRITE(iUnit,fmtL) 'ThSI:','stepFwd_oceMxL =',stepFwd_oceMxL
WRITE(iUnit,fmtR) 'ThSI:','tauRelax_MxL =',tauRelax_MxL
WRITE(iUnit,fmtR) 'ThSI:','tauRelax_MxL_salt =',tauRelax_MxL_salt
WRITE(iUnit,fmtR) 'ThSI:','hMxL_default =',hMxL_default
WRITE(iUnit,fmtR) 'ThSI:','sMxL_default =',sMxL_default
WRITE(iUnit,fmtR) 'ThSI:','vMxL_default =',vMxL_default
WRITE(iUnit,fmtR) 'ThSI:','thSIce_taveFreq =',thSIce_taveFreq
WRITE(iUnit,fmtR) 'ThSI:','thSIce_diagFreq =',thSIce_diagFreq
WRITE(iUnit,fmtR) 'ThSI:','thSIce_monFreq =',thSIce_monFreq
WRITE(iUnit,fmtI) 'ThSI:','startIceModel =',startIceModel
c WRITE(iUnit,fmtC) 'thSIceFract_InitFile =',thSIceFract_InitFile
IF (iUnit.NE.standardMessageUnit) CLOSE(iUnit)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
_END_MASTER(myThid)
C-- Everyone else must wait for the parameters to be loaded
_BARRIER
#endif /* ALLOW_THSICE */
RETURN
END