phys/module_fdda_spnudging.F

   1 !wrf:model_layer:physics
   2 !
   3 ! Code contributed by Gonzalo Miguez-Macho, Univ. de Santiago de Compostela, Galicia, Spain
   4 !
   5 ! Added capability to spectrally nudge water vapor mixing ratio, and added
   6 ! user-definable lid for nudging potential temperature and water vapor mixing
   7 ! ratio. (Tanya Spero, U.S. Environmental Protection Agency -- October 2017)
   8
   9 MODULE module_fdda_spnudging
  10
  11 #ifdef DM_PARALLEL
  12   USE module_dm , ONLY : ntasks_x, ntasks_y, local_communicator_x, local_communicator_y, data_order_xzy
  13 #endif
  14   USE module_wrf_error , ONLY : wrf_err_message
  15
  16 CONTAINS
  17 !
  18 !-------------------------------------------------------------------
  19 !
  20    SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fdda_hour, &
  21                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph,if_no_pbl_nudging_q,&
  22                if_zfac_uv, k_zfac_uv, dk_zfac_uv, &
  23                if_zfac_t, k_zfac_t, dk_zfac_t, &
  24                if_zfac_ph, k_zfac_ph, dk_zfac_ph, &
  25                if_zfac_q, k_zfac_q, dk_zfac_q, ktrop, &
  26                guv, gt, gph, gq, if_ramping, dtramp_min, xwavenum, ywavenum, &
  27                u3d,v3d,th3d,ph3d,qv3d,                 &
  28                u_ndg_old,v_ndg_old,t_ndg_old,ph_ndg_old,q_ndg_old,       &
  29                u_ndg_new,v_ndg_new,t_ndg_new,ph_ndg_new,q_ndg_new,       &
  30                RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RPHNDGDTEN,RQVNDGDTEN,&
  31                pblh, ht, z, z_at_w,                             &
  32                ids,ide, jds,jde, kds,kde,                           &
  33                ims,ime, jms,jme, kms,kme,                           &
  34                i_start,i_end, j_start,j_end, kts,kte, num_tiles, &
  35                ips,ipe,jps,jpe,kps,kpe,                   &
  36                imsx,imex,jmsx,jmex,kmsx,kmex,                    &
  37                ipsx,ipex,jpsx,jpex,kpsx,kpex,                   &
  38                imsy,imey,jmsy,jmey,kmsy,kmey,                    &
  39                ipsy,ipey,jpsy,jpey,kpsy,kpey                   )
  40
  41
  42
  43    USE module_state_description
  44    USE module_domain, ONLY : domain
  45
  46 !-------------------------------------------------------------------
  47    implicit none
  48 !-------------------------------------------------------------------
  49 !-- u3d         3d u-velocity staggered on u points
  50 !-- v3d         3d v-velocity staggered on v points
  51 !-- th3d        3d potential temperature (k)
  52 !-- ph3d        3d perturbation geopotential
  53 !-- qv3d        3d water vapor mixing ratio
  54 !-- rundgdten   staggered u tendency due to
  55 !               spectral nudging (m/s/s)
  56 !-- rvndgdten   staggered v tendency due to
  57 !               spectral nudging (m/s/s)
  58 !-- rthndgdten  theta tendency due to
  59 !               spectral nudging (K/s)
  60 !-- phndgdten  ph tendency due to
  61 !               spectral nudging (m2/s2/s)
  62 !-- qvndgdten   qv tendency due to
  63 !               spectral nudging (kg/kg/s)
  64 !-- ids         start index for i in domain
  65 !-- ide         end index for i in domain
  66 !-- jds         start index for j in domain
  67 !-- jde         end index for j in domain
  68 !-- kds         start index for k in domain
  69 !-- kde         end index for k in domain
  70 !-- ims         start index for i in memory
  71 !-- ime         end index for i in memory
  72 !-- jms         start index for j in memory
  73 !-- jme         end index for j in memory
  74 !-- kms         start index for k in memory
  75 !-- kme         end index for k in memory
  76 !-- its         start index for i in tile
  77 !-- ite         end index for i in tile
  78 !-- jts         start index for j in tile
  79 !-- jte         end index for j in tile
  80 !-- kts         start index for k in tile
  81 !-- kte         end index for k in tile
  82 !-------------------------------------------------------------------
  83 !
  84    TYPE(domain) , TARGET          :: grid
  85
  86    INTEGER,  INTENT(IN)   ::      itimestep, analysis_interval, end_fdda_hour
  87
  88    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
  89                                   if_no_pbl_nudging_ph, if_no_pbl_nudging_q
  90    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_ph, if_zfac_q
  91    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t, k_zfac_ph, k_zfac_q
  92    INTEGER,  INTENT(IN)   ::      dk_zfac_uv,  dk_zfac_t, dk_zfac_ph, dk_zfac_q
  93    INTEGER,  INTENT(IN)   ::      ktrop
  94    INTEGER,  INTENT(IN)   ::      if_ramping
  95    INTEGER,  INTENT(IN)   ::      xwavenum,ywavenum
  96
  97    INTEGER , INTENT(IN)   ::      id
  98    REAL,     INTENT(IN)   ::      DT, xtime, dtramp_min
  99
 100    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
 101                                   ims,ime, jms,jme, kms,kme, &
 102                                   kts,kte, num_tiles,        &
 103                                   ips,ipe,jps,jpe,kps,kpe,   &
 104                                   imsx,imex,jmsx,jmex,kmsx,kmex,   &
 105                                   ipsx,ipex,jpsx,jpex,kpsx,kpex,   &
 106                                   imsy,imey,jmsy,jmey,kmsy,kmey,   &
 107                                   ipsy,ipey,jpsy,jpey,kpsy,kpey
 108
 109    INTEGER, DIMENSION(num_tiles), INTENT(IN) ::                   &
 110   &                                    i_start,i_end,j_start,j_end
 111
 112    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
 113              INTENT(IN)   ::                  ph3d, &
 114                                               qv3d, &
 115                                               th3d, &
 116                                                  z, &
 117                                             z_at_w
 118
 119    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
 120              INTENT(INOUT)   ::           rundgdten, &
 121                                           rvndgdten, &
 122                                          rthndgdten, &
 123                                          rphndgdten, &
 124                                          rqvndgdten
 125
 126
 127    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
 128              INTENT(INOUT)   ::           u_ndg_old, &
 129                                           v_ndg_old, &
 130                                           t_ndg_old, &
 131                                           ph_ndg_old, &
 132                                           q_ndg_old, &
 133                                           u_ndg_new, &
 134                                           v_ndg_new, &
 135                                           t_ndg_new, &
 136                                           ph_ndg_new, &
 137                                           q_ndg_new
 138
 139    REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ), &
 140              INTENT(IN)   ::                    u3d, &
 141                                                 v3d
 142
 143    REAL,  DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
 144                                                          ht
 145
 146    REAL, INTENT(IN)    :: guv, gt ,gph, gq
 147
 148    INTEGER             :: its,ite, jts,jte,ij
 149    INTEGER             :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
 150    REAL                :: xtime_old, xtime_new, coef, val_analysis
 151    INTEGER             :: kpbl, dbg_level
 152
 153    REAL                :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
 154    REAL, DIMENSION( ips:ipe, kps:kpe, jps:jpe, 5 ) :: wpbl  ! 1: u, 2: v, 3: t, 4: ph, 5: q
 155    REAL, DIMENSION( kps:kpe, 5 )                   :: wzfac ! 1: u, 2: v, 3: t, 4: ph, 5: q
 156
 157    LOGICAL , EXTERNAL  :: wrf_dm_on_monitor
 158
 159    CHARACTER (LEN=256) :: wrf_err_message
 160
 161 #if  ! ( NMM_CORE == 1 )
 162
 163       DO ij = 1 , num_tiles
 164           its=i_start(ij)
 165           ite=i_end(ij)
 166           jts=j_start(ij)
 167           jte=j_end(ij)
 168       ENDDO
 169 !GMM default values for vertical coefficients, set here
 170    wpbl(:,:,:,:) = 1.0
 171    wzfac(:,:) = 1.0
 172
 173  !$OMP PARALLEL DO   &
 174  !$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation.
 175  DO ij = 1 , num_tiles
 176      DO j = jts, jte
 177      DO k = kts, kte
 178      DO i = its, ite
 179        RUNDGDTEN(i,k,j) = 0.0
 180        RVNDGDTEN(i,k,j) = 0.0
 181        RTHNDGDTEN(i,k,j) = 0.0
 182        RPHNDGDTEN(i,k,j) = 0.0
 183        RQVNDGDTEN(i,k,j) = 0.0
 184      ENDDO
 185      ENDDO
 186      ENDDO
 187  ENDDO
 188  !$OMP END PARALLEL DO
 189
 190  actual_end_fdda_min = end_fdda_hour*60.0
 191  IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
 192        actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
 193  IF( xtime > actual_end_fdda_min ) RETURN
 194
 195  !$OMP PARALLEL DO   &
 196  !$OMP PRIVATE ( ij, i,j,k ) ! only for part of the calculation.
 197  !GMM Transpose and filtering needs to be done on patch
 198
 199  DO ij = 1 , num_tiles
 200
 201 !  actual_end_fdda_min = end_fdda_hour*60.0
 202 !  IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
 203 !      actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
 204
 205    xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
 206    xtime_new = xtime_old + analysis_interval
 207    coef = (xtime-xtime_old)/(xtime_new-xtime_old)
 208
 209    IF ( wrf_dm_on_monitor()) THEN
 210
 211      CALL get_wrf_debug_level( dbg_level )
 212
 213      IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
 214
 215        IF( xtime < end_fdda_hour*60.0 ) THEN
 216          WRITE(wrf_err_message,FMT='(a,i2.2,a,f15.3,a)') &
 217           'D',id,' Spectral nudging read in new data at time = ', xtime, ' min.'
 218          CALL wrf_message( TRIM(wrf_err_message) )
 219          WRITE(wrf_err_message,FMT='(a,i2.2,a,2(f15.3,1x),a)') &
 220           'D',id,' Spectral nudging bracketing times = ', xtime_old, xtime_new, ' min.'
 221          CALL wrf_message( TRIM(wrf_err_message) )
 222        ENDIF
 223
 224        actual_end_fdda_min = end_fdda_hour*60.0
 225        IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
 226            actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
 227
 228        IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
 229 !        Find the mid point of the tile and print out the sample values
 230          i0 = (ite-its)/2+its
 231          j0 = (jte-jts)/2+jts
 232
 233          IF( guv > 0.0 ) THEN
 234            DO k = kts, kte
 235              WRITE(wrf_err_message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
 236                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
 237                ' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
 238              CALL wrf_message( TRIM(wrf_err_message) )
 239            ENDDO
 240            DO k = kts, kte
 241              WRITE(wrf_err_message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
 242                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
 243                ' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
 244              CALL wrf_message( TRIM(wrf_err_message) )
 245            ENDDO
 246          ENDIF
 247
 248          IF( gt > 0.0 ) THEN
 249            DO k = kts, kte
 250              WRITE(wrf_err_message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
 251                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
 252                ' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
 253              CALL wrf_message( TRIM(wrf_err_message) )
 254            ENDDO
 255          ENDIF
 256
 257          IF( gph > 0.0 ) THEN
 258            DO k = kts, kte
 259              WRITE(wrf_err_message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
 260                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
 261                ' ph_ndg_old=', ph_ndg_old(i0,k,j0), ' ph_ndg_new=', ph_ndg_new(i0,k,j0)
 262              CALL wrf_message( TRIM(wrf_err_message) )
 263            ENDDO
 264          ENDIF
 265
 266          IF( gq > 0.0 ) THEN
 267            DO k = kts, kte
 268              WRITE(wrf_err_message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
 269                '    D0',id,' sample analysis values at i,k,j=', i0, k, j0, &
 270                ' q_ndg_old=', q_ndg_old(i0,k,j0), ' q_ndg_new=', q_ndg_new(i0,k,j0)
 271              CALL wrf_message( TRIM(wrf_err_message) )
 272            ENDDO
 273          ENDIF
 274
 275        ENDIF
 276      ENDIF
 277    ENDIF
 278
 279    jtsv=MAX0(jts,jds+1)
 280    itsu=MAX0(its,ids+1)
 281
 282    jtf=MIN0(jte,jde-1)
 283    ktf=MIN0(kte,kde-1)
 284    itf=MIN0(ite,ide-1)
 285 !
 286 ! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t,
 287 ! if_no_pbl_nudging_ph swithes) are set to 1, compute the weighting function, wpbl(:,k,:,:),
 288 ! based on the PBL depth.  wpbl = 1 above the PBL and wpbl = 0 in the PBL.  If all
 289 ! the switche are set to zero, wpbl = 1 (default value).
 290 !
 291 !GMM If those are set to zero, then check if the user-defined namelist switches (if_zfac_uv, if_zfac_t,
 292 ! if_zfac_ph swithes) are set to 1, compute the weighting function, wzfac(k,:),
 293 ! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_ph)
 294 ! below which analysis nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x).
 295 !  The strength of nudging increases linearly from k_zfac to kzfac + dk_zfac
 296 ! (dk_zfac_uv, dk_zfac_t and kd_zfac_ph are also set in the namelist, default value is 1).
 297 !If all switches are set to zero, wzfac = 1 (default value).
 298 !
 299 ! TLS: The above conditions also apply to Qv. (October 2017)
 300 ! TLS: Allow user-definable lid (KTROP, via the namelist) that can be used to
 301 !      restrict nudging of potential temperature and water vapor mixing ratio
 302 !      above a layer.  Nudging of both theta and Qv will be at 0.5 strength
 303 !      at the user-defined layer and will be 0.0 above that layer.  Setting
 304 !      KTROP=0 will allow nudging to the top of the model.  U, V, andd PHI are
 305 !      not affected by KTROP.  KTROP is separate from the PBL and ZFAC options.
 306 !      Using KTROP with the PBL restriction can effectively limit nudging of
 307 !      theta and Qv to the free troposphere.  Ideally KTROP will be refined or
 308 !      supplemented with the space and time varying tropopause.  The lower limit
 309 !      for KTROP is arbitrarily set to the index of the middle layer of the
 310 !      number of layers in a given simulation.  (October 2017)
 311
 312    IF( if_no_pbl_nudging_uv == 1 ) THEN
 313
 314      DO j=jts,jtf
 315      DO i=itsu,itf
 316
 317        kpbl = 1
 318        zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
 319
 320        loop_ku: DO k=kts,ktf
 321          zagl_bot = 0.5 * ( z_at_w(i-1,k,  j)-ht(i-1,j) + z_at_w(i,k,  j)-ht(i,j) )
 322          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) )
 323          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
 324            kpbl = k
 325            EXIT loop_ku
 326          ENDIF
 327        ENDDO loop_ku
 328
 329        DO k=kts,ktf
 330           wpbl(i, k, j, 1) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
 331        ENDDO
 332
 333      ENDDO
 334      ENDDO
 335
 336      DO i=its,itf
 337      DO j=jtsv,jtf
 338
 339        kpbl = 1
 340        zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
 341
 342        loop_kv: DO k=kts,ktf
 343          zagl_bot = 0.5 * ( z_at_w(i,k,  j-1)-ht(i,j-1) + z_at_w(i,k,  j)-ht(i,j) )
 344          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) )
 345          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
 346            kpbl = k
 347            EXIT loop_kv
 348          ENDIF
 349        ENDDO loop_kv
 350
 351        DO k=kts,ktf
 352           wpbl(i, k, j, 2) = max ( min ( float(k-kpbl) / float(dk_zfac_uv) , 1. ) , 0.)
 353        ENDDO
 354
 355      ENDDO
 356      ENDDO
 357
 358    ELSEIF( if_zfac_uv == 1 ) THEN
 359
 360      DO j=jts,jte
 361      DO k=kts,ktf
 362      DO i=its,ite
 363           wpbl(i, k, j, 1) = max ( min ( float(k-k_zfac_uv) / float(dk_zfac_uv) , 1. ) , 0.)
 364           wpbl(i, k, j, 2) = wpbl(i, k, j, 1)
 365      ENDDO
 366      ENDDO
 367      ENDDO
 368
 369    ENDIF
 370
 371
 372    IF( if_no_pbl_nudging_t == 1 ) THEN
 373
 374      DO j=jts,jtf
 375      DO i=its,itf
 376
 377        kpbl = 1
 378        zpbl = pblh(i,j)
 379
 380        loop_kt: DO k=kts,ktf
 381          zagl_bot = z_at_w(i,k,  j)-ht(i,j)
 382          zagl_top = z_at_w(i,k+1,j)-ht(i,j)
 383          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
 384            kpbl = k
 385            EXIT loop_kt
 386          ENDIF
 387        ENDDO loop_kt
 388
 389        DO k=kts,ktf
 390           wpbl(i, k, j, 3) = max ( min ( float(k-kpbl) / float(dk_zfac_t) , 1. ) , 0.)
 391        ENDDO
 392
 393      ENDDO
 394      ENDDO
 395
 396    ELSEIF( if_zfac_t == 1 ) THEN
 397
 398      DO j=jts,jtf
 399      DO k=kts,ktf
 400      DO i=its,itf
 401           wpbl(i, k, j, 3) = max ( min ( float(k-k_zfac_t) / float(dk_zfac_t) , 1. ) , 0.)
 402      ENDDO
 403      ENDDO
 404      ENDDO
 405
 406    ENDIF
 407
 408
 409 !~~~Restricting nudging of temperature above a level KTROP to consider
 410 !~~~deficiencies and gradients in temperature above tropopause.
 411      IF ( ktrop > 0 ) THEN
 412
 413        DO j=jts,jtf
 414        DO i=its,itf
 415            wpbl(i, ktrop, j, 3) = wpbl(i, ktrop, j, 3) * 0.5
 416        ENDDO
 417        ENDDO
 418
 419        DO j=jts,jtf
 420        DO k=ktrop+1,ktf
 421        DO i=its,itf
 422             wpbl(i, k, j, 3) = 0.0
 423        ENDDO
 424        ENDDO
 425        ENDDO
 426
 427      ENDIF
 428 !~~~End KTROP mods for potential temperature
 429
 430
 431    IF( if_no_pbl_nudging_ph == 1 ) THEN
 432
 433      DO j=jts,jtf
 434      DO i=its,itf
 435
 436        kpbl = 1
 437        zpbl = pblh(i,j)
 438
 439        loop_kph: DO k=kts,kte
 440          zagl = z_at_w(i,k,  j)-ht(i,j)
 441          IF( zpbl >= zagl) THEN
 442            kpbl = k
 443            EXIT loop_kph
 444          ENDIF
 445        ENDDO loop_kph
 446
 447        DO k=kts,kte
 448           wpbl(i, k, j, 4) = max ( min ( float(k-kpbl) / float(dk_zfac_ph) , 1. ) , 0.)
 449        ENDDO
 450
 451      ENDDO
 452      ENDDO
 453
 454    ELSEIF( if_zfac_ph == 1 ) THEN
 455
 456      DO j=jts,jtf
 457      DO k=kts,kte
 458      DO i=its,itf
 459           wpbl(i, k, j, 4) = max ( min ( float(k-k_zfac_ph) / float(dk_zfac_ph) , 1. ) , 0.)
 460      ENDDO
 461      ENDDO
 462      ENDDO
 463
 464    ENDIF
 465
 466
 467    IF( if_no_pbl_nudging_q == 1 ) THEN
 468
 469      DO j=jts,jtf
 470      DO i=its,itf
 471
 472        kpbl = 1
 473        zpbl = pblh(i,j)
 474
 475        loop_kq: DO k=kts,ktf
 476          zagl_bot = z_at_w(i,k,  j)-ht(i,j)
 477          zagl_top = z_at_w(i,k+1,j)-ht(i,j)
 478          IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
 479            kpbl = k
 480            EXIT loop_kq
 481          ENDIF
 482        ENDDO loop_kq
 483
 484        DO k=kts,ktf
 485           wpbl(i, k, j, 5) = max ( min ( float(k-kpbl) / float(dk_zfac_q) , 1. ) , 0.)
 486        ENDDO
 487
 488      ENDDO
 489      ENDDO
 490
 491    ELSEIF( if_zfac_q == 1 ) THEN
 492
 493      DO j=jts,jtf
 494      DO k=kts,ktf
 495      DO i=its,itf
 496           wpbl(i, k, j, 5) = max ( min ( float(k-k_zfac_q) / float(dk_zfac_q) , 1. ) , 0.)
 497      ENDDO
 498      ENDDO
 499      ENDDO
 500
 501    ENDIF
 502
 503 !~~~Restricting nudging of moisture above a level KTROP to consider
 504 !~~~deficiencies in moisture above tropopause.
 505      IF ( ktrop > 0 ) THEN
 506
 507        DO j=jts,jtf
 508        DO i=its,itf
 509            wpbl(i, ktrop, j, 5) = wpbl(i, ktrop, j, 5) * 0.5
 510        ENDDO
 511        ENDDO
 512
 513        DO j=jts,jtf
 514        DO k=ktrop+1,ktf
 515        DO i=its,itf
 516             wpbl(i, k, j, 5) = 0.0
 517        ENDDO
 518        ENDDO
 519        ENDDO
 520
 521      ENDIF
 522 !~~~End KTROP mods for water vapor mixing ratio
 523
 524
 525
 526 ! If if_ramping and dtramp_min are defined by user, comput a time weighting function, tfac,
 527 ! for analysis nudging so that at the end of the nudging period (which has to be at a
 528 ! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
 529 !
 530 ! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at
 531 ! end_fdda_hour-ABS(dtramp_min).
 532 !
 533 ! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at
 534 ! end_fdda_hour+ABS(dtramp_min). In this case, the obs values are extrapolated using
 535 ! the obs tendency saved from the previous FDDA wondow.  More specifically for extrapolation,
 536 ! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
 537 !
 538    tfac = 1.0
 539
 540    IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
 541
 542      IF( dtramp_min <= 0.0 ) THEN
 543        actual_end_fdda_min = end_fdda_hour*60.0
 544      ELSE
 545        actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
 546      ENDIF
 547
 548      IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN
 549        tfac = 1.0
 550      ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
 551        tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
 552        IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/analysis_interval
 553      ELSE
 554        tfac = 0.0
 555      ENDIF
 556
 557    ENDIF
 558
 559    ENDDO
 560    !$OMP END PARALLEL DO
 561
 562 !
 563 ! Compute 3-D nudging tendencies for u, v
 564 !
 565 !
 566 !GMM Fist calculate differences between model variables and analysis values,
 567 !then filter in the x and y direction all wave numbers higher than
 568 ! xwavenum and ywavenum, as specified in the namelist.
 569 !If either xwavenum or ywavenum are not specified,
 570 ! default values are zero, and spectral nudging is turned off
 571 !Then use the filtered differences to calculate the spectral nudging tendencies
 572
 573 IF(guv > 0. ) then
 574
 575  !$OMP PARALLEL DO   &
 576  !$OMP PRIVATE ( ij, i,j,k )
 577  DO ij = 1 , num_tiles
 578
 579    DO j=jts,jte
 580    DO k=kts,ktf
 581    DO i=its,ite
 582      val_analysis = u_ndg_old(i,k,j) *( 1.0 - coef ) + u_ndg_new(i,k,j) * coef
 583
 584      grid%dif_analysis(i,k,j) = val_analysis - u3d(i,k,j)
 585
 586    ENDDO
 587    ENDDO
 588    ENDDO
 589
 590  ENDDO
 591  !$OMP END PARALLEL DO
 592
 593 !Filter
 594
 595 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 596 # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
 597
 598      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
 599                                 ids, ide, jds, jde, kds, kde,         &
 600                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 601                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
 602
 603 # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
 604
 605      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
 606                                 ids, ide, jds, jde, kds, kde,         &
 607                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
 608                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
 609
 610 # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
 611
 612
 613 #else
 614      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
 615                                 ids, ide, jds, jde, kds, kde,       &
 616                                 ims, ime, jms, jme, kms, kme,       &
 617                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 618
 619      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
 620                                 ids, ide, jds, jde, kds, kde,       &
 621                                 ims, ime, jms, jme, kms, kme,       &
 622                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 623 #endif
 624
 625 ! Calculate tendency
 626
 627  !$OMP PARALLEL DO   &
 628  !$OMP PRIVATE ( ij, i,j,k )
 629  DO ij = 1 , num_tiles
 630
 631    DO j=jts,jte
 632    DO k=kts,ktf
 633    DO i=its,ite
 634      RUNDGDTEN(i,k,j) = guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
 635                         grid%dif_analysis(i,k,j)
 636    ENDDO
 637    ENDDO
 638    ENDDO
 639
 640
 641 !
 642 ! Now V component
 643 !
 644
 645    DO j=jts,jte
 646    DO k=kts,ktf
 647    DO i=its,ite
 648      val_analysis = v_ndg_old(i,k,j) *( 1.0 - coef ) + v_ndg_new(i,k,j) * coef
 649
 650      grid%dif_analysis(i,k,j) = val_analysis - v3d(i,k,j)
 651
 652    ENDDO
 653    ENDDO
 654    ENDDO
 655
 656  ENDDO
 657  !$OMP END PARALLEL DO
 658
 659 !
 660 ! Filter
 661 !
 662 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 663 # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
 664      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,     &
 665                                 ids, ide, jds, jde, kds, kde,         &
 666                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 667                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
 668
 669 # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
 670      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,   &
 671                                 ids, ide, jds, jde, kds, kde-1,         &
 672                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
 673                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
 674
 675 # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
 676
 677
 678 #else
 679      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,    &
 680                                 ids, ide, jds, jde, kds, kde,       &
 681                                 ims, ime, jms, jme, kms, kme,       &
 682                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 683
 684      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,    &
 685                                 ids, ide, jds, jde, kds, kde,       &
 686                                 ims, ime, jms, jme, kms, kme,       &
 687                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 688 #endif
 689
 690 ! Calculate tendency
 691
 692  !$OMP PARALLEL DO   &
 693  !$OMP PRIVATE ( ij, i,j,k )
 694  DO ij = 1 , num_tiles
 695
 696    DO j=jts,jte
 697    DO k=kts,ktf
 698    DO i=its,ite
 699      RVNDGDTEN(i,k,j) = guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
 700                         grid%dif_analysis(i,k,j)
 701    ENDDO
 702    ENDDO
 703    ENDDO
 704
 705  ENDDO
 706  !$OMP END PARALLEL DO
 707
 708 ENDIF
 709
 710 IF(gt > 0. ) then
 711
 712  !$OMP PARALLEL DO   &
 713  !$OMP PRIVATE ( ij, i,j,k )
 714  DO ij = 1 , num_tiles
 715
 716    DO j=jts,jte
 717    DO k=kts,ktf
 718    DO i=its,ite
 719      val_analysis = t_ndg_old(i,k,j) *( 1.0 - coef ) + t_ndg_new(i,k,j) * coef
 720
 721      grid%dif_analysis(i,k,j) = val_analysis - th3d(i,k,j) + 300.
 722
 723    ENDDO
 724    ENDDO
 725    ENDDO
 726
 727  ENDDO
 728  !$OMP END PARALLEL DO
 729
 730 !Filter
 731
 732 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 733 # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
 734
 735      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
 736                                 ids, ide, jds, jde, kds, kde,         &
 737                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 738                                 ipsx, ipex, jpsx, jpex, kpsx, MIN(kde-1,kpex) )
 739
 740 # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
 741
 742      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
 743                                 ids, ide, jds, jde, kds, kde,         &
 744                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
 745                                 ipsy, ipey, jpsy, jpey, kpsy, MIN(kde-1,kpey) )
 746
 747 # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
 748
 749
 750 #else
 751      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
 752                                 ids, ide, jds, jde, kds, kde,       &
 753                                 ims, ime, jms, jme, kms, kme,       &
 754                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 755
 756      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
 757                                 ids, ide, jds, jde, kds, kde,       &
 758                                 ims, ime, jms, jme, kms, kme,       &
 759                                 ips, ipe, jps, jpe, kps, MIN(kde-1,kpe) )
 760 #endif
 761
 762 ! Calculate tendency
 763
 764  !$OMP PARALLEL DO   &
 765  !$OMP PRIVATE ( ij, i,j,k )
 766  DO ij = 1 , num_tiles
 767
 768    DO j=jts,jte
 769    DO k=kts,ktf
 770    DO i=its,ite
 771      RTHNDGDTEN(i,k,j) = gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
 772                         grid%dif_analysis(i,k,j)
 773    ENDDO
 774    ENDDO
 775    ENDDO
 776
 777  ENDDO
 778  !$OMP END PARALLEL DO
 779
 780 ENDIF
 781
 782 IF(gph > 0. ) then
 783
 784  !$OMP PARALLEL DO   &
 785  !$OMP PRIVATE ( ij, i,j,k )
 786  DO ij = 1 , num_tiles
 787
 788    DO j=jts,jte
 789    DO k=kts,kte
 790    DO i=its,ite
 791      val_analysis = ph_ndg_old(i,k,j) *( 1.0 - coef ) + ph_ndg_new(i,k,j) * coef
 792
 793      grid%dif_analysis(i,k,j) = val_analysis - ph3d(i,k,j)
 794
 795    ENDDO
 796    ENDDO
 797    ENDDO
 798
 799  ENDDO
 800  !$OMP END PARALLEL DO
 801
 802 !Filter
 803
 804 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 805 # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
 806
 807      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
 808                                 ids, ide, jds, jde, kds, kde,         &
 809                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 810                                 ipsx, ipex, jpsx, jpex, kpsx, kpex  )
 811
 812 # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
 813
 814      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
 815                                 ids, ide, jds, jde, kds, kde,         &
 816                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
 817                                 ipsy, ipey, jpsy, jpey, kpsy, kpey  )
 818
 819 # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
 820
 821
 822 #else
 823      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
 824                                 ids, ide, jds, jde, kds, kde,       &
 825                                 ims, ime, jms, jme, kms, kme,       &
 826                                 ips, ipe, jps, jpe, kps, kpe )
 827
 828      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
 829                                 ids, ide, jds, jde, kds, kde,       &
 830                                 ims, ime, jms, jme, kms, kme,       &
 831                                 ips, ipe, jps, jpe, kps, kpe  )
 832 #endif
 833
 834 ! Calculate tendency
 835
 836  !$OMP PARALLEL DO   &
 837  !$OMP PRIVATE ( ij, i,j,k )
 838  DO ij = 1 , num_tiles
 839
 840    DO j=jts,jte
 841    DO k=kts,kte
 842    DO i=its,ite
 843      RPHNDGDTEN(i,k,j) = gph * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
 844                         grid%dif_analysis(i,k,j)
 845    ENDDO
 846    ENDDO
 847    ENDDO
 848
 849  ENDDO
 850  !$OMP END PARALLEL DO
 851
 852 ENDIF
 853
 854
 855 IF(gq > 0. ) then
 856
 857  !$OMP PARALLEL DO   &
 858  !$OMP PRIVATE ( ij, i,j,k )
 859  DO ij = 1 , num_tiles
 860
 861    DO j=jts,jte
 862    DO k=kts,kte
 863    DO i=its,ite
 864      val_analysis = q_ndg_old(i,k,j) *( 1.0 - coef ) + q_ndg_new(i,k,j) * coef
 865
 866      grid%dif_analysis(i,k,j) = val_analysis - qv3d(i,k,j)
 867
 868    ENDDO
 869    ENDDO
 870    ENDDO
 871
 872  ENDDO
 873  !$OMP END PARALLEL DO
 874
 875 !Filter
 876
 877 #ifdef DM_PARALLEL
 878 # include "XPOSE_SPECTRAL_NUDGING_z2x.inc"
 879
 880      CALL spectral_nudging_filter_3dx( grid%dif_xxx, xwavenum,    &
 881                                 ids, ide, jds, jde, kds, kde,         &
 882                                 imsx, imex, jmsx, jmex, kmsx, kmex,     &
 883                                 ipsx, ipex, jpsx, jpex, kpsx, kpex  )
 884
 885 # include "XPOSE_SPECTRAL_NUDGING_x2y.inc"
 886
 887      CALL spectral_nudging_filter_3dy( grid%dif_yyy, ywavenum,    &
 888                                 ids, ide, jds, jde, kds, kde,         &
 889                                 imsy, imey, jmsy, jmey, kmsy, kmey,     &
 890                                 ipsy, ipey, jpsy, jpey, kpsy, kpey  )
 891
 892 # include "XPOSE_SPECTRAL_NUDGING_y2z.inc"
 893
 894
 895 #else
 896      CALL spectral_nudging_filter_3dx( grid%dif_analysis, xwavenum,     &
 897                                 ids, ide, jds, jde, kds, kde,       &
 898                                 ims, ime, jms, jme, kms, kme,       &
 899                                 ips, ipe, jps, jpe, kps, kpe )
 900
 901      CALL spectral_nudging_filter_3dy( grid%dif_analysis, ywavenum,      &
 902                                 ids, ide, jds, jde, kds, kde,       &
 903                                 ims, ime, jms, jme, kms, kme,       &
 904                                 ips, ipe, jps, jpe, kps, kpe  )
 905 #endif
 906
 907 ! Calculate tendency
 908
 909  !$OMP PARALLEL DO   &
 910  !$OMP PRIVATE ( ij, i,j,k )
 911  DO ij = 1 , num_tiles
 912
 913    DO j=jts,jte
 914    DO k=kts,kte
 915    DO i=its,ite
 916      RQVNDGDTEN(i,k,j) = gq * wpbl(i,k,j,5) * wzfac(k,5) * tfac * &
 917                         grid%dif_analysis(i,k,j)
 918    ENDDO
 919    ENDDO
 920    ENDDO
 921
 922  ENDDO
 923  !$OMP END PARALLEL DO
 924
 925 ENDIF
 926
 927
 928 #endif
 929
 930    END SUBROUTINE spectral_nudging
 931
 932 !------------------------------------------------------------------------------
 933
 934 SUBROUTINE spectral_nudging_filter_3dx( f, nwave,            &
 935                             ids, ide, jds, jde, kds, kde,    &
 936                             ims, ime, jms, jme, kms, kme,    &
 937                             its, ite, jts, jte, kts, kte )
 938
 939   IMPLICIT NONE
 940
 941   INTEGER ,       INTENT(IN   ) :: nwave
 942   INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 943                                    ims, ime, jms, jme, kms, kme, &
 944                                    its, ite, jts, jte, kts, kte
 945
 946   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) ::  f
 947
 948   REAL , DIMENSION(1:ide-ids+1,1:kte-kts+1) :: sheet
 949   INTEGER ::  i, j, j_end, k, nx, ny
 950
 951   ! Variables will stay in domain form since this routine is meaningless
 952   ! unless tile extent is the same as domain extent in E/W direction, i.e.,
 953   ! the processor has access to all grid points in E/W direction.
 954   ! There may be other ways of doing FFTs, but we haven't learned them yet...
 955
 956   ! Check to make sure we have full access to all E/W points
 957   IF ((its /= ids) .OR. (ite /= ide)) THEN
 958      WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (its /= ids) or (ite /= ide)',its,ids,ite,ide
 959      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
 960   END IF
 961
 962
 963   nx = ide-ids+1
 964   ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
 965   j_end = MIN(jte, jde-1)
 966   IF (j_end == jde-1) j_end = jde
 967   DO j = jts, j_end
 968
 969         DO k=kts,kte
 970         DO i=ids,ide-1
 971            sheet(i-ids+1,k-kts+1) = f(i,k,j)
 972         END DO
 973            sheet(ide,k-kts+1) = 0.
 974         END DO
 975
 976         CALL spectralnudgingfilterfft2dncar(nx,ny,nwave,sheet)
 977
 978         DO k=kts,kte
 979            DO i=ids,ide
 980               f(i,k,j) = sheet(i-ids+1,k-kts+1)
 981            END DO
 982         END DO
 983   END DO ! outer j (latitude) loop
 984
 985 END SUBROUTINE spectral_nudging_filter_3dx
 986 !------------------------------------------------------------------------------
 987
 988 SUBROUTINE spectral_nudging_filter_3dy( f, nwave,   &
 989                             ids, ide, jds, jde, kds, kde,    &
 990                             ims, ime, jms, jme, kms, kme,    &
 991                             its, ite, jts, jte, kts, kte )
 992
 993   IMPLICIT NONE
 994
 995   INTEGER ,       INTENT(IN   ) :: nwave
 996   INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 997                                    ims, ime, jms, jme, kms, kme, &
 998                                    its, ite, jts, jte, kts, kte
 999
1000   REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) ::  f
1001
1002   REAL , DIMENSION(1:jde-jds+1,1:kte-kts+1) :: sheet
1003   INTEGER ::  i, j, i_end, k, nx, ny
1004
1005   ! Variables will stay in domain form since this routine is meaningless
1006   ! unless tile extent is the same as domain extent in S/N direction, i.e.,
1007   ! the processor has access to all grid points in S/N direction.
1008   ! There may be other ways of doing FFTs, but we haven't learned them yet...
1009
1010   ! Check to make sure we have full access to all S/N points
1011   IF ((jts /= jds) .OR. (jte /= jde)) THEN
1012      WRITE ( wrf_err_message , * ) 'module_spectral_nudging: 3d: (jts /= jds) or (jte /= jde)',jts,jds,jte,jde
1013      CALL wrf_error_fatal ( TRIM( wrf_err_message ) )
1014   END IF
1015
1016
1017   nx = jde-jds+1
1018   ny = kte-kts+1 ! we can filter extra level for variables that are non-Z-staggered
1019   i_end = MIN(ite, ide-1)
1020   IF (i_end == ide-1) i_end = ide
1021   DO i = its, i_end
1022
1023         DO k=kts,kte
1024         DO j=jds,jde
1025            sheet(j-jds+1,k-kts+1) = f(i,k,j)
1026         END DO
1027            sheet(jde,k-kts+1) = 0.
1028         END DO
1029
1030         CALL spectralnudgingfilterfft2dncar(nx,ny,nwave,sheet)
1031
1032         DO k=kts,kte
1033            DO j=jds,jde
1034               f(i,k,j) = sheet(j-jds+1,k-kts+1)
1035            END DO
1036         END DO
1037   END DO ! outer i (longitude) loop
1038
1039 END SUBROUTINE spectral_nudging_filter_3dy
1040
1041 !------------------------------------------------------------------------------
1042
1043 SUBROUTINE spectralnudgingfilterfft2dncar(nx,ny,nwave,fin)
1044   IMPLICIT NONE
1045   INTEGER , INTENT(IN) :: nx, ny, nwave
1046   REAL , DIMENSION(nx,ny), INTENT(INOUT) :: fin
1047
1048   INTEGER :: i, j
1049   REAL, dimension(nx,ny) :: fp
1050
1051   INTEGER :: lensave, ier, nh, n1
1052   INTEGER :: lot, jump, n, inc, lenr, lensav, lenwrk
1053   REAL, DIMENSION(nx+15) :: wsave
1054   REAL, DIMENSION(nx,ny) :: work
1055
1056
1057
1058
1059 !  we are following the naming convention of the fftpack5 routines
1060
1061   n = nx
1062   lot = ny
1063   lensav = n+15
1064   inc = 1
1065   lenr = nx*ny
1066   jump = nx
1067   lenwrk = lenr
1068
1069 !  forward transform
1070 !  initialize coefficients, place in wsave
1071 !   (should place this in init and save wsave at program start)
1072
1073   call rfftmi(n,wsave,lensav,ier)
1074   IF(ier /= 0) THEN
1075     write(wrf_err_message,*) ' error in rfftmi ',ier
1076     CALL wrf_message(TRIM(wrf_err_message))
1077   END IF
1078
1079 !  do the forward transform
1080
1081   call rfftmf( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
1082   IF(ier /= 0) THEN
1083     write(wrf_err_message,*) ' error in rfftmf ',ier
1084     CALL wrf_message(TRIM(wrf_err_message))
1085   END IF
1086
1087       nh = min(max(1 + 2*nwave,0),n)
1088
1089
1090 ! filter all waves with wavenumber larger than nwave
1091
1092   fp = 1.
1093
1094   DO j=1,ny
1095      DO i=nh+1,n
1096          fp(i,j) = 0.
1097      ENDDO
1098   ENDDO
1099
1100   DO j=1,ny
1101     DO i=1,nx
1102       fin(i,j) = fp(i,j)*fin(i,j)
1103     ENDDO
1104   ENDDO
1105
1106 !  do the backward transform
1107
1108   call rfftmb( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
1109   IF(ier /= 0) THEN
1110     write(wrf_err_message,*) ' error in rfftmb ',ier
1111     CALL wrf_message(TRIM(wrf_err_message))
1112   END IF
1113
1114 END SUBROUTINE spectralnudgingfilterfft2dncar
1115
1116 !-------------------------------------------------------------------
1117
1118    SUBROUTINE fddaspnudginginit(id,rundgdten,rvndgdten,rthndgdten,rphndgdten,rqvndgdten, &
1119                run_hours,  &
1120                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph, &
1121                if_no_pbl_nudging_q, &
1122                if_zfac_uv, k_zfac_uv, dk_zfac_uv, &
1123                if_zfac_t, k_zfac_t, dk_zfac_t, &
1124                if_zfac_ph, k_zfac_ph, dk_zfac_ph, &
1125                if_zfac_q, k_zfac_q, dk_zfac_q, ktrop, &
1126                guv, gt, gph, gq, if_ramping, dtramp_min, end_fdda_hour, &
1127                xwavenum,ywavenum,                          &
1128                       restart, allowed_to_read,                    &
1129                       ids, ide, jds, jde, kds, kde,                &
1130                       ims, ime, jms, jme, kms, kme,                &
1131                       its, ite, jts, jte, kts, kte                 )
1132 !-------------------------------------------------------------------
1133    IMPLICIT NONE
1134 !-------------------------------------------------------------------
1135 !
1136    INTEGER , INTENT(IN)         ::  id
1137    LOGICAL, INTENT(IN)          ::  restart, allowed_to_read
1138    INTEGER, INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
1139                                     ims, ime, jms, jme, kms, kme, &
1140                                     its, ite, jts, jte, kts, kte
1141    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
1142                                                        rundgdten, &
1143                                                        rvndgdten, &
1144                                                       rthndgdten, &
1145                                                       rphndgdten, &
1146                                                       rqvndgdten
1147    INTEGER,  INTENT(IN)   ::      run_hours
1148    INTEGER,  INTENT(IN)   ::      if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
1149                                   if_no_pbl_nudging_ph, if_no_pbl_nudging_q,  end_fdda_hour
1150    INTEGER,  INTENT(IN)   ::      if_zfac_uv, if_zfac_t, if_zfac_ph, if_zfac_q
1151    INTEGER,  INTENT(IN)   ::      k_zfac_uv,  k_zfac_t,  k_zfac_ph, k_zfac_q
1152    INTEGER,  INTENT(IN)   ::      dk_zfac_uv,  dk_zfac_t,  dk_zfac_ph, dk_zfac_q
1153    INTEGER,  INTENT(IN)   ::      ktrop
1154    INTEGER,  INTENT(IN)   ::      if_ramping
1155    INTEGER,  INTENT(IN)   ::      xwavenum,ywavenum
1156    REAL,     INTENT(IN)   ::      dtramp_min
1157    REAL, INTENT(IN)       ::      guv, gt, gph, gq
1158    REAL                   ::      actual_end_fdda_min
1159
1160    INTEGER :: i, j, k
1161
1162    LOGICAL , EXTERNAL     ::      wrf_dm_on_monitor
1163
1164    CHARACTER (LEN=256) :: wrf_err_message
1165
1166    IF ( wrf_dm_on_monitor() ) THEN
1167
1168      IF( guv > 0.0) THEN
1169        WRITE(wrf_err_message,'(a,i1,a,e12.4,a,i4,a,i4)') &
1170            'D0',id,' Spectral nudging for wind is turned on and Guv= ', guv,' xwave= ',xwavenum,' ywavenum= ',ywavenum
1171        CALL wrf_message(TRIM(wrf_err_message))
1172      ELSE IF( guv < 0.0 ) THEN
1173        CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
1174      ELSE
1175        WRITE(wrf_err_message,'(a,i1,a,e12.4)') &
1176            'D0',id,' Spectral nudging for wind is turned off and Guv= ', guv
1177        CALL wrf_message(TRIM(wrf_err_message))
1178      ENDIF
1179
1180      IF( gt > 0.0) THEN
1181        WRITE(wrf_err_message,'(a,i1,a,e12.4,a,i4,a,i4)') &
1182            'D0',id,' Spectral nudging for temperature is turned on and Gt= ', gt,' xwave= ',xwavenum,' ywavenum= ',ywavenum
1183        CALL wrf_message(TRIM(wrf_err_message))
1184      ELSE IF( gt < 0.0 ) THEN
1185        CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
1186      ELSE
1187        WRITE(wrf_err_message,'(a,i1,a,e12.4)') &
1188            'D0',id,' Spectral nudging for temperature is turned off and Gt= ', gt
1189        CALL wrf_message(TRIM(wrf_err_message))
1190      ENDIF
1191
1192      IF( gph > 0.0) THEN
1193        WRITE(wrf_err_message,'(a,i1,a,e12.4,a,i4,a,i4)') &
1194          'D0',id,' Spectral nudging for geopotential is turned on and Gph= ', gph,' xwave= ',xwavenum,' ywavenum= ',ywavenum
1195        CALL wrf_message(TRIM(wrf_err_message))
1196      ELSE IF( gph < 0.0 ) THEN
1197        CALL wrf_error_fatal('In grid FDDA, Gph must be positive.')
1198      ELSE
1199        WRITE(wrf_err_message,'(a,i1,a,e12.4)') &
1200          'D0',id,' Spectral nudging for geopotential is turned off and Gph= ', gph
1201        CALL wrf_message(TRIM(wrf_err_message))
1202      ENDIF
1203
1204      IF( gq > 0.0) THEN
1205        WRITE(wrf_err_message,'(a,i1,a,e12.4,a,i4,a,i4)') &
1206          'D0',id,' Spectral nudging for water vapor mixing ratio is turned on and Gq= ', gq,' xwave= ',xwavenum,' ywavenum= ',ywavenum
1207        CALL wrf_message(TRIM(wrf_err_message))
1208      ELSE IF( gq < 0.0 ) THEN
1209        CALL wrf_error_fatal('In grid FDDA, Gq must be positive.')
1210      ELSE
1211        WRITE(wrf_err_message,'(a,i1,a,e12.4)') &
1212          'D0',id,' Spectral nudging for water vapor mixing ratio is turned off and Gq= ', gq
1213        CALL wrf_message(TRIM(wrf_err_message))
1214      ENDIF
1215
1216      IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
1217         WRITE(wrf_err_message,'(a,i1,a)') &
1218            'D0',id,' Spectral nudging for wind is turned off within the PBL.'
1219         CALL wrf_message(TRIM(wrf_err_message))
1220              IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must be greater or equal than 1.')
1221      ELSEIF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
1222         WRITE(wrf_err_message,'(a,i1,a,i3)') &
1223            'D0',id,' Spectral nudging for wind is turned off below layer', k_zfac_uv
1224         CALL wrf_message(TRIM(wrf_err_message))
1225              IF( dk_zfac_uv < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_uv must  be greater or equal than 1.')
1226      ENDIF
1227
1228
1229      IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
1230         WRITE(wrf_err_message,'(a,i1,a)') &
1231            'D0',id,' Spectral nudging for temperature is turned off within the PBL.'
1232         CALL wrf_message(TRIM(wrf_err_message))
1233              IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
1234      ELSEIF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
1235         WRITE(wrf_err_message,'(a,i1,a,i3)') &
1236            'D0',id,' Spectral nudging for temperature is turned off below layer', k_zfac_t
1237         CALL wrf_message(TRIM(wrf_err_message))
1238             IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
1239      ENDIF
1240
1241
1242      IF( gph > 0.0 .AND. if_no_pbl_nudging_ph == 1 ) THEN
1243         WRITE(wrf_err_message,'(a,i1,a)') &
1244          'D0',id,' Spectral nudging for geopotential is turned off within the PBL.'
1245         CALL wrf_message(TRIM(wrf_err_message))
1246             IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
1247      ELSEIF( gph > 0.0 .AND. if_zfac_ph == 1 ) THEN
1248         WRITE(wrf_err_message,'(a,i1,a,i3)') &
1249           'D0',id,' Spectral nudging for geopotential is turned off below layer', &
1250            k_zfac_ph
1251         CALL wrf_message(TRIM(wrf_err_message))
1252             IF( dk_zfac_ph < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_ph must be greater or equal than 1.')
1253      ENDIF
1254
1255      IF( gq > 0.0 .AND. if_no_pbl_nudging_q == 1 ) THEN
1256         WRITE(wrf_err_message,'(a,i1,a)') &
1257            'D0',id,' Spectral nudging for water vapor mixing ratio is turned off within the PBL.'
1258         CALL wrf_message(TRIM(wrf_err_message))
1259              IF( dk_zfac_t < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_t must be greater or equal than 1.')
1260      ELSEIF( gq > 0.0 .AND. if_zfac_q == 1 ) THEN
1261         WRITE(wrf_err_message,'(a,i1,a,i3)') &
1262            'D0',id,' Spectral nudging for water vapor mixing ratio is turned off below layer', k_zfac_q
1263         CALL wrf_message(TRIM(wrf_err_message))
1264             IF( dk_zfac_q < 1 ) CALL wrf_error_fatal('In spectral nudging, dk_zfac_q must be greater or equal than 1.')
1265      ENDIF
1266
1267
1268      IF( ktrop > 0 .AND. ktrop < NINT(FLOAT(kde) * 0.5) ) THEN
1269         WRITE(wrf_err_message,'(a,i3)') &
1270           'Tropopause layer for nudging restriction seems too low. KTROP set to:', ktrop
1271         CALL wrf_message(TRIM(wrf_err_message))
1272      ENDIF
1273
1274
1275      IF ( ktrop > FLOAT(kde) ) THEN
1276        CALL wrf_error_fatal('KTROP is set above the model top.')
1277      ENDIF
1278
1279
1280      IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
1281        IF( dtramp_min <= 0.0 ) THEN
1282          actual_end_fdda_min = end_fdda_hour*60.0
1283        ELSE
1284          actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
1285        ENDIF
1286
1287        IF( actual_end_fdda_min <= run_hours*60. ) THEN
1288           WRITE(wrf_err_message,'(a,i1,a)') &
1289             'D0',id,' Spectral nudging is ramped down near the end of the nudging period,'
1290           CALL wrf_message(TRIM(wrf_err_message))
1291
1292           WRITE(wrf_err_message,'(a,f6.2,a,f6.2,a)') &
1293              '      starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0,&
1294              'h, ending at ', actual_end_fdda_min/60.0,'h.'
1295           CALL wrf_message(TRIM(wrf_err_message))
1296        ENDIF
1297      ENDIF
1298
1299    ENDIF
1300
1301    IF(.not.restart) THEN
1302      DO j = jts,jte
1303      DO k = kts,kte
1304      DO i = its,ite
1305         rundgdten(i,k,j) = 0.
1306         rvndgdten(i,k,j) = 0.
1307         rthndgdten(i,k,j) = 0.
1308         rphndgdten(i,k,j) = 0.
1309         rqvndgdten(i,k,j) = 0.
1310      ENDDO
1311      ENDDO
1312      ENDDO
1313    ENDIF
1314
1315    END SUBROUTINE fddaspnudginginit
1316 !-------------------------------------------------------------------
1317
1318 END MODULE module_fdda_spnudging