PageRenderTime 67ms CodeModel.GetById 23ms 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

Large files files are truncated, but you can click here to view the full file

    

  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.             …
    1186.

Large files files are truncated, but you can click here to view the full file