PageRenderTime 52ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_fdda_spnudging.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1087 lines | 658 code | 230 blank | 199 comment | 18 complexity | 325b380e0394c020783b2b679289484d MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !wrf:model_layer:physics
    2. 

  2. !
    3. 

  3. ! Code contributed by Gonzalo Miguez-Macho, Univ. de Santiago de Compostela, Galicia, Spain
    4. 

  4. !
    5. 

  5. MODULE module_fdda_spnudging
    6. 

  6. 

  7. #ifdef DM_PARALLEL
    8. 

  8.   USE module_dm , ONLY : ntasks_x, ntasks_y, local_communicator_x, local_communicator_y, data_order_xzy
    9. 

  9. #endif
    10. 

  10.   USE module_wrf_error , ONLY : wrf_err_message
    11. 

  11. 

  12. CONTAINS
    13. 

  13. !
    14. 

  14. !-------------------------------------------------------------------
    15. 

  15. !
    16. 

  16.    SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fdda_hour, &
    17. 

  17.                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph,&
    18. 

  18.                if_zfac_uv, k_zfac_uv, dk_zfac_uv, & 
    19. 

  19.                if_zfac_t, k_zfac_t, dk_zfac_t, &
    20. 

  20.                if_zfac_ph, k_zfac_ph, dk_zfac_ph, &
    21. 

  21.                guv, gt, gph, if_ramping, dtramp_min, xwavenum, ywavenum, &
    22. 

  22.                u3d,v3d,th3d,ph3d,                 &
    23. 

  23.                u_ndg_old,v_ndg_old,t_ndg_old,ph_ndg_old,       &
    24. 

  24.                u_ndg_new,v_ndg_new,t_ndg_new,ph_ndg_new,       &
    25. 

  25.                RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RPHNDGDTEN,&
    26. 

  26.                pblh, ht, z, z_at_w,                             &
    27. 

  27.                ids,ide, jds,jde, kds,kde,                           &
    28. 

  28.                ims,ime, jms,jme, kms,kme,                           &
    29. 

  29.                i_start,i_end, j_start,j_end, kts,kte, num_tiles, &
    30. 

  30.                ips,ipe,jps,jpe,kps,kpe,                   &
    31. 

  31.                imsx,imex,jmsx,jmex,kmsx,kmex,                    &
    32. 

  32.                ipsx,ipex,jpsx,jpex,kpsx,kpex,                   &
    33. 

  33.                imsy,imey,jmsy,jmey,kmsy,kmey,                    &
    34. 

  34.                ipsy,ipey,jpsy,jpey,kpsy,kpey                   )
    35. 

  35. 

  36. 

  37. 

  38.    USE module_state_description
    39. 

  39.    USE module_domain, ONLY : domain
    40. 

  40. 

  41. !-------------------------------------------------------------------
    42. 

  42.    implicit none
    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. !---ph3d        3d perturbation geopotential
    48. 

  48. !-- rundgdten   staggered u tendency due to
    49. 

  49. !               spectral nudging (m/s/s)
    50. 

  50. !-- rvndgdten   staggered v tendency due to
    51. 

  51. !               spectral nudging (m/s/s)
    52. 

  52. !-- rthndgdten  theta tendency due to
    53. 

  53. !               spectral nudging (K/s)
    54. 

  54. !-- phndgdten  ph tendency due to
    55. 

  55. !               spectral nudging (m2/s2/s)
    56. 

  56. !-- ids         start index for i in domain
    57. 

  57. !-- ide         end index for i in domain
    58. 

  58. !-- jds         start index for j in domain
    59. 

  59. !-- jde         end index for j in domain
    60. 

  60. !-- kds         start index for k in domain
    61. 

  61. !-- kde         end index for k in domain
    62. 

  62. !-- ims         start index for i in memory
    63. 

  63. !-- ime         end index for i in memory
    64. 

  64. !-- jms         start index for j in memory
    65. 

  65. !-- jme         end index for j in memory
    66. 

  66. !-- kms         start index for k in memory
    67. 

  67. !-- kme         end index for k in memory
    68. 

  68. !-- its         start index for i in tile
    69. 

  69. !-- ite         end index for i in tile
    70. 

  70. !-- jts         start index for j in tile
    71. 

  71. !-- jte         end index for j in tile
    72. 

  72. !-- kts         start index for k in tile
    73. 

  73. !-- kte         end index for k in tile
    74. 

  74. !-------------------------------------------------------------------
    75. 

  75. !
    76. 

  76.    TYPE(domain) , TARGET          :: grid
    77. 

  77. 

  78.    INTEGER,  INTENT(IN)   ::      itimestep, analysis_interval, end_fdda_hour
    79. 

  79. 

  80.    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
    81. 

  81.                                   if_no_pbl_nudging_ph
    82. 

  82.    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_ph
    83. 

  83.    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t, k_zfac_ph
    84. 

  84.    INTEGER,  INTENT(IN)   ::      dk_zfac_uv,  dk_zfac_t, dk_zfac_ph
    85. 

  85.    INTEGER,  INTENT(IN)   ::      if_ramping
    86. 

  86.    INTEGER,  INTENT(IN)   ::      xwavenum,ywavenum
    87. 

  87. 

  88.    INTEGER , INTENT(IN)   ::      id
    89. 

  89.    REAL,     INTENT(IN)   ::      DT, xtime, dtramp_min
    90. 

  90. 

  91.    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
    92. 

  92.                                   ims,ime, jms,jme, kms,kme, &
    93. 

  93.                                   kts,kte, num_tiles,        &
    94. 

  94.                                   ips,ipe,jps,jpe,kps,kpe,   &
    95. 

  95.                                   imsx,imex,jmsx,jmex,kmsx,kmex,   &
    96. 

  96.                                   ipsx,ipex,jpsx,jpex,kpsx,kpex,   &
    97. 

  97.                                   imsy,imey,jmsy,jmey,kmsy,kmey,   &
    98. 

  98.                                   ipsy,ipey,jpsy,jpey,kpsy,kpey 
    99. 

  99. 

  100.    INTEGER, DIMENSION(num_tiles), INTENT(IN) ::                   &
    101. 

  101.   &                                    i_start,i_end,j_start,j_end
    102. 

  102.  
    103. 

  103.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    104. 

  104.              INTENT(IN)   ::                  ph3d, &
    105. 

  105.                                               th3d, &
    106. 

  106.                                                  z, &
    107. 

  107.                                             z_at_w
    108. 

  108. 

  109.    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
    110. 

  110.              INTENT(INOUT)   ::           rundgdten, &
    111. 

  111.                                           rvndgdten, &
    112. 

  112.                                          rthndgdten, &
    113. 

  113.                                          rphndgdten
    114. 

  114. 

  115. 

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

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

  118.                                           v_ndg_old, &
    119. 

  119.                                           t_ndg_old, &
    120. 

  120.                                           ph_ndg_old, &
    121. 

  121.                                           u_ndg_new, &
    122. 

  122.                                           v_ndg_new, &
    123. 

  123.                                           t_ndg_new, &
    124. 

  124.                                           ph_ndg_new
    125. 

  125. 

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

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

  128.                                                 v3d
    129. 

  129. 

  130.    REAL,  DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
    131. 

  131.                                                          ht
    132. 

  132. 

  133.    REAL, INTENT(IN)    :: guv, gt ,gph
    134. 

  134. 

  135.    INTEGER             :: its,ite, jts,jte,ij
    136. 

  136.    INTEGER             :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
    137. 

  137.    REAL                :: xtime_old, xtime_new, coef, val_analysis
    138. 

  138.    INTEGER             :: kpbl, dbg_level
    139. 

  139. 

  140.    REAL                :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
    141. 

  141.    REAL, DIMENSION( ips:ipe, kps:kpe, jps:jpe, 4 ) :: wpbl  ! 1: u, 2: v, 3: t, 4: ph
    142. 

  142.    REAL, DIMENSION( kps:kpe, 4 )                   :: wzfac ! 1: u, 2: v, 3: t, 4: ph
    143. 

  143. 

  144.    LOGICAL , EXTERNAL  :: wrf_dm_on_monitor
    145. 

  145. 

  146.    CHARACTER (LEN=256) :: message
    147. 

  147. 

  148. #if  ! ( NMM_CORE == 1 )
    149. 

  149. 

  150.       DO ij = 1 , num_tiles
    151. 

  151.           its=i_start(ij)
    152. 

  152.           ite=i_end(ij)
    153. 

  153.           jts=j_start(ij)
    154. 

  154.           jte=j_end(ij)
    155. 

  155.       ENDDO
    156. 

  156. !GMM default values for vertical coefficients, set here
    157. 

  157.    wpbl(:,:,:,:) = 1.0
    158. 

  158.    wzfac(:,:) = 1.0
    159. 

  159. 

  160.  !$OMP PARALLEL DO   &
    161. 

  161.  !$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation.
    162. 

  162.  DO ij = 1 , num_tiles
    163. 

  163.      DO j = jts, jte
    164. 

  164.      DO k = kts, kte
    165. 

  165.      DO i = its, ite
    166. 

  166.        RUNDGDTEN(i,k,j) = 0.0
    167. 

  167.        RVNDGDTEN(i,k,j) = 0.0
    168. 

  168.        RTHNDGDTEN(i,k,j) = 0.0
    169. 

  169.        RPHNDGDTEN(i,k,j) = 0.0
    170. 

  170.      ENDDO
    171. 

  171.      ENDDO
    172. 

  172.      ENDDO
    173. 

  173.  ENDDO
    174. 

  174.  !$OMP END PARALLEL DO
    175. 

  175. 

  176.  actual_end_fdda_min = end_fdda_hour*60.0
    177. 

  177.  IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
    178. 

  178.        actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    179. 

  179.  IF( xtime > actual_end_fdda_min ) RETURN
    180. 

  180. 

  181.  !$OMP PARALLEL DO   &
    182. 

  182.  !$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation. 
    183. 

  183.  !GMM Transpose and filtering needs to be done on patch
    184. 

  184.  
    185. 

  185.  DO ij = 1 , num_tiles
    186. 

  186. 

  187. !  actual_end_fdda_min = end_fdda_hour*60.0
    188. 

  188. !  IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
    189. 

  189. !      actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    190. 

  190. 

  191.    xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
    192. 

  192.    xtime_new = xtime_old + analysis_interval
    193. 

  193.    coef = (xtime-xtime_old)/(xtime_new-xtime_old)
    194. 

  194. 

  195.    IF ( wrf_dm_on_monitor()) THEN
    196. 

  196. 

  197.      CALL get_wrf_debug_level( dbg_level )
    198. 

  198. 

  199.      IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
    200. 

  200. 

  201.        IF( xtime < end_fdda_hour*60.0 ) THEN
    202. 

  202.          WRITE(message,FMT='(a,i2.2,a,f15.3,a)') &
    203. 

  203.           'D',id,' Spectral nudging read in new data at time = ', xtime, ' min.'
    204. 

  204.          CALL wrf_message( TRIM(message) )
    205. 

  205.          WRITE(message,FMT='(a,i2.2,a,2(f15.3,1x),a)') &
    206. 

  206.           'D',id,' Spectral nudging bracketing times = ', xtime_old, xtime_new, ' min.'
    207. 

  207.          CALL wrf_message( TRIM(message) )
    208. 

  208.        ENDIF
    209. 

  209. 

  210.        actual_end_fdda_min = end_fdda_hour*60.0
    211. 

  211.        IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
    212. 

  212.            actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    213. 

  213. 

  214.        IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
    215. 

  215. !        Find the mid point of the tile and print out the sample values
    216. 

  216.          i0 = (ite-its)/2+its
    217. 

  217.          j0 = (jte-jts)/2+jts 
    218. 

  218. 

  219.          IF( guv > 0.0 ) THEN
    220. 

  220.            DO k = kts, kte
    221. 

  221.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    222. 

  222.                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
    223. 

  223.                ' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
    224. 

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

  225.            ENDDO
    226. 

  226.            DO k = kts, kte
    227. 

  227.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    228. 

  228.                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
    229. 

  229.                ' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
    230. 

  230.              CALL wrf_message( TRIM(message) )
    231. 

  231.            ENDDO
    232. 

  232.          ENDIF
    233. 

  233. 

  234.          IF( gt > 0.0 ) THEN
    235. 

  235.            DO k = kts, kte
    236. 

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

  237.                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
    238. 

  238.                ' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
    239. 

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

  240.            ENDDO
    241. 

  241.          ENDIF
    242. 

  242. 

  243.          IF( gph > 0.0 ) THEN
    244. 

  244.            DO k = kts, kte
    245. 

  245.              WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
    246. 

  246.                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
    247. 

  247.                ' ph_ndg_old=', ph_ndg_old(i0,k,j0), ' ph_ndg_new=', ph_ndg_new(i0,k,j0)
    248. 

  248.              CALL wrf_message( TRIM(message) )
    249. 

  249.            ENDDO
    250. 

  250.          ENDIF
    251. 

  251. 

  252.        ENDIF
    253. 

  253.      ENDIF
    254. 

  254.    ENDIF
    255. 

  255. 

  256.    jtsv=MAX0(jts,jds+1)
    257. 

  257.    itsu=MAX0(its,ids+1)
    258. 

  258. 

  259.    jtf=MIN0(jte,jde-1)
    260. 

  260.    ktf=MIN0(kte,kde-1)
    261. 

  261.    itf=MIN0(ite,ide-1)
    262. 

  262. !
    263. 

  263. ! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t, 
    264. 

  264. ! if_no_pbl_nudging_ph swithes) are set to 1, compute the weighting function, wpbl(:,k,:,:), 
    265. 

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

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

  267. !
    268. 

  268. !GMM If those are set to zero, then check if the user-defined namelist switches (if_zfac_uv, if_zfac_t,
    269. 

  269. ! if_zfac_ph swithes) are set to 1, compute the weighting function, wzfac(k,:),
    270. 

  270. ! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_ph) 
    271. 

  271. ! below which analysis nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x).
    272. 

  272. !  The strength of nudging increases linearly from k_zfac to kzfac + dk_zfac 
    273. 

  273. ! (dk_zfac_uv, dk_zfac_t and kd_zfac_ph are also set in the namelist, default value is 1).
    274. 

  274. !If all switches are set to zero, wzfac = 1 (default value).
    275. 

  275. !
    276. 

  276. 

  277.    IF( if_no_pbl_nudging_uv == 1 ) THEN
    278. 

  278. 

  279.      DO j=jts,jtf 
    280. 

  280.      DO i=itsu,itf
    281. 

  281. 

  282.        kpbl = 1
    283. 

  283.        zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
    284. 

  284. 

  285.        loop_ku: DO k=kts,ktf 
    286. 

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

  287.          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) )
    288. 

  288.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    289. 

  289.            kpbl = k
    290. 

  290.            EXIT loop_ku
    291. 

  291.          ENDIF
    292. 

  292.        ENDDO loop_ku
    293. 

  293. 

  294.        DO k=kts,ktf
    295. 

  295.           wpbl(i, k, j, 1) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
    296. 

  296.        ENDDO
    297. 

  297. 

  298.      ENDDO
    299. 

  299.      ENDDO
    300. 

  300. 

  301.      DO i=its,itf
    302. 

  302.      DO j=jtsv,jtf
    303. 

  303. 

  304.        kpbl = 1
    305. 

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

  306. 

  307.        loop_kv: DO k=kts,ktf
    308. 

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

  309.          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) )
    310. 

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

  311.            kpbl = k
    312. 

  312.            EXIT loop_kv
    313. 

  313.          ENDIF
    314. 

  314.        ENDDO loop_kv
    315. 

  315. 

  316.        DO k=kts,ktf
    317. 

  317.           wpbl(i, k, j, 2) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
    318. 

  318.        ENDDO
    319. 

  319. 

  320.      ENDDO
    321. 

  321.      ENDDO
    322. 

  322. 

  323.    ELSEIF( if_zfac_uv == 1 ) THEN
    324. 

  324. 

  325.      DO j=jts,jte
    326. 

  326.      DO k=kts,ktf
    327. 

  327.      DO i=its,ite
    328. 

  328.           wpbl(i, k, j, 1) = max ( min ( float(k-k_zfac_uv) / float(dk_zfac_uv) , 1. ) , 0.)
    329. 

  329.           wpbl(i, k, j, 2) = wpbl(i, k, j, 1)
    330. 

  330.      ENDDO
    331. 

  331.      ENDDO
    332. 

  332.      ENDDO
    333. 

  333. 

  334.    ENDIF
    335. 

  335. 

  336. 

  337.    IF( if_no_pbl_nudging_t == 1 ) THEN
    338. 

  338.    
    339. 

  339.      DO j=jts,jtf
    340. 

  340.      DO i=its,itf
    341. 

  341. 

  342.        kpbl = 1
    343. 

  343.        zpbl = pblh(i,j)
    344. 

  344.         
    345. 

  345.        loop_kt: DO k=kts,ktf
    346. 

  346.          zagl_bot = z_at_w(i,k,  j)-ht(i,j)
    347. 

  347.          zagl_top = z_at_w(i,k+1,j)-ht(i,j)
    348. 

  348.          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
    349. 

  349.            kpbl = k
    350. 

  350.            EXIT loop_kt
    351. 

  351.          ENDIF
    352. 

  352.        ENDDO loop_kt
    353. 

  353. 

  354.        DO k=kts,ktf
    355. 

  355.           wpbl(i, k, j, 3) = max ( min ( float(k-kpbl) / float(dk_zfac_t) , 1. ) , 0.)
    356. 

  356.        ENDDO 
    357. 

  357.         
    358. 

  358.      ENDDO
    359. 

  359.      ENDDO
    360. 

  360. 

  361.    ELSEIF( if_zfac_t == 1 ) THEN
    362. 

  362. 

  363.      DO j=jts,jtf
    364. 

  364.      DO k=kts,ktf
    365. 

  365.      DO i=its,itf
    366. 

  366.           wpbl(i, k, j, 3) = max ( min ( float(k-k_zfac_t) / float(dk_zfac_t) , 1. ) , 0.)
    367. 

  367.      ENDDO
    368. 

  368.      ENDDO
    369. 

  369.      ENDDO
    370. 

  370. 

  371.    ENDIF
    372. 

  372. 

  373. 

  374.    IF( if_no_pbl_nudging_ph == 1 ) THEN
    375. 

  375.    
    376. 

  376.      DO j=jts,jtf
    377. 

  377.      DO i=its,itf
    378. 

  378. 

  379.        kpbl = 1
    380. 

  380.        zpbl = pblh(i,j)
    381. 

  381.           
    382. 

  382.        loop_kph: DO k=kts,kte
    383. 

  383.          zagl = z_at_w(i,k,  j)-ht(i,j)
    384. 

  384.          IF( zpbl >= zagl) THEN
    385. 

  385.            kpbl = k
    386. 

  386.            EXIT loop_kph
    387. 

  387.          ENDIF
    388. 

  388.        ENDDO loop_kph
    389. 

  389. 

  390.        DO k=kts,kte
    391. 

  391.           wpbl(i, k, j, 4) = max ( min ( float(k-kpbl) / float(dk_zfac_ph) , 1. ) , 0.)
    392. 

  392.        ENDDO 
    393. 

  393.             
    394. 

  394.      ENDDO  
    395. 

  395.      ENDDO
    396. 

  396.         
    397. 

  397.    ELSEIF( if_zfac_ph == 1 ) THEN
    398. 

  398. 

  399.      DO j=jts,jtf
    400. 

  400.      DO k=kts,kte
    401. 

  401.      DO i=its,itf
    402. 

  402.           wpbl(i, k, j, 4) = max ( min ( float(k-k_zfac_ph) / float(dk_zfac_ph) , 1. ) , 0.)
    403. 

  403.      ENDDO
    404. 

  404.      ENDDO
    405. 

  405.      ENDDO
    406. 

  406. 

  407.    ENDIF
    408. 

  408. 

  409. 

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

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

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

  413. !
    414. 

  414. ! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at 
    415. 

  415. ! end_fdda_hour-ABS(dtramp_min).  
    416. 

  416. !
    417. 

  417. ! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at 
    418. 

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

  419. ! the obs tendency saved from the previous FDDA wondow.  More specifically for extrapolation, 
    420. 

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

  421. !
    422. 

  422.    tfac = 1.0
    423. 

  423. 

  424.    IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
    425. 

  425.  
    426. 

  426.      IF( dtramp_min <= 0.0 ) THEN
    427. 

  427.        actual_end_fdda_min = end_fdda_hour*60.0
    428. 

  428.      ELSE
    429. 

  429.        actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
    430. 

  430.      ENDIF
    431. 

  431. 

  432.      IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN 
    433. 

  433.        tfac = 1.0
    434. 

  434.      ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
    435. 

  435.        tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
    436. 

  436.        IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/analysis_interval
    437. 

  437.      ELSE                                                     
    438. 

  438.        tfac = 0.0
    439. 

  439.      ENDIF
    440. 

  440. 

  441.    ENDIF                                                  
    442. 

  442. 

  443.    ENDDO
    444. 

  444.    !$OMP END PARALLEL DO
    445. 

  445. 

  446. !
    447. 

  447. ! Compute 3-D nudging tendencies for u, v
    448. 

  448. !
    449. 

  449. !
    450. 

  450. !GMM Fist calculate differences between model variables and analysis values,
    451. 

  451. !then filter in the x and y direction all wave numbers higher than
    452. 

  452. ! xwavenum and ywavenum, as specified in the namelist.
    453. 

  453. !If either xwavenum or ywavenum are not specified,
    454. 

  454. ! default values are zero, and spectral nudging is turned off
    455. 

  455. !Then use the filtered differences to calculate the spectral nudging tendencies
    456. 

  456. 

  457. IF(guv > 0. ) then
    458. 

  458. 

  459.  !$OMP PARALLEL DO   &
    460. 

  460.  !$OMP PRIVATE ( ij, i,j,k )
    461. 

  461.  DO ij = 1 , num_tiles
    462. 

  462. 

  463.    DO j=jts,jte
    464. 

  464.    DO k=kts,ktf
    465. 

  465.    DO i=its,ite
    466. 

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

  467.      
    468. 

  468.      grid%dif_analysis(i,k,j) = val_analysis - u3d(i,k,j)
    469. 

  469. 

  470.    ENDDO
    471. 

  471.    ENDDO
    472. 

  472.    ENDDO
    473. 

  473.  
    474. 

  474.  ENDDO
    475. 

  475.  !$OMP END PARALLEL DO
    476. 

  476. 

  477. !Filter
    478. 

  478. 

  479. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    480. 

  480. # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
    481. 

  481. 

  482.      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
    483. 

  483.                                 ids, ide, jds, jde, kds, kde,         &
    484. 

  484.                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
    485. 

  485.                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) ) 
    486. 

  486. 

  487. # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
    488. 

  488. 

  489.      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
    490. 

  490.                                 ids, ide, jds, jde, kds, kde,         &
    491. 

  491.                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
    492. 

  492.                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
    493. 

  493. 

  494. # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
    495. 

  495. 

  496. 

  497. #else
    498. 

  498.      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
    499. 

  499.                                 ids, ide, jds, jde, kds, kde,       &
    500. 

  500.                                 ims, ime, jms, jme, kms, kme,       &
    501. 

  501.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    502. 

  502. 

  503.      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
    504. 

  504.                                 ids, ide, jds, jde, kds, kde,       &
    505. 

  505.                                 ims, ime, jms, jme, kms, kme,       &
    506. 

  506.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    507. 

  507. #endif
    508. 

  508. 

  509. ! Calculate tendency
    510. 

  510. 

  511.  !$OMP PARALLEL DO   &
    512. 

  512.  !$OMP PRIVATE ( ij, i,j,k )
    513. 

  513.  DO ij = 1 , num_tiles
    514. 

  514. 

  515.    DO j=jts,jte
    516. 

  516.    DO k=kts,ktf
    517. 

  517.    DO i=its,ite
    518. 

  518.      RUNDGDTEN(i,k,j) = guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
    519. 

  519.                         grid%dif_analysis(i,k,j)
    520. 

  520.    ENDDO
    521. 

  521.    ENDDO
    522. 

  522.    ENDDO
    523. 

  523. 

  524. 

  525. !
    526. 

  526. ! Now V component
    527. 

  527. !
    528. 

  528. 

  529.    DO j=jts,jte
    530. 

  530.    DO k=kts,ktf
    531. 

  531.    DO i=its,ite
    532. 

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

  533. 

  534.      grid%dif_analysis(i,k,j) = val_analysis - v3d(i,k,j)
    535. 

  535. 

  536.    ENDDO
    537. 

  537.    ENDDO
    538. 

  538.    ENDDO
    539. 

  539. 

  540.  ENDDO
    541. 

  541.  !$OMP END PARALLEL DO
    542. 

  542. 

  543. !
    544. 

  544. ! Filter
    545. 

  545. !
    546. 

  546. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    547. 

  547. # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
    548. 

  548.      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,     &
    549. 

  549.                                 ids, ide, jds, jde, kds, kde,         &
    550. 

  550.                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
    551. 

  551.                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
    552. 

  552. 

  553. # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
    554. 

  554.      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,   &
    555. 

  555.                                 ids, ide, jds, jde, kds, kde-1,         &
    556. 

  556.                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
    557. 

  557.                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
    558. 

  558. 

  559. # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
    560. 

  560. 

  561. 

  562. #else
    563. 

  563.      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,    &
    564. 

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

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

  566.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    567. 

  567. 

  568.      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,    &
    569. 

  569.                                 ids, ide, jds, jde, kds, kde,       &
    570. 

  570.                                 ims, ime, jms, jme, kms, kme,       &
    571. 

  571.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    572. 

  572. #endif
    573. 

  573. 

  574. ! Calculate tendency
    575. 

  575. 

  576.  !$OMP PARALLEL DO   &
    577. 

  577.  !$OMP PRIVATE ( ij, i,j,k )
    578. 

  578.  DO ij = 1 , num_tiles
    579. 

  579. 

  580.    DO j=jts,jte
    581. 

  581.    DO k=kts,ktf
    582. 

  582.    DO i=its,ite
    583. 

  583.      RVNDGDTEN(i,k,j) = guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
    584. 

  584.                         grid%dif_analysis(i,k,j)
    585. 

  585.    ENDDO
    586. 

  586.    ENDDO
    587. 

  587.    ENDDO
    588. 

  588. 

  589.  ENDDO
    590. 

  590.  !$OMP END PARALLEL DO
    591. 

  591. 

  592. ENDIF
    593. 

  593. 

  594. IF(gt > 0. ) then
    595. 

  595. 

  596.  !$OMP PARALLEL DO   &
    597. 

  597.  !$OMP PRIVATE ( ij, i,j,k )
    598. 

  598.  DO ij = 1 , num_tiles
    599. 

  599. 

  600.    DO j=jts,jte
    601. 

  601.    DO k=kts,ktf
    602. 

  602.    DO i=its,ite
    603. 

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

  604. 

  605.      grid%dif_analysis(i,k,j) = val_analysis - th3d(i,k,j) + 300.
    606. 

  606. 

  607.    ENDDO
    608. 

  608.    ENDDO
    609. 

  609.    ENDDO
    610. 

  610. 

  611.  ENDDO
    612. 

  612.  !$OMP END PARALLEL DO
    613. 

  613. 

  614. !Filter
    615. 

  615. 

  616. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    617. 

  617. # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
    618. 

  618. 

  619.      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
    620. 

  620.                                 ids, ide, jds, jde, kds, kde,         &
    621. 

  621.                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
    622. 

  622.                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
    623. 

  623. 

  624. # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
    625. 

  625. 

  626.      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
    627. 

  627.                                 ids, ide, jds, jde, kds, kde,         &
    628. 

  628.                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
    629. 

  629.                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
    630. 

  630. 

  631. # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
    632. 

  632. 

  633. 

  634. #else
    635. 

  635.      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
    636. 

  636.                                 ids, ide, jds, jde, kds, kde,       &
    637. 

  637.                                 ims, ime, jms, jme, kms, kme,       &
    638. 

  638.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    639. 

  639. 

  640.      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
    641. 

  641.                                 ids, ide, jds, jde, kds, kde,       &
    642. 

  642.                                 ims, ime, jms, jme, kms, kme,       &
    643. 

  643.                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
    644. 

  644. #endif
    645. 

  645. 

  646. ! Calculate tendency
    647. 

  647. 

  648.  !$OMP PARALLEL DO   &
    649. 

  649.  !$OMP PRIVATE ( ij, i,j,k )
    650. 

  650.  DO ij = 1 , num_tiles
    651. 

  651. 

  652.    DO j=jts,jte
    653. 

  653.    DO k=kts,ktf
    654. 

  654.    DO i=its,ite
    655. 

  655.      RTHNDGDTEN(i,k,j) = gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
    656. 

  656.                         grid%dif_analysis(i,k,j)
    657. 

  657.    ENDDO
    658. 

  658.    ENDDO
    659. 

  659.    ENDDO
    660. 

  660. 

  661.  ENDDO
    662. 

  662.  !$OMP END PARALLEL DO
    663. 

  663. 

  664. ENDIF
    665. 

  665. 

  666. IF(gph > 0. ) then
    667. 

  667. 

  668.  !$OMP PARALLEL DO   &
    669. 

  669.  !$OMP PRIVATE ( ij, i,j,k )
    670. 

  670.  DO ij = 1 , num_tiles
    671. 

  671. 

  672.    DO j=jts,jte
    673. 

  673.    DO k=kts,kte
    674. 

  674.    DO i=its,ite
    675. 

  675.      val_analysis = ph_ndg_old(i,k,j) *( 1.0 - coef ) + ph_ndg_new(i,k,j) * coef
    676. 

  676. 

  677.      grid%dif_analysis(i,k,j) = val_analysis - ph3d(i,k,j)
    678. 

  678. 

  679.    ENDDO
    680. 

  680.    ENDDO
    681. 

  681.    ENDDO
    682. 

  682. 

  683.  ENDDO
    684. 

  684.  !$OMP END PARALLEL DO
    685. 

  685. 

  686. !Filter
    687. 

  687. 

  688. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
    689. 

  689. # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
    690. 

  690. 

  691.      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
    692. 

  692.                                 ids, ide, jds, jde, kds, kde,         &
    693. 

  693.                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
    694. 

  694.                                 ipsx, ipex, jpsx, jpex, kpsx, kpex  )
    695. 

  695. 

  696. # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
    697. 

  697. 

  698.      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
    699. 

  699.                                 ids, ide, jds, jde, kds, kde,         &
    700. 

  700.                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
    701. 

  701.                                 ipsy, ipey, jpsy, jpey, kpsy, kpey  )
    702. 

  702. 

  703. # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
    704. 

  704. 

  705. 

  706. #else
    707. 

  707.      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
    708. 

  708.                                 ids, ide, jds, jde, kds, kde,       &
    709. 

  709.                                 ims, ime, jms, jme, kms, kme,       &
    710. 

  710.                                 ips, ipe, jps, jpe, kps, kpe )
    711. 

  711. 

  712.      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
    713. 

  713.                                 ids, ide, jds, jde, kds, kde,       &
    714. 

  714.                                 ims, ime, jms, jme, kms, kme,       &
    715. 

  715.                                 ips, ipe, jps, jpe, kps, kpe  )
    716. 

  716. #endif
    717. 

  717. 

  718. ! Calculate tendency
    719. 

  719. 

  720.  !$OMP PARALLEL DO   &
    721. 

  721.  !$OMP PRIVATE ( ij, i,j,k )
    722. 

  722.  DO ij = 1 , num_tiles
    723. 

  723. 

  724.    DO j=jts,jte
    725. 

  725.    DO k=kts,kte
    726. 

  726.    DO i=its,ite
    727. 

  727.      RPHNDGDTEN(i,k,j) = gph * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
    728. 

  728.                         grid%dif_analysis(i,k,j)
    729. 

  729.    ENDDO
    730. 

  730.    ENDDO
    731. 

  731.    ENDDO
    732. 

  732. 

  733.  ENDDO
    734. 

  734.  !$OMP END PARALLEL DO
    735. 

  735. 

  736. ENDIF
    737. 

  737. 

  738. #endif
    739. 

  739. 

  740.    END SUBROUTINE spectral_nudging
    741. 

  741. 

  742. !------------------------------------------------------------------------------
    743. 

  743. 

  744. SUBROUTINE spectral_nudging_filter_3dx( f, nwave,            &
    745. 

  745.                             ids, ide, jds, jde, kds, kde,    &
    746. 

  746.                             ims, ime, jms, jme, kms, kme,    &
    747. 

  747.                             its, ite, jts, jte, kts, kte )
    748. 

  748. 

  749.   IMPLICIT NONE
    750. 

  750. 

  751.   INTEGER ,       INTENT(IN   ) :: nwave
    752. 

  752.   INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    753. 

  753.                                    ims, ime, jms, jme, kms, kme, &
    754. 

  754.                                    its, ite, jts, jte, kts, kte
    755. 

  755. 

  756.   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) ::  f
    757. 

  757. 

  758.   REAL , DIMENSION(1:ide-ids+1,1:kte-kts+1) :: sheet
    759. 

  759.   INTEGER ::  i, j, j_end, k, nx, ny
    760. 

  760. 

  761.   ! Variables will stay in domain form since this routine is meaningless
    762. 

  762.   ! unless tile extent is the same as domain extent in E/W direction, i.e.,
    763. 

  763.   ! the processor has access to all grid points in E/W direction.
    764. 

  764.   ! There may be other ways of doing FFTs, but we haven't learned them yet...
    765. 

    766. 

  766.   ! Check to make sure we have full access to all E/W points
    767. 

  767.   IF ((its /= ids) .OR. (ite /= ide)) THEN
    768. 

  768.      WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (its /= ids) or (ite /= ide)',its,ids,ite,ide
    769. 

  769.      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
    770. 

  770.   END IF
    771. 

  771. 

  772. 

  773.   nx = ide-ids+1 
    774. 

  774.   ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
    775. 

  775.   j_end = MIN(jte, jde-1)
    776. 

  776.   IF (j_end == jde-1) j_end = jde
    777. 

  777.   DO j = jts, j_end
    778. 

  778. 

  779.         DO k=kts,kte
    780. 

  780.         DO i=ids,ide-1
    781. 

  781.            sheet(i-ids+1,k-kts+1) = f(i,k,j)
    782. 

  782.         END DO
    783. 

  783.            sheet(ide,k-kts+1) = 0.
    784. 

  784.         END DO
    785. 

  785. 

  786.         CALL spectralnudgingfilterfft2dncar(nx,ny,nwave,sheet)
    787. 

  787. 

  788.         DO k=kts,kte
    789. 

  789.            DO i=ids,ide
    790. 

  790.               f(i,k,j) = sheet(i-ids+1,k-kts+1)
    791. 

  791.            END DO
    792. 

  792.         END DO
    793. 

  793.   END DO ! outer j (latitude) loop
    794. 

  794. 

  795. END SUBROUTINE spectral_nudging_filter_3dx
    796. 

  796. !------------------------------------------------------------------------------
    797. 

  797. 

  798. SUBROUTINE spectral_nudging_filter_3dy( f, nwave,   &
    799. 

  799.                             ids, ide, jds, jde, kds, kde,    &
    800. 

  800.                             ims, ime, jms, jme, kms, kme,    &
    801. 

  801.                             its, ite, jts, jte, kts, kte )
    802. 

  802. 

  803.   IMPLICIT NONE
    804. 

  804. 

  805.   INTEGER ,       INTENT(IN   ) :: nwave
    806. 

  806.   INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
    807. 

  807.                                    ims, ime, jms, jme, kms, kme, &
    808. 

  808.                                    its, ite, jts, jte, kts, kte
    809. 

  809. 

  810.   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) ::  f
    811. 

  811. 

  812.   REAL , DIMENSION(1:jde-jds+1,1:kte-kts+1) :: sheet
    813. 

  813.   INTEGER ::  i, j, i_end, k, nx, ny
    814. 

  814. 

  815.   ! Variables will stay in domain form since this routine is meaningless
    816. 

  816.   ! unless tile extent is the same as domain extent in S/N direction, i.e.,
    817. 

  817.   ! the processor has access to all grid points in S/N direction.
    818. 

  818.   ! There may be other ways of doing FFTs, but we haven't learned them yet...
    819. 

    820. 

  820.   ! Check to make sure we have full access to all S/N points
    821. 

  821.   IF ((jts /= jds) .OR. (jte /= jde)) THEN
    822. 

  822.      WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (jts /= jds) or (jte /= jde)',jts,ids,ite,ide
    823. 

  823.      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
    824. 

  824.   END IF
    825. 

  825. 

  826. 

  827.   nx = jde-jds+1
    828. 

  828.   ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
    829. 

  829.   i_end = MIN(ite, ide-1)
    830. 

  830.   IF (i_end == ide-1) i_end = ide
    831. 

  831.   DO i = its, i_end
    832. 

  832. 

  833.         DO k=kts,kte
    834. 

  834.         DO j=jds,jde
    835. 

  835.            sheet(j-jds+1,k-kts+1) = f(i,k,j)
    836. 

  836.         END DO
    837. 

  837.            sheet(jde,k-kts+1) = 0.
    838. 

  838.         END DO
    839. 

  839. 

  840.         CALL spectralnudgingfilterfft2dncar(nx,ny,nwave,sheet)
    841. 

  841. 

  842.         DO k=kts,kte
    843. 

  843.            DO j=jds,jde
    844. 

  844.               f(i,k,j) = sheet(j-jds+1,k-kts+1)
    845. 

  845.            END DO
    846. 

  846.         END DO
    847. 

  847.   END DO ! outer i (longitude) loop
    848. 

  848. 

  849. END SUBROUTINE spectral_nudging_filter_3dy
    850. 

  850. 

  851. !------------------------------------------------------------------------------
    852. 

  852. 

  853. SUBROUTINE spectralnudgingfilterfft2dncar(nx,ny,nwave,fin)
    854. 

  854.   IMPLICIT NONE
    855. 

  855.   INTEGER , INTENT(IN) :: nx, ny, nwave
    856. 

  856.   REAL , DIMENSION(nx,ny), INTENT(INOUT) :: fin
    857. 

  857. 

  858.   INTEGER :: i, j
    859. 

  859.   REAL, dimension(nx,ny) :: fp
    860. 

  860. 

  861.   INTEGER :: lensave, ier, nh, n1
    862. 

  862.   INTEGER :: lot, jump, n, inc, lenr, lensav, lenwrk
    863. 

  863.   REAL, DIMENSION(nx+15) :: wsave
    864. 

  864.   REAL, DIMENSION(nx,ny) :: work
    865. 

  865. 

  866. 

  867. 

  868. !  we are following the naming convention of the fftpack5 routines
    869. 

  869. 

  870.   n = nx
    871. 

  871.   lot = ny
    872. 

  872.   lensav = n+15
    873. 

  873.   inc = 1
    874. 

  874.   lenr = nx*ny
    875. 

  875.   jump = nx
    876. 

  876.   lenwrk = lenr
    877. 

  877. 

  878. !  forward transform
    879. 

  879. !  initialize coefficients, place in wsave
    880. 

  880. !   (should place this in init and save wsave at program start)
    881. 

  881. 

  882.   call rfftmi(n,wsave,lensav,ier)
    883. 

  883.   IF(ier /= 0) THEN
    884. 

  884.     write(0,*) ' error in rfftmi ',ier
    885. 

  885.   END IF
    886. 

  886. 

  887. !  do the forward transform
    888. 

  888. 

  889.   call rfftmf( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
    890. 

  890.   IF(ier /= 0) THEN
    891. 

  891.     write(0,*) ' error in rfftmf ',ier
    892. 

  892.   END IF
    893. 

  893. 

  894.   if(MOD(n,2) == 0) then
    895. 

  895.     nh = n/2 - 1
    896. 

  896.   else
    897. 

  897.     nh = (n-1)/2
    898. 

  898.   end if
    899. 

  899. 

  900. 

  901. ! filter all waves with wavenumber larger than nwave
    902. 

  902. 

  903.   fp = 1.
    904. 

  904. 

  905.   DO j=1,ny
    906. 

  906.      DO i=nwave+1,nh
    907. 

  907.          fp(2*i-1,j) = 0.
    908. 

  908.          fp(2*i,j) = 0.
    909. 

  909.      ENDDO
    910. 

  910.   ENDDO
    911. 

  911. 

  912.   DO j=1,ny
    913. 

  913.     DO i=1,nx
    914. 

  914.       fin(i,j) = fp(i,j)*fin(i,j)
    915. 

  915.     ENDDO
    916. 

  916.   ENDDO
    917. 

  917. 

  918. !  do the backward transform
    919. 

  919. 

  920.   call rfftmb( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
    921. 

  921.   IF(ier /= 0) THEN
    922. 

  922.     write(0,*) ' error in rfftmb ',ier
    923. 

  923.   END IF
    924. 

  924. 

  925. END SUBROUTINE spectralnudgingfilterfft2dncar
    926. 

  926. 

  927. !-------------------------------------------------------------------
    928. 

  928. 

  929.    SUBROUTINE fddaspnudginginit(id,rundgdten,rvndgdten,rthndgdten,rphndgdten, &
    930. 

  930.                run_hours,  &
    931. 

  931.                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph, &
    932. 

  932.                if_zfac_uv, k_zfac_uv, dk_zfac_uv, &
    933. 

  933.                if_zfac_t, k_zfac_t, dk_zfac_t, &
    934. 

  934.                if_zfac_ph, k_zfac_ph, dk_zfac_ph, &               
    935. 

  935.                guv, gt, gph, if_ramping, dtramp_min, end_fdda_hour, &
    936. 

  936.                xwavenum,ywavenum,                          &
    937. 

  937.                       restart, allowed_to_read,                    &
    938. 

  938.                       ids, ide, jds, jde, kds, kde,                &
    939. 

  939.                       ims, ime, jms, jme, kms, kme,                &
    940. 

  940.                       its, ite, jts, jte, kts, kte                 )
    941. 

  941. !-------------------------------------------------------------------
    942. 

  942.    IMPLICIT NONE
    943. 

  943. !-------------------------------------------------------------------
    944. 

  944. !
    945. 

  945.    INTEGER , INTENT(IN)         ::  id
    946. 

  946.    LOGICAL, INTENT(IN)          ::  restart, allowed_to_read
    947. 

  947.    INTEGER, INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
    948. 

  948.                                     ims, ime, jms, jme, kms, kme, &
    949. 

  949.                                     its, ite, jts, jte, kts, kte
    950. 

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

  951.                                                        rundgdten, &
    952. 

  952.                                                        rvndgdten, &
    953. 

  953.                                                       rthndgdten, &
    954. 

  954.                                                       rphndgdten
    955. 

  955.    INTEGER,  INTENT(IN)   ::      run_hours
    956. 

  956.    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
    957. 

  957.                                   if_no_pbl_nudging_ph, end_fdda_hour
    958. 

  958.    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_ph
    959. 

  959.    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t,  k_zfac_ph
    960. 

  960.    INTEGER,  INTENT(IN)   ::      dk_zfac_uv,  dk_zfac_t,  dk_zfac_ph
    961. 

  961.    INTEGER,  INTENT(IN)   ::      if_ramping
    962. 

  962.    INTEGER,  INTENT(IN)   ::      xwavenum,ywavenum
    963. 

  963.    REAL,     INTENT(IN)   ::      dtramp_min
    964. 

  964.    REAL, INTENT(IN)       ::      guv, gt, gph
    965. 

  965.    REAL                   ::      actual_end_fdda_min
    966. 

  966. 

  967.    INTEGER :: i, j, k
    968. 

  968. 

  969.    LOGICAL , EXTERNAL     ::      wrf_dm_on_monitor
    970. 

  970. 

  971.    CHARACTER (LEN=256) :: message
    972. 

  972. 

  973.    IF ( wrf_dm_on_monitor() ) THEN
    974. 

  974. 

  975.      IF( guv > 0.0) THEN
    976. 

  976.        WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
    977. 

  977.            'D0',id,' Spectral nudging for wind is turned on and Guv= ', guv,' xwave= ',xwavenum,' ywavenum= ',ywavenum
    978. 

  978.        CALL wrf_message(TRIM(message))
    979. 

  979.      ELSE IF( guv < 0.0 ) THEN
    980. 

  980.        CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
    981. 

  981.      ELSE
    982. 

  982.        WRITE(message,'(a,i1,a,e12.4)') &
    983. 

  983.            'D0',id,' Spectral nudging for wind is turned off and Guv= ', guv
    984. 

  984.        CALL wrf_message(TRIM(message))
    985. 

  985.      ENDIF
    986. 

  986. 

  987.      IF( gt > 0.0) THEN
    988. 

  988.        WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
    989. 

  989.            'D0',id,' Spectral nudging for temperature is turned on and Gt= ', gt,' xwave= ',xwavenum,' ywavenum= ',ywavenum
    990. 

  990.        CALL wrf_message(TRIM(message))
    991. 

  991.      ELSE IF( gt < 0.0 ) THEN
    992. 

  992.        CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
    993. 

  993.      ELSE
    994. 

  994.        WRITE(message,'(a,i1,a,e12.4)') &
    995. 

  995.            'D0',id,' Spectral nudging for temperature is turned off and Gt= ', gt
    996. 

  996.        CALL wrf_message(TRIM(message))
    997. 

  997.      ENDIF
    998. 

  998. 

  999.      IF( gph > 0.0) THEN
    1000. 

  1000.        WRITE(message,'(a,i1,a,e12.4,a,i4,a,i4)') &
    1001. 

  1001.          'D0',id,' Spectral nudging for geopotential is turned on and Gph= ', gph,' xwave= ',xwavenum,' ywavenum= ',ywavenum
    1002. 

  1002.        CALL wrf_message(TRIM(message))
    1003. 

  1003.      ELSE IF( gph < 0.0 ) THEN
    1004. 

  1004.        CALL wrf_error_fatal('In grid FDDA, Gph must be positive.')
    1005. 

  1005.      ELSE
    1006. 

  1006.        WRITE(message,'(a,i1,a,e12.4)') &
    1007. 

  1007.          'D0',id,' Spectral nudging for geopotential is turned off and Gph= ', gph
    1008. 

  1008.        CALL wrf_message(TRIM(message))
    1009. 

  1009.      ENDIF
    1010. 

  1010. 

  1011.      IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
    1012. 

  1012.         WRITE(message,'(a,i1,a)') &
    1013. 

  1013.            'D0',id,' Spectral nudging for wind is turned off within the PBL.'
    1014. 

  1014.         CALL wrf_message(TRIM(message))
    1015. 

  1015.              IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must be greater or equal than 1.')
    1016. 

  1016.      ELSEIF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
    1017. 

  1017.         WRITE(message,'(a,i1,a,i3)') &
    1018. 

  1018.            'D0',id,' Spectral nudging for wind is turned off below layer', k_zfac_uv
    1019. 

  1019.         CALL wrf_message(TRIM(message))
    1020. 

  1020.              IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must  be greater or equal than 1.')       
    1021. 

  1021.      ENDIF
    1022. 

  1022. 

  1023. 

  1024.      IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
    1025. 

  1025.         WRITE(message,'(a,i1,a)') &
    1026. 

  1026.            'D0',id,' Spectral nudging for temperature is turned off within the PBL.'
    1027. 

  1027.         CALL wrf_message(TRIM(message))
    1028. 

  1028.              IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
    1029. 

  1029.      ELSEIF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
    1030. 

  1030.         WRITE(message,'(a,i1,a,i3)') &
    1031. 

  1031.            'D0',id,' Spectral nudging for temperature is turned off below layer', k_zfac_t
    1032. 

  1032.         CALL wrf_message(TRIM(message))
    1033. 

  1033.             IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
    1034. 

  1034.      ENDIF
    1035. 

  1035. 

  1036. 

  1037.      IF( gph > 0.0 .AND. if_no_pbl_nudging_ph == 1 ) THEN
    1038. 

  1038.         WRITE(message,'(a,i1,a)') &
    1039. 

  1039.          'D0',id,' Spectral nudging for geopotential is turned off within the PBL.'
    1040. 

  1040.         CALL wrf_message(TRIM(message))
    1041. 

  1041.             IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
    1042. 

  1042.      ELSEIF( gph > 0.0 .AND. if_zfac_ph == 1 ) THEN
    1043. 

  1043.         WRITE(message,'(a,i1,a,i3)') &
    1044. 

  1044.           'D0',id,' Spectral nudging for geopotential is turned off below layer', &
    1045. 

  1045.            k_zfac_ph
    1046. 

  1046.         CALL wrf_message(TRIM(message))
    1047. 

  1047.             IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
    1048. 

  1048.      ENDIF
    1049. 

  1049. 

  1050.      IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
    1051. 

  1051.        IF( dtramp_min <= 0.0 ) THEN
    1052. 

  1052.          actual_end_fdda_min = end_fdda_hour*60.0
    1053. 

  1053.        ELSE
    1054. 

  1054.          actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
    1055. 

  1055.        ENDIF
    1056. 

  1056. 

  1057.        IF( actual_end_fdda_min <= run_hours*60. ) THEN
    1058. 

  1058.           WRITE(message,'(a,i1,a)') &
    1059. 

  1059.             'D0',id,' Spectral nudging is ramped down near the end of the nudging period,'
    1060. 

  1060.           CALL wrf_message(TRIM(message))
    1061. 

  1061. 

  1062.           WRITE(message,'(a,f6.2,a,f6.2,a)') &
    1063. 

  1063.              '      starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0,&
    1064. 

  1064.              'h, ending at ', actual_end_fdda_min/60.0,'h.'
    1065. 

  1065.           CALL wrf_message(TRIM(message))
    1066. 

  1066.        ENDIF
    1067. 

  1067.      ENDIF
    1068. 

  1068. 

  1069.    ENDIF
    1070. 

  1070. 

  1071.    IF(.not.restart) THEN
    1072. 

  1072.      DO j = jts,jte
    1073. 

  1073.      DO k = kts,kte
    1074. 

  1074.      DO i = its,ite
    1075. 

  1075.         rundgdten(i,k,j) = 0.
    1076. 

  1076.         rvndgdten(i,k,j) = 0.
    1077. 

  1077.         rthndgdten(i,k,j) = 0.
    1078. 

  1078.         rphndgdten(i,k,j) = 0.
    1079. 

  1079.      ENDDO
    1080. 

  1080.      ENDDO
    1081. 

  1081.      ENDDO
    1082. 

  1082.    ENDIF
    1083. 

  1083. 

  1084.    END SUBROUTINE fddaspnudginginit
    1085. 

  1085. !-------------------------------------------------------------------
    1086. 

  1086. 

  1087. END MODULE module_fdda_spnudging
    1088. 

  1088.