| blob | 524734338956ff283257254351aeb556542dcac5 |
1 subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, &
2 add_lvls, new_plvl, interp_type, &
3 debug_level, out_format, prefix)
4 ! !
5 ! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite" !
6 ! !
7 !*****************************************************************************!
8 ! !
10 use filelist
11 use gridinfo
12 use storage_module
13 use table
14 use module_debug
15 use misc_definitions_module
16 use stringutil
18 implicit none
20 !------------------------------------------------------------------------------
21 ! Arguments:
23 ! HSTART: Starting date of times to process
24 character (LEN=19) :: hstart
26 ! NTIMES: Number of time periods to process
27 integer :: ntimes
29 ! INTERVAL: Time inteval (seconds) of time periods to process.
30 integer :: interval
32 ! NLVL: The number of levels in the stored data.
33 integer :: nlvl
35 ! MAXLVL: The parameterized maximum number of levels to allow.
36 integer :: maxlvl
38 ! PLVL: Array of pressure levels (Pa) in the dataset
39 real , dimension(maxlvl) :: plvl
41 ! NEW_PLVL: Array of the additional pressure levels (Pa) to interpolate to
42 real , dimension(maxlvl) :: new_plvl
44 ! TLVL: Array combining pressure levels (Pa) in PLVL and NEW_PLVL
45 real , dimension(maxlvl) :: tlvl
47 ! ADD_LVLS: Should we add levels via interpolation?
48 logical :: add_lvls
50 ! INTERP_TYPE: vertical Interpolation type
51 ! (1=log in pressure, anything else=linear in pressure)
52 integer :: interp_type
55 ! DEBUG_LEVEL: Integer level of debug printing (from namelist)
56 integer :: debug_level
58 !------------------------------------------------------------------------------
60 character (LEN=25) :: units
61 character (LEN=46) :: Desc
62 real, allocatable, dimension(:,:) :: scr2d, tmp2d
63 real, pointer, dimension(:,:) :: ptr2d
65 integer :: k, kk, mm, n, ierr, ifv
66 integer :: itest, nn, nl, lvls, tvls
67 integer :: iunit=13
69 character(LEN=19) :: hdate, hend
70 character(LEN=24) :: hdate_output
71 character(LEN=3) :: out_format
72 character(LEN=MAX_FILENAME_LEN) :: prefix
73 real :: xfcst, level
74 character(LEN=9) :: field
76 integer :: ntime, idts
78 logical :: found_level
79 real, dimension(maxlvl) :: new_plvl_to_sort
80 integer :: largest_number_loc
81 integer :: new_plvl_counter
83 ! DATELEN: length of date strings to use for our output file names.
84 integer :: datelen
86 ! Decide the length of date strings to use for output file names.
87 ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
89 if (mod(interval,3600) == 0) then
90 datelen = 13
91 else if (mod(interval, 60) == 0) then
92 datelen = 16
93 else
94 datelen = 19
95 endif
97 if ( debug_level .gt. 100 ) then
98 call mprintf(.true.,DEBUG,"Begin rrpr")
99 call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes)
100 do n = 1, nfiles
101 call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n))
102 enddo
103 endif
105 ! Compute the ending time:
107 call geth_newdate(hend, hstart, interval*ntimes)
109 call clear_storage
111 ! We want to do something for each of the requested times:
112 TIMELOOP : do ntime = 1, ntimes
113 idts = (ntime-1) * interval
114 call geth_newdate(hdate, hstart, idts)
115 call mprintf(.true.,DEBUG, &
116 "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts)
118 ! Loop over the output file dates, and do stuff if the file date matches
119 ! the requested time we are working on now.
121 FILELOOP : do n = 1, nfiles
122 if ( debug_level .gt. 100 ) then
123 call mprintf(.true.,DEBUG, &
124 "hstart = %s , hend = %s",s1=hstart,s2=hend)
125 call mprintf(.true.,DEBUG, &
126 "filedates(n) = %s",s1=filedates(n))
127 call mprintf(.true.,DEBUG, &
128 "filedates(n) = %s",s1=filedates(n)(1:datelen))
129 end if
130 if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP
131 if (debug_level .gt. 50 ) then
132 call mprintf(.true.,INFORM, &
133 "RRPR Processing : %s",s1=filedates(n)(1:datelen))
134 endif
135 open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), &
136 form='unformatted',status='old')
138 ! Read the file:
140 rdloop: do
141 read (iunit, iostat=ierr) ifv
142 if (ierr.ne.0) exit rdloop
143 if ( ifv .eq. 5) then ! WPS
144 read (iunit) hdate_output, xfcst, map%source, field, units, Desc, &
145 level, 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, map%r_earth
150 case (3)
151 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
152 map%truelat1, map%truelat2, map%r_earth
153 case (5)
154 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, &
155 map%truelat1, map%r_earth
156 case (1)
157 read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
158 map%truelat1, map%r_earth
159 case (6)
160 read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
161 map%lat0,map%lon0, map%r_earth
162 case default
163 call mprintf(.true.,ERROR, &
164 "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid)
165 end select
166 read (iunit) map%grid_wind
168 else if ( ifv .eq. 4 ) then ! SI
169 read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, &
170 map%nx, map%ny, map%igrid
171 hdate = hdate_output(1:19)
172 select case (map%igrid)
173 case (0, 4)
174 read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
175 case (3)
176 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
177 map%lov, map%truelat1, map%truelat2
178 case (5)
179 read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
180 map%lov, map%truelat1
181 case default
182 call mprintf(.true.,ERROR, &
183 "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid)
184 end select
186 else if ( ifv .eq. 3 ) then ! MM5
187 read(iunit) hdate_output, xfcst, field, units, desc, level,&
188 map%nx, map%ny, map%igrid
189 hdate = hdate_output(1:19)
190 select case (map%igrid)
191 case (3) ! lamcon
192 read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
193 map%truelat1, map%truelat2
194 case (5) ! Polar Stereographic
195 read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
196 map%truelat1
197 case (0, 4) ! lat/lon
198 read (iunit) map%lat1, map%lon1, map%dy, map%dx
199 case (1) ! Mercator
200 read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
201 case default
202 call mprintf(.true.,ERROR, &
203 "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid)
204 end select
205 else
206 call mprintf(.true.,ERROR, &
207 "unknown out_format, ifv = %i",i1=ifv)
208 endif
210 allocate(ptr2d(map%nx,map%ny))
211 read (iunit) ptr2d
212 call refw_storage(nint(level), field, ptr2d, map%nx, map%ny)
213 nullify (ptr2d)
214 enddo rdloop
216 write (0,*) 'Name of source model =>',map%source
217 !
218 ! We have reached the end of file, so time to close it.
219 !
220 close(iunit)
221 if (debug_level .gt. 100 ) call print_storage
222 !
223 ! By now the file has been read completely. Now, see if we need to fill in
224 ! missing fields:
225 !
227 ! Retrieve the number of levels in storage:
228 !
229 call get_plvls(plvl, maxlvl, nlvl)
231 ! Merge list of pressure levels in data with requested pressure levels
232 if ( add_lvls ) then
233 ! The merging code expects the user-defined pressure levels to be in
234 ! order from highest pressure to lowest pressure.
235 ! Sort the user-defined pressure levels accordingly
236 new_plvl_to_sort = new_plvl
237 ! Set array containing pressure levels to add to the default value set in read_namelist.F
238 new_plvl = -99999
239 DO new_plvl_counter = 1,maxlvl
240 largest_number_loc = MAXLOC(new_plvl_to_sort, DIM=1)
241 new_plvl(new_plvl_counter)=new_plvl_to_sort(largest_number_loc)
242 new_plvl_to_sort(largest_number_loc)=-99999.
243 END DO
245 tvls = 1
246 lvls = 1
247 loop_nvls : do nn=1,maxlvl
248 loop_lvls : do nl=lvls,maxlvl
249 if ( tvls > maxlvl ) then
250 call mprintf(.true.,ERROR, "Adding user-defined pressure levels resulted in too &
251 &many total pressure levels. Please increase maxlvl in ungrib.F")
252 endif
253 if ( plvl(nn) > 0.0 .AND. plvl(nn) >= new_plvl(nl) ) then
254 tlvl(tvls) = plvl(nn)
255 tvls = tvls + 1
256 if ( plvl(nn) == new_plvl(nl) ) lvls = lvls + 1
257 exit loop_lvls
258 endif
259 if ( plvl(nn) > 0.0 .AND. plvl(nn) < new_plvl(nl) ) then
260 tlvl(tvls) = new_plvl(nl)
261 tvls = tvls + 1
262 lvls = lvls + 1
263 endif
264 if ( plvl(nn) < 0.0 ) exit loop_nvls
265 enddo loop_lvls
266 enddo loop_nvls
267 plvl = tlvl
268 nlvl = tvls - 1
269 end if
271 !
272 ! Fill the surface level (code 200100) from higher 200100s, as necessary
273 !
274 do k = 1, nlvl
275 if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
276 ! We found a level between 200100 and 200200, now find the field
277 ! corresponding to that level.
278 MLOOP : do mm = 1, maxvar
279 if (is_there(nint(plvl(k)), namvar(mm))) then
280 INLOOP : do kk = 200101, nint(plvl(k))
281 if (is_there(kk, namvar(mm))) then
282 if ( debug_level .gt. 100 ) then
283 call mprintf(.true.,DEBUG, &
284 "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk)
285 end if
286 call get_dims(kk, namvar(mm))
287 allocate(scr2d(map%nx,map%ny))
288 call get_storage &
289 (kk, namvar(mm), scr2d, map%nx, map%ny)
290 call put_storage &
291 (200100,namvar(mm), scr2d,map%nx,map%ny)
292 deallocate(scr2d)
293 EXIT INLOOP
294 endif
295 enddo INLOOP
296 endif
297 enddo MLOOP
298 endif
299 enddo
301 !
302 ! If upper-air U is missing, see if we can interpolate from surrounding levels.
303 ! This is a simple vertical interpolation, linear or log in pressure.
304 ! Currently, this simply fills in missing levels between two present levels.
305 !
307 do k = 2, nlvl-1, 1
308 if (plvl(k-1) .lt. 200000.) then
309 if ( (.not. is_there(nint(plvl(k)),'UU')) .and. &
310 ( is_there(nint(plvl(k-1)), 'UU')) ) then
311 found_level = .false.
312 uu_loop : do itest = k+1,nlvl,1
313 if ( ( is_there(nint(plvl(itest)), 'UU')) ) then
314 found_level = .true.
315 exit uu_loop
316 endif
317 enddo uu_loop
318 if( found_level ) then
319 call get_dims(nint(plvl(itest)), 'UU')
320 call vntrp(plvl, maxlvl, k, itest, interp_type, "UU ", map%nx, map%ny)
321 else
322 PRINT *,'WARNING: Could not interpolate UU to level k=',k,' p=',plvl(k),&
323 'because could not find any level above this level.'
324 endif
325 endif
326 endif
327 enddo
329 !
330 ! If upper-air V is missing, see if we can interpolate from surrounding levels.
331 ! This is a simple vertical interpolation, linear or log in pressure.
332 ! Currently, this simply fills in missing levels between two present levels.
333 !
336 do k = 2, nlvl-1, 1
337 if (plvl(k-1) .lt. 200000.) then
338 if ( (.not. is_there(nint(plvl(k)),'VV')) .and. &
339 ( is_there(nint(plvl(k-1)), 'VV')) ) then
340 found_level = .false.
341 VV_loop : do itest = k+1,nlvl,1
342 if ( ( is_there(nint(plvl(itest)), 'VV')) ) then
343 found_level = .true.
344 exit VV_loop
345 endif
346 enddo VV_loop
347 if( found_level ) then
348 call get_dims(nint(plvl(itest)), 'VV')
349 call vntrp(plvl, maxlvl, k, itest, interp_type, "VV ", map%nx, map%ny)
350 else
351 PRINT *,'WARNING: Could not interpolate VV to level k=',k,' p=',plvl(k),&
352 'because could not find any level above this level.'
353 endif
354 endif
355 endif
356 enddo
358 !
359 ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR:
360 !--- Tanya's change for initializing WRF with RUC
362 do k = 1, nlvl
363 if (plvl(k).lt.200000.) then
364 if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. &
365 is_there(nint(plvl(k)), 'QV')) then
366 call get_dims(nint(plvl(k)), 'QV')
367 call compute_spechumd_qvapor(map%nx, map%ny, plvl(k))
368 endif
369 endif
370 enddo
372 !--- Tanya's change for initializing WRF with RUC
373 ! This allows for the ingestion for RUC isentropic data
374 !
375 do k = 1, nlvl
376 if (plvl(k).lt.200000.) then
377 if (.not. is_there(nint(plvl(k)), 'TT').and. &
378 is_there(nint(plvl(k)), 'VPTMP').and. &
379 is_there(nint(plvl(k)), 'SPECHUMD')) then
380 call get_dims(nint(plvl(k)), 'VPTMP')
381 call compute_t_vptmp(map%nx, map%ny, plvl(k))
382 endif
383 endif
384 enddo
385 !!!
386 !
387 ! If upper-air T is missing, see if we can interpolate from surrounding levels.
388 ! This is a simple vertical interpolation, linear or log in pressure.
389 ! Currently, this simply fills in missing levels between two present levels.
390 !
392 do k = 2, nlvl-1, 1
393 if (plvl(k-1) .lt. 200000.) then
394 if ( (.not. is_there(nint(plvl(k)),'TT')) .and. &
395 ( is_there(nint(plvl(k-1)), 'TT')) ) then
396 found_level = .false.
397 TT_loop : do itest = k+1,nlvl,1
398 if ( ( is_there(nint(plvl(itest)), 'TT')) ) then
399 found_level = .true.
400 exit TT_loop
401 endif
402 enddo TT_loop
403 if( found_level ) then
404 call get_dims(nint(plvl(itest)), 'TT')
405 call vntrp(plvl, maxlvl, k, itest, interp_type, "TT ", map%nx, map%ny)
406 else
407 PRINT *,'WARNING: Could not interpolate TT to level k=',k,' p=',plvl(k),&
408 'because could not find any level above this level.'
409 endif
410 endif
411 endif
412 enddo
414 ! Vertically interpolate to fill in other moisture variables
415 ! It seems that ultimately this should be wrapped in a function and probably loop over
416 ! all variables that can be vertically interpolated
418 ! QC
419 do k = 2, nlvl-1, 1
420 if (plvl(k-1) .lt. 200000.) then
421 if ( (.not. is_there(nint(plvl(k)),'QC')) .and. &
422 ( is_there(nint(plvl(k-1)), 'QC')) ) then
423 found_level = .false.
424 QC_loop : do itest = k+1,nlvl,1
425 if ( ( is_there(nint(plvl(itest)), 'QC')) ) then
426 found_level = .true.
427 exit QC_loop
428 endif
429 enddo QC_loop
430 if( found_level ) then
431 call get_dims(nint(plvl(itest)), 'QC')
432 call vntrp(plvl, maxlvl, k, itest, interp_type, "QC ", map%nx, map%ny)
433 else
434 PRINT *,'WARNING: Could not interpolate QC to level k=',k,' p=',plvl(k),&
435 'because could not find any level above this level.'
436 endif
437 endif
438 endif
439 enddo
441 ! QR
442 do k = 2, nlvl-1, 1
443 if (plvl(k-1) .lt. 200000.) then
444 if ( (.not. is_there(nint(plvl(k)),'QR')) .and. &
445 ( is_there(nint(plvl(k-1)), 'QR')) ) then
446 found_level = .false.
447 QR_loop : do itest = k+1,nlvl,1
448 if ( ( is_there(nint(plvl(itest)), 'QR')) ) then
449 found_level = .true.
450 exit QR_loop
451 endif
452 enddo QR_loop
453 if( found_level ) then
454 call get_dims(nint(plvl(itest)), 'QR')
455 call vntrp(plvl, maxlvl, k, itest, interp_type, "QR ", map%nx, map%ny)
456 else
457 PRINT *,'WARNING: Could not interpolate QR to level k=',k,' p=',plvl(k),&
458 'because could not find any level above this level.'
459 endif
460 endif
461 endif
462 enddo
464 ! QS
465 do k = 2, nlvl-1, 1
466 if (plvl(k-1) .lt. 200000.) then
467 if ( (.not. is_there(nint(plvl(k)),'QS')) .and. &
468 ( is_there(nint(plvl(k-1)), 'QS')) ) then
469 found_level = .false.
470 QS_loop : do itest = k+1,nlvl,1
471 if ( ( is_there(nint(plvl(itest)), 'QS')) ) then
472 found_level = .true.
473 exit QS_loop
474 endif
475 enddo QS_loop
476 if( found_level ) then
477 call get_dims(nint(plvl(itest)), 'QS')
478 call vntrp(plvl, maxlvl, k, itest, interp_type, "QS ", map%nx, map%ny)
479 else
480 PRINT *,'WARNING: Could not interpolate QS to level k=',k,' p=',plvl(k),&
481 'because could not find any level above this level.'
482 endif
483 endif
484 endif
485 enddo
487 ! QG
488 do k = 2, nlvl-1, 1
489 if (plvl(k-1) .lt. 200000.) then
490 if ( (.not. is_there(nint(plvl(k)),'QG')) .and. &
491 ( is_there(nint(plvl(k-1)), 'QG')) ) then
492 found_level = .false.
493 QG_loop : do itest = k+1,nlvl,1
494 if ( ( is_there(nint(plvl(itest)), 'QG')) ) then
495 found_level = .true.
496 exit QG_loop
497 endif
498 enddo QG_loop
499 if( found_level ) then
500 call get_dims(nint(plvl(itest)), 'QG')
501 call vntrp(plvl, maxlvl, k, itest, interp_type, "QG ", map%nx, map%ny)
502 else
503 PRINT *,'WARNING: Could not interpolate QG to level k=',k,' p=',plvl(k),&
504 'because could not find any level above this level.'
505 endif
506 endif
507 endif
508 enddo
511 !
512 ! Check to see if we need to fill HGT from GEOPT.
513 !
514 ! First make sure no GEOPT is missing
516 do k = 2, nlvl-1, 1
517 if (plvl(k-1) .lt. 200000.) then
518 if ( (.not. is_there(nint(plvl(k)),'GEOPT')) .and. &
519 ( is_there(nint(plvl(k-1)), 'GEOPT')) ) then
520 found_level = .false.
521 gg_loop : do itest = k+1,nlvl,1
522 if ( ( is_there(nint(plvl(itest)), 'GEOPT')) ) then
523 found_level = .true.
524 exit gg_loop
525 endif
526 enddo gg_loop
527 if( found_level ) then
528 call get_dims(nint(plvl(itest)), 'GEOPT')
529 call vntrp(plvl, maxlvl, k, itest, interp_type, "GEOPT ", map%nx, map%ny)
530 else
531 PRINT *,'WARNING: Could not interpolate GEOPT to level k=',k,' p=',plvl(k),&
532 'because could not find any level above this level.'
533 endif
534 endif
535 endif
536 enddo
538 do k = 1, nlvl
539 if (plvl(k).lt.200000.) then
540 if (.not. is_there(nint(plvl(k)), 'HGT').and. &
541 is_there(nint(plvl(k)), 'GEOPT')) then
542 call get_dims(nint(plvl(k)), 'GEOPT')
543 allocate(scr2d(map%nx,map%ny))
544 call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny)
545 scr2d = scr2d / 9.81
546 call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny)
547 call mprintf(.true.,DEBUG, &
548 "RRPR: Computing GHT from GEOPT ")
549 deallocate(scr2d)
550 endif
551 if ( (.not. is_there(nint(plvl(k)),'HGT')) .and. &
552 ( is_there(nint(plvl(k-1)), 'HGT')) ) then
553 found_level = .false.
554 hg_loop : do itest = k+1,nlvl,1
555 if ( ( is_there(nint(plvl(itest)), 'HGT')) ) THEN
556 found_level = .true.
557 exit hg_loop
558 ENDIF
559 enddo hg_loop
560 if( found_level ) then
561 call get_dims(nint(plvl(itest)), 'HGT')
562 call vntrp(plvl, maxlvl, k, itest, interp_type, "HGT ", map%nx, map%ny)
563 else
564 PRINT *,'WARNING: Could not interpolate HGT to level k=',k,' p=',plvl(k),&
565 'because could not find any level above this level.'
566 endif
567 endif
568 endif
569 enddo
571 !
572 ! If this is GFS data, we might have data at the level of max wind speed,
573 ! or the level of the tropopause. If so, we want to replicate the pressures
574 ! at those levels (new names). The replicated names are to allow the
575 ! metgrid program to interpolate the 2d pressure array with both a nearest
576 ! neighbor AND a 4-pt technique. Those two pressures are used in ARW real
577 ! for vertical interpolation of the trop and max wind level data.
578 !
580 if (index(map%source,'NCEP GFS') .ne. 0 ) then
581 call mprintf(.true.,DEBUG, &
582 "RRPR: Replicating GFS pressures for max wind and trop")
583 if ( is_there(200100,'PMAXW ') .or. &
584 is_there(200100,'PTROP ') ) then
585 call gfs_trop_maxw_pressures (map%nx, map%ny)
586 endif
587 endif
589 ! Repair GFS and ECMWF pressure-level RH
590 if (index(map%source,'NCEP GFS') .ne. 0 .or. &
591 index(map%source,'NCEP GEFS') .ne. 0 .or. &
592 index(map%source,'NCEP CDAS CFSV2') .ne. 0 .or. &
593 index(map%source,'ECMWF') .ne. 0 ) then
594 call mprintf(.true.,DEBUG, &
595 "RRPR: Adjusting RH values ")
596 do k = 1, nlvl
597 if ( is_there(nint(plvl(k)),'RH') .and. &
598 is_there(nint(plvl(k)),'TT') ) then
599 call fix_gfs_rh (map%nx, map%ny, plvl(k))
600 endif
601 enddo
602 endif
604 ! If upper-air RH is missing, see if we can compute RH from Specific Humidity:
606 do k = 1, nlvl
607 if (plvl(k).lt.200000.) then
608 if (.not. is_there(nint(plvl(k)), 'RH') .and. &
609 is_there(nint(plvl(k)), 'TT') .and. &
610 is_there(nint(plvl(k)), 'SPECHUMD')) then
611 call get_dims(nint(plvl(k)), 'TT')
612 call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k))
613 endif
614 endif
615 enddo
617 ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure:
618 ! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.)
620 do k = 1, nlvl
621 if (plvl(k).lt.200000.) then
622 if (.not. is_there(nint(plvl(k)),'RH').and. &
623 is_there(nint(plvl(k)), 'TT') .and. &
624 is_there(nint(plvl(k)),'VAPP')) then
625 call get_dims(nint(plvl(k)),'TT')
626 call compute_rh_vapp_upa(map%nx, map%ny, plvl(k))
627 endif
628 endif
629 enddo
631 ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression:
633 do k = 1, nlvl
634 if (plvl(k).lt.200000.) then
635 if (.not. is_there(nint(plvl(k)),'RH').and. &
636 is_there(nint(plvl(k)), 'TT') .and. &
637 is_there(nint(plvl(k)),'DEPR')) then
638 call get_dims(nint(plvl(k)),'TT')
639 call compute_rh_depr(map%nx, map%ny, plvl(k))
640 endif
641 endif
642 enddo
643 !
644 ! If upper-air RH is missing, see if we can interpolate from surrounding levels.
645 ! This is a simple vertical interpolation, linear or log in pressure.
646 ! Currently, this simply fills in missing levels between two present levels.
648 !
650 do k = 2, nlvl-1, 1
651 if (plvl(k-1) .lt. 200000.) then
652 if ( (.not. is_there(nint(plvl(k)),'RH')) .and. &
653 ( is_there(nint(plvl(k-1)), 'RH')) ) then
654 found_level = .false.
655 RH_loop : do itest = k+1,nlvl,1
656 if ( ( is_there(nint(plvl(itest)), 'RH')) ) then
657 found_level = .true.
658 exit RH_loop
659 endif
660 enddo RH_loop
661 if( found_level ) then
662 call get_dims(nint(plvl(itest)), 'RH')
663 call vntrp(plvl, maxlvl, k, itest, interp_type, "RH ", map%nx, map%ny)
664 else
665 PRINT *,'WARNING: Could not interpolate RH to level k=',k,' p=',plvl(k),&
666 'because could not find any level above this level.'
667 endif
668 endif
669 endif
670 enddo
672 !
673 ! Check to see if we need to fill RH above 300 mb:
674 !
675 if (is_there(30000, 'RH')) then
676 call get_dims(30000, 'RH')
677 allocate(scr2d(map%nx,map%ny))
679 do k = 1, nlvl
680 ! Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa.
681 ! The stratospheric RH will be adjusted further in real.
682 if (plvl(k).le.7000.) then
683 scr2d = 0.
684 else if (plvl(k).lt.30000.) then
685 scr2d = 5.
686 endif
687 if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then
688 ! levels higher than .1 mb are special - do not fill
689 if (.not. is_there(nint(plvl(k)), 'RH')) then
690 call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny)
691 call mprintf(.true.,DEBUG, &
692 "RRPR: RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.))
693 endif
694 endif
695 enddo
696 deallocate(scr2d)
697 endif
698 !
699 ! If surface RH is missing, see if we can compute RH from Specific Humidity
700 ! or Dewpoint or Dewpoint depression:
701 !
702 if (.not. is_there (200100, 'RH')) then
703 if (is_there(200100, 'TT').and. &
704 is_there(200100, 'PSFC' ) .and. &
705 is_there(200100, 'SPECHUMD')) then
706 call get_dims(200100, 'TT')
707 call compute_rh_spechumd(map%nx, map%ny)
708 call mprintf(.true.,DEBUG, &
709 "RRPR: SURFACE RH is computed")
710 elseif (is_there(200100, 'TT' ).and. &
711 is_there(200100, 'DEWPT')) then
712 call get_dims(200100, 'TT')
713 call compute_rh_dewpt(map%nx, map%ny)
714 elseif (is_there(200100, 'TT').and. &
715 is_there(200100, 'DEPR')) then
716 call get_dims(200100, 'TT')
717 call compute_rh_depr(map%nx, map%ny, 200100.)
718 endif
719 endif
721 !
722 ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC
723 ! (From Wei Wang, 2007 June 21, modified 12/28/2007)
724 !
725 if (.not. is_there(200100, 'SNOW') .and. &
726 is_there(200100, 'SNOWRUC')) then
727 call get_dims(200100, 'SNOWRUC')
728 allocate(scr2d(map%nx,map%ny))
729 call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny)
730 scr2d = scr2d * 1000.
731 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
732 deallocate(scr2d)
733 endif
735 ! compute snow water equivalent (SNOW) for NCEP RUC models
736 ! As of Sept. 14 2011
737 if ( index(map%source,'NCEP RUC Model') .ne. 0) then
738 if (is_there(200100, 'SNOWH') .and. .not. is_there(200100, 'SNOW')) then
739 call get_dims(200100, 'SNOWH')
740 allocate(scr2d(map%nx,map%ny))
741 call get_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
742 call mprintf(.true.,DEBUG, &
743 "RRPR: Computing SNOWH from SNOW")
744 if (is_there(200100, 'RHOSN')) then ! If we have snow density, use it to compute snowh
745 call get_dims(200100, 'RHOSN')
746 allocate(tmp2d(map%nx,map%ny))
747 call get_storage(200100, 'RHOSN', tmp2d, map%nx, map%ny)
748 scr2d = scr2d * tmp2d
749 deallocate(tmp2d)
750 else
751 scr2d = scr2d * 200.0 ! Assume 200:1 ratio
752 endif
753 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
754 deallocate(scr2d)
755 endif
756 endif
758 ! Modify the 2017 GFS masked fields
759 if (index(map%source,'NCEP GFS') .ne. 0 .or. &
760 index(map%source,'NCEP GEFS') .ne. 0 ) then
761 call mprintf(.true.,DEBUG, &
762 "RRPR: Adjusting GFS masked fields ")
763 if ( is_there(200100, 'ST000010')) then
764 call get_dims(200100, 'ST000010')
765 call fix_gfs_miss (map%nx, map%ny, 200100.)
766 endif
767 endif
769 ! Fix the 23 September 2020 GEFS landmask
770 if (index(map%source,'NCEP GEFS') .ne. 0 ) then
771 call mprintf(.true.,DEBUG, &
772 "RRPR: Adjusting GEFS landmask")
773 if ( is_there(200100, 'ST000010') .and. &
774 is_there(200100, 'LANDSEA')) then
775 call get_dims(200100, 'LANDSEA')
776 call fix_gefs_landmask (map%nx, map%ny, 200100.)
777 endif
778 endif
780 ! Add residual soil moisture to SOILM* if initialized from the GSD RUC model or from NCEP RUC
781 if (index(map%source,'NOAA GSD') .ne. 0 .or. &
782 index(map%source,'NCEP RUC Model') .ne. 0) then
783 if ( .not. is_there(200100, 'SOILM000') .and.&
784 is_there(200100, 'SM000ruc') ) then
785 call get_dims(200100, 'SM000ruc')
786 print *,'Adjust RUC soil moisture'
787 call mprintf(.true.,DEBUG, &
788 "RRPR: Adjusting RUC soil moisture ")
789 call fix_ruc_soilm (map%nx, map%ny)
790 endif
791 endif
793 ! Convert soil moisture in ICON input from kg m-2 to m3 m-3
794 if (index(map%source,'DWD') .ne. 0) then
795 if (is_there(200100, 'SOILM001')) then
796 call get_dims(200100, 'SOILM001')
797 print *,'Adjusting ICON soil moisture'
798 call mprintf(.true.,DEBUG, &
799 "RRPR: Adjusting ICON soil moisture ")
800 call fix_icon_soilm (map%nx, map%ny)
801 endif
802 endif
804 !
805 ! Check to see if we need to fill SOILHGT from SOILGEO.
806 ! (From Wei Wang, 2007 June 21)
807 !
808 if (.not. is_there(200100, 'SOILHGT') .and. &
809 is_there(200100, 'SOILGEO')) then
810 call get_dims(200100, 'SOILGEO')
811 allocate(scr2d(map%nx,map%ny))
812 call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny)
813 scr2d = scr2d / 9.81
814 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
815 call mprintf(.true.,DEBUG, &
816 "RRPR: Computing SOILGHT from SOILGEO ")
817 deallocate(scr2d)
818 endif
820 ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40)
821 if (.not. is_there(200100, 'SOILHGT') .and. &
822 is_there(1, 'SOILGEO')) then
823 call get_dims(1, 'SOILGEO')
824 allocate(scr2d(map%nx,map%ny))
825 call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny)
826 scr2d = scr2d / 9.81
827 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny)
828 deallocate(scr2d)
829 endif
831 ! For NCEP RR (using the same ID as for RUC) native-level input,
832 ! may need to move PSFC from level 1 to 2001.
833 ! From TGS 8 Sept. 2011
834 if ( index(map%source,'NCEP RUC Model') .ne. 0) then
835 if (.not. is_there(200100, 'PSFC') .and. &
836 is_there(1, 'PRESSURE')) then
837 print *,'Process PSFC for NCEP RR'
838 call get_dims(1, 'PRESSURE')
839 allocate(scr2d(map%nx,map%ny))
840 call get_storage(1, 'PRESSURE', scr2d, map%nx, map%ny)
841 call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
842 deallocate(scr2d)
843 endif
844 endif
846 ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001.
847 if ( index(map%source,'ECMWF') .ne. 0) then
848 if (.not. is_there(200100, 'PSFC') .and. &
849 is_there(1, 'PSFCH')) then
850 call get_dims(1, 'PSFCH')
851 allocate(scr2d(map%nx,map%ny))
852 call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny)
853 call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny)
854 deallocate(scr2d)
855 endif
856 endif
858 ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2
859 !
860 if (is_there(200100, 'SNOW_EC')) then
861 call get_dims(200100, 'SNOW_EC')
862 allocate(scr2d(map%nx,map%ny))
863 call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny)
864 scr2d = scr2d * 1000.
865 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
866 deallocate(scr2d)
867 endif
869 ! Convert the ECMWF LANDSEA mask from a fraction to a flag
871 if ( index(map%source,'ECMWF') .ne. 0) then
872 if (is_there(200100, 'LANDSEA')) then
873 call get_dims(200100, 'LANDSEA')
874 call make_zero_or_one(map%nx, map%ny, 'LANDSEA')
875 endif
876 endif
878 ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF.
879 ! The GFS-based reanalyses values should be OK as is.
880 if ((index(map%source,'NCEP GFS') .ne. 0 .or. &
881 index(map%source,'NCEP GEFS') .ne. 0) .and. &
882 is_there(200100, 'SNOW')) then
883 call mprintf(.true.,DEBUG, &
884 "RRPR: Recomputing SNOW for NCEP GFS")
885 call get_dims(200100, 'SNOW')
886 allocate(scr2d(map%nx,map%ny))
887 call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
888 scr2d = scr2d * 2.
889 call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
890 deallocate(scr2d)
891 endif
893 ! compute physical snow depth (SNOWH) for various models
894 ! As of March 2011, this is done here instead of real because we have model
895 ! source information.
896 if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then
897 call get_dims(200100, 'SNOW')
898 allocate(scr2d(map%nx,map%ny))
899 call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny)
900 call mprintf(.true.,DEBUG, &
901 "RRPR: Computing SNOWH from SNOW")
902 if ( index(map%source,'NCEP ') .ne. 0) then
903 scr2d = scr2d * 0.005 ! Assume 200:1 ratio as used at NCEP and in NOAH
904 else if (index(map%source,'ECMWF') .ne. 0) then
905 if (is_there(200100, 'SNOW_DEN')) then ! If we have snow density, use it to compute snowh
906 call get_dims(200100, 'SNOW_DEN')
907 allocate(tmp2d(map%nx,map%ny))
908 call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny)
909 scr2d = scr2d / tmp2d
910 deallocate(tmp2d)
911 else
912 scr2d = scr2d * 0.004 ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio).
913 endif
914 else ! Other models
915 scr2d = scr2d * 0.005 ! Use real's default method (200:1)
916 endif
917 call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny)
918 deallocate(scr2d)
919 endif
921 ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted
922 ! to the appropriate values in real depending on whether or not the polar mods are used.
924 !! If we've got a SEAICE field, make sure that it is all Zeros and Ones:
926 ! if (is_there(200100, 'SEAICE')) then
927 ! call get_dims(200100, 'SEAICE')
928 ! call make_zero_or_one(map%nx, map%ny, 'SEAICE')
929 ! endif
931 ! If we've got an ICEMASK field, re-flag it for output to met_em and real:
932 ! Field | GRIB In | Out
933 ! -------------------------
934 ! water | 0 | 0
935 ! land | -1 | 1
936 ! ice | 1 | 0
938 if (is_there(200100, 'ICEMASK')) then
939 call get_dims(200100, 'ICEMASK')
940 call re_flag_ice_mask(map%nx, map%ny)
941 endif
943 ! If we have an ICEFRAC field, convert from % to fraction
944 if (is_there(200100, 'ICEFRAC')) then
945 call get_dims(200100, 'ICEFRAC')
946 allocate(scr2d(map%nx,map%ny))
947 call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
948 scr2d = scr2d / 100.
949 call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny)
950 deallocate(scr2d)
951 endif
954 call mprintf(.true.,INFORM, &
955 "RRPR: hdate = %s ",s1=hdate)
956 call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level)
957 call clear_storage
958 exit FILELOOP
959 enddo FILELOOP
960 enddo TIMELOOP
961 end subroutine rrpr
963 subroutine make_zero_or_one(ix, jx, infield)
964 ! Make sure the input field (SEAICE or LANDSEA) is zero or one.
965 use storage_module
966 implicit none
967 integer :: ix, jx
968 real, dimension(ix,jx) :: seaice
969 character(len=*) :: infield
971 call get_storage(200100, infield, seaice, ix, jx)
972 where(seaice > 0.5)
973 seaice = 1.0
974 elsewhere
975 seaice = 0.0
976 end where
977 call put_storage(200100, infield, seaice, ix, jx)
978 end subroutine make_zero_or_one
980 subroutine re_flag_ice_mask(ix, jx)
981 !
982 ! Change land points from -1 to 1
983 ! Change ice points from 1 to 0
984 ! Water points stay at 0
985 !
986 use storage_module
987 implicit none
988 integer :: ix, jx
989 real, dimension(ix,jx) :: iceflag
991 call get_storage(200100, 'ICEMASK',iceflag, ix, jx)
992 where(iceflag > 0.5) ! Ice points, set to water value
993 iceflag = 0.0
994 end where
995 where(iceflag < -0.5) ! Land points
996 iceflag = 1.0
997 end where
998 call put_storage(200100, 'ICEMASK',iceflag, ix, jx)
999 end subroutine re_flag_ice_mask
1001 subroutine compute_spechumd_qvapor(ix, jx, plvl)
1002 ! Compute specific humidity from water vapor mixing ratio.
1003 use storage_module
1004 implicit none
1005 integer :: ix, jx
1006 real :: plvl
1007 real, dimension(ix,jx) :: QVAPOR, SPECHUMD
1009 call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx)
1011 SPECHUMD = QVAPOR/(1.+QVAPOR)
1013 call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx)
1014 if(nint(plvl).eq.1) then
1015 call put_storage(200100,'SPECHUMD', spechumd, ix, jx)
1016 endif
1018 end subroutine compute_spechumd_qvapor
1020 subroutine compute_t_vptmp(ix, jx, plvl)
1021 ! Compute temperature from virtual potential temperature
1022 use storage_module
1023 implicit none
1024 integer :: ix, jx
1025 real :: plvl
1026 real, dimension(ix,jx) :: T, VPTMP, P, Q
1028 real, parameter :: rovcp=0.28571
1030 call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx)
1031 IF (nint(plvl) .LT. 200) THEN
1032 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
1033 ELSE
1034 p = plvl
1035 ENDIF
1036 call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
1038 t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q))
1040 call put_storage(nint(plvl), 'TT', t, ix, jx)
1041 if(nint(plvl).eq.1) then
1042 call put_storage(200100, 'PSFC', p, ix, jx)
1043 endif
1045 end subroutine compute_t_vptmp
1048 subroutine compute_rh_spechumd(ix, jx)
1049 ! Compute relative humidity from specific humidity.
1050 use storage_module
1051 implicit none
1052 integer :: ix, jx
1053 real, dimension(ix,jx) :: T, P, RH, Q
1055 real, parameter :: svp1=611.2
1056 real, parameter :: svp2=17.67
1057 real, parameter :: svp3=29.65
1058 real, parameter :: svpt0=273.15
1059 real, parameter :: eps = 0.622
1061 call get_storage(200100, 'TT', T, ix, jx)
1062 call get_storage(200100, 'PSFC', P, ix, jx)
1063 call get_storage(200100, 'SPECHUMD', Q, ix, jx)
1065 rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1067 call put_storage(200100, 'RH', rh, ix, jx)
1069 end subroutine compute_rh_spechumd
1071 subroutine compute_rh_spechumd_upa(ix, jx, plvl)
1072 ! Compute relative humidity from specific humidity.
1073 use storage_module
1074 implicit none
1075 integer :: ix, jx
1076 real :: plvl
1077 real, dimension(ix,jx) :: T, P, RH, Q
1079 real, parameter :: svp1=611.2
1080 real, parameter :: svp2=17.67
1081 real, parameter :: svp3=29.65
1082 real, parameter :: svpt0=273.15
1083 real, parameter :: eps = 0.622
1085 IF ( nint(plvl).LT. 200) THEN
1086 if (is_there(nint(plvl), 'PRESSURE')) then
1087 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
1088 else
1089 return ! if we don't have pressure on model levels, return
1090 endif
1091 ELSE
1092 P = plvl
1093 ENDIF
1094 call get_storage(nint(plvl), 'TT', T, ix, jx)
1095 call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx)
1096 Q=MAX(1.E-10,Q)
1098 rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1100 call put_storage(nint(plvl), 'RH', rh, ix, jx)
1102 end subroutine compute_rh_spechumd_upa
1104 subroutine compute_rh_vapp_upa(ix, jx, plvl)
1105 ! Compute relative humidity from vapor pressure.
1106 ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27.
1107 use storage_module
1108 implicit none
1109 integer :: ix, jx
1110 real :: plvl
1111 real, dimension(ix,jx) :: P, ES
1112 real, pointer, dimension(:,:) :: T, E, RH
1114 real, parameter :: svp1=611.2
1115 real, parameter :: svp2=17.67
1116 real, parameter :: svp3=29.65
1117 real, parameter :: svpt0=273.15
1119 allocate(RH(ix,jx))
1121 IF ( nint(plvl).LT. 200) THEN
1122 if (is_there(nint(plvl), 'PRESSURE')) then
1123 call get_storage(nint(plvl), 'PRESSURE', P, ix, jx)
1124 else
1125 return ! if we don't have pressure on model levels, return
1126 endif
1127 ELSE
1128 P = plvl
1129 ENDIF
1130 call refr_storage(nint(plvl), 'TT', T, ix, jx)
1131 call refr_storage(nint(plvl), 'VAPP', E, ix, jx)
1133 ES=svp1*exp(svp2*(T-svpt0)/(T-svp3))
1134 rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2)
1136 call refw_storage(nint(plvl), 'RH', rh, ix, jx)
1138 nullify(T,E)
1140 end subroutine compute_rh_vapp_upa
1142 subroutine compute_rh_depr(ix, jx, plvl)
1143 ! Compute relative humidity from Dewpoint Depression
1144 use storage_module
1145 implicit none
1146 integer :: ix, jx
1147 real :: plvl
1148 real, dimension(ix,jx) :: t, depr, rh
1150 real, parameter :: Xlv = 2.5e6
1151 real, parameter :: Rv = 461.5
1153 integer :: i, j
1155 call get_storage(nint(plvl), 'TT', T, ix, jx)
1156 call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx)
1158 where(DEPR < 100.)
1159 rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2
1160 elsewhere
1161 rh = 0.0
1162 endwhere
1164 call put_storage(nint(plvl),'RH ', rh, ix, jx)
1166 end subroutine compute_rh_depr
1168 subroutine compute_rh_dewpt(ix,jx)
1169 ! Compute relative humidity from Dewpoint
1170 use storage_module
1171 implicit none
1172 integer :: ix, jx
1173 real, dimension(ix,jx) :: t, dp, rh
1175 real, parameter :: Xlv = 2.5e6
1176 real, parameter :: Rv = 461.5
1178 call get_storage(200100, 'TT ', T, ix, jx)
1179 call get_storage(200100, 'DEWPT ', DP, ix, jx)
1181 rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2
1183 call put_storage(200100,'RH ', rh, ix, jx)
1185 end subroutine compute_rh_dewpt
1187 subroutine gfs_trop_maxw_pressures(ix,jx)
1188 ! These are duplicate pressure values from the GFS, for
1189 ! the level of max wind speed and for the trop level.
1190 ! The duplicates are saved with a different name, so that
1191 ! the metgrid program can horizontally interpolate them to
1192 ! the model domain with a nearest neighbor method.
1193 use storage_module
1194 implicit none
1195 integer :: ix, jx
1196 real, dimension(ix,jx) :: pmaxw, pmaxwnn, ptrop, ptropnn
1198 if ( is_there(200100, 'PMAXW ') ) then
1199 call get_storage(200100, 'PMAXW ', pmaxw , ix, jx)
1200 pmaxwnn = pmaxw
1201 call put_storage(200100, 'PMAXWNN ', pmaxwnn, ix, jx)
1202 end if
1204 if ( is_there(200100, 'PTROP ') ) then
1205 call get_storage(200100, 'PTROP ', ptrop , ix, jx)
1206 ptropnn = ptrop
1207 call put_storage(200100, 'PTROPNN ', ptropnn, ix, jx)
1208 end if
1210 end subroutine gfs_trop_maxw_pressures
1212 subroutine vntrp(plvl, maxlvl, k, k2, interp_type, name, ix, jx)
1213 use storage_module
1214 use module_debug
1215 implicit none
1216 integer :: ix, jx, k, k2, maxlvl
1217 real, dimension(maxlvl) :: plvl
1218 character(len=8) :: name
1219 real, dimension(ix,jx) :: a, b, c
1220 real :: frc
1221 integer :: interp_type
1223 write(*,'("Interpolating to fill in ", A, " at level ", F8.2, " hPa using levels ", F8.2," hPa and ",&
1224 & F8.2," hPa")') trim(name), plvl(k)/100.0,plvl(k-1)/100.0, plvl(k2)/100.0
1225 call mprintf(.true.,INFORM, &
1226 "RRPR: Interpolating to fill in %s at %i Pa using levels %i Pa and %i Pa",&
1227 s1=trim(name), i1=int(plvl(k)), i2=int(plvl(k-1)), i3=int(plvl(k2)))
1229 call get_storage(nint(plvl(k-1)), name, a, ix, jx)
1230 call get_storage(nint(plvl(k2)), name, c, ix, jx)
1232 if ( interp_type == 1 ) then
1233 frc = log(plvl(k)/plvl(k2)) / log(plvl(k-1)/plvl(k2))
1234 else
1235 frc = (plvl(k) - plvl(k2)) / (plvl(k-1)-plvl(k2))
1236 endif
1238 b = (1.-frc)*c + frc*a
1240 !KWM b = 0.5 * (a + c)
1241 call put_storage(nint(plvl(k)), name, b, ix, jx)
1243 end subroutine vntrp
1246 subroutine fix_gfs_miss (ix, jx, plvl)
1247 ! This routine replaces July 2017 GFS missing values with the WPS one.
1248 ! Earlier GFS files are unmodified.
1249 ! As of 2017, NCEP changed 'ocean' values for masked fields (ST, SM, SNOWH, etc.)
1250 ! from 0 to something other than missing. We will assume any large value
1251 ! (greater than 10^^18) is a missing code.
1252 ! We reset it to the WPS missing value (as set in METGRID.TBL)
1253 ! Changes described in http://www.nws.noaa.gov/os/notification/scn17-67gfsupgradeaaa.htm
1254 ! plvl must always be 200100.
1255 ! While, technically, any 3-d field could have a missing value in it we only deal with
1256 ! the surface fields which are known to have missing values.
1257 !
1258 use storage_module
1259 implicit none
1260 integer :: ix, jx, i, j, k
1261 real :: plvl
1262 real, allocatable, dimension(:,:) :: f, sea
1263 integer, parameter :: nvar = 10
1264 character, dimension(nvar) :: flist*8
1265 data flist/'ST000010','ST010040','ST040100','ST100200','ST010200', &
1266 'SM000010','SM010040','SM040100','SM100200','SM010200' /
1267 allocate(sea(ix,jx))
1268 allocate(f(ix,jx))
1269 ! If LANDN is present (July 2017 and later GFS output), use it for LANDSEA
1270 if ( is_there(200100, 'LANDN') ) then
1271 call get_storage(200100, 'LANDN', sea, ix, jx)
1272 call put_storage(200100, 'LANDSEA', sea, ix, jx)
1273 endif
1274 do k = 1, nvar
1275 if (is_there(200100, flist(k) )) then
1276 call get_storage(nint(plvl), flist(k), f, ix, jx)
1277 do j = 1, jx
1278 do i = 1, ix
1279 if (abs(f(i,j)) .gt. 1.e18) then
1280 f(i,j) = -1.e30
1281 endif
1282 enddo
1283 enddo
1284 ! Limit soil moisture to .468 (should only occur for permanent land ice points
1285 ! according to NCEP documentation.)
1286 if (flist(k)(1:2) .eq. 'SM' ) then
1287 do j = 1, jx
1288 do i = 1, ix
1289 if ((f(i,j)) .gt. 0.468) then
1290 f(i,j) = 0.468
1291 endif
1292 enddo
1293 enddo
1294 endif
1295 call put_storage(200100, flist(k), f, ix, jx)
1296 endif
1297 enddo
1298 ! Snow fields have a different WPS missing value and must be adjusted
1299 ! separately from the soil fields (unless we want to have an array of missing values,too).
1300 if (is_there(200100, 'SNOW' )) then
1301 call get_storage(200100, 'SNOW', f, ix, jx)
1302 do j = 1, jx
1303 do i = 1, ix
1304 if (abs(f(i,j)) .gt. 1.e18) then
1305 f(i,j) = 0
1306 endif
1307 enddo
1308 enddo
1309 call put_storage(200100, 'SNOW', f, ix, jx)
1310 endif
1311 if (is_there(200100, 'SNOWH' )) then
1312 call get_storage(200100, 'SNOWH', f, ix, jx)
1313 do j = 1, jx
1314 do i = 1, ix
1315 if (abs(f(i,j)) .gt. 1.e18) then
1316 f(i,j) = 0
1317 endif
1318 enddo
1319 enddo
1320 call put_storage(200100, 'SNOWH', f, ix, jx)
1321 endif
1322 deallocate (f)
1323 deallocate (sea)
1324 end subroutine fix_gfs_miss
1326 subroutine fix_gfs_rh (ix, jx, plvl)
1327 ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe).
1328 use storage_module
1329 implicit none
1330 integer :: ix, jx, i, j
1331 real :: plvl, eis, ews, r
1332 real, allocatable, dimension(:,:) :: rh, tt
1334 allocate(rh(ix,jx))
1335 allocate(tt(ix,jx))
1336 call get_storage(nint(plvl), 'RH', rh, ix, jx)
1337 call get_storage(nint(plvl), 'TT', tt, ix, jx)
1338 do j = 1, jx
1339 do i = 1, ix
1340 if ( tt(i,j) .le. 273.15 ) then
1341 ! Murphy and Koop 2005 ice saturation vapor pressure.
1342 ! eis and ews in hPA, tt is in K
1343 eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) &
1344 - (0.00728332 * tt(i,j)))
1345 ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most
1346 ! formulae are very similar from 0 to -20, so we don't need a more exact formula.
1348 ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5))
1349 if ( tt(i,j) .gt. 253.15 ) then
1350 ! A linear approximation to the GFS blending region ( -20 > T < 0 )
1351 r = ((273.15 - tt(i,j)) / 20.)
1352 r = (r * eis) + ((1-r)*ews)
1353 else
1354 r = eis
1355 endif
1356 rh(i,j) = rh(i,j) * (r / ews)
1357 endif
1358 enddo
1359 enddo
1360 call put_storage(nint(plvl), 'RH', rh, ix, jx)
1361 deallocate (rh)
1362 deallocate (tt)
1363 end subroutine fix_gfs_rh
1366 subroutine fix_ruc_soilm (ix, jx)
1367 ! This routine adds residual soil moisture if initialized fron RUC
1368 use storage_module
1369 implicit none
1370 integer :: ix, jx, i, j
1371 REAL , DIMENSION(100) :: lqmi
1372 real, allocatable, dimension(:,:) :: soilm000, soilm005, soilm020, &
1373 soilm040, soilm160, soilm300,soilcat
1374 allocate(soilm000(ix,jx))
1375 allocate(soilm005(ix,jx))
1376 allocate(soilm020(ix,jx))
1377 allocate(soilm040(ix,jx))
1378 allocate(soilm160(ix,jx))
1379 allocate(soilm300(ix,jx))
1380 allocate(soilcat(ix,jx))
1381 call get_storage(200100, 'SM000ruc', soilm000, ix, jx)
1382 call get_storage(200100, 'SM005ruc', soilm005, ix, jx)
1383 call get_storage(200100, 'SM020ruc', soilm020, ix, jx)
1384 call get_storage(200100, 'SM040ruc', soilm040, ix, jx)
1385 call get_storage(200100, 'SM160ruc', soilm160, ix, jx)
1386 call get_storage(200100, 'SM300ruc', soilm300, ix, jx)
1388 call get_storage(200100, 'SOILCAT', soilcat, ix, jx)
1390 lqmi(1:16) = &
1391 (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1392 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1393 0.004, 0.065 /)
1395 do j = 1, jx
1396 do i = 1, ix
1398 SOILM000(i,j)=SOILM000(i,j) + lqmi(nint(soilcat(i,j)))
1399 SOILM005(i,j)=SOILM005(i,j) + lqmi(nint(soilcat(i,j)))
1400 SOILM020(i,j)=SOILM020(i,j) + lqmi(nint(soilcat(i,j)))
1401 SOILM040(i,j)=SOILM040(i,j) + lqmi(nint(soilcat(i,j)))
1402 SOILM160(i,j)=SOILM160(i,j) + lqmi(nint(soilcat(i,j)))
1403 SOILM300(i,j)=SOILM300(i,j) + lqmi(nint(soilcat(i,j)))
1404 enddo
1405 enddo
1406 call put_storage(200100, 'SOILM000', soilm000, ix, jx)
1407 call put_storage(200100, 'SOILM005', soilm005, ix, jx)
1408 call put_storage(200100, 'SOILM020', soilm020, ix, jx)
1409 call put_storage(200100, 'SOILM040', soilm040, ix, jx)
1410 call put_storage(200100, 'SOILM160', soilm160, ix, jx)
1411 call put_storage(200100, 'SOILM300', soilm300, ix, jx)
1413 print *,'fix_ruc_soilm is done!'
1415 deallocate(soilm000)
1416 deallocate(soilm005)
1417 deallocate(soilm020)
1418 deallocate(soilm040)
1419 deallocate(soilm160)
1420 deallocate(soilm300)
1421 deallocate(soilcat)
1423 end subroutine fix_ruc_soilm
1425 subroutine fix_icon_soilm (ix, jx)
1426 use storage_module
1427 implicit none
1428 integer :: ix, jx
1429 real, allocatable, dimension(:,:) :: soilm001, soilm002, soilm006, &
1430 soilm018, soilm054, soilm162, soilm486, soilm999
1431 allocate(soilm001(ix,jx))
1432 allocate(soilm002(ix,jx))
1433 allocate(soilm006(ix,jx))
1434 allocate(soilm018(ix,jx))
1435 allocate(soilm054(ix,jx))
1436 allocate(soilm162(ix,jx))
1437 allocate(soilm486(ix,jx))
1438 allocate(soilm999(ix,jx))
1439 call get_storage(200100, 'SOILM001', soilm001, ix, jx)
1440 call get_storage(200100, 'SOILM002', soilm002, ix, jx)
1441 call get_storage(200100, 'SOILM006', soilm006, ix, jx)
1442 call get_storage(200100, 'SOILM018', soilm018, ix, jx)
1443 call get_storage(200100, 'SOILM054', soilm054, ix, jx)
1444 call get_storage(200100, 'SOILM162', soilm162, ix, jx)
1445 call get_storage(200100, 'SOILM486', soilm486, ix, jx)
1446 call get_storage(200100, 'SOILM999', soilm999, ix, jx)
1447 !convert kg m-2 to m3 m-3
1448 soilm001 = soilm001 / 0.01 / 1000.
1449 soilm002 = soilm002 / 0.02 / 1000.
1450 soilm006 = soilm006 / 0.06 / 1000.
1451 soilm018 = soilm018 / 0.18 / 1000.
1452 soilm054 = soilm054 / 0.54 / 1000.
1453 soilm162 = soilm162 / 1.62 / 1000.
1454 soilm486 = soilm486 / 4.86 / 1000.
1455 soilm999 = soilm999 / 14.58 / 1000.
1456 call put_storage(200100, 'SOILM001', soilm001, ix, jx)
1457 call put_storage(200100, 'SOILM002', soilm002, ix, jx)
1458 call put_storage(200100, 'SOILM006', soilm006, ix, jx)
1459 call put_storage(200100, 'SOILM018', soilm018, ix, jx)
1460 call put_storage(200100, 'SOILM054', soilm054, ix, jx)
1461 call put_storage(200100, 'SOILM162', soilm162, ix, jx)
1462 call put_storage(200100, 'SOILM486', soilm486, ix, jx)
1463 call put_storage(200100, 'SOILM999', soilm999, ix, jx)
1465 print *,'fix_icon_soilm is done!'
1467 deallocate(soilm001)
1468 deallocate(soilm002)
1469 deallocate(soilm006)
1470 deallocate(soilm018)
1471 deallocate(soilm054)
1472 deallocate(soilm162)
1473 deallocate(soilm486)
1474 deallocate(soilm999)
1476 end subroutine fix_icon_soilm
1478 subroutine fix_gefs_landmask (ix, jx, plvl)
1479 ! Set LANDSEA to values based on soil temp. Assumes fix_gfs_miss has already been called.
1480 ! Needed for September 23, 2020 changes to GEFS.
1481 ! Earlier GEFS files are unmodified. Must be called after fix_gfs_miss
1482 ! plvl must always be 200100.
1483 !
1484 use storage_module
1485 use module_debug
1486 implicit none
1487 integer :: ix, jx, i, j, k
1488 real :: plvl
1489 real, allocatable, dimension(:,:) :: f, sea
1490 !
1491 ! If LANDN is present (not currently the case) then exit
1492 if ( is_there(200100, 'LANDN') ) then
1493 return
1494 else
1495 allocate(sea(ix,jx))
1496 allocate(f(ix,jx))
1497 if (is_there(nint(plvl), 'ST000010' )) then
1498 call get_storage(nint(plvl), 'ST000010', f, ix, jx)
1499 call get_storage(nint(plvl), 'LANDSEA', sea, ix, jx)
1500 do j = 1, jx
1501 do i = 1, ix
1502 if (f(i,j) .le. 0.) then
1503 sea(i,j) = 0.
1504 else
1505 sea(i,j) = 1.
1506 endif
1507 enddo
1508 enddo
1509 call put_storage(200100, 'LANDSEA', sea, ix, jx)
1510 else
1511 call mprintf(.true.,INFORM, &
1512 "RRPR: fix_gefs_landmask: ST000010 not found, LANDSEA is not modified")
1513 endif
1514 deallocate (f)
1515 deallocate (sea)
1516 endif
1517 end subroutine fix_gefs_landmask