forked from NOAA-PSL/stochastic_physics
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcompns_stochy.F90
201 lines (191 loc) · 7.59 KB
/
compns_stochy.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
!-----------------------------------------------------------------------
subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
!$$$ Subprogram Documentation Block
!
! Subprogram: compns Check and compute namelist frequencies
! Prgmmr: Iredell Org: NP23 Date: 1999-01-26
!
! Abstract: This subprogram checks global spectral model namelist
! frequencies in hour units for validity. If they are valid,
! then the frequencies are computed in timestep units.
! The following rules are applied:
! 1. the timestep must be positive;
!
! Program History Log:
! 2016-10-11 Phil Pegion make the stochastic physics stand alone
!
! Usage: call compns_stochy (me,deltim,nlunit, stochy_namelist,iret)
! Input Arguments:
! deltim - real timestep in seconds
! Output Arguments:
! iret - integer return code (0 if successful or
! between 1 and 8 for which rule above was broken)
! stochy_namelist
!
! Attributes:
! Language: Fortran 90
!
!$$$
use stochy_namelist_def
implicit none
integer, intent(out) :: iret
integer, intent(in) :: nlunit,me,sz_nml
character(len=*), intent(in) :: input_nml_file(sz_nml)
character(len=64), intent(in) :: fn_nml
real, intent(in) :: deltim
real tol
integer k,ios
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
namelist /nam_stochy/ntrunc,lon_s,lat_s,sppt,sppt_tau,sppt_lscale,sppt_logit, &
iseed_shum,iseed_sppt,shum,shum_tau,&
shum_lscale,fhstoch,stochini,skeb_varspect_opt,sppt_sfclimit, &
skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, &
skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,&
shum_sigefold,skebint,skeb_npass,use_zmtnblck
namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts
pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land
tol=0.01 ! tolerance for calculations
! spectral resolution defintion
ntrunc=-999
lon_s=-999
lat_s=-999
! can specify up to 5 values for the stochastic physics parameters
! (each is an array of length 5)
sppt = -999. ! stochastic physics tendency amplitude
shum = -999. ! stochastic boundary layer spf hum amp
skeb = -999. ! stochastic KE backscatter amplitude
! mg, sfcperts
pertz0 = -999. ! momentum roughness length amplitude
pertshc = -999. ! soil hydraulic conductivity amp
pertzt = -999. ! mom/heat roughness length amplitude
pertlai = -999. ! leaf area index amplitude
pertvegf = -999. ! vegetation fraction amplitude
pertalb = -999. ! albedo perturbations amplitude
! logicals
do_sppt = .false.
use_zmtnblck = .false.
do_shum = .false.
do_skeb = .false.
! mg, sfcperts
do_sfcperts = .false.
sppt_land = .false.
nsfcpert = 0
! for sfcperts random patterns
sfc_lscale = -999. ! length scales
sfc_tau = -999. ! time scales
iseed_sfc = 0 ! random seeds (if 0 use system clock)
! for SKEB random patterns.
skeb_vfilt = 0
skebint = 0
skeb_npass = 11 ! number of passes of smoother for dissipation estiamte
sppt_tau = -999. ! time scales
shum_tau = -999.
skeb_tau = -999.
skeb_vdof = 5 ! proxy for vertical correlation, 5 is close to 40 passes of the 1-2-1 filter in the GFS
skebnorm = 0 ! 0 - random pattern is stream function, 1- pattern is kenorm, 2- pattern is vorticity
sppt_lscale = -999. ! length scales
shum_lscale = -999.
skeb_lscale = -999.
iseed_sppt = 0 ! random seeds (if 0 use system clock)
iseed_shum = 0
iseed_skeb = 0
! parameters to control vertical tapering of stochastic physics with
! height
sppt_sigtop1 = 0.1
sppt_sigtop2 = 0.025
skeb_sigtop1 = 0.1
skeb_sigtop2 = 0.025
shum_sigefold = 0.2
! reduce amplitude of sppt near surface (lowest 2 levels)
sppt_sfclimit = .false.
! gaussian or power law variance spectrum for skeb (0: gaussian, 1:
! power law). If power law, skeb_lscale interpreted as a power not a
! length scale.
skeb_varspect_opt = 0
sppt_logit = .false. ! logit transform for sppt to bounded interval [-1,+1]
fhstoch = -999.0 ! forecast hour to dump random patterns
stochini = .false. ! true= read in pattern, false=initialize from seed
#ifdef INTERNAL_FILE_NML
read(input_nml_file, nml=nam_stochy)
#else
rewind (nlunit)
open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios)
read(nlunit,nam_stochy)
#endif
#ifdef INTERNAL_FILE_NML
read(input_nml_file, nml=nam_sfcperts)
#else
rewind (nlunit)
open (unit=nlunit, file=fn_nml, READONLY, status='OLD', iostat=ios)
read(nlunit,nam_sfcperts)
#endif
if (me == 0) then
print *,' in compns_stochy'
print*,'skeb=',skeb
endif
! PJP stochastic physics additions
IF (sppt(1) > 0 ) THEN
do_sppt=.true.
ENDIF
IF (shum(1) > 0 ) THEN
do_shum=.true.
! shum parameter has units of 1/hour, to remove time step
! dependence.
! change shum parameter units from per hour to per timestep
DO k=1,5
IF (shum(k) .gt. 0.0) shum(k)=shum(k)*deltim/3600.0
ENDDO
ENDIF
IF (skeb(1) > 0 ) THEN
do_skeb=.true.
if (skebnorm==0) then ! stream function norm
skeb=skeb*1.111e3*sqrt(deltim)
!skeb=skeb*5.0e5/sqrt(deltim)
endif
if (skebnorm==1) then ! stream function norm
skeb=skeb*0.00222*sqrt(deltim)
!skeb=skeb*1/sqrt(deltim)
endif
if (skebnorm==2) then ! vorticty function norm
skeb=skeb*1.111e-9*sqrt(deltim)
!skeb=skeb*5.0e-7/sqrt(deltim)
endif
! adjust skeb values for resolution.
! scaling is such that a value of 1.0 at T574 with a 900 second
! timestep produces well-calibrated values of forecast spread.
! DO k=1,5
! IF (skeb(k) .gt. 0.0) THEN
! skeb(k)=skeb(k)*deltim/(ntrunc*(ntrunc+1))*365765.0 ! 365765 is a scale factor so the base SKEB value in the namelist is 1.0
! skeb(k)=skeb(k)*deltim/(ntrunc*(ntrunc+1))*2000.0 ! 2000 is new scale factor so the base SKEB value in the namelist is 1.0
! ENDIF
! ENDDO
ENDIF
! compute frequencty to estimate dissipation timescale
IF (skebint == 0.) skebint=deltim
nsskeb=nint(skebint/deltim) ! skebint in seconds
IF(nsskeb<=0 .or. abs(nsskeb-skebint/deltim)>tol) THEN
WRITE(0,*) "SKEB interval is invalid",skebint
iret=9
return
ENDIF
! mg, sfcperts
IF (pertz0(1) > 0 .OR. pertshc(1) > 0 .OR. pertzt(1) > 0 .OR. &
pertlai(1) > 0 .OR. pertvegf(1) > 0 .OR. pertalb(1) > 0) THEN
do_sfcperts=.true.
ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! All checks are successful.
!
if (me == 0) then
print *, 'stochastic physics'
print *, ' do_sppt : ', do_sppt
print *, ' do_shum : ', do_shum
print *, ' do_skeb : ', do_skeb
print *, ' do_sfcperts : ', do_sfcperts
endif
iret = 0
!
return
end subroutine compns_stochy