PageRenderTime 74ms CodeModel.GetById 18ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_fddaobs_rtfdda.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2978 lines | 1716 code | 378 blank | 884 comment | 137 complexity | 951b09af3ef210354f9fd397a7ef6ab5 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. MODULE module_fddaobs_rtfdda
    4. 

  4. 

  5. ! This obs-nudging FDDA module (RTFDDA) is developed by the 
    6. 

  6. ! NCAR/RAL/NSAP (National Security Application Programs), under the 
    7. 

  7. ! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is 
    8. 

  8. ! acknowledged for releasing this capability for WRF community 
    9. 

  9. ! research applications.
    10. 

  10. !
    11. 

  11. ! The NCAR/RAL RTFDDA module was adapted, and significantly modified 
    12. 

  12. ! from the obs-nudging module in the standard MM5V3.1 which was originally 
    13. 

  13. ! developed by PSU (Stauffer and Seaman, 1994). 
    14. 

  14. ! 
    15. 

  15. ! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module 
    16. 

  16. ! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
    17. 

  17. ! Nov. 2006
    18. 

  18. ! 
    19. 

  19. ! References:
    20. 

  20. !   
    21. 

  21. !   Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
    22. 

  22. !     implementation of obs-nudging-based FDDA into WRF for supporting 
    23. 

  23. !     ATEC test operations. 2005 WRF user workshop. Paper 10.7.
    24. 

  24. !
    25. 

  25. !   Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update 
    26. 

  26. !     on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE 
    27. 

  27. !     and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7. 
    28. 

  28. 

  29. !   
    30. 

  30. !   Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data 
    31. 

  31. !     assimilation. J. Appl. Meteor., 33, 416-434.
    32. 

  32. !
    33. 

  33. !   http://www.rap.ucar.edu/projects/armyrange/references.html
    34. 

  34. !
    35. 

  35. ! Modification History:
    36. 

  36. !   03212007 Modified fddaobs_init to compute Lambert cone factor. -Al Bourgeois
    37. 

  37. 

  38. CONTAINS
    39. 

  39. 

  40. !------------------------------------------------------------------------------
    41. 

  41.   SUBROUTINE fddaobs_init(nudge_opt, maxdom, inest, parid,             &
    42. 

  42.                           idynin, dtramp, fdaend, restart,             &
    43. 

  43.                           twindo_cg, twindo, itimestep,                &
    44. 

  44.                           no_pbl_nudge_uv,                             &
    45. 

  45.                           no_pbl_nudge_t,                              &
    46. 

  46.                           no_pbl_nudge_q,                              &
    47. 

  47.                           sfc_scheme_horiz, sfc_scheme_vert,           &
    48. 

  48.                           maxsnd_gap,                                  &
    49. 

  49.                           sfcfact, sfcfacr, dpsmx,                     &
    50. 

  50.                           nudge_wind, nudge_temp, nudge_mois,          &
    51. 

  51.                           nudgezfullr1_uv, nudgezrampr1_uv,            &
    52. 

  52.                           nudgezfullr2_uv, nudgezrampr2_uv,            &
    53. 

  53.                           nudgezfullr4_uv, nudgezrampr4_uv,            &
    54. 

  54.                           nudgezfullr1_t,  nudgezrampr1_t,             &
    55. 

  55.                           nudgezfullr2_t,  nudgezrampr2_t,             &
    56. 

  56.                           nudgezfullr4_t,  nudgezrampr4_t,             &
    57. 

  57.                           nudgezfullr1_q,  nudgezrampr1_q,             &
    58. 

  58.                           nudgezfullr2_q,  nudgezrampr2_q,             &
    59. 

  59.                           nudgezfullr4_q,  nudgezrampr4_q,             &
    60. 

  60.                           nudgezfullmin,   nudgezrampmin, nudgezmax,   &
    61. 

  61.                           xlat, xlong,                                 &
    62. 

  62.                           start_year, start_month, start_day,          &
    63. 

  63.                           start_hour, start_minute, start_second,      &
    64. 

  64.                           p00, t00, tlp,                               &
    65. 

  65.                           znu, p_top,                                  &
    66. 

  66. #if ( EM_CORE == 1 )
    67. 

  67.                           fdob,                                        &
    68. 

  68. #endif
    69. 

  69.                           iprt,                                        &
    70. 

  70.                           ids,ide, jds,jde, kds,kde,                   &
    71. 

  71.                           ims,ime, jms,jme, kms,kme,                   &
    72. 

  72.                           its,ite, jts,jte, kts,kte)     
    73. 

  73. !-----------------------------------------------------------------------
    74. 

  74. !  This routine does initialization for real time fdda obs-nudging.
    75. 

  75. !
    76. 

  76. !-----------------------------------------------------------------------
    77. 

  77.   USE module_model_constants, ONLY : g, r_d
    78. 

  78.   USE module_domain
    79. 

  79.   USE module_dm, ONLY : wrf_dm_min_real
    80. 

  80. !-----------------------------------------------------------------------
    81. 

  81.   IMPLICIT NONE
    82. 

  82. !-----------------------------------------------------------------------
    83. 

  83. 

  84. !=======================================================================
    85. 

  85. ! Definitions
    86. 

  86. !-----------
    87. 

  87.   INTEGER, intent(in)  :: maxdom
    88. 

  88.   INTEGER, intent(in)  :: nudge_opt(maxdom)
    89. 

  89.   INTEGER, intent(in)  :: ids,ide, jds,jde, kds,kde,                 &
    90. 

  90.                           ims,ime, jms,jme, kms,kme,                 &
    91. 

  91.                           its,ite, jts,jte, kts,kte
    92. 

  92.   INTEGER, intent(in)  :: inest
    93. 

  93.   INTEGER, intent(in)  :: parid(maxdom)
    94. 

  94.   INTEGER, intent(in)  :: idynin         ! flag for dynamic initialization
    95. 

  95.   REAL,    intent(in)  :: dtramp         ! time period for idynin ramping (min)
    96. 

  96.   REAL,    intent(in)  :: fdaend(maxdom) ! nudging end time for domain (min)
    97. 

  97.   LOGICAL, intent(in)  :: restart
    98. 

  98.   REAL, intent(in)     :: twindo_cg      ! time window on coarse grid
    99. 

  99.   REAL, intent(in)     :: twindo
    100. 

  100.   INTEGER, intent(in)  :: itimestep
    101. 

  101.   INTEGER , INTENT(IN) :: no_pbl_nudge_uv(maxdom)  ! flags for no wind nudging in pbl 
    102. 

  102.   INTEGER , INTENT(IN) :: no_pbl_nudge_t(maxdom)   ! flags for no temperature nudging in pbl
    103. 

  103.   INTEGER , INTENT(IN) :: no_pbl_nudge_q(maxdom)   ! flags for no moisture nudging in pbl
    104. 

  104.   INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
    105. 

  105.   INTEGER , INTENT(IN) :: sfc_scheme_vert  ! vertical spreading scheme for surf obs (orig or regime vif)
    106. 

  106.   REAL    , INTENT(IN) :: maxsnd_gap      ! max allowed pressure gap in soundings, for interp (centibars)
    107. 

  107.   REAL, intent(in)     :: sfcfact      ! scale factor applied to time window for surface obs  
    108. 

  108.   REAL, intent(in)     :: sfcfacr      ! scale fac applied to horiz rad of infl for sfc obs
    109. 

  109.   REAL, intent(in)     :: dpsmx        ! max press change allowed within hor rad of infl
    110. 

  110.   INTEGER , INTENT(IN) :: nudge_wind(maxdom)       ! wind-nudging flag
    111. 

  111.   INTEGER , INTENT(IN) :: nudge_temp(maxdom)       ! temperature-nudging flag
    112. 

  112.   INTEGER , INTENT(IN) :: nudge_mois(maxdom)       ! moisture-nudging flag
    113. 

  113.   REAL, INTENT(IN)     :: nudgezfullr1_uv  ! vert infl fcn, regime=1 full-wt   hght, winds
    114. 

  114.   REAL, INTENT(IN)     :: nudgezrampr1_uv  ! vert infl fcn, regime=1 ramp down hght, winds
    115. 

  115.   REAL, INTENT(IN)     :: nudgezfullr2_uv  ! vert infl fcn, regime=2 full-wt   hght, winds
    116. 

  116.   REAL, INTENT(IN)     :: nudgezrampr2_uv  ! vert infl fcn, regime=2 ramp down hght, winds
    117. 

  117.   REAL, INTENT(IN)     :: nudgezfullr4_uv  ! vert infl fcn, regime=4 full-wt   hght, winds
    118. 

  118.   REAL, INTENT(IN)     :: nudgezrampr4_uv  ! vert infl fcn, regime=4 ramp down hght, winds
    119. 

  119.   REAL, INTENT(IN)     :: nudgezfullr1_t   ! vert infl fcn, regime=1 full-wt   hght, temp
    120. 

  120.   REAL, INTENT(IN)     :: nudgezrampr1_t   ! vert infl fcn, regime=1 ramp down hght, temp
    121. 

  121.   REAL, INTENT(IN)     :: nudgezfullr2_t   ! vert infl fcn, regime=2 full-wt   hght, temp
    122. 

  122.   REAL, INTENT(IN)     :: nudgezrampr2_t   ! vert infl fcn, regime=2 ramp down hght, temp
    123. 

  123.   REAL, INTENT(IN)     :: nudgezfullr4_t   ! vert infl fcn, regime=4 full-wt   hght, temp
    124. 

  124.   REAL, INTENT(IN)     :: nudgezrampr4_t   ! vert infl fcn, regime=4 ramp down hght, temp
    125. 

  125.   REAL, INTENT(IN)     :: nudgezfullr1_q   ! vert infl fcn, regime=1 full-wt   hght, mois
    126. 

  126.   REAL, INTENT(IN)     :: nudgezrampr1_q   ! vert infl fcn, regime=1 ramp down hght, mois
    127. 

  127.   REAL, INTENT(IN)     :: nudgezfullr2_q   ! vert infl fcn, regime=2 full-wt   hght, mois
    128. 

  128.   REAL, INTENT(IN)     :: nudgezrampr2_q   ! vert infl fcn, regime=2 ramp down hght, mois
    129. 

  129.   REAL, INTENT(IN)     :: nudgezfullr4_q   ! vert infl fcn, regime=4 full-wt   hght, mois
    130. 

  130.   REAL, INTENT(IN)     :: nudgezrampr4_q   ! vert infl fcn, regime=4 ramp down hght, mois
    131. 

  131.   REAL, INTENT(IN)     :: nudgezfullmin    ! min dpth thru which vert infl fcn remains 1.0 (m)
    132. 

  132.   REAL, INTENT(IN)     :: nudgezrampmin    ! min dpth thru which vif decreases 1.0 to 0.0 (m)
    133. 

  133.   REAL, INTENT(IN)     :: nudgezmax        ! max dpth in which vif is nonzero (m)
    134. 

  134.   REAL, INTENT(IN)     :: xlat ( ims:ime, jms:jme )        ! latitudes on mass-point grid
    135. 

  135.   REAL, INTENT(IN)     :: xlong( ims:ime, jms:jme )        ! longitudes on mass-point grid
    136. 

  136.   INTEGER, intent(in)  :: start_year   ! Model start year
    137. 

  137.   INTEGER, intent(in)  :: start_month  ! Model start month
    138. 

  138.   INTEGER, intent(in)  :: start_day    ! Model start day
    139. 

  139.   INTEGER, intent(in)  :: start_hour   ! Model start hour
    140. 

  140.   INTEGER, intent(in)  :: start_minute ! Model start minute
    141. 

  141.   INTEGER, intent(in)  :: start_second ! Model start second
    142. 

  142.   REAL, INTENT(IN)     :: p00          ! base state pressure
    143. 

  143.   REAL, INTENT(IN)     :: t00          ! base state temperature
    144. 

  144.   REAL, INTENT(IN)     :: tlp          ! base state lapse rate
    145. 

  145.   REAL, INTENT(IN)     :: p_top        ! pressure at top of model
    146. 

  146.   REAL, INTENT(IN)     :: znu( kms:kme )      ! eta values on half (mass) levels
    147. 

  147. #if ( EM_CORE == 1 ) 
    148. 

  148.   TYPE(fdob_type), intent(inout)  :: fdob
    149. 

  149. #endif
    150. 

  150.   LOGICAL, intent(in)  :: iprt         ! Flag enabling printing warning messages
    151. 

  151. 

  152. ! Local variables
    153. 

  153.   logical            :: nudge_flag      ! nudging flag for this nest 
    154. 

  154.   integer            :: ktau            ! current timestep
    155. 

  155.   integer            :: nest            ! loop counter
    156. 

  156.   integer            :: idom            ! domain id
    157. 

  157.   integer            :: parent          ! parent domain
    158. 

  158.   real               :: conv            ! 180/pi
    159. 

  159.   real               :: tl1             ! Lambert standard parallel 1
    160. 

  160.   real               :: tl2             ! Lambert standard parallel 2
    161. 

  161.   real               :: xn1
    162. 

  162.   real               :: known_lat       ! Latitude of domain point (i,j)=(1,1)
    163. 

  163.   real               :: known_lon       ! Longitude of domain point (i,j)=(1,1) 
    164. 

  164.   character(len=200) :: msg             ! Argument to wrf_message
    165. 

  165.   real               :: z_at_p( kms:kme )  ! height at p levels
    166. 

  166.   integer            :: i,j,k           ! loop counters
    167. 

  167. 

  168. #if ( EM_CORE == 1 ) 
    169. 

  169. 

  170. ! Check to see if the nudging flag has been set. If not,
    171. 

  171. ! simply RETURN.
    172. 

  172.   nudge_flag = (nudge_opt(inest) .eq. 1)
    173. 

  173.   if (.not. nudge_flag) return
    174. 

  174. 

  175.   call wrf_message("")
    176. 

  176.   write(msg,fmt='(a,i2)') ' OBSERVATION NUDGING IS ACTIVATED FOR MESH ',inest
    177. 

  177.   call wrf_message(msg)
    178. 

  178. 

  179.   ktau  = itimestep
    180. 

  180.   if(restart) then
    181. 

  181.     fdob%ktaur = ktau
    182. 

  182.   else
    183. 

  183.     fdob%ktaur = 0 
    184. 

  184.   endif
    185. 

  185. 

  186. ! Create character string containing model starting-date
    187. 

  187.   CALL date_string(start_year, start_month, start_day, start_hour,        &
    188. 

  188.                    start_minute, start_second, fdob%sdate)
    189. 

  189. 

  190. ! Set flag for nudging on pressure (not sigma) surfaces.
    191. 

  191.   fdob%iwtsig = 0
    192. 

  192. 

  193. !**************************************************************************
    194. 

  194. ! *** Initialize datend for dynamic initialization (ajb added 08132008) ***
    195. 

  195. !**************************************************************************
    196. 

  196. ! Set ending nudging date (used with dynamic ramp-down) to zero.
    197. 

  197.   fdob%datend = 0.
    198. 

  198.   fdob%ieodi = 0
    199. 

  199. 

  200. ! Check for dynamic initialization flag
    201. 

  201.   if(idynin.eq.1)then
    202. 

  202. !    Set datend to time in minutes after which data are assumed to have ended.
    203. 

  203.      if(dtramp.gt.0.)then
    204. 

  204.         fdob%datend = fdaend(inest) - dtramp
    205. 

  205.      else
    206. 

  206.         fdob%datend = fdaend(inest)
    207. 

  207.      endif
    208. 

  208.      if(iprt) then
    209. 

  209.         call wrf_message("")
    210. 

  210.         write(msg,fmt='(a,i3,a)')                                              &
    211. 

  211.           ' *** DYNAMIC-INITIALIZATION OPTION FOR INEST = ', inest, ' ***'
    212. 

  212.         call wrf_message(msg)
    213. 

  213.         write(msg,*) ' FDAEND,DATEND,DTRAMP: ',fdaend(inest),fdob%datend,dtramp
    214. 

  214.         call wrf_message(msg)
    215. 

  215.         call wrf_message("")
    216. 

  216.      endif
    217. 

  217.   endif
    218. 

  218. ! *** end dynamic initialization section ***
    219. 

  219. !**************************************************************************
    220. 

  220. 

  221. ! Store flags for surface obs spreading schemes
    222. 

  222.   if(sfc_scheme_horiz.eq.1) then
    223. 

  223.      call wrf_message('MM5 scheme selected for horizontal spreading of surface obs')
    224. 

  224.   elseif (sfc_scheme_horiz.eq.0) then
    225. 

  225.      call wrf_message('WRF scheme selected for horizontal spreading of surface obs')
    226. 

  226.   else
    227. 

  227.      write(msg,fmt='(a,i3)') 'Unknown h-spreading scheme for surface obs: ',sfc_scheme_horiz
    228. 

  228.      call wrf_message(msg)
    229. 

  229.      call wrf_message("Valid selections: 0=WRF scheme, 1=Original MM5 scheme")
    230. 

  230.      call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
    231. 

  231.   endif
    232. 

  232. 

  233.   if(sfc_scheme_vert.eq.1) then
    234. 

  234.      call wrf_message('Original simple scheme selected for vertical spreading of surface obs')
    235. 

  235.   elseif (sfc_scheme_vert.eq.0) then
    236. 

  236.      call wrf_message("Regime-based VIF scheme selected for vertical spreading of surface obs")
    237. 

  237.   else
    238. 

  238.      write(msg,fmt='(a,i3)') 'Unknown v-spreading scheme for surface obs: ',sfc_scheme_vert
    239. 

  239.      call wrf_message(msg)
    240. 

  240.      call wrf_message("Valid selections: 0=Regime-based VIF scheme, 1=Original simple scheme")
    241. 

  241.      call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
    242. 

  242.   endif
    243. 

  243.   fdob%sfc_scheme_horiz = sfc_scheme_horiz 
    244. 

  244.   fdob%sfc_scheme_vert  = sfc_scheme_vert
    245. 

  245. 

  246. 

  247. ! Store temporal and spatial scales
    248. 

  248.   fdob%sfcfact = sfcfact
    249. 

  249.   fdob%sfcfacr = sfcfacr
    250. 

  250. 

  251. ! Set time window.
    252. 

  252.   fdob%window = twindo
    253. 

  253.   call wrf_message("")
    254. 

  254.   write(msg,fmt='(a,i3)') '*** TIME WINDOW SETTINGS FOR NEST ',inest
    255. 

  255.   call wrf_message(msg)
    256. 

  256.   write(msg,fmt='(a,f6.3,2(a,f5.3))') '    TWINDO (hrs) = ',twindo,      &
    257. 

  257.             '  SFCFACT = ',sfcfact,'  SFCFACR = ',sfcfacr
    258. 

  258.   call wrf_message(msg)
    259. 

  259.   call wrf_message("")
    260. 

  260. 

  261.   if(inest.eq.1) then
    262. 

  262.     if(twindo .eq. 0.) then
    263. 

  263.       if(iprt) then
    264. 

  264.         call wrf_message("")
    265. 

  265.         write(msg,*) '*** WARNING: TWINDO=0 on the coarse domain.'
    266. 

  266.         call wrf_message(msg)
    267. 

  267.         write(msg,*) '*** Did you forget to set twindo in the fdda namelist?'
    268. 

  268.         call wrf_message(msg)
    269. 

  269.         call wrf_message("")
    270. 

  270.       endif
    271. 

  271.     endif
    272. 

  272.   else        ! nest
    273. 

  273.     if(twindo .eq. 0.) then
    274. 

  274.       fdob%window = twindo_cg
    275. 

  275.       if(iprt) then
    276. 

  276.         call wrf_message("")
    277. 

  277.         write(msg,fmt='(a,i2)') 'WARNING: TWINDO=0. for nest ',inest
    278. 

  278.         call wrf_message(msg)
    279. 

  279.         write(msg,fmt='(a,f12.5,a)') 'Default to coarse-grid value of ', twindo_cg,' hrs'
    280. 

  280.         call wrf_message(msg)
    281. 

  281.         call wrf_message("")
    282. 

  282.       endif
    283. 

  283.     endif
    284. 

  284.   endif
    285. 

  285. 

  286. ! Initialize flags.
    287. 

  287. 

  288.   fdob%domain_tot=0
    289. 

  289.   do nest=1,maxdom
    290. 

  290.     fdob%domain_tot = fdob%domain_tot + nudge_opt(nest)
    291. 

  291.   end do
    292. 

  292. 

  293. ! Calculate and store dcon from dpsmx
    294. 

  294.   if(dpsmx.gt.0.) then
    295. 

  295.        fdob%dpsmx = dpsmx
    296. 

  296.        fdob%dcon = 1.0/fdob%dpsmx
    297. 

  297.   else 
    298. 

  298.        call wrf_error_fatal('fddaobs_init: Namelist variable dpsmx must be greater than zero!')
    299. 

  299.   endif
    300. 

  300. 

  301. ! Calculate and store base-state heights at half (mass) levels.
    302. 

  302.   CALL get_base_state_height_column( p_top, p00, t00, tlp, g, r_d, znu,   &
    303. 

  303.                                      fdob%base_state,  kts, kte, kds,kde, kms,kme )
    304. 

  304. 

  305. ! Initialize flags for nudging within PBL.
    306. 

  306.   fdob%nudge_uv_pbl  = .true.
    307. 

  307.   fdob%nudge_t_pbl   = .true.
    308. 

  308.   fdob%nudge_q_pbl   = .true.
    309. 

  309.   if(no_pbl_nudge_uv(inest) .eq. 1) fdob%nudge_uv_pbl  = .false.
    310. 

  310.   if(no_pbl_nudge_t(inest) .eq. 1)  fdob%nudge_t_pbl   = .false.
    311. 

  311.   if(no_pbl_nudge_q(inest) .eq. 1)  fdob%nudge_q_pbl   = .false.
    312. 

  312. 

  313.   if(no_pbl_nudge_uv(inest) .eq. 1) then
    314. 

  314.     fdob%nudge_uv_pbl  = .false.
    315. 

  315.     write(msg,*) '   --> Obs nudging for U/V is turned off in PBL'
    316. 

  316.     call wrf_message(msg)
    317. 

  317.   endif
    318. 

  318.   if(no_pbl_nudge_t(inest) .eq. 1)  then
    319. 

  319.     fdob%nudge_t_pbl   = .false.
    320. 

  320.     write(msg,*) '   --> Obs nudging for T is turned off in PBL'
    321. 

  321.     call wrf_message(msg)
    322. 

  322.   endif
    323. 

  323.   if(no_pbl_nudge_q(inest) .eq. 1)  then
    324. 

  324.     fdob%nudge_q_pbl   = .false.
    325. 

  325.     write(msg,*) '   --> Obs nudging for Q is turned off in PBL'
    326. 

  326.     call wrf_message(msg)
    327. 

  327.   endif
    328. 

  328. 

  329. ! Store max allowed pressure gap for interpolating between soundings
    330. 

  330.   fdob%max_sndng_gap = maxsnd_gap
    331. 

  331.   write(msg,fmt='(a,f6.1)')  &
    332. 

  332.   '*** MAX PRESSURE GAP (cb) for interpolation between soundings = ',maxsnd_gap
    333. 

  333.   call wrf_message(msg)
    334. 

  334.   call wrf_message("")
    335. 

  335. 

  336. ! Initialize vertical influence fcn for LML obs 
    337. 

  337.   if(sfc_scheme_vert.eq.0) then
    338. 

  338.     fdob%vif_uv(1) = nudgezfullr1_uv
    339. 

  339.     fdob%vif_uv(2) = nudgezrampr1_uv
    340. 

  340.     fdob%vif_uv(3) = nudgezfullr2_uv
    341. 

  341.     fdob%vif_uv(4) = nudgezrampr2_uv
    342. 

  342.     fdob%vif_uv(5) = nudgezfullr4_uv
    343. 

  343.     fdob%vif_uv(6) = nudgezrampr4_uv
    344. 

  344.     fdob%vif_t (1) = nudgezfullr1_t
    345. 

  345.     fdob%vif_t (2) = nudgezrampr1_t
    346. 

  346.     fdob%vif_t (3) = nudgezfullr2_t
    347. 

  347.     fdob%vif_t (4) = nudgezrampr2_t
    348. 

  348.     fdob%vif_t (5) = nudgezfullr4_t
    349. 

  349.     fdob%vif_t (6) = nudgezrampr4_t
    350. 

  350.     fdob%vif_q (1) = nudgezfullr1_q
    351. 

  351.     fdob%vif_q (2) = nudgezrampr1_q
    352. 

  352.     fdob%vif_q (3) = nudgezfullr2_q
    353. 

  353.     fdob%vif_q (4) = nudgezrampr2_q
    354. 

  354.     fdob%vif_q (5) = nudgezfullr4_q
    355. 

  355.     fdob%vif_q (6) = nudgezrampr4_q
    356. 

  356. 

  357. ! Sanity checks
    358. 

  358.     if(nudgezmax.le.0.) then
    359. 

  359.       write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezmax MUST BE GREATER THAN ZERO.'
    360. 

  360.       call wrf_message(msg)
    361. 

  361.       write(msg,*) 'THE NAMELIST VALUE IS',nudgezmax
    362. 

  362.       call wrf_message(msg)
    363. 

  363.       call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgemax value' )
    364. 

  364.     endif
    365. 

  365.     if(nudgezfullmin.lt.0.) then
    366. 

  366.       write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezfullmin MUST BE NONNEGATIVE.'
    367. 

  367.       call wrf_message(msg)
    368. 

  368.       write(msg,*) 'THE NAMELIST VALUE IS',nudgezfullmin
    369. 

  369.       call wrf_message(msg)
    370. 

  370.       call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgefullmin value' )
    371. 

  371.     endif
    372. 

  372.     if(nudgezrampmin.lt.0.) then
    373. 

  373.       write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezrampmin MUST BE NONNEGATIVE.'
    374. 

  374.       call wrf_message(msg)
    375. 

  375.       write(msg,*) 'THE NAMELIST VALUE IS',nudgezrampmin
    376. 

  376.       call wrf_message(msg)
    377. 

  377.       call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgerampmin value' )
    378. 

  378.     endif
    379. 

  379.     if(nudgezmax.lt.nudgezfullmin+nudgezrampmin) then
    380. 

  380.       write(msg,*) 'STOP! INCONSISTENT OBS NAMELIST INPUTS.'
    381. 

  381.       call wrf_message(msg)
    382. 

  382.       write(msg,fmt='(3(a,f12.3))') 'obs_nudgezmax = ',nudgezmax,                &
    383. 

  383.                               ' obs_nudgezfullmin = ',nudgezfullmin,       &
    384. 

  384.                               ' obs_nudgezrampmin = ',nudgezrampmin
    385. 

  385.       call wrf_message(msg)
    386. 

  386.       write(msg,*) 'REQUIRE NUDGEZMAX >= NUDGEZFULLMIN + NUDGEZRAMPMIN'
    387. 

  387.       call wrf_message(msg)
    388. 

  388.       call wrf_error_fatal ( 'fddaobs_init: STOP on inconsistent namelist values' )
    389. 

  389.     endif
    390. 

  390.  
    391. 

  391.     fdob%vif_fullmin = nudgezfullmin
    392. 

  392.     fdob%vif_rampmin = nudgezrampmin
    393. 

  393.     fdob%vif_max     = nudgezmax
    394. 

  394.  
    395. 

  395. ! Check to make sure that if nudgzfullmin > 0, then it must be at least as large as the
    396. 

  396. ! first model half-level will be anywhere in the domain at any time within the simulation.
    397. 

  397. ! We use 1.1 times the base-state value fdob%base_state(1) for this purpose. 
    398. 

  398. 

  399.     if(nudgezfullmin.gt.0.0) then
    400. 

  400.         if(nudgezfullmin .lt. 1.1*fdob%base_state(1)) then
    401. 

  401.            fdob%vif_fullmin = 1.1*fdob%base_state(1)
    402. 

  402.         endif
    403. 

  403.     endif
    404. 

  404. 

  405. ! Print vertical weight info only if wind, temperature, or moisture nudging is requested.
    406. 

  406.     if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1)                             &
    407. 

  407.                                                .or. (nudge_mois(inest).eq.1) ) then
    408. 

  408.     call wrf_message("")
    409. 

  409.     write(msg,fmt='(a,i2,a)') ' *** SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest,' :'
    410. 

  410.     call wrf_message(msg)
    411. 

  411. 

  412.     call wrf_message("")
    413. 

  413.     write(msg,fmt='(a,i5,a)') '  NUDGEZMAX: The maximum height at which nudging will be'//     &
    414. 

  414.                         ' applied from surface obs is ', nint(nudgezmax),' m AGL.'
    415. 

  415.     call wrf_message(msg)
    416. 

  416.     call wrf_message("")
    417. 

  417.     write(msg,fmt='(a,i3,a)') '  NUDGEZFULLMIN: The minimum height of full nudging weight'//   &
    418. 

  418.                           ' for surface obs is ', nint(fdob%vif_fullmin),' m.'
    419. 

  419.     call wrf_message(msg)
    420. 

  420.     if(nudgezfullmin.lt.fdob%vif_fullmin) then
    421. 

  421.         write(msg,fmt='(a,i3,a)') '  ***WARNING***: NUDGEZFULLMIN has been increased from'//   &
    422. 

  422.                               ' the user-input value of ',nint(nudgezfullmin),' m.'
    423. 

  423.         call wrf_message(msg)
    424. 

  424.         write(msg,fmt='(a,i3,a)') '  to ensure that at least the bottom model level is'//      &
    425. 

  425.                               ' included in full nudging.'
    426. 

  426.         call wrf_message(msg)
    427. 

  427.     endif
    428. 

  428.     call wrf_message("")
    429. 

  429.     write(msg,fmt='(a,i3,a)') '  NUDGEZRAMPMIN: The minimum height to ramp from full to no'//  &
    430. 

  430.                           ' nudging for surface obs is ', nint(nudgezrampmin),' m.'
    431. 

  431.     call wrf_message(msg)
    432. 

  432.     call wrf_message("")
    433. 

  433.     endif   !endif either wind, temperature, or moisture nudging is requested
    434. 

  434. 

  435. ! Print vif settings
    436. 

  436.     if(nudge_wind(inest) .eq. 1) then
    437. 

  437.       call print_vif_var('wind', fdob%vif_uv, nudgezfullmin, nudgezrampmin)
    438. 

  438.       call wrf_message("")
    439. 

  439.     endif
    440. 

  440.     if(nudge_temp(inest) .eq. 1) then
    441. 

  441.       call print_vif_var('temp', fdob%vif_t,  nudgezfullmin, nudgezrampmin)
    442. 

  442.       call wrf_message("")
    443. 

  443.     endif
    444. 

  444.     if(nudge_mois(inest) .eq. 1) then
    445. 

  445.       call print_vif_var('mois', fdob%vif_q,  nudgezfullmin, nudgezrampmin)
    446. 

  446.       call wrf_message("")
    447. 

  447.     endif
    448. 

  448. 

  449.     if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1)                             &
    450. 

  450.                                                  .or. (nudge_mois(inest).eq.1) ) then
    451. 

  451.     write(msg,fmt='(a,i2)') ' *** END SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest
    452. 

  452.     call wrf_message(msg)
    453. 

  453.     call wrf_message("")
    454. 

  454.     endif
    455. 

  455.   endif   !endif(sfc_scheme_vert.EQ.0)
    456. 

  456. 

  457. ! Set parameters.
    458. 

  458.   fdob%pfree = 50.0
    459. 

  459.   fdob%rinfmn = 1.0
    460. 

  460.   fdob%rinfmx = 2.0
    461. 

  461. 

  462. ! Get known lat and lon and store these on all processors in fdob structure, for
    463. 

  463. ! later use in projection x-formation to map (lat,lon) to (x,y) for each obs.
    464. 

  464.       IF (its .eq. 1 .AND. jts .eq. 1) then
    465. 

  465.          known_lat = xlat(1,1)
    466. 

  466.          known_lon = xlong(1,1)
    467. 

  467.       ELSE
    468. 

  468.          known_lat = 9999.
    469. 

  469.          known_lon = 9999.
    470. 

  470.       END IF
    471. 

  471.       fdob%known_lat = wrf_dm_min_real(known_lat)
    472. 

  472.       fdob%known_lon = wrf_dm_min_real(known_lon)
    473. 

  473. 

  474. ! Calculate the nest levels, levidn. Note that each nest
    475. 

  475. ! must know the nest levels levidn(maxdom) of each domain.
    476. 

  476.   do nest=1,maxdom
    477. 

  477. 

  478. ! Initialize nest level for each domain.
    479. 

  479.     if (nest .eq. 1) then
    480. 

  480.       fdob%levidn(nest) = 0  ! Mother domain has nest level 0
    481. 

  481.     else
    482. 

  482.       fdob%levidn(nest) = 1  ! All other domains start with 1
    483. 

  483.     endif
    484. 

  484.     idom = nest
    485. 

  485. 100 parent = parid(idom)      ! Go up the parent tree
    486. 

  486. 

  487.       if (parent .gt. 1) then   ! If not up to mother domain
    488. 

  488.         fdob%levidn(nest) = fdob%levidn(nest) + 1
    489. 

  489.         idom = parent
    490. 

  490.         goto 100
    491. 

  491.       endif
    492. 

  492.   enddo
    493. 

  493. 

  494. 

  495. ! fdob%LML_OBS_HT1_LEV = kte 
    496. 

  496. ! HT1: do k = kte, kts, -1
    497. 

  497. !        if( LML_HT1 .gt. z_at_p(k) ) then
    498. 

  498. !          fdob%LML_OBS_HT1_LEV = k
    499. 

  499. !          EXIT HT1 
    500. 

  500. !        endif
    501. 

  501. !      enddo  HT1
    502. 

  502. 

  503. ! fdob%LML_OBS_HT2_LEV = kte 
    504. 

  504. ! HT2: do k = kte, kts, -1
    505. 

  505. !        if( LML_HT2 .gt. z_at_p(k) ) then
    506. 

  506. !          fdob%LML_OBS_HT2_LEV = k
    507. 

  507. !          EXIT HT2 
    508. 

  508. !        endif
    509. 

  509. !      enddo HT2 
    510. 

  510.   RETURN
    511. 

  511. #endif
    512. 

  512.   END SUBROUTINE fddaobs_init
    513. 

  513. 

  514. #if ( EM_CORE == 1 )
    515. 

  515. !-----------------------------------------------------------------------
    516. 

  516. SUBROUTINE errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp,  &
    517. 

  517.                  z,                                             &
    518. 

  518.                  uratx, vratx, tratx, kpbl,                     &
    519. 

  519.                  nndgv, nerrf, niobf, maxdom,                   &
    520. 

  520.                  levidn, parid, nstat, nstaw,                   &
    521. 

  521.                  iswind, istemp, ismois, ispstr,                &
    522. 

  522.                  timeob, rio, rjo, rko,                         &
    523. 

  523.                  varobs, errf, ktau, xtime,                     &
    524. 

  524.                  iratio, npfi,                                  &
    525. 

  525.                  prt_max, prt_freq, iprt,                       &
    526. 

  526.                  obs_prt, stnid_prt, lat_prt, lon_prt,          &
    527. 

  527.                  mlat_prt, mlon_prt,                            & 
    528. 

  528.                  ids,ide, jds,jde, kds,kde,                     &
    529. 

  529.                  ims,ime, jms,jme, kms,kme,                     &
    530. 

  530.                  its,ite, jts,jte, kts,kte  )
    531. 

  531. 

  532. !-----------------------------------------------------------------------
    533. 

  533. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    534. 

  534.   USE module_dm, ONLY : get_full_obs_vector, wrf_dm_sum_real
    535. 

  535. #else
    536. 

  536.   USE module_dm, ONLY :                      wrf_dm_sum_real
    537. 

  537. #endif
    538. 

  538.   USE module_model_constants, ONLY : rcp
    539. 

  539. 

  540. !-----------------------------------------------------------------------
    541. 

  541.   IMPLICIT NONE
    542. 

  542. !-----------------------------------------------------------------------
    543. 

  543. !
    544. 

  544. ! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
    545. 

  545. !     OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
    546. 

  546. !     POINTS.  THE OBSERVED VALUES CLOSEST TO THE CURRENT
    547. 

  547. !     FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
    548. 

  548. !     IN4DOB AND STORED IN ARRAY VAROBS.  THE DIFFERENCES
    549. 

  549. !     CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
    550. 

  550. !     ERRF FOR THE NSTA OBSERVATION LOCATIONS.  MISSING
    551. 

  551. !     OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
    552. 

  552. !
    553. 

  553. !     HISTORY: Original author: MM5 version???
    554. 

  554. !              02/04/2004 - Creation of WRF version.           Al Bourgeois
    555. 

  555. !              08/28/2006 - Conversion from F77 to F90         Al Bourgeois
    556. 

  556. !------------------------------------------------------------------------------
    557. 

  557. 

  558. ! THE STORAGE ORDER IN ERRF IS AS FOLLOWS:
    559. 

  559. !        IVAR                VARIABLE TYPE(TAU-1)
    560. 

  560. !        ----                --------------------
    561. 

  561. !         1                    U error at obs loc
    562. 

  562. !         2                    V error at obs loc
    563. 

  563. !         3                    T error at obs loc
    564. 

  564. !         4                    Q error at obs loc
    565. 

  565. !         5                    Vertical layer index for PBL top at IOB, JOB
    566. 

  566. !         6                    Model surface press at obs loc (T-points)
    567. 

  567. !         7                    Model surface press at obs loc (U-points)
    568. 

  568. !         8                    Model surface press at obs loc (V-points)
    569. 

  569. !         9                    RKO at U-points
    570. 

  570. 

  571. !-----------------------------------------------------------------------
    572. 

  572. !
    573. 

  573. !     Description of input arguments.
    574. 

  574. !
    575. 

  575. !-----------------------------------------------------------------------
    576. 

  576. 

  577.   INTEGER, INTENT(IN)  :: inest                  ! Domain index.
    578. 

  578.   INTEGER, INTENT(IN)  :: nndgv                  ! Number of nudge variables.
    579. 

  579.   INTEGER, INTENT(IN)  :: nerrf                  ! Number of error fields.
    580. 

  580.   INTEGER, INTENT(IN)  :: niobf                  ! Number of observations.
    581. 

  581.   INTEGER, INTENT(IN)  :: maxdom                 ! Maximum number of domains.
    582. 

  582.   INTEGER, INTENT(IN)  :: levidn(maxdom)         ! Level of nest.
    583. 

  583.   INTEGER, INTENT(IN)  :: parid(maxdom)          ! Id of parent grid.
    584. 

  584.   INTEGER, INTENT(IN)  :: ktau                   ! Model time step index
    585. 

  585.   REAL, INTENT(IN)     :: xtime                  ! Forecast time in minutes
    586. 

  586.   INTEGER, INTENT(IN)  :: iratio                 ! Nest to parent gridsize ratio.
    587. 

  587.   INTEGER, INTENT(IN)  :: npfi                   ! Coarse-grid diagnostics freq.
    588. 

  588.   INTEGER, INTENT(IN)  :: prt_max                ! Max number of obs to print.
    589. 

  589.   INTEGER, INTENT(IN)  :: prt_freq               ! Frequency of obs to print.
    590. 

  590.   LOGICAL, INTENT(IN)  :: iprt                   ! Print flag
    591. 

  591.   INTEGER, INTENT(IN)  :: obs_prt(prt_max)       ! Print obs indices
    592. 

  592.   INTEGER, INTENT(IN)  :: stnid_prt(40,prt_max)  ! Print obs station ids
    593. 

  593.   REAL, INTENT(IN)     :: lat_prt(prt_max)       ! Print obs latitude
    594. 

  594.   REAL, INTENT(IN)     :: lon_prt(prt_max)       ! Print obs longitude
    595. 

  595.   REAL, INTENT(IN)     :: mlat_prt(prt_max)      ! Print model lat at obs loc
    596. 

  596.   REAL, INTENT(IN)     :: mlon_prt(prt_max)      ! Print model lon at obs loc
    597. 

  597.   INTEGER, INTENT(IN)  :: nstat                  ! # stations held for use
    598. 

  598.   INTEGER, INTENT(IN)  :: nstaw                  ! # stations in current window
    599. 

  599.   INTEGER, intent(in)  :: iswind
    600. 

  600.   INTEGER, intent(in)  :: istemp
    601. 

  601.   INTEGER, intent(in)  :: ismois
    602. 

  602.   INTEGER, intent(in)  :: ispstr
    603. 

  603.   INTEGER, INTENT(IN)  :: ids,ide, jds,jde, kds,kde  ! domain dims.
    604. 

  604.   INTEGER, INTENT(IN)  :: ims,ime, jms,jme, kms,kme  ! memory dims.
    605. 

  605.   INTEGER, INTENT(IN)  :: its,ite, jts,jte, kts,kte  ! tile   dims.
    606. 

  606. 

  607.   REAL,   INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
    608. 

  608.   REAL,   INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
    609. 

  609.   REAL,   INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
    610. 

  610.   REAL,   INTENT(IN) :: t0
    611. 

  611.   REAL,   INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
    612. 

  612.   REAL,   INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
    613. 

  613.   REAL,   INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
    614. 

  614.   REAL,   INTENT(IN)  :: rovcp
    615. 

  615.   REAL,    INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! Ht above sl on half-levels
    616. 

  616.   REAL,   INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
    617. 

  617.   REAL,   INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
    618. 

  618.   REAL,   INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
    619. 

  619.   INTEGER,INTENT(IN) :: kpbl( ims:ime, jms:jme )  ! Vertical layer index for PBL top
    620. 

  620.   REAL,   INTENT(IN) :: timeob(niobf)             ! Obs time (hrs)
    621. 

  621.   REAL,   INTENT(IN) :: rio(niobf)                ! Obs west-east coordinate (non-stag grid).
    622. 

  622.   REAL,   INTENT(IN) :: rjo(niobf)                ! Obs south-north coordinate (non-stag grid).
    623. 

  623.   REAL,   INTENT(INOUT) :: rko(niobf)             ! Obs bottom-top coordinate
    624. 

  624.   REAL,   INTENT(INOUT) :: varobs(nndgv, niobf)
    625. 

  625.   REAL,   INTENT(INOUT) :: errf(nerrf, niobf)
    626. 

  626. 

  627. ! Local variables
    628. 

  628.   INTEGER :: iobmg(niobf)   ! Obs i-coord on mass grid
    629. 

  629.   INTEGER :: jobmg(niobf)   ! Obs j-coord on mass grid
    630. 

  630.   INTEGER :: ia(niobf)
    631. 

  631.   INTEGER :: ib(niobf)
    632. 

  632.   INTEGER :: ic(niobf)
    633. 

  633.   REAL :: pbbo(kds:kde)    ! column base pressure (cb) at obs loc.
    634. 

  634.   REAL :: ppbo(kds:kde)    ! column pressure perturbation (cb) at obs loc.
    635. 

  635. 

  636.   REAL :: ra(niobf)
    637. 

  637.   REAL :: rb(niobf)
    638. 

  638.   REAL :: rc(niobf)
    639. 

  639.   REAL :: dxobmg(niobf)     ! Interp. fraction (x dir) referenced to mass-grid
    640. 

  640.   REAL :: dyobmg(niobf)     ! Interp. fraction (y dir) referenced to mass-grid
    641. 

  641.   INTEGER MM(MAXDOM)
    642. 

  642.   INTEGER NNL
    643. 

  643.   real :: uratio( ims:ime, jms:jme )   ! U to U10 ratio on momentum points.
    644. 

  644.   real :: vratio( ims:ime, jms:jme )   ! V to V10 ratio on momentum points.
    645. 

  645.   real :: pug1,pug2,pvg1,pvg2
    646. 

  646.   character(len=200) :: msg            ! Argument to wrf_message
    647. 

  647. 

  648. ! Define staggers for U, V, and T grids, referenced from non-staggered grid.
    649. 

  649.   real, parameter :: gridx_t = 0.5     ! Mass-point x stagger
    650. 

  650.   real, parameter :: gridy_t = 0.5     ! Mass-point y stagger
    651. 

  651.   real, parameter :: gridx_u = 0.0     ! U-point x stagger
    652. 

  652.   real, parameter :: gridy_u = 0.5     ! U-point y stagger
    653. 

  653.   real, parameter :: gridx_v = 0.5     ! V-point x stagger
    654. 

  654.   real, parameter :: gridy_v = 0.0     ! V-point y stagger
    655. 

  655. 

  656.   real :: dummy = 99999.
    657. 

  657. 

  658.   real :: pbhi, pphi
    659. 

  659.   real :: obs_pottemp                  ! Potential temperature at observation
    660. 

  660. 

  661. !***  DECLARATIONS FOR IMPLICIT NONE
    662. 

  662.   integer nsta,ivar,n,ityp
    663. 

  663.   integer iob,job,kob,iob_ms,job_ms
    664. 

  664.   integer k,kbot,nml,nlb,nle
    665. 

  665.   integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
    666. 

  666.   integer i_start,i_end,j_start,j_end    ! loop ranges for uratio,vratio calc.
    667. 

  667.   integer k_start,k_end
    668. 

  668.   integer ips                            ! For printing obs information
    669. 

  669.   integer pnx                            ! obs index for printout
    670. 

  670. 

  671.   real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
    672. 

  672.   real pob
    673. 

  673.   real hob
    674. 

  674.   real uratiob,vratiob,tratiob,tratxob,fnpf
    675. 

  675. 

  676. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    677. 

  677.   LOGICAL MP_LOCAL_DUMMASK(NIOBF)  ! Mask for work to be done on this processor
    678. 

  678.   LOGICAL MP_LOCAL_UOBMASK(NIOBF)  ! Dot-point mask
    679. 

  679.   LOGICAL MP_LOCAL_VOBMASK(NIOBF)  ! Dot-point mask
    680. 

  680.   LOGICAL MP_LOCAL_COBMASK(NIOBF)  ! Cross-point mask
    681. 

  681. #endif
    682. 

  682. 

  683. ! LOGICAL, EXTERNAL :: TILE_MASK
    684. 

  684. 

  685.   NSTA=NSTAT
    686. 

  686. 

  687. ! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
    688. 

  688. ! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
    689. 

  689. ! TO THE FINE MESH INDEX EQUIVALENTS
    690. 

  690. 

  691. ! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
    692. 

  692. 

  693.   if (iprt) then
    694. 

  694.     write(msg,fmt='(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
    695. 

  695.             KTAU,' AND INEST = ',INEST,':  NSTA = ',NSTAW,' ++++++'
    696. 

  696.     call wrf_message(msg)
    697. 

  697.   endif
    698. 

  698. 

  699.   ERRF = 0.0    ! Zero out errf array
    700. 

  700. 

  701. ! Set up loop bounds for this grid's boundary conditions
    702. 

  702.   i_start = max( its-1,ids )
    703. 

  703.   i_end   = min( ite+1,ide-1 )
    704. 

  704.   j_start = max( jts-1,jds )
    705. 

  705.   j_end   = min( jte+1,jde-1 )
    706. 

  706.   k_start = kts
    707. 

  707.   k_end = min( kte, kde-1 )
    708. 

  708. 

  709.   DO ITYP=1,3   ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP 
    710. 

  710. 

  711. !   Set grid stagger
    712. 

  712.     IF(ITYP.EQ.1) THEN        ! U-POINT CASE
    713. 

  713.        GRIDX = gridx_u
    714. 

  714.        GRIDY = gridy_u
    715. 

  715.     ELSE IF(ITYP.EQ.2) THEN   ! V-POINT CASE
    716. 

  716.        GRIDX = gridx_v
    717. 

  717.        GRIDY = gridy_v
    718. 

  718.     ELSE                      ! MASS-POINT CASE
    719. 

  719.        GRIDX = gridx_t
    720. 

  720.        GRIDY = gridy_t
    721. 

  721.     ENDIF
    722. 

  722. 

  723. !   Compute URATIO and VRATIO fields on momentum (u,v) points.
    724. 

  724.     IF(ityp.eq.1)THEN
    725. 

  725.       call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
    726. 

  726.     ELSE IF (ityp.eq.2) THEN
    727. 

  727.       call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
    728. 

  728.     ENDIF
    729. 

  729. 

  730.     IF(INEST.EQ.1) THEN       ! COARSE MESH CASE...
    731. 

  731.       DO N=1,NSTA
    732. 

  732.         RA(N)=RIO(N)-GRIDX
    733. 

  733.         RB(N)=RJO(N)-GRIDY
    734. 

  734.         IA(N)=RA(N)
    735. 

  735.         IB(N)=RB(N)
    736. 

  736.         IOB=MAX0(1,IA(N))
    737. 

  737.         IOB=MIN0(IOB,ide-1)
    738. 

  738.         JOB=MAX0(1,IB(N))
    739. 

  739.         JOB=MIN0(JOB,jde-1)
    740. 

  740.         DXOB=RA(N)-FLOAT(IA(N))
    741. 

  741.         DYOB=RB(N)-FLOAT(IB(N))
    742. 

  742. 

  743. !       Save mass-point arrays for computing rko for all var types
    744. 

  744.         if(ityp.eq.1) then
    745. 

  745.           iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
    746. 

  746.           jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
    747. 

  747.           dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
    748. 

  748.           dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
    749. 

  749.         endif
    750. 

  750.         iob_ms = iobmg(n)
    751. 

  751.         job_ms = jobmg(n)
    752. 

  752.         dxob_ms = dxobmg(n)
    753. 

  753.         dyob_ms = dyobmg(n)
    754. 

  754. 

  755. 

  756. ! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION
    757. 

  757. ! This is tricky: Initialize pob to zero in all procs. Only one proc actually
    758. 

  758. ! calculates pob. If this is an obs to be converted from height-to-pressure, then
    759. 

  759. ! by definition, varobs(5,n) will initially have the missing value -888888. After
    760. 

  760. ! the pob calculation, pob will be zero in all procs except the one that calculated
    761. 

  761. ! it, and so pob is dm_summed over all procs and stored into varobs(5,n). So on
    762. 

  762. ! subsequent passes, the dm_sum will not occur because varobs(5,n) will not have a
    763. 

  763. ! missing value. If this is not an obs-in-height,  
    764. 

  764. 

  765.         pob = 0.0
    766. 

  766. 

  767. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    768. 

  768. !       Set mask for obs to be handled by this processor
    769. 

  769.         MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
    770. 

  770.   
    771. 

  771.         IF ( MP_LOCAL_DUMMASK(N) ) THEN
    772. 

  772. #endif
    773. 

  773. 

  774. !         Interpolate pressure to obs location column and convert from Pa to cb.
    775. 

  775. 

  776.           do k = kds, kde
    777. 

  777.             pbbo(k) = .001*(                                            &
    778. 

  778.                (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) +     &
    779. 

  779.                                   DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
    780. 

  780.                    DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) +   &
    781. 

  781.                                   DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )  
    782. 

  782.             ppbo(k) = .001*(                                            &
    783. 

  783.                (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) +        &
    784. 

  784.                                   DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) +    &
    785. 

  785.                    DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) +      &
    786. 

  786.                                   DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )  
    787. 

  787.           enddo
    788. 

  788. 

  789. !ajb      20040119: Note based on bugfix for dot/cross points split across processors,
    790. 

  790. !ajb                which was actually a serial code fix: The ityp=2 (v-points) and
    791. 

  791. !ajb                ityp=3 (mass-points) cases should not use the ityp=1 (u-points)
    792. 

  792. !ajb                case rko! This is necessary for bit-for-bit reproducability
    793. 

  793. !ajb                with the parallel run.   (coarse mesh)
    794. 

  794. 

  795.           if(abs(rko(n)+99).lt.1.)then
    796. 

  796.             pob = varobs(5,n)
    797. 

  797. 

  798.             if(pob .eq.-888888.) then
    799. 

  799.                hob = varobs(6,n)
    800. 

  800.                if(hob .gt. -800000. ) then
    801. 

  801.                  pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms,          &
    802. 

  802.                                 dxob_ms, dyob_ms, k_start, k_end, kds,kde,   &
    803. 

  803.                                 ims,ime, jms,jme, kms,kme )
    804. 

  804.                endif
    805. 

  805.             endif
    806. 

  806. 

  807.             if(pob .gt.-800000.)then
    808. 

  808.               do k=k_end-1,1,-1
    809. 

  809.                 kbot = k
    810. 

  810.                 if(pob .le. pbbo(k)+ppbo(k)) then
    811. 

  811.                   goto 199
    812. 

  812.                 endif
    813. 

  813.               enddo
    814. 

  814.  199          continue
    815. 

  815. 

  816.               pphi = ppbo(kbot+1)
    817. 

  817.               pbhi = pbbo(kbot+1)
    818. 

  818. 

  819.               rko(n) = real(kbot+1)-                                    &
    820. 

  820.                  ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
    821. 

  821. 

  822.               rko(n)=max(rko(n),1.0)
    823. 

  823.             endif
    824. 

  824.           endif
    825. 

  825. 

  826. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    827. 

  827.         ENDIF       !end IF( MP_LOCAL_DUMMASK(N) )                                 !ajb
    828. 

  828. #endif
    829. 

  829. 

  830. ! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
    831. 

  831.         if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
    832. 

  832.            varobs(5,n) = wrf_dm_sum_real ( pob )
    833. 

  833.         endif
    834. 

  834. 

  835.         RC(N)=RKO(N)
    836. 

  836. 

  837.       ENDDO      ! END COARSE MESH LOOP OVER NSTA 
    838. 

  838. 

  839.     ELSE       ! FINE MESH CASE
    840. 

  840. 

  841. !**********************************************************************
    842. 

  842. !ajb_07012008: Conversion of obs coordinates to the fine mesh here
    843. 

  843. !ajb   is no longer necessary, since implementation of the WRF map pro-
    844. 

  844. !ajb   jections (to each nest directly) is implemented in sub in4dob. 
    845. 

  845. !**********************************************************************
    846. 

  846. !ajb
    847. 

  847. !ajb  GET (I,J,K) OF OBSERVATIONS ON FINE MESH VALUES.
    848. 

  848.       DO N=1,NSTA
    849. 

  849.           RA(N)=RIO(N)-GRIDX           ! ajb added 07012008
    850. 

  850.           RB(N)=RJO(N)-GRIDY           ! ajb added 07012008
    851. 

  851.           IA(N)=RA(N)
    852. 

  852.           IB(N)=RB(N)
    853. 

  853.           IOB=MAX0(1,IA(N))
    854. 

  854.           IOB=MIN0(IOB,ide-1)
    855. 

  855.           JOB=MAX0(1,IB(N))
    856. 

  856.           JOB=MIN0(JOB,jde-1)
    857. 

  857.           DXOB=RA(N)-FLOAT(IA(N))
    858. 

  858.           DYOB=RB(N)-FLOAT(IB(N))
    859. 

  859. 

  860. !         Save mass-point arrays for computing rko for all var types
    861. 

  861.           if(ityp.eq.1) then
    862. 

  862.             iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
    863. 

  863.             jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
    864. 

  864.             dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
    865. 

  865.             dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
    866. 

  866.           endif 
    867. 

  867.           iob_ms = iobmg(n)
    868. 

  868.           job_ms = jobmg(n)
    869. 

  869.           dxob_ms = dxobmg(n)
    870. 

  870.           dyob_ms = dyobmg(n)
    871. 

  871. 

  872. ! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION (see COARSE MESH comments)
    873. 

  873.         pob = 0.0
    874. 

  874. 

  875. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    876. 

  876. !       Set mask for obs to be handled by this processor
    877. 

  877.         MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
    878. 

  878. 

  879.         IF ( MP_LOCAL_DUMMASK(N) ) THEN
    880. 

  880. #endif
    881. 

  881. 

  882. !         Interpolate pressure to obs location column and convert from Pa to cb.
    883. 

  883. 

  884.           do k = kds, kde
    885. 

  885.             pbbo(k) = .001*(                                            &
    886. 

  886.                (1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) +     &
    887. 

  887.                                   DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
    888. 

  888.                    DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) +   &
    889. 

  889.                                   DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
    890. 

  890.             ppbo(k) = .001*(                                            &
    891. 

  891.                (1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) +        &
    892. 

  892.                                   DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) +    &
    893. 

  893.                    DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) +      &
    894. 

  894.                                   DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
    895. 

  895.           enddo
    896. 

  896. 

  897. !ajb      20040119: Note based on bugfix for dot/cross points split across processors,
    898. 

  898. !ajb                which was actually a serial code fix: The ityp=2 (v-points) and
    899. 

  899. !ajb                itype=3 (mass-points) cases should not use the ityp=1 (u-points)
    900. 

  900. !ajb                case) rko! This is necessary for bit-for-bit reproducability
    901. 

  901. !ajb                with parallel run.   (fine mesh)
    902. 

  902. 

  903.           if(abs(rko(n)+99).lt.1.)then
    904. 

  904.             pob = varobs(5,n)
    905. 

  905. 

  906.             if(pob .eq.-888888.) then
    907. 

  907.                hob = varobs(6,n)
    908. 

  908.                if(hob .gt. -800000. ) then
    909. 

  909.                  pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms,          &
    910. 

  910.                                 dxob_ms, dyob_ms, k_start, k_end, kds,kde,   &
    911. 

  911.                                 ims,ime, jms,jme, kms,kme )
    912. 

  912.                endif
    913. 

  913.             endif
    914. 

  914. 

  915.             if(pob .gt.-800000.)then
    916. 

  916.               do k=k_end-1,1,-1
    917. 

  917.                 kbot = k
    918. 

  918.                 if(pob .le. pbbo(k)+ppbo(k)) then
    919. 

  919.                   goto 198
    920. 

  920.                 endif
    921. 

  921.               enddo
    922. 

  922.  198          continue
    923. 

  923. 

  924.               pphi = ppbo(kbot+1)
    925. 

  925.               pbhi = pbbo(kbot+1)
    926. 

  926. 

  927.               rko(n) = real(kbot+1)-                                    &
    928. 

  928.                  ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
    929. 

  929.               rko(n)=max(rko(n),1.0)
    930. 

  930.             endif
    931. 

  931.           endif
    932. 

  932. 

  933. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    934. 

  934.         ENDIF       !end IF( MP_LOCAL_DUMMASK(N) )                                 !ajb
    935. 

  935. #endif
    936. 

  936. 

  937. ! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
    938. 

  938.         if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
    939. 

  939.            varobs(5,n) = wrf_dm_sum_real ( pob )
    940. 

  940.         endif
    941. 

  941. 

  942.         RC(N)=RKO(N)
    943. 

  943. 

  944.       ENDDO      ! END FINE MESH LOOP OVER NSTA
    945. 

  945.     
    946. 

  946.     ENDIF      ! end if(inest.eq.1)
    947. 

  947. 

  948. !   INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
    949. 

  949. !   AND THE LOCAL FORECAST VALUES TO ZERO.  FOR THE FINE MESH
    950. 

  950. !   ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
    951. 

  951. !   OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
    952. 

  952. 

  953. !   SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
    954. 

  954. !   CURRENT DOMAIN
    955. 

  955.     IF(ITYP.EQ.1) THEN
    956. 

  956.       NLB=1
    957. 

  957.       NLE=1
    958. 

  958.     ELSE IF(ITYP.EQ.2) THEN
    959. 

  959.       NLB=2
    960. 

  960.       NLE=2
    961. 

  961.     ELSE
    962. 

  962.       NLB=3
    963. 

  963.       NLE=5
    964. 

  964.     ENDIF
    965. 

  965.     DO IVAR=NLB,NLE
    966. 

  966.       DO N=1,NSTA
    967. 

  967.         IF((RA(N)-1.).LT.0)THEN
    968. 

  968.            ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
    969. 

  969.         ENDIF
    970. 

  970.         IF((RB(N)-1.).LT.0)THEN
    971. 

  971.            ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
    972. 

  972.         ENDIF
    973. 

  973.         IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
    974. 

  974.            ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
    975. 

  975.         ENDIF
    976. 

  976.         IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
    977. 

  977.            ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
    978. 

  978.         ENDIF
    979. 

  979.         if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
    980. 

  980.       ENDDO
    981. 

  981.     ENDDO
    982. 

  982. 

  983. !   NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
    984. 

  984. !   GRID POINT TOWARD THE LOWER LEFT
    985. 

  985.     DO N=1,NSTA
    986. 

  986.         IA(N)=RA(N)
    987. 

  987.         IB(N)=RB(N)
    988. 

  988.         IC(N)=RC(N)
    989. 

  989.     ENDDO
    990. 

  990.     DO N=1,NSTA
    991. 

  991.         RA(N)=RA(N)-FLOAT(IA(N))
    992. 

  992.         RB(N)=RB(N)-FLOAT(IB(N))
    993. 

  993.         RC(N)=RC(N)-FLOAT(IC(N))
    994. 

  994.     ENDDO
    995. 

  995. !   PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
    996. 

  996. !   TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
    997. 

  997. !   POINTS FOR U, V, T, AND Q.
    998. 

  998. 

  999. !   Compute local masks for dot and cross points.
    1000. 

  1000.     if(ityp.eq.1) then
    1001. 

  1001.       DO N=1,NSTA
    1002. 

  1002.         IOB=MAX0(1,IA(N))
    1003. 

  1003.         IOB=MIN0(IOB,ide-1)
    1004. 

  1004.         JOB=MAX0(1,IB(N))
    1005. 

  1005.         JOB=MIN0(JOB,jde-1)
    1006. 

  1006. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1007. 

  1007. !       Set mask for U-momemtum points to be handled by this processor
    1008. 

  1008.         MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
    1009. 

  1009. #endif
    1010. 

  1010.       ENDDO
    1011. 

  1011.     endif
    1012. 

  1012.     if(ityp.eq.2) then
    1013. 

  1013.       DO N=1,NSTA
    1014. 

  1014.         IOB=MAX0(1,IA(N))
    1015. 

  1015.         IOB=MIN0(IOB,ide-1)
    1016. 

  1016.         JOB=MAX0(1,IB(N))
    1017. 

  1017.         JOB=MIN0(JOB,jde-1)
    1018. 

  1018. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1019. 

  1019. !       Set mask for V-momentum points to be handled by this processor
    1020. 

  1020.         MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
    1021. 

  1021. #endif
    1022. 

  1022.       ENDDO
    1023. 

  1023.     endif
    1024. 

  1024.     if(ityp.eq.3) then
    1025. 

  1025.       DO N=1,NSTA
    1026. 

  1026.         IOB=MAX0(1,IA(N))
    1027. 

  1027.         IOB=MIN0(IOB,ide-1)
    1028. 

  1028.         JOB=MAX0(1,IB(N))
    1029. 

  1029.         JOB=MIN0(JOB,jde-1)
    1030. 

  1030. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1031. 

  1031. !       Set mask for cross (mass) points to be handled by this processor
    1032. 

  1032.         MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
    1033. 

  1033. #endif
    1034. 

  1034.       ENDDO
    1035. 

  1035.     endif
    1036. 

  1036. 

  1037. !**********************************************************
    1038. 

  1038. !   PROCESS U VARIABLE (IVAR=1)
    1039. 

  1039. !**********************************************************
    1040. 

  1040.     IF(ITYP.EQ.1) THEN
    1041. 

  1041. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1042. 

  1042.       DO N=1,NSTA
    1043. 

  1043.         IF(MP_LOCAL_UOBMASK(N)) THEN
    1044. 

  1044.           ERRF(9,N)=rko(n)       !RKO is needed by neighboring processors     !ajb
    1045. 

  1045.         ENDIF
    1046. 

  1046.       ENDDO
    1047. 

  1047. #endif
    1048. 

  1048.       IF(ISWIND.EQ.1) THEN
    1049. 

  1049.         DO N=1,NSTA
    1050. 

  1050.           IOB=MAX0(2,IA(N))
    1051. 

  1051.           IOB=MIN0(IOB,ide-1)
    1052. 

  1052.           IOBM=MAX0(1,IOB-1)
    1053. 

  1053.           IOBP=MIN0(ide-1,IOB+1)
    1054. 

  1054.           JOB=MAX0(2,IB(N))
    1055. 

  1055.           JOB=MIN0(JOB,jde-1)
    1056. 

  1056.           JOBM=MAX0(1,JOB-1)
    1057. 

  1057.           JOBP=MIN0(jde-1,JOB+1)
    1058. 

  1058.           KOB=MIN0(K_END,IC(N))
    1059. 

  1059. 

  1060. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1061. 

  1061.           IF(MP_LOCAL_UOBMASK(N))THEN     ! Do if obs on this processor
    1062. 

  1062. #endif
    1063. 

  1063.             KOBP=MIN0(KOB+1,k_end)
    1064. 

  1064.             DXOB=RA(N)
    1065. 

  1065.             DYOB=RB(N)
    1066. 

  1066.             DZOB=RC(N)
    1067. 

  1067. 

  1068. !           Compute surface pressure values at surrounding U and V points
    1069. 

  1069.             PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
    1070. 

  1070.             PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
    1071. 

  1071. 

  1072. !           This is to correct obs to model sigma level using reverse similarity theory
    1073. 

  1073.             if(rko(n).eq.1.0)then
    1074. 

  1074.               uratiob=((1.-DYOB)*((1.-DXOB)*uratio(IOB,JOB)+     &
    1075. 

  1075.                     DXOB*uratio(IOBP,JOB)                        &
    1076. 

  1076.                   )+DYOB*((1.-DXOB)*uratio(IOB,JOBP)+            &
    1077. 

  1077.                   DXOB*uratio(IOBP,JOBP)))
    1078. 

  1078.             else
    1079. 

  1079.               uratiob=1.
    1080. 

  1080.             endif
    1081. 

  1081. !YLIU       Some PBL scheme do not define the vratio/uratio
    1082. 

  1082.             if(abs(uratiob).lt.1.0e-3) then
    1083. 

  1083.               uratiob=1.
    1084. 

  1084.             endif
    1085. 

  1085. 

  1086. !           INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
    1087. 

  1087. !           WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
    1088. 

  1088. !           OUTSIDE THE CURRENT DOMAIN
    1089. 

  1089. 

  1090. !           U COMPONENT WIND ERROR
    1091. 

  1091.             ERRF(1,N)=ERRF(1,N)+uratiob*VAROBS(1,N)-((1.-DZOB)*        &
    1092. 

  1092.                       ((1.-DyOB)*((1.-                                 &
    1093. 

  1093.                       DxOB)*UB(IOB,KOB,JOB)+DxOB*UB(IOB+1,KOB,JOB)     &
    1094. 

  1094.                       )+DyOB*((1.-DxOB)*UB(IOB,KOB,JOB+1)+DxOB*        &
    1095. 

  1095.                       UB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
    1096. 

  1096.                       *UB(IOB,KOBP,JOB)+DxOB*UB(IOB+1,KOBP,JOB))+      &
    1097. 

  1097.                       DyOB*((1.-DxOB)*UB(IOB,KOBP,JOB+1)+DxOB*         &
    1098. 

  1098.                       UB(IOB+1,KOBP,JOB+1))))
    1099. 

  1099. 

  1100. !      if(n.le.10) then
    1101. 

  1101. !        write(6,*)
    1102. 

  1102. !        write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob,   &
    1103. 

  1103. !                                     ' N = ',n,' inest = ',inest
    1104. 

  1104. !        write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
    1105. 

  1105. !        write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
    1106. 

  1106. !        write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
    1107. 

  1107. !        write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
    1108. 

  1108. !        write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
    1109. 

  1109. !        write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
    1110. 

  1110. !        write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
    1111. 

  1111. !        write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
    1112. 

  1112. !        write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
    1113. 

  1113. !        write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
    1114. 

  1114. !        write(6,*) 'uratiob = ',uratiob
    1115. 

  1115. !        write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
    1116. 

  1116. !        write(6,*) 'ERRF(1,N) = ',errf(1,n)
    1117. 

  1117. !        write(6,*)
    1118. 

  1118. !      endif
    1119. 

  1119. 

  1120. 

  1121. !           Store model surface pressure (not the error!) at U point.
    1122. 

  1122.             ERRF(7,N)=.001*( (1.-DXOB)*PUG1 + DXOB*PUG2 )
    1123. 

  1123.   
    1124. 

  1124. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1125. 

  1125.           ENDIF       ! end IF( MP_LOCAL_UOBMASK(N) )
    1126. 

  1126. #endif
    1127. 

  1127.         ENDDO    ! END U-POINT LOOP OVER OBS
    1128. 

  1128. 

  1129.       ENDIF    ! end if(iswind.eq.1)
    1130. 

  1130. 

  1131.     ENDIF   ! ITYP=1: PROCESS U
    1132. 

  1132. 

  1133. !**********************************************************
    1134. 

  1134. !   PROCESS V VARIABLE (IVAR=2)
    1135. 

  1135. !**********************************************************
    1136. 

  1136.     IF(ITYP.EQ.2) THEN
    1137. 

  1137. 

  1138.       IF(ISWIND.EQ.1) THEN
    1139. 

  1139.         DO N=1,NSTA
    1140. 

  1140.           IOB=MAX0(2,IA(N))
    1141. 

  1141.           IOB=MIN0(IOB,ide-1)
    1142. 

  1142.           IOBM=MAX0(1,IOB-1)
    1143. 

  1143.           IOBP=MIN0(ide-1,IOB+1)
    1144. 

  1144.           JOB=MAX0(2,IB(N))
    1145. 

  1145.           JOB=MIN0(JOB,jde-1)
    1146. 

  1146.           JOBM=MAX0(1,JOB-1)
    1147. 

  1147.           JOBP=MIN0(jde-1,JOB+1)
    1148. 

  1148.           KOB=MIN0(K_END,IC(N))
    1149. 

  1149. 

  1150. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1151. 

  1151.           IF(MP_LOCAL_VOBMASK(N))THEN     ! Do if obs on this processor
    1152. 

  1152. #endif
    1153. 

  1153.             KOBP=MIN0(KOB+1,k_end)
    1154. 

  1154.             DXOB=RA(N)
    1155. 

  1155.             DYOB=RB(N)
    1156. 

  1156.             DZOB=RC(N)
    1157. 

  1157. 

  1158. !           Compute surface pressure values at surrounding U and V points
    1159. 

  1159.             PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
    1160. 

  1160.             PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
    1161. 

  1161. 

  1162. !           This is to correct obs to model sigma level using reverse similarity theory
    1163. 

  1163.             if(rko(n).eq.1.0)then
    1164. 

  1164.               vratiob=((1.-DYOB)*((1.-DXOB)*vratio(IOB,JOB)+     &
    1165. 

  1165.                     DXOB*vratio(IOBP,JOB)                        &
    1166. 

  1166.                   )+DYOB*((1.-DXOB)*vratio(IOB,JOBP)+            &
    1167. 

  1167.                   DXOB*vratio(IOBP,JOBP)))
    1168. 

  1168.             else
    1169. 

  1169.               vratiob=1.
    1170. 

  1170.             endif
    1171. 

  1171. !YLIU       Some PBL scheme do not define the vratio/uratio
    1172. 

  1172.             if(abs(vratiob).lt.1.0e-3) then
    1173. 

  1173.               vratiob=1.
    1174. 

  1174.             endif
    1175. 

  1175. 

  1176. !           INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
    1177. 

  1177. !           WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
    1178. 

  1178. !           OUTSIDE THE CURRENT DOMAIN
    1179. 

  1179.   
    1180. 

  1180. !           V COMPONENT WIND ERROR
    1181. 

  1181.             ERRF(2,N)=ERRF(2,N)+vratiob*VAROBS(2,N)-((1.-DZOB)*        &
    1182. 

  1182.                      ((1.-DyOB)*((1.-                                  &
    1183. 

  1183.                       DxOB)*VB(IOB,KOB,JOB)+DxOB*VB(IOB+1,KOB,JOB)     &
    1184. 

  1184.                       )+DyOB*((1.-DxOB)*VB(IOB,KOB,JOB+1)+DxOB*        &
    1185. 

  1185.                       VB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
    1186. 

  1186.                       *VB(IOB,KOBP,JOB)+DxOB*VB(IOB+1,KOBP,JOB))+      &
    1187. 

  1187.                       DyOB*((1.-DxOB)*VB(IOB,KOBP,JOB+1)+DxOB*         &
    1188. 

  1188.                       VB(IOB+1,KOBP,JOB+1))))
    1189. 

  1189. 

  1190. !      if(n.le.10) then
    1191. 

  1191. !        write(6,*) 
    1192. 

  1192. !        write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF2 at ',iob,job,kob,   &
    1193. 

  1193. !                                     ' N = ',n,' inest = ',inest
    1194. 

  1194. !        write(6,*) 'VAROBS(2,N) = ',varobs(2,n)
    1195. 

  1195. !        write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
    1196. 

  1196. !        write(6,*) 'VB(IOB,KOB,JOB) = ',VB(IOB,KOB,JOB)
    1197. 

  1197. !        write(6,*) 'ERRF(2,N) = ',errf(2,n)
    1198. 

  1198. !        write(6,*) 'vratiob = ',vratiob
    1199. 

  1199. !        write(6,*)
    1200. 

  1200. !      endif
    1201. 

  1201. 

  1202.   
    1203. 

  1203. !           Store model surface pressure (not the error!) at V point.
    1204. 

  1204.             ERRF(8,N)=.001*( (1.-DYOB)*PVG1 + DYOB*PVG2 )
    1205. 

  1205.   
    1206. 

  1206. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1207. 

  1207.           ENDIF       ! end IF( MP_LOCAL_VOBMASK(N) )
    1208. 

  1208. #endif
    1209. 

  1209.         ENDDO    ! END V-POINT LOOP OVER OBS
    1210. 

  1210. 

  1211.       ENDIF    ! end if(iswind.eq.1)
    1212. 

  1212. 

  1213.     ENDIF   ! ITYP=1: PROCESS V
    1214. 

  1214. 

  1215. !**********************************************************
    1216. 

  1216. !   PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
    1217. 

  1217. !**********************************************************
    1218. 

  1218.     IF(ITYP.EQ.3) THEN
    1219. 

  1219. 

  1220.       IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
    1221. 

  1221.         DO N=1,NSTA
    1222. 

  1222.           IOB=MAX0(1,IA(N))
    1223. 

  1223.           IOB=MIN0(IOB,ide-1)
    1224. 

  1224.           JOB=MAX0(1,IB(N))
    1225. 

  1225.           JOB=MIN0(JOB,jde-1)
    1226. 

  1226. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1227. 

  1227.           IF(MP_LOCAL_COBMASK(N)) THEN     ! Do if obs on this processor
    1228. 

  1228. #endif
    1229. 

  1229.             KOB=MIN0(k_end,IC(N))
    1230. 

  1230.             KOBP=MIN0(KOB+1,K_END)
    1231. 

  1231.             DXOB=RA(N)
    1232. 

  1232.             DYOB=RB(N)
    1233. 

  1233.             DZOB=RC(N)
    1234. 

  1234. 

  1235. !           This is to correct obs to model sigma level using reverse similarity theory
    1236. 

  1236.             if(rko(n).eq.1.0)then
    1237. 

  1237.               tratxob=((1.-DYOB)*((1.-DXOB)*tratx(IOB,JOB)+        &
    1238. 

  1238.                     DXOB*tratx(IOB+1,JOB)                          &
    1239. 

  1239.                   )+DYOB*((1.-DXOB)*tratx(IOB,JOB+1)+              &
    1240. 

  1240.                   DXOB*tratx(IOB+1,JOB+1)))
    1241. 

  1241.             else
    1242. 

  1242.               tratxob=1.
    1243. 

  1243.             endif
    1244. 

  1244. 

  1245. !yliu
    1246. 

  1246.             if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
    1247. 

  1247. 

  1248. !ajb We must convert temperature to potential temperature
    1249. 

  1249.             obs_pottemp = -888888.
    1250. 

  1250.             if(varobs(3,n).gt.-800000. .and. varobs(5,n).gt.-800000) then
    1251. 

  1251.               obs_pottemp = varobs(3,n)*(100./varobs(5,n))**RCP - t0
    1252. 

  1252.             endif
    1253. 

  1253. 

  1254.             ERRF(3,N)=ERRF(3,N)+tratxob*obs_pottemp-((1.-DZOB)*     &
    1255. 

  1255.                       ((1.-DyOB)*((1.-                              &
    1256. 

  1256.                       DxOB)*(TB(IOB,KOB,JOB))+DxOB*                 &
    1257. 

  1257.                       (TB(IOB+1,KOB,JOB)))+DyOB*((1.-DxOB)*         &
    1258. 

  1258.                       (TB(IOB,KOB,JOB+1))+DxOB*                     &
    1259. 

  1259.                       (TB(IOB+1,KOB,JOB+1))))+DZOB*((1.-            &
    1260. 

  1260.                       DyOB)*((1.-DxOB)*(TB(IOB,KOBP,JOB))+DxOB*     &
    1261. 

  1261.                       (TB(IOB+1,KOBP,JOB)))+DyOB*((1.-DxOB)*        &
    1262. 

  1262.                       (TB(IOB,KOBP,JOB+1))+DxOB*                    &
    1263. 

  1263.                       (TB(IOB+1,KOBP,JOB+1)))))
    1264. 

  1264. 

  1265. !       if(n.le.10) then
    1266. 

  1266. !         write(6,*)
    1267. 

  1267. !         write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob,   &
    1268. 

  1268. !                                      ' N = ',n,' inest = ',inest
    1269. 

  1269. !         write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
    1270. 

  1270. !         write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
    1271. 

  1271. !         write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
    1272. 

  1272. !         write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
    1273. 

  1273. !         write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
    1274. 

  1274. !         write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
    1275. 

  1275. !         write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
    1276. 

  1276. !         write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
    1277. 

  1277. !         write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
    1278. 

  1278. !         write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
    1279. 

  1279. !         write(6,*) 'tratxob = ',tratxob
    1280. 

  1280. !         write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
    1281. 

  1281. !         write(6,*) 'ERRF(3,N) = ',errf(3,n)
    1282. 

  1282. !         write(6,*)
    1283. 

  1283. !       endif
    1284. 

  1284. 

  1285. !           MOISTURE ERROR
    1286. 

  1286.             ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)*((1.-DyOB)*((1.- &
    1287. 

  1287.                       DxOB)*QVB(IOB,KOB,JOB)+DxOB*                      &
    1288. 

  1288.                       QVB(IOB+1,KOB,JOB))+DyOB*((1.-DxOB)*              &
    1289. 

  1289.                       QVB(IOB,KOB,JOB+1)+DxOB*                          &
    1290. 

  1290.                       QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.-                 &
    1291. 

  1291.                       DyOB)*((1.-DxOB)*QVB(IOB,KOBP,JOB)+DxOB           &
    1292. 

  1292.                       *QVB(IOB+1,KOBP,JOB))+DyOB*((1.-DxOB              &
    1293. 

  1293.                       )*QVB(IOB,KOBP,JOB+1)+DxOB*                       &
    1294. 

  1294.                       QVB(IOB+1,KOBP,JOB+1))))
    1295. 

  1295. 

  1296. !           Store model surface pressure (not the error!) at T-point
    1297. 

  1297.             ERRF(6,N)= .001*                                            &
    1298. 

  1298.                       ((1.-DyOB)*((1.-DxOB)*pbase(IOB,1,JOB)+DxOB*      &
    1299. 

  1299.                       pbase(IOB+1,1,JOB))+DyOB*((1.-DxOB)*              &
    1300. 

  1300.                       pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
    1301. 

  1301. 

  1302. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1303. 

  1303.           ENDIF       ! end IF( MP_LOCAL_COBMASK(N) )
    1304. 

  1304. #endif
    1305. 

  1305.         ENDDO     ! END T and QV LOOP OVER OBS
    1306. 

  1306. 

  1307.       ENDIF   !end if(istemp.eq.1 .or. ismois.eq.1)
    1308. 

  1308. 

  1309. !*******************************************************
    1310. 

  1310. !   USE ERRF(5,N) TO STORE KPBL AT (I,J) NEAREST THE OBS
    1311. 

  1311. !*******************************************************
    1312. 

  1312.       DO N=1,NSTA
    1313. 

  1313.         IOB=MAX0(1,IA(N))
    1314. 

  1314.         IOB=MIN0(IOB,ide-1)
    1315. 

  1315.         JOB=MAX0(1,IB(N))
    1316. 

  1316.         JOB=MIN0(JOB,jde-1)
    1317. 

  1317. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1318. 

  1318.         IF(MP_LOCAL_COBMASK(N)) THEN    ! Do if obs on this processor
    1319. 

  1319. #endif
    1320. 

  1320.           DXOB=RA(N)
    1321. 

  1321.           DYOB=RB(N)
    1322. 

  1322.           ERRF(5,N) = kpbl(iob+nint(dxob),job+nint(dyob))
    1323. 

  1323. 

  1324. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1325. 

  1325.         ENDIF       ! end IF( MP_LOCAL_COBMASK(N) )
    1326. 

  1326. #endif
    1327. 

  1327.       ENDDO
    1328. 

  1328. !!**********************************************************
    1329. 

  1329.     ENDIF   ! end if(ityp.eq.3)
    1330. 

  1330. 

  1331.   ENDDO   ! END BIG LOOP
    1332. 

  1332. 

  1333. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1334. 

  1334. ! Gather the errf values calculated by the processors with the obs.
    1335. 

  1335.   CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask,     &
    1336. 

  1336.                            mp_local_vobmask, mp_local_cobmask, errf)
    1337. 

  1337. 

  1338. ! 02252010: Go ahead and assign rko for "obs-off" processors here, to
    1339. 

  1339. !           fix the problem where duplicate obs can be dropped from
    1340. 

  1340. !           the "obs-on" processor, but not from the others, due to
    1341. 

  1341. !           rko being -99 on the "obs-off" processors.
    1342. 

  1342.   do n=1,nsta
    1343. 

  1343.     rko(n) = errf(9,n)
    1344. 

  1344.   enddo
    1345. 

  1345. ! 02252010: End bugfix.
    1346. 

  1346. #endif
    1347. 

  1347. 

  1348. ! Print obs information.
    1349. 

  1349.   CALL print_obs_info(iprt,inest,niobf,rio,rjo,rko,                  &
    1350. 

  1350.                       prt_max,prt_freq,obs_prt,stnid_prt,            &
    1351. 

  1351.                       lat_prt,lon_prt,mlat_prt,mlon_prt,             &
    1352. 

  1352.                       timeob,xtime)
    1353. 

  1353. 

  1354. ! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
    1355. 

  1355.   IF(INEST.EQ.1)THEN
    1356. 

  1356.     INPF=NPFI
    1357. 

  1357.   ELSE
    1358. 

  1358.     FNPF=IRATIO**LEVIDN(INEST)
    1359. 

  1359.     INPF=FNPF*NPFI
    1360. 

  1360.   ENDIF
    1361. 

  1361. ! Gross error check for temperature. Set all vars bad.
    1362. 

  1362.   do n=1,nsta
    1363. 

  1363.     if((abs(errf(3,n)).gt.20.).and.           &
    1364. 

  1364.            (errf(3,n).gt.-800000.))then
    1365. 

  1365. 

  1366.        errf(1,n)=-888888.
    1367. 

  1367.        errf(2,n)=-888888.
    1368. 

  1368.        errf(3,n)=-888888.
    1369. 

  1369.        errf(4,n)=-888888.
    1370. 

  1370.        varobs(1,n)=-888888.
    1371. 

  1371.        varobs(2,n)=-888888.
    1372. 

  1372.        varobs(3,n)=-888888.
    1373. 

  1373.        varobs(4,n)=-888888.
    1374. 

  1374.     endif
    1375. 

  1375.   enddo
    1376. 

  1376. 

  1377. ! For printout
    1378. 

  1378. !  IF(MOD(KTAU,INPF).NE.0) THEN
    1379. 

  1379. !      RETURN
    1380. 

  1380. !  ENDIF
    1381. 

  1381. 

  1382.   RETURN
    1383. 

  1383.   END SUBROUTINE errob
    1384. 

  1384. 

  1385.   SUBROUTINE upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme,  &
    1386. 

  1386.                     arrin, arrout)
    1387. 

  1387. !------------------------------------------------------------------------------
    1388. 

  1388. !     PURPOSE: This subroutine interpolates a real 2D array defined over mass
    1389. 

  1389. !              coordinate points, to U (momentum) points.
    1390. 

  1390. !
    1391. 

  1391. !------------------------------------------------------------------------------
    1392. 

  1392.   IMPLICIT NONE
    1393. 

  1393. 

  1394.   INTEGER, INTENT(IN) :: i_start     ! Starting i index for this model tile
    1395. 

  1395.   INTEGER, INTENT(IN) :: i_end       ! Ending   i index for this model tile
    1396. 

  1396.   INTEGER, INTENT(IN) :: j_start     ! Starting j index for this model tile
    1397. 

  1397.   INTEGER, INTENT(IN) :: j_end       ! Ending   j index for this model tile
    1398. 

  1398.   INTEGER, INTENT(IN) :: ids         ! Starting i index for entire model domain
    1399. 

  1399.   INTEGER, INTENT(IN) :: ide         ! Ending   i index for entire model domain
    1400. 

  1400.   INTEGER, INTENT(IN) :: ims         ! Starting i index for model patch 
    1401. 

  1401.   INTEGER, INTENT(IN) :: ime         ! Ending   i index for model patch 
    1402. 

  1402.   INTEGER, INTENT(IN) :: jms         ! Starting j index for model patch 
    1403. 

  1403.   INTEGER, INTENT(IN) :: jme         ! Ending   j index for model patch 
    1404. 

  1404.   REAL,   INTENT(IN)  :: arrin ( ims:ime, jms:jme )  ! input array on mass points
    1405. 

  1405.   REAL,   INTENT(OUT) :: arrout( ims:ime, jms:jme )  ! output array on U points 
    1406. 

  1406. 

  1407. ! Local variables
    1408. 

  1408.   integer :: i, j
    1409. 

  1409. 

  1410. ! Do domain interior first
    1411. 

  1411.   do j = j_start, j_end
    1412. 

  1412.     do i = max(2,i_start), i_end
    1413. 

  1413.        arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
    1414. 

  1414.     enddo
    1415. 

  1415.   enddo
    1416. 

  1416. 

  1417. ! Do west-east boundaries
    1418. 

  1418.   if(i_start .eq. ids) then
    1419. 

  1419.     do j = j_start, j_end
    1420. 

  1420.       arrout(i_start,j) = arrin(i_start,j)
    1421. 

  1421.     enddo
    1422. 

  1422.   endif
    1423. 

  1423.   if(i_end .eq. ide-1) then
    1424. 

  1424.     do j = j_start, j_end
    1425. 

  1425.       arrout(i_end+1,j) = arrin(i_end,j)
    1426. 

  1426.     enddo
    1427. 

  1427.   endif
    1428. 

  1428. 

  1429.   RETURN
    1430. 

  1430.   END SUBROUTINE upoint
    1431. 

  1431. 

  1432.   SUBROUTINE vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme,  &
    1433. 

  1433.                     arrin, arrout)
    1434. 

  1434. !------------------------------------------------------------------------------
    1435. 

  1435. !     PURPOSE: This subroutine interpolates a real 2D array defined over mass
    1436. 

  1436. !              coordinate points, to V (momentum) points.
    1437. 

  1437. !
    1438. 

  1438. !------------------------------------------------------------------------------
    1439. 

  1439.   IMPLICIT NONE
    1440. 

  1440. 

  1441.   INTEGER, INTENT(IN) :: i_start     ! Starting i index for this model tile
    1442. 

  1442.   INTEGER, INTENT(IN) :: i_end       ! Ending   i index for this model tile
    1443. 

  1443.   INTEGER, INTENT(IN) :: j_start     ! Starting j index for this model tile
    1444. 

  1444.   INTEGER, INTENT(IN) :: j_end       ! Ending   j index for this model tile
    1445. 

  1445.   INTEGER, INTENT(IN) :: jds         ! Starting j index for entire model domain
    1446. 

  1446.   INTEGER, INTENT(IN) :: jde         ! Ending   j index for entire model domain
    1447. 

  1447.   INTEGER, INTENT(IN) :: ims         ! Starting i index for model patch 
    1448. 

  1448.   INTEGER, INTENT(IN) :: ime         ! Ending   i index for model patch 
    1449. 

  1449.   INTEGER, INTENT(IN) :: jms         ! Starting j index for model patch 
    1450. 

  1450.   INTEGER, INTENT(IN) :: jme         ! Ending   j index for model patch 
    1451. 

  1451.   REAL,   INTENT(IN)  :: arrin ( ims:ime, jms:jme )  ! input array on mass points
    1452. 

  1452.   REAL,   INTENT(OUT) :: arrout( ims:ime, jms:jme )  ! output array on V points 
    1453. 

  1453. 

  1454. ! Local variables
    1455. 

  1455.   integer :: i, j
    1456. 

  1456. 

  1457. ! Do domain interior first
    1458. 

  1458.   do j = max(2,j_start), j_end
    1459. 

  1459.     do i = i_start, i_end
    1460. 

  1460.       arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
    1461. 

  1461.     enddo
    1462. 

  1462.   enddo
    1463. 

  1463. 

  1464. ! Do south-north boundaries
    1465. 

  1465.   if(j_start .eq. jds) then
    1466. 

  1466.     do i = i_start, i_end
    1467. 

  1467.       arrout(i,j_start) = arrin(i,j_start)
    1468. 

  1468.     enddo
    1469. 

  1469.   endif
    1470. 

  1470.   if(j_end .eq. jde-1) then
    1471. 

  1471.     do i = i_start, i_end
    1472. 

  1472.       arrout(i,j_end+1) = arrin(i,j_end)
    1473. 

  1473.     enddo
    1474. 

  1474.   endif
    1475. 

  1475. 

  1476.   RETURN
    1477. 

  1477.   END SUBROUTINE vpoint
    1478. 

  1478. 

  1479.   LOGICAL FUNCTION TILE_MASK(iloc, jloc, its, ite, jts, jte)
    1480. 

  1480. !------------------------------------------------------------------------------
    1481. 

  1481. ! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
    1482. 

  1482. !      
    1483. 

  1483. ! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
    1484. 

  1484. !                  tile-range indices (its,jts) and (ite,jte)
    1485. 

  1485. !          FALSE otherwise.
    1486. 

  1486. !
    1487. 

  1487. !------------------------------------------------------------------------------
    1488. 

  1488.   IMPLICIT NONE
    1489. 

  1489. 

  1490.   INTEGER, INTENT(IN) :: iloc
    1491. 

  1491.   INTEGER, INTENT(IN) :: jloc
    1492. 

  1492.   INTEGER, INTENT(IN) :: its
    1493. 

  1493.   INTEGER, INTENT(IN) :: ite
    1494. 

  1494.   INTEGER, INTENT(IN) :: jts
    1495. 

  1495.   INTEGER, INTENT(IN) :: jte
    1496. 

  1496. 

  1497. ! Local variables
    1498. 

  1498.   LOGICAL :: retval
    1499. 

  1499. 

  1500.   TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND.    &
    1501. 

  1501.                jloc .LE. jte .AND. jloc .GE. jts )
    1502. 

  1502. 

  1503.   RETURN
    1504. 

  1504.   END FUNCTION TILE_MASK
    1505. 

  1505. 

  1506. !-----------------------------------------------------------------------
    1507. 

  1507.   SUBROUTINE nudob(j, ivar, aten, inest, ifrest, ktau, ktaur,         &
    1508. 

  1508.                        xtime, mu, msfx, msfy, nndgv, nerrf, niobf, maxdom,   &
    1509. 

  1509.                        npfi, ionf, rinxy, twindo,                     &
    1510. 

  1510.                        nudge_pbl,                                     &
    1511. 

  1511.                        sfcfact, sfcfacr,                              &
    1512. 

  1512.                        levidn,                                        &
    1513. 

  1513.                        parid, nstat,                                  &
    1514. 

  1514.                        rinfmn, rinfmx, pfree, dcon, tfaci,            &
    1515. 

  1515.                        sfc_scheme_horiz, sfc_scheme_vert, maxsnd_gap, &
    1516. 

  1516.                        lev_in_ob, plfo, nlevs_ob,                     &
    1517. 

  1517.                        iratio, dx, dtmin, rio, rjo, rko,              &
    1518. 

  1518.                        timeob, varobs, errf, pbase, ptop, pp,         &
    1519. 

  1519.                        iswind, istemp, ismois, giv, git, giq,         &
    1520. 

  1520.                        savwt, kpblt, nscan,                           &
    1521. 

  1521.                        vih1, vih2, terrh, zslab,                      &
    1522. 

  1522.                        iprt,                                          &
    1523. 

  1523.                        ids,ide, jds,jde, kds,kde,                     &  ! domain dims
    1524. 

  1524.                        ims,ime, jms,jme, kms,kme,                     &  ! memory dims
    1525. 

  1525.                        its,ite, jts,jte, kts,kte )                       ! tile   dims
    1526. 

  1526. 

  1527. !-----------------------------------------------------------------------
    1528. 

  1528.   USE module_model_constants
    1529. 

  1529.   USE module_domain
    1530. 

  1530. !-----------------------------------------------------------------------
    1531. 

  1531.   IMPLICIT NONE
    1532. 

  1532. !-----------------------------------------------------------------------
    1533. 

  1533. !
    1534. 

  1534. ! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
    1535. 

  1535. !     VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
    1536. 

  1536. !     ASSIMILATION FROM INDIVIDUAL OBSERVATIONS.  THE NUDGING
    1537. 

  1537. !     TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
    1538. 

  1538. !     WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
    1539. 

  1539. !     ANALYSIS.  THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
    1540. 

  1540. !     AND MINIMAL STORAGE REQUIREMENTS.  ALGORITHMS SHOULD BE
    1541. 

  1541. !     VECTORIZED WHEREVER POSSIBLE.
    1542. 

  1542. !
    1543. 

  1543. !     HISTORY: Original author: MM5 version???
    1544. 

  1544. !              02/04/2004 - Creation of WRF version.           Al Bourgeois
    1545. 

  1545. !              08/28/2006 - Conversion from F77 to F90         Al Bourgeois
    1546. 

  1546. !------------------------------------------------------------------------------
    1547. 

  1547. !
    1548. 

  1548. ! NOTE: This routine was originally designed for MM5, which uses
    1549. 

  1549. !       a nonstandard (I,J) coordinate system. For WRF, I is the 
    1550. 

  1550. !       east-west running coordinate, and J is the south-north
    1551. 

  1551. !       running coordinate. So a "J-slab" here is west-east in
    1552. 

  1552. !       extent, not south-north as for MM5.      -ajb 06/10/2004
    1553. 

  1553. !
    1554. 

  1554. !     NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
    1555. 

  1555. !     AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
    1556. 

  1556. !     TYPES OF FACTORS:
    1557. 

  1557. !       1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
    1558. 

  1558. !          TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
    1559. 

  1559. !          TIME (XTIME) ARE USED.  OBSERVATIONS CLOSEST TO
    1560. 

  1560. !          XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
    1561. 

  1561. !       2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
    1562. 

  1562. !          CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
    1563. 

  1563. !          (RINSIG).
    1564. 

  1564. !       3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
    1565. 

  1565. !          CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY).  THE
    1566. 

  1566. !          VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
    1567. 

  1567. !          TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
    1568. 

  1568. !
    1569. 

  1569. !     THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
    1570. 

  1570. !     VALUE OF IVAR AS FOLLOWS:
    1571. 

  1571. !             IVAR                     VARIABLE(TAU-1)
    1572. 

  1572. !             ----                     ---------------
    1573. 

  1573. !               1                             U
    1574. 

  1574. !               2                             V
    1575. 

  1575. !               3                             T
    1576. 

  1576. !               4                             QV
    1577. 

  1577. !               5                          PSB(CROSS)   REMOVED IN V3
    1578. 

  1578. !              (6)                           PSB(DOT)
    1579. 

  1579. !
    1580. 

  1580. !-----------------------------------------------------------------------
    1581. 

  1581. !
    1582. 

  1582. !     Description of input arguments.
    1583. 

  1583. !
    1584. 

  1584. !-----------------------------------------------------------------------
    1585. 

  1585. 

  1586.   INTEGER, INTENT(IN)  :: ids,ide, jds,jde, kds,kde  ! domain dims.
    1587. 

  1587.   INTEGER, INTENT(IN)  :: ims,ime, jms,jme, kms,kme  ! memory dims.
    1588. 

  1588.   INTEGER, INTENT(IN)  :: its,ite, jts,jte, kts,kte  ! tile   dims.
    1589. 

  1589.   INTEGER, INTENT(IN)  :: j                          ! south-north running coordinate.
    1590. 

  1590.   INTEGER, INTENT(IN)  :: ivar
    1591. 

  1591.   INTEGER, INTENT(IN)  :: inest                      ! domain index
    1592. 

  1592.   LOGICAL, INTENT(IN)  :: ifrest
    1593. 

  1593.   INTEGER, INTENT(IN)  :: ktau
    1594. 

  1594.   INTEGER, INTENT(IN)  :: ktaur
    1595. 

  1595.   REAL, INTENT(IN)     :: xtime                      ! forecast time in minutes
    1596. 

  1596.   INTEGER, INTENT(IN)  :: nndgv                      ! number of nudge variables
    1597. 

  1597.   INTEGER, INTENT(IN)  :: nerrf                      ! number of error fields
    1598. 

  1598.   INTEGER, INTENT(IN)  :: niobf                      ! number of observations
    1599. 

  1599.   INTEGER, INTENT(IN)  :: maxdom                     ! maximum number of domains
    1600. 

  1600.   INTEGER, INTENT(IN)  :: npfi 
    1601. 

  1601.   INTEGER, INTENT(IN)  :: ionf
    1602. 

  1602.   REAL, INTENT(IN)     :: rinxy
    1603. 

  1603.   REAL, INTENT(IN)     :: twindo
    1604. 

  1604.   REAL, intent(in)     :: sfcfact                    ! scale for time window (surface obs) 
    1605. 

  1605.   REAL, intent(in)     :: sfcfacr                    ! scale for hor rad inf (surface obs)
    1606. 

  1606.   LOGICAL, intent(in)  :: nudge_pbl                  ! flag for nudging in pbl
    1607. 

  1607.   INTEGER, INTENT(IN)  :: levidn(maxdom)             ! level of nest
    1608. 

  1608.   INTEGER, INTENT(IN)  :: parid(maxdom)              ! parent domain id
    1609. 

  1609.   INTEGER, INTENT(IN)  :: nstat                      ! number of obs stations
    1610. 

  1610.   REAL,    INTENT(IN)  :: rinfmn          ! minimum radius of influence
    1611. 

  1611.   REAL,    INTENT(IN)  :: rinfmx          ! maximum radius of influence
    1612. 

  1612.   REAL,    INTENT(IN)  :: pfree           ! pressure level (cb) where terrain effect becomes small
    1613. 

  1613.   REAL,    INTENT(IN)  :: dcon            ! 1/DPSMX
    1614. 

  1614.   REAL,    INTENT(IN)  :: tfaci           ! scale factor used for ramp-down in dynamic initialization
    1615. 

  1615.   INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
    1616. 

  1616.   INTEGER , INTENT(IN) :: sfc_scheme_vert  ! vertical   spreading scheme for surf obs (orig or regime vif)
    1617. 

  1617.   REAL    , INTENT(IN) :: maxsnd_gap      ! max allowed pressure gap in soundings, for interp (centibars)
    1618. 

  1618.   REAL, INTENT(IN)     :: lev_in_ob(niobf)           ! Level in sounding-type obs.
    1619. 

  1619.   REAL, intent(IN)     :: plfo(niobf)                ! Index for type of obs platform
    1620. 

  1620.   REAL, INTENT(IN)     :: nlevs_ob(niobf)            ! Number of levels in sounding.
    1621. 

  1621.   INTEGER, INTENT(IN)  :: iratio                     ! Nest to parent gridsize ratio.
    1622. 

  1622.   REAL, INTENT(IN)     :: dx                         ! This domain grid cell-size (m)
    1623. 

  1623.   REAL, INTENT(IN)     :: dtmin                      ! Model time step in minutes
    1624. 

  1624.   REAL, INTENT(IN)     :: rio(niobf)                 ! Obs west-east coordinate (non-stag grid).
    1625. 

  1625.   REAL, INTENT(IN)     :: rjo(niobf)                 ! Obs south-north coordinate (non-stag grid).
    1626. 

  1626.   REAL, INTENT(INOUT)  :: rko(niobf)                 ! Obs vertical coordinate.
    1627. 

  1627.   REAL, INTENT(IN)     :: timeob(niobf)
    1628. 

  1628.   REAL, INTENT(IN)     :: varobs(nndgv,niobf)
    1629. 

  1629.   REAL, INTENT(IN)     :: errf(nerrf, niobf)
    1630. 

  1630.   REAL, INTENT(IN)     :: pbase( ims:ime, kms:kme )  ! Base pressure.
    1631. 

  1631.   REAL, INTENT(IN)     :: ptop
    1632. 

  1632.   REAL, INTENT(IN)     :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
    1633. 

  1633.   REAL, INTENT(IN)     :: mu(ims:ime)   ! Air mass on u, v, or mass-grid
    1634. 

  1634.   REAL, INTENT(IN)     :: msfx(ims:ime)  ! Map scale (only used for vars u & v)
    1635. 

  1635.   REAL, INTENT(IN)     :: msfy(ims:ime)  ! Map scale (only used for vars u & v)
    1636. 

  1636.   INTEGER, intent(in)  :: iswind        ! Nudge flag for wind
    1637. 

  1637.   INTEGER, intent(in)  :: istemp        ! Nudge flag for temperature
    1638. 

  1638.   INTEGER, intent(in)  :: ismois        ! Nudge flag for moisture
    1639. 

  1639.   REAL, intent(in)     :: giv           ! Coefficient for wind
    1640. 

  1640.   REAL, intent(in)     :: git           ! Coefficient for temperature
    1641. 

  1641.   REAL, intent(in)     :: giq           ! Coefficient for moisture
    1642. 

  1642.   REAL, INTENT(INOUT)  :: aten( ims:ime, kms:kme)
    1643. 

  1643.   REAL, INTENT(INOUT)  :: savwt( nndgv, ims:ime, kms:kme )
    1644. 

  1644.   INTEGER, INTENT(IN)  :: kpblt(ims:ime)
    1645. 

  1645.   INTEGER, INTENT(IN)  :: nscan                      ! number of scans
    1646. 

  1646.   REAL, INTENT(IN)     :: vih1(its:ite) ! Vert infl ht (m) abv grd for full wts
    1647. 

  1647.   REAL, INTENT(IN)     :: vih2(its:ite) ! Vert infl ht (m) abv grd for ramp
    1648. 

  1648.   REAL, INTENT(IN)     :: terrh(ims:ime) ! Terrain height (m)
    1649. 

  1649. ! INTEGER, INTENT(IN)  :: vik1(its:ite) ! Vertical infl k-level for full wts
    1650. 

  1650. ! INTEGER, INTENT(IN)  :: vik2(its:ite) ! Vertical infl k-level for ramp
    1651. 

  1651.   REAL, INTENT(IN)     :: zslab(ims:ime, kms:kme)    ! model ht above ground (m)
    1652. 

  1652.   LOGICAL, INTENT(IN)  :: iprt                       ! print flag
    1653. 

  1653. 

  1654. ! Local variables
    1655. 

  1655.   integer :: mm(maxdom)
    1656. 

  1656.   integer :: kobs                  ! k-lev below obs (for obs straddling pblt)
    1657. 

  1657.   integer :: kpbl_obs(nstat)       ! kpbl at grid point (IOB,JOB) (ajb 20090519)
    1658. 

  1658.   real :: ra(niobf)
    1659. 

  1659.   real :: rb(niobf)
    1660. 

  1660.   real :: psurf(niobf)
    1661. 

  1661.   real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
    1662. 

  1662.   real :: rscale(ims:ime)           ! For converting to rho-coupled units.
    1663. 

  1663.   real :: wtij(ims:ime)             ! For holding weights in i-loop
    1664. 

  1664.   real :: reserf(100)
    1665. 

  1665.   character*40 name
    1666. 

  1666.   character*3 chr_hr
    1667. 

  1667.   character(len=200) :: msg            ! Argument to wrf_message
    1668. 

  1668. 

  1669. !***  DECLARATIONS FOR IMPLICIT NONE
    1670. 

  1670.   integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
    1671. 

  1671.   integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
    1672. 

  1672.   integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
    1673. 

  1673.   integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
    1674. 

  1674.   integer :: komin,komax,nn,nhi,nlo,nnjc
    1675. 

  1675.   integer :: i_s,i_e
    1676. 

  1676.   integer :: istq
    1677. 

  1677.   real :: gfactor,rfactor,gridx,gridy,rindx,ris
    1678. 

  1678.   real :: grfacx,grfacy
    1679. 

  1679.   real :: timewt,pob
    1680. 

  1680.   real :: ri,rj,rx,ry,rsq,pdfac,erfivr,dk,slope,rinfac
    1681. 

  1681.   real :: dprim,dsq,d     ! vars for mm5 surface-ob weighting
    1682. 

  1682.   real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
    1683. 

  1683.   real :: dz_ramp         ! For ramping weights for surface obs 
    1684. 

  1684. 

  1685.   real :: scratch
    1686. 

  1686.   integer :: kk !ajb temp
    1687. 

  1687. 

  1688. !      print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
    1689. 

  1689. !      if(ivar.ne.4)return
    1690. 

  1690. !yliu start -- for multi-scans: NSCAN=0: original
    1691. 

  1691. !                               NSCAN=1: added a scan with a larger Ri and smaller G
    1692. 

  1692. !       if(NSCAN.ne.0 .and. NSCAN.ne.1)  stop
    1693. 

  1693. ! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
    1694. 

  1694.   if(NSCAN.ne.0) then
    1695. 

  1695.     IF (iprt) then
    1696. 

  1696.         write(msg,*) 'SAVWT must be resized for NSCAN=1'
    1697. 

  1697.         call wrf_message(msg)
    1698. 

  1698.     ENDIF
    1699. 

  1699.     call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
    1700. 

  1700.   endif
    1701. 

  1701.   IPLO=0  + NSCAN*4
    1702. 

  1702.   GFACTOR=1. +  NSCAN*(-1. + 0.33333) 
    1703. 

  1703.   RFACTOR=1. +  NSCAN*(-1. + 3.0)
    1704. 

  1704. !yliu end
    1705. 

  1705. ! jc
    1706. 

  1706. 

  1707. ! return if too close to j boundary
    1708. 

  1708.   if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
    1709. 

  1709. !       write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
    1710. 

  1710. !    $             ' too close to boundary.'
    1711. 

  1711.     return
    1712. 

  1712.   endif
    1713. 

  1713.   if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
    1714. 

  1714. !       write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
    1715. 

  1715. !    $             ' too close to boundary.'
    1716. 

  1716.     return
    1717. 

  1717.   endif
    1718. 

  1718. 

  1719. ! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
    1720. 

  1720.   ICUT=0
    1721. 

  1721.   IF(INEST.GT.1)ICUT=1
    1722. 

  1722.   i_s = max0(2+icut,its)
    1723. 

  1723.   i_e = min0(ide-1-icut,ite)
    1724. 

  1724. 

  1725.   IPL=IVAR    + IPLO     !yliu +IPLO
    1726. 

  1726. 

  1727. ! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
    1728. 

  1728. 

  1729.   INPF=(IRATIO**LEVIDN(INEST))*NPFI
    1730. 

  1730.   INFR=(IRATIO**LEVIDN(INEST))*IONF
    1731. 

  1731. 

  1732.   GRIDX=0.0
    1733. 

  1733.   GRIDY=0.0
    1734. 

  1734.   IGRID=0
    1735. 

  1735.   IF(IVAR.GE.3)THEN
    1736. 

  1736.     GRIDX=0.5
    1737. 

  1737.     GRIDY=0.5
    1738. 

  1738.     IGRID=1
    1739. 

  1739.   ELSEIF(IVAR.eq.1) THEN
    1740. 

  1740.     GRIDY=0.5
    1741. 

  1741.     GRIDX=0.0
    1742. 

  1742.     IGRID=1
    1743. 

  1743.   ELSEIF(IVAR.eq.2) THEN
    1744. 

  1744.     GRIDX=0.5
    1745. 

  1745.     GRIDY=0.0
    1746. 

  1746.     IGRID=1
    1747. 

  1747.   ENDIF
    1748. 

  1748. 

  1749. ! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
    1750. 

  1750. ! KILOMETERS TO GRID LENGTHS, RINDX
    1751. 

  1751. 

  1752.   RINDX=RINXY*1000./DX          * RFACTOR   !yliu *RFACTOR
    1753. 

  1753.   RIS=RINDX*RINDX
    1754. 

  1754.   IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
    1755. 

  1755.   IF(MOD(KTAU,INFR).NE.0)GOTO 126
    1756. 

  1756. 5 CONTINUE
    1757. 

  1757.   IF (iprt) THEN
    1758. 

  1758.    IF(J.EQ.10) then
    1759. 

  1759.        write(msg,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
    1760. 

  1760.        call wrf_message(msg)
    1761. 

  1761.    ENDIF
    1762. 

  1762.   ENDIF
    1763. 

  1763. 6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,',                    &
    1764. 

  1764.             'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2,       &
    1765. 

  1765.             ' rindx=',f4.1)
    1766. 

  1766. 

  1767. !********************************************************************
    1768. 

  1768. !ajb_07012008  Setting ra and rb for the fine mesh her is now simple.
    1769. 

  1769. !              Values are no longer calculated here based on the
    1770. 

  1770. !              coarse mesh, since direct use of WRF map projections
    1771. 

  1771. !              on each nest was implemented in subroutine in4dob.
    1772. 

  1772. !********************************************************************
    1773. 

  1773. ! SET RA AND RB
    1774. 

  1774.   DO N=1,NSTAT
    1775. 

  1775.     RA(N)=RIO(N)-GRIDX
    1776. 

  1776.     RB(N)=RJO(N)-GRIDY
    1777. 

  1777.   ENDDO
    1778. 

  1778. 

  1779. ! INITIALIZE WEIGHTING ARRAYS TO ZERO
    1780. 

  1780.   DO I=its,ite
    1781. 

  1781.     DO K=1,kte
    1782. 

  1782.       WT(I,K)=0.0
    1783. 

  1783.       WT2ERR(I,K)=0.0
    1784. 

  1784.     ENDDO
    1785. 

  1785.   ENDDO
    1786. 

  1786. 

  1787. ! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
    1788. 

  1788. ! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
    1789. 

  1789. !
    1790. 

  1790. ! COMPUTE P* AT OBS LOCATION (RA,RB).  DO THIS AS SEPARATE VECTOR LOOP H
    1791. 

  1791. ! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
    1792. 

  1792. 

  1793. ! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
    1794. 

  1794. ! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
    1795. 

  1795. ! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
    1796. 

  1796. ! ajb 05052009: psurf is actually pbase(k=1) interpolated to obs (i,j).
    1797. 

  1797.   DO N=1,NSTAT
    1798. 

  1798.     IF(IVAR.GE.3)THEN
    1799. 

  1799.       PSURF(N)=ERRF(6,N)
    1800. 

  1800.     ELSE
    1801. 

  1801.       IF(IVAR.EQ.1)THEN
    1802. 

  1802.         PSURF(N)=ERRF(7,N)        ! U-points
    1803. 

  1803.       ELSE
    1804. 

  1804.         PSURF(N)=ERRF(8,N)        ! V-points
    1805. 

  1805.       ENDIF
    1806. 

  1806.     ENDIF
    1807. 

  1807.   ENDDO
    1808. 

  1808. 

  1809. ! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
    1810. 

  1810. ! J-STRIP
    1811. 

  1811. 

  1812.   MAXJ=J+IFIX(RINDX*RINFMX+0.99)                                             !ajb
    1813. 

  1813.   MINJ=J-IFIX(RINDX*RINFMX+0.99)                                             !ajb
    1814. 

  1814. 

  1815. ! jc comment out this? want to use obs beyond the domain?
    1816. 

  1816. !      MAXJ=MIN0(JL-IGRID,MAXJ)           !yliu
    1817. 

  1817. !      MINJ=MAX0(1,MINJ)                  !yliu
    1818. 

  1818. 

  1819.   n=1
    1820. 

  1820. 

  1821. !***********************************************************************
    1822. 

  1822.   DO nnn=1,NSTAT   ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
    1823. 

  1823. !***********************************************************************
    1824. 

  1824. ! Soundings are consecutive obs, but they need to be treated as a single 
    1825. 

  1825. ! entity. Thus change the looping to nnn, with n defined separately.
    1826. 

  1826. 

  1827. 

  1828. !yliu 
    1829. 

  1829. !  note for sfc data: nsndlev=1 and njcsnd=1
    1830. 

  1830.     nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1   
    1831. 

  1831. 

  1832. ! yliu start -- set together with the other parts
    1833. 

  1833. ! test: do the sounding levels as individual obs
    1834. 

  1834. !   nsndlev=1
    1835. 

  1835. ! yliu end
    1836. 

  1836.     njcsnd=nsndlev
    1837. 

  1837. ! set pob here, to be used later
    1838. 

  1838.     pob=varobs(5,n)
    1839. 

  1839. 

  1840. !     print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
    1841. 

  1841. !     print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
    1842. 

  1842. ! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
    1843. 

  1843. ! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP.  THIS
    1844. 

  1844. ! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
    1845. 

  1845. ! ATION.
    1846. 

  1846. 

  1847. !yliu: Skip bad obs if it is sfc or single level sounding.
    1848. 

  1848. !yliu: Before this (020918), a snd will be skipped if its first 
    1849. 

  1849. !yliu        level has bad data- often true due to elevation
    1850. 

  1850. 

  1851.     IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
    1852. 

  1852. !     print *, " bad obs skipped"
    1853. 

  1853. 

  1854.     ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
    1855. 

  1855. !     print *, " skipped obs far away from this J-slice"
    1856. 

  1856. 

  1857. !----------------------------------------------------------------------
    1858. 

  1858.     ELSE    ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
    1859. 

  1859. !----------------------------------------------------------------------
    1860. 

  1860. 

  1861. ! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
    1862. 

  1862. ! FOR THE VERTICAL WEIGHTING, WTSIG
    1863. 

  1863. 

  1864. ! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
    1865. 

  1865. !ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
    1866. 

  1866. !ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
    1867. 

  1867. 

  1868. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    1869. 

  1869.       rko(n) = errf(9,n)        !ajb 20021210
    1870. 

  1870.       kpbl_obs(n) = errf(5,n)   !ajb 20090519
    1871. 

  1871. #endif
    1872. 

  1872. !ajb 20090427 added .45 to rko so KOB is equal to 1 only if RKO > 1.05
    1873. 

  1873.       KOB=nint(RKO(N)+0.45)
    1874. 

  1874.       KOB=MIN0(kte,KOB)
    1875. 

  1875.       KOB=MAX0(1,KOB)
    1876. 

  1876. 

  1877. ! ASSIMILATE SURFACE LAYER DATA ON SIGMA
    1878. 

  1878.       IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
    1879. 

  1879. 

  1880. ! Compute temporal weight
    1881. 

  1881.         timewt = get_timewt(xtime,dtmin,twindo,sfcfact,timeob(n))
    1882. 

  1882. 

  1883.         DO K=1,kte
    1884. 

  1884.           WTSIG(K)=0.0
    1885. 

  1885.         ENDDO
    1886. 

  1886. ! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
    1887. 

  1887. !       WTSIG(1)=1.0
    1888. 

  1888. !       WTSIG(2)=0.67
    1889. 

  1889. !       WTSIG(3)=0.33
    1890. 

  1890. !       KOMIN=3
    1891. 

  1891. !       KOMAX=1
    1892. 

  1892. ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
    1893. 

  1893. ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
    1894. 

  1894. ! fix this because kpblt at 1 and il is 0
    1895. 

  1895.         MAXI=IFIX(RA(N)+0.99+RINDX*sfcfacr)
    1896. 

  1896.         MAXI=MIN0(ide-1,MAXI)
    1897. 

  1897.         MINI=IFIX(RA(N)-RINDX*sfcfacr-0.99)
    1898. 

  1898.         MINI=MAX0(2,MINI)
    1899. 

  1899. !yliu start
    1900. 

  1900. ! use also obs outside of this domain  -- surface obs
    1901. 

  1901. !     if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
    1902. 

  1902. !    &   RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
    1903. 

  1903. !        print *, " skipped obs far away from this domain"
    1904. 

  1904. ! currently can use obs within this domain or ones very close to (1/3
    1905. 

  1905. !   influence of radius in the coarse domain) this
    1906. 

  1906. !   domain. In later case, use BC column value to approximate the model value
    1907. 

  1907. !   at obs point -- ERRF need model field in errob.F !!
    1908. 

  1908.         if (  RA(N).GE.(0.-RINDX*sfcfacr/3)                        &
    1909. 

  1909.         .and. RA(N).LE.float(ide)+RINDX*sfcfacr/3                  &
    1910. 

  1910.         .and. RB(N).GE.(0.-RINDX*sfcfacr/3)                        &
    1911. 

  1911.         .and. RB(N).LE.float(jde)+RINDX*sfcfacr/3) then
    1912. 

  1912. 

  1913. ! or use obs within this domain only
    1914. 

  1914. !     if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
    1915. 

  1915. !    &   RB(N).LT.1 .or. RB(N).GT.float(JL)) then
    1916. 

  1916. !        print *, " skipped obs far outside of this domain"
    1917. 

  1917. !        if(j.eq.3 .and. ivar.eq.3) then
    1918. 

  1918. !           write(6,*) 'N = ',n,' exit 120 3'
    1919. 

  1919. !        endif
    1920. 

  1920. !yliu end
    1921. 

  1921. !
    1922. 

  1922. ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
    1923. 

  1923. ! OBSERVATION N.  COMPUTE THE HORIZONTAL DISTANCE TO
    1924. 

  1924. ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
    1925. 

  1925.           RJ=FLOAT(J)
    1926. 

  1926.           RX=RJ-RB(N)
    1927. 

  1927. ! WEIGHTS FOR THE 3-D VARIABLES
    1928. 

  1928.           ERFIVR=ERRF(IVAR,N)
    1929. 

  1929.  
    1930. 

  1930. !ajb Compute and add weights to sum only if nudge_pbl switch is on.
    1931. 

  1931.           if(nudge_pbl) then
    1932. 

  1932. 

  1933. ! Apply selected horizontal spreading scheme.
    1934. 

  1934.             if(SFC_SCHEME_HORIZ.eq.1) then        ! original mm5 scheme 
    1935. 

  1935.               DO I=max0(its,MINI),min0(ite,MAXI)
    1936. 

  1936.                 RI=FLOAT(I)
    1937. 

  1937.                 RY=RI-RA(N)
    1938. 

  1938.                 RIS=RINDX*RINDX*sfcfacr*sfcfacr
    1939. 

  1939.                 RSQ=RX*RX+RY*RY
    1940. 

  1940. ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
    1941. 

  1941.                 DPRIM=SQRT(RSQ)
    1942. 

  1942.                 D=DPRIM+RINDX*DCON*ABS(psurf(n)-.001*pbase(i,1))
    1943. 

  1943.                 DSQ=D*D
    1944. 

  1944.                 WTIJ(i)=(RIS-DSQ)/(RIS+DSQ)
    1945. 

  1945.                 WTIJ(i)=AMAX1(0.0,WTIJ(i))
    1946. 

  1946.               ENDDO
    1947. 

  1947.             else                                 ! wrf scheme 
    1948. 

  1948.               DO I=max0(its,MINI),min0(ite,MAXI)
    1949. 

  1949.                 RI=FLOAT(I)
    1950. 

  1950.                 RY=RI-RA(N)
    1951. 

  1951.                 RIS=RINDX*RINDX*sfcfacr*sfcfacr
    1952. 

  1952.                 RSQ=RX*RX+RY*RY
    1953. 

  1953. ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
    1954. 

  1954.                 wtij(i)=(ris-rsq)/(ris+rsq)
    1955. 

  1955.                 scratch = (abs (psurf(n)-.001*pbase(i,1))*DCON)
    1956. 

  1956.                 pdfac=1.-AMIN1(1.0,scratch)
    1957. 

  1957.                 wtij(i)=wtij(i)*pdfac
    1958. 

  1958.                 WTIJ(i)=AMAX1(0.0,WTIJ(i))
    1959. 

  1959.               ENDDO
    1960. 

  1960.             endif
    1961. 

  1961. 

  1962. ! Apply selected vertical spreading scheme.
    1963. 

  1963.             if(SFC_SCHEME_VERT.eq.1) then         ! original simple scheme 
    1964. 

  1964. 

  1965.               DO I=max0(its,MINI),min0(ite,MAXI)
    1966. 

  1966. 

  1967. ! try making sfc obs weighting go thru pbl
    1968. 

  1968.                 komax=max0(3,kpblt(i))
    1969. 

  1969.                 IF (iprt) THEN
    1970. 

  1970.                   if (kpblt(i).gt.25 .and. ktau.ne.0)                         &
    1971. 

  1971.                                      write(6,552)inest,i,j,kpblt(i)
    1972. 

  1972. 552               FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
    1973. 

  1973.                 ENDIF
    1974. 

  1974. 

  1975.                 if(kpblt(i).gt.25) komax=3
    1976. 

  1976.                 komin=1
    1977. 

  1977.                 dk=float(komax)
    1978. 

  1978. 

  1979.                 do k=komin,komax
    1980. 

  1980. 

  1981.                   wtsig(k)=float(komax-k+1)/dk
    1982. 

  1982.                   WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
    1983. 

  1983. 

  1984.                   WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K)    &
    1985. 

  1985.                               *WTSIG(K)*ERFIVR
    1986. 

  1986.                 enddo
    1987. 

  1987.               ENDDO
    1988. 

  1988. 

  1989.             else                                ! regime-based vif scheme
    1990. 

  1990. 

  1991. ! Here we calculate weights in the vertical coordinate, based on vih1 and vih2.
    1992. 

  1992. ! In the equation for wtsig(k), Z=zslab(i,k)-terrh(i) contains the model Z-values
    1993. 

  1993. ! (height above ground in meters) on a J-slab. The equation produces wtsig = 1.0 at
    1994. 

  1994. ! levels 1 to K, where z(K) < vih1 < z(K+1). For the example below, in the equation
    1995. 

  1995. ! for wtsig(k), the expression vih1(i)-Z(i,k) is strictly positive for k=1,2,3 since
    1996. 

  1996. ! levels 1, 2, and 3 are below vih1. So xtsig(k)=min(1.0, 1.0-x) where x > 0 ==>
    1997. 

  1997. ! wtsig(k)=1 for k=1,2,3.
    1998. 

  1998. !
    1999. 

  1999. !    For levels K+1 and up, wtsig will decrease linearly with height, with values
    2000. 

  2000. ! along the ramp that has value 1.0 at vih1 and 0.0 at vih2. In the example:
    2001. 

  2001. !
    2002. 

  2002. !   dz_ramp  = 1/(200-150) = 1/50
    2003. 

  2003. !   xtsig(4) = 1 + (150-175)/50 = 1 - 1/2 = 1/2
    2004. 

  2004. !
    2005. 

  2005. !       WTSIG
    2006. 

  2006. !         1 -|*    *     *     *     *     *
    2007. 

  2007. !            |
    2008. 

  2008. !            |                                *
    2009. 

  2009. !            |
    2010. 

  2010. !            |                                   *
    2011. 

  2011. !            |
    2012. 

  2012. !            |                                      *
    2013. 

  2013. !         0 -|--|-------|-----------|------|----|----|---------|---->  Z = HT ABOVE
    2014. 

  2014. !              15      55          115    150  175  200       250          GROUND
    2015. 

  2015. !              k=1     k=2         k=3    vih1 k=4  vih2      k=5
    2016. 

  2016.  
    2017. 

  2017.               DO I=max0(its,MINI),min0(ite,MAXI)
    2018. 

  2018. 

  2019.                 dz_ramp = 1.0 / max( 1.0, vih2(i)-vih1(i) )   ! vih2 >= vih1 by construct
    2020. 

  2020. 

  2021.            LML: do k = kts, kte 
    2022. 

  2022.                   wtsig(k) = min( 1.0, 1.0 + ( vih1(i)-zslab(i,k)+terrh(i) ) * dz_ramp )
    2023. 

  2023.                   wtsig(k) = max( 0.0, wtsig(k))
    2024. 

  2024. 

  2025.                   if(wtsig(k).le.0.0) EXIT LML
    2026. 

  2026.                     WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
    2027. 

  2027.                     WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K)    &
    2028. 

  2028.                                 *WTSIG(K)*ERFIVR
    2029. 

  2029.                   enddo LML
    2030. 

  2030.               ENDDO
    2031. 

  2031.             endif
    2032. 

  2032. 

  2033.           endif ! end if nudge-pbl switch is on
    2034. 

  2034. 

  2035.         endif   ! end check for obs in domain
    2036. 

  2036. ! END SURFACE-LAYER OBS NUDGING
    2037. 

  2037. 

  2038.       ELSE    
    2039. 

  2039. 

  2040. ! Compute temporal weight
    2041. 

  2041.         timewt = get_timewt(xtime,dtmin,twindo,1.,timeob(n))
    2042. 

  2042. 

  2043. 

  2044. ! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
    2045. 

  2045. !
    2046. 

  2046. !       print *,'in upper air section'
    2047. 

  2047. ! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
    2048. 

  2048. ! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
    2049. 

  2049. ! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
    2050. 

  2050. ! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
    2051. 

  2051. !ajb          SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
    2052. 

  2052. 

  2053.         slope = (RINFMN-RINFMX)/(psurf(n)-PFREE)
    2054. 

  2054. 

  2055.         RINFAC=SLOPE*POB+RINFMX-SLOPE*pfree
    2056. 

  2056.         RINFAC=AMAX1(RINFAC,RINFMN)
    2057. 

  2057.         RINFAC=AMIN1(RINFAC,RINFMX)
    2058. 

  2058. !yliu: for multilevel upper-air data, take the maximum
    2059. 

  2059. !      for the I loop.
    2060. 

  2060.         if(nsndlev.gt.1) RINFAC = RINFMX 
    2061. 

  2061. !yliu end
    2062. 

  2062. 

  2063.         MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
    2064. 

  2064.         MAXI=MIN0(ide-IGRID,MAXI)
    2065. 

  2065.         MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
    2066. 

  2066.         MINI=MAX0(1,MINI)
    2067. 

  2067. 

  2068. ! yliu start
    2069. 

  2069. ! use also obs outside of but close to this domain  -- upr data   
    2070. 

  2070. !     if(   RA(N).LT.(0.-RINFAC*RINDX)
    2071. 

  2071. !    & .or. RA(N).GT.float(IL)+RINFAC*RINDX
    2072. 

  2072. !    & .or. RB(N).LT.(0.-RINFAC*RINDX)
    2073. 

  2073. !    & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then          
    2074. 

  2074. !        print *, " skipped obs far away from this I-range"
    2075. 

  2075. ! currently can use obs within this domain or ones very close to (1/3
    2076. 

  2076. !   influence of radius in the coarse domain) this 
    2077. 

  2077. !   domain. In later case, use BC column value to approximate the model value 
    2078. 

  2078. !   at obs point -- ERRF need model field in errob.F !!
    2079. 

  2079. 

  2080. !cc         if (i.eq.39 .and. j.eq.34) then
    2081. 

  2081. !cc            write(6,*) 'RA(N) = ',ra(n) 
    2082. 

  2082. !cc            write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
    2083. 

  2083. !cc         endif
    2084. 

  2084.         if(   RA(N).GE.(0.-RINFAC*RINDX/3)                      &
    2085. 

  2085.         .and. RA(N).LE.float(ide)+RINFAC*RINDX/3                &
    2086. 

  2086.         .and. RB(N).GE.(0.-RINFAC*RINDX/3)                      &
    2087. 

  2087.         .and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
    2088. 

  2088. ! or use obs within this domain only
    2089. 

  2089. !     if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
    2090. 

  2090. !    &   RB(N).LT.1 .or. RB(N).GT.float(JL)) then
    2091. 

  2091. !        print *, " skipped obs far outside of this domain"
    2092. 

  2092. 

  2093. ! yliu end
    2094. 

  2094. ! is this 2 needed here - kpbl not used?
    2095. 

  2095. !          MINI=MAX0(2,MINI)
    2096. 

  2096. 

  2097. ! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
    2098. 

  2098. ! OBSERVATION N.  COMPUTE THE HORIZONTAL DISTANCE TO
    2099. 

  2099. ! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
    2100. 

  2100.           RJ=FLOAT(J)
    2101. 

  2101.           RX=RJ-RB(N)
    2102. 

  2102. ! WEIGHTS FOR THE 3-D VARIABLES
    2103. 

  2103. !
    2104. 

  2104.           ERFIVR=ERRF(IVAR,N)
    2105. 

  2105. ! jc
    2106. 

  2106.           nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
    2107. 

  2107. ! yliu start
    2108. 

  2108. ! test: do the sounding levels as individual obs
    2109. 

  2109. !        nsndlev=1
    2110. 

  2110. ! yliu end
    2111. 

  2111.           njcsnd=nsndlev
    2112. 

  2112. !
    2113. 

  2113.           DO I=max0(its,MINI),min0(ite,MAXI)
    2114. 

  2114. ! jc
    2115. 

  2115.             RI=FLOAT(I)
    2116. 

  2116.             RY=RI-RA(N)
    2117. 

  2117.             RIS=RINDX*RINFAC*RINDX*RINFAC
    2118. 

  2118.             RSQ=RX*RX+RY*RY
    2119. 

  2119. ! yliu test: for upper-air data, keep D1 influence radii
    2120. 

  2120.             WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
    2121. 

  2121. 

  2122.             WTIJ=AMAX1(0.0,WTIJ(i))
    2123. 

  2123. ! weight ob in vertical with +- 50 mb
    2124. 

  2124. ! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
    2125. 

  2125.             if(nsndlev.eq.1) then
    2126. 

  2126.               rinprs=7.5
    2127. 

  2127. ! 
    2128. 

  2128.             else
    2129. 

  2129.              rinprs=3.0
    2130. 

  2130.             endif
    2131. 

  2131. ! yliu end
    2132. 

  2132. 

  2133. !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    2134. 

  2134. !   ---  HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY  ---
    2135. 

  2135. !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    2136. 

  2136. 

  2137.             if(nsndlev.eq.1)then 
    2138. 

  2138. !----------------------------------------------------------------------
    2139. 

  2139. !         ---   HANDLE 1-LEVEL OBSERVATIONS  ---
    2140. 

  2140. !----------------------------------------------------------------------
    2141. 

  2141. 

  2142. !         if(I.eq.MINI) print *, "  Single snd "
    2143. 

  2143. ! ERFIVR is the residual (difference) between the ob and the model
    2144. 

  2144. ! at that point. We can analyze that residual up and down.
    2145. 

  2145. ! First find komin for ob.
    2146. 

  2146. !yliu start -- in the old code, komax and komin were reversed!
    2147. 

  2147.               do k=kte,1,-1
    2148. 

  2148.                 pijk = .001*(pbase(i,k)+pp(i,k))
    2149. 

  2149. !               print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
    2150. 

  2150.                 if(pijk.ge.(pob+rinprs)) then
    2151. 

  2151.                   komin=k
    2152. 

  2152.                   go to 325
    2153. 

  2153.                 endif
    2154. 

  2154.               enddo
    2155. 

  2155.               komin=1
    2156. 

  2156.  325          continue
    2157. 

  2157. ! now find komax for ob
    2158. 

  2158.               do k=3,kte
    2159. 

  2159.                 pijk = .001*(pbase(i,k)+pp(i,k))
    2160. 

  2160.                 if(pijk.le.(pob-rinprs)) then
    2161. 

  2161.                   komax=k
    2162. 

  2162.                   go to 326
    2163. 

  2163.                 endif
    2164. 

  2164.               enddo
    2165. 

  2165.               komax=kte   ! yliu 20050706
    2166. 

  2166.  326          continue
    2167. 

  2167. 

  2168. ! yliu: single-level upper-air data will act either above or below the PBL top
    2169. 

  2169. ! ajb: Reset komin or komax. if kobs>kpbl_obs, komin=kpbl_obs+1, else komax=kpbl_obs  
    2170. 

  2170. 

  2171.               if( (kpblt(i).le.komax) .and. (kpblt(i).ge.komin) ) then
    2172. 

  2172.                  kobs = 1
    2173. 

  2173.                  OBS_K: do k = komin, komax
    2174. 

  2174.                      if( pob .gt. .001*(pbase(i,k)+pp(i,k)) ) then
    2175. 

  2175.                         kobs = k
    2176. 

  2176.                         EXIT OBS_K
    2177. 

  2177.                      endif
    2178. 

  2178.                  enddo OBS_K
    2179. 

  2179. 

  2180.                  if(kobs.gt.kpbl_obs(n)) then
    2181. 

  2181. ! Obs will act only above the PBL top
    2182. 

  2182.                      komin=max0(kobs, komin)   ! kobs here is kpblt(i)+1
    2183. 

  2183.                  else                          ! Obs acts below PBL top
    2184. 

  2184. ! Obs will act only below the PBL top
    2185. 

  2185.                      komax=min0(kpblt(i), komax)
    2186. 

  2186.                  endif
    2187. 

  2187.               endif
    2188. 

  2188. ! yliu end
    2189. 

  2189. !
    2190. 

  2190. !           print *,'1 level, komin,komax=',komin,komax
    2191. 

  2191. !           if(i.eq.MINI) then
    2192. 

  2192. !             print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
    2193. 

  2193. !             ERFIVR=0
    2194. 

  2194. !           endif
    2195. 

  2195.               do k=1,kte
    2196. 

  2196.                 reserf(k)=0.0
    2197. 

  2197.                 wtsig(k)=0.0
    2198. 

  2198.               enddo
    2199. 

  2199. !yliu end
    2200. 

  2200. 

  2201. !ajb Add weights to sum only if nudge_pbl switch is on OR obs is above pbl top.
    2202. 

  2202.               if(nudge_pbl .or. komin.ge.kpblt(i)) then
    2203. 

  2203.                 do k=komin,komax
    2204. 

  2204.                   pijk = .001*(pbase(i,k)+pp(i,k))
    2205. 

  2205.                   reserf(k)=erfivr
    2206. 

  2206.                   wtsig(k)=1.-abs(pijk-pob)/rinprs
    2207. 

  2207.                   wtsig(k)=amax1(wtsig(k),0.0)
    2208. 

  2208. !             print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
    2209. 

  2209. ! Now calculate WT and WT2ERR for each i,j,k point                      cajb
    2210. 

  2210.                   WT(I,K)=WT(I,K)+TIMEWT*WTIJ(i)*wtsig(k)
    2211. 

  2211. 

  2212.                   WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*        &
    2213. 

  2213.                               reserf(k)*wtsig(k)*wtsig(k)
    2214. 

  2214.                 enddo
    2215. 

  2215.               endif
    2216. 

  2216. 

  2217.             else
    2218. 

  2218. 

  2219. !----------------------------------------------------------------------
    2220. 

  2220. !         ---   HANDLE MULTI-LEVEL OBSERVATIONS  ---
    2221. 

  2221. !----------------------------------------------------------------------
    2222. 

  2222. !yliu start 
    2223. 

  2223. !           if(I.eq.MINI) print *, "  Multi-level  snd "
    2224. 

  2224. !             print *, "   n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n)  &
    2225. 

  2225. !                  ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1) 
    2226. 

  2226.               if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
    2227. 

  2227.                 IF (iprt) THEN
    2228. 

  2228.                   write(msg,*) "n = ",n,"nsndlev = ",nsndlev 
    2229. 

  2229.                   call wrf_message(msg)
    2230. 

  2230.                   write(msg,*) "nlevs_ob,lev_in_ob",                          &
    2231. 

  2231.                            nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
    2232. 

  2232.                   call wrf_message(msg)
    2233. 

  2233.                   call wrf_message("in nudobs.F: sounding level messed up, stopping")
    2234. 

  2234.                 ENDIF
    2235. 

  2235.                 call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
    2236. 

  2236.              endif       
    2237. 

  2237. !yliu end
    2238. 

  2238. ! This is for a multi-level observation
    2239. 

  2239. ! The trick here is that the sounding is "one ob". You don't
    2240. 

  2240. !    want multiple levels to each be treated like separate
    2241. 

  2241. !    and independent observations.
    2242. 

  2242. ! At each i,j want to interpolate sounding to the model levels at that
    2243. 

  2243. !    particular point.
    2244. 

  2244.               komin=1
    2245. 

  2245.               komax=kte-2
    2246. 

  2246. ! this loop goes to 1501
    2247. 

  2247. ! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
    2248. 

  2248. !yliu start
    2249. 

  2249.               do k=1,kte
    2250. 

  2250.                 reserf(k)=0.0
    2251. 

  2251.                 wtsig(k)=0.0
    2252. 

  2252.               enddo
    2253. 

  2253. !yliu end
    2254. 

  2254. 

  2255.               do k=komax,komin,-1
    2256. 

  2256.   
    2257. 

  2257.                 pijk = .001*(pbase(i,k)+pp(i,k))
    2258. 

  2258. 

  2259. ! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
    2260. 

  2260.                 if(pijk.gt.varobs(5,n)) then
    2261. 

  2261.                   go to 1501
    2262. 

  2262.                 endif
    2263. 

  2263. 

  2264. ! if sigma level pressure is .lt. than the highest ob level, don't interpolate
    2265. 

  2265.                 if(pijk.le.varobs(5,n+nsndlev-1)) then 
    2266. 

  2266.                   go to 1501
    2267. 

  2267.                 endif
    2268. 

  2268. 

  2269. ! now interpolate sounding to this k
    2270. 

  2270. ! yliu start-- recalculate WTij for each k-level
    2271. 

  2271. !ajb      SLOPE = (RINFMN-RINFMX)/(pdoc(i,j)+PTOP-PFREE)        
    2272. 

  2272.                 slope = (RINFMN-RINFMX)/ (.001*pbase(i,1)-PFREE)
    2273. 

  2273.                 RINFAC=SLOPE*pijk+RINFMX-SLOPE*PFREE              
    2274. 

  2274.                 RINFAC=AMAX1(RINFAC,RINFMN)      
    2275. 

  2275.                 RINFAC=AMIN1(RINFAC,RINFMX)
    2276. 

  2276.                 RIS=RINDX*RINFAC*RINDX*RINFAC  
    2277. 

  2277.                 RSQ=RX*RX+RY*RY               
    2278. 

  2278. 

  2279. ! for upper-air data, keep D1 influence radii
    2280. 

  2280.                 WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)      
    2281. 

  2281.                 WTIJ(i)=AMAX1(0.0,WTIJ(i))
    2282. 

  2282. ! yliu end
    2283. 

  2283. 

  2284. ! this loop goes to 1503
    2285. 

  2285.                 do nn=2,nsndlev
    2286. 

  2286. 

  2287. ! only set pobhi if varobs(ivar) is ok
    2288. 

  2288.                   pobhi=-888888.
    2289. 

  2289. 

  2290.                   if(varobs(ivar,n+nn-1).gt.-800000.                           &
    2291. 

  2291.                   .and. varobs(5,n+nn-1).gt.-800000.) then
    2292. 

  2292.                     pobhi=varobs(5,n+nn-1)
    2293. 

  2293.                     nhi=n+nn-1
    2294. 

  2294.                     if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.0.5*maxsnd_gap) then
    2295. 

  2295.                       go to 1502        ! within maxsnd_gap/2 of obs height
    2296. 

  2296.                     endif
    2297. 

  2297.                   endif
    2298. 

  2298. 

  2299.                 enddo
    2300. 

  2300. 

  2301. ! did not find any ob above within maxsnd_gap/2, so jump out 
    2302. 

  2302.                 go to 1501
    2303. 

  2303.  1502           continue
    2304. 

  2304. 

  2305.                 nlo=nhi-1
    2306. 

  2306.                 do nnjc=nhi-1,n,-1 
    2307. 

  2307.                   if(varobs(ivar,nnjc).gt.-800000.                             &
    2308. 

  2308.                   .and. varobs(5,nnjc).gt.-800000.) then
    2309. 

  2309.                     poblo=varobs(5,nnjc)
    2310. 

  2310.                     nlo=nnjc
    2311. 

  2311.                     if(poblo.gt.pijk .and. abs(poblo-pijk).lt.0.5*maxsnd_gap) then
    2312. 

  2312.                       go to 1505        ! within maxsnd_gap/2 of obs height
    2313. 

  2313.                     endif
    2314. 

  2314.                   endif
    2315. 

  2315.                 enddo
    2316. 

  2316. !yliu end --
    2317. 

  2317. 

  2318. ! did not find any ob below within maxsnd_gap/2, so jump out 
    2319. 

  2319.                 go to 1501
    2320. 

  2320.  1505           continue
    2321. 

  2321. 

  2322. ! interpolate to model level
    2323. 

  2323.                 pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
    2324. 

  2324.                 reserf(k)=errf(ivar,nlo)+                               &
    2325. 

  2325.                             (errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
    2326. 

  2326.                 wtsig(k)=1.
    2327. 

  2327.   
    2328. 

  2328.  1501           continue
    2329. 

  2329. 

  2330. ! now calculate WT and WT2ERR for each i,j,k point                               cajb
    2331. 

  2331. !ajb Add weights to sum only if nudge_pbl switch is on OR k > kpblt.
    2332. 

  2332.                 if(nudge_pbl .or. k.gt.kpblt(i)) then
    2333. 

  2333. 

  2334.                   WT(I,K)=WT(I,K)+TIMEWT*WTIJ(i)*wtsig(k)
    2335. 

  2335.   
    2336. 

  2336.                   WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*        &
    2337. 

  2337.                               reserf(k)*wtsig(k)*wtsig(k)
    2338. 

  2338.                 endif
    2339. 

  2339. 

  2340.               enddo   ! enddo k levels
    2341. 

  2341. 

  2342. ! end multi-levels
    2343. 

  2343.             endif  ! end if(nsndlev.eq.1)
    2344. 

  2344. !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    2345. 

  2345. !   END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
    2346. 

  2346. !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    2347. 

  2347. !
    2348. 

  2348.           ENDDO ! END DO MINI,MAXI LOOP
    2349. 

  2349. 

  2350.         endif ! check for obs in domain
    2351. 

  2351. 

  2352. ! END OF NUDGING TO OBS ON PRESSURE LEVELS
    2353. 

  2353. 

  2354.       ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
    2355. 

  2355. 

  2356. !----------------------------------------------------------------------
    2357. 

  2357.     ENDIF  ! END SECTION FOR PROCESSING OF OBSERVATION
    2358. 

  2358. !----------------------------------------------------------------------
    2359. 

  2359. 

  2360. !   n=n+1
    2361. 

  2361.     n=n+njcsnd
    2362. 

  2362. 

  2363. !yliu  1202 continue
    2364. 

  2364.     if(n.gt.nstat)then
    2365. 

  2365. !     print *,'n,nstat=',n,nstat,ivar,j
    2366. 

  2366.       go to 1203
    2367. 

  2367.     endif
    2368. 

  2368. !   print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n) 
    2369. 

  2369. 

  2370. !***********************************************************************
    2371. 

  2371.   ENDDO  ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
    2372. 

  2372. !***********************************************************************
    2373. 

  2373. 

  2374.  1203 continue
    2375. 

  2375. 

  2376. ! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED.  NOW
    2377. 

  2377. ! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
    2378. 

  2378. ! THE ATEN ARRAY
    2379. 

  2379. ! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
    2380. 

  2380. ! THEY ARE USED BELOW IN THE DENOMINATOR.
    2381. 

  2381.   DO K=kts,kte
    2382. 

  2382.     DO I=its,ite
    2383. 

  2383.       IF(WT(I,K).EQ.0)THEN
    2384. 

  2384.         WT2ERR(I,K)=0.0
    2385. 

  2385.       ENDIF
    2386. 

  2386.       IF(WT(I,K).EQ.0)THEN
    2387. 

  2387.         WT(I,K)=1.0
    2388. 

  2388.       ENDIF
    2389. 

  2389.     ENDDO
    2390. 

  2390.   ENDDO
    2391. 

  2391. 

  2392. 126 CONTINUE
    2393. 

  2393. 

  2394.   IF(IVAR.GE.3)GOTO 170
    2395. 

  2395. ! this is for u,v
    2396. 

  2396. ! 3-D DOT POINT TENDENCIES
    2397. 

  2397.  
    2398. 

  2398. ! Calculate scales for converting nudge factor from u (v)
    2399. 

  2399. ! to rho_u (or rho_v) units.
    2400. 

  2400. 

  2401.   IF (IVAR == 1) THEN
    2402. 

  2402.      call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
    2403. 

  2403.   ELSE IF (IVAR == 2) THEN
    2404. 

  2404.      call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
    2405. 

  2405.   END IF
    2406. 

  2406.  
    2407. 

  2407.   DO K=1,kte
    2408. 

  2408. 

  2409.     DO I=i_s,i_e
    2410. 

  2410. 

  2411.       IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
    2412. 

  2412.         W2EOWT=WT2ERR(I,K)/WT(I,K)
    2413. 

  2413.       ELSE
    2414. 

  2414.         W2EOWT=SAVWT(IPL,I,K)
    2415. 

  2415. 

  2416. !        write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,W2EOWT = ',i,j,k,ipl,W2EOWT
    2417. 

  2417. 

  2418.       ENDIF
    2419. 

  2419. 

  2420. !      if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
    2421. 

  2421. !           scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
    2422. 

  2422. !           write(6,*) 'ATEN calc: k = ',k
    2423. 

  2423. !           write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
    2424. 

  2424. !           write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
    2425. 

  2425. !     $               ' W2EOWT = ',w2eowt
    2426. 

  2426. !           write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
    2427. 

  2427. !     $               ' GFACTOR = ',gfactor
    2428. 

  2428. !       endif
    2429. 

  2429. !
    2430. 

  2430. !      if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
    2431. 

  2431. !           scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
    2432. 

  2432. !           write(6,*) 'ATEN calc: k = ',k
    2433. 

  2433. !           write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
    2434. 

  2434. !           write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
    2435. 

  2435. !     $                ' W2EOWT = ',w2eowt
    2436. 

  2436. !           write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
    2437. 

  2437. !     $                ' GFACTOR = ',gfactor
    2438. 

  2438. !       endif
    2439. 

  2439. 

  2440. !      if(ivar .eq. 1 .and. i.eq.38 .and. (j.ge.36.and.j.le.62) .and. k.eq.7) then
    2441. 

  2441. !           scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
    2442. 

  2442. !           write(6,*)
    2443. 

  2443. !           write(6,*) 'ATEN calc: j = ',j,' k = ',k
    2444. 

  2444. !           write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
    2445. 

  2445. !           write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),' W2EOWT = ',w2eowt
    2446. 

  2446. !           write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,' GFACTOR = ',gfactor
    2447. 

  2447. !       endif
    2448. 

  2448. 

  2449.         ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I)                        &
    2450. 

  2450.                     *W2EOWT*TFACI                           &
    2451. 

  2451.                     *ISWIND       *GFACTOR   !yliu *GFACTOR 
    2452. 

  2452. 

  2453. !      if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
    2454. 

  2454. !           write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
    2455. 

  2455. !      endif
    2456. 

  2456. !      if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
    2457. 

  2457. !           write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
    2458. 

  2458. !      endif
    2459. 

  2459. 

  2460.     ENDDO
    2461. 

  2461.   ENDDO
    2462. 

  2462. 

  2463.   IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
    2464. 

  2464.     DO K=1,kte
    2465. 

  2465.       DO I=its,ite
    2466. 

  2466.         SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
    2467. 

  2467. 

  2468. !        write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,savwt = ',i,j,k,ipl,savwt(ipl,i,k)
    2469. 

  2469.       ENDDO
    2470. 

  2470.     ENDDO
    2471. 

  2471.   ENDIF
    2472. 

  2472. 

  2473.   RETURN
    2474. 

  2474. 

  2475. 170 CONTINUE
    2476. 

  2476. 

  2477. ! 3-D CROSS-POINT TENDENCIES
    2478. 

  2478. ! this is for t (ivar=3) and q (ivsr=4)
    2479. 

  2479.   IF(3-IVAR.LT.0)THEN
    2480. 

  2480.     GITQ=GIQ
    2481. 

  2481.   ELSE
    2482. 

  2482.     GITQ=GIT
    2483. 

  2483.   ENDIF
    2484. 

  2484.   IF(3-IVAR.LT.0)THEN
    2485. 

  2485.     ISTQ=ISMOIS
    2486. 

  2486.   ELSE
    2487. 

  2487.     ISTQ=ISTEMP
    2488. 

  2488.   ENDIF
    2489. 

  2489. 

  2490.   DO K=1,kte
    2491. 

  2491.     DO I=i_s,i_e
    2492. 

  2492.       IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
    2493. 

  2493.         W2EOWT=WT2ERR(I,K)/WT(I,K)
    2494. 

  2494.       ELSE
    2495. 

  2495.         W2EOWT=SAVWT(IPL,I,K)
    2496. 

  2496.       ENDIF
    2497. 

  2497. 

  2498. !        if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
    2499. 

  2499. !            scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
    2500. 

  2500. !            write(6,*) 'ATEN calc: k = ',k
    2501. 

  2501. !            write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
    2502. 

  2502. !            write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),                  &
    2503. 

  2503. !                       ' W2EOWT = ',w2eowt
    2504. 

  2504. !            write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq,              &
    2505. 

  2505. !                       ' GFACTOR = ',gfactor
    2506. 

  2506. !        endif
    2507. 

  2507.  
    2508. 

  2508. !        if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
    2509. 

  2509. !            scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
    2510. 

  2510. !            write(6,*) 'ATEN calc: k = ',k
    2511. 

  2511. !            write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
    2512. 

  2512. !            write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
    2513. 

  2513. !     $                 ' W2EOWT = ',w2eowt
    2514. 

  2514. !            write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq,
    2515. 

  2515. !     $                 ' GFACTOR = ',gfactor
    2516. 

  2516. !        endif
    2517. 

  2517. 

  2518.       ATEN(i,k)=ATEN(i,k)+GITQ*MU(I)                       &
    2519. 

  2519.                   *W2EOWT*TFACI*ISTQ       *GFACTOR   !yliu *GFACTOR
    2520. 

  2520. 

  2521. !        if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
    2522. 

  2522. !            write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
    2523. 

  2523. !        endif
    2524. 

  2524. !        if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
    2525. 

  2525. !            write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
    2526. 

  2526. !        endif
    2527. 

  2527. 

  2528.     ENDDO
    2529. 

  2529.   ENDDO
    2530. 

  2530. 

  2531.   IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
    2532. 

  2532.     DO K=1,kte
    2533. 

  2533.       DO I=its,ite
    2534. 

  2534.         SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
    2535. 

  2535.       ENDDO
    2536. 

  2536.     ENDDO
    2537. 

  2537.   ENDIF
    2538. 

  2538. 

  2539.   RETURN
    2540. 

  2540.   END SUBROUTINE nudob
    2541. 

  2541. 

  2542.   SUBROUTINE date_string(year, month, day, hour, minute, second, cdate)
    2543. 

  2543. !-----------------------------------------------------------------------
    2544. 

  2544. !  PURPOSE: Form a date string (YYYY-MM-DD_hh:mm:ss) from integer
    2545. 

  2545. !           components.
    2546. 

  2546. !-----------------------------------------------------------------------
    2547. 

  2547.   IMPLICIT NONE
    2548. 

  2548. !-----------------------------------------------------------------------
    2549. 

  2549. 

  2550.   INTEGER, INTENT(IN)  :: year
    2551. 

  2551.   INTEGER, INTENT(IN)  :: month
    2552. 

  2552.   INTEGER, INTENT(IN)  :: day
    2553. 

  2553.   INTEGER, INTENT(IN)  :: hour
    2554. 

  2554.   INTEGER, INTENT(IN)  :: minute
    2555. 

  2555.   INTEGER, INTENT(IN)  :: second
    2556. 

  2556.   CHARACTER*19, INTENT(INOUT) :: cdate
    2557. 

  2557. 

  2558. ! Local variables
    2559. 

  2559.   integer   :: ic                    ! loop counter
    2560. 

  2560. 

  2561.       cdate(1:19)  = "0000-00-00_00:00:00"
    2562. 

  2562.       write(cdate( 1: 4),'(i4)') year
    2563. 

  2563.       write(cdate( 6: 7),'(i2)') month
    2564. 

  2564.       write(cdate( 9:10),'(i2)') day
    2565. 

  2565.       write(cdate(12:13),'(i2)') hour
    2566. 

  2566.       write(cdate(15:16),'(i2)') minute
    2567. 

  2567.       write(cdate(18:19),'(i2)') second
    2568. 

  2568.       do ic = 1,19
    2569. 

  2569.         if(cdate(ic:ic) .eq. " ") cdate(ic:ic) = "0"
    2570. 

  2570.       enddo
    2571. 

  2571. 

  2572.   RETURN
    2573. 

  2573.   END SUBROUTINE date_string
    2574. 

  2574. 

  2575.   SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
    2576. 

  2576. !-----------------------------------------------------------------------
    2577. 

  2577.   IMPLICIT NONE
    2578. 

  2578. !-----------------------------------------------------------------------
    2579. 

  2579. 

  2580.   INTEGER, INTENT(IN)  :: ims,ime           ! Memory dimensions
    2581. 

  2581.   INTEGER, INTENT(IN)  :: its,ite           ! Tile   dimensions
    2582. 

  2582.   REAL, INTENT(IN)     :: a( ims:ime )      ! Air mass array 
    2583. 

  2583.   REAL, INTENT(IN)     :: msf( ims:ime )    ! Map scale factor array
    2584. 

  2584.   REAL, INTENT(OUT)    :: rscale( ims:ime ) ! Scales for rho-coupling
    2585. 

  2585. 

  2586. ! Local variables
    2587. 

  2587.   integer :: i
    2588. 

  2588. 

  2589. ! Calculate scales to be used for producing rho-coupled nudging factors.
    2590. 

  2590.   do i = its,ite
    2591. 

  2591.     rscale(i) = a(i)/msf(i)
    2592. 

  2592.   enddo
    2593. 

  2593. 

  2594.   RETURN
    2595. 

  2595.   END SUBROUTINE calc_rcouple_scales
    2596. 

  2596. 

  2597.   SUBROUTINE print_obs_info(iprt,inest,niobf,rio,rjo,rko,                &
    2598. 

  2598.                             prt_max,prt_freq,obs,stnid,lat,lon,          &
    2599. 

  2599.                             mlat,mlon,timeob,xtime)
    2600. 

  2600. !*************************************************************************
    2601. 

  2601. ! Purpose: Print obs information.
    2602. 

  2602. !*************************************************************************
    2603. 

  2603. 

  2604.   IMPLICIT NONE
    2605. 

  2605. 

  2606.   LOGICAL, intent(in)    :: iprt          ! Print flag
    2607. 

  2607.   INTEGER, intent(in)    :: inest         ! Nest level
    2608. 

  2608.   INTEGER, intent(in)    :: niobf         ! Maximum number of observations
    2609. 

  2609.   REAL,    intent(in)    :: rio(niobf)    ! West-east coord (non-stagger)
    2610. 

  2610.   REAL,    intent(in)    :: rjo(niobf)    ! South-north coord (non-stagger)
    2611. 

  2611.   REAL,    intent(in)    :: rko(niobf)    ! Bottom-top north coord (non-stagger)
    2612. 

  2612.   INTEGER, intent(in)    :: prt_max        ! Max no. of obs for diagnostic printout
    2613. 

  2613.   INTEGER, intent(in)    :: prt_freq       ! Frequency for diagnostic printout
    2614. 

  2614.   INTEGER, intent(in)    :: obs(prt_max)  ! Saved obs indices to print
    2615. 

  2615.   INTEGER, intent(in)    :: stnid(40,prt_max) ! Saved station ids
    2616. 

  2616.   REAL,    intent(in)    :: lat(prt_max)  ! Saved latitudes
    2617. 

  2617.   REAL,    intent(in)    :: lon(prt_max)  ! Saved longitudes
    2618. 

  2618.   REAL,    intent(in)    :: mlat(prt_max) ! Saved model latitudes
    2619. 

  2619.   REAL,    intent(in)    :: mlon(prt_max) ! Saved longitudes
    2620. 

  2620.   REAL,    intent(in)    :: timeob(niobf) ! Times of each observation (hours)
    2621. 

  2621.   REAL,    intent(in)    :: xtime         ! Model time in minutes
    2622. 

  2622. 

  2623. ! Local variables
    2624. 

  2624.   integer :: i                    ! Loop counter over obs station chars
    2625. 

  2625.   integer :: n                    ! Loop counter over obs
    2626. 

  2626.   integer :: pnx                  ! Obs index for printout
    2627. 

  2627.   character(len=200) :: msg       ! Argument to wrf_message
    2628. 

  2628.   character(len=20)  :: station_id ! Station id of observation
    2629. 

  2629. 

  2630.   if(iprt) then
    2631. 

  2631.     if(prt_max.gt.0) then
    2632. 

  2632. 

  2633.       if(obs(1).ne.-999) then
    2634. 

  2634. 

  2635.         call wrf_message("")
    2636. 

  2636.         write(msg,fmt='(a,i4,a,f8.1,a)') 'REPORTING OBS MASS-PT LOCS FOR NEST ',  &
    2637. 

  2637.                                      inest,' AT XTIME=',xtime,' MINUTES'
    2638. 

  2638.         call wrf_message(msg)
    2639. 

  2639. 

  2640.         write(msg,fmt='(a,i4,a,i5,a)') 'FREQ=',prt_freq,', MAX=',prt_max,         &
    2641. 

  2641.                            ' LOCS, NEWLY READ OBS ONLY, -999 => OBS OFF PROC'
    2642. 

  2642.         call wrf_message(msg)
    2643. 

  2643.         call wrf_message("")
    2644. 

  2644. 

  2645.         write(msg,fmt='(3a)') '    OBS#     I       J       K     OBS LAT',       &
    2646. 

  2646.                           '  OBS LON   XLAT(I,J)  XLONG(I,J)  TIME(hrs)',     &
    2647. 

  2647.                           '  OBS STATION ID'
    2648. 

  2648.         call wrf_message(msg)
    2649. 

  2649. 

  2650.       endif
    2651. 

  2651.     endif
    2652. 

  2652. 

  2653. ! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
    2654. 

  2654. !       Hence subtract .5 from each to get mass-point coords.
    2655. 

  2655.     do n=1,prt_max
    2656. 

  2656.        pnx = obs(n)
    2657. 

  2657.        if(pnx.ne.-999) then
    2658. 

  2658. !          Retrieve 15 chars of station id
    2659. 

  2659.            do i = 1,15
    2660. 

  2660.               station_id(i:i) = char(stnid(i,n))
    2661. 

  2661.            enddo
    2662. 

  2662.            write(msg,fmt='(2x,i7,3f8.3,2f9.3,2x,f9.3,2x,f9.3,3x,f6.2,7x,a15)')    &
    2663. 

  2663.                pnx,rio(pnx)-.5,rjo(pnx)-.5,rko(pnx),lat(n),lon(n),            &
    2664. 

  2664.                mlat(n),mlon(n),timeob(pnx),station_id
    2665. 

  2665.         call wrf_message(msg)
    2666. 

  2666.        endif
    2667. 

  2667.     enddo
    2668. 

  2668.     if(obs(1).ne.-999) call wrf_message("")
    2669. 

  2669.   endif
    2670. 

  2670.   END SUBROUTINE print_obs_info
    2671. 

  2671. 

  2672.   REAL FUNCTION ht_to_p( h, pbbc, ppbc, z, ic, jc, dx, dy,                    &
    2673. 

  2673.                          k_start, k_end, kds,kde, ims,ime, jms,jme, kms,kme )
    2674. 

  2674. 

  2675. !******************************************************************************
    2676. 

  2676. ! Purpose: Interpolate pressure at a specified x (ic), y (jc), and height (h).
    2677. 

  2677. !          The input pressure column pbbc+ppbc (base and perturbn) must already
    2678. 

  2678. !          be horizontally interpolated to the x, y position. The subroutine
    2679. 

  2679. !          get_height_column is called here to horizontally interpolated the
    2680. 

  2680. !          3D height field z to get a height column at (iob, job).
    2681. 

  2681. !******************************************************************************
    2682. 

  2682. 

  2683.   IMPLICIT NONE
    2684. 

  2684. 

  2685.   REAL,    INTENT(IN)  :: h                                ! height value (m)
    2686. 

  2686.   INTEGER, INTENT(IN)  :: k_start, k_end                   ! loop bounds  
    2687. 

  2687.   INTEGER, INTENT(IN)  :: kds,kde                          ! vertical dim.
    2688. 

  2688.   INTEGER, INTENT(IN)  :: ims,ime, jms,jme, kms,kme        ! memory dims.
    2689. 

  2689.   REAL,    INTENT(IN)  :: pbbc(kds:kde)                    ! column base pressure (cb)
    2690. 

  2690.   REAL,    INTENT(IN)  :: ppbc(kds:kde)                    ! column pressure perturbation (cb)
    2691. 

  2691.   REAL,    INTENT(IN)  :: z( ims:ime, kms:kme, jms:jme )   ! ht (m) above sl on half-levels
    2692. 

  2692.   INTEGER, INTENT(IN)  :: ic                               ! i-coord of desired p
    2693. 

  2693.   INTEGER, INTENT(IN)  :: jc                               ! j-coord of desired p
    2694. 

  2694.   REAL,    INTENT(IN)  :: dx                               ! interp. fraction (x dir)
    2695. 

  2695.   REAL,    INTENT(IN)  :: dy                               ! interp. fraction (y dir)
    2696. 

  2696. 

  2697. ! Local variables
    2698. 

  2698.   INTEGER :: k               ! loop counter
    2699. 

  2699.   INTEGER :: klo             ! lower k bound
    2700. 

  2700.   REAL :: zlo                ! lower z bound for h
    2701. 

  2701.   REAL :: zhi                ! upper z bound for h
    2702. 

  2702.   REAL :: p                  ! interpolated pressure value
    2703. 

  2703.   REAL :: ln_p               ! log p
    2704. 

  2704.   REAL :: ln_plo             ! log plo
    2705. 

  2705.   REAL :: ln_phi             ! log phi
    2706. 

  2706.   REAL :: z_at_p( kms:kme )  ! height at p levels
    2707. 

  2707. 

  2708. ! Get interpolated z column on pressure (half-) levels at (ic,jc)
    2709. 

  2709.   call get_height_column(z, ic, jc, dx, dy, z_at_p,                   &
    2710. 

  2710.                          k_start, k_end, kds,kde,                     &
    2711. 

  2711.                          ims,ime, jms,jme, kms,kme )
    2712. 

  2712. 

  2713. ! Now we have pbbc, ppbc, z_at_p, so compute p at h. First, find
    2714. 

  2714. ! bounding layers klo and khi so that z_at_p(klo) <= h <= z_at_p(khi)
    2715. 

  2715. 

  2716.   ZLEVS: do k = k_start+1, k_end
    2717. 

  2717.     klo = k-1
    2718. 

  2718.     if(h .le. z_at_p(k)) then
    2719. 

  2719.       EXIT ZLEVS
    2720. 

  2720.     endif
    2721. 

  2721.   enddo ZLEVS
    2722. 

  2722. 

  2723.   zlo = z_at_p(klo)
    2724. 

  2724.   zhi = z_at_p(klo+1)
    2725. 

  2725. 

  2726. ! Interpolate natural log of pressure
    2727. 

  2727.   ln_plo = log( pbbc(klo+1) + ppbc(klo+1) )
    2728. 

  2728.   ln_phi = log( pbbc(klo) + ppbc(klo) )
    2729. 

  2729.   if(h.le.zlo) then
    2730. 

  2730.     ln_p = ln_phi     ! set to k=1 pressure 
    2731. 

  2731.   else if (h.ge.zhi) then
    2732. 

  2732.     ln_p = ln_plo     ! set to k=k_end pressure 
    2733. 

  2733.   else
    2734. 

  2734.     ln_p = ln_plo + (ln_phi-ln_plo)*((zhi-h)/(zhi-zlo)) 
    2735. 

  2735.   endif
    2736. 

  2736. 

  2737. ! Return pressure
    2738. 

  2738.   p = exp(ln_p)
    2739. 

  2739.   ht_to_p = p
    2740. 

  2740.   RETURN
    2741. 

  2741.   END FUNCTION ht_to_p
    2742. 

  2742. 

  2743.   SUBROUTINE get_height_column( z, ic, jc, dx, dy, z_at_p,                  &
    2744. 

  2744.                                 k_start, k_end, kds,kde,                    &
    2745. 

  2745.                                 ims,ime, jms,jme, kms,kme )
    2746. 

  2746. !*************************************************************************
    2747. 

  2747. ! Purpose: Compute column height at ic, jc location on p levels
    2748. 

  2748. !*************************************************************************
    2749. 

  2749. 

  2750.   IMPLICIT NONE
    2751. 

  2751. 

  2752.   INTEGER, INTENT(IN)  :: k_start, k_end                   ! Loop bounds  
    2753. 

  2753.   INTEGER, INTENT(IN)  :: kds,kde                          ! vertical dim.
    2754. 

  2754.   INTEGER, INTENT(IN)  :: ims,ime, jms,jme, kms,kme        ! memory dims.
    2755. 

  2755.   REAL,    INTENT(IN)  :: z( ims:ime, kms:kme, jms:jme )   ! ht (m) on half-levels
    2756. 

  2756.   INTEGER, INTENT(IN)  :: ic                               ! i-coord of desired p
    2757. 

  2757.   INTEGER, INTENT(IN)  :: jc                               ! j-coord of desired p
    2758. 

  2758.   REAL,    INTENT(IN)  :: dx                               ! interp. fraction (x dir)
    2759. 

  2759.   REAL,    INTENT(IN)  :: dy                               ! interp. fraction (y dir)
    2760. 

  2760.   REAL,    INTENT(OUT) :: z_at_p( kms:kme )                ! column height at p levels
    2761. 

  2761. 

  2762. ! Local variables
    2763. 

  2763.   INTEGER :: k             ! loop counter
    2764. 

  2764. 

  2765. 

  2766.   do k = kds, kde
    2767. 

  2767.       z_at_p(k) =                                     & 
    2768. 

  2768.          (1.-DY)*( (1.-DX)*z(IC,K,JC) +               &
    2769. 

  2769.                             DX *z(IC+1,K,JC) ) +      &
    2770. 

  2770.              DY* ( (1.-DX)*z(IC,K,JC+1) +             &
    2771. 

  2771.                             DX *z(IC+1,K,JC+1) )
    2772. 

  2772.   enddo
    2773. 

  2773. 

  2774.   END SUBROUTINE get_height_column
    2775. 

  2775. 

  2776.   SUBROUTINE get_base_state_height_column( p_top, p00, t00, a, g, r_d,    &
    2777. 

  2777.                                znu, z_at_p,  k_start, k_end, kds,kde, kms,kme )
    2778. 

  2778. !********************************************************
    2779. 

  2779. ! Purpose: Compute base-state column height on p levels.
    2780. 

  2780. !    See, e.g., eqns 1.3-1.5 in MM5 User's Guide.
    2781. 

  2781. !    Height is a function of pressure plus the constants:
    2782. 

  2782. !    p00  - sea level pressure (Pa)
    2783. 

  2783. !    t00  - sea level temperature (K)
    2784. 

  2784. !    a    - base state lapse rate (k/m)
    2785. 

  2785. !    r_d  - gas constant (J/Kg/K)   (Joule=1kg/m**2/s**2)
    2786. 

  2786. !    g    - gravitational constant (m/s**2)
    2787. 

  2787. !********************************************************
    2788. 

  2788. 

  2789.   IMPLICIT NONE
    2790. 

  2790. 

  2791.   REAL, INTENT(IN)     :: p_top        ! pressure at top of model
    2792. 

  2792.   REAL, INTENT(IN)     :: p00          ! base state pressure
    2793. 

  2793.   REAL, INTENT(IN)     :: t00          ! base state temperature
    2794. 

  2794.   REAL, INTENT(IN)     :: a            ! base state lapse rate
    2795. 

  2795.   REAL, INTENT(IN)     :: g                ! gravity constant
    2796. 

  2796.   REAL, INTENT(IN)     :: r_d              ! gas constant
    2797. 

  2797.   INTEGER, INTENT(IN)  :: k_start, k_end   ! Loop bounds  
    2798. 

  2798.   INTEGER, INTENT(IN)  :: kds,kde          ! vertical dim.
    2799. 

  2799.   INTEGER, INTENT(IN)  :: kms,kme          ! vertical memory dim.
    2800. 

  2800.   REAL, INTENT(IN)  :: znu( kms:kme )      ! eta values on half (mass) levels
    2801. 

  2801.   REAL, INTENT(OUT) :: z_at_p( kms:kme )   ! column height at p levels
    2802. 

  2802. 

  2803. ! Local variables
    2804. 

  2804.   integer :: k             ! loop counter
    2805. 

  2805.   real    :: ps0           ! base state pressure at surface
    2806. 

  2806.   real    :: pb(kms:kme)   ! pressure on half eta levels
    2807. 

  2807.   real    :: logterm       ! temporary for Z calculation 
    2808. 

  2808.   real    :: ginv          ! 1/g
    2809. 

  2809.   
    2810. 

  2810.   ginv = 1/g
    2811. 

  2811. 

  2812. ! Compute base state pressure on half eta levels.
    2813. 

  2813.    do k = k_start, k_end
    2814. 

  2814.      pb(k) = znu(k)*(p00 - p_top) + p_top
    2815. 

  2815.    enddo
    2816. 

  2816. 

  2817. ! Use hydrostatic relation to compute height at pressure levels.
    2818. 

  2818.    do k = k_start, k_end
    2819. 

  2819.      logterm = log(pb(k)/p00)
    2820. 

  2820.      z_at_p(k) = .5*r_d*a*ginv*logterm*logterm - r_d*t00*ginv*logterm
    2821. 

  2821.    enddo
    2822. 

  2822. 

  2823.   END SUBROUTINE get_base_state_height_column
    2824. 

  2824. 

  2825.   REAL FUNCTION get_timewt(xtime,dtmin,twindo,scalef,obtime)
    2826. 

  2826. !*************************************************************************
    2827. 

  2827. ! Purpose: Compute the temporal weight factor for an observation
    2828. 

  2828. !*************************************************************************
    2829. 

  2829. 

  2830.   IMPLICIT NONE
    2831. 

  2831. 

  2832.   REAL, INTENT(IN)  :: xtime              ! model time (minutes)
    2833. 

  2833.   REAL, INTENT(IN)  :: dtmin              ! model timestep (minutes)
    2834. 

  2834.   REAL, INTENT(IN)  :: twindo             ! half window (hours)
    2835. 

  2835.   REAL, INTENT(IN)  :: scalef             ! window scale factor
    2836. 

  2836.   REAL, INTENT(IN)  :: obtime             ! observation time (hours)
    2837. 

  2837. 

  2838. ! Local variables
    2839. 

  2839.   real :: fdtim            ! reference time (minutes)
    2840. 

  2840.   real :: tw1              ! half of twindo, scaled, in minutes
    2841. 

  2841.   real :: tw2              ! twindo, scaled, in minutes
    2842. 

  2842.   real :: tconst           ! reciprical of tw1
    2843. 

  2843.   real :: ttim             ! obtime in minutes
    2844. 

  2844.   real :: dift             ! | fdtim-ttim |
    2845. 

  2845.   real :: timewt           ! returned weight
    2846. 

  2846. 

  2847. ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
    2848. 

  2848.   FDTIM=XTIME-DTMIN
    2849. 

  2849. ! TWINDO IS IN MINUTES:
    2850. 

  2850.   TW1=TWINDO/2.*60.*scalef
    2851. 

  2851.   TW2=TWINDO*60.*scalef
    2852. 

  2852.   TCONST=1./TW1
    2853. 

  2853.   TIMEWT=0.0
    2854. 

  2854.   TTIM=obtime*60.
    2855. 

  2855. !***********TTIM=TARGET TIME IN MINUTES
    2856. 

  2856.   DIFT=ABS(FDTIM-TTIM)
    2857. 

  2857.   IF(DIFT.LE.TW1)TIMEWT=1.0
    2858. 

  2858.   IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
    2859. 

  2859.      IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
    2860. 

  2860.      IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
    2861. 

  2861.   ENDIF
    2862. 

  2862.   get_timewt = timewt
    2863. 

  2863.   END FUNCTION get_timewt
    2864. 

  2864. 

  2865.   SUBROUTINE print_vif_var(var, vif, nfullmin, nrampmin )
    2866. 

  2866. !********************************************************
    2867. 

  2867. ! Purpose: Print a description of the vertical influence
    2868. 

  2868. !          function for a given variable.
    2869. 

  2869. !********************************************************
    2870. 

  2870.   IMPLICIT NONE
    2871. 

  2871. 

  2872.   character(len=4), intent(in)  :: var      ! Variable (wind, temp, mois)
    2873. 

  2873.   real,             intent(in)  :: vif(6)   ! Vertical influence function
    2874. 

  2874.   real,             intent(in)  :: nfullmin ! Vert infl fcn full nudge min
    2875. 

  2875.   real,             intent(in)  :: nrampmin ! Vert infl fcn ramp decr min
    2876. 

  2876. 

  2877. ! Local variables
    2878. 

  2878.   character(len=200) :: msg1, msg2
    2879. 

  2879.   character(len=8) :: regime
    2880. 

  2880.   real :: nfullr1, nrampr1
    2881. 

  2881.   real :: nfullr2, nrampr2
    2882. 

  2882.   real :: nfullr4, nrampr4
    2883. 

  2883. 

  2884.   nfullr1 = vif(1)
    2885. 

  2885.   nrampr1 = vif(2)
    2886. 

  2886.   nfullr2 = vif(3)
    2887. 

  2887.   nrampr2 = vif(4)
    2888. 

  2888.   nfullr4 = vif(5)
    2889. 

  2889.   nrampr4 = vif(6)
    2890. 

  2890. 

  2891.   if(var.eq.'wind') then
    2892. 

  2892.     write(msg1,fmt='(a)') '  For winds:'
    2893. 

  2893.   elseif (var.eq.'temp') then
    2894. 

  2894.     write(msg1,fmt='(a)') '  For temperature:'
    2895. 

  2895.   elseif (var.eq.'mois') then
    2896. 

  2896.     write(msg1,fmt='(a)') '  For moisture:'
    2897. 

  2897.   else
    2898. 

  2898.     write(msg1,fmt='(a,a4)') 'Unknown variable type: ',var
    2899. 

  2899.     call wrf_message(msg1)
    2900. 

  2900.     call wrf_error_fatal ( 'print_vif_var: module_fddaobs_rtfdda STOP' )
    2901. 

  2901.   endif
    2902. 

  2902.       
    2903. 

  2903.   call wrf_message(msg1)
    2904. 

  2904. 

  2905. ! For this variable, print a description of the vif for each regime 
    2906. 

  2906.   call print_vif_regime(1, nfullr1, nrampr1, nfullmin, nrampmin) 
    2907. 

  2907.   call print_vif_regime(2, nfullr2, nrampr2, nfullmin, nrampmin) 
    2908. 

  2908.   call print_vif_regime(4, nfullr4, nrampr4, nfullmin, nrampmin) 
    2909. 

  2909. 

  2910.   END SUBROUTINE print_vif_var
    2911. 

  2911. 

  2912.   SUBROUTINE print_vif_regime(reg, nfullr, nrampr, nfullmin, nrampmin )
    2913. 

  2913. !********************************************************
    2914. 

  2914. ! Purpose: Print a description of the vertical influence
    2915. 

  2915. !          function for a given regime.
    2916. 

  2916. !********************************************************
    2917. 

  2917.   IMPLICIT NONE
    2918. 

  2918. 

  2919.   integer, intent(in)  :: reg          ! Regime number (1, 2, 4)
    2920. 

  2920.   real,    intent(in)  :: nfullr       ! Full nudge range for regime
    2921. 

  2921.   real,    intent(in)  :: nrampr       ! Rampdown range for regime
    2922. 

  2922.   real,    intent(in)  :: nfullmin     ! Vert infl fcn full nudge min
    2923. 

  2923.   real,    intent(in)  :: nrampmin     ! Vert infl fcn ramp decr min
    2924. 

  2924. 

  2925. ! Local variables
    2926. 

  2926.   character(len=200) :: msg1, msg2
    2927. 

  2927.   character(len=8) :: regime
    2928. 

  2928. 

  2929.   if(reg.eq.1) then
    2930. 

  2930.      write(regime,fmt='(a)') 'Regime 1'
    2931. 

  2931.   elseif (reg.eq.2) then
    2932. 

  2932.      write(regime,fmt='(a)') 'Regime 2'
    2933. 

  2933.   elseif (reg.eq.4) then
    2934. 

  2934.      write(regime,fmt='(a)') 'Regime 4'
    2935. 

  2935.   else
    2936. 

  2936.      write(msg1,fmt='(a,i3)') 'Unknown regime number: ',reg
    2937. 

  2937.      call wrf_message(msg1)
    2938. 

  2938.      call wrf_error_fatal ( 'print_vif_regime: module_fddaobs_rtfdda STOP' )
    2939. 

  2939.   endif
    2940. 

  2940. 

  2941. !Set msg1 for description of full weighting range
    2942. 

  2942.   if(nfullr.lt.0) then
    2943. 

  2943.      if(nfullr.eq.-5000) then
    2944. 

  2944.        write(msg1,fmt='(2x,a8,a)') regime, ': Full weighting to the PBL top'
    2945. 

  2945.      elseif (nfullr.lt.-5000) then
    2946. 

  2946.        write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(-5000.-nfullr), &
    2947. 

  2947.                                           ' m above the PBL top'
    2948. 

  2948.      else
    2949. 

  2949.        write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(nfullr+5000.),  &
    2950. 

  2950.                                           ' m below the PBL top'
    2951. 

  2951.      endif
    2952. 

  2952.   else
    2953. 

  2953.      write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting through ',                 &
    2954. 

  2954.                                      int(max(nfullr,nfullmin)),' m'
    2955. 

  2955.   endif
    2956. 

  2956. 

  2957. !Set msg2 for description of rampdown range
    2958. 

  2958.   if(nrampr.lt.0) then
    2959. 

  2959.      if(nrampr.eq.-5000) then
    2960. 

  2960.        write(msg2,fmt='(a)') ' and a vertical rampdown up to the PBL top.'
    2961. 

  2961.      elseif (nrampr.lt.-5000) then
    2962. 

  2962.        write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(-5000.-nrampr),    &
    2963. 

  2963.                             ' m above the PBL top.'
    2964. 

  2964.      else
    2965. 

  2965.        write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(nrampr+5000.),     &
    2966. 

  2966.                             ' m below the PBL top.'
    2967. 

  2967.      endif
    2968. 

  2968.   else
    2969. 

  2969.      write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown in the next ',                &
    2970. 

  2970.                           int(max(nrampr,nrampmin)),' m.'
    2971. 

  2971.   endif
    2972. 

  2972.   call wrf_message(TRIM(msg1)//msg2)
    2973. 

  2973. 

  2974.   END SUBROUTINE print_vif_regime
    2975. 

  2975. #endif
    2976. 

  2976. 

  2977. END MODULE module_fddaobs_rtfdda
    2978. 

  2978. 

  2979.