/wrfv2_fire/share/wrf_fddaobs_in.F
FORTRAN Legacy | 1917 lines | 1186 code | 202 blank | 529 comment | 93 complexity | e351a9244cba5754755f2c8cd21ad44d MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !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)
- 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) ) 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, &
- 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/
- ! 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.
- INQUIRE (NVOLA+INEST-1,OPENED=OPENED)
- IF (.NOT. OPENED) THEN
- ifon(inest)=1
- write(fonc(1:2),'(i2)')ifon(inest)
- if(fonc(1:1).eq.' ')fonc(1:1)='0'
- INQUIRE (file='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2) &
- ,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(NVOLA+INEST-1, &
- FILE='OBS_DOMAIN'//CHAR(INEST+ICHAR('0'))//fonc(1:2), &
- 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
- 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
- TI…
Large files files are truncated, but you can click here to view the full file