PageRenderTime 62ms CodeModel.GetById 28ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_fdda_psufddagd.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1288 lines | 885 code | 172 blank | 231 comment | 38 complexity | a954ca769b64aad37a55e85a99436c37 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !wrf:model_layer:physics
    2. 

  2. !
    3. 

  3. !
    4. 

  4. !
    5. 

  5. MODULE module_fdda_psufddagd
    6. 

  6. 

  7. CONTAINS
    8. 

  8. !
    9. 

  9. !-------------------------------------------------------------------
    10. 

  10. !
    11. 

  11.    SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
    12. 

  12.                id,analysis_interval, end_fdda_hour, &
    13. 

  13.                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
    14. 

  14.                if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
    15. 

  15.                guv, gt, gq, if_ramping, dtramp_min,  &
    16. 

  16. 

  17.                grid_sfdda, &
    18. 

  18.                analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
    19. 

  19.                rinblw, &
    20. 

  20. 

  21.                u3d,v3d,th3d,t3d,                 &
    22. 

  22.                qv3d,     &
    23. 

  23.                p3d,pi3d,                &
    24. 

  24.                u_ndg_old,v_ndg_old,t_ndg_old,q_ndg_old,mu_ndg_old,       &
    25. 

  25.                u_ndg_new,v_ndg_new,t_ndg_new,q_ndg_new,mu_ndg_new,       &
    26. 

  26.      u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
    27. 

  27.      rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
    28. 

  28.      u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
    29. 

  29.      rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
    30. 

  30.                RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
    31. 

  31.                pblh, ht, regime, znt, z, z_at_w,                             &
    32. 

  32.                ids,ide, jds,jde, kds,kde,                           &
    33. 

  33.                ims,ime, jms,jme, kms,kme,                           &
    34. 

  34.                its,ite, jts,jte, kts,kte                        )
    35. 

  35. 

  36. !-------------------------------------------------------------------
    37. 

  37.    implicit none
    38. 

  38. !-------------------------------------------------------------------
    39. 

  39. !
    40. 

  40. !   This code was implemented by Aijun Deng (Penn State).  The 3-D analysis nudging was comleted 
    41. 

  41. !   and released in December 2006.  The surface analysis nudging capability was added and 
    42. 

  42. !   released in March 2009 with WRFV3.1.
    43. 

  43. !
    44. 

  44. !-- u3d         3d u-velocity staggered on u points
    45. 

  45. !-- v3d         3d v-velocity staggered on v points
    46. 

  46. !-- th3d        3d potential temperature (k)
    47. 

  47. !-- t3d         temperature (k)
    48. 

  48. !-- qv3d        3d water vapor mixing ratio (kg/kg)
    49. 

  49. !-- p3d         3d pressure (pa)
    50. 

  50. !-- pi3d        3d exner function (dimensionless)
    51. 

  51. !-- rundgdten   staggered u tendency due to
    52. 

  52. !               fdda grid nudging (m/s/s)
    53. 

  53. !-- rvndgdten   staggered v tendency due to
    54. 

  54. !               fdda grid nudging (m/s/s)
    55. 

  55. !-- rthndgdten  theta tendency due to
    56. 

  56. !               fdda grid nudging (K/s)
    57. 

  57. !-- rqvndgdten  qv tendency due to
    58. 

  58. !               fdda grid nudging (kg/kg/s)
    59. 

  59. !-- rmundgdten  mu tendency due to
    60. 

  60. !               fdda grid nudging (Pa/s)
    61. 

  61. !-- ids         start index for i in domain
    62. 

  62. !-- ide         end index for i in domain
    63. 

  63. !-- jds         start index for j in domain
    64. 

  64. !-- jde         end index for j in domain
    65. 

  65. !-- kds         start index for k in domain
    66. 

  66. !-- kde         end index for k in domain
    67. 

  67. !-- ims         start index for i in memory
    68. 

  68. !-- ime         end index for i in memory
    69. 

  69. !-- jms         start index for j in memory
    70. 

  70. !-- jme         end index for j in memory
    71. 

  71. !-- kms         start index for k in memory
    72. 

  72. !-- kme         end index for k in memory
    73. 

  73. !-- its         start index for i in tile
    74. 

  74. !-- ite         end index for i in tile
    75. 

  75. !-- jts         start index for j in tile
    76. 

  76. !-- jte         end index for j in tile
    77. 

  77. !-- kts         start index for k in tile
    78. 

  78. !-- kte         end index for k in tile
    79. 

  79. !-------------------------------------------------------------------
    80. 

  80. !
    81. 

  81.    INTEGER,  INTENT(IN)   ::      itimestep, analysis_interval, end_fdda_hour
    82. 

  82.    INTEGER,  INTENT(IN)   ::      analysis_interval_sfc, end_fdda_hour_sfc
    83. 

  83.    INTEGER,  INTENT(IN)   ::      grid_sfdda
    84. 

  84. 

  85.    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
    86. 

  86.                                   if_no_pbl_nudging_q
    87. 

  87.    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_q
    88. 

  88.    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t,  k_zfac_q
    89. 

  89.    INTEGER,  INTENT(IN)   ::      if_ramping
    90. 

  90. 

  91.    INTEGER , INTENT(IN)   ::      id
    92. 

  92.    REAL,     INTENT(IN)   ::      DT, dx, xtime, dtramp_min
    93. 

  93. 

  94.    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
    95. 

  95.                                   ims,ime, jms,jme, kms,kme, &
    96. 

  96.                                   its,ite, jts,jte, kts,kte
    97. 

  97.  
    98. 

  98.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    99. 

  99.              INTENT(IN)   ::                   qv3d, &
    100. 

  100.                                                p3d, &
    101. 

  101.                                               pi3d, &
    102. 

  102.                                               th3d, &
    103. 

  103.                                                t3d, &
    104. 

  104.                                                  z, &
    105. 

  105.                                             z_at_w
    106. 

  106. 

  107.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    108. 

  108.              INTENT(INOUT)   ::           rundgdten, &
    109. 

  109.                                           rvndgdten, &
    110. 

  110.                                          rthndgdten, &
    111. 

  111.                                          rqvndgdten
    112. 

  112. 

  113.    REAL,     DIMENSION( ims:ime, jms:jme ), &
    114. 

  114.              INTENT(INOUT)   ::          rmundgdten
    115. 

  115. 

  116.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    117. 

  117.              INTENT(IN)      ::           u_ndg_old, &
    118. 

  118.                                           v_ndg_old, &
    119. 

  119.                                           t_ndg_old, &
    120. 

  120.                                           q_ndg_old, &
    121. 

  121.                                           u_ndg_new, &
    122. 

  122.                                           v_ndg_new, &
    123. 

  123.                                           t_ndg_new, &
    124. 

  124.                                           q_ndg_new
    125. 

  125.                                                            
    126. 

  126.    REAL,       DIMENSION( ims:ime, jms:jme ),            &   
    127. 

  127.                INTENT(IN)       ::                       u10_ndg_old,  &
    128. 

  128.                                                          v10_ndg_old,  &
    129. 

  129.                                                          t2_ndg_old,   &
    130. 

  130.                                                          th2_ndg_old,  &
    131. 

  131.                                                          q2_ndg_old,   &
    132. 

  132.                                                          rh_ndg_old,   &
    133. 

  133.                                                          psl_ndg_old,  &
    134. 

  134.                                                          ps_ndg_old,   &
    135. 

  135.                                                          u10_ndg_new,  &
    136. 

  136.                                                          v10_ndg_new,  &
    137. 

  137.                                                          t2_ndg_new,   &
    138. 

  138.                                                          th2_ndg_new,  &
    139. 

  139.                                                          q2_ndg_new,   &
    140. 

  140.                                                          rh_ndg_new,   &
    141. 

  141.                                                          psl_ndg_new,  &
    142. 

  142.                                                          ps_ndg_new
    143. 

  143. 

  144.    REAL,       DIMENSION( ims:ime, jms:jme ),            &   
    145. 

  145.                INTENT(IN)       ::                       tob_ndg_old,  &
    146. 

  146.                                                          tob_ndg_new
    147. 

  147. 

  148.    REAL,     DIMENSION( ims:ime, jms:jme ), &
    149. 

  149.              INTENT(INOUT) ::   mu_ndg_old, &
    150. 

  150.                                 mu_ndg_new
    151. 

  151. 

  152.    REAL,       DIMENSION( ims:ime, jms:jme ),            &   
    153. 

  153.                INTENT(IN)       ::                       odis_ndg_old, odis_ndg_new
    154. 

  154. 

  155.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    156. 

  156.              INTENT(IN)   ::                    u3d, &
    157. 

  157.                                                 v3d
    158. 

  158. 

  159.    REAL,  DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
    160. 

  160.                                                        ht, &
    161. 

  161.                                                        regime, &
    162. 

  162.                                                        znt
    163. 

  163. 

  164.    REAL, INTENT(IN)    :: guv, gt, gq
    165. 

  165.    REAL, INTENT(IN)    :: guv_sfc, gt_sfc, gq_sfc, rinblw
    166. 

  166. 

  167.    INTEGER             :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
    168. 

  168.    REAL                :: xtime_old, xtime_new, coef, val_analysis
    169. 

  169.    INTEGER             :: kpbl, dbg_level
    170. 

  170. 

  171.    REAL                :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
    172. 

  172.    REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ) :: wpbl  ! 1: u, 2: v, 3: t, 4: q
    173. 

  173.    REAL, DIMENSION( kts:kte, 4 )                   :: wzfac ! 1: u, 2: v, 3: t, 4: q
    174. 

  174. 

  175.    LOGICAL , EXTERNAL  :: wrf_dm_on_monitor
    176. 

  176. 

  177.    CHARACTER (LEN=256) :: message
    178. 

  178.    INTEGER :: int4
    179. 

  179. 

  180.    int4 = 1  ! 1: temporal interpolation. else: target nudging toward *_ndg_new values
    181. 

  181. 

  182.    actual_end_fdda_min = end_fdda_hour*60.0
    183. 

  183.    IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
    184. 

  184.        actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    185. 

  185.    IF( xtime > actual_end_fdda_min ) THEN
    186. 

  186. !  If xtime is greater than the end time, no need to calculate tendencies. Just set the tnedencies 
    187. 

  187. !  to zero to turn off nudging and return.
    188. 

  188.      DO j = jts, jte
    189. 

  189.      DO k = kts, kte
    190. 

  190.      DO i = its, ite
    191. 

  191.        RUNDGDTEN(i,k,j) = 0.0
    192. 

  192.        RVNDGDTEN(i,k,j) = 0.0
    193. 

  193.        RTHNDGDTEN(i,k,j) = 0.0
    194. 

  194.        RQVNDGDTEN(i,k,j) = 0.0
    195. 

  195.        IF( k .EQ. kts ) RMUNDGDTEN(i,j) = 0.
    196. 

  196.      ENDDO
    197. 

  197.      ENDDO
    198. 

  198.      ENDDO
    199. 

  199.      RETURN
    200. 

  200.    ENDIF
    201. 

  201. 

  202.    IF( analysis_interval <= 0 )CALL wrf_error_fatal('In grid FDDA, gfdda_interval_m must be > 0')
    203. 

  203.    xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
    204. 

  204.    xtime_new = xtime_old + analysis_interval
    205. 

  205.    IF( int4 == 1 ) THEN
    206. 

  206.      coef = (xtime-xtime_old)/(xtime_new-xtime_old)
    207. 

  207.    ELSE
    208. 

  208.      coef = 1.0          ! Nudging toward a target value (*_ndg_new values)
    209. 

  209.    ENDIF
    210. 

  210. 

  211.    IF ( wrf_dm_on_monitor()) THEN
    212. 

  212. 

  213.      CALL get_wrf_debug_level( dbg_level )
    214. 

  214. 

  215.      IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
    216. 

  216. 

  217.        IF( xtime < end_fdda_hour*60.0 ) THEN
    218. 

  218.          WRITE(message,'(a,i1,a,f10.3,a)') &
    219. 

  219.           'D0',id,' 3-D analysis nudging reads new data at time = ', xtime, ' min.'
    220. 

  220.          CALL wrf_message( TRIM(message) )
    221. 

  221.          WRITE(message,'(a,i1,a,2f8.2,a)') &
    222. 

  222.           'D0',id,' 3-D analysis nudging bracketing times = ', xtime_old, xtime_new, ' min.'
    223. 

  223.          CALL wrf_message( TRIM(message) )
    224. 

  224.        ENDIF
    225. 

  225. 

  226.        actual_end_fdda_min = end_fdda_hour*60.0
    227. 

  227.        IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
    228. 

  228.            actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    229. 

  229. 

  230.        IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    231. 

  231. !        Find the mid point of the tile and print out the sample values
    232. 

  232.          i0 = (ite-its)/2+its
    233. 

  233.          j0 = (jte-jts)/2+jts 
    234. 

  234. 

  235.          IF( guv > 0.0 ) THEN
    236. 

  236.            DO k = kts, kte
    237. 

  237.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    238. 

  238.                '    D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
    239. 

  239.                ' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
    240. 

  240.              CALL wrf_message( TRIM(message) )
    241. 

  241.            ENDDO
    242. 

  242.            WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    243. 

  243.              '    D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
    244. 

  244.              ' mu_ndg_old=', mu_ndg_old(i0,j0), ' mu_ndg_new=', mu_ndg_new(i0,j0)
    245. 

  245.            CALL wrf_message( TRIM(message) )
    246. 

  246.            DO k = kts, kte
    247. 

  247.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    248. 

  248.                '    D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
    249. 

  249.                ' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
    250. 

  250.              CALL wrf_message( TRIM(message) )
    251. 

  251.            ENDDO
    252. 

  252.          ENDIF
    253. 

  253. 

  254.          IF( gt > 0.0 ) THEN
    255. 

  255.            DO k = kts, kte
    256. 

  256.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    257. 

  257.                '    D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
    258. 

  258.                ' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
    259. 

  259.              CALL wrf_message( TRIM(message) )
    260. 

  260.            ENDDO
    261. 

  261.          ENDIF
    262. 

  262. 

  263.          IF( gq > 0.0 ) THEN
    264. 

  264.            DO k = kts, kte
    265. 

  265.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    266. 

  266.                '    D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
    267. 

  267.                ' q_ndg_old=', q_ndg_old(i0,k,j0), ' q_ndg_new=', q_ndg_new(i0,k,j0)
    268. 

  268.              CALL wrf_message( TRIM(message) )
    269. 

  269.            ENDDO
    270. 

  270.          ENDIF
    271. 

  271. 

  272.         IF( int4 == 1 ) then
    273. 

  273.            WRITE(message,'(a,i1,a)') '    D0',id, &
    274. 

  274.               ' 3-D nudging towards the temporally interpolated analysis'
    275. 

  275.          ELSE
    276. 

  276.            WRITE(message,'(a,i1,a)') '    D0',id, &
    277. 

  277.               ' 3-D nudging towards the target analysis'
    278. 

  278.          ENDIF
    279. 

  279. 

  280.        ENDIF
    281. 

  281.      ENDIF
    282. 

  282.    ENDIF
    283. 

  283. 

  284.    jtsv=MAX0(jts,jds+1)
    285. 

  285.    itsu=MAX0(its,ids+1)
    286. 

  286. 

  287.    jtf=MIN0(jte,jde-1)
    288. 

  288.    ktf=MIN0(kte,kde-1)
    289. 

  289.    itf=MIN0(ite,ide-1)
    290. 

  290. !
    291. 

  291. ! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t, 
    292. 

  292. ! if_no_pbl_nudging_q swithes) are set to 1, compute the weighting function, wpbl(:,k,:,:), 
    293. 

  293. ! based on the PBL depth.  wpbl = 1 above the PBL and wpbl = 0 in the PBL.  If all 
    294. 

  294. ! the switche are set to zero, wpbl = 1 (default value).
    295. 

  295. !
    296. 

  296.    wpbl(:,:,:,:) = 1.0
    297. 

  297. 

  298.    IF( if_no_pbl_nudging_uv == 1 .OR. grid_sfdda == 1 ) THEN
    299. 

  299. 

  300.      DO j=jts,jtf 
    301. 

  301.      DO i=itsu,itf
    302. 

  302. 

  303.        kpbl = 1
    304. 

  304.        zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
    305. 

  305. 

  306.        loop_ku: DO k=kts,ktf 
    307. 

  307. !        zagl     = 0.5 * ( z(i-1,k,j)-ht(i-1,j) + z(i,k,j)-ht(i,j) )
    308. 

  308.          zagl_bot = 0.5 * ( z_at_w(i-1,k,  j)-ht(i-1,j) + z_at_w(i,k,  j)-ht(i,j) )
    309. 

  309.          zagl_top = 0.5 * ( z_at_w(i-1,k+1,j)-ht(i-1,j) + z_at_w(i,k+1,j)-ht(i,j) )
    310. 

  310.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    311. 

  311.            kpbl = k
    312. 

  312.            EXIT loop_ku
    313. 

  313.          ENDIF
    314. 

  314.        ENDDO loop_ku
    315. 

  315. 

  316.        DO k=kts,ktf 
    317. 

  317.          IF( k <= kpbl   ) wpbl(i, k, j, 1) = 0.0
    318. 

  318.          IF( k == kpbl+1 ) wpbl(i, k, j, 1) = 0.1
    319. 

  319.          IF( k >  kpbl+1 ) wpbl(i, k, j, 1) = 1.0
    320. 

  320.        ENDDO
    321. 

  321. 

  322.      ENDDO
    323. 

  323.      ENDDO
    324. 

  324. 

  325.      DO i=its,itf
    326. 

  326.      DO j=jtsv,jtf
    327. 

  327. 

  328.        kpbl = 1
    329. 

  329.        zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
    330. 

  330. 

  331.        loop_kv: DO k=kts,ktf
    332. 

  332. !        zagl     = 0.5 * ( z(i,k,j-1)-ht(i,j-1) + z(i,k,j)-ht(i,j) )
    333. 

  333.          zagl_bot = 0.5 * ( z_at_w(i,k,  j-1)-ht(i,j-1) + z_at_w(i,k,  j)-ht(i,j) )
    334. 

  334.          zagl_top = 0.5 * ( z_at_w(i,k+1,j-1)-ht(i,j-1) + z_at_w(i,k+1,j)-ht(i,j) )
    335. 

  335.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    336. 

  336.            kpbl = k
    337. 

  337.            EXIT loop_kv
    338. 

  338.          ENDIF
    339. 

  339.        ENDDO loop_kv
    340. 

  340. 

  341.        DO k=kts,ktf
    342. 

  342.          IF( k <= kpbl   ) wpbl(i, k, j, 2) = 0.0
    343. 

  343.          IF( k == kpbl+1 ) wpbl(i, k, j, 2) = 0.1
    344. 

  344.          IF( k >  kpbl+1 ) wpbl(i, k, j, 2) = 1.0
    345. 

  345.        ENDDO
    346. 

  346. 

  347.      ENDDO
    348. 

  348.      ENDDO
    349. 

  349. 

  350.    ENDIF
    351. 

  351. 

  352.    IF( if_no_pbl_nudging_t == 1 .OR. grid_sfdda == 1 ) THEN
    353. 

  353.    
    354. 

  354.      DO j=jts,jtf
    355. 

  355.      DO i=its,itf
    356. 

  356. 

  357.        kpbl = 1
    358. 

  358.        zpbl = pblh(i,j)
    359. 

  359.         
    360. 

  360.        loop_kt: DO k=kts,ktf
    361. 

  361. !        zagl     = z(i,k,j)-ht(i,j)
    362. 

  362.          zagl_bot = z_at_w(i,k,  j)-ht(i,j)
    363. 

  363.          zagl_top = z_at_w(i,k+1,j)-ht(i,j)
    364. 

  364.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    365. 

  365.            kpbl = k
    366. 

  366.            EXIT loop_kt
    367. 

  367.          ENDIF
    368. 

  368.        ENDDO loop_kt
    369. 

  369. 

  370.        DO k=kts,ktf
    371. 

  371.          IF( k <= kpbl   ) wpbl(i, k, j, 3) = 0.0
    372. 

  372.          IF( k == kpbl+1 ) wpbl(i, k, j, 3) = 0.1
    373. 

  373.          IF( k >  kpbl+1 ) wpbl(i, k, j, 3) = 1.0
    374. 

  374.        ENDDO 
    375. 

  375.         
    376. 

  376.      ENDDO
    377. 

  377.      ENDDO
    378. 

  378. 

  379.    ENDIF
    380. 

  380. 

  381.    IF( if_no_pbl_nudging_q == 1 .OR. grid_sfdda == 1 ) THEN
    382. 

  382.    
    383. 

  383.      DO j=jts,jtf
    384. 

  384.      DO i=its,itf
    385. 

  385. 

  386.        kpbl = 1
    387. 

  387.        zpbl = pblh(i,j)
    388. 

  388.           
    389. 

  389.        loop_kq: DO k=kts,ktf
    390. 

  390. !        zagl     = z(i,k,j)-ht(i,j)
    391. 

  391.          zagl_bot = z_at_w(i,k,  j)-ht(i,j)
    392. 

  392.          zagl_top = z_at_w(i,k+1,j)-ht(i,j)
    393. 

  393.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    394. 

  394.            kpbl = k
    395. 

  395.            EXIT loop_kq
    396. 

  396.          ENDIF
    397. 

  397.        ENDDO loop_kq
    398. 

  398. 

  399.        DO k=kts,ktf
    400. 

  400.          IF( k <= kpbl   ) wpbl(i, k, j, 4) = 0.0
    401. 

  401.          IF( k == kpbl+1 ) wpbl(i, k, j, 4) = 0.1
    402. 

  402.          IF( k >  kpbl+1 ) wpbl(i, k, j, 4) = 1.0
    403. 

  403.        ENDDO 
    404. 

  404.             
    405. 

  405.      ENDDO  
    406. 

  406.      ENDDO
    407. 

  407.         
    408. 

  408.    ENDIF
    409. 

  409. !
    410. 

  410. ! If the user-defined namelist switches (if_zfac_uv, if_zfac_t,
    411. 

  411. ! if_zfac_q swithes) are set to 1, compute the weighting function, wzfac(k,:),
    412. 

  412. ! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_q) below which analysis 
    413. 

  413. ! nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x).  If all
    414. 

  414. ! the switche are set to zero, wzfac = 1 (default value).
    415. 

  415. !
    416. 

  416.    wzfac(:,:) = 1.0
    417. 

  417. 

  418.    IF( if_zfac_uv == 1 ) THEN
    419. 

  419. 

  420.      DO j=jts,jtf
    421. 

  421.      DO i=itsu,itf
    422. 

  422.      DO k=kts,ktf
    423. 

  423.        IF( k <= k_zfac_uv   ) wzfac(k, 1:2) = 0.0
    424. 

  424.        IF( k == k_zfac_uv+1 ) wzfac(k, 1:2) = 0.1
    425. 

  425.        IF( k >  k_zfac_uv+1 ) wzfac(k, 1:2) = 1.0
    426. 

  426.      ENDDO
    427. 

  427.      ENDDO
    428. 

  428.      ENDDO
    429. 

  429. 

  430.    ENDIF
    431. 

  431. 

  432.    IF( if_zfac_t == 1 ) THEN
    433. 

  433. 

  434.      DO j=jts,jtf
    435. 

  435.      DO i=itsu,itf
    436. 

  436.      DO k=kts,ktf
    437. 

  437.        IF( k <= k_zfac_t   ) wzfac(k, 3) = 0.0
    438. 

  438.        IF( k == k_zfac_t+1 ) wzfac(k, 3) = 0.1
    439. 

  439.        IF( k >  k_zfac_t+1 ) wzfac(k, 3) = 1.0
    440. 

  440.      ENDDO
    441. 

  441.      ENDDO
    442. 

  442.      ENDDO
    443. 

  443. 

  444.    ENDIF
    445. 

  445. 

  446.    IF( if_zfac_q == 1 ) THEN
    447. 

  447.        
    448. 

  448.      DO j=jts,jtf
    449. 

  449.      DO i=itsu,itf 
    450. 

  450.      DO k=kts,ktf
    451. 

  451.        IF( k <= k_zfac_q   ) wzfac(k, 4) = 0.0
    452. 

  452.        IF( k == k_zfac_q+1 ) wzfac(k, 4) = 0.1 
    453. 

  453.        IF( k >  k_zfac_q+1 ) wzfac(k, 4) = 1.0
    454. 

  454.      ENDDO  
    455. 

  455.      ENDDO
    456. 

  456.      ENDDO
    457. 

  457. 

  458.    ENDIF
    459. 

  459. !
    460. 

  460. ! If if_ramping and dtramp_min are defined by user, comput a time weighting function, tfac, 
    461. 

  461. ! for analysis nudging so that at the end of the nudging period (which has to be at a 
    462. 

  462. ! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
    463. 

  463. !
    464. 

  464. ! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at 
    465. 

  465. ! end_fdda_hour-ABS(dtramp_min).  
    466. 

  466. !
    467. 

  467. ! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at 
    468. 

  468. ! end_fdda_hour+ABS(dtramp_min). In this case, the obs values are extrapolated using 
    469. 

  469. ! the obs tendency saved from the previous FDDA wondow.  More specifically for extrapolation, 
    470. 

  470. ! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
    471. 

  471. !
    472. 

  472.    tfac = 1.0
    473. 

  473. 

  474.    IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
    475. 

  475.  
    476. 

  476.      IF( dtramp_min <= 0.0 ) THEN
    477. 

  477.        actual_end_fdda_min = end_fdda_hour*60.0
    478. 

  478.      ELSE
    479. 

  479.        actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
    480. 

  480.      ENDIF
    481. 

  481. 

  482.      IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN 
    483. 

  483.        tfac = 1.0
    484. 

  484.      ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
    485. 

  485.        tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
    486. 

  486.        IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/(analysis_interval*1.0)
    487. 

  487.      ELSE                                                     
    488. 

  488.        tfac = 0.0
    489. 

  489.      ENDIF
    490. 

  490. 

  491.    ENDIF                                                  
    492. 

  492. !
    493. 

  493. ! Surface Analysis Nudging
    494. 

  494. !
    495. 

  495.    IF( grid_sfdda == 1 ) THEN
    496. 

  496.      CALL SFDDAGD(itimestep,dx,dt,xtime, id, &
    497. 

  497.      analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
    498. 

  498.      rinblw, &
    499. 

  499.                u3d,v3d,th3d,t3d,                 &
    500. 

  500.                qv3d,     &
    501. 

  501.                p3d,pi3d,        &
    502. 

  502.      u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
    503. 

  503.      rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old,  &
    504. 

  504.      u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
    505. 

  505.      rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new,  &
    506. 

  506.      RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
    507. 

  507.      pblh, ht, regime, znt, z, z_at_w,                             &
    508. 

  508.      ids,ide, jds,jde, kds,kde,                           &
    509. 

  509.      ims,ime, jms,jme, kms,kme,                           &
    510. 

  510.      its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
    511. 

  511.      actual_end_fdda_min, tfac )
    512. 

  512.    ENDIF
    513. 

  513. !
    514. 

  514. ! Compute 3-D nudging tendencies for u, v, t and q
    515. 

  515. !
    516. 

  516.    DO j=jts,jtf
    517. 

  517.    DO k=kts,ktf
    518. 

  518.    DO i=itsu,itf
    519. 

  519.      val_analysis = u_ndg_old(i,k,j) *( 1.0 - coef ) + u_ndg_new(i,k,j) * coef
    520. 

  520.      RUNDGDTEN(i,k,j) = RUNDGDTEN(i,k,j) + guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
    521. 

  521.                          ( val_analysis - u3d(i,k,j) )
    522. 

  522.    ENDDO
    523. 

  523.    ENDDO
    524. 

  524.    ENDDO
    525. 

  525. 

  526.    DO j=jtsv,jtf
    527. 

  527.    DO k=kts,ktf
    528. 

  528.    DO i=its,itf
    529. 

  529.      val_analysis = v_ndg_old(i,k,j) *( 1.0 - coef ) + v_ndg_new(i,k,j) * coef
    530. 

  530.      RVNDGDTEN(i,k,j) = RVNDGDTEN(i,k,j) + guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
    531. 

  531.                        ( val_analysis - v3d(i,k,j) )
    532. 

  532.    ENDDO
    533. 

  533.    ENDDO
    534. 

  534.    ENDDO
    535. 

  535. 

  536.    DO j=jts,jtf
    537. 

  537.    DO k=kts,ktf
    538. 

  538.    DO i=its,itf
    539. 

  539.      val_analysis = t_ndg_old(i,k,j) *( 1.0 - coef ) + t_ndg_new(i,k,j) * coef
    540. 

  540.      RTHNDGDTEN(i,k,j) = RTHNDGDTEN(i,k,j) +   gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
    541. 

  541.                           ( val_analysis - th3d(i,k,j) + 300.0 )
    542. 

  542. 

  543.      val_analysis = q_ndg_old(i,k,j) *( 1.0 - coef ) + q_ndg_new(i,k,j) * coef
    544. 

  544.      RQVNDGDTEN(i,k,j) = RQVNDGDTEN(i,k,j) + gq * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
    545. 

  545.                           ( val_analysis - qv3d(i,k,j) )
    546. 

  546.    ENDDO
    547. 

  547.    ENDDO
    548. 

  548.    ENDDO
    549. 

  549. 

  550.    END SUBROUTINE fddagd
    551. 

  551. 

  552. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    553. 

  553.    SUBROUTINE sfddagd(itimestep,dx,dt,xtime,  &
    554. 

  554.                id, analysis_interval_sfc, end_fdda_hour_sfc, &
    555. 

  555.                guv_sfc, gt_sfc, gq_sfc, rinblw,  &
    556. 

  556.                u3d,v3d,th3d,t3d,                 &
    557. 

  557.                qv3d,     &
    558. 

  558.                p3d,pi3d,                &
    559. 

  559.                u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
    560. 

  560.                rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old,  &
    561. 

  561.                u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
    562. 

  562.                rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new,  &
    563. 

  563.                RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN, &
    564. 

  564.                pblh, ht, regime, znt, z, z_at_w,                             &
    565. 

  565.                ids,ide, jds,jde, kds,kde,                           &
    566. 

  566.                ims,ime, jms,jme, kms,kme,                           &
    567. 

  567.                its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
    568. 

  568.                actual_end_fdda_min, tfac)
    569. 

  569. 

  570. !-------------------------------------------------------------------
    571. 

  571.    USE module_model_constants
    572. 

  572.  
    573. 

  573.    implicit none
    574. 

  574. 

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

  576. !  
    577. 

  577. !   This code was implemented by Aijun Deng (Penn State).  The 3-D analysis nudging was comleted
    578. 

  578. !   and released in December 2006.  The surface analysis nudging capability was added and
    579. 

  579. !   released in March 2009 with WRFV3.1.
    580. 

  580. !
    581. 

  581. !-- u3d         3d u-velocity staggered on u points
    582. 

  582. !-- v3d         3d v-velocity staggered on v points
    583. 

  583. !-- th3d        3d potential temperature (k)
    584. 

  584. !-- t3d         temperature (k)
    585. 

  585. !-- qv3d        3d water vapor mixing ratio (kg/kg)
    586. 

  586. !-- p3d         3d pressure (pa)
    587. 

  587. !-- pi3d        3d exner function (dimensionless)
    588. 

  588. !-- rundgdten   staggered u tendency due to
    589. 

  589. !               fdda grid nudging (m/s/s)
    590. 

  590. !-- rvndgdten   staggered v tendency due to
    591. 

  591. !               fdda grid nudging (m/s/s)
    592. 

  592. !-- rthndgdten  theta tendency due to
    593. 

  593. !               fdda grid nudging (K/s)
    594. 

  594. !-- rqvndgdten  qv tendency due to
    595. 

  595. !               fdda grid nudging (kg/kg/s)
    596. 

  596. !-- rmundgdten  mu tendency due to
    597. 

  597. !               fdda grid nudging (Pa/s)
    598. 

  598. !-- ids         start index for i in domain
    599. 

  599. !-- ide         end index for i in domain
    600. 

  600. !-- jds         start index for j in domain
    601. 

  601. !-- jde         end index for j in domain
    602. 

  602. !-- kds         start index for k in domain
    603. 

  603. !-- kde         end index for k in domain
    604. 

  604. !-- ims         start index for i in memory
    605. 

  605. !-- ime         end index for i in memory
    606. 

  606. !-- jms         start index for j in memory
    607. 

  607. !-- jme         end index for j in memory
    608. 

  608. !-- kms         start index for k in memory
    609. 

  609. !-- kme         end index for k in memory
    610. 

  610. !-- its         start index for i in tile
    611. 

  611. !-- ite         end index for i in tile
    612. 

  612. !-- jts         start index for j in tile
    613. 

  613. !-- jte         end index for j in tile
    614. 

  614. !-- kts         start index for k in tile
    615. 

  615. !-- kte         end index for k in tile
    616. 

  616. !-------------------------------------------------------------------
    617. 

  617. !
    618. 

  618.    INTEGER,  INTENT(IN)   ::      itimestep, analysis_interval_sfc, end_fdda_hour_sfc
    619. 

  619. 

  620.    INTEGER , INTENT(IN)   ::      id
    621. 

  621.    REAL,     INTENT(IN)   ::      dx,DT, xtime
    622. 

  622. 

  623.    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
    624. 

  624.                                   ims,ime, jms,jme, kms,kme, &
    625. 

  625.                                   its,ite, jts,jte, kts,kte
    626. 

  626.  
    627. 

  627.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    628. 

  628.              INTENT(IN)   ::                   qv3d, &
    629. 

  629.                                                p3d, &
    630. 

  630.                                               pi3d, &
    631. 

  631.                                               th3d, &
    632. 

  632.                                                t3d, &
    633. 

  633.                                                  z, &
    634. 

  634.                                             z_at_w
    635. 

  635. 

  636.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    637. 

  637.              INTENT(INOUT)   ::           rundgdten, &
    638. 

  638.                                           rvndgdten, &
    639. 

  639.                                          rthndgdten, &
    640. 

  640.                                          rqvndgdten
    641. 

  641. 

  642.    REAL,     DIMENSION( ims:ime, jms:jme ), &
    643. 

  643.              INTENT(INOUT)   ::          rmundgdten
    644. 

  644. 

  645.    REAL,       DIMENSION( ims:ime, jms:jme ),            &
    646. 

  646.                INTENT(IN)       ::                       u10_ndg_old,  &
    647. 

  647.                                                          v10_ndg_old,  &
    648. 

  648.                                                          t2_ndg_old,   &
    649. 

  649.                                                          th2_ndg_old,  &
    650. 

  650.                                                          q2_ndg_old,   &
    651. 

  651.                                                          rh_ndg_old,   &
    652. 

  652.                                                          psl_ndg_old,  &
    653. 

  653.                                                          ps_ndg_old,   &
    654. 

  654.                                                          u10_ndg_new,  &
    655. 

  655.                                                          v10_ndg_new,  &
    656. 

  656.                                                          t2_ndg_new,   &
    657. 

  657.                                                          th2_ndg_new,  &
    658. 

  658.                                                          q2_ndg_new,   &
    659. 

  659.                                                          rh_ndg_new,   &
    660. 

  660.                                                          psl_ndg_new,  &
    661. 

  661.                                                          ps_ndg_new
    662. 

  662. 

  663.    REAL,       DIMENSION( ims:ime, jms:jme ),            &
    664. 

  664.                INTENT(IN)       ::                       tob_ndg_old,  &
    665. 

  665.                                                          tob_ndg_new
    666. 

  666. 

  667.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    668. 

  668.              INTENT(IN)   ::                    u3d, &
    669. 

  669.                                                 v3d
    670. 

  670. 

  671.    REAL,       DIMENSION( ims:ime, jms:jme ),            &
    672. 

  672.                INTENT(IN)         ::                    odis_ndg_old, odis_ndg_new
    673. 

  673. 

  674.    REAL,  DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
    675. 

  675.                                                        ht,   &
    676. 

  676.                                                        regime, &
    677. 

  677.                                                        znt
    678. 

  678. 

  679.    REAL, INTENT(IN)    :: guv_sfc, gt_sfc, gq_sfc, rinblw
    680. 

  680. 

  681.    INTEGER             :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, j0
    682. 

  682.    REAL                :: xtime_old_sfc, xtime_new_sfc, coef, val_analysis, es
    683. 

  683.    INTEGER             :: kpbl, dbg_level
    684. 

  684. 

  685.    REAL                :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
    686. 

  686. 

  687.    REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ), &
    688. 

  688.          INTENT(IN)                                   :: wpbl  ! 1: u, 2: v, 3: t, 4: q
    689. 

  689.    REAL, DIMENSION( kts:kte, 4 ),  &
    690. 

  690.          INTENT(IN)                                   :: wzfac ! 1: u, 2: v, 3: t, 4: q
    691. 

  691.    REAL, DIMENSION( its:ite, jts:jte)                 :: wndcor_u, wndcor_v
    692. 

  692.    REAL, DIMENSION( its-2:ite+2, jts-2:jte+2)         :: blw_old, blw_new
    693. 

  693.    REAL, DIMENSION( its:ite, kts:kte, jts:jte)        :: qsat
    694. 

  694. 

  695.    REAL                                               :: m, b=1.8, blw, rindx, x
    696. 

  696. 

  697.    REAL :: difz, wr14, wr1z, wr24, wr2z, wndfac, reg, znt0
    698. 

  698.    INTEGER,  INTENT(IN)   :: if_ramping
    699. 

  699.    REAL, INTENT(IN)       :: dtramp_min
    700. 

  700. 

  701.    LOGICAL , EXTERNAL     :: wrf_dm_on_monitor
    702. 

  702. 

  703.    CHARACTER (LEN=256)    :: message
    704. 

  704.    INTEGER :: iwinds, idd, iqsat, int4
    705. 

  705. 

  706.    iwinds = 1    ! 1: Scale the surface wind analysis to the lowest model level, 
    707. 

  707.                  ! if the first model half-layer is greater than 10 meters 
    708. 

  708.                  ! and we are in the free convection regime (REGIME=4.0). else: no
    709. 

  709.    idd = 1       ! 1: Obs data density correction is applied. else: no
    710. 

  710.    iqsat = 1     ! 1: Remove super saturation. eles: no
    711. 

  711.    int4 = 1      ! 1: temporal ionterpolation. else: target nudging toward *_ndg_new values
    712. 

  712. 

  713.    IF( analysis_interval_sfc <= 0 )CALL wrf_error_fatal('In grid sfc FDDA, sgfdda_interval_m must be > 0')
    714. 

  714.    xtime_old_sfc = FLOOR(xtime/analysis_interval_sfc) * analysis_interval_sfc * 1.0
    715. 

  715.    xtime_new_sfc = xtime_old_sfc + analysis_interval_sfc
    716. 

  716.    IF( int4 == 1 ) THEN
    717. 

  717.      coef = (xtime-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc)  ! Temporal interpolation
    718. 

  718.    ELSE
    719. 

  719.      coef = 1.0          ! Nudging toward a target value (*_ndg_new values)
    720. 

  720.    ENDIF
    721. 

  721. 

  722.    IF ( wrf_dm_on_monitor()) THEN
    723. 

  723. 

  724.      CALL get_wrf_debug_level( dbg_level )
    725. 

  725. 

  726.      IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
    727. 

  727. 

  728.        IF( xtime < end_fdda_hour_sfc*60.0 ) THEN
    729. 

  729.          WRITE(message,'(a,i1,a,f10.3,a)') &
    730. 

  730.           'D0',id,' surface analysis nudging reads new data at time = ', xtime, ' min.'
    731. 

  731.          CALL wrf_message( TRIM(message) )
    732. 

  732.          WRITE(message,'(a,i1,a,2f8.2,a)') &
    733. 

  733.           'D0',id,' surface analysis nudging bracketing times = ', xtime_old_sfc, xtime_new_sfc, ' min.'
    734. 

  734.          CALL wrf_message( TRIM(message) )
    735. 

  735.        ENDIF
    736. 

  736. 

  737.        IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    738. 

  738. !        Find the mid point of the tile and print out the sample values
    739. 

  739.          i0 = (ite-its)/2+its
    740. 

  740.          j0 = (jte-jts)/2+jts 
    741. 

  741. 

  742.          IF( guv_sfc > 0.0 ) THEN
    743. 

  743.            WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
    744. 

  744.              '    D0',id,' sample surface analysis values at i,j=', i0, j0, &
    745. 

  745.              ' u10_ndg_old=', u10_ndg_old(i0,j0), ' u10_ndg_new=', u10_ndg_new(i0,j0)
    746. 

  746.            CALL wrf_message( TRIM(message) )
    747. 

  747.            WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
    748. 

  748.              '    D0',id,' sample surface analysis values at i,j=', i0, j0, &
    749. 

  749.              ' v10_ndg_old=', v10_ndg_old(i0,j0), ' v10_ndg_new=', v10_ndg_new(i0,j0)
    750. 

  750.            CALL wrf_message( TRIM(message) )
    751. 

  751.          ENDIF
    752. 

  752. 

  753.          IF( gt_sfc > 0.0 ) THEN
    754. 

  754.            WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
    755. 

  755.              '    D0',id,' sample surface analysis values at i,j=', i0, j0, &
    756. 

  756.              ' th2_ndg_old=', th2_ndg_old(i0,j0), ' th2_ndg_new=', th2_ndg_new(i0,j0)
    757. 

  757.            CALL wrf_message( TRIM(message) )
    758. 

  758.          ENDIF
    759. 

  759. 

  760.          IF( gq_sfc > 0.0 ) THEN
    761. 

  761.            WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
    762. 

  762.              '    D0',id,' sample surface analysis values at i,j=', i0, j0, &
    763. 

  763.              ' q2_ndg_old=', q2_ndg_old(i0,j0), ' q2_ndg_new=', q2_ndg_new(i0,j0)
    764. 

  764.            CALL wrf_message( TRIM(message) )
    765. 

  765.          ENDIF
    766. 

  766. 

  767.          IF( iwinds ==  1 ) &
    768. 

  768.            WRITE(message,'(a,i1,a)') '    D0',id, &
    769. 

  769.               ' surface wind analysis s scaled to the lowest model level, if dz1 > 10m and REGIME=4.'
    770. 

  770. 

  771.          IF( idd ==  1 ) &
    772. 

  772.            WRITE(message,'(a,i1,a)') '    D0',id, &
    773. 

  773.               ' obs data density is used for additional weighting function'
    774. 

  774. 

  775.          IF( iqsat ==  1 ) &
    776. 

  776.            WRITE(message,'(a,i1,a)') '    D0',id, &
    777. 

  777.               ' super saturation is not allowed for q analysis'
    778. 

  778. 

  779.          IF( int4 ==  1 ) then
    780. 

  780.            WRITE(message,'(a,i1,a)') '    D0',id, &
    781. 

  781.               ' surface nudging towards the temporally interpolated analysis'
    782. 

  782.          ELSE
    783. 

  783.            WRITE(message,'(a,i1,a)') '    D0',id, &
    784. 

  784.               ' surface nudging towards the target analysis'
    785. 

  785.          ENDIF
    786. 

  786. 

  787.        ENDIF
    788. 

  788.      ENDIF
    789. 

  789.    ENDIF
    790. 

  790. 

  791. 

  792.    jtsv=MAX0(jts,jds+1)
    793. 

  793.    itsu=MAX0(its,ids+1)
    794. 

  794. 

  795.    jtf=MIN0(jte,jde-1)
    796. 

  796.    ktf=MIN0(kte,kde-1)
    797. 

  797.    itf=MIN0(ite,ide-1)
    798. 

  798. 

  799. !
    800. 

  800. ! Compute the vertical weighting function to scale the surface wind analysis to 
    801. 

  801. ! the lowest model level, if the first model half-layer is greater
    802. 

  802. ! than 10 meters and we are in the free convection regime (REGIME=4.0). 
    803. 

  803. !
    804. 

  804.    IF( iwinds == 1 ) THEN
    805. 

  805.      wndcor_u(:,:) = 1.0
    806. 

  806.      DO j=jts,jtf
    807. 

  807.      DO i=itsu,itf
    808. 

  808.        reg =  0.5 * ( regime(i-1,  j) + regime(i,  j) )
    809. 

  809.        difz = 0.5 * ( z(i-1,1,j) - ht(i-1,j)  &
    810. 

  810.                     + z(i,  1,j) - ht(i,  j)    ) 
    811. 

  811.        IF( reg > 3.5 .AND. difz > 10.0 ) THEN
    812. 

  812.          znt0 = 0.5 * (    znt(i-1,  j) +    znt(i,  j) )
    813. 

  813.          IF( znt0 <= 0.2) THEN
    814. 

  814.            wndcor_u(i,j) = 1.0+0.320*znt0**0.2
    815. 

  815.          ELSE
    816. 

  816.            wndcor_u(i,j) = 1.169+0.315*znt0
    817. 

  817.          ENDIF
    818. 

  818. 

  819.          wr14 = log(40.0/0.05)
    820. 

  820.          wr1z = log(difz/0.05)
    821. 

  821.          wr24 = log(40.0/1.0)
    822. 

  822.          wr2z = log(difz/1.0)
    823. 

  823.          wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24) 
    824. 

  824.          wndcor_u(i,j) = wndfac*wndcor_u(i,j)
    825. 

  825.        ENDIF
    826. 

  826.      ENDDO
    827. 

  827.      ENDDO
    828. 

  828. 

  829.      IF ( wrf_dm_on_monitor()) THEN
    830. 

  830.        IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
    831. 

  831.          IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    832. 

  832.            i0 = (ite-its)/2+its
    833. 

  833.            j0 = (jte-jts)/2+jts
    834. 

  834.            WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
    835. 

  835.              '    D0',id,' sample wndcor_u values at i,j=', i0, j0, &
    836. 

  836.              ' wndcor_u=', wndcor_u(i0,j0)
    837. 

  837.            CALL wrf_message( TRIM(message) )
    838. 

  838.          ENDIF
    839. 

  839.        ENDIF
    840. 

  840.      ENDIF
    841. 

  841.    ELSE 
    842. 

  842.      wndcor_u(:,:) = 1.0
    843. 

  843.    ENDIF
    844. 

  844. 

  845.    IF( iwinds == 1 ) THEN
    846. 

  846.      wndcor_v(:,:) = 1.0
    847. 

  847.      DO j=jtsv,jtf
    848. 

  848.      DO i=its,itf
    849. 

  849.        reg =  0.5 * ( regime(i,  j-1) + regime(i,  j) )
    850. 

  850.        difz = 0.5 * ( z(i,1,j-1) - ht(i,j-1)  &
    851. 

  851.                   +   z(i,1,j  ) - ht(i,j  )    )
    852. 

  852.        IF( reg > 3.5 .AND. difz > 10.0 ) THEN 
    853. 

  853.          znt0 = 0.5 * (    znt(i,  j-1) +    znt(i,  j) )
    854. 

  854.          IF( znt0 <= 0.2) THEN
    855. 

  855.            wndcor_v(i,j) = 1.0+0.320*znt0**0.2
    856. 

  856.          ELSE
    857. 

  857.            wndcor_v(i,j) = 1.169+0.315*znt0
    858. 

  858.          ENDIF
    859. 

  859.        
    860. 

  860.          wr14 = log(40.0/0.05)
    861. 

  861.          wr1z = log(difz/0.05)
    862. 

  862.          wr24 = log(40.0/1.0)
    863. 

  863.          wr2z = log(difz/1.0)
    864. 

  864.          wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
    865. 

  865.            wndcor_v(i,j) = wndfac*wndcor_v(i,j)
    866. 

  866.        ENDIF
    867. 

  867.      ENDDO
    868. 

  868.      ENDDO
    869. 

  869.      IF ( wrf_dm_on_monitor()) THEN
    870. 

  870.        IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
    871. 

  871.          IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    872. 

  872.            i0 = (ite-its)/2+its
    873. 

  873.            j0 = (jte-jts)/2+jts
    874. 

  874.            WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
    875. 

  875.              '    D0',id,' sample wndcor_v values at i,j=', i0, j0, &
    876. 

  876.              ' wndcor_v=', wndcor_v(i0,j0)
    877. 

  877.            CALL wrf_message( TRIM(message) )
    878. 

  878.          ENDIF
    879. 

  879.        ENDIF
    880. 

  880.      ENDIF
    881. 

  881.    ELSE 
    882. 

  882.      wndcor_v(:,:) = 1.0
    883. 

  883.    ENDIF
    884. 

  884. !
    885. 

  885. ! Compute saturation mixing ratio so that nudging to a super-saturated state
    886. 

  886. ! is not allowed.
    887. 

  887. !
    888. 

  888.    IF( iqsat == 1 ) THEN
    889. 

  889.      DO j=jts,jtf
    890. 

  890.      DO k=kts,ktf
    891. 

  891.      DO i=its,itf
    892. 

  892.        es = SVP1*EXP(SVP2*(t3d(i,k,j)-SVPT0)/(t3d(i,k,j)-SVP3))  * 10.0    ! mb
    893. 

  893.        qsat(i,k,j) = EP_2*es/(p3d(i,k,j)/100.0-es)
    894. 

  894.      ENDDO
    895. 

  895.      ENDDO
    896. 

  896.      ENDDO
    897. 

  897. 

  898.      IF ( wrf_dm_on_monitor()) THEN
    899. 

  899.        IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
    900. 

  900.          IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    901. 

  901.            i0 = (ite-its)/2+its
    902. 

  902.            j0 = (jte-jts)/2+jts 
    903. 

  903.            DO k = kts, kte
    904. 

  904.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    905. 

  905.                '    D0',id,' sample moisture values (g/kg) at i,k,j=', i0, k, j0, &
    906. 

  906.                ' qv3d=', qv3d(i0,k,j0)*1000.0, ' qsat=', qsat(i0,k,j0)*1000.0
    907. 

  907.              CALL wrf_message( TRIM(message) )
    908. 

  908.            ENDDO
    909. 

  909.          ENDIF
    910. 

  910.        ENDIF
    911. 

  911.      ENDIF
    912. 

  912.    ENDIF
    913. 

  913. !
    914. 

  914. ! Obs data density weighting.
    915. 

  915. !
    916. 

  916.    IF( idd == 1 ) THEN
    917. 

  917. 

  918.      IF( rinblw < 0.001 ) THEN
    919. 

  919.        IF ( wrf_dm_on_monitor()) THEN
    920. 

  920.          WRITE(message,'(a)') 'Error in rinblw, please specify a reasonable value ***'
    921. 

  921.          CALL wrf_message( TRIM(message) )
    922. 

  922.        ENDIF
    923. 

  923.        CALL wrf_error_fatal('In grid FDDA')
    924. 

  924.      ENDIF
    925. 

  925. 

  926.      rindx = rinblw*1000.0/dx
    927. 

  927.      m = -0.8*2.0/rindx
    928. 

  928. 

  929.      DO j=MAX(jts-2,jds),MIN(jte+2,jde-1)
    930. 

  930.      DO i=MAX(its-2,ids),MIN(ite+2,ide-1)
    931. 

  931.        IF( odis_ndg_old(i,j) < 0.5*rinblw ) THEN
    932. 

  932.          blw_old(i,j) = 1.0
    933. 

  933.        ELSE
    934. 

  934.          x = min( odis_ndg_old(i,j)*1000./dx, rindx )
    935. 

  935.          blw_old(i,j) = m * x + b
    936. 

  936.        ENDIF
    937. 

  937. 

  938.        IF( odis_ndg_new(i,j) < 0.5*rinblw ) THEN
    939. 

  939.          blw_new(i,j) = 1.0
    940. 

  940.        ELSE
    941. 

  941.          x = min( odis_ndg_new(i,j)*1000./dx, rindx )
    942. 

  942.          blw_new(i,j) = m * x + b
    943. 

  943.        ENDIF
    944. 

  944.      ENDDO
    945. 

  945.      ENDDO
    946. 

  946. 

  947. ! Smoother applies one point outside the tile, but one point in from boundaries
    948. 

  948.      CALL smther(blw_old, its-2,ite+2, jts-2,jte+2, 1, &
    949. 

  949.                  MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
    950. 

  950.      CALL smther(blw_new, its-2,ite+2, jts-2,jte+2, 1, &
    951. 

  951.                  MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
    952. 

  952. 

  953.      WHERE ( blw_old > 1.0)
    954. 

  954.         blw_old = 1.0
    955. 

  955.      END WHERE
    956. 

  956.      WHERE ( blw_new > 1.0)
    957. 

  957.         blw_new = 1.0
    958. 

  958.      END WHERE
    959. 

  959.      WHERE ( blw_old < 0.0)
    960. 

  960.         blw_old = 0.0
    961. 

  961.      END WHERE
    962. 

  962.      WHERE ( blw_new < 0.0)
    963. 

  963.         blw_new = 0.0
    964. 

  964.      END WHERE
    965. 

  965. 

  966.      IF ( wrf_dm_on_monitor()) THEN
    967. 

  967.        IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
    968. 

  968.          IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    969. 

  969.            i0 = (ite-its)/2+its
    970. 

  970.            j0 = (jte-jts)/2+jts
    971. 

  971.            WRITE(message,'(a,i1,a,2i4,4(a,f10.4))') &
    972. 

  972.              '    D0',id,' sample blw values at i,j=', i0, j0, &
    973. 

  973.              ' odis_ndg_old=', odis_ndg_old(i0,j0), ' km odis_ndg_new=', odis_ndg_new(i0,j0), &
    974. 

  974.              ' km  blw_old=', blw_old(i0,j0), ' blw_new=', blw_new(i0,j0)
    975. 

  975.            CALL wrf_message( TRIM(message) )
    976. 

  976.          ENDIF
    977. 

  977.        ENDIF
    978. 

  978.      ENDIF
    979. 

  979. 

  980.    ENDIF
    981. 

  981. !
    982. 

  982. ! TFAC for surface analysis nudging
    983. 

  983. !
    984. 

  984.    IF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min &
    985. 

  985.        .AND. dtramp_min > 0.0 .AND. if_ramping == 1 )  &
    986. 

  986.     coef = (xtime-xtime_old_sfc+analysis_interval_sfc)/(analysis_interval_sfc*1.0)
    987. 

  987. !  print*, 'coef =', xtime_old_sfc, xtime, xtime_new_sfc, coef
    988. 

  988. !              
    989. 

  989. ! Compute surface analysis nudging tendencies for u, v, t and q
    990. 

  990. !              
    991. 

  991.    DO j=jts,jtf
    992. 

  992.    DO k=kts,ktf
    993. 

  993.    DO i=itsu,itf
    994. 

  994.      IF( idd == 1 ) THEN
    995. 

  995.        blw = 0.5* (blw_old(i-1,j)+blw_old(i,j)) * ( 1.0 - coef ) &
    996. 

  996.            + 0.5* (blw_new(i-1,j)+blw_new(i,j)) * coef
    997. 

  997.      ELSE
    998. 

  998.        blw = 1.0
    999. 

  999.      ENDIF
    1000. 

  1000.      val_analysis = u10_ndg_old(i,j) *( 1.0 - coef ) + u10_ndg_new(i,j) * coef
    1001. 

  1001.      val_analysis = val_analysis * wndcor_u(i,j)
    1002. 

  1002.      RUNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,1)) * wzfac(k,1) * tfac * blw * &
    1003. 

  1003.                          ( val_analysis - u3d(i,1,j) )
    1004. 

  1004.    ENDDO
    1005. 

  1005.    ENDDO
    1006. 

  1006.    ENDDO
    1007. 

  1007. 

  1008.    DO j=jtsv,jtf
    1009. 

  1009.    DO k=kts,ktf 
    1010. 

  1010.    DO i=its,itf
    1011. 

  1011.      IF( idd == 1 ) THEN
    1012. 

  1012.        blw = 0.5* (blw_old(i,j-1)+blw_old(i,j)) * ( 1.0 - coef ) &
    1013. 

  1013.            + 0.5* (blw_new(i,j-1)+blw_new(i,j)) * coef
    1014. 

  1014.      ELSE
    1015. 

  1015.        blw = 1.0
    1016. 

  1016.      ENDIF
    1017. 

  1017.      val_analysis = v10_ndg_old(i,j) *( 1.0 - coef ) + v10_ndg_new(i,j) * coef
    1018. 

  1018.      val_analysis = val_analysis * wndcor_v(i,j)
    1019. 

  1019.      RVNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,2)) * wzfac(k,2) * tfac * blw * &
    1020. 

  1020.                        ( val_analysis - v3d(i,1,j) )
    1021. 

  1021.    ENDDO
    1022. 

  1022.    ENDDO
    1023. 

  1023.    ENDDO
    1024. 

  1024. 

  1025.    DO j=jts,jtf
    1026. 

  1026.    DO k=kts,ktf
    1027. 

  1027.    DO i=its,itf
    1028. 

  1028.      IF( idd == 1 ) THEN
    1029. 

  1029.        blw = blw_old(i,j) * ( 1.0 - coef ) + blw_new(i,j) * coef
    1030. 

  1030.      ELSE
    1031. 

  1031.        blw = 1.0
    1032. 

  1032.      ENDIF
    1033. 

  1033.      val_analysis = th2_ndg_old(i,j) *( 1.0 - coef ) + th2_ndg_new(i,j) * coef
    1034. 

  1034.      RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac * blw * &
    1035. 

  1035.                           ( val_analysis - th3d(i,1,j))
    1036. 

  1036. 

  1037.      val_analysis = q2_ndg_old(i,j) *( 1.0 - coef ) + q2_ndg_new(i,j) * coef
    1038. 

  1038.      IF( iqsat ==  1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
    1039. 

  1039.      RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac * blw * &
    1040. 

  1040.                           ( val_analysis - qv3d(i,k,j) )
    1041. 

  1041.    ENDDO
    1042. 

  1042.    ENDDO
    1043. 

  1043.    ENDDO
    1044. 

  1044. 

  1045.    END SUBROUTINE sfddagd
    1046. 

  1046. 

  1047. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1048. 

  1048. 

  1049.    SUBROUTINE fddagdinit(id,rundgdten,rvndgdten,rthndgdten,rqvndgdten,rmundgdten,&
    1050. 

  1050.                run_hours,  &
    1051. 

  1051.                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
    1052. 

  1052.                if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
    1053. 

  1053.                guv, gt, gq, if_ramping, dtramp_min, end_fdda_hour, &
    1054. 

  1054.                grid_sfdda, guv_sfc, gt_sfc, gq_sfc,                &
    1055. 

  1055.                       restart, allowed_to_read,                    &
    1056. 

  1056.                       ids, ide, jds, jde, kds, kde,                &
    1057. 

  1057.                       ims, ime, jms, jme, kms, kme,                &
    1058. 

  1058.                       its, ite, jts, jte, kts, kte                 )
    1059. 

  1059. !-------------------------------------------------------------------
    1060. 

  1060.    IMPLICIT NONE
    1061. 

  1061. !-------------------------------------------------------------------
    1062. 

  1062. !
    1063. 

  1063.    INTEGER , INTENT(IN)         ::  id
    1064. 

  1064.    LOGICAL, INTENT(IN)          ::  restart, allowed_to_read
    1065. 

  1065.    INTEGER, INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
    1066. 

  1066.                                     ims, ime, jms, jme, kms, kme, &
    1067. 

  1067.                                     its, ite, jts, jte, kts, kte
    1068. 

  1068.    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
    1069. 

  1069.                                                        rundgdten, &
    1070. 

  1070.                                                        rvndgdten, &
    1071. 

  1071.                                                       rthndgdten, &
    1072. 

  1072.                                                       rqvndgdten
    1073. 

  1073.    INTEGER,  INTENT(IN)   ::      run_hours
    1074. 

  1074.    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
    1075. 

  1075.                                   if_no_pbl_nudging_q, end_fdda_hour
    1076. 

  1076.    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_q
    1077. 

  1077.    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t,  k_zfac_q
    1078. 

  1078.    INTEGER,  INTENT(IN)   ::      if_ramping, grid_sfdda
    1079. 

  1079.    REAL,     INTENT(IN)   ::      dtramp_min
    1080. 

  1080.    REAL, INTENT(IN)       ::      guv, gt, gq
    1081. 

  1081.    REAL, INTENT(IN)       ::      guv_sfc, gt_sfc, gq_sfc
    1082. 

  1082.    REAL                   ::      actual_end_fdda_min
    1083. 

  1083. 

  1084.    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: rmundgdten
    1085. 

  1085.    INTEGER :: i, j, k
    1086. 

  1086. 

  1087.    LOGICAL , EXTERNAL     ::      wrf_dm_on_monitor
    1088. 

  1088. 

  1089.    CHARACTER (LEN=256) :: message
    1090. 

  1090. 

  1091.    IF ( wrf_dm_on_monitor() ) THEN  
    1092. 

  1092. 

  1093.      IF( guv > 0.0 ) THEN
    1094. 

  1094.        WRITE(message,'(a,i1,a,e12.4)') &
    1095. 

  1095.            'D0',id,' 3-D analysis nudging for wind is applied and Guv= ', guv
    1096. 

  1096.        CALL wrf_message(TRIM(message))
    1097. 

  1097.      ELSE IF( guv < 0.0 ) THEN
    1098. 

  1098.        CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
    1099. 

  1099.      ELSE 
    1100. 

  1100.        WRITE(message,'(a,i1,a,e12.4)') &
    1101. 

  1101.            'D0',id,' 3-D analysis nudging for wind is not applied and Guv= ', guv
    1102. 

  1102.        CALL wrf_message(TRIM(message))
    1103. 

  1103.      ENDIF
    1104. 

  1104. 

  1105.      IF( gt > 0.0 ) THEN
    1106. 

  1106.        WRITE(message,'(a,i1,a,e12.4)') &
    1107. 

  1107.            'D0',id,' 3-D analysis nudging for temperature is applied and Gt= ', gt
    1108. 

  1108.        CALL wrf_message(TRIM(message))
    1109. 

  1109.      ELSE IF( gt < 0.0 ) THEN
    1110. 

  1110.        CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
    1111. 

  1111.      ELSE 
    1112. 

  1112.        WRITE(message,'(a,i1,a,e12.4)') &
    1113. 

  1113.            'D0',id,' 3-D analysis nudging for temperature is not applied and Gt= ', gt
    1114. 

  1114.        CALL wrf_message(TRIM(message))
    1115. 

  1115.      ENDIF
    1116. 

  1116. 

  1117.      IF( gq > 0.0 ) THEN
    1118. 

  1118.        WRITE(message,'(a,i1,a,e12.4)') &
    1119. 

  1119.          'D0',id,' 3-D analysis nudging for water vapor mixing ratio is applied and Gq= ', gq
    1120. 

  1120.        CALL wrf_message(TRIM(message))
    1121. 

  1121.      ELSE IF( gq < 0.0 ) THEN
    1122. 

  1122.        CALL wrf_error_fatal('In grid FDDA, Gq must be positive.')
    1123. 

  1123.      ELSE
    1124. 

  1124.        WRITE(message,'(a,i1,a,e12.4)') &
    1125. 

  1125.          'D0',id,' 3-D analysis nudging for water vapor mixing ratio is not applied and Gq= ', gq
    1126. 

  1126.        CALL wrf_message(TRIM(message))
    1127. 

  1127.      ENDIF
    1128. 

  1128. 

  1129.      IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
    1130. 

  1130.         WRITE(message,'(a,i1,a)') &
    1131. 

  1131.            'D0',id,' 3-D analysis nudging for wind is turned off within the PBL.'
    1132. 

  1132.         CALL wrf_message(TRIM(message))
    1133. 

  1133.      ENDIF
    1134. 

  1134. 

  1135.      IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
    1136. 

  1136.         WRITE(message,'(a,i1,a)') &
    1137. 

  1137.            'D0',id,' 3-D analysis nudging for temperature is turned off within the PBL.'
    1138. 

  1138.         CALL wrf_message(TRIM(message))
    1139. 

  1139.      ENDIF
    1140. 

  1140. 

  1141.      IF( gq > 0.0 .AND. if_no_pbl_nudging_q == 1 ) THEN
    1142. 

  1142.         WRITE(message,'(a,i1,a)') &
    1143. 

  1143.          'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off within the PBL.'
    1144. 

  1144.         CALL wrf_message(TRIM(message))
    1145. 

  1145.      ENDIF
    1146. 

  1146. 

  1147.      IF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
    1148. 

  1148.         WRITE(message,'(a,i1,a,i3)') &
    1149. 

  1149.            'D0',id,' 3-D analysis nudging for wind is turned off below layer', k_zfac_uv
    1150. 

  1150.         CALL wrf_message(TRIM(message))
    1151. 

  1151.      ENDIF
    1152. 

  1152. 

  1153.      IF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
    1154. 

  1154.         WRITE(message,'(a,i1,a,i3)') &
    1155. 

  1155.            'D0',id,' 3-D analysis nudging for temperature is turned off below layer', k_zfac_t
    1156. 

  1156.         CALL wrf_message(TRIM(message))
    1157. 

  1157.      ENDIF
    1158. 

  1158. 

  1159.      IF( gq > 0.0 .AND. if_zfac_q == 1 ) THEN
    1160. 

  1160.         WRITE(message,'(a,i1,a,i3)') &
    1161. 

  1161.           'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off below layer', &
    1162. 

  1162.            k_zfac_q
    1163. 

  1163.         CALL wrf_message(TRIM(message))
    1164. 

  1164.      ENDIF
    1165. 

  1165. 

  1166.      IF( grid_sfdda ==1 ) THEN
    1167. 

  1167.        IF( guv_sfc > 0.0 ) THEN
    1168. 

  1168.          WRITE(message,'(a,i1,a,e12.4)') &
    1169. 

  1169.              'D0',id,' surface analysis nudging for wind is applied and Guv_sfc= ', guv_sfc
    1170. 

  1170.          CALL wrf_message(TRIM(message))
    1171. 

  1171.        ELSE IF( guv_sfc < 0.0 ) THEN
    1172. 

  1172.          CALL wrf_error_fatal('In grid FDDA, Guv_sfc must be positive.')
    1173. 

  1173.        ELSE
    1174. 

  1174.          WRITE(message,'(a,i1,a,e12.4)') &
    1175. 

  1175.              'D0',id,' surface analysis nudging for wind is not applied and Guv_sfc= ', guv_sfc
    1176. 

  1176.          CALL wrf_message(TRIM(message))
    1177. 

  1177.        ENDIF
    1178. 

  1178. 

  1179.        IF( gt_sfc > 0.0 ) THEN
    1180. 

  1180.          WRITE(message,'(a,i1,a,e12.4)') &
    1181. 

  1181.              'D0',id,' surface analysis nudging for temperature is applied and Gt_sfc= ', gt_sfc
    1182. 

  1182.          CALL wrf_message(TRIM(message))
    1183. 

  1183.        ELSE IF( gt_sfc < 0.0 ) THEN
    1184. 

  1184.          CALL wrf_error_fatal('In grid FDDA, Gt_sfc must be positive.')
    1185. 

  1185.        ELSE
    1186. 

  1186.          WRITE(message,'(a,i1,a,e12.4)') &
    1187. 

  1187.              'D0',id,' surafce analysis nudging for temperature is not applied and Gt_sfc= ', gt_sfc
    1188. 

  1188.          CALL wrf_message(TRIM(message))
    1189. 

  1189.        ENDIF
    1190. 

  1190. 

  1191.        IF( gq_sfc > 0.0 ) THEN
    1192. 

  1192.          WRITE(message,'(a,i1,a,e12.4)') &
    1193. 

  1193.            'D0',id,' surface analysis nudging for water vapor mixing ratio is applied and Gq_sfc= ', gq_sfc
    1194. 

  1194.          CALL wrf_message(TRIM(message))
    1195. 

  1195.        ELSE IF( gq_sfc < 0.0 ) THEN
    1196. 

  1196.          CALL wrf_error_fatal('In grid FDDA, Gq_sfc must be positive.')
    1197. 

  1197.        ELSE
    1198. 

  1198.          WRITE(message,'(a,i1,a,e12.4)') &
    1199. 

  1199.            'D0',id,' surface analysis nudging for water vapor mixing ratio is not applied and Gq_sfc= ', gq_sfc
    1200. 

  1200.          CALL wrf_message(TRIM(message))
    1201. 

  1201.        ENDIF
    1202. 

  1202. 

  1203.      ENDIF
    1204. 

  1204. 

  1205.      IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
    1206. 

  1206.        IF( dtramp_min <= 0.0 ) THEN
    1207. 

  1207.          actual_end_fdda_min = end_fdda_hour*60.0
    1208. 

  1208.        ELSE
    1209. 

  1209.          actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    1210. 

  1210.        ENDIF
    1211. 

  1211. 

  1212.        IF( actual_end_fdda_min <= run_hours*60. ) THEN
    1213. 

  1213.           WRITE(message,'(a,i1,a)') &
    1214. 

  1214.             'D0',id,' analysis nudging is ramped down near the end of the nudging period,'
    1215. 

  1215.           CALL wrf_message(TRIM(message))
    1216. 

  1216. 

  1217.           WRITE(message,'(a,f6.2,a,f6.2,a)') &
    1218. 

  1218.              '      starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0, &
    1219. 

  1219.              'h, ending at ', actual_end_fdda_min/60.0,'h.'
    1220. 

  1220.           CALL wrf_message(TRIM(message))
    1221. 

  1221.        ENDIF
    1222. 

  1222.      ENDIF
    1223. 

  1223. 

  1224.    ENDIF
    1225. 

  1225. 

  1226.    IF(.not.restart) THEN
    1227. 

  1227.      DO j = jts,jte
    1228. 

  1228.      DO k = kts,kte
    1229. 

  1229.      DO i = its,ite
    1230. 

  1230.         rundgdten(i,k,j) = 0.
    1231. 

  1231.         rvndgdten(i,k,j) = 0.
    1232. 

  1232.         rthndgdten(i,k,j) = 0.
    1233. 

  1233.         rqvndgdten(i,k,j) = 0.
    1234. 

  1234.         if(k.eq.kts) rmundgdten(i,j) = 0.
    1235. 

  1235.      ENDDO
    1236. 

  1236.      ENDDO
    1237. 

  1237.      ENDDO
    1238. 

  1238.    ENDIF
    1239. 

  1239. 

  1240.    END SUBROUTINE fddagdinit
    1241. 

  1241. 

  1242. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1243. 

  1243.   SUBROUTINE smther(slab, idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd)
    1244. 

  1244. !
    1245. 

  1245. !    PURPOSE: SPATIALLY SMOOTH DATA IN SLAB TO DAMPEN SHORT
    1246. 

  1246. !             WAVELENGTH COMPONENTS.
    1247. 

  1247. !
    1248. 

  1248. !    Implemented based on the same smoothing subroutine used in MM5, with modifications to
    1249. 

  1249. !    remove the extra loop that causes unneccessary desmoothing.  Aijun Deng (Penn State),
    1250. 

  1250. !    December 2008 
    1251. 

  1251. !
    1252. 

  1252.     IMPLICIT NONE
    1253. 

  1253. 

  1254.     INTEGER :: idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd
    1255. 

  1255.     INTEGER :: i, j, k, kk
    1256. 

  1256.     REAL    :: avg
    1257. 

  1257.     REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB
    1258. 

  1258.     REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB_ORIG
    1259. 

  1259.     REAL, DIMENSION(2)                            :: XNU
    1260. 

  1260. 

  1261.     IF(NPASS.EQ.0)RETURN
    1262. 

  1262.     XNU(1)=0.50
    1263. 

  1263.     XNU(2)=-0.52
    1264. 

  1264.     DO K=1,NPASS
    1265. 

  1265.       KK = 2 - MOD(K,2)
    1266. 

  1266. 

  1267.       DO J=JDIMST,JDIMND
    1268. 

  1268.         DO I=IDIMST,IDIMND
    1269. 

  1269.            SLAB_ORIG(I,J) = SLAB(I,J)
    1270. 

  1270.         END DO
    1271. 

  1271.       END DO
    1272. 

  1272. 

  1273.       DO J=JST,JND
    1274. 

  1274.         DO I=IST,IND
    1275. 

  1275.           AVG = ( SLAB_ORIG(I+1,J  ) + &
    1276. 

  1276.                   SLAB_ORIG(I-1,J  ) + &
    1277. 

  1277.                   SLAB_ORIG(I  ,J+1) + &
    1278. 

  1278.                   SLAB_ORIG(I  ,J-1) ) * 0.25
    1279. 

  1279.           SLAB(I,J)=SLAB_ORIG(I,J)+XNU(KK)*(AVG - SLAB_ORIG(I,J))
    1280. 

  1280.         ENDDO
    1281. 

  1281.       ENDDO
    1282. 

  1282. 

  1283.     ENDDO
    1284. 

  1284. 

  1285.   END SUBROUTINE smther
    1286. 

  1286. 

  1287. !-------------------------------------------------------------------
    1288. 

  1288. END MODULE module_fdda_psufddagd
    1289. 

  1289.