From 500527b3f23d1b9b66346b08e18da855ebc34c89 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 14 Apr 2017 18:43:28 -0600 Subject: [PATCH] Add ungrib capability to interpolate to new pressure levels (#14) This merge adds changes that allow for the interpolation to a set of additional, runtime-defined pressure levels with the ungrib program. Now, the &ungrib namelist record contains three new options: * add_lvls: logical, whether to interpolate to additional levels * interp_type: integer, 0=linear in pressure 1=linear in log pressure * new_plvl: additional pressure levels (in Pa) to interpolate to These modifications were provided by Brian Reen (US Army Research Lab) and were based on earlier code provided to Penn State by Cindy Bruyere (NCAR). Brian notes that: "This capability is beneficial when quality controlling observations via Obsgrid (e.g., for observation nudging or verification). Obsgrid interpolates multi-level observations to the vertical levels of the coarse grid model used to drive WRF (e.g., GFS) and then compares the interpolated observed profile to the coarse grid model for quality control purposes. Thus, increased vertical resolution in the coarse grid model data results in increased vertical resolution in the interpolated observed profile used for obs nudging or verification; this allows the observed vertical structures to be better represented in the interpolated profile used for verification of obs nudging. The additional vertical resolution in the coarse grid model data also improves processing of single-level above-surface observations (e.g., aircraft obs) by providing model data closer to the level of the aircraft ob." --- ungrib/src/read_namelist.F | 47 ++++++- ungrib/src/rrpr.F | 325 ++++++++++++++++++++++++++++++++++++++++----- ungrib/src/ungrib.F | 39 +++++- 3 files changed, 368 insertions(+), 43 deletions(-) diff --git a/ungrib/src/read_namelist.F b/ungrib/src/read_namelist.F index aeef36c..b931c7c 100644 --- a/ungrib/src/read_namelist.F +++ b/ungrib/src/read_namelist.F @@ -1,5 +1,6 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& - ordered_by_date, debug_level, out_format, prefix) + ordered_by_date, debug_level, out_format, prefix, & + add_lvls, new_plvl_in, interp_type) use misc_definitions_module use module_debug @@ -12,6 +13,9 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& integer :: ntimes logical :: ordered_by_date integer :: debug_level + real, dimension(:) :: new_plvl_in + logical :: add_lvls + integer :: interp_type integer :: ierr integer :: idts @@ -46,6 +50,11 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& character(len=MAX_FILENAME_LEN) :: prefix logical :: nocolons + real :: target_end, incr + integer :: il + + real, dimension(:), allocatable :: new_plvl + namelist /share/ wrf_core, max_dom, & start_year, start_month, start_day, start_hour, & start_minute, start_second, & @@ -59,7 +68,10 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& nocolons namelist /ungrib/ out_format, & - ordered_by_date, prefix + ordered_by_date, prefix, & + add_lvls, new_plvl, interp_type + + allocate(new_plvl(size(new_plvl_in))) start_year = 0 start_month = 0 @@ -82,6 +94,10 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& debug_level = 0 nocolons = .false. + add_lvls = .false. + new_plvl = -99999. + interp_type = 0 + ! Start routine: ! Build the namelist file name: @@ -249,8 +265,35 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& 'read_namelist: I do not recognize the output format, %s , stopping.',s1=out_format) endif +! Check to see if I should create my own set of new_plvl's + if ( add_lvls .AND. new_plvl(2) > -99999. .AND. new_plvl(2) < 0.0 ) then + target_end = abs(new_plvl(2)) + incr = new_plvl(3) + il = 2 + make_plvl : do + if(il.gt.size(new_plvl)) then + call mprintf(.true.,ERROR,& + 'Too many new levels specified via new_plvl. Increase maxlvl in ungrib.F') + end if + new_plvl(il) = new_plvl(il-1) - incr + ! If we are past the end of the range of pressures over which new levels + ! are to be created, then discard the pressure we just calculated. + ! This occurs when the user-chosen increment did not evenly divide + ! the range of pressures over which new pressures were to be added. + if ( new_plvl(il) < target_end ) then + new_plvl(il) = -99999. + exit make_plvl + end if + if ( new_plvl(il) == target_end ) exit make_plvl + il = il + 1 + end do make_plvl + endif + ! Close the namelist file: close(10) + + new_plvl_in(:) = new_plvl(:) + deallocate(new_plvl) end subroutine read_namelist diff --git a/ungrib/src/rrpr.F b/ungrib/src/rrpr.F index 50ba18b..d4d89b2 100644 --- a/ungrib/src/rrpr.F +++ b/ungrib/src/rrpr.F @@ -1,4 +1,6 @@ -subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix) +subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, & + add_lvls, new_plvl, interp_type, & + debug_level, out_format, prefix) ! ! ! In case you are wondering, RRPR stands for "Read, ReProcess, and wRite" ! ! ! @@ -36,6 +38,20 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f ! PLVL: Array of pressure levels (Pa) in the dataset real , dimension(maxlvl) :: plvl +! NEW_PLVL: Array of the additional pressure levels (Pa) to interpolate to + real , dimension(maxlvl) :: new_plvl + +! TLVL: Array combining pressure levels (Pa) in PLVL and NEW_PLVL + real , dimension(maxlvl) :: tlvl + +! ADD_LVLS: Should we add levels via interpolation? + logical :: add_lvls + +! INTERP_TYPE: vertical Interpolation type +! (1=log in pressure, anything else=linear in pressure) + integer :: interp_type + + ! DEBUG_LEVEL: Integer level of debug printing (from namelist) integer :: debug_level @@ -47,6 +63,7 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f real, pointer, dimension(:,:) :: ptr2d integer :: k, kk, mm, n, ierr, ifv + integer :: itest, nn, nl, lvls, tvls integer :: iunit=13 character(LEN=19) :: hdate, hend @@ -58,6 +75,11 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f integer :: ntime, idts + logical :: found_level + real, dimension(maxlvl) :: new_plvl_to_sort + integer :: largest_number_loc + integer :: new_plvl_counter + ! DATELEN: length of date strings to use for our output file names. integer :: datelen @@ -205,6 +227,47 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f ! Retrieve the number of levels in storage: ! call get_plvls(plvl, maxlvl, nlvl) + +! Merge list of pressure levels in data with requested pressure levels + if ( add_lvls ) then + ! The merging code expects the user-defined pressure levels to be in + ! order from highest pressure to lowest pressure. + ! Sort the user-defined pressure levels accordingly + new_plvl_to_sort = new_plvl + ! Set array containing pressure levels to add to the default value set in read_namelist.F + new_plvl = -99999 + DO new_plvl_counter = 1,maxlvl + largest_number_loc = MAXLOC(new_plvl_to_sort, DIM=1) + new_plvl(new_plvl_counter)=new_plvl_to_sort(largest_number_loc) + new_plvl_to_sort(largest_number_loc)=-99999. + END DO + + tvls = 1 + lvls = 1 + loop_nvls : do nn=1,maxlvl + loop_lvls : do nl=lvls,maxlvl + if ( tvls > maxlvl ) then + call mprintf(.true.,ERROR, "Adding user-defined pressure levels resulted in too & + &many total pressure levels. Please increase maxlvl in ungrib.F") + endif + if ( plvl(nn) > 0.0 .AND. plvl(nn) >= new_plvl(nl) ) then + tlvl(tvls) = plvl(nn) + tvls = tvls + 1 + if ( plvl(nn) == new_plvl(nl) ) lvls = lvls + 1 + exit loop_lvls + endif + if ( plvl(nn) > 0.0 .AND. plvl(nn) < new_plvl(nl) ) then + tlvl(tvls) = new_plvl(nl) + tvls = tvls + 1 + lvls = lvls + 1 + endif + if ( plvl(nn) < 0.0 ) exit loop_nvls + enddo loop_lvls + enddo loop_nvls + plvl = tlvl + nlvl = tvls - 1 + end if + ! ! Fill the surface level (code 200100) from higher 200100s, as necessary ! @@ -237,34 +300,57 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f ! ! If upper-air U is missing, see if we can interpolate from surrounding levels. -! This is a simple vertical interpolation, linear in pressure. -! Currently, this simply fills in one missing level between two present levels. +! This is a simple vertical interpolation, linear or log in pressure. +! Currently, this simply fills in missing levels between two present levels. ! do k = 2, nlvl-1, 1 if (plvl(k-1) .lt. 200000.) then if ( (.not. is_there(nint(plvl(k)),'UU')) .and. & - ( is_there(nint(plvl(k-1)), 'UU')) .and.& - ( is_there(nint(plvl(k+1)), 'UU')) ) then - call get_dims(nint(plvl(k+1)), 'UU') - call vntrp(plvl, maxlvl, k, "UU ", map%nx, map%ny) + ( is_there(nint(plvl(k-1)), 'UU')) ) then + found_level = .false. + uu_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'UU')) ) then + found_level = .true. + exit uu_loop + endif + enddo uu_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'UU') + call vntrp(plvl, maxlvl, k, itest, interp_type, "UU ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate UU to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif endif endif enddo ! ! If upper-air V is missing, see if we can interpolate from surrounding levels. -! This is a simple vertical interpolation, linear in pressure. -! Currently, this simply fills in one missing level between two present levels. +! This is a simple vertical interpolation, linear or log in pressure. +! Currently, this simply fills in missing levels between two present levels. ! + do k = 2, nlvl-1, 1 if (plvl(k-1) .lt. 200000.) then if ( (.not. is_there(nint(plvl(k)),'VV')) .and. & - ( is_there(nint(plvl(k-1)), 'VV')) .and.& - ( is_there(nint(plvl(k+1)), 'VV')) ) then - call get_dims(nint(plvl(k+1)), 'VV') - call vntrp(plvl, maxlvl, k, "VV ", map%nx, map%ny) + ( is_there(nint(plvl(k-1)), 'VV')) ) then + found_level = .false. + VV_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'VV')) ) then + found_level = .true. + exit VV_loop + endif + enddo VV_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'VV') + call vntrp(plvl, maxlvl, k, itest, interp_type, "VV ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate VV to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif endif endif enddo @@ -299,24 +385,156 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f !!! ! ! If upper-air T is missing, see if we can interpolate from surrounding levels. -! This is a simple vertical interpolation, linear in pressure. -! Currently, this simply fills in one missing level between two present levels. +! This is a simple vertical interpolation, linear or log in pressure. +! Currently, this simply fills in missing levels between two present levels. ! do k = 2, nlvl-1, 1 if (plvl(k-1) .lt. 200000.) then if ( (.not. is_there(nint(plvl(k)),'TT')) .and. & - ( is_there(nint(plvl(k-1)), 'TT')) .and.& - ( is_there(nint(plvl(k+1)), 'TT')) ) then - call get_dims(nint(plvl(k+1)), 'TT') - call vntrp(plvl, maxlvl, k, "TT ", map%nx, map%ny) + ( is_there(nint(plvl(k-1)), 'TT')) ) then + found_level = .false. + TT_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'TT')) ) then + found_level = .true. + exit TT_loop + endif + enddo TT_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'TT') + call vntrp(plvl, maxlvl, k, itest, interp_type, "TT ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate TT to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif + endif + enddo + +! Vertically interpolate to fill in other moisture variables +! It seems that ultimately this should be wrapped in a function and probably loop over +! all variables that can be vertically interpolated + +! QC + do k = 2, nlvl-1, 1 + if (plvl(k-1) .lt. 200000.) then + if ( (.not. is_there(nint(plvl(k)),'QC')) .and. & + ( is_there(nint(plvl(k-1)), 'QC')) ) then + found_level = .false. + QC_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'QC')) ) then + found_level = .true. + exit QC_loop + endif + enddo QC_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'QC') + call vntrp(plvl, maxlvl, k, itest, interp_type, "QC ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate QC to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif endif endif enddo +! QR + do k = 2, nlvl-1, 1 + if (plvl(k-1) .lt. 200000.) then + if ( (.not. is_there(nint(plvl(k)),'QR')) .and. & + ( is_there(nint(plvl(k-1)), 'QR')) ) then + found_level = .false. + QR_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'QR')) ) then + found_level = .true. + exit QR_loop + endif + enddo QR_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'QR') + call vntrp(plvl, maxlvl, k, itest, interp_type, "QR ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate QR to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif + endif + enddo + +! QS + do k = 2, nlvl-1, 1 + if (plvl(k-1) .lt. 200000.) then + if ( (.not. is_there(nint(plvl(k)),'QS')) .and. & + ( is_there(nint(plvl(k-1)), 'QS')) ) then + found_level = .false. + QS_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'QS')) ) then + found_level = .true. + exit QS_loop + endif + enddo QS_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'QS') + call vntrp(plvl, maxlvl, k, itest, interp_type, "QS ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate QS to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif + endif + enddo + +! QG + do k = 2, nlvl-1, 1 + if (plvl(k-1) .lt. 200000.) then + if ( (.not. is_there(nint(plvl(k)),'QG')) .and. & + ( is_there(nint(plvl(k-1)), 'QG')) ) then + found_level = .false. + QG_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'QG')) ) then + found_level = .true. + exit QG_loop + endif + enddo QG_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'QG') + call vntrp(plvl, maxlvl, k, itest, interp_type, "QG ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate QG to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif + endif + enddo + + ! ! Check to see if we need to fill HGT from GEOPT. ! +! First make sure no GEOPT is missing + + do k = 2, nlvl-1, 1 + if (plvl(k-1) .lt. 200000.) then + if ( (.not. is_there(nint(plvl(k)),'GEOPT')) .and. & + ( is_there(nint(plvl(k-1)), 'GEOPT')) ) then + found_level = .false. + gg_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'GEOPT')) ) then + found_level = .true. + exit gg_loop + endif + enddo gg_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'GEOPT') + call vntrp(plvl, maxlvl, k, itest, interp_type, "GEOPT ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate GEOPT to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif + endif + enddo + do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)), 'HGT').and. & @@ -327,14 +545,31 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f scr2d = scr2d / 9.81 call put_storage(nint(plvl(k)), 'HGT', scr2d, map%nx, map%ny) call mprintf(.true.,DEBUG, & - "RRPR: Computing GHT from GEOPT ") + "RRPR: Computing GHT from GEOPT ") deallocate(scr2d) endif + if ( (.not. is_there(nint(plvl(k)),'HGT')) .and. & + ( is_there(nint(plvl(k-1)), 'HGT')) ) then + found_level = .false. + hg_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'HGT')) ) THEN + found_level = .true. + exit hg_loop + ENDIF + enddo hg_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'HGT') + call vntrp(plvl, maxlvl, k, itest, interp_type, "HGT ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate HGT to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif + endif endif enddo ! -! If this is GFS data, we might have data at the level of max wind speed, +! If this is GFS data, we might have data at the level of max wind speed, ! or the level of the tropopause. If so, we want to replicate the pressures ! at those levels (new names). The replicated names are to allow the ! metgrid program to interpolate the 2d pressure array with both a nearest @@ -406,19 +641,29 @@ subroutine rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_f enddo ! ! If upper-air RH is missing, see if we can interpolate from surrounding levels. -! This is a simple vertical interpolation, linear in pressure. -! Currently, this simply fills in one missing level between two present levels. -! May expand this in the future to fill in additional levels. May also expand -! this in the future to vertically interpolate other variables. +! This is a simple vertical interpolation, linear or log in pressure. +! Currently, this simply fills in missing levels between two present levels. + ! do k = 2, nlvl-1, 1 if (plvl(k-1) .lt. 200000.) then if ( (.not. is_there(nint(plvl(k)),'RH')) .and. & - ( is_there(nint(plvl(k-1)), 'RH')) .and.& - ( is_there(nint(plvl(k+1)), 'RH')) ) then - call get_dims(nint(plvl(k+1)), 'RH') - call vntrp(plvl, maxlvl, k, "RH ", map%nx, map%ny) + ( is_there(nint(plvl(k-1)), 'RH')) ) then + found_level = .false. + RH_loop : do itest = k+1,nlvl,1 + if ( ( is_there(nint(plvl(itest)), 'RH')) ) then + found_level = .true. + exit RH_loop + endif + enddo RH_loop + if( found_level ) then + call get_dims(nint(plvl(itest)), 'RH') + call vntrp(plvl, maxlvl, k, itest, interp_type, "RH ", map%nx, map%ny) + else + PRINT *,'WARNING: Could not interpolate RH to level k=',k,' p=',plvl(k),& + 'because could not find any level above this level.' + endif endif endif enddo @@ -930,28 +1175,40 @@ subroutine gfs_trop_maxw_pressures(ix,jx) end subroutine gfs_trop_maxw_pressures -subroutine vntrp(plvl, maxlvl, k, name, ix, jx) +subroutine vntrp(plvl, maxlvl, k, k2, interp_type, name, ix, jx) use storage_module + use module_debug implicit none - integer :: ix, jx, k, maxlvl + integer :: ix, jx, k, k2, maxlvl real, dimension(maxlvl) :: plvl character(len=8) :: name real, dimension(ix,jx) :: a, b, c real :: frc + integer :: interp_type - write(*,'("Interpolating to fill in ", A, " at level ", I8)') trim(name), nint(plvl(k)) + write(*,'("Interpolating to fill in ", A, " at level ", F8.2, " hPa using levels ", F8.2," hPa and ", & + F8.2," hPa")') trim(name), plvl(k)/100.0,plvl(k-1)/100.0, plvl(k2)/100.0 + call mprintf(.true.,INFORM, & + "RRPR: Interpolating to fill in %s at %i Pa using levels %i Pa and %i Pa",& + s1=trim(name), i1=int(plvl(k)), i2=int(plvl(k-1)), i3=int(plvl(k2))) call get_storage(nint(plvl(k-1)), name, a, ix, jx) - call get_storage(nint(plvl(k+1)), name, c, ix, jx) + call get_storage(nint(plvl(k2)), name, c, ix, jx) - frc = (plvl(k) - plvl(k+1)) / ( plvl(k-1)-plvl(k+1)) + if ( interp_type == 1 ) then + frc = log(plvl(k)/plvl(k2)) / log(plvl(k-1)/plvl(k2)) + else + frc = (plvl(k) - plvl(k2)) / (plvl(k-1)-plvl(k2)) + endif b = (1.-frc)*c + frc*a + !KWM b = 0.5 * (a + c) call put_storage(nint(plvl(k)), name, b, ix, jx) end subroutine vntrp + subroutine fix_gfs_rh (ix, jx, plvl) ! This routine replaces GFS RH (wrt ice) with RH wrt liquid (which is what is assumed in real.exe). use storage_module diff --git a/ungrib/src/ungrib.F b/ungrib/src/ungrib.F index d18360d..9dd7ed7 100644 --- a/ungrib/src/ungrib.F +++ b/ungrib/src/ungrib.F @@ -54,15 +54,37 @@ program ungrib implicit none + interface + subroutine read_namelist(hstart, hend, delta_time, ntimes,& + ordered_by_date, debug_level, out_format, prefix, & + add_lvls, new_plvl, interp_type) + + use misc_definitions_module + + character(len=19) :: hstart, hend + integer :: delta_time + integer :: ntimes + logical :: ordered_by_date + integer :: debug_level + character(len=3) :: out_format + character(len=MAX_FILENAME_LEN) :: prefix + logical :: add_lvls + real, dimension(:) :: new_plvl + integer :: interp_type + end subroutine read_namelist + end interface + integer :: nunit1 = 12 character(LEN=132) :: gribflnm = 'GRIBFILE.AAA ' ! won't work with len=12 integer :: debug_level - integer , parameter :: maxlvl = 150 + integer , parameter :: maxlvl = 250 - real , dimension(maxlvl) :: plvl + real , dimension(maxlvl) :: plvl, new_plvl integer :: iplvl + logical :: add_lvls + integer :: interp_type integer :: nlvl @@ -93,7 +115,8 @@ program ungrib ! Read the namelist, and return the information we want: call read_namelist(hstart, hend, interval, ntimes, & - ordered_by_date, debug_level, out_format, prefix) + ordered_by_date, debug_level, out_format, prefix, & + add_lvls, new_plvl, interp_type) call mprintf(.true.,INFORM,"Interval value: %i seconds or %f hours", & i1=interval, f1=float(interval)/3600.) @@ -306,9 +329,9 @@ program ungrib if ((hsave(1:4).ne.'0000').and.(hsave.le.hend)) then if (debug_level .gt. 100) then - print*, 'Calling output: '//hsave(1:13) - call mprintf(.true.,DEBUG,"Calling output: %s ",s1=hsave(1:13)) - endif + print*, 'Calling output: '//hsave(1:13) + call mprintf(.true.,DEBUG,"Calling output: %s ",s1=hsave(1:13)) + endif call output(hsave, nlvl, maxlvl, plvl, interval, 1, out_format, prefix, debug_level) hsave=hdate @@ -343,7 +366,9 @@ program ungrib call mprintf(.true.,INFORM,"First pass done, doing a reprocess") - call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, debug_level, out_format, prefix) + call rrpr(hstart, ntimes, interval, nlvl, maxlvl, plvl, & + add_lvls, new_plvl, interp_type, & + debug_level, out_format, prefix) ! Make sure the filedates are in order, with an inefficient sort: -- 2.11.4.GIT