CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
maximum_hourly_diagnostics.F90
1
3
5
6 use machine, only: kind_phys
7
8 implicit none
9
10 private
11
12 public maximum_hourly_diagnostics_run
13
14 ! DH* TODO - cleanup use of constants
15 real(kind=kind_phys), parameter ::pq0=379.90516e0, a2a=17.2693882, a3=273.16, a4=35.86, rhmin=1.0e-6
16 ! *DH
17
18 ! Conversion from flashes per five minutes to flashes per minute.
19 real(kind=kind_phys), parameter :: scaling_factor = 0.2
20
21contains
22
23#if 0
24
27#endif
28 subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, &
29 imp_physics_gfdl, imp_physics_thompson, &
30 imp_physics_fer_hires, imp_physics_nssl, &
31 con_g, phil, &
32 gt0, refl_10cm, refdmax, refdmax263k, u10m, v10m, &
33 u10max, v10max, spd10max, pgr, t2m, q2m, t02max, &
34 t02min, rh02max, rh02min, dtp, rain, pratemax, &
35 lightning_threat, ltg1_max,ltg2_max,ltg3_max, &
36 wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,&
37 prsl, kdt, errmsg, errflg)
38
39 ! Interface variables
40 integer, intent(in) :: im, levs, kdt
41 logical, intent(in) :: reset, lradar, lightning_threat
42 integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, &
43 imp_physics_nssl
44 real(kind_phys), intent(in ) :: con_g
45 real(kind_phys), intent(in ) :: con_rd
46 real(kind_phys), intent(in ) :: phil(:,:)
47 real(kind_phys), intent(in ) :: gt0(:,:)
48 real(kind_phys), intent(in ) :: refl_10cm(:,:)
49 real(kind_phys), intent(inout) :: refdmax(:)
50 real(kind_phys), intent(inout) :: refdmax263k(:)
51 real(kind_phys), intent(in ) :: u10m(:)
52 real(kind_phys), intent(in ) :: v10m(:)
53 real(kind_phys), intent(inout) :: u10max(:)
54 real(kind_phys), intent(inout) :: v10max(:)
55 real(kind_phys), intent(inout) :: spd10max(:)
56 real(kind_phys), intent(in ) :: pgr(:)
57 real(kind_phys), intent(in ) :: t2m(:)
58 real(kind_phys), intent(in ) :: q2m(:)
59 real(kind_phys), intent(inout) :: t02max(:)
60 real(kind_phys), intent(inout) :: t02min(:)
61 real(kind_phys), intent(inout) :: rh02max(:)
62 real(kind_phys), intent(inout) :: rh02min(:)
63 real(kind_phys), intent(in ) :: dtp
64 real(kind_phys), intent(in ) :: rain(:)
65 real(kind_phys), intent(in ) :: tgrs(:,:)
66 real(kind_phys), intent(in ) :: prsl(:,:)
67 real(kind_phys), intent(inout) :: pratemax(:)
68
69 real(kind_phys), intent(in), dimension(:,:) :: prsi, qgraupel, qsnowwat, qicewat
70 real(kind_phys), intent(in), dimension(:,:), optional :: wgrs
71 real(kind_phys), intent(inout), dimension(:), optional :: ltg1_max, ltg2_max, ltg3_max
72 character(len=*), intent(out) :: errmsg
73 integer, intent(out) :: errflg
74
75 ! Local variables
76 real(kind_phys), dimension(:), allocatable :: refd, refd263k
77 real(kind_phys) :: tem, pshltr, qcq, rh02, dp, q
78 integer :: i
79
80 ! Initialize CCPP error handling variables
81 errmsg = ''
82 errflg = 0
83
84!Lightning threat indices
85 if (lightning_threat) then
86 call lightning_threat_indices
87 endif
88
89!Calculate hourly max 1-km agl and -10C reflectivity
90 if (lradar .and. (imp_physics == imp_physics_gfdl .or. &
91 imp_physics == imp_physics_thompson .or. &
92 imp_physics == imp_physics_fer_hires .or. &
93 imp_physics == imp_physics_nssl )) then
94 allocate(refd(im))
95 allocate(refd263k(im))
96 call max_fields(phil,refl_10cm,con_g,im,levs,refd,gt0,refd263k)
97 if (reset) then
98 IF ( imp_physics == imp_physics_nssl ) THEN ! ERM: might not need this as a separate assignment
99 do i=1,im
100 refdmax(i) = 0.
101 refdmax263k(i) = 0.
102 enddo
103 ELSE
104 do i=1,im
105 refdmax(i) = -35.
106 refdmax263k(i) = -35.
107 enddo
108 ENDIF
109 endif
110 do i=1,im
111 refdmax(i) = max(refdmax(i),refd(i))
112 refdmax263k(i) = max(refdmax263k(i),refd263k(i))
113 enddo
114 deallocate (refd)
115 deallocate (refd263k)
116 endif
117!
118 if (reset) then
119 do i=1,im
120 spd10max(i) = -999.
121 u10max(i) = -999.
122 v10max(i) = -999.
123 t02max(i) = -999.
124 t02min(i) = 999.
125 rh02max(i) = -999.
126 rh02min(i) = 999.
127 pratemax(i) = 0.
128 enddo
129 endif
130 do i=1,im
131! find max hourly wind speed then decompose
132 tem = sqrt(u10m(i)*u10m(i) + v10m(i)*v10m(i))
133 if (tem > spd10max(i)) then
134 spd10max(i) = tem
135 u10max(i) = u10m(i)
136 v10max(i) = v10m(i)
137 endif
138 pshltr=pgr(i)*exp(-0.068283/gt0(i,1))
139 qcq=pq0/pshltr*exp(a2a*(t2m(i)-a3)/(t2m(i)-a4))
140 rh02=q2m(i)/qcq
141 IF (rh02 > 1.0) THEN
142 rh02 = 1.0
143 ENDIF
144 IF (rh02 < rhmin) THEN !use smaller RH limit for stratosphere
145 rh02 = rhmin
146 ENDIF
147 rh02max(i) = max(rh02max(i),rh02)
148 rh02min(i) = min(rh02min(i),rh02)
149 t02max(i) = max(t02max(i),t2m(i))
150 t02min(i) = min(t02min(i),t2m(i))
151 pratemax(i) = max(pratemax(i),(3.6e6/dtp)*rain(i))
152 enddo
153
154 contains
155
156 subroutine lightning_threat_indices
157 implicit none
158 REAL(kind_phys), PARAMETER :: clim1=1.50
159 REAL(kind_phys), PARAMETER :: clim2=0.40*1.22
160 REAL(kind_phys), PARAMETER :: clim3=0.02*1.22
161 ! coef1 and coef2 are modified from the values given
162 ! in McCaul et al.
163 ! coef1 is x 1000 x 1.22
164 ! coef2 is x 1.22
165 ! are these tuning factors, scale factors??
166 ! McCaul et al. used a 2-km WRF simulation
167 REAL(kind_phys), PARAMETER :: coef1=0.042*1000.*1.22
168 REAL(kind_phys), PARAMETER :: coef2=0.20*1.22
169
170 REAL(kind_phys) :: totice_colint(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho
171 LOGICAL :: ltg1_calc(im)
172 integer :: k, i, count
173
174 count = 0
175 high_ltg1 = 0
176 high_wgrs = 0
177 high_graupel = 0
178
179 totice_colint = 0
180 ltg1_calc = .false.
181 do k=1,levs-1
182 do i=1,im
183 dp = prsi(i,k) - prsi(i,k+1)
184 q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k)
185 rho = prsl(i,k) / (con_rd * tgrs(i,k))
186 totice_colint(i) = totice_colint(i) + q * rho * dp / con_g
187
188 IF ( .not.ltg1_calc(i) ) THEN
189 IF ( 0.5*(tgrs(i,k+1) + tgrs(i,k)) < 258.15 ) THEN
190 count = count + 1
191 ltg1_calc(i) = .true.
192
193 ltg1 = coef1*wgrs(i,k)* &
194 (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 )
195 if(ltg1 > high_ltg1) then
196 high_ltg1 = ltg1
197 high_graupel = qgraupel(i,k)
198 high_wgrs = wgrs(i,k)
199 endif
200
201 IF ( ltg1 .LT. clim1 ) ltg1 = 0.
202
203 ! Scale to flashes per minue
204 ltg1 = ltg1 * scaling_factor
205
206 IF ( ltg1 .GT. ltg1_max(i) ) THEN
207 ltg1_max(i) = ltg1
208 ENDIF
209 ENDIF
210 ENDIF
211 enddo
212 enddo
213
214 do i=1,im
215 ltg2 = coef2 * totice_colint(i)
216
217 IF ( ltg2 .LT. clim2 ) ltg2 = 0.
218
219 ! Scale to flashes per minute
220 ltg2 = ltg2 * scaling_factor
221
222 IF ( ltg2 .GT. ltg2_max(i) ) THEN
223 ltg2_max(i) = ltg2
224 ENDIF
225
226 ! This calculation is already in flashes per minute.
227 ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i)
228
229 ! Thus, we must scale clim3. The compiler will optimize this away.
230 IF ( ltg3_max(i) .LT. clim3 * scaling_factor ) ltg3_max(i) = 0.
231 enddo
232
233 end subroutine lightning_threat_indices
234
235 end subroutine maximum_hourly_diagnostics_run
236
237 subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k)
238 integer, intent(in) :: im,levs
239 real (kind=kind_phys), intent(in) :: grav
240 real (kind=kind_phys), intent(in),dimension(:,:) :: phil,ref3d,tk
241 integer :: i,k,ll,ipt,kpt
242 real :: dbz1avg,zmidp1,zmidloc,refl,fact
243 real, dimension(im,levs) :: z
244 real, dimension(im) :: zintsfc
245 real, dimension(:), intent(inout) :: refd,refd263k
246 REAL :: dbz1(2),dbzk,dbzk1
247 logical :: counter
248 do i=1,im
249 do k=1,levs
250 z(i,k)=phil(i,k)/grav
251 enddo
252 enddo
253 do i=1,im
254 refd(i) = -35.
255 vloop: do k=1,levs-1
256 if ( z(i,k+1) >= 1000. .and. z(i,k) <= 1000.) then
257 zmidp1=z(i,k+1)
258 zmidloc=z(i,k)
259 dbz1(1)=ref3d(i,k+1) !- dBZ (not Z) values
260 dbz1(2)=ref3d(i,k) !- dBZ values
261 exit vloop
262 endif
263 enddo vloop
264
265!!! Initial curefl value without reduction above freezing level
266!
267! curefl=0.
268! if (cprate(i,j)>0.) then
269! cuprate=rdtphs*cprate(i,j)
270! curefl=cu_a*cuprate**cu_b
271! endif
272 do ll=1,2
273 refl=0.
274 if (dbz1(ll)>-35.) refl=10.**(0.1*dbz1(ll))
275! dbz1(l)=curefl+refl !- in Z units
276 dbz1(ll)=refl
277 enddo
278!-- Vertical interpolation of Z (units of mm**6/m**3)
279 fact=(1000.-zmidloc)/(zmidloc-zmidp1)
280 dbz1avg=dbz1(2)+(dbz1(2)-dbz1(1))*fact
281!-- Convert to dBZ (10*logZ) as the last step
282 if (dbz1avg>0.01) then
283 dbz1avg=10.*alog10(dbz1avg)
284 else
285 dbz1avg=-35.
286 endif
287 refd(i)=max(refd(i),dbz1avg)
288 enddo
289
290!-- refl at -10C
291 do i=1,im
292 dbz1(1) = -35.
293 dbz1(2) = -35.
294 vloopm10: do k=1,levs-1
295 if (tk(i,k+1) .le. 263.15 .and. tk(i,k) .ge. 263.15) then
296 dbz1(1)=ref3d(i,k+1) !- dBZ (not Z) values
297 dbz1(2)=ref3d(i,k) !- dBZ values
298 exit vloopm10
299 endif
300 enddo vloopm10
301
302 do ll=1,2
303 refl=0.
304 if (dbz1(ll)>-35.) refl=10.**(0.1*dbz1(ll))
305! dbz1(l)=curefl+refl !- in Z units
306 dbz1(ll)=refl
307 enddo
308!-- Take max of bounding reflectivity values
309 dbz1avg=maxval(dbz1)
310!-- Convert to dBZ (10*logZ) as the last step
311 if (dbz1avg>0.01) then
312 dbz1avg=10.*alog10(dbz1avg)
313 else
314 dbz1avg=-35.
315 endif
316 refd263k(i)=dbz1avg
317 enddo
318 end subroutine max_fields
319
subroutine count(slimsk, sno, ijmax)
This subroutine counts number of points for the four surface conditions.
Definition sfcsub.F:2631