/wrfv2_fire/phys/module_fddaobs_rtfdda.F
FORTRAN Legacy | 2978 lines | 1716 code | 378 blank | 884 comment | 137 complexity | 951b09af3ef210354f9fd397a7ef6ab5 MD5 | raw file
Possible License(s): AGPL-1.0
- !WRF:MODEL_LAYER:PHYSICS
- !
- MODULE module_fddaobs_rtfdda
- ! 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
- !
- ! Modification History:
- ! 03212007 Modified fddaobs_init to compute Lambert cone factor. -Al Bourgeois
- CONTAINS
- !------------------------------------------------------------------------------
- SUBROUTINE fddaobs_init(nudge_opt, maxdom, inest, parid, &
- idynin, dtramp, fdaend, restart, &
- twindo_cg, twindo, itimestep, &
- no_pbl_nudge_uv, &
- no_pbl_nudge_t, &
- no_pbl_nudge_q, &
- sfc_scheme_horiz, sfc_scheme_vert, &
- maxsnd_gap, &
- sfcfact, sfcfacr, dpsmx, &
- nudge_wind, nudge_temp, nudge_mois, &
- nudgezfullr1_uv, nudgezrampr1_uv, &
- nudgezfullr2_uv, nudgezrampr2_uv, &
- nudgezfullr4_uv, nudgezrampr4_uv, &
- nudgezfullr1_t, nudgezrampr1_t, &
- nudgezfullr2_t, nudgezrampr2_t, &
- nudgezfullr4_t, nudgezrampr4_t, &
- nudgezfullr1_q, nudgezrampr1_q, &
- nudgezfullr2_q, nudgezrampr2_q, &
- nudgezfullr4_q, nudgezrampr4_q, &
- nudgezfullmin, nudgezrampmin, nudgezmax, &
- xlat, xlong, &
- start_year, start_month, start_day, &
- start_hour, start_minute, start_second, &
- p00, t00, tlp, &
- znu, p_top, &
- #if ( EM_CORE == 1 )
- fdob, &
- #endif
- iprt, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte)
- !-----------------------------------------------------------------------
- ! This routine does initialization for real time fdda obs-nudging.
- !
- !-----------------------------------------------------------------------
- USE module_model_constants, ONLY : g, r_d
- USE module_domain
- USE module_dm, ONLY : wrf_dm_min_real
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !=======================================================================
- ! Definitions
- !-----------
- INTEGER, intent(in) :: maxdom
- INTEGER, intent(in) :: nudge_opt(maxdom)
- INTEGER, intent(in) :: ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- INTEGER, intent(in) :: inest
- INTEGER, intent(in) :: parid(maxdom)
- INTEGER, intent(in) :: idynin ! flag for dynamic initialization
- REAL, intent(in) :: dtramp ! time period for idynin ramping (min)
- REAL, intent(in) :: fdaend(maxdom) ! nudging end time for domain (min)
- LOGICAL, intent(in) :: restart
- REAL, intent(in) :: twindo_cg ! time window on coarse grid
- REAL, intent(in) :: twindo
- INTEGER, intent(in) :: itimestep
- INTEGER , INTENT(IN) :: no_pbl_nudge_uv(maxdom) ! flags for no wind nudging in pbl
- INTEGER , INTENT(IN) :: no_pbl_nudge_t(maxdom) ! flags for no temperature nudging in pbl
- INTEGER , INTENT(IN) :: no_pbl_nudge_q(maxdom) ! flags for no moisture nudging in pbl
- INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
- INTEGER , INTENT(IN) :: sfc_scheme_vert ! vertical spreading scheme for surf obs (orig or regime vif)
- REAL , INTENT(IN) :: maxsnd_gap ! max allowed pressure gap in soundings, for interp (centibars)
- REAL, intent(in) :: sfcfact ! scale factor applied to time window for surface obs
- REAL, intent(in) :: sfcfacr ! scale fac applied to horiz rad of infl for sfc obs
- REAL, intent(in) :: dpsmx ! max press change allowed within hor rad of infl
- INTEGER , INTENT(IN) :: nudge_wind(maxdom) ! wind-nudging flag
- INTEGER , INTENT(IN) :: nudge_temp(maxdom) ! temperature-nudging flag
- INTEGER , INTENT(IN) :: nudge_mois(maxdom) ! moisture-nudging flag
- REAL, INTENT(IN) :: nudgezfullr1_uv ! vert infl fcn, regime=1 full-wt hght, winds
- REAL, INTENT(IN) :: nudgezrampr1_uv ! vert infl fcn, regime=1 ramp down hght, winds
- REAL, INTENT(IN) :: nudgezfullr2_uv ! vert infl fcn, regime=2 full-wt hght, winds
- REAL, INTENT(IN) :: nudgezrampr2_uv ! vert infl fcn, regime=2 ramp down hght, winds
- REAL, INTENT(IN) :: nudgezfullr4_uv ! vert infl fcn, regime=4 full-wt hght, winds
- REAL, INTENT(IN) :: nudgezrampr4_uv ! vert infl fcn, regime=4 ramp down hght, winds
- REAL, INTENT(IN) :: nudgezfullr1_t ! vert infl fcn, regime=1 full-wt hght, temp
- REAL, INTENT(IN) :: nudgezrampr1_t ! vert infl fcn, regime=1 ramp down hght, temp
- REAL, INTENT(IN) :: nudgezfullr2_t ! vert infl fcn, regime=2 full-wt hght, temp
- REAL, INTENT(IN) :: nudgezrampr2_t ! vert infl fcn, regime=2 ramp down hght, temp
- REAL, INTENT(IN) :: nudgezfullr4_t ! vert infl fcn, regime=4 full-wt hght, temp
- REAL, INTENT(IN) :: nudgezrampr4_t ! vert infl fcn, regime=4 ramp down hght, temp
- REAL, INTENT(IN) :: nudgezfullr1_q ! vert infl fcn, regime=1 full-wt hght, mois
- REAL, INTENT(IN) :: nudgezrampr1_q ! vert infl fcn, regime=1 ramp down hght, mois
- REAL, INTENT(IN) :: nudgezfullr2_q ! vert infl fcn, regime=2 full-wt hght, mois
- REAL, INTENT(IN) :: nudgezrampr2_q ! vert infl fcn, regime=2 ramp down hght, mois
- REAL, INTENT(IN) :: nudgezfullr4_q ! vert infl fcn, regime=4 full-wt hght, mois
- REAL, INTENT(IN) :: nudgezrampr4_q ! vert infl fcn, regime=4 ramp down hght, mois
- REAL, INTENT(IN) :: nudgezfullmin ! min dpth thru which vert infl fcn remains 1.0 (m)
- REAL, INTENT(IN) :: nudgezrampmin ! min dpth thru which vif decreases 1.0 to 0.0 (m)
- REAL, INTENT(IN) :: nudgezmax ! max dpth in which vif is nonzero (m)
- REAL, INTENT(IN) :: xlat ( ims:ime, jms:jme ) ! latitudes on mass-point grid
- REAL, INTENT(IN) :: xlong( ims:ime, jms:jme ) ! longitudes on mass-point grid
- INTEGER, intent(in) :: start_year ! Model start year
- INTEGER, intent(in) :: start_month ! Model start month
- INTEGER, intent(in) :: start_day ! Model start day
- INTEGER, intent(in) :: start_hour ! Model start hour
- INTEGER, intent(in) :: start_minute ! Model start minute
- INTEGER, intent(in) :: start_second ! Model start second
- REAL, INTENT(IN) :: p00 ! base state pressure
- REAL, INTENT(IN) :: t00 ! base state temperature
- REAL, INTENT(IN) :: tlp ! base state lapse rate
- REAL, INTENT(IN) :: p_top ! pressure at top of model
- REAL, INTENT(IN) :: znu( kms:kme ) ! eta values on half (mass) levels
- #if ( EM_CORE == 1 )
- TYPE(fdob_type), intent(inout) :: fdob
- #endif
- LOGICAL, intent(in) :: iprt ! Flag enabling printing warning messages
- ! Local variables
- logical :: nudge_flag ! nudging flag for this nest
- integer :: ktau ! current timestep
- integer :: nest ! loop counter
- integer :: idom ! domain id
- integer :: parent ! parent domain
- real :: conv ! 180/pi
- real :: tl1 ! Lambert standard parallel 1
- real :: tl2 ! Lambert standard parallel 2
- real :: xn1
- real :: known_lat ! Latitude of domain point (i,j)=(1,1)
- real :: known_lon ! Longitude of domain point (i,j)=(1,1)
- character(len=200) :: msg ! Argument to wrf_message
- real :: z_at_p( kms:kme ) ! height at p levels
- integer :: i,j,k ! loop counters
- #if ( EM_CORE == 1 )
- ! Check to see if the nudging flag has been set. If not,
- ! simply RETURN.
- nudge_flag = (nudge_opt(inest) .eq. 1)
- if (.not. nudge_flag) return
- call wrf_message("")
- write(msg,fmt='(a,i2)') ' OBSERVATION NUDGING IS ACTIVATED FOR MESH ',inest
- call wrf_message(msg)
- ktau = itimestep
- if(restart) then
- fdob%ktaur = ktau
- else
- fdob%ktaur = 0
- endif
- ! Create character string containing model starting-date
- CALL date_string(start_year, start_month, start_day, start_hour, &
- start_minute, start_second, fdob%sdate)
- ! Set flag for nudging on pressure (not sigma) surfaces.
- fdob%iwtsig = 0
- !**************************************************************************
- ! *** Initialize datend for dynamic initialization (ajb added 08132008) ***
- !**************************************************************************
- ! Set ending nudging date (used with dynamic ramp-down) to zero.
- fdob%datend = 0.
- fdob%ieodi = 0
- ! Check for dynamic initialization flag
- if(idynin.eq.1)then
- ! Set datend to time in minutes after which data are assumed to have ended.
- if(dtramp.gt.0.)then
- fdob%datend = fdaend(inest) - dtramp
- else
- fdob%datend = fdaend(inest)
- endif
- if(iprt) then
- call wrf_message("")
- write(msg,fmt='(a,i3,a)') &
- ' *** DYNAMIC-INITIALIZATION OPTION FOR INEST = ', inest, ' ***'
- call wrf_message(msg)
- write(msg,*) ' FDAEND,DATEND,DTRAMP: ',fdaend(inest),fdob%datend,dtramp
- call wrf_message(msg)
- call wrf_message("")
- endif
- endif
- ! *** end dynamic initialization section ***
- !**************************************************************************
- ! Store flags for surface obs spreading schemes
- if(sfc_scheme_horiz.eq.1) then
- call wrf_message('MM5 scheme selected for horizontal spreading of surface obs')
- elseif (sfc_scheme_horiz.eq.0) then
- call wrf_message('WRF scheme selected for horizontal spreading of surface obs')
- else
- write(msg,fmt='(a,i3)') 'Unknown h-spreading scheme for surface obs: ',sfc_scheme_horiz
- call wrf_message(msg)
- call wrf_message("Valid selections: 0=WRF scheme, 1=Original MM5 scheme")
- call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
- endif
- if(sfc_scheme_vert.eq.1) then
- call wrf_message('Original simple scheme selected for vertical spreading of surface obs')
- elseif (sfc_scheme_vert.eq.0) then
- call wrf_message("Regime-based VIF scheme selected for vertical spreading of surface obs")
- else
- write(msg,fmt='(a,i3)') 'Unknown v-spreading scheme for surface obs: ',sfc_scheme_vert
- call wrf_message(msg)
- call wrf_message("Valid selections: 0=Regime-based VIF scheme, 1=Original simple scheme")
- call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
- endif
- fdob%sfc_scheme_horiz = sfc_scheme_horiz
- fdob%sfc_scheme_vert = sfc_scheme_vert
- ! Store temporal and spatial scales
- fdob%sfcfact = sfcfact
- fdob%sfcfacr = sfcfacr
- ! Set time window.
- fdob%window = twindo
- call wrf_message("")
- write(msg,fmt='(a,i3)') '*** TIME WINDOW SETTINGS FOR NEST ',inest
- call wrf_message(msg)
- write(msg,fmt='(a,f6.3,2(a,f5.3))') ' TWINDO (hrs) = ',twindo, &
- ' SFCFACT = ',sfcfact,' SFCFACR = ',sfcfacr
- call wrf_message(msg)
- call wrf_message("")
- if(inest.eq.1) then
- if(twindo .eq. 0.) then
- if(iprt) then
- call wrf_message("")
- write(msg,*) '*** WARNING: TWINDO=0 on the coarse domain.'
- call wrf_message(msg)
- write(msg,*) '*** Did you forget to set twindo in the fdda namelist?'
- call wrf_message(msg)
- call wrf_message("")
- endif
- endif
- else ! nest
- if(twindo .eq. 0.) then
- fdob%window = twindo_cg
- if(iprt) then
- call wrf_message("")
- write(msg,fmt='(a,i2)') 'WARNING: TWINDO=0. for nest ',inest
- call wrf_message(msg)
- write(msg,fmt='(a,f12.5,a)') 'Default to coarse-grid value of ', twindo_cg,' hrs'
- call wrf_message(msg)
- call wrf_message("")
- endif
- endif
- endif
- ! Initialize flags.
- fdob%domain_tot=0
- do nest=1,maxdom
- fdob%domain_tot = fdob%domain_tot + nudge_opt(nest)
- end do
- ! Calculate and store dcon from dpsmx
- if(dpsmx.gt.0.) then
- fdob%dpsmx = dpsmx
- fdob%dcon = 1.0/fdob%dpsmx
- else
- call wrf_error_fatal('fddaobs_init: Namelist variable dpsmx must be greater than zero!')
- endif
- ! Calculate and store base-state heights at half (mass) levels.
- CALL get_base_state_height_column( p_top, p00, t00, tlp, g, r_d, znu, &
- fdob%base_state, kts, kte, kds,kde, kms,kme )
- ! Initialize flags for nudging within PBL.
- fdob%nudge_uv_pbl = .true.
- fdob%nudge_t_pbl = .true.
- fdob%nudge_q_pbl = .true.
- if(no_pbl_nudge_uv(inest) .eq. 1) fdob%nudge_uv_pbl = .false.
- if(no_pbl_nudge_t(inest) .eq. 1) fdob%nudge_t_pbl = .false.
- if(no_pbl_nudge_q(inest) .eq. 1) fdob%nudge_q_pbl = .false.
- if(no_pbl_nudge_uv(inest) .eq. 1) then
- fdob%nudge_uv_pbl = .false.
- write(msg,*) ' --> Obs nudging for U/V is turned off in PBL'
- call wrf_message(msg)
- endif
- if(no_pbl_nudge_t(inest) .eq. 1) then
- fdob%nudge_t_pbl = .false.
- write(msg,*) ' --> Obs nudging for T is turned off in PBL'
- call wrf_message(msg)
- endif
- if(no_pbl_nudge_q(inest) .eq. 1) then
- fdob%nudge_q_pbl = .false.
- write(msg,*) ' --> Obs nudging for Q is turned off in PBL'
- call wrf_message(msg)
- endif
- ! Store max allowed pressure gap for interpolating between soundings
- fdob%max_sndng_gap = maxsnd_gap
- write(msg,fmt='(a,f6.1)') &
- '*** MAX PRESSURE GAP (cb) for interpolation between soundings = ',maxsnd_gap
- call wrf_message(msg)
- call wrf_message("")
- ! Initialize vertical influence fcn for LML obs
- if(sfc_scheme_vert.eq.0) then
- fdob%vif_uv(1) = nudgezfullr1_uv
- fdob%vif_uv(2) = nudgezrampr1_uv
- fdob%vif_uv(3) = nudgezfullr2_uv
- fdob%vif_uv(4) = nudgezrampr2_uv
- fdob%vif_uv(5) = nudgezfullr4_uv
- fdob%vif_uv(6) = nudgezrampr4_uv
- fdob%vif_t (1) = nudgezfullr1_t
- fdob%vif_t (2) = nudgezrampr1_t
- fdob%vif_t (3) = nudgezfullr2_t
- fdob%vif_t (4) = nudgezrampr2_t
- fdob%vif_t (5) = nudgezfullr4_t
- fdob%vif_t (6) = nudgezrampr4_t
- fdob%vif_q (1) = nudgezfullr1_q
- fdob%vif_q (2) = nudgezrampr1_q
- fdob%vif_q (3) = nudgezfullr2_q
- fdob%vif_q (4) = nudgezrampr2_q
- fdob%vif_q (5) = nudgezfullr4_q
- fdob%vif_q (6) = nudgezrampr4_q
- ! Sanity checks
- if(nudgezmax.le.0.) then
- write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezmax MUST BE GREATER THAN ZERO.'
- call wrf_message(msg)
- write(msg,*) 'THE NAMELIST VALUE IS',nudgezmax
- call wrf_message(msg)
- call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgemax value' )
- endif
- if(nudgezfullmin.lt.0.) then
- write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezfullmin MUST BE NONNEGATIVE.'
- call wrf_message(msg)
- write(msg,*) 'THE NAMELIST VALUE IS',nudgezfullmin
- call wrf_message(msg)
- call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgefullmin value' )
- endif
- if(nudgezrampmin.lt.0.) then
- write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezrampmin MUST BE NONNEGATIVE.'
- call wrf_message(msg)
- write(msg,*) 'THE NAMELIST VALUE IS',nudgezrampmin
- call wrf_message(msg)
- call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgerampmin value' )
- endif
- if(nudgezmax.lt.nudgezfullmin+nudgezrampmin) then
- write(msg,*) 'STOP! INCONSISTENT OBS NAMELIST INPUTS.'
- call wrf_message(msg)
- write(msg,fmt='(3(a,f12.3))') 'obs_nudgezmax = ',nudgezmax, &
- ' obs_nudgezfullmin = ',nudgezfullmin, &
- ' obs_nudgezrampmin = ',nudgezrampmin
- call wrf_message(msg)
- write(msg,*) 'REQUIRE NUDGEZMAX >= NUDGEZFULLMIN + NUDGEZRAMPMIN'
- call wrf_message(msg)
- call wrf_error_fatal ( 'fddaobs_init: STOP on inconsistent namelist values' )
- endif
-
- fdob%vif_fullmin = nudgezfullmin
- fdob%vif_rampmin = nudgezrampmin
- fdob%vif_max = nudgezmax
-
- ! Check to make sure that if nudgzfullmin > 0, then it must be at least as large as the
- ! first model half-level will be anywhere in the domain at any time within the simulation.
- ! We use 1.1 times the base-state value fdob%base_state(1) for this purpose.
- if(nudgezfullmin.gt.0.0) then
- if(nudgezfullmin .lt. 1.1*fdob%base_state(1)) then
- fdob%vif_fullmin = 1.1*fdob%base_state(1)
- endif
- endif
- ! Print vertical weight info only if wind, temperature, or moisture nudging is requested.
- if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1) &
- .or. (nudge_mois(inest).eq.1) ) then
- call wrf_message("")
- write(msg,fmt='(a,i2,a)') ' *** SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest,' :'
- call wrf_message(msg)
- call wrf_message("")
- write(msg,fmt='(a,i5,a)') ' NUDGEZMAX: The maximum height at which nudging will be'// &
- ' applied from surface obs is ', nint(nudgezmax),' m AGL.'
- call wrf_message(msg)
- call wrf_message("")
- write(msg,fmt='(a,i3,a)') ' NUDGEZFULLMIN: The minimum height of full nudging weight'// &
- ' for surface obs is ', nint(fdob%vif_fullmin),' m.'
- call wrf_message(msg)
- if(nudgezfullmin.lt.fdob%vif_fullmin) then
- write(msg,fmt='(a,i3,a)') ' ***WARNING***: NUDGEZFULLMIN has been increased from'// &
- ' the user-input value of ',nint(nudgezfullmin),' m.'
- call wrf_message(msg)
- write(msg,fmt='(a,i3,a)') ' to ensure that at least the bottom model level is'// &
- ' included in full nudging.'
- call wrf_message(msg)
- endif
- call wrf_message("")
- write(msg,fmt='(a,i3,a)') ' NUDGEZRAMPMIN: The minimum height to ramp from full to no'// &
- ' nudging for surface obs is ', nint(nudgezrampmin),' m.'
- call wrf_message(msg)
- call wrf_message("")
- endif !endif either wind, temperature, or moisture nudging is requested
- ! Print vif settings
- if(nudge_wind(inest) .eq. 1) then
- call print_vif_var('wind', fdob%vif_uv, nudgezfullmin, nudgezrampmin)
- call wrf_message("")
- endif
- if(nudge_temp(inest) .eq. 1) then
- call print_vif_var('temp', fdob%vif_t, nudgezfullmin, nudgezrampmin)
- call wrf_message("")
- endif
- if(nudge_mois(inest) .eq. 1) then
- call print_vif_var('mois', fdob%vif_q, nudgezfullmin, nudgezrampmin)
- call wrf_message("")
- endif
- if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1) &
- .or. (nudge_mois(inest).eq.1) ) then
- write(msg,fmt='(a,i2)') ' *** END SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest
- call wrf_message(msg)
- call wrf_message("")
- endif
- endif !endif(sfc_scheme_vert.EQ.0)
- ! Set parameters.
- fdob%pfree = 50.0
- fdob%rinfmn = 1.0
- fdob%rinfmx = 2.0
- ! Get known lat and lon and store these on all processors in fdob structure, for
- ! later use in projection x-formation to map (lat,lon) to (x,y) for each obs.
- IF (its .eq. 1 .AND. jts .eq. 1) then
- known_lat = xlat(1,1)
- known_lon = xlong(1,1)
- ELSE
- known_lat = 9999.
- known_lon = 9999.
- END IF
- fdob%known_lat = wrf_dm_min_real(known_lat)
- fdob%known_lon = wrf_dm_min_real(known_lon)
- ! Calculate the nest levels, levidn. Note that each nest
- ! must know the nest levels levidn(maxdom) of each domain.
- do nest=1,maxdom
- ! Initialize nest level for each domain.
- if (nest .eq. 1) then
- fdob%levidn(nest) = 0 ! Mother domain has nest level 0
- else
- fdob%levidn(nest) = 1 ! All other domains start with 1
- endif
- idom = nest
- 100 parent = parid(idom) ! Go up the parent tree
- if (parent .gt. 1) then ! If not up to mother domain
- fdob%levidn(nest) = fdob%levidn(nest) + 1
- idom = parent
- goto 100
- endif
- enddo
- ! fdob%LML_OBS_HT1_LEV = kte
- ! HT1: do k = kte, kts, -1
- ! if( LML_HT1 .gt. z_at_p(k) ) then
- ! fdob%LML_OBS_HT1_LEV = k
- ! EXIT HT1
- ! endif
- ! enddo HT1
- ! fdob%LML_OBS_HT2_LEV = kte
- ! HT2: do k = kte, kts, -1
- ! if( LML_HT2 .gt. z_at_p(k) ) then
- ! fdob%LML_OBS_HT2_LEV = k
- ! EXIT HT2
- ! endif
- ! enddo HT2
- RETURN
- #endif
- END SUBROUTINE fddaobs_init
- #if ( EM_CORE == 1 )
- !-----------------------------------------------------------------------
- SUBROUTINE errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp, &
- z, &
- uratx, vratx, tratx, kpbl, &
- nndgv, nerrf, niobf, maxdom, &
- levidn, parid, nstat, nstaw, &
- iswind, istemp, ismois, ispstr, &
- timeob, rio, rjo, rko, &
- varobs, errf, ktau, xtime, &
- iratio, npfi, &
- prt_max, prt_freq, iprt, &
- obs_prt, stnid_prt, lat_prt, lon_prt, &
- mlat_prt, mlon_prt, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------------------
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- USE module_dm, ONLY : get_full_obs_vector, wrf_dm_sum_real
- #else
- USE module_dm, ONLY : wrf_dm_sum_real
- #endif
- USE module_model_constants, ONLY : rcp
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !
- ! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
- ! OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
- ! POINTS. THE OBSERVED VALUES CLOSEST TO THE CURRENT
- ! FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
- ! IN4DOB AND STORED IN ARRAY VAROBS. THE DIFFERENCES
- ! CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
- ! ERRF FOR THE NSTA OBSERVATION LOCATIONS. MISSING
- ! OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
- !
- ! HISTORY: Original author: MM5 version???
- ! 02/04/2004 - Creation of WRF version. Al Bourgeois
- ! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
- !------------------------------------------------------------------------------
- ! THE STORAGE ORDER IN ERRF IS AS FOLLOWS:
- ! IVAR VARIABLE TYPE(TAU-1)
- ! ---- --------------------
- ! 1 U error at obs loc
- ! 2 V error at obs loc
- ! 3 T error at obs loc
- ! 4 Q error at obs loc
- ! 5 Vertical layer index for PBL top at IOB, JOB
- ! 6 Model surface press at obs loc (T-points)
- ! 7 Model surface press at obs loc (U-points)
- ! 8 Model surface press at obs loc (V-points)
- ! 9 RKO at U-points
- !-----------------------------------------------------------------------
- !
- ! Description of input arguments.
- !
- !-----------------------------------------------------------------------
- INTEGER, INTENT(IN) :: inest ! Domain index.
- INTEGER, INTENT(IN) :: nndgv ! Number of nudge variables.
- INTEGER, INTENT(IN) :: nerrf ! Number of error fields.
- INTEGER, INTENT(IN) :: niobf ! Number of observations.
- INTEGER, INTENT(IN) :: maxdom ! Maximum number of domains.
- INTEGER, INTENT(IN) :: levidn(maxdom) ! Level of nest.
- INTEGER, INTENT(IN) :: parid(maxdom) ! Id of parent grid.
- INTEGER, INTENT(IN) :: ktau ! Model time step index
- REAL, INTENT(IN) :: xtime ! Forecast time in minutes
- INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
- INTEGER, INTENT(IN) :: npfi ! Coarse-grid diagnostics freq.
- INTEGER, INTENT(IN) :: prt_max ! Max number of obs to print.
- INTEGER, INTENT(IN) :: prt_freq ! Frequency of obs to print.
- LOGICAL, INTENT(IN) :: iprt ! Print flag
- INTEGER, INTENT(IN) :: obs_prt(prt_max) ! Print obs indices
- INTEGER, INTENT(IN) :: stnid_prt(40,prt_max) ! Print obs station ids
- REAL, INTENT(IN) :: lat_prt(prt_max) ! Print obs latitude
- REAL, INTENT(IN) :: lon_prt(prt_max) ! Print obs longitude
- REAL, INTENT(IN) :: mlat_prt(prt_max) ! Print model lat at obs loc
- REAL, INTENT(IN) :: mlon_prt(prt_max) ! Print model lon at obs loc
- INTEGER, INTENT(IN) :: nstat ! # stations held for use
- INTEGER, INTENT(IN) :: nstaw ! # stations in current window
- INTEGER, intent(in) :: iswind
- INTEGER, intent(in) :: istemp
- INTEGER, intent(in) :: ismois
- INTEGER, intent(in) :: ispstr
- INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
- INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
- INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
- REAL, INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
- REAL, INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
- REAL, INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
- REAL, INTENT(IN) :: t0
- REAL, INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
- REAL, INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
- REAL, INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
- REAL, INTENT(IN) :: rovcp
- REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! Ht above sl on half-levels
- REAL, INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
- REAL, INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
- REAL, INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
- INTEGER,INTENT(IN) :: kpbl( ims:ime, jms:jme ) ! Vertical layer index for PBL top
- REAL, INTENT(IN) :: timeob(niobf) ! Obs time (hrs)
- REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
- REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
- REAL, INTENT(INOUT) :: rko(niobf) ! Obs bottom-top coordinate
- REAL, INTENT(INOUT) :: varobs(nndgv, niobf)
- REAL, INTENT(INOUT) :: errf(nerrf, niobf)
- ! Local variables
- INTEGER :: iobmg(niobf) ! Obs i-coord on mass grid
- INTEGER :: jobmg(niobf) ! Obs j-coord on mass grid
- INTEGER :: ia(niobf)
- INTEGER :: ib(niobf)
- INTEGER :: ic(niobf)
- REAL :: pbbo(kds:kde) ! column base pressure (cb) at obs loc.
- REAL :: ppbo(kds:kde) ! column pressure perturbation (cb) at obs loc.
- REAL :: ra(niobf)
- REAL :: rb(niobf)
- REAL :: rc(niobf)
- REAL :: dxobmg(niobf) ! Interp. fraction (x dir) referenced to mass-grid
- REAL :: dyobmg(niobf) ! Interp. fraction (y dir) referenced to mass-grid
- INTEGER MM(MAXDOM)
- INTEGER NNL
- real :: uratio( ims:ime, jms:jme ) ! U to U10 ratio on momentum points.
- real :: vratio( ims:ime, jms:jme ) ! V to V10 ratio on momentum points.
- real :: pug1,pug2,pvg1,pvg2
- character(len=200) :: msg ! Argument to wrf_message
- ! Define staggers for U, V, and T grids, referenced from non-staggered grid.
- real, parameter :: gridx_t = 0.5 ! Mass-point x stagger
- real, parameter :: gridy_t = 0.5 ! Mass-point y stagger
- real, parameter :: gridx_u = 0.0 ! U-point x stagger
- real, parameter :: gridy_u = 0.5 ! U-point y stagger
- real, parameter :: gridx_v = 0.5 ! V-point x stagger
- real, parameter :: gridy_v = 0.0 ! V-point y stagger
- real :: dummy = 99999.
- real :: pbhi, pphi
- real :: obs_pottemp ! Potential temperature at observation
- !*** DECLARATIONS FOR IMPLICIT NONE
- integer nsta,ivar,n,ityp
- integer iob,job,kob,iob_ms,job_ms
- integer k,kbot,nml,nlb,nle
- integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
- integer i_start,i_end,j_start,j_end ! loop ranges for uratio,vratio calc.
- integer k_start,k_end
- integer ips ! For printing obs information
- integer pnx ! obs index for printout
- real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
- real pob
- real hob
- real uratiob,vratiob,tratiob,tratxob,fnpf
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- LOGICAL MP_LOCAL_DUMMASK(NIOBF) ! Mask for work to be done on this processor
- LOGICAL MP_LOCAL_UOBMASK(NIOBF) ! Dot-point mask
- LOGICAL MP_LOCAL_VOBMASK(NIOBF) ! Dot-point mask
- LOGICAL MP_LOCAL_COBMASK(NIOBF) ! Cross-point mask
- #endif
- ! LOGICAL, EXTERNAL :: TILE_MASK
- NSTA=NSTAT
- ! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
- ! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
- ! TO THE FINE MESH INDEX EQUIVALENTS
- ! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
- if (iprt) then
- write(msg,fmt='(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
- KTAU,' AND INEST = ',INEST,': NSTA = ',NSTAW,' ++++++'
- call wrf_message(msg)
- endif
- ERRF = 0.0 ! Zero out errf array
- ! Set up loop bounds for this grid's boundary conditions
- i_start = max( its-1,ids )
- i_end = min( ite+1,ide-1 )
- j_start = max( jts-1,jds )
- j_end = min( jte+1,jde-1 )
- k_start = kts
- k_end = min( kte, kde-1 )
- DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
- ! Set grid stagger
- IF(ITYP.EQ.1) THEN ! U-POINT CASE
- GRIDX = gridx_u
- GRIDY = gridy_u
- ELSE IF(ITYP.EQ.2) THEN ! V-POINT CASE
- GRIDX = gridx_v
- GRIDY = gridy_v
- ELSE ! MASS-POINT CASE
- GRIDX = gridx_t
- GRIDY = gridy_t
- ENDIF
- ! Compute URATIO and VRATIO fields on momentum (u,v) points.
- IF(ityp.eq.1)THEN
- call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
- ELSE IF (ityp.eq.2) THEN
- call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
- ENDIF
- IF(INEST.EQ.1) THEN ! COARSE MESH CASE...
- DO N=1,NSTA
- RA(N)=RIO(N)-GRIDX
- RB(N)=RJO(N)-GRIDY
- IA(N)=RA(N)
- IB(N)=RB(N)
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- DXOB=RA(N)-FLOAT(IA(N))
- DYOB=RB(N)-FLOAT(IB(N))
- ! Save mass-point arrays for computing rko for all var types
- if(ityp.eq.1) then
- iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
- jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
- dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
- dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
- endif
- iob_ms = iobmg(n)
- job_ms = jobmg(n)
- dxob_ms = dxobmg(n)
- dyob_ms = dyobmg(n)
- ! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION
- ! This is tricky: Initialize pob to zero in all procs. Only one proc actually
- ! calculates pob. If this is an obs to be converted from height-to-pressure, then
- ! by definition, varobs(5,n) will initially have the missing value -888888. After
- ! the pob calculation, pob will be zero in all procs except the one that calculated
- ! it, and so pob is dm_summed over all procs and stored into varobs(5,n). So on
- ! subsequent passes, the dm_sum will not occur because varobs(5,n) will not have a
- ! missing value. If this is not an obs-in-height,
- pob = 0.0
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Set mask for obs to be handled by this processor
- MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
-
- IF ( MP_LOCAL_DUMMASK(N) ) THEN
- #endif
- ! Interpolate pressure to obs location column and convert from Pa to cb.
- do k = kds, kde
- pbbo(k) = .001*( &
- (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
- DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
- DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
- DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
- ppbo(k) = .001*( &
- (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
- DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
- DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
- DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
- enddo
- !ajb 20040119: Note based on bugfix for dot/cross points split across processors,
- !ajb which was actually a serial code fix: The ityp=2 (v-points) and
- !ajb ityp=3 (mass-points) cases should not use the ityp=1 (u-points)
- !ajb case rko! This is necessary for bit-for-bit reproducability
- !ajb with the parallel run. (coarse mesh)
- if(abs(rko(n)+99).lt.1.)then
- pob = varobs(5,n)
- if(pob .eq.-888888.) then
- hob = varobs(6,n)
- if(hob .gt. -800000. ) then
- pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms, &
- dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
- ims,ime, jms,jme, kms,kme )
- endif
- endif
- if(pob .gt.-800000.)then
- do k=k_end-1,1,-1
- kbot = k
- if(pob .le. pbbo(k)+ppbo(k)) then
- goto 199
- endif
- enddo
- 199 continue
- pphi = ppbo(kbot+1)
- pbhi = pbbo(kbot+1)
- rko(n) = real(kbot+1)- &
- ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
- rko(n)=max(rko(n),1.0)
- endif
- endif
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
- #endif
- ! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
- if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
- varobs(5,n) = wrf_dm_sum_real ( pob )
- endif
- RC(N)=RKO(N)
- ENDDO ! END COARSE MESH LOOP OVER NSTA
- ELSE ! FINE MESH CASE
- !**********************************************************************
- !ajb_07012008: Conversion of obs coordinates to the fine mesh here
- !ajb is no longer necessary, since implementation of the WRF map pro-
- !ajb jections (to each nest directly) is implemented in sub in4dob.
- !**********************************************************************
- !ajb
- !ajb GET (I,J,K) OF OBSERVATIONS ON FINE MESH VALUES.
- DO N=1,NSTA
- RA(N)=RIO(N)-GRIDX ! ajb added 07012008
- RB(N)=RJO(N)-GRIDY ! ajb added 07012008
- IA(N)=RA(N)
- IB(N)=RB(N)
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- DXOB=RA(N)-FLOAT(IA(N))
- DYOB=RB(N)-FLOAT(IB(N))
- ! Save mass-point arrays for computing rko for all var types
- if(ityp.eq.1) then
- iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
- jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
- dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
- dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
- endif
- iob_ms = iobmg(n)
- job_ms = jobmg(n)
- dxob_ms = dxobmg(n)
- dyob_ms = dyobmg(n)
- ! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION (see COARSE MESH comments)
- pob = 0.0
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Set mask for obs to be handled by this processor
- MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
- IF ( MP_LOCAL_DUMMASK(N) ) THEN
- #endif
- ! Interpolate pressure to obs location column and convert from Pa to cb.
- do k = kds, kde
- pbbo(k) = .001*( &
- (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
- DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
- DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
- DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
- ppbo(k) = .001*( &
- (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
- DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
- DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
- DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
- enddo
- !ajb 20040119: Note based on bugfix for dot/cross points split across processors,
- !ajb which was actually a serial code fix: The ityp=2 (v-points) and
- !ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
- !ajb case) rko! This is necessary for bit-for-bit reproducability
- !ajb with parallel run. (fine mesh)
- if(abs(rko(n)+99).lt.1.)then
- pob = varobs(5,n)
- if(pob .eq.-888888.) then
- hob = varobs(6,n)
- if(hob .gt. -800000. ) then
- pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms, &
- dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
- ims,ime, jms,jme, kms,kme )
- endif
- endif
- if(pob .gt.-800000.)then
- do k=k_end-1,1,-1
- kbot = k
- if(pob .le. pbbo(k)+ppbo(k)) then
- goto 198
- endif
- enddo
- 198 continue
- pphi = ppbo(kbot+1)
- pbhi = pbbo(kbot+1)
- rko(n) = real(kbot+1)- &
- ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
- rko(n)=max(rko(n),1.0)
- endif
- endif
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
- #endif
- ! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
- if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
- varobs(5,n) = wrf_dm_sum_real ( pob )
- endif
- RC(N)=RKO(N)
- ENDDO ! END FINE MESH LOOP OVER NSTA
-
- ENDIF ! end if(inest.eq.1)
- ! INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
- ! AND THE LOCAL FORECAST VALUES TO ZERO. FOR THE FINE MESH
- ! ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
- ! OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
- ! SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
- ! CURRENT DOMAIN
- IF(ITYP.EQ.1) THEN
- NLB=1
- NLE=1
- ELSE IF(ITYP.EQ.2) THEN
- NLB=2
- NLE=2
- ELSE
- NLB=3
- NLE=5
- ENDIF
- DO IVAR=NLB,NLE
- DO N=1,NSTA
- IF((RA(N)-1.).LT.0)THEN
- ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
- ENDIF
- IF((RB(N)-1.).LT.0)THEN
- ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
- ENDIF
- IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
- ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
- ENDIF
- IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
- ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
- ENDIF
- if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
- ENDDO
- ENDDO
- ! NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
- ! GRID POINT TOWARD THE LOWER LEFT
- DO N=1,NSTA
- IA(N)=RA(N)
- IB(N)=RB(N)
- IC(N)=RC(N)
- ENDDO
- DO N=1,NSTA
- RA(N)=RA(N)-FLOAT(IA(N))
- RB(N)=RB(N)-FLOAT(IB(N))
- RC(N)=RC(N)-FLOAT(IC(N))
- ENDDO
- ! PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
- ! TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
- ! POINTS FOR U, V, T, AND Q.
- ! Compute local masks for dot and cross points.
- if(ityp.eq.1) then
- DO N=1,NSTA
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Set mask for U-momemtum points to be handled by this processor
- MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
- #endif
- ENDDO
- endif
- if(ityp.eq.2) then
- DO N=1,NSTA
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Set mask for V-momentum points to be handled by this processor
- MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
- #endif
- ENDDO
- endif
- if(ityp.eq.3) then
- DO N=1,NSTA
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Set mask for cross (mass) points to be handled by this processor
- MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
- #endif
- ENDDO
- endif
- !**********************************************************
- ! PROCESS U VARIABLE (IVAR=1)
- !**********************************************************
- IF(ITYP.EQ.1) THEN
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- DO N=1,NSTA
- IF(MP_LOCAL_UOBMASK(N)) THEN
- ERRF(9,N)=rko(n) !RKO is needed by neighboring processors !ajb
- ENDIF
- ENDDO
- #endif
- IF(ISWIND.EQ.1) THEN
- DO N=1,NSTA
- IOB=MAX0(2,IA(N))
- IOB=MIN0(IOB,ide-1)
- IOBM=MAX0(1,IOB-1)
- IOBP=MIN0(ide-1,IOB+1)
- JOB=MAX0(2,IB(N))
- JOB=MIN0(JOB,jde-1)
- JOBM=MAX0(1,JOB-1)
- JOBP=MIN0(jde-1,JOB+1)
- KOB=MIN0(K_END,IC(N))
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- IF(MP_LOCAL_UOBMASK(N))THEN ! Do if obs on this processor
- #endif
- KOBP=MIN0(KOB+1,k_end)
- DXOB=RA(N)
- DYOB=RB(N)
- DZOB=RC(N)
- ! Compute surface pressure values at surrounding U and V points
- PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
- PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
- ! This is to correct obs to model sigma level using reverse similarity theory
- if(rko(n).eq.1.0)then
- uratiob=((1.-DYOB)*((1.-DXOB)*uratio(IOB,JOB)+ &
- DXOB*uratio(IOBP,JOB) &
- )+DYOB*((1.-DXOB)*uratio(IOB,JOBP)+ &
- DXOB*uratio(IOBP,JOBP)))
- else
- uratiob=1.
- endif
- !YLIU Some PBL scheme do not define the vratio/uratio
- if(abs(uratiob).lt.1.0e-3) then
- uratiob=1.
- endif
- ! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
- ! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
- ! OUTSIDE THE CURRENT DOMAIN
- ! U COMPONENT WIND ERROR
- ERRF(1,N)=ERRF(1,N)+uratiob*VAROBS(1,N)-((1.-DZOB)* &
- ((1.-DyOB)*((1.- &
- DxOB)*UB(IOB,KOB,JOB)+DxOB*UB(IOB+1,KOB,JOB) &
- )+DyOB*((1.-DxOB)*UB(IOB,KOB,JOB+1)+DxOB* &
- UB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
- *UB(IOB,KOBP,JOB)+DxOB*UB(IOB+1,KOBP,JOB))+ &
- DyOB*((1.-DxOB)*UB(IOB,KOBP,JOB+1)+DxOB* &
- UB(IOB+1,KOBP,JOB+1))))
- ! if(n.le.10) then
- ! write(6,*)
- ! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob, &
- ! ' N = ',n,' inest = ',inest
- ! write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
- ! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
- ! write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
- ! write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
- ! write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
- ! write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
- ! write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
- ! write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
- ! write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
- ! write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
- ! write(6,*) 'uratiob = ',uratiob
- ! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
- ! write(6,*) 'ERRF(1,N) = ',errf(1,n)
- ! write(6,*)
- ! endif
- ! Store model surface pressure (not the error!) at U point.
- ERRF(7,N)=.001*( (1.-DXOB)*PUG1 + DXOB*PUG2 )
-
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF ! end IF( MP_LOCAL_UOBMASK(N) )
- #endif
- ENDDO ! END U-POINT LOOP OVER OBS
- ENDIF ! end if(iswind.eq.1)
- ENDIF ! ITYP=1: PROCESS U
- !**********************************************************
- ! PROCESS V VARIABLE (IVAR=2)
- !**********************************************************
- IF(ITYP.EQ.2) THEN
- IF(ISWIND.EQ.1) THEN
- DO N=1,NSTA
- IOB=MAX0(2,IA(N))
- IOB=MIN0(IOB,ide-1)
- IOBM=MAX0(1,IOB-1)
- IOBP=MIN0(ide-1,IOB+1)
- JOB=MAX0(2,IB(N))
- JOB=MIN0(JOB,jde-1)
- JOBM=MAX0(1,JOB-1)
- JOBP=MIN0(jde-1,JOB+1)
- KOB=MIN0(K_END,IC(N))
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- IF(MP_LOCAL_VOBMASK(N))THEN ! Do if obs on this processor
- #endif
- KOBP=MIN0(KOB+1,k_end)
- DXOB=RA(N)
- DYOB=RB(N)
- DZOB=RC(N)
- ! Compute surface pressure values at surrounding U and V points
- PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
- PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
- ! This is to correct obs to model sigma level using reverse similarity theory
- if(rko(n).eq.1.0)then
- vratiob=((1.-DYOB)*((1.-DXOB)*vratio(IOB,JOB)+ &
- DXOB*vratio(IOBP,JOB) &
- )+DYOB*((1.-DXOB)*vratio(IOB,JOBP)+ &
- DXOB*vratio(IOBP,JOBP)))
- else
- vratiob=1.
- endif
- !YLIU Some PBL scheme do not define the vratio/uratio
- if(abs(vratiob).lt.1.0e-3) then
- vratiob=1.
- endif
- ! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
- ! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
- ! OUTSIDE THE CURRENT DOMAIN
-
- ! V COMPONENT WIND ERROR
- ERRF(2,N)=ERRF(2,N)+vratiob*VAROBS(2,N)-((1.-DZOB)* &
- ((1.-DyOB)*((1.- &
- DxOB)*VB(IOB,KOB,JOB)+DxOB*VB(IOB+1,KOB,JOB) &
- )+DyOB*((1.-DxOB)*VB(IOB,KOB,JOB+1)+DxOB* &
- VB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
- *VB(IOB,KOBP,JOB)+DxOB*VB(IOB+1,KOBP,JOB))+ &
- DyOB*((1.-DxOB)*VB(IOB,KOBP,JOB+1)+DxOB* &
- VB(IOB+1,KOBP,JOB+1))))
- ! if(n.le.10) then
- ! write(6,*)
- ! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF2 at ',iob,job,kob, &
- ! ' N = ',n,' inest = ',inest
- ! write(6,*) 'VAROBS(2,N) = ',varobs(2,n)
- ! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
- ! write(6,*) 'VB(IOB,KOB,JOB) = ',VB(IOB,KOB,JOB)
- ! write(6,*) 'ERRF(2,N) = ',errf(2,n)
- ! write(6,*) 'vratiob = ',vratiob
- ! write(6,*)
- ! endif
-
- ! Store model surface pressure (not the error!) at V point.
- ERRF(8,N)=.001*( (1.-DYOB)*PVG1 + DYOB*PVG2 )
-
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF ! end IF( MP_LOCAL_VOBMASK(N) )
- #endif
- ENDDO ! END V-POINT LOOP OVER OBS
- ENDIF ! end if(iswind.eq.1)
- ENDIF ! ITYP=1: PROCESS V
- !**********************************************************
- ! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
- !**********************************************************
- IF(ITYP.EQ.3) THEN
- IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
- DO N=1,NSTA
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
- #endif
- KOB=MIN0(k_end,IC(N))
- KOBP=MIN0(KOB+1,K_END)
- DXOB=RA(N)
- DYOB=RB(N)
- DZOB=RC(N)
- ! This is to correct obs to model sigma level using reverse similarity theory
- if(rko(n).eq.1.0)then
- tratxob=((1.-DYOB)*((1.-DXOB)*tratx(IOB,JOB)+ &
- DXOB*tratx(IOB+1,JOB) &
- )+DYOB*((1.-DXOB)*tratx(IOB,JOB+1)+ &
- DXOB*tratx(IOB+1,JOB+1)))
- else
- tratxob=1.
- endif
- !yliu
- if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
- !ajb We must convert temperature to potential temperature
- obs_pottemp = -888888.
- if(varobs(3,n).gt.-800000. .and. varobs(5,n).gt.-800000) then
- obs_pottemp = varobs(3,n)*(100./varobs(5,n))**RCP - t0
- endif
- ERRF(3,N)=ERRF(3,N)+tratxob*obs_pottemp-((1.-DZOB)* &
- ((1.-DyOB)*((1.- &
- DxOB)*(TB(IOB,KOB,JOB))+DxOB* &
- (TB(IOB+1,KOB,JOB)))+DyOB*((1.-DxOB)* &
- (TB(IOB,KOB,JOB+1))+DxOB* &
- (TB(IOB+1,KOB,JOB+1))))+DZOB*((1.- &
- DyOB)*((1.-DxOB)*(TB(IOB,KOBP,JOB))+DxOB* &
- (TB(IOB+1,KOBP,JOB)))+DyOB*((1.-DxOB)* &
- (TB(IOB,KOBP,JOB+1))+DxOB* &
- (TB(IOB+1,KOBP,JOB+1)))))
- ! if(n.le.10) then
- ! write(6,*)
- ! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob, &
- ! ' N = ',n,' inest = ',inest
- ! write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
- ! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
- ! write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
- ! write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
- ! write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
- ! write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
- ! write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
- ! write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
- ! write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
- ! write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
- ! write(6,*) 'tratxob = ',tratxob
- ! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
- ! write(6,*) 'ERRF(3,N) = ',errf(3,n)
- ! write(6,*)
- ! endif
- ! MOISTURE ERROR
- ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)*((1.-DyOB)*((1.- &
- DxOB)*QVB(IOB,KOB,JOB)+DxOB* &
- QVB(IOB+1,KOB,JOB))+DyOB*((1.-DxOB)* &
- QVB(IOB,KOB,JOB+1)+DxOB* &
- QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.- &
- DyOB)*((1.-DxOB)*QVB(IOB,KOBP,JOB)+DxOB &
- *QVB(IOB+1,KOBP,JOB))+DyOB*((1.-DxOB &
- )*QVB(IOB,KOBP,JOB+1)+DxOB* &
- QVB(IOB+1,KOBP,JOB+1))))
- ! Store model surface pressure (not the error!) at T-point
- ERRF(6,N)= .001* &
- ((1.-DyOB)*((1.-DxOB)*pbase(IOB,1,JOB)+DxOB* &
- pbase(IOB+1,1,JOB))+DyOB*((1.-DxOB)* &
- pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
- #endif
- ENDDO ! END T and QV LOOP OVER OBS
- ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
- !*******************************************************
- ! USE ERRF(5,N) TO STORE KPBL AT (I,J) NEAREST THE OBS
- !*******************************************************
- DO N=1,NSTA
- IOB=MAX0(1,IA(N))
- IOB=MIN0(IOB,ide-1)
- JOB=MAX0(1,IB(N))
- JOB=MIN0(JOB,jde-1)
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
- #endif
- DXOB=RA(N)
- DYOB=RB(N)
- ERRF(5,N) = kpbl(iob+nint(dxob),job+nint(dyob))
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
- #endif
- ENDDO
- !!**********************************************************
- ENDIF ! end if(ityp.eq.3)
- ENDDO ! END BIG LOOP
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! Gather the errf values calculated by the processors with the obs.
- CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask, &
- mp_local_vobmask, mp_local_cobmask, errf)
- ! 02252010: Go ahead and assign rko for "obs-off" processors here, to
- ! fix the problem where duplicate obs can be dropped from
- ! the "obs-on" processor, but not from the others, due to
- ! rko being -99 on the "obs-off" processors.
- do n=1,nsta
- rko(n) = errf(9,n)
- enddo
- ! 02252010: End bugfix.
- #endif
- ! Print obs information.
- CALL print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
- prt_max,prt_freq,obs_prt,stnid_prt, &
- lat_prt,lon_prt,mlat_prt,mlon_prt, &
- timeob,xtime)
- ! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
- IF(INEST.EQ.1)THEN
- INPF=NPFI
- ELSE
- FNPF=IRATIO**LEVIDN(INEST)
- INPF=FNPF*NPFI
- ENDIF
- ! Gross error check for temperature. Set all vars bad.
- do n=1,nsta
- if((abs(errf(3,n)).gt.20.).and. &
- (errf(3,n).gt.-800000.))then
- errf(1,n)=-888888.
- errf(2,n)=-888888.
- errf(3,n)=-888888.
- errf(4,n)=-888888.
- varobs(1,n)=-888888.
- varobs(2,n)=-888888.
- varobs(3,n)=-888888.
- varobs(4,n)=-888888.
- endif
- enddo
- ! For printout
- ! IF(MOD(KTAU,INPF).NE.0) THEN
- ! RETURN
- ! ENDIF
- RETURN
- END SUBROUTINE errob
- SUBROUTINE upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, &
- arrin, arrout)
- !------------------------------------------------------------------------------
- ! PURPOSE: This subroutine interpolates a real 2D array defined over mass
- ! coordinate points, to U (momentum) points.
- !
- !------------------------------------------------------------------------------
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
- INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
- INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
- INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
- INTEGER, INTENT(IN) :: ids ! Starting i index for entire model domain
- INTEGER, INTENT(IN) :: ide ! Ending i index for entire model domain
- INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
- INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
- INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
- INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
- REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
- REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on U points
- ! Local variables
- integer :: i, j
- ! Do domain interior first
- do j = j_start, j_end
- do i = max(2,i_start), i_end
- arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
- enddo
- enddo
- ! Do west-east boundaries
- if(i_start .eq. ids) then
- do j = j_start, j_end
- arrout(i_start,j) = arrin(i_start,j)
- enddo
- endif
- if(i_end .eq. ide-1) then
- do j = j_start, j_end
- arrout(i_end+1,j) = arrin(i_end,j)
- enddo
- endif
- RETURN
- END SUBROUTINE upoint
- SUBROUTINE vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, &
- arrin, arrout)
- !------------------------------------------------------------------------------
- ! PURPOSE: This subroutine interpolates a real 2D array defined over mass
- ! coordinate points, to V (momentum) points.
- !
- !------------------------------------------------------------------------------
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
- INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
- INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
- INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
- INTEGER, INTENT(IN) :: jds ! Starting j index for entire model domain
- INTEGER, INTENT(IN) :: jde ! Ending j index for entire model domain
- INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
- INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
- INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
- INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
- REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
- REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on V points
- ! Local variables
- integer :: i, j
- ! Do domain interior first
- do j = max(2,j_start), j_end
- do i = i_start, i_end
- arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
- enddo
- enddo
- ! Do south-north boundaries
- if(j_start .eq. jds) then
- do i = i_start, i_end
- arrout(i,j_start) = arrin(i,j_start)
- enddo
- endif
- if(j_end .eq. jde-1) then
- do i = i_start, i_end
- arrout(i,j_end+1) = arrin(i,j_end)
- enddo
- endif
- RETURN
- END SUBROUTINE vpoint
- LOGICAL FUNCTION TILE_MASK(iloc, jloc, its, ite, jts, jte)
- !------------------------------------------------------------------------------
- ! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
- !
- ! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
- ! tile-range indices (its,jts) and (ite,jte)
- ! FALSE otherwise.
- !
- !------------------------------------------------------------------------------
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: iloc
- INTEGER, INTENT(IN) :: jloc
- INTEGER, INTENT(IN) :: its
- INTEGER, INTENT(IN) :: ite
- INTEGER, INTENT(IN) :: jts
- INTEGER, INTENT(IN) :: jte
- ! Local variables
- LOGICAL :: retval
- TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
- jloc .LE. jte .AND. jloc .GE. jts )
- RETURN
- END FUNCTION TILE_MASK
- !-----------------------------------------------------------------------
- SUBROUTINE nudob(j, ivar, aten, inest, ifrest, ktau, ktaur, &
- xtime, mu, msfx, msfy, nndgv, nerrf, niobf, maxdom, &
- npfi, ionf, rinxy, twindo, &
- nudge_pbl, &
- sfcfact, sfcfacr, &
- levidn, &
- parid, nstat, &
- rinfmn, rinfmx, pfree, dcon, tfaci, &
- sfc_scheme_horiz, sfc_scheme_vert, maxsnd_gap, &
- lev_in_ob, plfo, nlevs_ob, &
- iratio, dx, dtmin, rio, rjo, rko, &
- timeob, varobs, errf, pbase, ptop, pp, &
- iswind, istemp, ismois, giv, git, giq, &
- savwt, kpblt, nscan, &
- vih1, vih2, terrh, zslab, &
- iprt, &
- ids,ide, jds,jde, kds,kde, & ! domain dims
- ims,ime, jms,jme, kms,kme, & ! memory dims
- its,ite, jts,jte, kts,kte ) ! tile dims
- !-----------------------------------------------------------------------
- USE module_model_constants
- USE module_domain
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !
- ! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
- ! VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
- ! ASSIMILATION FROM INDIVIDUAL OBSERVATIONS. THE NUDGING
- ! TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
- ! WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
- ! ANALYSIS. THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
- ! AND MINIMAL STORAGE REQUIREMENTS. ALGORITHMS SHOULD BE
- ! VECTORIZED WHEREVER POSSIBLE.
- !
- ! HISTORY: Original author: MM5 version???
- ! 02/04/2004 - Creation of WRF version. Al Bourgeois
- ! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
- !------------------------------------------------------------------------------
- !
- ! 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 a "J-slab" here is west-east in
- ! extent, not south-north as for MM5. -ajb 06/10/2004
- !
- ! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
- ! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
- ! TYPES OF FACTORS:
- ! 1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
- ! TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
- ! TIME (XTIME) ARE USED. OBSERVATIONS CLOSEST TO
- ! XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
- ! 2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
- ! CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
- ! (RINSIG).
- ! 3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
- ! CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY). THE
- ! VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
- ! TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
- !
- ! THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
- ! VALUE OF IVAR AS FOLLOWS:
- ! IVAR VARIABLE(TAU-1)
- ! ---- ---------------
- ! 1 U
- ! 2 V
- ! 3 T
- ! 4 QV
- ! 5 PSB(CROSS) REMOVED IN V3
- ! (6) PSB(DOT)
- !
- !-----------------------------------------------------------------------
- !
- ! Description of input arguments.
- !
- !-----------------------------------------------------------------------
- INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
- INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
- INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
- INTEGER, INTENT(IN) :: j ! south-north running coordinate.
- INTEGER, INTENT(IN) :: ivar
- INTEGER, INTENT(IN) :: inest ! domain index
- LOGICAL, INTENT(IN) :: ifrest
- INTEGER, INTENT(IN) :: ktau
- INTEGER, INTENT(IN) :: ktaur
- REAL, INTENT(IN) :: xtime ! forecast time in minutes
- INTEGER, INTENT(IN) :: nndgv ! number of nudge variables
- INTEGER, INTENT(IN) :: nerrf ! number of error fields
- INTEGER, INTENT(IN) :: niobf ! number of observations
- INTEGER, INTENT(IN) :: maxdom ! maximum number of domains
- INTEGER, INTENT(IN) :: npfi
- INTEGER, INTENT(IN) :: ionf
- REAL, INTENT(IN) :: rinxy
- REAL, INTENT(IN) :: twindo
- REAL, intent(in) :: sfcfact ! scale for time window (surface obs)
- REAL, intent(in) :: sfcfacr ! scale for hor rad inf (surface obs)
- LOGICAL, intent(in) :: nudge_pbl ! flag for nudging in pbl
- INTEGER, INTENT(IN) :: levidn(maxdom) ! level of nest
- INTEGER, INTENT(IN) :: parid(maxdom) ! parent domain id
- INTEGER, INTENT(IN) :: nstat ! number of obs stations
- REAL, INTENT(IN) :: rinfmn ! minimum radius of influence
- REAL, INTENT(IN) :: rinfmx ! maximum radius of influence
- REAL, INTENT(IN) :: pfree ! pressure level (cb) where terrain effect becomes small
- REAL, INTENT(IN) :: dcon ! 1/DPSMX
- REAL, INTENT(IN) :: tfaci ! scale factor used for ramp-down in dynamic initialization
- INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
- INTEGER , INTENT(IN) :: sfc_scheme_vert ! vertical spreading scheme for surf obs (orig or regime vif)
- REAL , INTENT(IN) :: maxsnd_gap ! max allowed pressure gap in soundings, for interp (centibars)
- REAL, INTENT(IN) :: lev_in_ob(niobf) ! Level in sounding-type obs.
- REAL, intent(IN) :: plfo(niobf) ! Index for type of obs platform
- REAL, INTENT(IN) :: nlevs_ob(niobf) ! Number of levels in sounding.
- INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
- REAL, INTENT(IN) :: dx ! This domain grid cell-size (m)
- REAL, INTENT(IN) :: dtmin ! Model time step in minutes
- REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
- REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
- REAL, INTENT(INOUT) :: rko(niobf) ! Obs vertical coordinate.
- REAL, INTENT(IN) :: timeob(niobf)
- REAL, INTENT(IN) :: varobs(nndgv,niobf)
- REAL, INTENT(IN) :: errf(nerrf, niobf)
- REAL, INTENT(IN) :: pbase( ims:ime, kms:kme ) ! Base pressure.
- REAL, INTENT(IN) :: ptop
- REAL, INTENT(IN) :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
- REAL, INTENT(IN) :: mu(ims:ime) ! Air mass on u, v, or mass-grid
- REAL, INTENT(IN) :: msfx(ims:ime) ! Map scale (only used for vars u & v)
- REAL, INTENT(IN) :: msfy(ims:ime) ! Map scale (only used for vars u & v)
- INTEGER, intent(in) :: iswind ! Nudge flag for wind
- INTEGER, intent(in) :: istemp ! Nudge flag for temperature
- INTEGER, intent(in) :: ismois ! Nudge flag for moisture
- REAL, intent(in) :: giv ! Coefficient for wind
- REAL, intent(in) :: git ! Coefficient for temperature
- REAL, intent(in) :: giq ! Coefficient for moisture
- REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
- REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
- INTEGER, INTENT(IN) :: kpblt(ims:ime)
- INTEGER, INTENT(IN) :: nscan ! number of scans
- REAL, INTENT(IN) :: vih1(its:ite) ! Vert infl ht (m) abv grd for full wts
- REAL, INTENT(IN) :: vih2(its:ite) ! Vert infl ht (m) abv grd for ramp
- REAL, INTENT(IN) :: terrh(ims:ime) ! Terrain height (m)
- ! INTEGER, INTENT(IN) :: vik1(its:ite) ! Vertical infl k-level for full wts
- ! INTEGER, INTENT(IN) :: vik2(its:ite) ! Vertical infl k-level for ramp
- REAL, INTENT(IN) :: zslab(ims:ime, kms:kme) ! model ht above ground (m)
- LOGICAL, INTENT(IN) :: iprt ! print flag
- ! Local variables
- integer :: mm(maxdom)
- integer :: kobs ! k-lev below obs (for obs straddling pblt)
- integer :: kpbl_obs(nstat) ! kpbl at grid point (IOB,JOB) (ajb 20090519)
- real :: ra(niobf)
- real :: rb(niobf)
- real :: psurf(niobf)
- real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
- real :: rscale(ims:ime) ! For converting to rho-coupled units.
- real :: wtij(ims:ime) ! For holding weights in i-loop
- real :: reserf(100)
- character*40 name
- character*3 chr_hr
- character(len=200) :: msg ! Argument to wrf_message
- !*** DECLARATIONS FOR IMPLICIT NONE
- integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
- integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
- integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
- integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
- integer :: komin,komax,nn,nhi,nlo,nnjc
- integer :: i_s,i_e
- integer :: istq
- real :: gfactor,rfactor,gridx,gridy,rindx,ris
- real :: grfacx,grfacy
- real :: timewt,pob
- real :: ri,rj,rx,ry,rsq,pdfac,erfivr,dk,slope,rinfac
- real :: dprim,dsq,d ! vars for mm5 surface-ob weighting
- real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
- real :: dz_ramp ! For ramping weights for surface obs
- real :: scratch
- integer :: kk !ajb temp
- ! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
- ! if(ivar.ne.4)return
- !yliu start -- for multi-scans: NSCAN=0: original
- ! NSCAN=1: added a scan with a larger Ri and smaller G
- ! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
- ! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
- if(NSCAN.ne.0) then
- IF (iprt) then
- write(msg,*) 'SAVWT must be resized for NSCAN=1'
- call wrf_message(msg)
- ENDIF
- call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
- endif
- IPLO=0 + NSCAN*4
- GFACTOR=1. + NSCAN*(-1. + 0.33333)
- RFACTOR=1. + NSCAN*(-1. + 3.0)
- !yliu end
- ! jc
- ! return if too close to j boundary
- if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
- ! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
- ! $ ' too close to boundary.'
- return
- endif
- if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
- ! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
- ! $ ' too close to boundary.'
- return
- endif
- ! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
- ICUT=0
- IF(INEST.GT.1)ICUT=1
- i_s = max0(2+icut,its)
- i_e = min0(ide-1-icut,ite)
- IPL=IVAR + IPLO !yliu +IPLO
- ! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
- INPF=(IRATIO**LEVIDN(INEST))*NPFI
- INFR=(IRATIO**LEVIDN(INEST))*IONF
- GRIDX=0.0
- GRIDY=0.0
- IGRID=0
- IF(IVAR.GE.3)THEN
- GRIDX=0.5
- GRIDY=0.5
- IGRID=1
- ELSEIF(IVAR.eq.1) THEN
- GRIDY=0.5
- GRIDX=0.0
- IGRID=1
- ELSEIF(IVAR.eq.2) THEN
- GRIDX=0.5
- GRIDY=0.0
- IGRID=1
- ENDIF
- ! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
- ! KILOMETERS TO GRID LENGTHS, RINDX
- RINDX=RINXY*1000./DX * RFACTOR !yliu *RFACTOR
- RIS=RINDX*RINDX
- IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
- IF(MOD(KTAU,INFR).NE.0)GOTO 126
- 5 CONTINUE
- IF (iprt) THEN
- IF(J.EQ.10) then
- write(msg,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
- call wrf_message(msg)
- ENDIF
- ENDIF
- 6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
- 'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
- ' rindx=',f4.1)
- !********************************************************************
- !ajb_07012008 Setting ra and rb for the fine mesh her is now simple.
- ! Values are no longer calculated here based on the
- ! coarse mesh, since direct use of WRF map projections
- ! on each nest was implemented in subroutine in4dob.
- !********************************************************************
- ! SET RA AND RB
- DO N=1,NSTAT
- RA(N)=RIO(N)-GRIDX
- RB(N)=RJO(N)-GRIDY
- ENDDO
- ! INITIALIZE WEIGHTING ARRAYS TO ZERO
- DO I=its,ite
- DO K=1,kte
- WT(I,K)=0.0
- WT2ERR(I,K)=0.0
- ENDDO
- ENDDO
- ! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
- ! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
- !
- ! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
- ! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
- ! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
- ! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
- ! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
- ! ajb 05052009: psurf is actually pbase(k=1) interpolated to obs (i,j).
- DO N=1,NSTAT
- IF(IVAR.GE.3)THEN
- PSURF(N)=ERRF(6,N)
- ELSE
- IF(IVAR.EQ.1)THEN
- PSURF(N)=ERRF(7,N) ! U-points
- ELSE
- PSURF(N)=ERRF(8,N) ! V-points
- ENDIF
- ENDIF
- ENDDO
- ! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
- ! J-STRIP
- MAXJ=J+IFIX(RINDX*RINFMX+0.99) !ajb
- MINJ=J-IFIX(RINDX*RINFMX+0.99) !ajb
- ! jc comment out this? want to use obs beyond the domain?
- ! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
- ! MINJ=MAX0(1,MINJ) !yliu
- n=1
- !***********************************************************************
- DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
- !***********************************************************************
- ! Soundings are consecutive obs, but they need to be treated as a single
- ! entity. Thus change the looping to nnn, with n defined separately.
- !yliu
- ! note for sfc data: nsndlev=1 and njcsnd=1
- nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
- ! yliu start -- set together with the other parts
- ! test: do the sounding levels as individual obs
- ! nsndlev=1
- ! yliu end
- njcsnd=nsndlev
- ! set pob here, to be used later
- pob=varobs(5,n)
- ! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
- ! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
- ! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
- ! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
- ! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
- ! ATION.
- !yliu: Skip bad obs if it is sfc or single level sounding.
- !yliu: Before this (020918), a snd will be skipped if its first
- !yliu level has bad data- often true due to elevation
- IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
- ! print *, " bad obs skipped"
- ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
- ! print *, " skipped obs far away from this J-slice"
- !----------------------------------------------------------------------
- ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
- !----------------------------------------------------------------------
- ! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
- ! FOR THE VERTICAL WEIGHTING, WTSIG
- ! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
- !ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
- !ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- rko(n) = errf(9,n) !ajb 20021210
- kpbl_obs(n) = errf(5,n) !ajb 20090519
- #endif
- !ajb 20090427 added .45 to rko so KOB is equal to 1 only if RKO > 1.05
- KOB=nint(RKO(N)+0.45)
- KOB=MIN0(kte,KOB)
- KOB=MAX0(1,KOB)
- ! ASSIMILATE SURFACE LAYER DATA ON SIGMA
- IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
- ! Compute temporal weight
- timewt = get_timewt(xtime,dtmin,twindo,sfcfact,timeob(n))
- DO K=1,kte
- WTSIG(K)=0.0
- ENDDO
- ! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
- ! WTSIG(1)=1.0
- ! WTSIG(2)=0.67
- ! WTSIG(3)=0.33
- ! KOMIN=3
- ! KOMAX=1
- ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
- ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
- ! fix this because kpblt at 1 and il is 0
- MAXI=IFIX(RA(N)+0.99+RINDX*sfcfacr)
- MAXI=MIN0(ide-1,MAXI)
- MINI=IFIX(RA(N)-RINDX*sfcfacr-0.99)
- MINI=MAX0(2,MINI)
- !yliu start
- ! use also obs outside of this domain -- surface obs
- ! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
- ! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
- ! print *, " skipped obs far away from this domain"
- ! currently can use obs within this domain or ones very close to (1/3
- ! influence of radius in the coarse domain) this
- ! domain. In later case, use BC column value to approximate the model value
- ! at obs point -- ERRF need model field in errob.F !!
- if ( RA(N).GE.(0.-RINDX*sfcfacr/3) &
- .and. RA(N).LE.float(ide)+RINDX*sfcfacr/3 &
- .and. RB(N).GE.(0.-RINDX*sfcfacr/3) &
- .and. RB(N).LE.float(jde)+RINDX*sfcfacr/3) then
- ! or use obs within this domain only
- ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
- ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
- ! print *, " skipped obs far outside of this domain"
- ! if(j.eq.3 .and. ivar.eq.3) then
- ! write(6,*) 'N = ',n,' exit 120 3'
- ! endif
- !yliu end
- !
- ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
- ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
- ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
- RJ=FLOAT(J)
- RX=RJ-RB(N)
- ! WEIGHTS FOR THE 3-D VARIABLES
- ERFIVR=ERRF(IVAR,N)
-
- !ajb Compute and add weights to sum only if nudge_pbl switch is on.
- if(nudge_pbl) then
- ! Apply selected horizontal spreading scheme.
- if(SFC_SCHEME_HORIZ.eq.1) then ! original mm5 scheme
- DO I=max0(its,MINI),min0(ite,MAXI)
- RI=FLOAT(I)
- RY=RI-RA(N)
- RIS=RINDX*RINDX*sfcfacr*sfcfacr
- RSQ=RX*RX+RY*RY
- ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
- DPRIM=SQRT(RSQ)
- D=DPRIM+RINDX*DCON*ABS(psurf(n)-.001*pbase(i,1))
- DSQ=D*D
- WTIJ(i)=(RIS-DSQ)/(RIS+DSQ)
- WTIJ(i)=AMAX1(0.0,WTIJ(i))
- ENDDO
- else ! wrf scheme
- DO I=max0(its,MINI),min0(ite,MAXI)
- RI=FLOAT(I)
- RY=RI-RA(N)
- RIS=RINDX*RINDX*sfcfacr*sfcfacr
- RSQ=RX*RX+RY*RY
- ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
- wtij(i)=(ris-rsq)/(ris+rsq)
- scratch = (abs (psurf(n)-.001*pbase(i,1))*DCON)
- pdfac=1.-AMIN1(1.0,scratch)
- wtij(i)=wtij(i)*pdfac
- WTIJ(i)=AMAX1(0.0,WTIJ(i))
- ENDDO
- endif
- ! Apply selected vertical spreading scheme.
- if(SFC_SCHEME_VERT.eq.1) then ! original simple scheme
- DO I=max0(its,MINI),min0(ite,MAXI)
- ! try making sfc obs weighting go thru pbl
- komax=max0(3,kpblt(i))
- IF (iprt) THEN
- if (kpblt(i).gt.25 .and. ktau.ne.0) &
- write(6,552)inest,i,j,kpblt(i)
- 552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
- ENDIF
- if(kpblt(i).gt.25) komax=3
- komin=1
- dk=float(komax)
- do k=komin,komax
- wtsig(k)=float(komax-k+1)/dk
- WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
- WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K) &
- *WTSIG(K)*ERFIVR
- enddo
- ENDDO
- else ! regime-based vif scheme
- ! Here we calculate weights in the vertical coordinate, based on vih1 and vih2.
- ! In the equation for wtsig(k), Z=zslab(i,k)-terrh(i) contains the model Z-values
- ! (height above ground in meters) on a J-slab. The equation produces wtsig = 1.0 at
- ! levels 1 to K, where z(K) < vih1 < z(K+1). For the example below, in the equation
- ! for wtsig(k), the expression vih1(i)-Z(i,k) is strictly positive for k=1,2,3 since
- ! levels 1, 2, and 3 are below vih1. So xtsig(k)=min(1.0, 1.0-x) where x > 0 ==>
- ! wtsig(k)=1 for k=1,2,3.
- !
- ! For levels K+1 and up, wtsig will decrease linearly with height, with values
- ! along the ramp that has value 1.0 at vih1 and 0.0 at vih2. In the example:
- !
- ! dz_ramp = 1/(200-150) = 1/50
- ! xtsig(4) = 1 + (150-175)/50 = 1 - 1/2 = 1/2
- !
- ! WTSIG
- ! 1 -|* * * * * *
- ! |
- ! | *
- ! |
- ! | *
- ! |
- ! | *
- ! 0 -|--|-------|-----------|------|----|----|---------|----> Z = HT ABOVE
- ! 15 55 115 150 175 200 250 GROUND
- ! k=1 k=2 k=3 vih1 k=4 vih2 k=5
-
- DO I=max0(its,MINI),min0(ite,MAXI)
- dz_ramp = 1.0 / max( 1.0, vih2(i)-vih1(i) ) ! vih2 >= vih1 by construct
- LML: do k = kts, kte
- wtsig(k) = min( 1.0, 1.0 + ( vih1(i)-zslab(i,k)+terrh(i) ) * dz_ramp )
- wtsig(k) = max( 0.0, wtsig(k))
- if(wtsig(k).le.0.0) EXIT LML
- WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
- WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K) &
- *WTSIG(K)*ERFIVR
- enddo LML
- ENDDO
- endif
- endif ! end if nudge-pbl switch is on
- endif ! end check for obs in domain
- ! END SURFACE-LAYER OBS NUDGING
- ELSE
- ! Compute temporal weight
- timewt = get_timewt(xtime,dtmin,twindo,1.,timeob(n))
- ! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
- !
- ! print *,'in upper air section'
- ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
- ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
- ! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
- ! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
- !ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
- slope = (RINFMN-RINFMX)/(psurf(n)-PFREE)
- RINFAC=SLOPE*POB+RINFMX-SLOPE*pfree
- RINFAC=AMAX1(RINFAC,RINFMN)
- RINFAC=AMIN1(RINFAC,RINFMX)
- !yliu: for multilevel upper-air data, take the maximum
- ! for the I loop.
- if(nsndlev.gt.1) RINFAC = RINFMX
- !yliu end
- MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
- MAXI=MIN0(ide-IGRID,MAXI)
- MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
- MINI=MAX0(1,MINI)
- ! yliu start
- ! use also obs outside of but close to this domain -- upr data
- ! if( RA(N).LT.(0.-RINFAC*RINDX)
- ! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
- ! & .or. RB(N).LT.(0.-RINFAC*RINDX)
- ! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
- ! print *, " skipped obs far away from this I-range"
- ! currently can use obs within this domain or ones very close to (1/3
- ! influence of radius in the coarse domain) this
- ! domain. In later case, use BC column value to approximate the model value
- ! at obs point -- ERRF need model field in errob.F !!
- !cc if (i.eq.39 .and. j.eq.34) then
- !cc write(6,*) 'RA(N) = ',ra(n)
- !cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
- !cc endif
- if( RA(N).GE.(0.-RINFAC*RINDX/3) &
- .and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
- .and. RB(N).GE.(0.-RINFAC*RINDX/3) &
- .and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
- ! or use obs within this domain only
- ! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
- ! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
- ! print *, " skipped obs far outside of this domain"
- ! yliu end
- ! is this 2 needed here - kpbl not used?
- ! MINI=MAX0(2,MINI)
- ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
- ! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
- ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
- RJ=FLOAT(J)
- RX=RJ-RB(N)
- ! WEIGHTS FOR THE 3-D VARIABLES
- !
- ERFIVR=ERRF(IVAR,N)
- ! jc
- nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
- ! yliu start
- ! test: do the sounding levels as individual obs
- ! nsndlev=1
- ! yliu end
- njcsnd=nsndlev
- !
- DO I=max0(its,MINI),min0(ite,MAXI)
- ! jc
- RI=FLOAT(I)
- RY=RI-RA(N)
- RIS=RINDX*RINFAC*RINDX*RINFAC
- RSQ=RX*RX+RY*RY
- ! yliu test: for upper-air data, keep D1 influence radii
- WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
- WTIJ=AMAX1(0.0,WTIJ(i))
- ! weight ob in vertical with +- 50 mb
- ! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
- if(nsndlev.eq.1) then
- rinprs=7.5
- !
- else
- rinprs=3.0
- endif
- ! yliu end
- !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
- ! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
- !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
- if(nsndlev.eq.1)then
- !----------------------------------------------------------------------
- ! --- HANDLE 1-LEVEL OBSERVATIONS ---
- !----------------------------------------------------------------------
- ! if(I.eq.MINI) print *, " Single snd "
- ! ERFIVR is the residual (difference) between the ob and the model
- ! at that point. We can analyze that residual up and down.
- ! First find komin for ob.
- !yliu start -- in the old code, komax and komin were reversed!
- do k=kte,1,-1
- pijk = .001*(pbase(i,k)+pp(i,k))
- ! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
- if(pijk.ge.(pob+rinprs)) then
- komin=k
- go to 325
- endif
- enddo
- komin=1
- 325 continue
- ! now find komax for ob
- do k=3,kte
- pijk = .001*(pbase(i,k)+pp(i,k))
- if(pijk.le.(pob-rinprs)) then
- komax=k
- go to 326
- endif
- enddo
- komax=kte ! yliu 20050706
- 326 continue
- ! yliu: single-level upper-air data will act either above or below the PBL top
- ! ajb: Reset komin or komax. if kobs>kpbl_obs, komin=kpbl_obs+1, else komax=kpbl_obs
- if( (kpblt(i).le.komax) .and. (kpblt(i).ge.komin) ) then
- kobs = 1
- OBS_K: do k = komin, komax
- if( pob .gt. .001*(pbase(i,k)+pp(i,k)) ) then
- kobs = k
- EXIT OBS_K
- endif
- enddo OBS_K
- if(kobs.gt.kpbl_obs(n)) then
- ! Obs will act only above the PBL top
- komin=max0(kobs, komin) ! kobs here is kpblt(i)+1
- else ! Obs acts below PBL top
- ! Obs will act only below the PBL top
- komax=min0(kpblt(i), komax)
- endif
- endif
- ! yliu end
- !
- ! print *,'1 level, komin,komax=',komin,komax
- ! if(i.eq.MINI) then
- ! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
- ! ERFIVR=0
- ! endif
- do k=1,kte
- reserf(k)=0.0
- wtsig(k)=0.0
- enddo
- !yliu end
- !ajb Add weights to sum only if nudge_pbl switch is on OR obs is above pbl top.
- if(nudge_pbl .or. komin.ge.kpblt(i)) then
- do k=komin,komax
- pijk = .001*(pbase(i,k)+pp(i,k))
- reserf(k)=erfivr
- wtsig(k)=1.-abs(pijk-pob)/rinprs
- wtsig(k)=amax1(wtsig(k),0.0)
- ! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
- ! Now calculate WT and WT2ERR for each i,j,k point cajb
- WT(I,K)=WT(I,K)+TIMEWT*WTIJ(i)*wtsig(k)
- WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)* &
- reserf(k)*wtsig(k)*wtsig(k)
- enddo
- endif
- else
- !----------------------------------------------------------------------
- ! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
- !----------------------------------------------------------------------
- !yliu start
- ! if(I.eq.MINI) print *, " Multi-level snd "
- ! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
- ! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
- if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
- IF (iprt) THEN
- write(msg,*) "n = ",n,"nsndlev = ",nsndlev
- call wrf_message(msg)
- write(msg,*) "nlevs_ob,lev_in_ob", &
- nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
- call wrf_message(msg)
- call wrf_message("in nudobs.F: sounding level messed up, stopping")
- ENDIF
- call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
- endif
- !yliu end
- ! This is for a multi-level observation
- ! The trick here is that the sounding is "one ob". You don't
- ! want multiple levels to each be treated like separate
- ! and independent observations.
- ! At each i,j want to interpolate sounding to the model levels at that
- ! particular point.
- komin=1
- komax=kte-2
- ! this loop goes to 1501
- ! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
- !yliu start
- do k=1,kte
- reserf(k)=0.0
- wtsig(k)=0.0
- enddo
- !yliu end
- do k=komax,komin,-1
-
- pijk = .001*(pbase(i,k)+pp(i,k))
- ! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
- if(pijk.gt.varobs(5,n)) then
- go to 1501
- endif
- ! if sigma level pressure is .lt. than the highest ob level, don't interpolate
- if(pijk.le.varobs(5,n+nsndlev-1)) then
- go to 1501
- endif
- ! now interpolate sounding to this k
- ! yliu start-- recalculate WTij for each k-level
- !ajb SLOPE = (RINFMN-RINFMX)/(pdoc(i,j)+PTOP-PFREE)
- slope = (RINFMN-RINFMX)/ (.001*pbase(i,1)-PFREE)
- RINFAC=SLOPE*pijk+RINFMX-SLOPE*PFREE
- RINFAC=AMAX1(RINFAC,RINFMN)
- RINFAC=AMIN1(RINFAC,RINFMX)
- RIS=RINDX*RINFAC*RINDX*RINFAC
- RSQ=RX*RX+RY*RY
- ! for upper-air data, keep D1 influence radii
- WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
- WTIJ(i)=AMAX1(0.0,WTIJ(i))
- ! yliu end
- ! this loop goes to 1503
- do nn=2,nsndlev
- ! only set pobhi if varobs(ivar) is ok
- pobhi=-888888.
- if(varobs(ivar,n+nn-1).gt.-800000. &
- .and. varobs(5,n+nn-1).gt.-800000.) then
- pobhi=varobs(5,n+nn-1)
- nhi=n+nn-1
- if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.0.5*maxsnd_gap) then
- go to 1502 ! within maxsnd_gap/2 of obs height
- endif
- endif
- enddo
- ! did not find any ob above within maxsnd_gap/2, so jump out
- go to 1501
- 1502 continue
- nlo=nhi-1
- do nnjc=nhi-1,n,-1
- if(varobs(ivar,nnjc).gt.-800000. &
- .and. varobs(5,nnjc).gt.-800000.) then
- poblo=varobs(5,nnjc)
- nlo=nnjc
- if(poblo.gt.pijk .and. abs(poblo-pijk).lt.0.5*maxsnd_gap) then
- go to 1505 ! within maxsnd_gap/2 of obs height
- endif
- endif
- enddo
- !yliu end --
- ! did not find any ob below within maxsnd_gap/2, so jump out
- go to 1501
- 1505 continue
- ! interpolate to model level
- pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
- reserf(k)=errf(ivar,nlo)+ &
- (errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
- wtsig(k)=1.
-
- 1501 continue
- ! now calculate WT and WT2ERR for each i,j,k point cajb
- !ajb Add weights to sum only if nudge_pbl switch is on OR k > kpblt.
- if(nudge_pbl .or. k.gt.kpblt(i)) then
- WT(I,K)=WT(I,K)+TIMEWT*WTIJ(i)*wtsig(k)
-
- WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)* &
- reserf(k)*wtsig(k)*wtsig(k)
- endif
- enddo ! enddo k levels
- ! end multi-levels
- endif ! end if(nsndlev.eq.1)
- !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
- ! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
- !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
- !
- ENDDO ! END DO MINI,MAXI LOOP
- endif ! check for obs in domain
- ! END OF NUDGING TO OBS ON PRESSURE LEVELS
- ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
- !----------------------------------------------------------------------
- ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
- !----------------------------------------------------------------------
- ! n=n+1
- n=n+njcsnd
- !yliu 1202 continue
- if(n.gt.nstat)then
- ! print *,'n,nstat=',n,nstat,ivar,j
- go to 1203
- endif
- ! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
- !***********************************************************************
- ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
- !***********************************************************************
- 1203 continue
- ! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
- ! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
- ! THE ATEN ARRAY
- ! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
- ! THEY ARE USED BELOW IN THE DENOMINATOR.
- DO K=kts,kte
- DO I=its,ite
- IF(WT(I,K).EQ.0)THEN
- WT2ERR(I,K)=0.0
- ENDIF
- IF(WT(I,K).EQ.0)THEN
- WT(I,K)=1.0
- ENDIF
- ENDDO
- ENDDO
- 126 CONTINUE
- IF(IVAR.GE.3)GOTO 170
- ! this is for u,v
- ! 3-D DOT POINT TENDENCIES
-
- ! Calculate scales for converting nudge factor from u (v)
- ! to rho_u (or rho_v) units.
- IF (IVAR == 1) THEN
- call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
- ELSE IF (IVAR == 2) THEN
- call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
- END IF
-
- DO K=1,kte
- DO I=i_s,i_e
- IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
- W2EOWT=WT2ERR(I,K)/WT(I,K)
- ELSE
- W2EOWT=SAVWT(IPL,I,K)
- ! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,W2EOWT = ',i,j,k,ipl,W2EOWT
- ENDIF
- ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
- ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
- ! write(6,*) 'ATEN calc: k = ',k
- ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
- ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
- ! $ ' W2EOWT = ',w2eowt
- ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
- ! $ ' GFACTOR = ',gfactor
- ! endif
- !
- ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
- ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
- ! write(6,*) 'ATEN calc: k = ',k
- ! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
- ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
- ! $ ' W2EOWT = ',w2eowt
- ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
- ! $ ' GFACTOR = ',gfactor
- ! endif
- ! if(ivar .eq. 1 .and. i.eq.38 .and. (j.ge.36.and.j.le.62) .and. k.eq.7) then
- ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
- ! write(6,*)
- ! write(6,*) 'ATEN calc: j = ',j,' k = ',k
- ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
- ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),' W2EOWT = ',w2eowt
- ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,' GFACTOR = ',gfactor
- ! endif
- ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
- *W2EOWT*TFACI &
- *ISWIND *GFACTOR !yliu *GFACTOR
- ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
- ! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
- ! endif
- ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
- ! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
- ! endif
- ENDDO
- ENDDO
- IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
- DO K=1,kte
- DO I=its,ite
- SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
- ! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,savwt = ',i,j,k,ipl,savwt(ipl,i,k)
- ENDDO
- ENDDO
- ENDIF
- RETURN
- 170 CONTINUE
- ! 3-D CROSS-POINT TENDENCIES
- ! this is for t (ivar=3) and q (ivsr=4)
- IF(3-IVAR.LT.0)THEN
- GITQ=GIQ
- ELSE
- GITQ=GIT
- ENDIF
- IF(3-IVAR.LT.0)THEN
- ISTQ=ISMOIS
- ELSE
- ISTQ=ISTEMP
- ENDIF
- DO K=1,kte
- DO I=i_s,i_e
- IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
- W2EOWT=WT2ERR(I,K)/WT(I,K)
- ELSE
- W2EOWT=SAVWT(IPL,I,K)
- ENDIF
- ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
- ! scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
- ! write(6,*) 'ATEN calc: k = ',k
- ! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
- ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i), &
- ! ' W2EOWT = ',w2eowt
- ! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq, &
- ! ' GFACTOR = ',gfactor
- ! endif
-
- ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
- ! scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
- ! write(6,*) 'ATEN calc: k = ',k
- ! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
- ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
- ! $ ' W2EOWT = ',w2eowt
- ! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq,
- ! $ ' GFACTOR = ',gfactor
- ! endif
- ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
- *W2EOWT*TFACI*ISTQ *GFACTOR !yliu *GFACTOR
- ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
- ! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
- ! endif
- ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
- ! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
- ! endif
- ENDDO
- ENDDO
- IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
- DO K=1,kte
- DO I=its,ite
- SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
- ENDDO
- ENDDO
- ENDIF
- RETURN
- END SUBROUTINE nudob
- SUBROUTINE date_string(year, month, day, hour, minute, second, cdate)
- !-----------------------------------------------------------------------
- ! PURPOSE: Form a date string (YYYY-MM-DD_hh:mm:ss) from integer
- ! components.
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- INTEGER, INTENT(IN) :: year
- INTEGER, INTENT(IN) :: month
- INTEGER, INTENT(IN) :: day
- INTEGER, INTENT(IN) :: hour
- INTEGER, INTENT(IN) :: minute
- INTEGER, INTENT(IN) :: second
- CHARACTER*19, INTENT(INOUT) :: cdate
- ! Local variables
- integer :: ic ! loop counter
- cdate(1:19) = "0000-00-00_00:00:00"
- write(cdate( 1: 4),'(i4)') year
- write(cdate( 6: 7),'(i2)') month
- write(cdate( 9:10),'(i2)') day
- write(cdate(12:13),'(i2)') hour
- write(cdate(15:16),'(i2)') minute
- write(cdate(18:19),'(i2)') second
- do ic = 1,19
- if(cdate(ic:ic) .eq. " ") cdate(ic:ic) = "0"
- enddo
- RETURN
- END SUBROUTINE date_string
- SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
- INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
- REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
- REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
- REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
- ! Local variables
- integer :: i
- ! Calculate scales to be used for producing rho-coupled nudging factors.
- do i = its,ite
- rscale(i) = a(i)/msf(i)
- enddo
- RETURN
- END SUBROUTINE calc_rcouple_scales
- SUBROUTINE print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
- prt_max,prt_freq,obs,stnid,lat,lon, &
- mlat,mlon,timeob,xtime)
- !*************************************************************************
- ! Purpose: Print obs information.
- !*************************************************************************
- IMPLICIT NONE
- LOGICAL, intent(in) :: iprt ! Print flag
- INTEGER, intent(in) :: inest ! Nest level
- 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, intent(in) :: rko(niobf) ! Bottom-top north coord (non-stagger)
- INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
- INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
- INTEGER, intent(in) :: obs(prt_max) ! Saved obs indices to print
- INTEGER, intent(in) :: stnid(40,prt_max) ! Saved station ids
- REAL, intent(in) :: lat(prt_max) ! Saved latitudes
- REAL, intent(in) :: lon(prt_max) ! Saved longitudes
- REAL, intent(in) :: mlat(prt_max) ! Saved model latitudes
- REAL, intent(in) :: mlon(prt_max) ! Saved longitudes
- REAL, intent(in) :: timeob(niobf) ! Times of each observation (hours)
- REAL, intent(in) :: xtime ! Model time in minutes
- ! Local variables
- integer :: i ! Loop counter over obs station chars
- integer :: n ! Loop counter over obs
- integer :: pnx ! Obs index for printout
- character(len=200) :: msg ! Argument to wrf_message
- character(len=20) :: station_id ! Station id of observation
- if(iprt) then
- if(prt_max.gt.0) then
- if(obs(1).ne.-999) then
- call wrf_message("")
- write(msg,fmt='(a,i4,a,f8.1,a)') 'REPORTING OBS MASS-PT LOCS FOR NEST ', &
- inest,' AT XTIME=',xtime,' MINUTES'
- call wrf_message(msg)
- write(msg,fmt='(a,i4,a,i5,a)') 'FREQ=',prt_freq,', MAX=',prt_max, &
- ' LOCS, NEWLY READ OBS ONLY, -999 => OBS OFF PROC'
- call wrf_message(msg)
- call wrf_message("")
- write(msg,fmt='(3a)') ' OBS# I J K OBS LAT', &
- ' OBS LON XLAT(I,J) XLONG(I,J) TIME(hrs)', &
- ' OBS STATION ID'
- call wrf_message(msg)
- endif
- endif
- ! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
- ! Hence subtract .5 from each to get mass-point coords.
- do n=1,prt_max
- pnx = obs(n)
- if(pnx.ne.-999) then
- ! Retrieve 15 chars of station id
- do i = 1,15
- station_id(i:i) = char(stnid(i,n))
- enddo
- write(msg,fmt='(2x,i7,3f8.3,2f9.3,2x,f9.3,2x,f9.3,3x,f6.2,7x,a15)') &
- pnx,rio(pnx)-.5,rjo(pnx)-.5,rko(pnx),lat(n),lon(n), &
- mlat(n),mlon(n),timeob(pnx),station_id
- call wrf_message(msg)
- endif
- enddo
- if(obs(1).ne.-999) call wrf_message("")
- endif
- END SUBROUTINE print_obs_info
- REAL FUNCTION ht_to_p( h, pbbc, ppbc, z, ic, jc, dx, dy, &
- k_start, k_end, kds,kde, ims,ime, jms,jme, kms,kme )
- !******************************************************************************
- ! Purpose: Interpolate pressure at a specified x (ic), y (jc), and height (h).
- ! The input pressure column pbbc+ppbc (base and perturbn) must already
- ! be horizontally interpolated to the x, y position. The subroutine
- ! get_height_column is called here to horizontally interpolated the
- ! 3D height field z to get a height column at (iob, job).
- !******************************************************************************
- IMPLICIT NONE
- REAL, INTENT(IN) :: h ! height value (m)
- INTEGER, INTENT(IN) :: k_start, k_end ! loop bounds
- INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
- INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
- REAL, INTENT(IN) :: pbbc(kds:kde) ! column base pressure (cb)
- REAL, INTENT(IN) :: ppbc(kds:kde) ! column pressure perturbation (cb)
- REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) above sl on half-levels
- INTEGER, INTENT(IN) :: ic ! i-coord of desired p
- INTEGER, INTENT(IN) :: jc ! j-coord of desired p
- REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
- REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
- ! Local variables
- INTEGER :: k ! loop counter
- INTEGER :: klo ! lower k bound
- REAL :: zlo ! lower z bound for h
- REAL :: zhi ! upper z bound for h
- REAL :: p ! interpolated pressure value
- REAL :: ln_p ! log p
- REAL :: ln_plo ! log plo
- REAL :: ln_phi ! log phi
- REAL :: z_at_p( kms:kme ) ! height at p levels
- ! Get interpolated z column on pressure (half-) levels at (ic,jc)
- call get_height_column(z, ic, jc, dx, dy, z_at_p, &
- k_start, k_end, kds,kde, &
- ims,ime, jms,jme, kms,kme )
- ! Now we have pbbc, ppbc, z_at_p, so compute p at h. First, find
- ! bounding layers klo and khi so that z_at_p(klo) <= h <= z_at_p(khi)
- ZLEVS: do k = k_start+1, k_end
- klo = k-1
- if(h .le. z_at_p(k)) then
- EXIT ZLEVS
- endif
- enddo ZLEVS
- zlo = z_at_p(klo)
- zhi = z_at_p(klo+1)
- ! Interpolate natural log of pressure
- ln_plo = log( pbbc(klo+1) + ppbc(klo+1) )
- ln_phi = log( pbbc(klo) + ppbc(klo) )
- if(h.le.zlo) then
- ln_p = ln_phi ! set to k=1 pressure
- else if (h.ge.zhi) then
- ln_p = ln_plo ! set to k=k_end pressure
- else
- ln_p = ln_plo + (ln_phi-ln_plo)*((zhi-h)/(zhi-zlo))
- endif
- ! Return pressure
- p = exp(ln_p)
- ht_to_p = p
- RETURN
- END FUNCTION ht_to_p
- SUBROUTINE get_height_column( z, ic, jc, dx, dy, z_at_p, &
- k_start, k_end, kds,kde, &
- ims,ime, jms,jme, kms,kme )
- !*************************************************************************
- ! Purpose: Compute column height at ic, jc location on p levels
- !*************************************************************************
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
- INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
- INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
- REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) on half-levels
- INTEGER, INTENT(IN) :: ic ! i-coord of desired p
- INTEGER, INTENT(IN) :: jc ! j-coord of desired p
- REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
- REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
- REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
- ! Local variables
- INTEGER :: k ! loop counter
- do k = kds, kde
- z_at_p(k) = &
- (1.-DY)*( (1.-DX)*z(IC,K,JC) + &
- DX *z(IC+1,K,JC) ) + &
- DY* ( (1.-DX)*z(IC,K,JC+1) + &
- DX *z(IC+1,K,JC+1) )
- enddo
- END SUBROUTINE get_height_column
- SUBROUTINE get_base_state_height_column( p_top, p00, t00, a, g, r_d, &
- znu, z_at_p, k_start, k_end, kds,kde, kms,kme )
- !********************************************************
- ! Purpose: Compute base-state column height on p levels.
- ! See, e.g., eqns 1.3-1.5 in MM5 User's Guide.
- ! Height is a function of pressure plus the constants:
- ! p00 - sea level pressure (Pa)
- ! t00 - sea level temperature (K)
- ! a - base state lapse rate (k/m)
- ! r_d - gas constant (J/Kg/K) (Joule=1kg/m**2/s**2)
- ! g - gravitational constant (m/s**2)
- !********************************************************
- IMPLICIT NONE
- REAL, INTENT(IN) :: p_top ! pressure at top of model
- REAL, INTENT(IN) :: p00 ! base state pressure
- REAL, INTENT(IN) :: t00 ! base state temperature
- REAL, INTENT(IN) :: a ! base state lapse rate
- REAL, INTENT(IN) :: g ! gravity constant
- REAL, INTENT(IN) :: r_d ! gas constant
- INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
- INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
- INTEGER, INTENT(IN) :: kms,kme ! vertical memory dim.
- REAL, INTENT(IN) :: znu( kms:kme ) ! eta values on half (mass) levels
- REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
- ! Local variables
- integer :: k ! loop counter
- real :: ps0 ! base state pressure at surface
- real :: pb(kms:kme) ! pressure on half eta levels
- real :: logterm ! temporary for Z calculation
- real :: ginv ! 1/g
-
- ginv = 1/g
- ! Compute base state pressure on half eta levels.
- do k = k_start, k_end
- pb(k) = znu(k)*(p00 - p_top) + p_top
- enddo
- ! Use hydrostatic relation to compute height at pressure levels.
- do k = k_start, k_end
- logterm = log(pb(k)/p00)
- z_at_p(k) = .5*r_d*a*ginv*logterm*logterm - r_d*t00*ginv*logterm
- enddo
- END SUBROUTINE get_base_state_height_column
- REAL FUNCTION get_timewt(xtime,dtmin,twindo,scalef,obtime)
- !*************************************************************************
- ! Purpose: Compute the temporal weight factor for an observation
- !*************************************************************************
- IMPLICIT NONE
- REAL, INTENT(IN) :: xtime ! model time (minutes)
- REAL, INTENT(IN) :: dtmin ! model timestep (minutes)
- REAL, INTENT(IN) :: twindo ! half window (hours)
- REAL, INTENT(IN) :: scalef ! window scale factor
- REAL, INTENT(IN) :: obtime ! observation time (hours)
- ! Local variables
- real :: fdtim ! reference time (minutes)
- real :: tw1 ! half of twindo, scaled, in minutes
- real :: tw2 ! twindo, scaled, in minutes
- real :: tconst ! reciprical of tw1
- real :: ttim ! obtime in minutes
- real :: dift ! | fdtim-ttim |
- real :: timewt ! returned weight
- ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
- FDTIM=XTIME-DTMIN
- ! TWINDO IS IN MINUTES:
- TW1=TWINDO/2.*60.*scalef
- TW2=TWINDO*60.*scalef
- TCONST=1./TW1
- TIMEWT=0.0
- TTIM=obtime*60.
- !***********TTIM=TARGET TIME IN MINUTES
- DIFT=ABS(FDTIM-TTIM)
- IF(DIFT.LE.TW1)TIMEWT=1.0
- IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
- IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
- IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
- ENDIF
- get_timewt = timewt
- END FUNCTION get_timewt
- SUBROUTINE print_vif_var(var, vif, nfullmin, nrampmin )
- !********************************************************
- ! Purpose: Print a description of the vertical influence
- ! function for a given variable.
- !********************************************************
- IMPLICIT NONE
- character(len=4), intent(in) :: var ! Variable (wind, temp, mois)
- real, intent(in) :: vif(6) ! Vertical influence function
- real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
- real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
- ! Local variables
- character(len=200) :: msg1, msg2
- character(len=8) :: regime
- real :: nfullr1, nrampr1
- real :: nfullr2, nrampr2
- real :: nfullr4, nrampr4
- nfullr1 = vif(1)
- nrampr1 = vif(2)
- nfullr2 = vif(3)
- nrampr2 = vif(4)
- nfullr4 = vif(5)
- nrampr4 = vif(6)
- if(var.eq.'wind') then
- write(msg1,fmt='(a)') ' For winds:'
- elseif (var.eq.'temp') then
- write(msg1,fmt='(a)') ' For temperature:'
- elseif (var.eq.'mois') then
- write(msg1,fmt='(a)') ' For moisture:'
- else
- write(msg1,fmt='(a,a4)') 'Unknown variable type: ',var
- call wrf_message(msg1)
- call wrf_error_fatal ( 'print_vif_var: module_fddaobs_rtfdda STOP' )
- endif
-
- call wrf_message(msg1)
- ! For this variable, print a description of the vif for each regime
- call print_vif_regime(1, nfullr1, nrampr1, nfullmin, nrampmin)
- call print_vif_regime(2, nfullr2, nrampr2, nfullmin, nrampmin)
- call print_vif_regime(4, nfullr4, nrampr4, nfullmin, nrampmin)
- END SUBROUTINE print_vif_var
- SUBROUTINE print_vif_regime(reg, nfullr, nrampr, nfullmin, nrampmin )
- !********************************************************
- ! Purpose: Print a description of the vertical influence
- ! function for a given regime.
- !********************************************************
- IMPLICIT NONE
- integer, intent(in) :: reg ! Regime number (1, 2, 4)
- real, intent(in) :: nfullr ! Full nudge range for regime
- real, intent(in) :: nrampr ! Rampdown range for regime
- real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
- real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
- ! Local variables
- character(len=200) :: msg1, msg2
- character(len=8) :: regime
- if(reg.eq.1) then
- write(regime,fmt='(a)') 'Regime 1'
- elseif (reg.eq.2) then
- write(regime,fmt='(a)') 'Regime 2'
- elseif (reg.eq.4) then
- write(regime,fmt='(a)') 'Regime 4'
- else
- write(msg1,fmt='(a,i3)') 'Unknown regime number: ',reg
- call wrf_message(msg1)
- call wrf_error_fatal ( 'print_vif_regime: module_fddaobs_rtfdda STOP' )
- endif
- !Set msg1 for description of full weighting range
- if(nfullr.lt.0) then
- if(nfullr.eq.-5000) then
- write(msg1,fmt='(2x,a8,a)') regime, ': Full weighting to the PBL top'
- elseif (nfullr.lt.-5000) then
- write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(-5000.-nfullr), &
- ' m above the PBL top'
- else
- write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(nfullr+5000.), &
- ' m below the PBL top'
- endif
- else
- write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting through ', &
- int(max(nfullr,nfullmin)),' m'
- endif
- !Set msg2 for description of rampdown range
- if(nrampr.lt.0) then
- if(nrampr.eq.-5000) then
- write(msg2,fmt='(a)') ' and a vertical rampdown up to the PBL top.'
- elseif (nrampr.lt.-5000) then
- write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(-5000.-nrampr), &
- ' m above the PBL top.'
- else
- write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(nrampr+5000.), &
- ' m below the PBL top.'
- endif
- else
- write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown in the next ', &
- int(max(nrampr,nrampmin)),' m.'
- endif
- call wrf_message(TRIM(msg1)//msg2)
- END SUBROUTINE print_vif_regime
- #endif
- END MODULE module_fddaobs_rtfdda