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" ! ! ! !*****************************************************************************! ! ! use filelist use gridinfo use storage_module use table use module_debug use misc_definitions_module use stringutil implicit none !------------------------------------------------------------------------------ ! Arguments: ! HSTART: Starting date of times to process character (LEN=19) :: hstart ! NTIMES: Number of time periods to process integer :: ntimes ! INTERVAL: Time inteval (seconds) of time periods to process. integer :: interval ! NLVL: The number of levels in the stored data. integer :: nlvl ! MAXLVL: The parameterized maximum number of levels to allow. integer :: maxlvl ! 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 !------------------------------------------------------------------------------ character (LEN=25) :: units character (LEN=46) :: Desc real, allocatable, dimension(:,:) :: scr2d, tmp2d 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 character(LEN=24) :: hdate_output character(LEN=3) :: out_format character(LEN=MAX_FILENAME_LEN) :: prefix real :: xfcst, level character(LEN=9) :: field 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 ! Decide the length of date strings to use for output file names. ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds. if (mod(interval,3600) == 0) then datelen = 13 else if (mod(interval, 60) == 0) then datelen = 16 else datelen = 19 endif if ( debug_level .gt. 100 ) then call mprintf(.true.,DEBUG,"Begin rrpr") call mprintf(.true.,DEBUG,"nfiles = %i , ntimes = %i )",i1=nfiles,i2=ntimes) do n = 1, nfiles call mprintf(.true.,DEBUG,"filedates(%i) = %s",i1=n,s1=filedates(n)) enddo endif ! Compute the ending time: call geth_newdate(hend, hstart, interval*ntimes) call clear_storage ! We want to do something for each of the requested times: TIMELOOP : do ntime = 1, ntimes idts = (ntime-1) * interval call geth_newdate(hdate, hstart, idts) call mprintf(.true.,DEBUG, & "RRPR: hstart = %s , hdate = %s , idts = %i",s1=hstart,s2=hdate,i1=idts) ! Loop over the output file dates, and do stuff if the file date matches ! the requested time we are working on now. FILELOOP : do n = 1, nfiles if ( debug_level .gt. 100 ) then call mprintf(.true.,DEBUG, & "hstart = %s , hend = %s",s1=hstart,s2=hend) call mprintf(.true.,DEBUG, & "filedates(n) = %s",s1=filedates(n)) call mprintf(.true.,DEBUG, & "filedates(n) = %s",s1=filedates(n)(1:datelen)) end if if (filedates(n)(1:datelen).ne.hdate(1:datelen)) cycle FILELOOP if (debug_level .gt. 50 ) then call mprintf(.true.,INFORM, & "RRPR Processing : %s",s1=filedates(n)(1:datelen)) endif open(iunit, file=trim(get_path(prefix))//'PFILE:'//filedates(n)(1:datelen), & form='unformatted',status='old') ! Read the file: rdloop: do read (iunit, iostat=ierr) ifv if (ierr.ne.0) exit rdloop if ( ifv .eq. 5) then ! WPS read (iunit) hdate_output, xfcst, map%source, field, units, Desc, & level, map%nx, map%ny, map%igrid hdate = hdate_output(1:19) select case (map%igrid) case (0, 4) read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth case (3) read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%truelat2, map%r_earth case (5) read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%r_earth case (1) read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, & map%truelat1, map%r_earth case (6) read (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, & map%lat0,map%lon0, map%r_earth case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in RRPR 1",i1=map%igrid) end select read (iunit) map%grid_wind else if ( ifv .eq. 4 ) then ! SI read (iunit) hdate_output, xfcst, map%source, field, units, desc, level, & map%nx, map%ny, map%igrid hdate = hdate_output(1:19) select case (map%igrid) case (0, 4) read(iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx case (3) read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2 case (5) read (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in RRPR 2",i1=map%igrid) end select else if ( ifv .eq. 3 ) then ! MM5 read(iunit) hdate_output, xfcst, field, units, desc, level,& map%nx, map%ny, map%igrid hdate = hdate_output(1:19) select case (map%igrid) case (3) ! lamcon read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%truelat2 case (5) ! Polar Stereographic read (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1 case (0, 4) ! lat/lon read (iunit) map%lat1, map%lon1, map%dy, map%dx case (1) ! Mercator read (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in RRPR 3",i1=map%igrid) end select else call mprintf(.true.,ERROR, & "unknown out_format, ifv = %i",i1=ifv) endif allocate(ptr2d(map%nx,map%ny)) read (iunit) ptr2d call refw_storage(nint(level), field, ptr2d, map%nx, map%ny) nullify (ptr2d) enddo rdloop write (0,*) 'Name of source model =>',map%source ! ! We have reached the end of file, so time to close it. ! close(iunit) if (debug_level .gt. 100 ) call print_storage ! ! By now the file has been read completely. Now, see if we need to fill in ! missing fields: ! ! 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 ! do k = 1, nlvl if ((plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then ! We found a level between 200100 and 200200, now find the field ! corresponding to that level. MLOOP : do mm = 1, maxvar if (is_there(nint(plvl(k)), namvar(mm))) then INLOOP : do kk = 200101, nint(plvl(k)) if (is_there(kk, namvar(mm))) then if ( debug_level .gt. 100 ) then call mprintf(.true.,DEBUG, & "Copying %s at level %i to level 200100.",s1=namvar(mm),i1=kk) end if call get_dims(kk, namvar(mm)) allocate(scr2d(map%nx,map%ny)) call get_storage & (kk, namvar(mm), scr2d, map%nx, map%ny) call put_storage & (200100,namvar(mm), scr2d,map%nx,map%ny) deallocate(scr2d) EXIT INLOOP endif enddo INLOOP endif enddo MLOOP endif enddo ! ! If upper-air U is missing, see if we can interpolate from surrounding 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')) ) 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 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')) ) 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 ! ! If upper-air SPECHUMD is missing, see if we can compute SPECHUMD from QVAPOR: !--- Tanya's change for initializing WRF with RUC do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)), 'SPECHUMD').and. & is_there(nint(plvl(k)), 'QV')) then call get_dims(nint(plvl(k)), 'QV') call compute_spechumd_qvapor(map%nx, map%ny, plvl(k)) endif endif enddo !--- Tanya's change for initializing WRF with RUC ! This allows for the ingestion for RUC isentropic data ! do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)), 'TT').and. & is_there(nint(plvl(k)), 'VPTMP').and. & is_there(nint(plvl(k)), 'SPECHUMD')) then call get_dims(nint(plvl(k)), 'VPTMP') call compute_t_vptmp(map%nx, map%ny, plvl(k)) endif endif enddo !!! ! ! If upper-air T is missing, see if we can interpolate from surrounding 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')) ) 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. & is_there(nint(plvl(k)), 'GEOPT')) then call get_dims(nint(plvl(k)), 'GEOPT') allocate(scr2d(map%nx,map%ny)) call get_storage(nint(plvl(k)), 'GEOPT', scr2d, map%nx, map%ny) 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 ") 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, ! 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 ! neighbor AND a 4-pt technique. Those two pressures are used in ARW real ! for vertical interpolation of the trop and max wind level data. ! if (index(map%source,'NCEP GFS') .ne. 0 ) then call mprintf(.true.,DEBUG, & "RRPR: Replicating GFS pressures for max wind and trop") if ( is_there(200100,'PMAXW ') .or. & is_there(200100,'PTROP ') ) then call gfs_trop_maxw_pressures (map%nx, map%ny) endif endif ! Repair GFS and ECMWF pressure-level RH if (index(map%source,'NCEP GFS') .ne. 0 .or. & index(map%source,'NCEP CDAS CFSV2') .ne. 0 .or. & index(map%source,'ECMWF') .ne. 0 ) then call mprintf(.true.,DEBUG, & "RRPR: Adjusting RH values ") do k = 1, nlvl if ( is_there(nint(plvl(k)),'RH') .and. & is_there(nint(plvl(k)),'TT') ) then call fix_gfs_rh (map%nx, map%ny, plvl(k)) endif enddo endif ! If upper-air RH is missing, see if we can compute RH from Specific Humidity: do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)), 'RH') .and. & is_there(nint(plvl(k)), 'TT') .and. & is_there(nint(plvl(k)), 'SPECHUMD')) then call get_dims(nint(plvl(k)), 'TT') call compute_rh_spechumd_upa(map%nx, map%ny, plvl(k)) endif endif enddo ! If upper-air RH is missing, see if we can compute RH from Vapor Pressure: ! (Thanks to Bob Hart of PSU ESSC -- 1999-05-27.) do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)),'RH').and. & is_there(nint(plvl(k)), 'TT') .and. & is_there(nint(plvl(k)),'VAPP')) then call get_dims(nint(plvl(k)),'TT') call compute_rh_vapp_upa(map%nx, map%ny, plvl(k)) endif endif enddo ! If upper-air RH is missing, see if we can compute RH from Dewpoint Depression: do k = 1, nlvl if (plvl(k).lt.200000.) then if (.not. is_there(nint(plvl(k)),'RH').and. & is_there(nint(plvl(k)), 'TT') .and. & is_there(nint(plvl(k)),'DEPR')) then call get_dims(nint(plvl(k)),'TT') call compute_rh_depr(map%nx, map%ny, plvl(k)) endif endif enddo ! ! If upper-air RH is missing, see if we can interpolate from surrounding 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)),'RH')) .and. & ( 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 ! ! Check to see if we need to fill RH above 300 mb: ! if (is_there(30000, 'RH')) then call get_dims(30000, 'RH') allocate(scr2d(map%nx,map%ny)) do k = 1, nlvl ! Set missing RH to 5% between 300 and 70 hPa. Set RH to 0 above 70 hPa. ! The stratospheric RH will be adjusted further in real. if (plvl(k).le.7000.) then scr2d = 0. else if (plvl(k).lt.30000.) then scr2d = 5. endif if (plvl(k).lt.30000. .and. plvl(k) .gt. 10. ) then ! levels higher than .1 mb are special - do not fill if (.not. is_there(nint(plvl(k)), 'RH')) then call put_storage(nint(plvl(k)),'RH',scr2d,map%nx,map%ny) call mprintf(.true.,DEBUG, & "RRPR: RH missing at %i hPa, inserting synthetic RH ",i1=nint(plvl(k)/100.)) endif endif enddo deallocate(scr2d) endif ! ! If surface RH is missing, see if we can compute RH from Specific Humidity ! or Dewpoint or Dewpoint depression: ! if (.not. is_there (200100, 'RH')) then if (is_there(200100, 'TT').and. & is_there(200100, 'PSFC' ) .and. & is_there(200100, 'SPECHUMD')) then call get_dims(200100, 'TT') call compute_rh_spechumd(map%nx, map%ny) call mprintf(.true.,DEBUG, & "RRPR: SURFACE RH is computed") elseif (is_there(200100, 'TT' ).and. & is_there(200100, 'DEWPT')) then call get_dims(200100, 'TT') call compute_rh_dewpt(map%nx, map%ny) elseif (is_there(200100, 'TT').and. & is_there(200100, 'DEPR')) then call get_dims(200100, 'TT') call compute_rh_depr(map%nx, map%ny, 200100.) endif endif ! ! If surface SNOW is missing, see if we can compute SNOW from SNOWRUC ! (From Wei Wang, 2007 June 21, modified 12/28/2007) ! if (.not. is_there(200100, 'SNOW') .and. & is_there(200100, 'SNOWRUC')) then call get_dims(200100, 'SNOWRUC') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SNOWRUC', scr2d, map%nx, map%ny) scr2d = scr2d * 1000. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny) deallocate(scr2d) endif ! compute snow water equivalent (SNOW) for NCEP RUC models ! As of Sept. 14 2011 if ( index(map%source,'NCEP RUC Model') .ne. 0) then if (is_there(200100, 'SNOWH') .and. .not. is_there(200100, 'SNOW')) then call get_dims(200100, 'SNOWH') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SNOWH', scr2d, map%nx, map%ny) call mprintf(.true.,DEBUG, & "RRPR: Computing SNOWH from SNOW") if (is_there(200100, 'RHOSN')) then ! If we have snow density, use it to compute snowh call get_dims(200100, 'RHOSN') allocate(tmp2d(map%nx,map%ny)) call get_storage(200100, 'RHOSN', tmp2d, map%nx, map%ny) scr2d = scr2d * tmp2d deallocate(tmp2d) else scr2d = scr2d * 200.0 ! Assume 200:1 ratio endif call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny) deallocate(scr2d) endif endif ! Modify the 2017 GFS masked fields if (index(map%source,'NCEP GFS') .ne. 0 ) then call mprintf(.true.,DEBUG, & "RRPR: Adjusting GFS masked fields ") if ( is_there(200100, 'ST000010')) then call get_dims(200100, 'ST000010') call fix_gfs_miss (map%nx, map%ny, 200100.) endif endif ! Add residual soil moisture to SOILM* if initialized from the GSD RUC model or from NCEP RUC if (index(map%source,'NOAA GSD') .ne. 0 .or. & index(map%source,'NCEP RUC Model') .ne. 0) then if ( .not. is_there(200100, 'SOILM000') .and.& is_there(200100, 'SM000ruc') ) then call get_dims(200100, 'SM000ruc') print *,'Adjust RUC soil moisture' call mprintf(.true.,DEBUG, & "RRPR: Adjusting RUC soil moisture ") call fix_ruc_soilm (map%nx, map%ny) endif endif ! ! Check to see if we need to fill SOILHGT from SOILGEO. ! (From Wei Wang, 2007 June 21) ! if (.not. is_there(200100, 'SOILHGT') .and. & is_there(200100, 'SOILGEO')) then call get_dims(200100, 'SOILGEO') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SOILGEO', scr2d, map%nx, map%ny) scr2d = scr2d / 9.81 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny) call mprintf(.true.,DEBUG, & "RRPR: Computing SOILGHT from SOILGEO ") deallocate(scr2d) endif ! For hybrid-level input, soilgeo is in level 1 (e.g. ERA40) if (.not. is_there(200100, 'SOILHGT') .and. & is_there(1, 'SOILGEO')) then call get_dims(1, 'SOILGEO') allocate(scr2d(map%nx,map%ny)) call get_storage(1, 'SOILGEO', scr2d, map%nx, map%ny) scr2d = scr2d / 9.81 call put_storage(200100, 'SOILHGT', scr2d, map%nx, map%ny) deallocate(scr2d) endif ! For NCEP RR (using the same ID as for RUC) native-level input, ! may need to move PSFC from level 1 to 2001. ! From TGS 8 Sept. 2011 if ( index(map%source,'NCEP RUC Model') .ne. 0) then if (.not. is_there(200100, 'PSFC') .and. & is_there(1, 'PRESSURE')) then print *,'Process PSFC for NCEP RR' call get_dims(1, 'PRESSURE') allocate(scr2d(map%nx,map%ny)) call get_storage(1, 'PRESSURE', scr2d, map%nx, map%ny) call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny) deallocate(scr2d) endif endif ! For ECMWF hybrid-level input, may need to move psfc from level 1 to 2001. if ( index(map%source,'ECMWF') .ne. 0) then if (.not. is_there(200100, 'PSFC') .and. & is_there(1, 'PSFCH')) then call get_dims(1, 'PSFCH') allocate(scr2d(map%nx,map%ny)) call get_storage(1, 'PSFCH', scr2d, map%nx, map%ny) call put_storage(200100, 'PSFC', scr2d, map%nx, map%ny) deallocate(scr2d) endif endif ! ECMWF snow depth in meters of water equivalent (Table 128). Convert to kg/m2 ! if (is_there(200100, 'SNOW_EC')) then call get_dims(200100, 'SNOW_EC') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SNOW_EC', scr2d, map%nx, map%ny) scr2d = scr2d * 1000. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny) deallocate(scr2d) endif ! Convert the ECMWF LANDSEA mask from a fraction to a flag if ( index(map%source,'ECMWF') .ne. 0) then if (is_there(200100, 'LANDSEA')) then call get_dims(200100, 'LANDSEA') call make_zero_or_one(map%nx, map%ny, 'LANDSEA') endif endif ! NCEP GFS weasd is one-half of the NAM value. Increase it for use in WRF. ! The GFS-based reanalyses values should be OK as is. if ((index(map%source,'NCEP GFS') .ne. 0 .or. & index(map%source,'NCEP GEFS') .ne. 0) .and. & is_there(200100, 'SNOW')) then call mprintf(.true.,DEBUG, & "RRPR: Recomputing SNOW for NCEP GFS") call get_dims(200100, 'SNOW') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny) scr2d = scr2d * 2. call put_storage(200100, 'SNOW', scr2d, map%nx, map%ny) deallocate(scr2d) endif ! compute physical snow depth (SNOWH) for various models ! As of March 2011, this is done here instead of real because we have model ! source information. if (is_there(200100, 'SNOW') .and. .not. is_there(200100, 'SNOWH')) then call get_dims(200100, 'SNOW') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'SNOW', scr2d, map%nx, map%ny) call mprintf(.true.,DEBUG, & "RRPR: Computing SNOWH from SNOW") if ( index(map%source,'NCEP ') .ne. 0) then scr2d = scr2d * 0.005 ! Assume 200:1 ratio as used at NCEP and in NOAH else if (index(map%source,'ECMWF') .ne. 0) then if (is_there(200100, 'SNOW_DEN')) then ! If we have snow density, use it to compute snowh call get_dims(200100, 'SNOW_DEN') allocate(tmp2d(map%nx,map%ny)) call get_storage(200100, 'SNOW_DEN', tmp2d, map%nx, map%ny) scr2d = scr2d / tmp2d deallocate(tmp2d) else scr2d = scr2d * 0.004 ! otherwise, assume a density of 250 mm/m (i.e. 250:1 ratio). endif else ! Other models scr2d = scr2d * 0.005 ! Use real's default method (200:1) endif call put_storage(200100, 'SNOWH', scr2d, map%nx, map%ny) deallocate(scr2d) endif ! As of March 2011, SEAICE can be a flag or a fraction. It will be converted ! to the appropriate values in real depending on whether or not the polar mods are used. !! If we've got a SEAICE field, make sure that it is all Zeros and Ones: ! if (is_there(200100, 'SEAICE')) then ! call get_dims(200100, 'SEAICE') ! call make_zero_or_one(map%nx, map%ny, 'SEAICE') ! endif ! If we've got an ICEMASK field, re-flag it for output to met_em and real: ! Field | GRIB In | Out ! ------------------------- ! water | 0 | 0 ! land | -1 | 1 ! ice | 1 | 0 if (is_there(200100, 'ICEMASK')) then call get_dims(200100, 'ICEMASK') call re_flag_ice_mask(map%nx, map%ny) endif ! If we have an ICEFRAC field, convert from % to fraction if (is_there(200100, 'ICEFRAC')) then call get_dims(200100, 'ICEFRAC') allocate(scr2d(map%nx,map%ny)) call get_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny) scr2d = scr2d / 100. call put_storage(200100, 'ICEFRAC', scr2d, map%nx, map%ny) deallocate(scr2d) endif call mprintf(.true.,INFORM, & "RRPR: hdate = %s ",s1=hdate) call output(hdate, nlvl, maxlvl, plvl, interval, 2, out_format, prefix, debug_level) call clear_storage exit FILELOOP enddo FILELOOP enddo TIMELOOP end subroutine rrpr subroutine make_zero_or_one(ix, jx, infield) ! Make sure the input field (SEAICE or LANDSEA) is zero or one. use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: seaice character(len=*) :: infield call get_storage(200100, infield, seaice, ix, jx) where(seaice > 0.5) seaice = 1.0 elsewhere seaice = 0.0 end where call put_storage(200100, infield, seaice, ix, jx) end subroutine make_zero_or_one subroutine re_flag_ice_mask(ix, jx) ! ! Change land points from -1 to 1 ! Change ice points from 1 to 0 ! Water points stay at 0 ! use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: iceflag call get_storage(200100, 'ICEMASK',iceflag, ix, jx) where(iceflag > 0.5) ! Ice points, set to water value iceflag = 0.0 end where where(iceflag < -0.5) ! Land points iceflag = 1.0 end where call put_storage(200100, 'ICEMASK',iceflag, ix, jx) end subroutine re_flag_ice_mask subroutine compute_spechumd_qvapor(ix, jx, plvl) ! Compute specific humidity from water vapor mixing ratio. use storage_module implicit none integer :: ix, jx real :: plvl real, dimension(ix,jx) :: QVAPOR, SPECHUMD call get_storage(nint(plvl), 'QV', QVAPOR, ix, jx) SPECHUMD = QVAPOR/(1.+QVAPOR) call put_storage(nint(plvl), 'SPECHUMD', spechumd, ix, jx) if(nint(plvl).eq.1) then call put_storage(200100,'SPECHUMD', spechumd, ix, jx) endif end subroutine compute_spechumd_qvapor subroutine compute_t_vptmp(ix, jx, plvl) ! Compute temperature from virtual potential temperature use storage_module implicit none integer :: ix, jx real :: plvl real, dimension(ix,jx) :: T, VPTMP, P, Q real, parameter :: rovcp=0.28571 call get_storage(nint(plvl), 'VPTMP', VPTMP, ix, jx) IF (nint(plvl) .LT. 200) THEN call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) ELSE p = plvl ENDIF call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx) t=vptmp * (p*1.e-5)**rovcp * (1./(1.+0.6078*Q)) call put_storage(nint(plvl), 'TT', t, ix, jx) if(nint(plvl).eq.1) then call put_storage(200100, 'PSFC', p, ix, jx) endif end subroutine compute_t_vptmp subroutine compute_rh_spechumd(ix, jx) ! Compute relative humidity from specific humidity. use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: T, P, RH, Q real, parameter :: svp1=611.2 real, parameter :: svp2=17.67 real, parameter :: svp3=29.65 real, parameter :: svpt0=273.15 real, parameter :: eps = 0.622 call get_storage(200100, 'TT', T, ix, jx) call get_storage(200100, 'PSFC', P, ix, jx) call get_storage(200100, 'SPECHUMD', Q, ix, jx) rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3))) call put_storage(200100, 'RH', rh, ix, jx) end subroutine compute_rh_spechumd subroutine compute_rh_spechumd_upa(ix, jx, plvl) ! Compute relative humidity from specific humidity. use storage_module implicit none integer :: ix, jx real :: plvl real, dimension(ix,jx) :: T, P, RH, Q real, parameter :: svp1=611.2 real, parameter :: svp2=17.67 real, parameter :: svp3=29.65 real, parameter :: svpt0=273.15 real, parameter :: eps = 0.622 IF ( nint(plvl).LT. 200) THEN if (is_there(nint(plvl), 'PRESSURE')) then call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) else return ! if we don't have pressure on model levels, return endif ELSE P = plvl ENDIF call get_storage(nint(plvl), 'TT', T, ix, jx) call get_storage(nint(plvl), 'SPECHUMD', Q, ix, jx) Q=MAX(1.E-10,Q) rh = 1.E2 * (p*q/(q*(1.-eps) + eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3))) call put_storage(nint(plvl), 'RH', rh, ix, jx) end subroutine compute_rh_spechumd_upa subroutine compute_rh_vapp_upa(ix, jx, plvl) ! Compute relative humidity from vapor pressure. ! Thanks to Bob Hart of PSU ESSC -- 1999-05-27. use storage_module implicit none integer :: ix, jx real :: plvl real, dimension(ix,jx) :: P, ES real, pointer, dimension(:,:) :: T, E, RH real, parameter :: svp1=611.2 real, parameter :: svp2=17.67 real, parameter :: svp3=29.65 real, parameter :: svpt0=273.15 allocate(RH(ix,jx)) IF ( nint(plvl).LT. 200) THEN if (is_there(nint(plvl), 'PRESSURE')) then call get_storage(nint(plvl), 'PRESSURE', P, ix, jx) else return ! if we don't have pressure on model levels, return endif ELSE P = plvl ENDIF call refr_storage(nint(plvl), 'TT', T, ix, jx) call refr_storage(nint(plvl), 'VAPP', E, ix, jx) ES=svp1*exp(svp2*(T-svpt0)/(T-svp3)) rh=min(1.E2*(P-ES)*E/((P-E)*ES), 1.E2) call refw_storage(nint(plvl), 'RH', rh, ix, jx) nullify(T,E) end subroutine compute_rh_vapp_upa subroutine compute_rh_depr(ix, jx, plvl) ! Compute relative humidity from Dewpoint Depression use storage_module implicit none integer :: ix, jx real :: plvl real, dimension(ix,jx) :: t, depr, rh real, parameter :: Xlv = 2.5e6 real, parameter :: Rv = 461.5 integer :: i, j call get_storage(nint(plvl), 'TT', T, ix, jx) call get_storage(nint(plvl), 'DEPR', DEPR, ix, jx) where(DEPR < 100.) rh = exp(Xlv/Rv*(1./T - 1./(T-depr))) * 1.E2 elsewhere rh = 0.0 endwhere call put_storage(nint(plvl),'RH ', rh, ix, jx) end subroutine compute_rh_depr subroutine compute_rh_dewpt(ix,jx) ! Compute relative humidity from Dewpoint use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: t, dp, rh real, parameter :: Xlv = 2.5e6 real, parameter :: Rv = 461.5 call get_storage(200100, 'TT ', T, ix, jx) call get_storage(200100, 'DEWPT ', DP, ix, jx) rh = exp(Xlv/Rv*(1./T - 1./dp)) * 1.E2 call put_storage(200100,'RH ', rh, ix, jx) end subroutine compute_rh_dewpt subroutine gfs_trop_maxw_pressures(ix,jx) ! These are duplicate pressure values from the GFS, for ! the level of max wind speed and for the trop level. ! The duplicates are saved with a different name, so that ! the metgrid program can horizontally interpolate them to ! the model domain with a nearest neighbor method. use storage_module implicit none integer :: ix, jx real, dimension(ix,jx) :: pmaxw, pmaxwnn, ptrop, ptropnn if ( is_there(200100, 'PMAXW ') ) then call get_storage(200100, 'PMAXW ', pmaxw , ix, jx) pmaxwnn = pmaxw call put_storage(200100, 'PMAXWNN ', pmaxwnn, ix, jx) end if if ( is_there(200100, 'PTROP ') ) then call get_storage(200100, 'PTROP ', ptrop , ix, jx) ptropnn = ptrop call put_storage(200100, 'PTROPNN ', ptropnn, ix, jx) end if end subroutine gfs_trop_maxw_pressures subroutine vntrp(plvl, maxlvl, k, k2, interp_type, name, ix, jx) use storage_module use module_debug implicit none 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 ", 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(k2)), name, c, ix, jx) 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_miss (ix, jx, plvl) ! This routine replaces July 2017 GFS missing values with the WPS one. ! Earlier GFS files are unmodified. ! As of 2017, NCEP changed 'ocean' values for masked fields (ST, SM, SNOWH, etc.) ! from 0 to something other than missing. We will assume any large value ! (greater than 10^^18) is a missing code. ! We reset it to the WPS missing value (as set in METGRID.TBL) ! Changes described in http://www.nws.noaa.gov/os/notification/scn17-67gfsupgradeaaa.htm ! plvl must always be 200100. ! While, technically, any 3-d field could have a missing value in it we only deal with ! the surface fields which are known to have missing values. ! use storage_module implicit none integer :: ix, jx, i, j, k real :: plvl real, allocatable, dimension(:,:) :: f, sea integer, parameter :: nvar = 10 character, dimension(nvar) :: flist*8 data flist/'ST000010','ST010040','ST040100','ST100200','ST010200', & 'SM000010','SM010040','SM040100','SM100200','SM010200' / allocate(sea(ix,jx)) allocate(f(ix,jx)) ! If LANDN is present (July 2017 and later GFS output), use it for LANDSEA if ( is_there(200100, 'LANDN') ) then call get_storage(200100, 'LANDN', sea, ix, jx) call put_storage(200100, 'LANDSEA', sea, ix, jx) endif do k = 1, nvar if (is_there(200100, flist(k) )) then call get_storage(nint(plvl), flist(k), f, ix, jx) do j = 1, jx do i = 1, ix if (abs(f(i,j)) .gt. 1.e18) then f(i,j) = -1.e30 endif enddo enddo ! Limit soil moisture to .468 (should only occur for permanent land ice points ! according to NCEP documentation.) if (flist(k)(1:2) .eq. 'SM' ) then do j = 1, jx do i = 1, ix if ((f(i,j)) .gt. 0.468) then f(i,j) = 0.468 endif enddo enddo endif call put_storage(200100, flist(k), f, ix, jx) endif enddo ! Snow fields have a different WPS missing value and must be adjusted ! separately from the soil fields (unless we want to have an array of missing values,too). if (is_there(200100, 'SNOW' )) then call get_storage(200100, 'SNOW', f, ix, jx) do j = 1, jx do i = 1, ix if (abs(f(i,j)) .gt. 1.e18) then f(i,j) = 0 endif enddo enddo call put_storage(200100, 'SNOW', f, ix, jx) endif if (is_there(200100, 'SNOWH' )) then call get_storage(200100, 'SNOWH', f, ix, jx) do j = 1, jx do i = 1, ix if (abs(f(i,j)) .gt. 1.e18) then f(i,j) = 0 endif enddo enddo call put_storage(200100, 'SNOWH', f, ix, jx) endif deallocate (f) deallocate (sea) end subroutine fix_gfs_miss 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 implicit none integer :: ix, jx, i, j real :: plvl, eis, ews, r real, allocatable, dimension(:,:) :: rh, tt allocate(rh(ix,jx)) allocate(tt(ix,jx)) call get_storage(nint(plvl), 'RH', rh, ix, jx) call get_storage(nint(plvl), 'TT', tt, ix, jx) do j = 1, jx do i = 1, ix if ( tt(i,j) .le. 273.15 ) then ! Murphy and Koop 2005 ice saturation vapor pressure. ! eis and ews in hPA, tt is in K eis = .01 * exp (9.550426 - (5723.265 / tt(i,j)) + (3.53068 * alog(tt(i,j))) & - (0.00728332 * tt(i,j))) ! Bolton 1980 liquid saturation vapor pressure. For water saturation, most ! formulae are very similar from 0 to -20, so we don't need a more exact formula. ews = 6.112 * exp(17.67 * (tt(i,j)-273.15) / ((tt(i,j)-273.15)+243.5)) if ( tt(i,j) .gt. 253.15 ) then ! A linear approximation to the GFS blending region ( -20 > T < 0 ) r = ((273.15 - tt(i,j)) / 20.) r = (r * eis) + ((1-r)*ews) else r = eis endif rh(i,j) = rh(i,j) * (r / ews) endif enddo enddo call put_storage(nint(plvl), 'RH', rh, ix, jx) deallocate (rh) deallocate (tt) end subroutine fix_gfs_rh subroutine fix_ruc_soilm (ix, jx) ! This routine adds residual soil moisture if initialized fron RUC use storage_module implicit none integer :: ix, jx, i, j REAL , DIMENSION(100) :: lqmi real, allocatable, dimension(:,:) :: soilm000, soilm005, soilm020, & soilm040, soilm160, soilm300,soilcat allocate(soilm000(ix,jx)) allocate(soilm005(ix,jx)) allocate(soilm020(ix,jx)) allocate(soilm040(ix,jx)) allocate(soilm160(ix,jx)) allocate(soilm300(ix,jx)) allocate(soilcat(ix,jx)) call get_storage(200100, 'SM000ruc', soilm000, ix, jx) call get_storage(200100, 'SM005ruc', soilm005, ix, jx) call get_storage(200100, 'SM020ruc', soilm020, ix, jx) call get_storage(200100, 'SM040ruc', soilm040, ix, jx) call get_storage(200100, 'SM160ruc', soilm160, ix, jx) call get_storage(200100, 'SM300ruc', soilm300, ix, jx) call get_storage(200100, 'SOILCAT', soilcat, ix, jx) lqmi(1:16) = & (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, & 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, & 0.004, 0.065 /) do j = 1, jx do i = 1, ix SOILM000(i,j)=SOILM000(i,j) + lqmi(nint(soilcat(i,j))) SOILM005(i,j)=SOILM005(i,j) + lqmi(nint(soilcat(i,j))) SOILM020(i,j)=SOILM020(i,j) + lqmi(nint(soilcat(i,j))) SOILM040(i,j)=SOILM040(i,j) + lqmi(nint(soilcat(i,j))) SOILM160(i,j)=SOILM160(i,j) + lqmi(nint(soilcat(i,j))) SOILM300(i,j)=SOILM300(i,j) + lqmi(nint(soilcat(i,j))) enddo enddo call put_storage(200100, 'SOILM000', soilm000, ix, jx) call put_storage(200100, 'SOILM005', soilm005, ix, jx) call put_storage(200100, 'SOILM020', soilm020, ix, jx) call put_storage(200100, 'SOILM040', soilm040, ix, jx) call put_storage(200100, 'SOILM160', soilm160, ix, jx) call put_storage(200100, 'SOILM300', soilm300, ix, jx) print *,'fix_ruc_soilm is done!' deallocate(soilm000) deallocate(soilm005) deallocate(soilm020) deallocate(soilm040) deallocate(soilm160) deallocate(soilm300) deallocate(soilcat) end subroutine fix_ruc_soilm