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