| blob | e286d811241f4e106a4b7ed70bbb791a3f350335 |
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 ! !
8 use filelist
9 use gridinfo
10 use storage_module
11 use table
12 use module_debug
13 use misc_definitions_module
14 use stringutil
16 implicit none
18 !------------------------------------------------------------------------------
19 ! Arguments:
21 ! HSTART: Starting date of times to process
22 character (LEN=19) :: hstart
24 ! NTIMES: Number of time periods to process
25 integer :: ntimes
27 ! INTERVAL: Time inteval (seconds) of time periods to process.
28 integer :: interval
30 ! NLVL: The number of levels in the stored data.
31 integer :: nlvl
33 ! MAXLVL: The parameterized maximum number of levels to allow.
34 integer :: maxlvl
36 ! PLVL: Array of pressure levels (Pa) in the dataset
37 real , dimension(maxlvl) :: plvl
39 ! DEBUG_LEVEL: Integer level of debug printing (from namelist)
40 integer :: debug_level
42 !------------------------------------------------------------------------------
44 character (LEN=25) :: units
45 character (LEN=46) :: Desc
46 real, allocatable, dimension(:,:) :: scr2d, tmp2d
47 real, pointer, dimension(:,:) :: ptr2d
49 integer :: k, kk, mm, n, ierr, ifv
50 integer :: iunit=13
52 character(LEN=19) :: hdate, hend
53 character(LEN=24) :: hdate_output
54 character(LEN=3) :: out_format
55 character(LEN=MAX_FILENAME_LEN) :: prefix
56 real :: xfcst, level
57 character(LEN=9) :: field
59 integer :: ntime, idts
61 ! DATELEN: length of date strings to use for our output file names.
62 integer :: datelen
64 ! Decide the length of date strings to use for output file names.
65 ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
67 if (mod(interval,3600) == 0) then
68 datelen = 13
69 else if (mod(interval, 60) == 0) then
70 datelen = 16
71 else
72 datelen = 19
73 endif
75 if ( debug_level .gt. 100 ) then
76 call mprintf(.true.,DEBUG,"Begin rrpr")
77 call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
78 do n = 1, nfiles
79 call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
80 enddo
81 endif
83 ! Compute the ending time:
85 call geth_newdate(hend, hstart, interval*ntimes)
87 call clear_storage
89 ! We want to do something for each of the requested times:
90 TIMELOOP : do ntime = 1, ntimes
91 idts = (ntime-1) * interval
92 call geth_newdate(hdate, hstart, idts)
93 call mprintf(.true.,DEBUG, &
94 "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
96 ! Loop over the output file dates, and do stuff if the file date matches
97 ! the requested time we are working on now.
99 FILELOOP : do n = 1, nfiles
100 if ( debug_level .gt. 100 ) then
101 call mprintf(.true.,DEBUG, &
102 "hstart = %s , hend = %s",s1=hstart,s2=hend)
103 call mprintf(.true.,DEBUG, &
104 "filedates(n) = %s",s1=filedates(n))
105 call mprintf(.true.,DEBUG, &
106 "filedates(n) = %s",s1=filedates(n)(1:datelen))
107 end if
108 if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
109 if (debug_level .gt. 50 ) then
110 call mprintf(.true.,INFORM, &
111 "RRPR Processing : %s",s1=filedates(n)(1:datelen))
112 endif
113 open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
114 form='unformatted',status='old')
116 ! Read the file:
118 rdloop: do
119 read (iunit, iostat=ierr) ifv
120 if (ierr.ne.0) exit rdloop
121 if ( ifv .eq. 5) then ! WPS
122 read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
123 level, map%nx, map%ny, map%igrid
124 hdate = hdate_output(1:19)
125 select case (map%igrid)
126 case (0, 4)
127 read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
128 case (3)
129 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
130 map%truelat1, map%truelat2, map%r_earth
131 case (5)
132 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
133 map%truelat1, map%r_earth
134 case (1)
135 read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
136 map%truelat1, map%r_earth
137 case default
138 call mprintf(.true.,ERROR, &
139 "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
140 end select
141 read (iunit) map%grid_wind
143 else if ( ifv .eq. 4 ) then ! SI
144 read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
145 map%nx, map%ny, map%igrid
146 hdate = hdate_output(1:19)
147 select case (map%igrid)
148 case (0, 4)
149 read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
150 case (3)
151 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
152 map%lov, map%truelat1, map%truelat2
153 case (5)
154 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
155 map%lov, map%truelat1
156 case default
157 call mprintf(.true.,ERROR, &
158 "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
159 end select
161 else if ( ifv .eq. 3 ) then ! MM5
162 read(iunit) hdate_output, xfcst, field, units, desc, level,&
163 map%nx, map%ny, map%igrid
164 hdate = hdate_output(1:19)
165 select case (map%igrid)
166 case (3) ! lamcon
167 read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
168 map%truelat1, map%truelat2
169 case (5) ! Polar Stereographic
170 read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
171 map%truelat1
172 case (0, 4) ! lat/lon
173 read (iunit) map%lat1, map%lon1, map%dy, map%dx
174 case (1) ! Mercator
175 read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
176 case default
177 call mprintf(.true.,ERROR, &
178 "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
179 end select
180 else
181 call mprintf(.true.,ERROR, &
182 "unknown out_format, ifv = %i",i1=ifv)
183 endif
185 allocate(ptr2d(map%nx,map%ny))
186 read (iunit) ptr2d
187 call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
188 nullify (ptr2d)
189 enddo rdloop
191 write (0,*) 'Name of source model =>',map%source
192 !
193 ! We have reached the end of file, so time to close it.
194 !
195 close(iunit)
196 if (debug_level .gt. 100 ) call print_storage
197 !
198 ! By now the file has been read completely. Now, see if we need to fill in
199 ! missing fields:
200 !
202 ! Retrieve the number of levels in storage:
203 !
204 call get_plvls(plvl, maxlvl, nlvl)
205 !
206 ! Fill the surface level (code 200100) from higher 200100s, as necessary
207 !
208 do k = 1, nlvl
209 if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
210 ! We found a level between 200100 and 200200, now find the field
211 ! corresponding to that level.
212 MLOOP : do mm = 1, maxvar
213 if (is_there(nint(plvl(k)), namvar(mm))) then
214 INLOOP : do kk = 200101, nint(plvl(k))
215 if (is_there(kk, namvar(mm))) then
216 if ( debug_level .gt. 100 ) then
217 call mprintf(.true.,DEBUG, &
218 "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
219 end if
220 call get_dims(kk, namvar(mm))
221 allocate(scr2d(map%nx,map%ny))
222 call get_storage &
223 (kk, namvar(mm), scr2d, map%nx, map%ny)
224 call put_storage &
225 (200100,namvar(mm), scr2d,map%nx,map%ny)
226 deallocate(scr2d)
227 EXIT INLOOP
228 endif
229 enddo INLOOP
230 endif
231 enddo MLOOP
232 endif
233 enddo
235 !
236 ! If upper-air U is missing, see if we can interpolate from surrounding levels.
237 ! This is a simple vertical interpolation, linear in pressure.
238 ! Currently, this simply fills in one missing level between two present levels.
239 !
241 do k = 2, nlvl-1, 1
242 if (plvl(k-1) .lt. 200000.) then
243 if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
244 ( is_there(nint(plvl(k-1)), 'UU')) .and.&
245 ( is_there(nint(plvl(k+1)), 'UU')) ) then
246 call get_dims(nint(plvl(k+1)), 'UU')
247 call vntrp(plvl, maxlvl, k, "UU ", map%nx, map%ny)
248 endif
249 endif
250 enddo
252 !
253 ! If upper-air V is missing, see if we can interpolate from surrounding levels.
254 ! This is a simple vertical interpolation, linear in pressure.
255 ! Currently, this simply fills in one missing level between two present levels.
256 !
258 do k = 2, nlvl-1, 1
259 if (plvl(k-1) .lt. 200000.) then
260 if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
261 ( is_there(nint(plvl(k-1)), 'VV')) .and.&
262 ( is_there(nint(plvl(k+1)), 'VV')) ) then
263 call get_dims(nint(plvl(k+1)), 'VV')
264 call vntrp(plvl, maxlvl, k, "VV ", map%nx, map%ny)
265 endif
266 endif
267 enddo
269 !
270 ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
271 !--- Tanya's change for initializing WRF with RUC
273 do k = 1, nlvl
274 if (plvl(k).lt.200000.) then
275 if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
276 is_there(nint(plvl(k)), 'QV')) then
277 call get_dims(nint(plvl(k)), 'QV')
278 call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
279 endif
280 endif
281 enddo
283 !--- Tanya's change for initializing WRF with RUC
284 ! This allows for the ingestion for RUC isentropic data
285 !
286 do k = 1, nlvl
287 if (plvl(k).lt.200000.) then
288 if (.not. is_there(nint(plvl(k)), 'TT').and. &
289 is_there(nint(plvl(k)), 'VPTMP').and. &
290 is_there(nint(plvl(k)), 'SPECHUMD')) then
291 call get_dims(nint(plvl(k)), 'VPTMP')
292 call compute_t_vptmp(map%nx, map%ny, plvl(k))
293 endif
294 endif
295 enddo
296 !!!
297 !
298 ! If upper-air T is missing, see if we can interpolate from surrounding levels.
299 ! This is a simple vertical interpolation, linear in pressure.
300 ! Currently, this simply fills in one missing level between two present levels.
301 !
303 do k = 2, nlvl-1, 1
304 if (plvl(k-1) .lt. 200000.) then
305 if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
306 ( is_there(nint(plvl(k-1)), 'TT')) .and.&
307 ( is_there(nint(plvl(k+1)), 'TT')) ) then
308 call get_dims(nint(plvl(k+1)), 'TT')
309 call vntrp(plvl, maxlvl, k, "TT ", map%nx, map%ny)
310 endif
311 endif
312 enddo
314 !
315 ! Check to see if we need to fill HGT from GEOPT.
316 !
317 do k = 1, nlvl
318 if (plvl(k).lt.200000.) then
319 if (.not. is_there(nint(plvl(k)), 'HGT').and. &
320 is_there(nint(plvl(k)), 'GEOPT')) then
321 call get_dims(nint(plvl(k)), 'GEOPT')
322 allocate(scr2d(map%nx,map%ny))
323 call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
324 scr2d = scr2d / 9.81
325 call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny)
326 deallocate(scr2d)
327 endif
328 endif
329 enddo
332 ! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
334 do k = 1, nlvl
335 if (plvl(k).lt.200000.) then
336 if (.not. is_there(nint(plvl(k)), 'RH') .and. &
337 is_there(nint(plvl(k)), 'TT') .and. &
338 is_there(nint(plvl(k)), 'SPECHUMD')) then
339 call get_dims(nint(plvl(k)), 'TT')
340 call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
341 endif
342 endif
343 enddo
345 ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
346 ! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
348 do k = 1, nlvl
349 if (plvl(k).lt.200000.) then
350 if (.not. is_there(nint(plvl(k)),'RH').and. &
351 is_there(nint(plvl(k)), 'TT') .and. &
352 is_there(nint(plvl(k)),'VAPP')) then
353 call get_dims(nint(plvl(k)),'TT')
354 call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
355 endif
356 endif
357 enddo
359 ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
361 do k = 1, nlvl
362 if (plvl(k).lt.200000.) then
363 if (.not. is_there(nint(plvl(k)),'RH').and. &
364 is_there(nint(plvl(k)), 'TT') .and. &
365 is_there(nint(plvl(k)),'DEPR')) then
366 call get_dims(nint(plvl(k)),'TT')
367 call compute_rh_depr(map%nx, map%ny, plvl(k))
368 endif
369 endif
370 enddo
371 !
372 ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
373 ! This is a simple vertical interpolation, linear in pressure.
374 ! Currently, this simply fills in one missing level between two present levels.
375 ! May expand this in the future to fill in additional levels. May also expand
376 ! this in the future to vertically interpolate other variables.
377 !
379 do k = 2, nlvl-1, 1
380 if (plvl(k-1) .lt. 200000.) then
381 if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
382 ( is_there(nint(plvl(k-1)), 'RH')) .and.&
383 ( is_there(nint(plvl(k+1)), 'RH')) ) then
384 call get_dims(nint(plvl(k+1)), 'RH')
385 call vntrp(plvl, maxlvl, k, "RH ", map%nx, map%ny)
386 endif
387 endif
388 enddo
390 ! Repair GFS and ECMWF pressure-level RH
391 if (index(map%source,'NCEP GFS') .ne. 0 .or. &
392 index(map%source,'ECMWF') .ne. 0 ) then
393 call mprintf(.true.,DEBUG, &
394 "RRPR: Adjusting GFS/ECMWF RH values ")
395 do k = 1, nlvl
396 if ( is_there(nint(plvl(k)),'RH') .and. &
397 is_there(nint(plvl(k)),'TT') ) then
398 call fix_gfs_rh (map%nx, map%ny, plvl(k))
399 endif
400 enddo
401 endif
402 !
403 ! Check to see if we need to fill RH above 300 mb:
404 !
405 if (is_there(30000, 'RH')) then
406 call get_dims(30000, 'RH')
407 allocate(scr2d(map%nx,map%ny))
409 do k = 1, nlvl
410 ! Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
411 ! The stratospheric RH will be adjusted further in real.
412 if (plvl(k).le.7000.) then
413 scr2d = 0.
414 else if (plvl(k).lt.30000.) then
415 scr2d = 5.
416 endif
417 if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
418 ! levels higher than .1 mb are special - do not fill
419 if (.not. is_there(nint(plvl(k)), 'RH')) then
420 call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
421 call mprintf(.true.,DEBUG, &
422 "RRPR: RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
423 endif
424 endif
425 enddo
426 deallocate(scr2d)
427 endif
428 !
429 ! If surface RH is missing, see if we can compute RH from Specific Humidity
430 ! or Dewpoint or Dewpoint depression:
431 !
432 if (.not. is_there (200100, 'RH')) then
433 if (is_there(200100, 'TT').and. &
434 is_there(200100, 'PSFC' ) .and. &
435 is_there(200100, 'SPECHUMD')) then
436 call get_dims(200100, 'TT')
437 call compute_rh_spechumd(map%nx, map%ny)
438 call mprintf(.true.,DEBUG, &
439 "RRPR: SURFACE RH is computed")
440 elseif (is_there(200100, 'TT' ).and. &
441 is_there(200100, 'DEWPT')) then
442 call get_dims(200100, 'TT')
443 call compute_rh_dewpt(map%nx, map%ny)
444 elseif (is_there(200100, 'TT').and. &
445 is_there(200100, 'DEPR')) then
446 call get_dims(200100, 'TT')
447 call compute_rh_depr(map%nx, map%ny, 200100.)
448 endif
449 endif
451 !
452 ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
453 ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
454 !
455 if (.not. is_there(200100, 'SNOW') .and. &
456 is_there(200100, 'SNOWRUC')) then
457 call get_dims(200100, 'SNOWRUC')
458 allocate(scr2d(map%nx,map%ny))
459 call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
460 scr2d = scr2d * 1000.
461 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
462 deallocate(scr2d)
463 endif
465 ! compute snow water equivalent (SNOW) for NCEP RUC models
466 ! As of Sept. 14 2011
467 if ( index(map%source,'NCEP RUC Model') .ne. 0) then
468 if (is_there(200100, 'SNOWH') .and. .not. is_there(200100, 'SNOW')) then
469 call get_dims(200100, 'SNOWH')
470 allocate(scr2d(map%nx,map%ny))
471 call get_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
472 call mprintf(.true.,DEBUG, &
473 "RRPR: Computing SNOWH from SNOW")
474 if (is_there(200100, 'RHOSN')) then ! If we have snow density, use it to compute snowh
475 call get_dims(200100, 'RHOSN')
476 allocate(tmp2d(map%nx,map%ny))
477 call get_storage(200100, 'RHOSN', tmp2d, map%nx, map%ny)
478 scr2d = scr2d * tmp2d
479 deallocate(tmp2d)
480 else
481 scr2d = scr2d * 200.0 ! Assume 200:1 ratio
482 endif
483 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
484 deallocate(scr2d)
485 endif
486 endif
488 ! Add residual soil moisture to SOILM* if initialized from the GSD RUC model or from NCEP RUC
489 if (index(map%source,'NOAA GSD') .ne. 0 .or. &
490 index(map%source,'NCEP RUC Model') .ne. 0) then
491 if ( .not. is_there(200100, 'SOILM000') .and.&
492 is_there(200100, 'SM000ruc') ) then
493 call get_dims(200100, 'SM000ruc')
494 print *,'Adjust RUC soil moisture'
495 call mprintf(.true.,DEBUG, &
496 "RRPR: Adjusting RUC soil moisture ")
497 call fix_ruc_soilm (map%nx, map%ny)
498 endif
499 endif
501 !
502 ! Check to see if we need to fill SOILHGT from SOILGEO.
503 ! (From Wei Wang, 2007 June 21)
504 !
505 if (.not. is_there(200100, 'SOILHGT') .and. &
506 is_there(200100, 'SOILGEO')) then
507 call get_dims(200100, 'SOILGEO')
508 allocate(scr2d(map%nx,map%ny))
509 call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
510 scr2d = scr2d / 9.81
511 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
512 deallocate(scr2d)
513 endif
515 ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
516 if (.not. is_there(200100, 'SOILHGT') .and. &
517 is_there(1, 'SOILGEO')) then
518 call get_dims(1, 'SOILGEO')
519 allocate(scr2d(map%nx,map%ny))
520 call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
521 scr2d = scr2d / 9.81
522 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
523 deallocate(scr2d)
524 endif
526 ! For NCEP RR (using the same ID as for RUC) native-level input,
527 ! may need to move PSFC from level 1 to 2001.
528 ! From TGS 8 Sept. 2011
529 if ( index(map%source,'NCEP RUC Model') .ne. 0) then
530 if (.not. is_there(200100, 'PSFC') .and. &
531 is_there(1, 'PRESSURE')) then
532 print *,'Process PSFC for NCEP RR'
533 call get_dims(1, 'PRESSURE')
534 allocate(scr2d(map%nx,map%ny))
535 call get_storage(1, 'PRESSURE', scr2d, map%nx, map%ny)
536 call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
537 deallocate(scr2d)
538 endif
539 endif
541 ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
542 if ( index(map%source,'ECMWF') .ne. 0) then
543 if (.not. is_there(200100, 'PSFC') .and. &
544 is_there(1, 'PSFCH')) then
545 call get_dims(1, 'PSFCH')
546 allocate(scr2d(map%nx,map%ny))
547 call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
548 call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
549 deallocate(scr2d)
550 endif
551 endif
553 ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2
554 !
555 if (is_there(200100, 'SNOW_EC')) then
556 call get_dims(200100, 'SNOW_EC')
557 allocate(scr2d(map%nx,map%ny))
558 call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
559 scr2d = scr2d * 1000.
560 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
561 deallocate(scr2d)
562 endif
564 ! Convert the ECMWF LANDSEA mask from a fraction to a flag
566 if ( index(map%source,'ECMWF') .ne. 0) then
567 if (is_there(200100, 'LANDSEA')) then
568 call get_dims(200100, 'LANDSEA')
569 call make_zero_or_one(map%nx, map%ny, 'LANDSEA')
570 endif
571 endif
573 ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF.
574 ! The GFS-based reanalyses values should be OK as is.
575 if ((index(map%source,'NCEP GFS') .ne. 0 .or. &
576 index(map%source,'NCEP GEFS') .ne. 0) .and. &
577 is_there(200100, 'SNOW')) then
578 call mprintf(.true.,DEBUG, &
579 "RRPR: Recomputing SNOW for NCEP GFS")
580 call get_dims(200100, 'SNOW')
581 allocate(scr2d(map%nx,map%ny))
582 call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
583 scr2d = scr2d * 2.
584 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
585 deallocate(scr2d)
586 endif
588 ! compute physical snow depth (SNOWH) for various models
589 ! As of March 2011, this is done here instead of real because we have model
590 ! source information.
591 if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then
592 call get_dims(200100, 'SNOW')
593 allocate(scr2d(map%nx,map%ny))
594 call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
595 call mprintf(.true.,DEBUG, &
596 "RRPR: Computing SNOWH from SNOW")
597 if ( index(map%source,'NCEP ') .ne. 0) then
598 scr2d = scr2d * 0.005 ! Assume 200:1 ratio as used at NCEP and in NOAH
599 else if (index(map%source,'ECMWF') .ne. 0) then
600 if (is_there(200100, 'SNOW_DEN')) then ! If we have snow density, use it to compute snowh
601 call get_dims(200100, 'SNOW_DEN')
602 allocate(tmp2d(map%nx,map%ny))
603 call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny)
604 scr2d = scr2d / tmp2d
605 deallocate(tmp2d)
606 else
607 scr2d = scr2d * 0.004 ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio).
608 endif
609 else ! Other models
610 scr2d = scr2d * 0.005 ! Use real's default method (200:1)
611 endif
612 call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
613 deallocate(scr2d)
614 endif
616 ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted
617 ! to the appropriate values in real depending on whether or not the polar mods are used.
619 !! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
621 ! if (is_there(200100, 'SEAICE')) then
622 ! call get_dims(200100, 'SEAICE')
623 ! call make_zero_or_one(map%nx, map%ny, 'SEAICE')
624 ! endif
626 ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
627 ! Field | GRIB In | Out
628 ! -------------------------
629 ! water | 0 | 0
630 ! land | -1 | 1
631 ! ice | 1 | 0
633 if (is_there(200100, 'ICEMASK')) then
634 call get_dims(200100, 'ICEMASK')
635 call re_flag_ice_mask(map%nx, map%ny)
636 endif
638 ! If we have an ICEFRAC field, convert from % to fraction
639 if (is_there(200100, 'ICEFRAC')) then
640 call get_dims(200100, 'ICEFRAC')
641 allocate(scr2d(map%nx,map%ny))
642 call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
643 scr2d = scr2d / 100.
644 call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
645 deallocate(scr2d)
646 endif
649 call mprintf(.true.,INFORM, &
650 "RRPR: hdate = %s ",s1=hdate)
651 call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
652 call clear_storage
653 exit FILELOOP
654 enddo FILELOOP
655 enddo TIMELOOP
656 end subroutine rrpr
658 subroutine make_zero_or_one(ix, jx, infield)
659 ! Make sure the input field (SEAICE or LANDSEA) is zero or one.
660 use storage_module
661 implicit none
662 integer :: ix, jx
663 real, dimension(ix,jx) :: seaice
664 character(len=*) :: infield
666 call get_storage(200100, infield, seaice, ix, jx)
667 where(seaice > 0.5)
668 seaice = 1.0
669 elsewhere
670 seaice = 0.0
671 end where
672 call put_storage(200100, infield, seaice, ix, jx)
673 end subroutine make_zero_or_one
675 subroutine re_flag_ice_mask(ix, jx)
676 !
677 ! Change land points from -1 to 1
678 ! Change ice points from 1 to 0
679 ! Water points stay at 0
680 !
681 use storage_module
682 implicit none
683 integer :: ix, jx
684 real, dimension(ix,jx) :: iceflag
686 call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
687 where(iceflag > 0.5) ! Ice points, set to water value
688 iceflag = 0.0
689 end where
690 where(iceflag < -0.5) ! Land points
691 iceflag = 1.0
692 end where
693 call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
694 end subroutine re_flag_ice_mask
696 subroutine compute_spechumd_qvapor(ix, jx, plvl)
697 ! Compute specific humidity from water vapor mixing ratio.
698 use storage_module
699 implicit none
700 integer :: ix, jx
701 real :: plvl
702 real, dimension(ix,jx) :: QVAPOR, SPECHUMD
704 call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
706 SPECHUMD = QVAPOR/(1.+QVAPOR)
708 call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
709 if(nint(plvl).eq.1) then
710 call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
711 endif
713 end subroutine compute_spechumd_qvapor
715 subroutine compute_t_vptmp(ix, jx, plvl)
716 ! Compute temperature from virtual potential temperature
717 use storage_module
718 implicit none
719 integer :: ix, jx
720 real :: plvl
721 real, dimension(ix,jx) :: T, VPTMP, P, Q
723 real, parameter :: rovcp=0.28571
725 call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx)
726 IF (nint(plvl) .LT. 200) THEN
727 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
728 ELSE
729 p = plvl
730 ENDIF
731 call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
733 t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
735 call put_storage(nint(plvl), 'TT', t, ix, jx)
736 if(nint(plvl).eq.1) then
737 call put_storage(200100, 'PSFC', p, ix, jx)
738 endif
740 end subroutine compute_t_vptmp
743 subroutine compute_rh_spechumd(ix, jx)
744 ! Compute relative humidity from specific humidity.
745 use storage_module
746 implicit none
747 integer :: ix, jx
748 real, dimension(ix,jx) :: T, P, RH, Q
750 real, parameter :: svp1=611.2
751 real, parameter :: svp2=17.67
752 real, parameter :: svp3=29.65
753 real, parameter :: svpt0=273.15
754 real, parameter :: eps = 0.622
756 call get_storage(200100, 'TT', T, ix, jx)
757 call get_storage(200100, 'PSFC', P, ix, jx)
758 call get_storage(200100, 'SPECHUMD', Q, ix, jx)
760 rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
762 call put_storage(200100, 'RH', rh, ix, jx)
764 end subroutine compute_rh_spechumd
766 subroutine compute_rh_spechumd_upa(ix, jx, plvl)
767 ! Compute relative humidity from specific humidity.
768 use storage_module
769 implicit none
770 integer :: ix, jx
771 real :: plvl
772 real, dimension(ix,jx) :: T, P, RH, Q
774 real, parameter :: svp1=611.2
775 real, parameter :: svp2=17.67
776 real, parameter :: svp3=29.65
777 real, parameter :: svpt0=273.15
778 real, parameter :: eps = 0.622
780 IF ( nint(plvl).LT. 200) THEN
781 if (is_there(nint(plvl), 'PRESSURE')) then
782 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
783 else
784 return ! if we don't have pressure on model levels, return
785 endif
786 ELSE
787 P = plvl
788 ENDIF
789 call get_storage(nint(plvl), 'TT', T, ix, jx)
790 call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
791 Q=MAX(1.E-10,Q)
793 rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
795 call put_storage(nint(plvl), 'RH', rh, ix, jx)
797 end subroutine compute_rh_spechumd_upa
799 subroutine compute_rh_vapp_upa(ix, jx, plvl)
800 ! Compute relative humidity from vapor pressure.
801 ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
802 use storage_module
803 implicit none
804 integer :: ix, jx
805 real :: plvl
806 real, dimension(ix,jx) :: P, ES
807 real, pointer, dimension(:,:) :: T, E, RH
809 real, parameter :: svp1=611.2
810 real, parameter :: svp2=17.67
811 real, parameter :: svp3=29.65
812 real, parameter :: svpt0=273.15
814 allocate(RH(ix,jx))
816 IF ( nint(plvl).LT. 200) THEN
817 if (is_there(nint(plvl), 'PRESSURE')) then
818 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
819 else
820 return ! if we don't have pressure on model levels, return
821 endif
822 ELSE
823 P = plvl
824 ENDIF
825 call refr_storage(nint(plvl), 'TT', T, ix, jx)
826 call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
828 ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
829 rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
831 call refw_storage(nint(plvl), 'RH', rh, ix, jx)
833 nullify(T,E)
835 end subroutine compute_rh_vapp_upa
837 subroutine compute_rh_depr(ix, jx, plvl)
838 ! Compute relative humidity from Dewpoint Depression
839 use storage_module
840 implicit none
841 integer :: ix, jx
842 real :: plvl
843 real, dimension(ix,jx) :: t, depr, rh
845 real, parameter :: Xlv = 2.5e6
846 real, parameter :: Rv = 461.5
848 integer :: i, j
850 call get_storage(nint(plvl), 'TT', T, ix, jx)
851 call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
853 where(DEPR < 100.)
854 rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
855 elsewhere
856 rh = 0.0
857 endwhere
859 call put_storage(nint(plvl),'RH ', rh, ix, jx)
861 end subroutine compute_rh_depr
863 subroutine compute_rh_dewpt(ix,jx)
864 ! Compute relative humidity from Dewpoint
865 use storage_module
866 implicit none
867 integer :: ix, jx
868 real, dimension(ix,jx) :: t, dp, rh
870 real, parameter :: Xlv = 2.5e6
871 real, parameter :: Rv = 461.5
873 call get_storage(200100, 'TT ', T, ix, jx)
874 call get_storage(200100, 'DEWPT ', DP, ix, jx)
876 rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
878 call put_storage(200100,'RH ', rh, ix, jx)
880 end subroutine compute_rh_dewpt
882 subroutine vntrp(plvl, maxlvl, k, name, ix, jx)
883 use storage_module
884 implicit none
885 integer :: ix, jx, k, maxlvl
886 real, dimension(maxlvl) :: plvl
887 character(len=8) :: name
888 real, dimension(ix,jx) :: a, b, c
889 real :: frc
891 write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k))
893 call get_storage(nint(plvl(k-1)), name, a, ix, jx)
894 call get_storage(nint(plvl(k+1)), name, c, ix, jx)
896 frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1))
898 b = (1.-frc)*a + frc*c
899 !KWM b = 0.5 * (a + c)
900 call put_storage(nint(plvl(k)), name, b, ix, jx)
902 end subroutine vntrp
904 subroutine fix_gfs_rh (ix, jx, plvl)
905 ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
906 use storage_module
907 implicit none
908 integer :: ix, jx, i, j
909 real :: plvl, eis, ews, r
910 real, allocatable, dimension(:,:) :: rh, tt
912 allocate(rh(ix,jx))
913 allocate(tt(ix,jx))
914 call get_storage(nint(plvl), 'RH', rh, ix, jx)
915 call get_storage(nint(plvl), 'TT', tt, ix, jx)
916 do j = 1, jx
917 do i = 1, ix
918 if ( tt(i,j) .le. 273.15 ) then
919 ! Murphy and Koop 2005 ice saturation vapor pressure.
920 ! eis and ews in hPA, tt is in K
921 eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
922 - (0.00728332 * tt(i,j)))
923 ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
924 ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
926 ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
927 if ( tt(i,j) .gt. 253.15 ) then
928 ! A linear approximation to the GFS blending region ( -20 > T < 0 )
929 r = ((273.15 - tt(i,j)) / 20.)
930 r = (r * eis) + ((1-r)*ews)
931 else
932 r = eis
933 endif
934 rh(i,j) = rh(i,j) * (r / ews)
935 endif
936 enddo
937 enddo
938 call put_storage(nint(plvl), 'RH', rh, ix, jx)
939 deallocate (rh)
940 deallocate (tt)
941 end subroutine fix_gfs_rh
944 subroutine fix_ruc_soilm (ix, jx)
945 ! This routine adds residual soil moisture if initialized fron RUC
946 use storage_module
947 implicit none
948 integer :: ix, jx, i, j
949 REAL , DIMENSION(100) :: lqmi
950 real, allocatable, dimension(:,:) :: soilm000, soilm005, soilm020, &
951 soilm040, soilm160, soilm300,soilcat
952 allocate(soilm000(ix,jx))
953 allocate(soilm005(ix,jx))
954 allocate(soilm020(ix,jx))
955 allocate(soilm040(ix,jx))
956 allocate(soilm160(ix,jx))
957 allocate(soilm300(ix,jx))
958 allocate(soilcat(ix,jx))
959 call get_storage(200100, 'SM000ruc', soilm000, ix, jx)
960 call get_storage(200100, 'SM005ruc', soilm005, ix, jx)
961 call get_storage(200100, 'SM020ruc', soilm020, ix, jx)
962 call get_storage(200100, 'SM040ruc', soilm040, ix, jx)
963 call get_storage(200100, 'SM160ruc', soilm160, ix, jx)
964 call get_storage(200100, 'SM300ruc', soilm300, ix, jx)
966 call get_storage(200100, 'SOILCAT', soilcat, ix, jx)
968 lqmi(1:16) = &
969 (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
970 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
971 0.004, 0.065 /)
973 do j = 1, jx
974 do i = 1, ix
976 SOILM000(i,j)=SOILM000(i,j) + lqmi(nint(soilcat(i,j)))
977 SOILM005(i,j)=SOILM005(i,j) + lqmi(nint(soilcat(i,j)))
978 SOILM020(i,j)=SOILM020(i,j) + lqmi(nint(soilcat(i,j)))
979 SOILM040(i,j)=SOILM040(i,j) + lqmi(nint(soilcat(i,j)))
980 SOILM160(i,j)=SOILM160(i,j) + lqmi(nint(soilcat(i,j)))
981 SOILM300(i,j)=SOILM300(i,j) + lqmi(nint(soilcat(i,j)))
982 enddo
983 enddo
984 call put_storage(200100, 'SOILM000', soilm000, ix, jx)
985 call put_storage(200100, 'SOILM005', soilm005, ix, jx)
986 call put_storage(200100, 'SOILM020', soilm020, ix, jx)
987 call put_storage(200100, 'SOILM040', soilm040, ix, jx)
988 call put_storage(200100, 'SOILM160', soilm160, ix, jx)
989 call put_storage(200100, 'SOILM300', soilm300, ix, jx)
991 print *,'fix_ruc_soilm is done!'
993 deallocate(soilm000)
994 deallocate(soilm005)
995 deallocate(soilm020)
996 deallocate(soilm040)
997 deallocate(soilm160)
998 deallocate(soilm300)
999 deallocate(soilcat)
1001 end subroutine fix_ruc_soilm