-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathphysconst.F90
executable file
·343 lines (274 loc) · 16.5 KB
/
physconst.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
module physconst
! Physical constants. Use CCSM shared values whenever available.
use shr_kind_mod, only: r8 => shr_kind_r8
use shr_const_mod, only: shr_const_g, shr_const_stebol, shr_const_tkfrz, &
shr_const_mwdair, shr_const_rdair, shr_const_mwwv, &
shr_const_latice, shr_const_latvap, shr_const_cpdair, &
shr_const_rhofw, shr_const_cpwv, shr_const_rgas, &
shr_const_karman, shr_const_pstd, shr_const_rhodair,&
shr_const_avogad, shr_const_boltz, shr_const_cpfw, &
shr_const_rwv, shr_const_zvir, shr_const_pi, &
shr_const_rearth, shr_const_sday, shr_const_cday, &
shr_const_spval, shr_const_omega, shr_const_cpvir, &
shr_const_tktrip
#ifdef COMMENT_OUT
use ppgrid, only: pcols, pver, pverp, begchunk, endchunk ! Dimensions and chunk bounds
#endif
implicit none
private
#ifdef COMMENT_OUT
public :: physconst_init
public :: physconst_readnl
public :: physconst_update
#endif
save
! Constants based off share code or defined in physconst
real(r8), public, parameter :: avogad = shr_const_avogad ! Avogadro s number (molecules/kmole)
real(r8), public, parameter :: boltz = shr_const_boltz ! Boltzman s constant (J/K/molecule)
real(r8), public, parameter :: cday = shr_const_cday ! sec in calendar day ~ sec
real(r8), public, parameter :: cpair = shr_const_cpdair ! specific heat of dry air (J/K/kg)
real(r8), public, parameter :: cpliq = shr_const_cpfw ! specific heat of fresh h2o (J/K/kg)
real(r8), public, parameter :: karman = shr_const_karman ! Von Karman constant
real(r8), public, parameter :: latice = shr_const_latice ! Latent heat of fusion (J/kg)
real(r8), public, parameter :: latvap = shr_const_latvap ! Latent heat of vaporization (J/kg)
real(r8), public, parameter :: pi = shr_const_pi ! 3.14...
real(r8), public, parameter :: pstd = shr_const_pstd ! Standard pressure (Pascals)
real(r8), public, parameter :: r_universal = shr_const_rgas ! Universal gas constant (J/K/kmol)
real(r8), public, parameter :: rhoh2o = shr_const_rhofw ! Density of liquid water (STP)
real(r8), public, parameter :: spval = shr_const_spval !special value
real(r8), public, parameter :: stebol = shr_const_stebol ! Stefan-Boltzmann s constant (W/m^2/K^4)
real(r8), public, parameter :: h2otrip = shr_const_tktrip ! Triple point temperature of water (K)
real(r8), public, parameter :: c0 = 2.99792458e8_r8 ! Speed of light in a vacuum (m/s)
real(r8), public, parameter :: planck = 6.6260755e-34_r8 ! Planck s constant (J.s)
! Molecular weights
real(r8), public, parameter :: mwco2 = 44._r8 ! molecular weight co2
real(r8), public, parameter :: mwn2o = 44._r8 ! molecular weight n2o
real(r8), public, parameter :: mwch4 = 16._r8 ! molecular weight ch4
real(r8), public, parameter :: mwf11 = 136._r8 ! molecular weight cfc11
real(r8), public, parameter :: mwf12 = 120._r8 ! molecular weight cfc12
real(r8), public, parameter :: mwo3 = 48._r8 ! molecular weight O3
real(r8), public, parameter :: mwso2 = 64._r8
real(r8), public, parameter :: mwso4 = 96._r8
real(r8), public, parameter :: mwh2o2 = 34._r8
real(r8), public, parameter :: mwdms = 62._r8
! modifiable physical constants for aquaplanet
real(r8), public :: gravit = shr_const_g ! gravitational acceleration (m/s**2)
real(r8), public :: sday = shr_const_sday ! sec in siderial day ~ sec
real(r8), public :: mwh2o = shr_const_mwwv ! molecular weight h2o
real(r8), public :: cpwv = shr_const_cpwv ! specific heat of water vapor (J/K/kg)
real(r8), public :: mwdry = shr_const_mwdair! molecular weight dry air
real(r8), public :: rearth = shr_const_rearth! radius of earth (m)
real(r8), public :: tmelt = shr_const_tkfrz ! Freezing point of water (K)
!--------------- Variables below here are derived from those above -----------------------
real(r8), public :: rga = 1._r8/shr_const_g ! reciprocal of gravit
real(r8), public :: ra = 1._r8/shr_const_rearth ! reciprocal of earth radius
real(r8), public :: omega = shr_const_omega ! earth rot ~ rad/sec
real(r8), public :: rh2o = shr_const_rwv ! Water vapor gas constant ~ J/K/kg
real(r8), public :: rair = shr_const_rdair ! Dry air gas constant ~ J/K/kg
real(r8), public :: epsilo = shr_const_mwwv/shr_const_mwdair ! ratio of h2o to dry air molecular weights
real(r8), public :: zvir = shr_const_zvir ! (rh2o/rair) - 1
real(r8), public :: cpvir = shr_const_cpvir ! CPWV/CPDAIR - 1.0
real(r8), public :: rhodair = shr_const_rhodair ! density of dry air at STP ~ kg/m^3
real(r8), public :: cappa = (shr_const_rgas/shr_const_mwdair)/shr_const_cpdair ! R/Cp
real(r8), public :: ez ! Coriolis expansion coeff -> omega/sqrt(0.375)
real(r8), public :: Cpd_on_Cpv = shr_const_cpdair/shr_const_cpwv
!--------------- Variables below here are for WACCM-X -----------------------
real(r8), public, dimension(:,:,:), pointer :: cpairv ! composition dependent specific heat at constant pressure
real(r8), public, dimension(:,:,:), pointer :: rairv ! composition dependent gas "constant"
real(r8), public, dimension(:,:,:), pointer :: cappav ! rairv/cpairv
real(r8), public, dimension(:,:,:), pointer :: mbarv ! composition dependent atmosphere mean mass
real(r8), public, dimension(:,:,:), pointer :: kmvis ! molecular viscosity kg/m/s
real(r8), public, dimension(:,:,:), pointer :: kmcnd ! molecular conductivity J/m/s/K
!--------------- Variables below here are for turbulent mountain stress -----------------------
real(r8), public :: tms_orocnst
real(r8), public :: tms_z0fac
#ifdef COMMENT_OUT
!================================================================================================
contains
!================================================================================================
subroutine physconst_init()
use abortutils, only: endrun
implicit none
integer :: ierr
!-------------------------------------------------------------------------------
! Allocate constituent dependent properties
!-------------------------------------------------------------------------------
allocate( cpairv(pcols,pver,begchunk:endchunk), &
rairv(pcols,pver,begchunk:endchunk), &
cappav(pcols,pver,begchunk:endchunk), &
mbarv(pcols,pver,begchunk:endchunk), &
kmvis(pcols,pverp,begchunk:endchunk), &
kmcnd(pcols,pverp,begchunk:endchunk), stat=ierr )
if ( ierr /= 0 ) call endrun('physconst: allocate failed in physconst_init')
!-------------------------------------------------------------------------------
! Initialize constituent dependent properties
!-------------------------------------------------------------------------------
cpairv(:pcols,:pver,begchunk:endchunk) = cpair
rairv(:pcols,:pver,begchunk:endchunk) = rair
cappav(:pcols,:pver,begchunk:endchunk) = rair/cpair
mbarv(:pcols,:pver,begchunk:endchunk) = mwdry
return
end subroutine physconst_init
!==============================================================================
! Read namelist variables.
subroutine physconst_readnl(nlfile)
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
use spmd_utils, only: masterproc
use abortutils, only: endrun
use cam_logfile, only: iulog
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'physconst_readnl'
logical newg, newsday, newmwh2o, newcpwv, newmwdry, newrearth, newtmelt
! Physical constants needing to be reset (ie. for aqua planet experiments)
namelist /physconst_nl/ cpwv, gravit, mwdry, mwh2o, rearth, sday, tmelt, tms_orocnst, tms_z0fac
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'physconst_nl', status=ierr)
if (ierr == 0) then
read(unitn, physconst_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(cpwv, 1, mpir8, 0, mpicom)
call mpibcast(gravit, 1, mpir8, 0, mpicom)
call mpibcast(mwdry, 1, mpir8, 0, mpicom)
call mpibcast(mwh2o, 1, mpir8, 0, mpicom)
call mpibcast(rearth, 1, mpir8, 0, mpicom)
call mpibcast(sday, 1, mpir8, 0, mpicom)
call mpibcast(tmelt, 1, mpir8, 0, mpicom)
call mpibcast(tms_orocnst, 1, mpir8, 0, mpicom)
call mpibcast(tms_z0fac, 1, mpir8, 0, mpicom)
#endif
newg = gravit .ne. shr_const_g
newsday = sday .ne. shr_const_sday
newmwh2o = mwh2o .ne. shr_const_mwwv
newcpwv = cpwv .ne. shr_const_cpwv
newmwdry = mwdry .ne. shr_const_mwdair
newrearth= rearth .ne. shr_const_rearth
newtmelt = tmelt .ne. shr_const_tkfrz
if (newg .or. newsday .or. newmwh2o .or. newcpwv .or. newmwdry .or. newrearth .or. newtmelt) then
if (masterproc) then
write(iulog,*)'****************************************************************************'
write(iulog,*)'*** New Physical Constant Values set via namelist ***'
write(iulog,*)'*** ***'
write(iulog,*)'*** Physical Constant Old Value New Value ***'
if (newg) write(iulog,*)'*** GRAVITY ',shr_const_g,gravit,'***'
if (newsday) write(iulog,*)'*** SDAY ',shr_const_sday,sday,'***'
if (newmwh2o) write(iulog,*)'*** MWH20 ',shr_const_mwwv,mwh2o,'***'
if (newcpwv) write(iulog,*)'*** CPWV ',shr_const_cpwv,cpwv,'***'
if (newmwdry) write(iulog,*)'*** MWDRY ',shr_const_mwdair,mwdry,'***'
if (newrearth) write(iulog,*)'*** REARTH ',shr_const_rearth,rearth,'***'
if (newtmelt) write(iulog,*)'*** TMELT ',shr_const_tkfrz,tmelt,'***'
write(iulog,*)'****************************************************************************'
end if
rga = 1._r8/gravit
ra = 1._r8/rearth
omega = 2.0_R8*pi/sday
cpvir = cpwv/cpair - 1._r8
epsilo = mwh2o/mwdry
! rair and rh2o have to be defined before any of the variables that use them
rair = r_universal/mwdry
rh2o = r_universal/mwh2o
cappa = rair/cpair
rhodair = pstd/(rair*tmelt)
zvir = (rh2o/rair)-1.0_R8
ez = omega / sqrt(0.375_r8)
Cpd_on_Cpv = cpair/cpwv
else
ez = omega / sqrt(0.375_r8)
end if
end subroutine physconst_readnl
!===============================================================================
subroutine physconst_update(mmr, t, cnst_mw_o, cnst_mw_o2, cnst_mw_h, cnst_mw_n, ixo, ixo2, ixh, ncnst, lchnk, ncol)
!-----------------------------------------------------------------------
! Update the physics "constants" that vary
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------------------------------------
integer, intent(in) :: ncnst ! number of constituents in mmr array
real(r8), intent(in) :: mmr(pcols,pver,ncnst) ! constituents q array from state structure
real(r8), intent(in) :: t(pcols,pver) ! temperature t array from state structure
real(r8), intent(in) :: cnst_mw_o, cnst_mw_o2, cnst_mw_h, cnst_mw_n ! Molecular weight of O, O2, H, and N
integer, intent(in) :: ixo, ixo2, ixh ! indices for O, O2, and H
integer, intent(in) :: lchnk ! Chunk number
integer, intent(in) :: ncol ! number of columns
!
!---------------------------Local storage-------------------------------------------------------------
integer :: i,k ! column,level,constituent indices
real(r8):: mmro, mmro2, mmrh, mmrn2 ! Mass mixing ratios of O, O2, H, and N
real(r8):: mbarvi, tint ! Mean mass, temperature, and specific heat on interface levels
real(r8):: dof1, dof2 ! Degress of freedom for cpairv calculation
real(r8):: kv1, kv2, kv3, kv4 ! Coefficients for kmvis calculation
real(r8):: kc1, kc2, kc3, kc4 ! Coefficients for kmcnd calculation
!--------------------------------------------
! Set constants needed for updates
!--------------------------------------------
dof1 = 5._r8
dof2 = 7._r8
kv1 = 4.03_r8
kv2 = 3.42_r8
kv3 = 3.9_r8
kv4 = 0.69_r8
kc1 = 56._r8
kc2 = 56._r8
kc3 = 75.9_r8
kc4 = 0.69_r8
!--------------------------------------------
! update cpairv, rairv, mbarv, and cappav
!--------------------------------------------
do k=1,pver
do i=1,ncol
mmro = mmr(i,k, ixo)
mmro2 = mmr(i,k, ixo2)
mmrh = mmr(i,k, ixh)
mmrn2 = 1._r8-mmro-mmro2-mmrh
mbarv(i,k,lchnk) = 1._r8/(mmro/cnst_mw_o + &
mmro2/cnst_mw_o2 + &
mmrn2/cnst_mw_n/2._r8 + &
mmrh/cnst_mw_h)
rairv(i,k,lchnk) = shr_const_rgas / mbarv(i,k,lchnk)
cpairv(i,k,lchnk) = 0.5_r8*shr_const_rgas* &
(dof1*mmro/cnst_mw_o+ &
dof2*mmro2/cnst_mw_o2+ &
dof2*mmrn2/cnst_mw_n/2._r8+ &
dof1*mmrh/cnst_mw_h)
cappav(i,k,lchnk) = rairv(i,k,lchnk)/cpairv(i,k,lchnk)
enddo
enddo
do k=2,pver
do i=1,ncol
mmro = .5_r8*(mmr(i,k-1, ixo)+mmr(i,k,ixo))
mmro2 = .5_r8*(mmr(i,k-1, ixo2)+mmr(i,k,ixo2))
mmrn2 = 1._r8-mmro-mmro2
mbarvi = .5_r8*(mbarv(i,k-1,lchnk)+mbarv(i,k,lchnk))
tint = .5_r8*(t(i,k-1)+t(i,k))
kmvis(i,k,lchnk) = (kv1*mmro2/cnst_mw_o2+ &
kv2*mmrn2/cnst_mw_n/2._r8+ &
kv3*mmro/cnst_mw_o)*mbarvi* &
tint**kv4 * 1.e-7_r8
kmcnd(i,k,lchnk) = (kc1*mmro2/cnst_mw_o2+ &
kc2*mmrn2/cnst_mw_n/2._r8+ &
kc3*mmro/cnst_mw_o)*mbarvi* &
tint**kc4 * 1.e-5_r8
enddo
enddo
do i=1,ncol
kmvis(i,1,lchnk) = 1.5_r8*kmvis(i,2,lchnk)-.5_r8*kmvis(i,3,lchnk)
kmcnd(i,1,lchnk) = 1.5_r8*kmcnd(i,2,lchnk)-.5_r8*kmcnd(i,3,lchnk)
kmvis(i,pverp,lchnk) = kmvis(i,pver,lchnk)
kmcnd(i,pverp,lchnk) = kmcnd(i,pver,lchnk)
enddo
end subroutine physconst_update
#endif
end module physconst