!WRF:MEDIATION_LAYER:IO
! ---
! This obs-nudging FDDA module (RTFDDA) is developed by the
! NCAR/RAL/NSAP (National Security Application Programs), under the
! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
! acknowledged for releasing this capability for WRF community
! research applications.
!
! The NCAR/RAL RTFDDA module was adapted, and significantly modified
! from the obs-nudging module in the standard MM5V3.1 which was originally
! developed by PSU (Stauffer and Seaman, 1994).
!
! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
! Nov. 2006
!
! References:
!
! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
! implementation of obs-nudging-based FDDA into WRF for supporting
! ATEC test operations. 2005 WRF user workshop. Paper 10.7.
!
! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
! on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
! and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
!
! Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
! assimilation. J. Appl. Meteor., 33, 416-434.
!
! http://www.rap.ucar.edu/projects/armyrange/references.html
!
SUBROUTINE wrf_fddaobs_in (grid ,config_flags) 1,5
USE module_domain
USE module_configure
USE module_model_constants !rovg
IMPLICIT NONE
TYPE(domain) :: grid
TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
#if ( EM_CORE == 1 )
! Local variables
integer :: ktau ! timestep index corresponding to xtime
integer :: krest ! restart timestep
integer :: inest ! nest level
integer :: infreq ! input frequency
integer :: nstlev ! nest level
real :: dtmin ! dt in minutes
real :: xtime ! forecast time in minutes
logical :: iprt_in4dob ! print flag
INTEGER ids , ide , jds , jde , kds , kde , &
ims , ime , jms , jme , kms , kme , &
ips , ipe , jps , jpe , kps , kpe
INTEGER ij, its, ite, jts, jte
! Modified to also call in4dob intially, since subr in4dob is no
! longer called from subr fddaobs_init. Note that itimestep is now
! the model step BEFORE the model integration step, because this
! routine is now called by med_before_solve_io.
ktau = grid%itimestep ! ktau corresponds to xtime
krest = grid%fdob%ktaur
inest = grid%grid_id
nstlev = grid%fdob%levidn(inest)
infreq = grid%obs_ionf*(grid%parent_grid_ratio**nstlev)
iprt_in4dob = grid%obs_ipf_in4dob
IF( ((ktau.GT.krest.AND.MOD(ktau,infreq).EQ.0) &
.OR.(ktau.EQ.krest)) .AND. grid%xtime <= grid%fdda_end ) then
! Calculate forecast time.
dtmin = grid%dt/60.
xtime = grid%xtime
CALL get_ijk_from_grid ( grid , &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe )
!$OMP PARALLEL DO &
!$OMP PRIVATE ( ij )
DO ij = 1 , grid%num_tiles
its = grid%i_start(ij)
ite = min(grid%i_end(ij),ide-1)
jts = grid%j_start(ij)
jte = min(grid%j_end(ij),jde-1)
CALL in4dob(inest, xtime, ktau, krest, dtmin, &
grid%julyr, grid%julday, grid%gmt, & !obsnypatch
grid%obs_nudge_opt, grid%obs_nudge_wind, grid%obs_nudge_temp, &
grid%obs_nudge_mois, grid%obs_nudge_pstr, grid%obs_coef_wind, &
grid%obs_coef_temp, grid%obs_coef_mois, grid%obs_coef_pstr, &
grid%obs_rinxy, grid%obs_rinsig, grid%fdob%window, &
grid%obs_npfi, grid%obs_ionf, &
grid%obs_prt_max, grid%obs_prt_freq, &
grid%obs_idynin, &
grid%obs_dtramp, grid%fdob, grid%fdob%varobs, &
grid%fdob%timeob, grid%fdob%nlevs_ob, grid%fdob%lev_in_ob, &
grid%fdob%plfo, grid%fdob%elevob, grid%fdob%rio, &
grid%fdob%rjo, grid%fdob%rko, &
grid%xlat, grid%xlong, &
config_flags%cen_lat, &
config_flags%cen_lon, &
config_flags%stand_lon, &
config_flags%truelat1, config_flags%truelat2, &
grid%fdob%known_lat, grid%fdob%known_lon, &
config_flags%dx, config_flags%dy, rovg, t0, &
grid%fdob%obsprt, &
grid%fdob%latprt, grid%fdob%lonprt, &
grid%fdob%mlatprt, grid%fdob%mlonprt, &
grid%fdob%stnidprt, &
ide, jde, &
ims, ime, jms, jme, &
its, ite, jts, jte, &
config_flags%map_proj, &
model_config_rec%parent_grid_ratio, &
model_config_rec%i_parent_start(inest), &
model_config_rec%j_parent_start(inest), &
model_config_rec%max_dom, &
model_config_rec%nobs_ndg_vars, grid%max_obs, iprt_in4dob)
ENDDO
!$OMP END PARALLEL DO
ENDIF
RETURN
#endif
END SUBROUTINE wrf_fddaobs_in
#if ( EM_CORE == 1 )
!------------------------------------------------------------------------------
! Begin subroutine in4dob and its subroutines
!------------------------------------------------------------------------------
SUBROUTINE in4dob(inest, xtime, ktau, ktaur, dtmin, & 1,68
myear, julday, gmt, & !obsnypatch
nudge_opt, iswind, istemp, &
ismois, ispstr, giv, &
git, giq, gip, &
rinxy, rinsig, twindo, &
npfi, ionf, prt_max, prt_freq, idynin, &
dtramp, fdob, varobs, &
timeob, nlevs_ob, lev_in_ob, &
plfo, elevob, rio, &
rjo, rko, &
xlat, xlong, &
cen_lat, &
cen_lon, &
stand_lon, &
true_lat1, true_lat2, &
known_lat, known_lon, &
dxm, dym, rovg, t0, &
obs_prt, &
lat_prt, lon_prt, &
mlat_prt, mlon_prt, &
stnid_prt, &
e_we, e_sn, &
ims, ime, jms, jme, &
its, ite, jts, jte, &
map_proj, &
parent_grid_ratio, &
i_parent_start, &
j_parent_start, &
maxdom, &
nndgv, niobf, iprt)
USE module_domain
USE module_model_constants, ONLY : rcp
USE module_date_time , ONLY : geth_idts
USE module_llxy
IMPLICIT NONE
! THIS IS SUBROUTINE READS AN OBSERVATION DATA FILE AND
! SELECTS ONLY THOSE VALUES OBSERVED AT TIMES THAT FALL
! WITHIN A TIME WINDOW (TWINDO) CENTERED ABOUT THE CURRENT
! FORECAST TIME (XTIME). THE INCOMING OBS FILES MUST BE
! IN CHRONOLOGICAL ORDER.
!
! NOTE: This routine was originally designed for MM5, which uses
! a nonstandard (I,J) coordinate system. For WRF, I is the
! east-west running coordinate, and J is the south-north
! running coordinate. So "J-slab" here is west-east in
! extent, not south-north as for MM5. RIO and RJO have
! the opposite orientation here as for MM5. -ajb 06/10/2004
! NOTE - IN4DOB IS CALLED ONLY FOR THE COARSE MESH TIMES
! - VAROBS(IVAR,N) HOLDS THE OBSERVATIONS.
! IVAR=1 OBS U
! IVAR=2 OBS V
! IVAR=3 OBS T
! IVAR=4 OBS Q
! IVAR=5 OBS Pressure
! IVAR=6 OBS Height
INTEGER, intent(in) :: niobf ! maximum number of observations
INTEGER, intent(in) :: nndgv ! number of nudge variables
INTEGER, intent(in) :: INEST ! nest level
REAL, intent(in) :: xtime ! model time in minutes
INTEGER, intent(in) :: KTAU ! current timestep
INTEGER, intent(in) :: KTAUR ! restart timestep
REAL, intent(in) :: dtmin ! dt in minutes
INTEGER, intent(in) :: myear ! model year !obsnypatch
INTEGER, intent(in) :: julday ! Julian day
REAL, intent(in) :: gmt ! Model GMT at start of run
INTEGER, intent(in) :: nudge_opt ! obs-nudge flag for this nest
INTEGER, intent(in) :: iswind ! nudge flag for wind
INTEGER, intent(in) :: istemp ! nudge flag for temperature
INTEGER, intent(in) :: ismois ! nudge flag for moisture
INTEGER, intent(in) :: ispstr ! nudge flag for pressure (obsolete)
REAL, intent(in) :: giv ! coefficient for wind
REAL, intent(in) :: git ! coefficient for temperature
REAL, intent(in) :: giq ! coefficient for moisture
REAL, intent(in) :: gip ! coefficient for pressure
REAL, intent(in) :: rinxy ! horizontal radius of influence (km)
REAL, intent(in) :: rinsig ! vertical radius of influence (on sigma)
REAL, intent(inout) :: twindo ! (time window)/2 (min) for nudging
INTEGER, intent(in) :: npfi ! coarse-grid time-step frequency for diagnostics
INTEGER, intent(in) :: ionf ! coarse-grid time-step frequency for obs-nudging calcs
INTEGER, intent(in) :: prt_max ! max number of entries of obs for diagnostic printout
INTEGER, intent(in) :: prt_freq ! frequency (in obs index) for diagnostic printout.
INTEGER, intent(in) :: idynin ! for dynamic initialization using a ramp-down function
REAL, intent(in) :: dtramp ! time period in minutes for ramping
TYPE(fdob_type), intent(inout) :: fdob ! derived data type for obs data
REAL, intent(inout) :: varobs(nndgv,niobf) ! observational values in each variable
REAL, intent(inout) :: timeob(niobf) ! model times for each observation (hours)
REAL, intent(inout) :: nlevs_ob(niobf) ! numbers of levels in sounding obs
REAL, intent(inout) :: lev_in_ob(niobf) ! level in sounding-type obs
REAL, intent(inout) :: plfo(niobf) ! index for type of obs-platform
REAL, intent(inout) :: elevob(niobf) ! elevations of observations (meters)
REAL, intent(inout) :: rio(niobf) ! west-east grid coordinate (non-staggered grid)
REAL, intent(inout) :: rjo(niobf) ! south-north grid coordinate (non-staggered grid)
REAL, intent(inout) :: rko(niobf) ! vertical grid coordinate
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: xlat, xlong ! lat/lon on mass-pt grid (for diagnostics only)
REAL, intent(in) :: cen_lat ! center latitude for map projection
REAL, intent(in) :: cen_lon ! center longitude for map projection
REAL, intent(in) :: stand_lon ! standard longitude for map projection
REAL, intent(in) :: true_lat1 ! truelat1 for map projection
REAL, intent(in) :: true_lat2 ! truelat2 for map projection
REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
REAL, intent(in) :: dxm ! grid size in x (meters)
REAL, intent(in) :: dym ! grid size in y (meters)
REAL, intent(in) :: rovg ! constant rho over g
REAL, intent(in) :: t0 ! background temperature
INTEGER, intent(inout) :: obs_prt(prt_max) ! For printout of obs index
REAL, intent(inout) :: lat_prt(prt_max) ! For printout of obs latitude
REAL, intent(inout) :: lon_prt(prt_max) ! For printout of obs longitude
REAL, intent(inout) :: mlat_prt(prt_max) ! For printout of model lat at obs (ri,rj)
REAL, intent(inout) :: mlon_prt(prt_max) ! For printout of model lon at obs (ri,rj)
INTEGER, intent(inout) :: stnid_prt(40,prt_max) ! For printout of model lon at obs (ri,rj)
INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
INTEGER, intent(in) :: ims ! grid memory start index (west-east dim)
INTEGER, intent(in) :: ime ! grid memory end index (west-east dim)
INTEGER, intent(in) :: jms ! grid memory start index (south-north dim)
INTEGER, intent(in) :: jme ! grid memory end index (south-north dim)
INTEGER, intent(in) :: its ! grid tile start index (west-east dim)
INTEGER, intent(in) :: ite ! grid tile end index (west-east dim)
INTEGER, intent(in) :: jts ! grid tile start index (south-north dim)
INTEGER, intent(in) :: jte ! grid tile end index (south-north dim)
INTEGER, intent(in) :: map_proj ! map projection index
INTEGER, intent(in) :: parent_grid_ratio ! parent to nest grid ration
INTEGER, intent(in) :: i_parent_start ! starting i coordinate in parent domain
INTEGER, intent(in) :: j_parent_start ! starting j coordinate in parent domain
INTEGER, intent(in) :: maxdom ! maximum number of domains
LOGICAL, intent(in) :: iprt ! print flag
!*** DECLARATIONS FOR IMPLICIT NONE
integer :: n, ndum, nopen, nvol, idate, imm, iss
integer :: nlast ! last obs in list of valid obs from prev window
integer :: nsta ! number of stations held in timeobs array
integer :: nstaw ! # stations within the actual time window
integer :: nprev ! previous n in obs read loop (for printout only)
integer :: meas_count, imc, njend, njc, njcc, julob, kn
real :: hourob, rjulob
real :: xhour, tback, tforwd, rjdate1, timanl1, rtimob
real :: rj, ri, elevation, pressure_data
real :: pressure_qc, height_data, height_qc, temperature_data
real :: temperature_qc, u_met_data, u_met_qc, v_met_data
real :: v_met_qc, rh_data, rh_qc, r_data, slp_data, slp_qc
real :: ref_pres_data, ref_pres_qc, psfc_data, psfc_qc
real :: precip_data, precip_qc, tbar, twdop
real*8 :: tempob
INTEGER, EXTERNAL :: nvals_le_limit ! for finding #obs with timeobs <= tforwd
! Local variables
TYPE (PROJ_INFO) :: obs_proj ! Structure for obs projection information.
character*14 date_char
character*19 obs_date !obsnypatch
integer idts !obsnypatch
character*40 platform,source,id,namef
character*2 fonc
character(len=200) :: msg ! Argument to wrf_message
real latitude,longitude
logical :: newpass ! New pass flag (used for printout)
logical is_sound,bogus
LOGICAL OPENED,exist
integer :: ieof(5),ifon(5)
data ieof/0,0,0,0,0/
data ifon/0,0,0,0,0/
integer :: nmove, nvola
integer :: iyear, itimob !obsnypatch
integer :: errcnt
DATA NMOVE/0/,NVOLA/61/
character*140 file_name_on_obs_unit, obs_domain_file_name
integer :: obs_domain_file_unit
! if(ieof(inest).eq.2.and.fdob%nstat.eq.0)then
! IF (iprt) print *,'returning from in4dob'
! return
! endif
! IF (iprt) print *,'start in4dob ',inest,xtime
IF(nudge_opt.NE.1)RETURN
! Initialize obs for printout
obs_prt = -999
newpass = .true.
errcnt = 0
! if start time, or restart time, set obs array to missing value
IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
DO N=1,NIOBF
TIMEOB(N)=99999.
ENDDO
fdob%xtime_at_rest = xtime !yliu 20080127
ENDIF
! set number of obs=0 if at start or restart
IF(KTAU.EQ.KTAUR)fdob%NSTAT=0
NSTA=fdob%NSTAT
XHOUR=XTIME/60.
XHOUR=AMAX1(XHOUR,0.0)
! DEFINE THE MAX LIMITS OF THE WINDOW
TBACK=XHOUR-TWINDO
TFORWD=XHOUR+TWINDO
IF (iprt) then
write(msg,fmt='(2(a,f8.3),a,i2)') &
'OBS NUDGING: Reading new obs for time window TBACK = ', &
tback,' TFORWD = ',tforwd,' for grid = ',inest
call wrf_message(msg)
ENDIF
! For obs that have become invalid because they are too old for the current time
! window, mark with 99999 to set up for removal from the rolling valid-obs list.
IF(NSTA.NE.0) THEN
NDUM=0
t_window : DO N=1,NSTA+1
IF((TIMEOB(N)-TBACK).LT.0) THEN
TIMEOB(N)=99999.
ENDIF
IF(TIMEOB(N).LT.9.E4) EXIT t_window
NDUM=N
ENDDO t_window
IF (iprt .and. ndum>0) THEN
write(msg,fmt='(a,i5,2a)') 'OBS NUDGING: ',ndum,' previously read obs ', &
'are now too old for the current window and have been removed.'
call wrf_message(msg)
ENDIF
! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
! IF (iprt) print *,'ndum at 20=',ndum
NDUM=ABS(NDUM)
NMOVE=NIOBF-NDUM
IF(NMOVE.GT.0 .AND. NDUM.NE.0 ) THEN
DO N=1,NMOVE
do KN = 1,nndgv
VAROBS(KN,N)=VAROBS(KN,N+NDUM)
enddo
! RIO is the west-east coordinate. RJO is south-north. (ajb)
RJO(N)=RJO(N+NDUM)
RIO(N)=RIO(N+NDUM)
RKO(N)=RKO(N+NDUM)
TIMEOB(N)=TIMEOB(N+NDUM)
nlevs_ob(n)=nlevs_ob(n+ndum)
lev_in_ob(n)=lev_in_ob(n+ndum)
plfo(n)=plfo(n+ndum)
elevob(n)=elevob(n+ndum)
ENDDO
ENDIF
NOPEN=NMOVE+1
IF(NOPEN.LE.NIOBF) THEN
DO N=NOPEN,NIOBF
do KN = 1,nndgv
VAROBS(KN,N)=99999.
enddo
RIO(N)=99999.
RJO(N)=99999.
RKO(N)=99999.
TIMEOB(N)=99999.
nlevs_ob(n)=99999.
lev_in_ob(n)=99999.
plfo(n)=99999.
elevob(n)=99999.
ENDDO
ENDIF
ENDIF
! Compute map projection info.
call set_projection(obs_proj, map_proj, cen_lat, cen_lon, &
true_lat1, true_lat2, stand_lon, &
known_lat, known_lon, &
e_we, e_sn, dxm, dym )
! FIND THE LAST OBS IN THE LIST
NLAST=0
last_ob : DO N=1,NIOBF
! print *,'nlast,n,timeob(n)=',nlast,n,timeob(n)
IF(TIMEOB(N).GT.9.E4) EXIT last_ob
NLAST=N
ENDDO last_ob
! print *,'in4dob, after 90 ',nlast,ktau,ktaur,nsta
! open file if at beginning or restart
IF(KTAU.EQ.0.OR.KTAU.EQ.KTAUR) THEN
fdob%RTLAST=-999.
obs_domain_file_unit = NVOLA+INEST-1
INQUIRE (obs_domain_file_unit,NAME=file_name_on_obs_unit,OPENED=OPENED)
!Build obs nudging file name (OBS_DOMAINX01 where X is the current domain number)
ifon(inest)=1
write(fonc(1:2),'(i2)')ifon(inest)
if(fonc(1:1).eq.' ')fonc(1:1)='0'
obs_domain_file_name='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2)
IF (.NOT. OPENED) THEN
INQUIRE (file=obs_domain_file_name,EXIST=exist)
if(exist)then
IF (iprt) THEN
write(msg,*) 'opening first fdda obs file, fonc=', &
fonc,' inest=',inest
call wrf_message(msg)
write(msg,*) 'ifon=',ifon(inest)
call wrf_message(msg)
ENDIF
OPEN(obs_domain_file_unit,FILE=obs_domain_file_name,FORM='FORMATTED',STATUS='OLD')
else
! no first file to open
IF (iprt) call wrf_message("there are no fdda obs files to open")
return
endif
ELSE
!If the unit for observation nudging file is already open make sure the file
!name matches the expected name to ensure that some other part of WRF has
!not opened a file using the same unit number
IF(file_name_on_obs_unit .ne. obs_domain_file_name) THEN
write(msg,*) 'File open on obs nudging unit (',obs_domain_file_unit,') with wrong name'
call wrf_message(msg)
write(msg,*) 'File open on obs nudging unit is named ',&
trim(adjustl(file_name_on_obs_unit)),' but it should be named ',&
trim(adjustl(obs_domain_file_name))
call wrf_message(msg)
write(msg,*) 'This likely means this unit number was opened elsewhere in WRF'
call wrf_message(msg)
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP Obs nudging file name mismatch' )
ENDIF
ENDIF
ENDIF !end if(KTAU.EQ.0.OR.KTAU.EQ.KTAUR)
! print *,'at jc check1'
!**********************************************************************
! -------------- BIG 100 LOOP OVER N --------------
!**********************************************************************
! NOW CHECK TO SEE IF EXTRA DATA MUST BE READ IN FROM THE
! DATA FILE. CONTINUE READING UNTIL THE REACHING THE EOF
! (DATA TIME IS NEGATIVE) OR FIRST TIME PAST TFORWD. THE
! LAST OBS CURRENTLY AVAILABLE IS IN N=NMOVE.
N=NLAST
IF(N.EQ.0)GOTO 110
1001 continue
! ieof=2 means no more files
IF(IEOF(inest).GT.1) then
GOTO 130
endif
100 CONTINUE
!ajb 20070116 bugfix for zero array index. N=0 if first obs is not in the domain.
IF(N.ne.0) THEN
IF(TIMEOB(N).GT.TFORWD.and.timeob(n).lt.99999.) THEN
GOTO 130
ENDIF
ENDIF
! OBSFILE: no more data in the obsfile
! AJB note: This is where we would implement multi-file reading.
if(ieof(inest).eq.1 )then
ieof(inest)=2
goto 130
endif
!**********************************************************************
! -------------- 110 SUBLOOP OVER N --------------
!**********************************************************************
110 continue
! THE TIME OF THE MOST RECENTLY ACQUIRED OBS IS .LE. TFORWD,
! SO CONTINUE READING
IF(N.GT.NIOBF-1)GOTO 120
! REPLACE NVOLA WITH LUN 70, AND USE NVOLA AS A FILE COUNTER
NVOL=NVOLA+INEST-1
IF(fdob%IEODI.EQ.1)GOTO 111
read(nvol,101,end=111,err=111)date_char
101 FORMAT(1x,a14)
n=n+1
! Convert the form of the observation date for geth_idts.
call fmt_date(date_char, obs_date)
! Compute the time period in seconds from the model reference
! date (fdob%sdate) until the observation date.
call geth_idts(obs_date, fdob%sdate(1:19), idts)
! Convert time in seconds to hours.
! In case of restart, correct for new sdate.
idts = idts + nint(fdob%xtime_at_rest*60.) ! yliu 20080127
rtimob =float(idts)/3600.
timeob(n)=rtimob
! print *,'read in ob ',n,timeob(n),rtimob
IF(IDYNIN.EQ.1.AND.TIMEOB(N)*60..GT.fdob%DATEND) THEN
IF (iprt) THEN
write(msg,*) ' IN4DOB: FOR INEST = ',INEST,' AT XTIME = ',XTIME, &
' TIMEOB = ',TIMEOB(N)*60.,' AND DATEND = ',fdob%DATEND,' :'
call wrf_message(msg)
write(msg,*) ' END-OF-DATA FLAG SET FOR OBS-NUDGING', &
' DYNAMIC INITIALIZATION'
call wrf_message(msg)
ENDIF
fdob%IEODI=1
TIMEOB(N)=99999.
rtimob=timeob(n)
ENDIF
read(nvol,102)latitude,longitude
102 FORMAT(2x,2(f9.4,1x))
! if(ifon.eq.4)print *,'ifon=4',latitude,longitude
! this works only for lc projection
! yliu: add llxy for all 3 projection
!ajb Arguments ri and rj have been switched from MM5 orientation.
CALL latlon_to_ij(obs_proj, latitude, longitude, ri, rj)
!ajb ri and rj are referenced to the non-staggered grid (not mass-pt!).
! (For MM5, they were referenced to the dot grid.)
ri = ri + .5 !ajb Adjust from mass-pt to non-staggered grid.
rj = rj + .5 !ajb Adjust from mass-pt to non-staggered grid.
rio(n)=ri
rjo(n)=rj
read(nvol,1021)id,namef
1021 FORMAT(2x,2(a40,3x))
read(nvol,103)platform,source,elevation,is_sound,bogus,meas_count
103 FORMAT( 2x,2(a16,2x),f8.0,2x,2(l4,2x),i5)
! write(6,*) '----- OBS description ----- '
! write(6,*) 'platform,source,elevation,is_sound,bogus,meas_count:'
! write(6,*) platform,source,elevation,is_sound,bogus,meas_count
! yliu
elevob(n)=elevation
! jc
! change platform from synop to profiler when needed
if(namef(2:9).eq.'PROFILER')platform(7:14)='PROFILER'
! yliu
if(namef(2:6).eq.'ACARS')platform(7:11)='ACARS'
if(namef(1:7).eq.'SATWNDS') platform(1:11)='SATWNDS '
if(namef(1:8).eq.'CLASS DA')platform(7:10)='TEMP'
! yliu end
rko(n)=-99.
!yliu 20050706
! if((platform(7:11).eq.'METAR').or.(platform(7:11).eq.'SPECI').or.
! 1 (platform(7:10).eq.'SHIP').or.(platform(7:11).eq.'SYNOP').or.
! 1 (platform(1:4).eq.'SAMS'))
! 1 rko(n)=1.0
if(.NOT. is_sound) rko(n)=1.0
!yliu 20050706 end
! plfo is inFORMATion on what platform. May use this later in adjusting weights
plfo(n)=99.
if(platform(7:11).eq.'METAR')plfo(n)=1.
if(platform(7:11).eq.'SPECI')plfo(n)=2.
if(platform(7:10).eq.'SHIP')plfo(n)=3.
if(platform(7:11).eq.'SYNOP')plfo(n)=4.
if(platform(7:10).eq.'TEMP')plfo(n)=5.
if(platform(7:11).eq.'PILOT')plfo(n)=6.
if(platform(1:7).eq.'SATWNDS')plfo(n)=7.
if(platform(7:11).eq.'SATWI')plfo(n)=7.
if(platform(1:4).eq.'SAMS')plfo(n)=8.
if(platform(7:14).eq.'PROFILER')plfo(n)=9.
! yliu: ACARS->SATWINDS
if(platform(7:11).eq.'ACARS')plfo(n)=7.
! yliu: end
if(plfo(n).eq.99.) then
IF (iprt) then
write(msg,*) 'n=',n,' unknown ob of type ',platform
call wrf_message(msg)
ENDIF
endif
!======================================================================
!======================================================================
! THIS PART READS SOUNDING INFO
IF(is_sound)THEN
nlevs_ob(n)=real(meas_count)
lev_in_ob(n)=1.
do imc=1,meas_count
! write(6,*) '0 inest = ',inest,' n = ',n
! the sounding has one header, many levels. This part puts it into
! "individual" observations. There's no other way for nudob to deal
! with it.
if(imc.gt.1)then ! sub-loop over N
n=n+1
if(n.gt.niobf)goto 120
nlevs_ob(n)=real(meas_count)
lev_in_ob(n)=real(imc)
timeob(n)=rtimob
rio(n)=ri
rjo(n)=rj
rko(n)=-99.
plfo(n)=plfo(n-imc+1)
elevob(n)=elevation
endif
read(nvol,104)pressure_data,pressure_qc, &
height_data,height_qc, &
temperature_data,temperature_qc, &
u_met_data,u_met_qc, &
v_met_data,v_met_qc, &
rh_data,rh_qc
104 FORMAT( 1x,6(f11.3,1x,f11.3,1x))
! yliu: Ensemble - add disturbance to upr obs
! if(plfo(n).eq.5.or.plfo(n).eq.6.or.plfo(n).eq.9) then FORE07E08
! if(imc.eq.1) then FORE07E08
! call srand(n)
! t_rand =- (rand(2)-0.5)*6
! call srand(n+100000)
! u_rand =- (rand(2)-0.5)*6
! call srand(n+200000)
! v_rand =- (rand(2)-0.5)*6
! endif FORE07E08
! if(temperature_qc.ge.0..and.temperature_qc.lt.30000..and.
! & temperature_data .gt. -88880.0 )
! & temperature_data = temperature_data + t_rand
! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
! make sure at least 1 of the components is .ne.0
! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
! u_met_data = u_met_data + u_rand
! v_met_data = v_met_data + v_rand
! endif
! endif FORE07E08
! yliu: Ens test - end
! jc
! hardwire to switch -777777. qc to 0. here temporarily
! -777777. is a sounding level that no qc was done on.
if(temperature_qc.eq.-777777.)temperature_qc=0.
if(pressure_qc.eq.-777777.)pressure_qc=0.
if(height_qc.eq.-777777.)height_qc=0.
if(u_met_qc.eq.-777777.)u_met_qc=0.
if(v_met_qc.eq.-777777.)v_met_qc=0.
if(rh_qc.eq.-777777.)rh_qc=0.
if(temperature_data.eq.-888888.)temperature_qc=-888888.
if(pressure_data.eq.-888888.)pressure_qc=-888888.
if(height_data.eq.-888888.)height_qc=-888888.
if(u_met_data.eq.-888888.)u_met_qc=-888888.
if(v_met_data.eq.-888888.)v_met_qc=-888888.
if(rh_data.eq.-888888.)rh_qc=-888888.
! jc
! Hardwire so that only use winds in pilot obs (no winds from temp) and
! only use temperatures and rh in temp obs (no temps from pilot obs)
! Do this because temp and pilot are treated as 2 platforms, but pilot
! has most of the winds, and temp has most of the temps. If use both,
! the second will smooth the effect of the first. Usually temps come in after
! pilots. pilots usually don't have any temps, but temp obs do have some
! winds usually.
! plfo=5 is TEMP ob, range sounding is an exception
!yliu start -- comment out to test with merged PILOT and TEMP and
! do not use obs interpolated by little_r
! if(plfo(n).eq.5. .and. namef(1:8).ne.'CLASS DA')then
! u_met_data=-888888.
! v_met_data=-888888.
! u_met_qc=-888888.
! v_met_qc=-888888.
! endif
if(plfo(n).eq.5..and.(u_met_qc.eq.256..or.v_met_qc.eq.256.))then
u_met_data=-888888.
v_met_data=-888888.
u_met_qc=-888888.
v_met_qc=-888888.
endif
!yliu end
! plfo=6 is PILOT ob
if(plfo(n).eq.6.)then
temperature_data=-888888.
rh_data=-888888.
temperature_qc=-888888.
rh_qc=-888888.
endif
!ajb Store temperature for WRF
! NOTE: The conversion to potential temperature, performed later in subroutine
! errob, requires good pressure data, either directly or via good height data.
! So here, in addition to checking for good temperature data, we must also
! do a check for good pressure or height.
if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
if( (pressure_qc.ge.0..and.pressure_qc.lt.30000.) .or. &
(height_qc .ge.0..and.height_qc .lt.30000.) ) then
varobs(3,n) = temperature_data
else
varobs(3,n)=-888888.
endif
else
varobs(3,n)=-888888.
endif
!ajb Store obs height
if(height_qc.ge.0..and.height_qc.lt.30000.)then
varobs(6,n)=height_data
else
varobs(6,n)=-888888.
endif
if(pressure_qc.ge.0..and.pressure_qc.lt.30000.)then
! if(pressure_qc.ge.0.)then
varobs(5,n)=pressure_data
else
varobs(5,n)=-888888.
IF (iprt) THEN
if(varobs(6,n).eq.-888888.000) then
if (errcnt.le.10) then
write(msg,*) '*** PROBLEM: sounding, p and ht undefined',latitude,longitude
call wrf_message(msg)
errcnt = errcnt + 1
if (errcnt.gt.10) call wrf_message("MAX of 10 warnings issued.")
endif
endif
ENDIF
endif
if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
! don't use data above 80 mb
if((varobs(5,n).gt.0.).and.(varobs(5,n).le.8.))then
u_met_data=-888888.
v_met_data=-888888.
u_met_qc=-888888.
v_met_qc=-888888.
temperature_data=-888888.
temperature_qc=-888888.
rh_data=-888888.
rh_qc=-888888.
endif
! Store horizontal wind components for WRF
if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
(v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
! make sure at least 1 of the components is .ne.0
(u_met_data.ne.0..or.v_met_data.ne.0.))then
! If Earth-relative wind vector, need to rotate it to grid-relative coords
if(u_met_qc.eq.129. .and. v_met_qc.eq.129.) then
CALL rotate_vector(longitude,u_met_data,v_met_data, &
obs_proj,map_proj)
endif
varobs(1,n)=u_met_data
varobs(2,n)=v_met_data
else
varobs(1,n)=-888888.
varobs(2,n)=-888888.
endif
r_data=-888888.
if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
if((pressure_qc.ge.0.).and.(temperature_qc.ge.0.).and. &
(pressure_qc.lt.30000.).and.(temperature_qc.lt.30000.))then
call rh2r(rh_data,temperature_data,pressure_data*.01, &
r_data,0) ! yliu, change last arg from 1 to 0
else
! print *,'rh, but no t or p to convert',temperature_qc, &
! pressure_qc,n
r_data=-888888.
endif
endif
varobs(4,n)=r_data
enddo ! end do imc=1,meas_count
! print *,'--- sdng n=',n,nlevs_ob(n),lev_in_ob(n),timeob(n)
! read in non-sounding obs
ELSEIF(.NOT.is_sound)THEN
nlevs_ob(n)=1.
lev_in_ob(n)=1.
read(nvol,105)slp_data,slp_qc, &
ref_pres_data,ref_pres_qc, &
height_data,height_qc, &
temperature_data,temperature_qc, &
u_met_data,u_met_qc, &
v_met_data,v_met_qc, &
rh_data,rh_qc, &
psfc_data,psfc_qc, &
precip_data,precip_qc
105 FORMAT( 1x,9(f11.3,1x,f11.3,1x))
! Ensemble: add disturbance to sfc obs
! call srand(n)
! t_rand =+ (rand(2)-0.5)*5
! call srand(n+100000)
! u_rand =+ (rand(2)-0.5)*5
! call srand(n+200000)
! v_rand =+ (rand(2)-0.5)*5
! if(temperature_qc.ge.0..and.temperature_qc.lt.30000. .and.
! & temperature_data .gt. -88880.0 )
! & temperature_data = temperature_data + t_rand
! if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and.
! & (v_met_qc.ge.0..and.v_met_qc.lt.30000.).and.
! make sure at least 1 of the components is .ne.0
! & (u_met_data.ne.0..or.v_met_data.ne.0.) .and.
! & (u_met_data.gt.-88880.0 .and. v_met_data.gt.-88880.0) )then
! u_met_data = u_met_data + u_rand
! v_met_data = v_met_data + v_rand
! endif
! yliu: Ens test - end
!Lilis
! calculate psfc if slp is there
if((psfc_qc.lt.0.).and.(slp_qc.ge.0..and.slp_qc.lt.30000.).and. &
(temperature_qc.ge.0..and.temperature_qc.lt.30000.).and. &
(slp_data.gt.90000.))then
tbar=temperature_data+0.5*elevation*.0065
psfc_data=slp_data*exp(-elevation/(rovg*tbar))
varobs(5,n)=psfc_data
psfc_qc=0.
endif
!c *No* **Very rough** estimate of psfc from sfc elevation if UUtah ob and elev>1000m
! estimate psfc from temp and elevation
! Do not know sfc pressure in model at this point.
! if((psfc_qc.lt.0.).and.(elevation.gt.1000.).and.
! 1 (temperature_qc.ge.0..and.temperature_qc.lt.30000.)
! 1 .and.(platform(7:16).eq.'SYNOP PRET'))then
if((psfc_qc.lt.0.).and. &
(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
tbar=temperature_data+0.5*elevation*.0065
psfc_data=100000.*exp(-elevation/(rovg*tbar))
varobs(5,n)=psfc_data
psfc_qc=0.
endif
if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and.(psfc_data.gt.70000. &
.and.psfc_data.lt.105000.))then
varobs(5,n)=psfc_data
else
varobs(5,n)=-888888.
endif
if(varobs(5,n).ge.0.)varobs(5,n)=varobs(5,n)*1.e-3
!Lilie
!ajb Store temperature for WRF
if(temperature_qc.ge.0..and.temperature_qc.lt.30000.)then
if((psfc_qc.ge.0..and.psfc_qc.lt.30000.).and. &
(psfc_data.gt.70000. .and.psfc_data.lt.105000.))then
varobs(3,n) = temperature_data
else
varobs(3,n)=-888888.
endif
else
varobs(3,n)=-888888.
endif
! Store horizontal wind components for WRF
if((u_met_qc.ge.0..and.u_met_qc.lt.30000.).and. &
(v_met_qc.ge.0..and.v_met_qc.lt.30000.).and. &
! make sure at least 1 of the components is .ne.0
(u_met_data.ne.0..or.v_met_data.ne.0.))then
! If Earth-relative wind vector, need to rotate it to grid-relative coords
if(u_met_qc.eq.129. .and. v_met_qc.eq.129.) then
CALL rotate_vector(longitude,u_met_data,v_met_data, &
obs_proj,map_proj)
endif
varobs(1,n)=u_met_data
varobs(2,n)=v_met_data
else
varobs(1,n)=-888888.
varobs(2,n)=-888888.
endif
! jc
! if a ship ob has rh<70%, then throw out
if(plfo(n).eq.3..and.rh_qc.ge.0..and.rh_data.lt.70.)then
rh_qc=-888888.
rh_data=-888888.
endif
!
r_data=-888888.
if(rh_qc.ge.0..and.rh_qc.lt.30000.)then
if((psfc_qc.ge.0..and.psfc_qc.lt.30000.) &
.and.(temperature_qc.ge.0..and.temperature_qc.lt.30000.))then
! rh_data=amin1(rh_data,96.) ! yliu: do not allow surface to be saturated
call rh2r(rh_data,temperature_data,psfc_data*.01, &
r_data,0) ! yliu, change last arg from 1 to 0
else
! print *,'rh, but no t or p',temperature_data,
! 1 psfc_data,n
r_data=-888888.
endif
endif
varobs(4,n)=r_data
ELSE
IF (iprt) THEN
call wrf_message(" ====== ")
call wrf_message(" NO Data Found ")
ENDIF
ENDIF !end if(is_sound)
! END OF SFC OBS INPUT SECTION
!======================================================================
!======================================================================
! check if ob time is too early (only applies to beginning)
IF(RTIMOB.LT.TBACK-TWINDO)then
IF (iprt) call wrf_message("ob too early")
n=n-1
GOTO 110
ENDIF
! check if this ob is a duplicate
! this check has to be before other checks
njend=n-1
if(is_sound)njend=n-meas_count
do njc=1,njend
! Check that time, lat, lon, and platform all match exactly.
! Platforms 1-4 (surface obs) can match with each other. Otherwise,
! platforms have to match exactly.
if( (timeob(n).eq.timeob(njc)) .and. &
(rio(n).eq.rio(njc)) .and. &
(rjo(n).eq.rjo(njc)) .and. &
(plfo(njc).ne.99.) ) then
!yliu: if two sfc obs are departed less than 1km, consider they are redundant
! (abs(rio(n)-rio(njc))*dscg.gt.1000.) &
! .or. (abs(rjo(n)-rjo(njc))*dscg.gt.1000.) &
! .or. (plfo(njc).eq.99.) )goto 801
!yliu end
! If platforms different, and either > 4, jump out
if( ( (plfo(n).le.4.).and.(plfo(njc).le.4.) ) .or. &
(plfo(n).eq.plfo(njc)) ) then
! if not a sounding, and levels are the same then replace first occurrence
if((.not.is_sound).and.(rko(njc).eq.rko(n))) then
! print *,'dup single ob-replace ',n,inest,
! plfo(n),plfo(njc)
! this is the sfc ob replacement part
do KN = 1,nndgv
VAROBS(KN,njc)=VAROBS(KN,n)
enddo
! don't need to switch these because they're the same
! RIO(njc)=RIO(n)
! RJO(njc)=RJO(n)
! RKO(njc)=RKO(n)
! TIMEOB(njc)=TIMEOB(n)
! nlevs_ob(njc)=nlevs_ob(n)
! lev_in_ob(njc)=lev_in_ob(n)
! plfo(njc)=plfo(n)
! end sfc ob replacement part
n=n-1
goto 100
! It's harder to fix the soundings, since the number of levels may be different
! The easiest thing to do is to just set the first occurrence to all missing, and
! keep the second occurrence, or vice versa.
! For temp or profiler keep the second, for pilot keep the one with more levs
! This is for a temp or prof sounding, equal to same
! also if a pilot, but second one has more obs
elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
( (plfo(njc).eq.5.).or.(plfo(njc).eq.9.).or. &
( (plfo(njc).eq.6.).and. &
(nlevs_ob(n).ge.nlevs_ob(njc)) ) ) )then
IF (iprt) THEN
write(msg,*) 'duplicate sounding - eliminate first occurrence', &
n,inest,meas_count,nlevs_ob(njc), &
latitude,longitude,plfo(njc)
call wrf_message(msg)
ENDIF
if(lev_in_ob(njc).ne.1.) then
IF (iprt) THEN
write(msg,*) 'problem ******* - dup sndg ', &
lev_in_ob(njc),nlevs_ob(njc)
call wrf_message(msg)
ENDIF
endif
! n=n-meas_count
! set the first sounding ob to missing
do njcc=njc,njc+nint(nlevs_ob(njc))-1
do KN = 1,nndgv
VAROBS(KN,njcc)=-888888.
enddo
plfo(njcc)=99.
enddo
goto 100
! if a pilot, but first one has more obs
elseif( (is_sound).and.(plfo(njc).eq.plfo(n)) .and. &
(plfo(njc).eq.6.).and. &
(nlevs_ob(n).lt.nlevs_ob(njc)) )then
IF (iprt) THEN
write(msg,*) &
'duplicate pilot sounding - eliminate second occurrence', &
n,inest,meas_count,nlevs_ob(njc), &
latitude,longitude,plfo(njc)
call wrf_message(msg)
ENDIF
if(lev_in_ob(njc).ne.1.) then
IF (iprt) THEN
write(msg,*) 'problem ******* - dup sndg ', &
lev_in_ob(njc),nlevs_ob(njc)
call wrf_message(msg)
ENDIF
endif
n=n-meas_count
!ajb Reset timeob for discarded indices.
do imc = n+1, n+meas_count
timeob(imc) = 99999.
enddo
goto 100
! This is for a single-level satellite upper air ob - replace first
elseif( (is_sound).and. &
(nlevs_ob(njc).eq.1.).and. &
(nlevs_ob(n).eq.1.).and. &
(varobs(5,njc).eq.varobs(5,n)).and. &
(plfo(njc).eq.7.).and.(plfo(n).eq.7.) ) then
IF (iprt) then
write(msg,*) &
'duplicate single lev sat-wind ob - replace first',n, &
inest,meas_count,varobs(5,n)
call wrf_message(msg)
ENDIF
! this is the single ua ob replacement part
do KN = 1,nndgv
VAROBS(KN,njc)=VAROBS(KN,n)
enddo
! don't need to switch these because they're the same
! RIO(njc)=RIO(n)
! RJO(njc)=RJO(n)
! RKO(njc)=RKO(n)
! TIMEOB(njc)=TIMEOB(n)
! nlevs_ob(njc)=nlevs_ob(n)
! lev_in_ob(njc)=lev_in_ob(n)
! plfo(njc)=plfo(n)
! end single ua ob replacement part
n=n-1
goto 100
else
! IF (iprt) THEN
! write(msg,*) 'duplicate location, but no match otherwise',n,njc, &
! plfo(n),varobs(5,n),nlevs_ob(n),lev_in_ob(n), &
! plfo(njc),varobs(5,njc),nlevs_ob(njc),lev_in_ob(njc)
! call wrf_message(msg)
! ENDIF
endif
endif
endif
! end of njc do loop
enddo
! check if ob is a sams ob that came in via UUtah - discard
if( plfo(n).eq.4..and.(platform(7:16).eq.'SYNOP PRET').and. &
(id(7:15).eq.'METNET= 3') )then
! print *,'elim metnet=3',latitude,longitude,rtimob
n=n-1
goto 100
endif
! check if ob is in the domain
if( (ri.lt.2.).or.(ri.gt.real(e_we-1)).or.(rj.lt.2.).or. &
(rj.gt.real(e_sn-1)) ) then
n=n-meas_count
!ajb Reset timeob for discarded indices.
do imc = n+1, n+meas_count
timeob(imc) = 99999.
enddo
goto 100
endif
IF(TIMEOB(N).LT.fdob%RTLAST) THEN
IF (iprt) THEN
call wrf_message("2 OBS ARE NOT IN CHRONOLOGICAL ORDER")
call wrf_message("NEW YEAR?")
write(msg,*) 'timeob,rtlast,n=',timeob(n),fdob%rtlast,n
call wrf_message(msg)
ENDIF
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 111' )
ELSE
fdob%RTLAST=TIMEOB(N)
ENDIF
! Save obs and model latitude and longitude for printout
CALL collect_obs_info(newpass,inest,n,latitude,longitude, &
nlast,nprev,niobf,id,stnid_prt, &
rio,rjo,prt_max,prt_freq,xlat,xlong, &
obs_prt,lat_prt,lon_prt,mlat_prt,mlon_prt, &
e_we,e_sn,ims,ime,jms,jme,its,ite,jts,jte)
GOTO 100
111 CONTINUE
!**********************************************************************
! -------------- END BIG 100 LOOP OVER N --------------
!**********************************************************************
if (iprt) then
write(msg,5403) NVOL,XTIME
call wrf_message(msg)
endif
IEOF(inest)=1
close(NVOLA+INEST-1)
IF (iprt) then
write(msg,*) 'closed fdda file for inest=',inest,nsta
call wrf_message(msg)
ENDIF
! AJB note: Go back and check for more files. (Multi-file implementation)
goto 1001
120 CONTINUE
! THE OBSERVATION ARRAYS ARE FULL AND THE MOST RECENTLY
! ACQUIRED OBS STILL HAS TIMEOB .LE. TFORWD. SO START
! DECREASING THE SIZE OF THE WINDOW
! get here if too many obs
IF (iprt) THEN
write(msg,121) N,NIOBF
call wrf_message(msg)
ENDIF
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 122' )
130 CONTINUE
! READ CYCLE IS COMPLETED. DETERMINE THE NUMBER OF OBS IN
! THE CURRENT WINDOW
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! BUT FIRST, WHEN KTAU.EQ.0 (OR IN GENERAL, KTAUR), DISCARD THE
! "OLD" OBS FIRST...
! Get here if at end of file, or if obs time is beyond what we need right now.
! On startup, we report the index of the last obs read.
! For restarts, we need to remove any old obs and then repack obs list.
IF(KTAU.EQ.KTAUR)THEN
NSTA=0
keep_obs : DO N=1,NIOBF
! try to keep all obs, but just don't use yet
! (don't want to throw away last obs read in - especially if
! its a sounding, in which case it looks like many obs)
IF(TIMEOB(N).GT.9.e4) EXIT keep_obs
if(timeob(n).gt.tforwd) then
if(iprt) then
write(msg,950) inest
call wrf_message(msg)
write(msg,951) n,timeob(n),tforwd
call wrf_message(msg)
endif
950 FORMAT('Saving index of first ob after end of current time window ', &
'for nest = ', i3,':')
951 FORMAT(' ob index = ',i8,', time of ob = ',f8.4, &
' hrs, end of time window = ',f8.4,' hrs')
endif
NSTA=N
ENDDO keep_obs
NDUM=0
! make time=99999. if ob is too old
! print *,'tback,nsta=',tback,nsta
old_obs : DO N=1,NSTA+1
IF((TIMEOB(N)-TBACK).LT.0)THEN
TIMEOB(N)=99999.
ENDIF
! print *,'n,ndum,timeob=',n,ndum,timeob(n)
IF(TIMEOB(N).LT.9.E4) EXIT old_obs
NDUM=N
ENDDO old_obs
! REMOVE OLD OBS DENOTED BY 99999. AT THE FRONT OF TIMEOB ARRAY
IF (iprt .and. ktaur > 0) THEN
write(msg,fmt='(a,i5,a)') 'OBS NUDGING: Upon restart, skipped over ',ndum, &
' obs that are now too old for the current obs window.'
call wrf_message(msg)
ENDIF
NDUM=ABS(NDUM)
NMOVE=NIOBF-NDUM
IF( NMOVE.GT.0 .AND. NDUM.NE.0) THEN
DO N=1,NMOVE
do KN = 1,nndgv
VAROBS(KN,N)=VAROBS(KN,N+NDUM)
enddo
RJO(N)=RJO(N+NDUM)
RIO(N)=RIO(N+NDUM)
RKO(N)=RKO(N+NDUM)
TIMEOB(N)=TIMEOB(N+NDUM)
nlevs_ob(n)=nlevs_ob(n+ndum)
lev_in_ob(n)=lev_in_ob(n+ndum)
plfo(n)=plfo(n+ndum)
ENDDO
ENDIF
! moved obs up. now fill remaining space with 99999.
NOPEN=NMOVE+1
IF(NOPEN.LE.NIOBF) THEN
DO N=NOPEN,NIOBF
do KN = 1,nndgv
VAROBS(KN,N)=99999.
enddo
RIO(N)=99999.
RJO(N)=99999.
RKO(N)=99999.
TIMEOB(N)=99999.
ENDDO
ENDIF
ENDIF
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
NSTA=0
! print *,'nsta at restart setting is ',nsta
! recalculate nsta after moving things around
recalc : DO N=1,NIOBF
! try to save all obs - don't throw away latest read in
IF(TIMEOB(N).GT.9.e4) EXIT recalc
NSTA=N
! nsta=n-1 ! yliu test
ENDDO recalc
! Find the number of stations that are actually within the time window.
nstaw = nvals_le_limit(nsta, timeob, tforwd)
IF (iprt) then
write(msg,160) KTAU,XTIME,NSTAW
call wrf_message(msg)
ENDIF
IF(KTAU.EQ.KTAUR)THEN
IF(nudge_opt.EQ.1)THEN
TWDOP=TWINDO*60.
IF (iprt) THEN
write(msg,1449) INEST,RINXY,RINSIG,TWDOP
call wrf_message(msg)
IF(ISWIND.EQ.1) then
write(msg,1450) GIV
call wrf_message(msg)
ELSE
write(msg,1455) INEST
call wrf_message("")
call wrf_message(msg)
call wrf_message("")
ENDIF
IF(ISTEMP.EQ.1) then
write(msg,1451) GIT
call wrf_message(msg)
ELSE
write(msg,1456) INEST
call wrf_message("")
call wrf_message(msg)
ENDIF
IF(ISMOIS.EQ.1) then
call wrf_message("")
write(msg,1452) GIQ
call wrf_message(msg)
ELSE
write(msg,1457) INEST
call wrf_message("")
call wrf_message(msg)
call wrf_message("")
ENDIF
ENDIF
ENDIF
ENDIF
IF(KTAU.EQ.KTAUR)THEN
IF(fdob%IWTSIG.NE.1)THEN
IF (iprt) THEN
write(msg,555)
call wrf_message(msg)
write(msg,556) fdob%RINFMN*RINXY,fdob%RINFMX*RINXY,fdob%PFREE*10.
call wrf_message(msg)
ENDIF
IF(fdob%RINFMN.GT.fdob%RINFMX) then
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 556' )
ENDIF
! IS MINIMUM GREATER THAN MAXIMUM?
IF (iprt) then
write(msg,557) fdob%DPSMX*10.,fdob%DCON
call wrf_message(msg)
ENDIF
IF(fdob%DPSMX.GT.10.) then
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP 557' )
ENDIF
ENDIF
ENDIF
IF(KTAU.EQ.KTAUR)THEN
IF (iprt) then
write(msg,601) INEST,IONF
call wrf_message(msg)
call wrf_message("")
ENDIF
ENDIF
fdob%NSTAT=NSTA
fdob%NSTAW=NSTAW
555 FORMAT(1X,' ABOVE THE SURFACE LAYER, OBS NUDGING IS PERFORMED', &
' ON PRESSURE LEVELS,')
556 FORMAT(1X,' WHERE RINXY VARIES LINEARLY FROM ',E11.3,' KM AT', &
' THE SURFACE TO ',E11.3,' KM AT ',F7.2,' MB AND ABOVE')
557 FORMAT(1X,' IN THE SURFACE LAYER, WXY IS A FUNCTION OF ', &
'DPSMX = ',F7.2,' MB WITH DCON = ',E11.3, &
' - SEE SUBROUTINE NUDOB')
601 FORMAT('FOR EFFICIENCY, THE OBS NUDGING FREQUENCY ', &
'FOR MESH #',I2,' IS ',1I2,' CGM TIMESTEPS ')
121 FORMAT(' WARNING: NOBS = ',I4,' IS GREATER THAN NIOBF = ', &
I4,': INCREASE PARAMETER NIOBF')
5403 FORMAT(1H0,'-------------EOF REACHED FOR NVOL = ',I3, &
' AND XTIME = ',F10.2,'-------------------')
160 FORMAT('****** CALL IN4DOB AT KTAU = ',I5,' AND XTIME = ', &
F10.2,': NSTA = ',I7,' ******')
1449 FORMAT('*****NUDGING INDIVIDUAL OBS ON MESH #',I2, &
' WITH RINXY = ', &
E11.3,' KM, RINSIG = ',E11.3,' AND TWINDO (HALF-PERIOD) = ', &
E11.3,' MIN')
1450 FORMAT(1X,'NUDGING IND. OBS WINDS WITH GIV = ',E11.3)
1451 FORMAT(1X,'NUDGING IND. OBS TEMPERATURE WITH GIT = ',E11.3)
1452 FORMAT(1X,'NUDGING IND. OBS MOISTURE WITH GIQ = ',E11.3)
1455 FORMAT(1X,'*** OBS WIND NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
1456 FORMAT(1X,'*** OBS TEMPERATURE NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
1457 FORMAT(1X,'*** OBS MOISTURE NUDGING FOR MESH ',I2,' IS TURNED OFF!!')
RETURN
END SUBROUTINE in4dob
SUBROUTINE julgmt(mdate,julgmtn,timanl,julday,gmt,ind)
! CONVERT MDATE YYMMDDHH TO JULGMT (JULIAN DAY * 100. +GMT)
! AND TO TIMANL (TIME IN MINUTES WITH RESPECT TO MODEL TIME)
! IF IND=0 INPUT MDATE, OUTPUT JULGMTN AND TIMANL
! IF IND=1 INPUT TIMANL, OUTPUT JULGMTN
! IF IND=2 INPUT JULGMTN, OUTPUT TIMANL
INTEGER, intent(in) :: MDATE
REAL, intent(out) :: JULGMTN
REAL, intent(out) :: TIMANL
INTEGER, intent(in) :: JULDAY
REAL, intent(in) :: GMT
INTEGER, intent(in) :: IND
!*** DECLARATIONS FOR IMPLICIT NONE
real :: MO(12), rjulanl, houranl, rhr
integer :: iyr, idate1, imo, idy, ihr, my1, my2, my3, ileap
integer :: juldayn, juldanl, idymax, mm
IF(IND.EQ.2)GOTO 150
IYR=INT(MDATE/1000000.+0.001)
IDATE1=MDATE-IYR*1000000
IMO=INT(IDATE1/10000.+0.001)
IDY=INT((IDATE1-IMO*10000.)/100.+0.001)
IHR=IDATE1-IMO*10000-IDY*100
MO(1)=31
MO(2)=28
! IS THE YEAR A LEAP YEAR? (IN THIS CENTURY)
IYR=IYR+1900
MY1=MOD(IYR,4)
MY2=MOD(IYR,100)
MY3=MOD(IYR,400)
ILEAP=0
! jc
! IF(MY1.EQ.0.AND.MY2.NE.0.OR.MY3.EQ.0)THEN
IF(MY1.EQ.0)THEN
ILEAP=1
MO(2)=29
ENDIF
IF(IND.EQ.1)GOTO 200
MO(3)=31
MO(4)=30
MO(5)=31
MO(6)=30
MO(7)=31
MO(8)=31
MO(9)=30
MO(10)=31
MO(11)=30
MO(12)=31
JULDAYN=0
DO 100 MM=1,IMO-1
JULDAYN=JULDAYN+MO(MM)
100 CONTINUE
IF(IHR.GE.24)THEN
IDY=IDY+1
IHR=IHR-24
ENDIF
JULGMTN=(JULDAYN+IDY)*100.+IHR
! CONVERT JULGMT TO TIMANL WRT MODEL TIME IN MINUTES (XTIME)
150 CONTINUE
JULDANL=INT(JULGMTN/100.+0.000001)
RJULANL=FLOAT(JULDANL)*100.
HOURANL=JULGMTN-RJULANL
TIMANL=(FLOAT(JULDANL-JULDAY)*24.-GMT+HOURANL)*60.
RETURN
200 CONTINUE
RHR=GMT+TIMANL/60.+0.000001
IDY=JULDAY
IDYMAX=365+ILEAP
300 IF(RHR.GE.24.0)THEN
RHR=RHR-24.0
IDY=IDY+1
GOTO 300
ENDIF
IF(IDY.GT.IDYMAX)IDY=IDY-IDYMAX
JULGMTN=FLOAT(IDY)*100.+RHR
RETURN
END SUBROUTINE julgmt
SUBROUTINE rh2r(rh,t,p,r,iice) 2
! convert rh to r
! if iice=1, use saturation with respect to ice
! rh is 0-100.
! r is g/g
! t is K
! p is mb
!
REAL, intent(in) :: rh
REAL, intent(in) :: t
REAL, intent(in) :: p
REAL, intent(out) :: r
INTEGER, intent(in) :: iice
!*** DECLARATIONS FOR IMPLICIT NONE
real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
real esat, rsat
eps=0.62197
e0=6.1078
eslcon1=17.2693882
eslcon2=35.86
esicon1=21.8745584
esicon2=7.66
t0=260.
! print *,'rh2r input=',rh,t,p
rh1=rh*.01
if(iice.eq.1.and.t.le.t0)then
esat=e0*exp(esicon1*(t-273.16)/(t-esicon2))
else
esat=e0*exp(eslcon1*(t-273.16)/(t-eslcon2))
endif
rsat=eps*esat/(p-esat)
! print *,'rsat,esat=',rsat,esat
r=rh1*rsat
! print *,'rh2r rh,t,p,r=',rh1,t,p,r
return
END SUBROUTINE rh2r
SUBROUTINE rh2rb(rh,t,p,r,iice),2
! convert rh to r
! if iice=1, use daturation with respect to ice
! rh is 0-100.
! r is g/g
! t is K
! p is mb
REAL, intent(in) :: rh
REAL, intent(in) :: t
REAL, intent(in) :: p
REAL, intent(out) :: r
INTEGER, intent(in) :: iice
!*** DECLARATIONS FOR IMPLICIT NONE
real eps, e0, eslcon1, eslcon2, esicon1, esicon2, t0, rh1
real esat, rsat
character(len=200) :: msg ! Argument to wrf_message
eps=0.622
e0=6.112
eslcon1=17.67
eslcon2=29.65
esicon1=22.514
esicon2=6.15e3
t0=273.15
write(msg,*) 'rh2r input=',rh,t,p
call wrf_message(msg)
rh1=rh*.01
if(iice.eq.1.and.t.le.t0)then
esat=e0*exp(esicon1-esicon2/t)
else
esat=e0*exp(eslcon1*(t-t0)/(t-eslcon2))
endif
rsat=eps*esat/(p-esat)
! print *,'rsat,esat=',rsat,esat
r=rh1*eps*rsat/(eps+rsat*(1.-rh1))
write(msg,*) 'rh2r rh,t,p,r=',rh1,t,p,r
call wrf_message(msg)
rh1=rh*.01
return
END SUBROUTINE rh2rb
SUBROUTINE set_projection (obs_proj, map_proj, cen_lat, cen_lon, & 1,7
true_lat1, true_lat2, stand_lon, &
known_lat, known_lon, &
e_we, e_sn, dxm, dym )
USE module_llxy
!*************************************************************************
! Purpose: Set map projection information which will be used to map the
! observation (lat,lon) location to its corresponding (x,y)
! location on the WRF (coarse) grid. using the selected map
! projection (e.g., Lambert, Mercator, Polar Stereo, etc).
!*************************************************************************
IMPLICIT NONE
TYPE(PROJ_INFO), intent(out) :: obs_proj ! structure for obs projection info.
INTEGER, intent(in) :: map_proj ! map projection index
REAL, intent(in) :: cen_lat ! center latitude for map projection
REAL, intent(in) :: cen_lon ! center longiture for map projection
REAL, intent(in) :: true_lat1 ! truelat1 for map projection
REAL, intent(in) :: true_lat2 ! truelat2 for map projection
REAL, intent(in) :: stand_lon ! standard longitude for map projection
INTEGER, intent(in) :: e_we ! max grid index in south-north coordinate
INTEGER, intent(in) :: e_sn ! max grid index in west-east coordinate
REAL, intent(in) :: known_lat ! latitude of domain origin point (i,j)=(1,1)
REAL, intent(in) :: known_lon ! longigude of domain origin point (i,j)=(1,1)
REAL, intent(in) :: dxm ! grid size in x (meters)
REAL, intent(in) :: dym ! grid size in y (meters)
! Set up map transformation structure
CALL map_init(obs_proj)
! Mercator
IF (map_proj == PROJ_MERC) THEN
CALL map_set(PROJ_MERC, obs_proj, &
truelat1 = true_lat1, &
lat1 = known_lat, &
lon1 = known_lon, &
knowni = 1., &
knownj = 1., &
dx = dxm)
! Lambert conformal
ELSE IF (map_proj == PROJ_LC) THEN
CALL map_set(PROJ_LC, obs_proj, &
truelat1 = true_lat1, &
truelat2 = true_lat2, &
stdlon = stand_lon, &
lat1 = known_lat, &
lon1 = known_lon, &
knowni = 1., &
knownj = 1., &
dx = dxm)
! Polar stereographic
ELSE IF (map_proj == PROJ_PS) THEN
CALL map_set(PROJ_PS, obs_proj, &
truelat1 = true_lat1, &
stdlon = stand_lon, &
lat1 = known_lat, &
lon1 = known_lon, &
knowni = 1., &
knownj = 1., &
dx = dxm)
! Cassini (global ARW)
ELSE IF (map_proj == PROJ_CASSINI) THEN
CALL map_set(PROJ_CASSINI, obs_proj, &
latinc = dym*360.0/(2.0*EARTH_RADIUS_M*PI), &
loninc = dxm*360.0/(2.0*EARTH_RADIUS_M*PI), &
lat1 = known_lat, &
lon1 = known_lon, &
! We still need to get POLE_LAT and POLE_LON metadata variables before
! this will work for rotated poles.
lat0 = 90.0, &
lon0 = 0.0, &
knowni = 1., &
knownj = 1., &
stdlon = stand_lon)
! Rotated latitude-longitude
ELSE IF (map_proj == PROJ_ROTLL) THEN
CALL map_set(PROJ_ROTLL, obs_proj, &
! I have no idea how this should work for NMM nested domains
ixdim = e_we-1, &
jydim = e_sn-1, &
phi = real(e_sn-2)*dym/2.0, &
lambda = real(e_we-2)*dxm, &
lat1 = cen_lat, &
lon1 = cen_lon, &
latinc = dym, &
loninc = dxm, &
stagger = HH)
END IF
END SUBROUTINE set_projection
SUBROUTINE fmt_date(idate,odate) !obsnypatch 1
!*************************************************************************
! Purpose: Re-format a character date string from YYYYMMDDHHmmss form
! to YYYY-MM-DD_HH:mm:ss form.
! INPUT:
! IDATE - Date string as YYYYMMDDHHmmss
! OUTPUT:
! ODATE - Date string as YYYY-MM-DD_HH:mm:ss
!*************************************************************************
IMPLICIT NONE
CHARACTER*14, intent(in) :: idate ! input date string
CHARACTER*19, intent(out) :: odate ! output date string
odate(1:19) = "0000-00-00_00:00:00"
odate(1:4) = idate(1:4) ! Year
odate(6:7) = idate(5:6) ! Month
odate(9:10) = idate(7:8) ! Day
odate(12:13) = idate(9:10) ! Hours
odate(15:16) = idate(11:12) ! Minutes
odate(18:19) = idate(13:14) ! Seconds
RETURN
END SUBROUTINE fmt_date
INTEGER FUNCTION nvals_le_limit(isize, values, limit) 1
!------------------------------------------------------------------------------
! PURPOSE: Return the number of values in a (real) non-decreasing array that
! are less than or equal to the specified upper limit.
! NOTE: It is important that the array is non-decreasing!
!
!------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: isize ! Size of input array
REAL, INTENT(IN) :: values(isize) ! Input array of reals
REAL, INTENT(IN) :: limit ! Upper limit
! Local variables
integer :: n
! Search the array from largest to smallest values.
find_nvals: DO n = isize, 1, -1
if(values(n).le.limit) EXIT find_nvals
ENDDO find_nvals
nvals_le_limit = n
RETURN
END FUNCTION nvals_le_limit
SUBROUTINE collect_obs_info(newpass,inest,n,latitude,longitude, & 1,4
nlast,nprev,niobf,station_id,stnid, &
rio,rjo,prt_max,prt_freq,xlat,xlong, &
obs, lat,lon, mlat,mlon, &
e_we,e_sn,ims,ime,jms,jme,its,ite,jts,jte)
!*************************************************************************
! Purpose: Collect the obs index, obs latitude, obs longitude, obs station
! id, and model latitude and longitude values for print
! diagnostics. Note that THIS SUBROUTINE IS CALLED INTERATIVELY
! FROM IN4DOB, WITHIN THE OBS READ LOOP that reads new obser-
! vations needed for the new time window. Flag newpass is true
! the first time collect_obs_info is called from the read-loop
! for a new time window. So for each pass of IN4DOB, newpass is
! true the first time IN4DOB calls collec_obs_info.
! OBS (soundings) contain multiple obs levels. So on each sub-
! sequent call of collect_obs_info for a specific pass of IN4DOB,
! n will jump by the number of levels in the sounding.
!
! Here, nlast refers to the index of the last valid-time obs
! that was read in during the last pass of IN4DOB (after the old
! obs were removed). This way we can properly start storing
! obs information for the new obs that are being read on this
! pass of IN4DOB, beginning with the first newly read obs for
! this pass of IN4DOB.
!
! Note that nprev is needed to properly handle soundings. On
! each pass, n is stored into nprev, and on each subsequent
! pass of collect_obs_info, a loop is performed from nprev+1 to n.
!*************************************************************************
IMPLICIT NONE
LOGICAL, intent(inout) :: newpass ! New pass flag
INTEGER, intent(in) :: inest ! nest index
INTEGER, intent(in) :: n ! Observation index
REAL, intent(in) :: latitude ! Latitude of obs
REAL, intent(in) :: longitude ! Latitude of obs
INTEGER, intent(in) :: nlast ! Last obs of valid obs, prev window
INTEGER, intent(inout) :: nprev ! Previous obs in new window read seq
INTEGER, intent(in) :: niobf ! Maximum number of observations
CHARACTER*15, intent(in) :: station_id ! First 15 chars of station id for obs n
INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
INTEGER, intent(inout) :: stnid(40,prt_max) ! Station ids for diagnostic printout
REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
REAL, DIMENSION( ims:ime, jms:jme ), &
intent(in ) :: xlat, xlong ! Lat/lon on mass-pt grid
INTEGER, intent(inout) :: obs(prt_max) ! Obs index for printout
REAL, intent(inout) :: lat(prt_max) ! Obs latitude for printout
REAL, intent(inout) :: lon(prt_max) ! Obs longitude for printout
REAL, intent(inout) :: mlat(prt_max) ! Model latitude at (rio,rjo) for printout
REAL, intent(inout) :: mlon(prt_max) ! Model longitude at (rio,rjo) for printout
INTEGER, intent(in) :: e_we ! Max grid index in south-north
INTEGER, intent(in) :: e_sn ! Max grid index in west-east
INTEGER, intent(in) :: ims ! Grid mem start (west-east)
INTEGER, intent(in) :: ime ! Grid mem end (west-east)
INTEGER, intent(in) :: jms ! Grid mem start (south-north)
INTEGER, intent(in) :: jme ! Grid mem end (south-north)
INTEGER, intent(in) :: its ! Grid tile start (west-east)
INTEGER, intent(in) :: ite ! Grid tile end (west-east)
INTEGER, intent(in) :: jts ! Grid tile start (south-north)
INTEGER, intent(in) :: jte ! Grid tile end (south-north)
! Local variables
integer i ! Loop counter over station id character
integer nn ! Loop counter over obs index
integer ndx,ndxp ! Index into printout arrays (ndx and prev ndx)
real :: ri, rj ! Mass-pt coord of obs
integer :: ril, rjl ! Mass-pt integer coord immed sw of obs
integer :: iend, jend ! Upper i, j index for interpolation
real :: dxob, dyob ! Grid fractions for interp
logical :: llsave ! Save lat/lon values if true
character(len=200) :: msg ! Argument to wrf_message
if(newpass) then
newpass = .false.
nprev = nlast ! Reset in case old obs have been discarded from prev window
endif
! Start iteration only if we have not yet stored prt_max number of obs for printing.
! Note: The loop below could represent multiple levels in a sounding, so we
! go ahead and start the loop if the beginning index (ndx) is prt_max or
! less, and then exit the loop if ndx exceeds prt_max.
if(prt_freq.gt.0) then
ndx = (n-nlast-1)/prt_freq + 1
ndxp = (nprev-nlast-1)/prt_freq + 1
else
write(msg,*) 'STOP! OBS NAMELIST INPUT obs_prt_freq MUST BE GREATER THAN ZERO.'
call wrf_message(msg)
write(msg,*) 'THE NAMELIST VALUE IS',prt_freq,' FOR NEST ',inest
call wrf_message(msg)
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob STOP' )
endif
! write(6,'5(a,i5),a,a15') 'n = ',n,' nlast = ',nlast,' ndx = ',ndx, &
! ' nprev = ',nprev,' ndxp = ',ndxp, &
! ' station id = ',station_id
if(ndxp .lt. prt_max) then
MODCHK : do nn = nprev+1, n
llsave = .false.
! if( mod(nn-1,prt_freq) .eq. 0 ) then
if( mod(nn-nlast-1,prt_freq) .eq. 0 ) then
ndx = (nn-nlast-1)/prt_freq + 1
if(ndx.gt.prt_max) EXIT MODCHK ! Limit printout to prt_max entries
llsave = .true.
endif
if(llsave) then
! Collect obs index and latitude and longitude.
obs(ndx) = nn
lat(ndx) = latitude
lon(ndx) = longitude
! Collect first 15 chars of obs station id (in integer format).
do i = 1,15
stnid(i,ndx) = ichar(station_id(i:i))
enddo
! Compute and collect the model latitude and longitude at the obs point.
CALL get_model_latlon(nn,niobf,rio,rjo,xlat,xlong,e_we,e_sn, &
ims,ime,jms,jme,its,ite,jts,jte, &
mlat(ndx),mlon(ndx))
endif !end if(llsave)
enddo MODCHK
endif !end if(ndx .le. prt_max)
! Save index of previous obs in read loop.
nprev = n
END SUBROUTINE collect_obs_info
SUBROUTINE get_model_latlon(n,niobf,rio,rjo,xlat,xlong,e_we,e_sn, & 1
ims,ime,jms,jme,its,ite,jts,jte, &
mlat,mlon)
!*************************************************************************
! Purpose: Use bilinear interpolation to compute the model latitude and
! longitude at the observation point.
!*************************************************************************
IMPLICIT NONE
INTEGER, intent(in) :: n ! Observation index
INTEGER, intent(in) :: niobf ! Maximum number of observations
REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
REAL, DIMENSION( ims:ime, jms:jme ), &
intent(in ) :: xlat, xlong ! Lat/lon on mass-pt grid
INTEGER, intent(in) :: e_we ! Max grid index in south-north
INTEGER, intent(in) :: e_sn ! Max grid index in west-east
INTEGER, intent(in) :: ims ! Grid mem start (west-east)
INTEGER, intent(in) :: ime ! Grid mem end (west-east)
INTEGER, intent(in) :: jms ! Grid mem start (south-north)
INTEGER, intent(in) :: jme ! Grid mem end (south-north)
INTEGER, intent(in) :: its ! Grid tile start (west-east)
INTEGER, intent(in) :: ite ! Grid tile end (west-east)
INTEGER, intent(in) :: jts ! Grid tile start (south-north)
INTEGER, intent(in) :: jte ! Grid tile end (south-north)
REAL, intent(out) :: mlat ! Model latitude at obs point
REAL, intent(out) :: mlon ! Model longitude at obs point
! Local variables
integer ndx ! Index into save arrays
real :: ri, rj ! Mass-pt coord of obs
integer :: ril, rjl ! Mass-pt integer coord immed sw of obs
integer :: iend, jend ! Upper i, j index for interpolation
real :: dxob, dyob ! Grid fractions for interp
! Compute model latitude and longitude if point on tile.
ri = rio(n) - .5 ! mass-pt west-east obs grid coord
rj = rjo(n) - .5 ! mass-pt south-north obs grid coord
ril = int(ri)
rjl = int(rj)
dxob = ri - float(ril)
dyob = rj - float(rjl)
iend = min(ite+1,e_we-2)
jend = min(jte+1,e_sn-2)
mlat = -999
mlon = -999
if(ri.ge.its .and. ri.lt.iend .and. rj.ge.jts .and. rj.lt.jend) then
! bilinear interpolation
mlat = ((1.-dyob)*((1.-dxob)*xlat(ril,rjl)+ &
dxob*xlat(ril+1,rjl) &
)+dyob*((1.-dxob)*xlat(ril,rjl+1)+ &
dxob*xlat(ril+1,rjl+1)))
mlon = ((1.-dyob)*((1.-dxob)*xlong(ril,rjl)+ &
dxob*xlong(ril+1,rjl) &
)+dyob*((1.-dxob)*xlong(ril,rjl+1)+ &
dxob*xlong(ril+1,rjl+1)))
endif
END SUBROUTINE get_model_latlon
SUBROUTINE rotate_vector(lon,u,v,obs_proj,map_proj) 2,1
USE module_llxy
!*************************************************************************
! Purpose: Rotate a single Earth-relative wind vector to a grid-relative
! wind vector.
!*************************************************************************
IMPLICIT NONE
REAL, intent(in) :: lon ! Longitude (deg)
REAL, intent(inout) :: u ! U-component of wind vector
REAL, intent(inout) :: v ! V-component of wind vector
TYPE(PROJ_INFO),intent(in) :: obs_proj ! Structure for obs projection
INTEGER, intent(in) :: map_proj ! Map projection index
! Local variables
real diff, alpha
double precision udbl, vdbl
! Only rotate winds for Lambert conformal or polar stereographic
if (map_proj == PROJ_LC .or. map_proj == PROJ_PS) then
diff = obs_proj%stdlon - lon
if (diff > 180.) then
diff = diff - 360.
else if (diff < -180.) then
diff = diff + 360.
end if
! Calculate the rotation angle, alpha, in radians
if (map_proj == PROJ_LC) then
alpha = diff * obs_proj%cone * rad_per_deg * obs_proj%hemi
else
alpha = diff * rad_per_deg * obs_proj%hemi
end if
udbl = v*sin(alpha) + u*cos(alpha)
vdbl = v*cos(alpha) - u*sin(alpha)
u = udbl
v = vdbl
endif
END SUBROUTINE rotate_vector
#endif
!-----------------------------------------------------------------------
! End subroutines for in4dob
!-----------------------------------------------------------------------