Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rrpr.F90 @ 87

Last change on this file since 87 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 24.4 KB

Line
1	subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix)
2	! !
3	! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite" !
4	! !
5	!*****************************************************************************!
6	! !
7	! Recent changes: !
8	! !
9	! 2004-10-29: !
10	! - Sync rrpr.F WRFSI v2.0.1 with MM5 v3.5 !
11	! Added DATELEN: length of date strings to use for !
12	! our output file names. !
13	! Added MM5 interpolation from surrounding levels !
14	! if upper-air U or V are missing !
15	! Added test to see if we've got a SEAICE field, !
16	! make sure that it is all Zeros and Ones: !
17	! !
18	! 2002-05-16: !
19	! - Handle the Mercator projection. !
20	! This change also required changes to output.F, rd_grib.F, !
21	! datint.F, gribcode.F !
22	! !
23	! 2002-02-13: !
24	! - Added vertical interpolation in pressure in case of missing !
25	! U, V, T (the check for RH was already there) !
26	! !
27	! 2001-02-14: !
28	! - Allow file names to have date stamps out to minutes or !
29	! seconds, if the user requests a time interval (in seconds) !
30	! that is not evenly divisible into hours or minutes. !
31	! INTERVAL is checked for divisibility into 3600 (for hours) !
32	! or 60 (for minutes). The local variable DATELEN is set !
33	! to be the number of characters to use in our character !
34	! dates. Valid values for DATELEN are 13 (for hours), !
35	! 16 (for minutes), and 19 (for seconds). !
36	! !
37	! This change also requires changes to pregrid_grib.F, !
38	! output.F, datint.F, file_delete.F !
39	! !
40	! - Do processing not just if the requested date matches one we !
41	! want, but if the requested date falls between the startdate !
42	! and the enddate. !
43	! !
44	!*****************************************************************************!
45
46	use filelist
47	use gridinfo
48	use storage_module
49	use table
50	use module_debug
51	use stringutil
52
53	implicit none
54
55	!------------------------------------------------------------------------------
56	! Arguments:
57
58	! HSTART: Starting date of times to process
59	character (LEN=19) :: hstart
60
61	! NTIMES: Number of time periods to process
62	integer :: ntimes
63
64	! INTERVAL: Time inteval (seconds) of time periods to process.
65	integer :: interval
66
67	! NLVL: The number of levels in the stored data.
68	integer :: nlvl
69
70	! MAXLVL: The parameterized maximum number of levels to allow.
71	integer :: maxlvl
72
73	! PLVL: Array of pressure levels (Pa) in the dataset
74	real , dimension(maxlvl) :: plvl
75
76	! DEBUG_LEVEL: Integer level of debug printing (from namelist)
77	integer :: debug_level
78
79	!------------------------------------------------------------------------------
80
81	character (LEN=25) :: units
82	character (LEN=46) :: Desc
83	real, allocatable, dimension(:,:) :: scr2d
84	real, pointer, dimension(:,:) :: ptr2d
85
86	integer :: k, kk, m, n, ierr, ifv
87	integer :: iunit=13
88
89	character(LEN=19) :: hdate, hend
90	character(LEN=24) :: hdate_output
91	character(LEN=3) :: out_format
92	character(LEN=256) :: prefix
93	real :: xfcst, level
94	character(LEN=9) :: field
95
96	integer :: ntime, idts
97
98	! DATELEN: length of date strings to use for our output file names.
99	integer :: datelen
100
101	! Decide the length of date strings to use for output file names.
102	! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
103	if (mod(interval,3600) == 0) then
104	datelen = 13
105	else if (mod(interval, 60) == 0) then
106	datelen = 16
107	else
108	datelen = 19
109	endif
110
111	if ( debug_level .gt. 100 ) then
112	call mprintf(.true.,DEBUG,"Begin rrpr")
113	call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
114	do n = 1, nfiles
115	call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
116	enddo
117	endif
118
119	! Compute the ending time:
120
121	call geth_newdate(hend, hstart, interval*ntimes)
122
123	call clear_storage
124
125	! We want to do something for each of the requested times:
126	TIMELOOP : do ntime = 1, ntimes
127	idts = (ntime-1) * interval
128	call geth_newdate(hdate, hstart, idts)
129	call mprintf(.true.,DEBUG, &
130	"RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
131
132	! Loop over the output file dates, and do stuff if the file date matches
133	! the requested time we are working on now.
134
135	FILELOOP : do n = 1, nfiles
136	if ( debug_level .gt. 100 ) then
137	call mprintf(.true.,DEBUG, &
138	"hstart = %s , hend = %s",s1=hstart,s2=hend)
139	call mprintf(.true.,DEBUG, &
140	"filedates(n) = %s",s1=filedates(n))
141	call mprintf(.true.,DEBUG, &
142	"filedates(n) = %s",s1=filedates(n)(1:datelen))
143	end if
144	if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
145	if (debug_level .gt. 50 ) then
146	call mprintf(.true.,INFORM, &
147	"RRPR Processing : %s",s1=filedates(n)(1:datelen))
148	endif
149	open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
150	form='unformatted',status='old')
151
152	! Read the file:
153
154	rdloop: do
155	read (iunit, iostat=ierr) ifv
156	if (ierr.ne.0) exit rdloop
157	if ( ifv .eq. 5) then ! WPS
158	read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
159	level, map%nx, map%ny, map%igrid
160	hdate = hdate_output(1:19)
161	select case (map%igrid)
162	case (0, 4)
163	read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
164	case (3)
165	read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
166	map%truelat1, map%truelat2, map%r_earth
167	case (5)
168	read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
169	map%truelat1, map%r_earth
170	case (1)
171	read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
172	map%truelat1, map%r_earth
173	case default
174	call mprintf(.true.,ERROR, &
175	"Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
176	end select
177	read (iunit) map%grid_wind
178
179	else if ( ifv .eq. 4 ) then ! SI
180	read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
181	map%nx, map%ny, map%igrid
182	hdate = hdate_output(1:19)
183	select case (map%igrid)
184	case (0, 4)
185	read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
186	case (3)
187	read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
188	map%lov, map%truelat1, map%truelat2
189	case (5)
190	read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
191	map%lov, map%truelat1
192	case default
193	call mprintf(.true.,ERROR, &
194	"Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
195	end select
196
197	else if ( ifv .eq. 3 ) then ! MM5
198	read(iunit) hdate_output, xfcst, field, units, desc, level,&
199	map%nx, map%ny, map%igrid
200	hdate = hdate_output(1:19)
201	select case (map%igrid)
202	case (3) ! lamcon
203	read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
204	map%truelat1, map%truelat2
205	case (5) ! Polar Stereographic
206	read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
207	map%truelat1
208	case (0, 4) ! lat/lon
209	read (iunit) map%lat1, map%lon1, map%dy, map%dx
210	case (1) ! Mercator
211	read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
212	case default
213	call mprintf(.true.,ERROR, &
214	"Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
215	end select
216	else
217	call mprintf(.true.,ERROR, &
218	"unknown out_format, ifv = %i",i1=ifv)
219	endif
220
221	allocate(ptr2d(map%nx,map%ny))
222	read (iunit) ptr2d
223	call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
224	nullify (ptr2d)
225	enddo rdloop
226	!
227	! We have reached the end of file, so time to close it.
228	!
229	close(iunit)
230	if (debug_level .gt. 100 ) call print_storage
231	!
232	! By now the file has been read completely. Now, see if we need to fill in
233	! missing fields:
234	!
235
236	! Retrieve the number of levels in storage:
237	!
238	call get_plvls(plvl, maxlvl, nlvl)
239	!
240	! Fill the surface level (code 200100) from higher 200100s, as necessary
241	!
242	do k = 1, nlvl
243	if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
244	! We found a level between 200100 and 200200, now find the field
245	! corresponding to that level.
246	MLOOP : do m = 1, maxvar
247	if (is_there(nint(plvl(k)), namvar(m))) then
248	INLOOP : do kk = 200101, nint(plvl(k))
249	if (is_there(kk, namvar(m))) then
250	if ( debug_level .gt. 100 ) then
251	call mprintf(.true.,DEBUG, &
252	"Copying %s at level %i to level 200100.",s1=namvar(m),i1=kk)
253	end if
254	call get_dims(kk, namvar(m))
255	allocate(scr2d(map%nx,map%ny))
256	call get_storage &
257	(kk, namvar(m), scr2d, map%nx, map%ny)
258	call put_storage &
259	(200100,namvar(m), scr2d,map%nx,map%ny)
260	deallocate(scr2d)
261	EXIT INLOOP
262	endif
263	enddo INLOOP
264	endif
265	enddo MLOOP
266	endif
267	enddo
268
269	!
270	! If upper-air U is missing, see if we can interpolate from surrounding levels.
271	! This is a simple vertical interpolation, linear in pressure.
272	! Currently, this simply fills in one missing level between two present levels.
273	!
274
275	do k = 2, nlvl-1, 1
276	if (plvl(k-1) .lt. 200000.) then
277	if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
278	( is_there(nint(plvl(k-1)), 'UU')) .and.&
279	( is_there(nint(plvl(k+1)), 'UU')) ) then
280	call get_dims(nint(plvl(k+1)), 'UU')
281	call vntrp(plvl, maxlvl, k, "UU ", map%nx, map%ny)
282	endif
283	endif
284	enddo
285
286	!
287	! If upper-air V is missing, see if we can interpolate from surrounding levels.
288	! This is a simple vertical interpolation, linear in pressure.
289	! Currently, this simply fills in one missing level between two present levels.
290	!
291
292	do k = 2, nlvl-1, 1
293	if (plvl(k-1) .lt. 200000.) then
294	if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
295	( is_there(nint(plvl(k-1)), 'VV')) .and.&
296	( is_there(nint(plvl(k+1)), 'VV')) ) then
297	call get_dims(nint(plvl(k+1)), 'VV')
298	call vntrp(plvl, maxlvl, k, "VV ", map%nx, map%ny)
299	endif
300	endif
301	enddo
302
303	!--- Tanya's change for initializing WRF with RUC
304	! This allows for the ingestion for RUC isentropic data
305	!
306	do k = 1, nlvl
307	if (plvl(k).lt.200000.) then
308	if (.not. is_there(nint(plvl(k)), 'TT').and. &
309	is_there(nint(plvl(k)), 'VPTMP')) then
310	call get_dims(nint(plvl(k)), 'VPTMP')
311	call compute_t_vptmp(map%nx, map%ny, plvl(k))
312	endif
313	endif
314	enddo
315	!!!
316	!
317	! If upper-air T is missing, see if we can interpolate from surrounding levels.
318	! This is a simple vertical interpolation, linear in pressure.
319	! Currently, this simply fills in one missing level between two present levels.
320	!
321
322	do k = 2, nlvl-1, 1
323	if (plvl(k-1) .lt. 200000.) then
324	if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
325	( is_there(nint(plvl(k-1)), 'TT')) .and.&
326	( is_there(nint(plvl(k+1)), 'TT')) ) then
327	call get_dims(nint(plvl(k+1)), 'TT')
328	call vntrp(plvl, maxlvl, k, "TT ", map%nx, map%ny)
329	endif
330	endif
331	enddo
332
333	!
334	! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
335	!--- Tanya's change for initializing WRF with RUC
336
337	do k = 1, nlvl
338	if (plvl(k).lt.200000.) then
339	if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
340	is_there(nint(plvl(k)), 'QV')) then
341	call get_dims(nint(plvl(k)), 'QV')
342	call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
343	endif
344	endif
345	enddo
346
347	!
348	! Check to see if we need to fill HGT from GEOPT.
349	!
350	do k = 1, nlvl
351	if (plvl(k).lt.200000.) then
352	if (.not. is_there(nint(plvl(k)), 'HGT').and. &
353	is_there(nint(plvl(k)), 'GEOPT')) then
354	call get_dims(nint(plvl(k)), 'GEOPT')
355	allocate(scr2d(map%nx,map%ny))
356	call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
357	scr2d = scr2d / 9.81
358	call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny)
359	deallocate(scr2d)
360	endif
361	endif
362	enddo
363
364
365	! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
366
367	do k = 1, nlvl
368	if (plvl(k).lt.200000.) then
369	if (.not. is_there(nint(plvl(k)), 'RH').and. &
370	is_there(nint(plvl(k)), 'SPECHUMD')) then
371	call get_dims(nint(plvl(k)), 'TT')
372	call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
373	endif
374	endif
375	enddo
376
377	! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
378	! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
379
380	do k = 1, nlvl
381	if (plvl(k).lt.200000.) then
382	if (.not. is_there(nint(plvl(k)),'RH').and. &
383	is_there(nint(plvl(k)),'VAPP')) then
384	call get_dims(nint(plvl(k)),'TT')
385	call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
386	endif
387	endif
388	enddo
389
390	! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
391
392	do k = 1, nlvl
393	if (plvl(k).lt.200000.) then
394	if (.not. is_there(nint(plvl(k)),'RH').and. &
395	is_there(nint(plvl(k)),'DEPR')) then
396	call get_dims(nint(plvl(k)),'TT')
397	call compute_rh_depr(map%nx, map%ny, plvl(k))
398	endif
399	endif
400	enddo
401	!
402	! If upper-air RH is missing, see if we can interpolate from surrounding levels.
403	! This is a simple vertical interpolation, linear in pressure.
404	! Currently, this simply fills in one missing level between two present levels.
405	! May expand this in the future to fill in additional levels. May also expand
406	! this in the future to vertically interpolate other variables.
407	!
408
409	do k = 2, nlvl-1, 1
410	if (plvl(k-1) .lt. 200000.) then
411	if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
412	( is_there(nint(plvl(k-1)), 'RH')) .and.&
413	( is_there(nint(plvl(k+1)), 'RH')) ) then
414	call get_dims(nint(plvl(k+1)), 'RH')
415	call vntrp(plvl, maxlvl, k, "RH ", map%nx, map%ny)
416	endif
417	endif
418	enddo
419
420	!
421	! Check to see if we need to fill RH above 300 mb to 10%:
422	!
423	if (is_there(30000, 'RH')) then
424	call get_dims(30000, 'RH')
425	allocate(scr2d(map%nx,map%ny))
426	scr2d = 10.
427
428	do k = 1, nlvl
429	if (plvl(k).lt.30000.) then
430	if (.not. is_there(nint(plvl(k)), 'RH')) then
431	call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
432	endif
433	endif
434	enddo
435	deallocate(scr2d)
436	endif
437	!
438	! If surface RH is missing, see if we can compute RH from Specific Humidity
439	! or Dewpoint or Dewpoint depression:
440	!
441	if (.not. is_there (200100, 'RH')) then
442	if (is_there(200100, 'TT').and. &
443	is_there(200100, 'PSFC' ) .and. &
444	is_there(200100, 'SPECHUMD')) then
445	call get_dims(200100, 'TT')
446	call compute_rh_spechumd(map%nx, map%ny)
447	call mprintf(.true.,DEBUG, &
448	"RRPR: SURFACE RH is computed")
449	elseif (is_there(200100, 'TT' ).and. &
450	is_there(200100, 'DEWPT')) then
451	call get_dims(200100, 'TT')
452	call compute_rh_dewpt(map%nx, map%ny)
453	elseif (is_there(200100, 'TT').and. &
454	is_there(200100, 'DEPR')) then
455	call get_dims(200100, 'TT')
456	call compute_rh_depr(map%nx, map%ny, 200100.)
457	endif
458	endif
459
460	! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
461
462	if (is_there(200100, 'SEAICE')) then
463	call get_dims(200100, 'SEAICE')
464	call make_zero_or_one(map%nx, map%ny)
465	endif
466
467	call mprintf(.true.,INFORM, &
468	"RRPR: hdate = %s ",s1=hdate)
469	call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
470	call clear_storage
471	exit FILELOOP
472	enddo FILELOOP
473	enddo TIMELOOP
474	end subroutine rrpr
475
476	subroutine make_zero_or_one(ix, jx)
477	! Make sure the SEAICE field is zero or one.
478	use storage_module
479	implicit none
480	integer :: ix, jx
481	real, dimension(ix,jx) :: seaice
482
483	call get_storage(200100, 'SEAICE',seaice, ix, jx)
484	where(seaice > 0.5)
485	seaice = 1.0
486	elsewhere
487	seaice = 0.0
488	end where
489	call put_storage(200100, 'SEAICE',seaice, ix, jx)
490	end subroutine make_zero_or_one
491
492
493	subroutine compute_spechumd_qvapor(ix, jx, plvl)
494	! Compute specific humidity from water vapor mixing ratio.
495	use storage_module
496	implicit none
497	integer :: ix, jx
498	real :: plvl
499	real :: lat1, lon1, dx, dy
500	real, dimension(ix,jx) :: QVAPOR, SPECHUMD
501
502	real startlat, startlon, deltalat, deltalon
503
504	call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
505
506	SPECHUMD = QVAPOR/(1.+QVAPOR)
507
508	call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
509	if(nint(plvl).eq.1) then
510	call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
511	endif
512
513	end subroutine compute_spechumd_qvapor
514
515	subroutine compute_t_vptmp(ix, jx, plvl)
516	! Compute relative humidity from specific humidity.
517	use storage_module
518	implicit none
519	integer :: ix, jx
520	real :: plvl
521	real :: lat1, lon1, dx, dy
522	real, dimension(ix,jx) :: T, VPTMP, P, Q
523
524	real, parameter :: rovcp=0.28571
525
526	real startlat, startlon, deltalat, deltalon
527
528	call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx)
529	IF (nint(plvl) .LT. 200) THEN
530	call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
531	ELSE
532	p = plvl
533	ENDIF
534	call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
535
536	t=vptmp * (p1.e-5)rovcp (1./(1.+0.6078*Q))
537
538	call put_storage(nint(plvl), 'TT', t, ix, jx)
539	if(nint(plvl).eq.1) then
540	call put_storage(200100, 'PSFC', p, ix, jx)
541	endif
542
543	end subroutine compute_t_vptmp
544
545
546	subroutine compute_rh_spechumd(ix, jx)
547	! Compute relative humidity from specific humidity.
548	use storage_module
549	implicit none
550	integer :: ix, jx
551	real :: lat1, lon1, dx, dy
552	real, dimension(ix,jx) :: T, P, RH, Q
553
554	real, parameter :: svp1=611.2
555	real, parameter :: svp2=17.67
556	real, parameter :: svp3=29.65
557	real, parameter :: svpt0=273.15
558	real, parameter :: eps = 0.622
559
560	real startlat, startlon, deltalat, deltalon
561
562	call get_storage(200100, 'TT', T, ix, jx)
563	call get_storage(200100, 'PSFC', P, ix, jx)
564	call get_storage(200100, 'SPECHUMD', Q, ix, jx)
565
566	rh = 1.E2 * (pq/(q(1.-eps) + eps))/(svp1exp(svp2(t-svpt0)/(T-svp3)))
567
568	call put_storage(200100, 'RH', rh, ix, jx)
569
570	end subroutine compute_rh_spechumd
571
572	subroutine compute_rh_spechumd_upa(ix, jx, plvl)
573	! Compute relative humidity from specific humidity.
574	use storage_module
575	implicit none
576	integer :: ix, jx
577	real :: plvl
578	real :: lat1, lon1, dx, dy
579	real, dimension(ix,jx) :: T, P, RH, Q
580
581	real, parameter :: svp1=611.2
582	real, parameter :: svp2=17.67
583	real, parameter :: svp3=29.65
584	real, parameter :: svpt0=273.15
585	real, parameter :: eps = 0.622
586
587	real startlat, startlon, deltalat, deltalon
588
589	IF ( nint(plvl).LT. 200) THEN
590	call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
591	ELSE
592	P = plvl
593	ENDIF
594	call get_storage(nint(plvl), 'TT', T, ix, jx)
595	call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
596
597	rh = 1.E2 * (pq/(q(1.-eps) + eps))/(svp1exp(svp2(t-svpt0)/(T-svp3)))
598
599	call put_storage(nint(plvl), 'RH', rh, ix, jx)
600
601	end subroutine compute_rh_spechumd_upa
602
603	subroutine compute_rh_vapp_upa(ix, jx, plvl)
604	! Compute relative humidity from vapor pressure.
605	! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
606	use storage_module
607	implicit none
608	integer :: ix, jx
609	real :: plvl
610	real :: lat1, lon1, dx, dy
611	real, dimension(ix,jx) :: P, ES
612	real, pointer, dimension(:,:) :: T, E, RH
613
614	real, parameter :: svp1=611.2
615	real, parameter :: svp2=17.67
616	real, parameter :: svp3=29.65
617	real, parameter :: svpt0=273.15
618	real, parameter :: eps = 0.622
619
620	real :: startlat, startlon, deltalat, deltalon
621
622	allocate(RH(ix,jx))
623
624	P = plvl
625	call refr_storage(nint(plvl), 'TT', T, ix, jx)
626	call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
627
628	ES=svp1exp(svp2(T-svpt0)/(T-svp3))
629	rh=min(1.E2(P-ES)E/((P-E)*ES), 1.E2)
630
631	call refw_storage(nint(plvl), 'RH', rh, ix, jx)
632
633	nullify(T,E)
634
635	end subroutine compute_rh_vapp_upa
636
637	subroutine compute_rh_depr(ix, jx, plvl)
638	! Compute relative humidity from Dewpoint Depression
639	use storage_module
640	implicit none
641	integer :: ix, jx
642	real :: plvl
643	real :: lat1, lon1, dx, dy
644	real, dimension(ix,jx) :: t, depr, rh
645
646	real, parameter :: Xlv = 2.5e6
647	real, parameter :: Rv = 461.5
648
649	integer :: i, j
650
651	call get_storage(nint(plvl), 'TT', T, ix, jx)
652	call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
653
654	where(DEPR < 100.)
655	rh = exp(Xlv/Rv(1./T - 1./(T-depr))) 1.E2
656	elsewhere
657	rh = 0.0
658	endwhere
659
660	call put_storage(nint(plvl),'RH ', rh, ix, jx)
661
662	end subroutine compute_rh_depr
663
664	subroutine compute_rh_dewpt(ix,jx)
665	! Compute relative humidity from Dewpoint
666	use storage_module
667	implicit none
668	integer :: ix, jx
669	real :: lat1, lon1, dx, dy
670	real, dimension(ix,jx) :: t, dp, rh
671
672	real, parameter :: Xlv = 2.5e6
673	real, parameter :: Rv = 461.5
674
675	call get_storage(200100, 'TT ', T, ix, jx)
676	call get_storage(200100, 'DEWPT ', DP, ix, jx)
677
678	rh = exp(Xlv/Rv(1./T - 1./dp)) 1.E2
679
680	call put_storage(200100,'RH ', rh, ix, jx)
681
682	end subroutine compute_rh_dewpt
683
684	subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
685	use storage_module
686	implicit none
687	integer :: ix, jx, k, maxlvl
688	real, dimension(maxlvl) :: plvl
689	character(len=8) :: name
690	real, dimension(ix,jx) :: a, b, c
691	real :: frc
692
693	write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
694
695	call get_storage(nint(plvl(k-1)), name, a, ix, jx)
696	call get_storage(nint(plvl(k+1)), name, c, ix, jx)
697
698	frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
699
700	b = (1.-frc)a + frcc
701	!KWM b = 0.5 * (a + c)
702	call put_storage(nint(plvl(k)), name, b, ix, jx)
703
704	end subroutine vntrp
705

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format