UPP (develop)
Loading...
Searching...
No Matches
CALWXT_BOURG.f
Go to the documentation of this file.
1
54!--------------------------------------------------------------------------------------
81!------------------------------------------------------------------------------------------------------------
82 subroutine calwxt_bourg_post(im,ista_2l,iend_2u,ista,iend,jm, &
83 & jsta_2l,jend_2u,jsta,jend,lm,lp1, &
84 & iseed,g,pthresh, &
85 & t,q,pmid,pint,lmh,prec,zint,ptype,me)
86 implicit none
87!
88! input:
89 integer,intent(in):: im,jm,jsta_2l,jend_2u,jsta,jend,lm,lp1,iseed,me,&
90 ista_2l,iend_2u,ista,iend
91 real,intent(in):: g,pthresh
92 real,intent(in), dimension(ista_2l:iend_2u,jsta_2l:jend_2u,lm) :: t, q, pmid
93 real,intent(in), dimension(ista_2l:iend_2u,jsta_2l:jend_2u,lp1) :: pint, zint
94 real,intent(in), dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: lmh, prec
95!
96! output:
97! real,intent(out) :: ptype(im,jm)
98 integer,intent(out) :: ptype(ista:iend,jsta:jend)
99!
100 integer i,j,ifrzl,iwrml,l,lhiwrm,lmhk,jlen
101 real pintk1,areane,tlmhk,areape,pintk2,surfw,area1,dzkl,psfck,r1,r2
102 real rn(im*jm*2)
103 integer :: rn_seed_size
104 integer, allocatable, dimension(:) :: rn_seed
105 logical, parameter :: debugprint = .false.
106!
107! initialize weather type array to zero (ie, off).
108! we do this since we want ptype to represent the
109! instantaneous weather type on return.
110 if (debugprint) then
111 print *,'in calwxtbg, jsta,jend=',jsta,jend,' im=',im
112 print *,'in calwxtbg,me=',me,'iseed=',iseed
113 endif
114!
115!$omp parallel do
116 do j=jsta,jend
117 do i=ista,iend
118 ptype(i,j) = 0
119 enddo
120 enddo
121!
122 jlen = jend - jsta + 1
123
124 call random_seed(size = rn_seed_size)
125 allocate(rn_seed(rn_seed_size))
126 rn_seed = iseed
127 call random_seed(put = rn_seed)
128 call random_number(rn)
129!
130!!$omp parallel do &
131! & private(a,lmhk,tlmhk,iwrml,psfck,lhiwrm,pintk1,pintk2,area1, &
132! & areape,dzkl,surfw,r1,r2)
133! print *,'incalwxtbg, rn',maxval(rn),minval(rn)
134
135 do j=jsta,jend
136! if(me==1)print *,'incalwxtbg, j=',j
137 do i=ista,iend
138 lmhk = min(nint(lmh(i,j)),lm)
139 psfck = pint(i,j,lmhk+1)
140!
141 if (prec(i,j) <= pthresh) cycle ! skip this point if no precip this time step
142
143! find the depth of the warm layer based at the surface
144! this will be the cut off point between computing
145! the surface based warm air and the warm air aloft
146!
147 tlmhk = t(i,j,lmhk) ! lowest layer t
148 iwrml = lmhk + 1
149 if (tlmhk >= 273.15) then
150 do l = lmhk, 2, -1
151 if (t(i,j,l) >= 273.15 .and. t(i,j,l-1) < 273.15 .and. &
152 & iwrml == lmhk+1) iwrml = l
153 end do
154 end if
155!
156! now find the highest above freezing level
157!
158! gsm added 250 mb check to prevent stratospheric warming situations
159! from counting as warm layers aloft
160 lhiwrm = lmhk + 1
161 do l = lmhk, 1, -1
162 if (t(i,j,l) >= 273.15 .and. pmid(i,j,l) > 25000.) lhiwrm = l
163 end do
164
165! energy variables
166
167! surfw is the positive energy between ground and the first sub-freezing layer above ground
168! areane is the negative energy between ground and the highest layer above ground
169! that is above freezing
170! areape is the positive energy "aloft" which is the warm energy not based at the ground
171! (the total warm energy = surfw + areape)
172!
173! pintk1 is the pressure at the bottom of the layer
174! pintk2 is the pressure at the top of the layer
175! dzkl is the thickness of the layer
176! ifrzl is a flag that tells us if we have hit a below freezing layer
177!
178 pintk1 = psfck
179 ifrzl = 0
180 areane = 0.0
181 areape = 0.0
182 surfw = 0.0
183
184 do l = lmhk, 1, -1
185 if (ifrzl == 0.and.t(i,j,l) <= 273.15) ifrzl = 1
186 pintk2 = pint(i,j,l)
187 dzkl = zint(i,j,l)-zint(i,j,l+1)
188 area1 = log(t(i,j,l)/273.15) * g * dzkl
189 if (t(i,j,l) >= 273.15.and. pmid(i,j,l) > 25000.) then
190 if (l < iwrml) areape = areape + area1
191 if (l >= iwrml) surfw = surfw + area1
192 else
193 if (l > lhiwrm) areane = areane + abs(area1)
194 end if
195 pintk1 = pintk2
196 end do
197!
198! decision tree time
199!
200 if (areape < 2.0) then
201! very little or no positive energy aloft, check for
202! positive energy just above the surface to determine rain vs. snow
203 if (surfw < 5.6) then
204! not enough positive energy just above the surface
205! snow = 1
206 ptype(i,j) = 1
207 else if (surfw > 13.2) then
208! enough positive energy just above the surface
209! rain = 8
210 ptype(i,j) = 8
211 else
212! transition zone, assume equally likely rain/snow
213! picking a random number, if <=0.5 snow
214 r1 = rn(i+im*(j-1))
215 if (r1 <= 0.5) then
216 ptype(i,j) = 1 ! snow = 1
217 else
218 ptype(i,j) = 8 ! rain = 8
219 end if
220 end if
221!
222 else
223! some positive energy aloft, check for enough negative energy
224! to freeze and make ice pellets to determine ip vs. zr
225 if (areane > 66.0+0.66*areape) then
226! enough negative area to make ip,
227! now need to check if there is enough positive energy
228! just above the surface to melt ip to make rain
229 if (surfw < 5.6) then
230! not enough energy at the surface to melt ip
231 ptype(i,j) = 2 ! ice pellets = 2
232 else if (surfw > 13.2) then
233! enough energy at the surface to melt ip
234 ptype(i,j) = 8 ! rain = 8
235 else
236! transition zone, assume equally likely ip/rain
237! picking a random number, if <=0.5 ip
238 r1 = rn(i+im*(j-1))
239 if (r1 <= 0.5) then
240 ptype(i,j) = 2 ! ice pellets = 2
241 else
242 ptype(i,j) = 8 ! rain = 8
243 end if
244 end if
245 else if (areane < 46.0+0.66*areape) then
246! not enough negative energy to refreeze, check surface temp
247! to determine rain vs. zr
248 if (tlmhk < 273.15) then
249 ptype(i,j) = 4 ! freezing rain = 4
250 else
251 ptype(i,j) = 8 ! rain = 8
252 end if
253 else
254! transition zone, assume equally likely ip/zr
255! picking a random number, if <=0.5 ip
256 r1 = rn(i+im*(j-1))
257 if (r1 <= 0.5) then
258! still need to check positive energy
259! just above the surface to melt ip vs. rain
260 if (surfw < 5.6) then
261 ptype(i,j) = 2 ! ice pellets = 2
262 else if (surfw > 13.2) then
263 ptype(i,j) = 8 ! rain = 8
264 else
265! transition zone, assume equally likely ip/rain
266! picking a random number, if <=0.5 ip
267 r2 = rn(i+im*(j-1)+im*jm)
268 if (r2 <= 0.5) then
269 ptype(i,j) = 2 ! ice pellets = 2
270 else
271 ptype(i,j) = 8 ! rain = 8
272 end if
273 end if
274 else
275! not enough negative energy to refreeze, check surface temp
276! to determine rain vs. zr
277 if (tlmhk < 273.15) then
278 ptype(i,j) = 4 ! freezing rain = 4
279 else
280 ptype(i,j) = 8 ! rain = 8
281 end if
282 end if
283 end if
284 end if
285! write(1000+me,*)' finished for i, j, from calbourge me=',me,i,j
286 end do
287! write(1000+me,*)' finished for j, from calbourge me=',me,j
288 end do
289! write(1000+me,*)' returning from calbourge me=',me
290 return
291 end
subroutine calwxt_bourg_post(im, ista_2l, iend_2u, ista, iend, jm, jsta_2l, jend_2u, jsta, jend, lm, lp1, iseed, g, pthresh, t, q, pmid, pint, lmh, prec, zint, ptype, me)
calwxt_bourg_post Subroutine that calculates precipitation type (Bourgouin).