phys/module_fddaobs_driver.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 MODULE module_fddaobs_driver
   3
   4 ! This obs-nudging FDDA module (RTFDDA) is developed by the
   5 ! NCAR/RAL/NSAP (National Security Application Programs), under the
   6 ! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
   7 ! acknowledged for releasing this capability for WRF community
   8 ! research applications.
   9 !
  10 ! The NCAR/RAL RTFDDA module was adapted, and significantly modified
  11 ! from the obs-nudging module in the standard MM5V3.1 which was originally
  12 ! developed by PSU (Stauffer and Seaman, 1994).
  13 !
  14 ! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
  15 ! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
  16 ! Nov. 2006
  17 !
  18 ! References:
  19 !
  20 !   Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
  21 !     implementation of obs-nudging-based FDDA into WRF for supporting
  22 !     ATEC test operations. 2005 WRF user workshop. Paper 10.7.
  23 !
  24 !   Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
  25 !     on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
  26 !     and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
  27
  28 !
  29 !   Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
  30 !     assimilation. J. Appl. Meteor., 33, 416-434.
  31 !
  32 !   http://www.rap.ucar.edu/projects/armyrange/references.html
  33 !
  34
  35 CONTAINS
  36
  37 !-----------------------------------------------------------------------
  38 SUBROUTINE fddaobs_driver( inest, domid, parid, restart,         &
  39                config_flags,                                     &
  40                nudge_opt, iprt_errob, iprt_nudob,                &
  41                fdasta, fdaend,                                   &
  42                nudge_wind, nudge_temp, nudge_mois,               &
  43                nudge_pstr,                                       &
  44                coef_wind, coef_temp, coef_mois,                  &
  45                coef_pstr, rinxy, rinsig,                         &
  46                npfi, ionf,                                       &
  47                obs_prt_max, obs_prt_freq, idynin, dtramp,        &
  48                parent_grid_ratio, maxdom, itimestep,             &
  49                xtime,                                            &
  50                dt, gmt, julday,                                  &
  51 #if ( EM_CORE == 1 )
  52                fdob,                                             &
  53 #endif
  54                max_obs, nobs_ndg_vars,                           &
  55                nobs_err_flds, nstat, varobs, errf, dx,           &
  56                KPBL, HT, mut, muu, muv, c1h, c2h,                &
  57                msftx, msfty, msfux, msfuy, msfvx, msfvy, p_phy, t_tendf, t0,             &
  58                ub, vb, tb, qvb, pbase, ptop, pp, phb, ph,        &
  59                uratx, vratx, tratx, ru_tendf, rv_tendf,          &
  60                moist_tend, savwt,                                &
  61                regime, pblh, z_at_w,                             &
  62                z,                                                &
  63                ids,ide, jds,jde, kds,kde,                        & ! domain dims
  64                ims,ime, jms,jme, kms,kme,                        & ! memory dims
  65                its,ite, jts,jte, kts,kte                         ) ! tile   dims
  66
  67 !-----------------------------------------------------------------------
  68   USE module_domain
  69   ! USE module_bc  seems to not be needed
  70   USE module_model_constants, ONLY : g, rcp
  71   USE module_fddaobs_rtfdda
  72
  73 ! This driver calls subroutines for fdda obs-nudging and
  74 ! returns computed tendencies
  75
  76 !
  77 !-----------------------------------------------------------------------
  78   IMPLICIT NONE
  79 !-----------------------------------------------------------------------
  80 ! taken from MM5 code - 03 Feb 2004.
  81 !-----------------------------------------------------------------------
  82
  83 !=======================================================================
  84 ! Definitions
  85 !-----------
  86 !-- KPBL          vertical layer index for PBL top
  87 !-- HT            terrain height (m)
  88 !-- p_phy         pressure (Pa)
  89 !-- t_tendf       temperature tendency
  90
  91   INTEGER, intent(in)  :: ids,ide, jds,jde, kds,kde  ! domain dims.
  92   INTEGER, intent(in)  :: ims,ime, jms,jme, kms,kme  ! memory dims.
  93   INTEGER, intent(in)  :: its,ite, jts,jte, kts,kte  ! tile   dims.
  94
  95   INTEGER, intent(in)  :: inest
  96   INTEGER, intent(in)  :: maxdom
  97   INTEGER, intent(in)  :: domid(maxdom)           ! Domain IDs
  98   INTEGER, intent(in)  :: parid(maxdom)           ! Parent domain IDs
  99   LOGICAL, intent(in)  :: restart
 100   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
 101   INTEGER, intent(in)  :: itimestep
 102   INTEGER, intent(in)  :: nudge_opt
 103   LOGICAL, intent(in)  :: iprt_errob
 104   LOGICAL, intent(in)  :: iprt_nudob
 105   REAL, intent(in)     :: fdasta
 106   REAL, intent(in)     :: fdaend
 107   INTEGER, intent(in)  :: nudge_wind
 108   INTEGER, intent(in)  :: nudge_temp
 109   INTEGER, intent(in)  :: nudge_mois
 110   INTEGER, intent(in)  :: nudge_pstr
 111   REAL, intent(in) :: coef_wind
 112   REAL, intent(in) :: coef_temp
 113   REAL, intent(in) :: coef_mois
 114   REAL, intent(in) :: coef_pstr
 115   REAL, intent(inout)  :: rinxy
 116   REAL, intent(inout)  :: rinsig
 117   INTEGER, intent(in) :: npfi
 118   INTEGER, intent(in) :: ionf
 119   INTEGER, intent(in) :: obs_prt_max      ! max number of obs in printout
 120   INTEGER, intent(in) :: obs_prt_freq     ! frequency (in obs index) printout
 121   INTEGER, intent(in) :: idynin
 122   REAL, intent(inout) :: dtramp
 123   INTEGER, intent(in) :: parent_grid_ratio
 124   REAL, intent(in)     :: xtime           ! forecast time in minutes
 125   REAL, intent(in)     :: dt
 126   REAL, intent(in)     :: gmt
 127   INTEGER, intent(in)  :: julday
 128   INTEGER, intent(in)  :: max_obs         ! max number of observations
 129   INTEGER, intent(in)  :: nobs_ndg_vars
 130   INTEGER, intent(in)  :: nobs_err_flds
 131   INTEGER, intent(in)  :: nstat
 132   REAL, intent(inout)  :: varobs(nobs_ndg_vars, max_obs)
 133   REAL, intent(inout)  :: errf(nobs_err_flds, max_obs)
 134   REAL, intent(in)     :: dx           ! this-domain grid cell-size (m)
 135   INTEGER, INTENT(IN) :: kpbl( ims:ime, jms:jme ) ! vertical layer index for PBL top
 136   REAL, INTENT(IN) :: ht( ims:ime, jms:jme )
 137   REAL, INTENT(IN) :: mut( ims:ime , jms:jme )   ! Air mass on t-grid
 138   REAL, INTENT(IN) :: muu( ims:ime , jms:jme )   ! Air mass on u-grid
 139   REAL, INTENT(IN) :: muv( ims:ime , jms:jme )   ! Air mass on v-grid
 140   REAL, INTENT(IN) :: msftx( ims:ime , jms:jme )  ! Map scale on t-grid
 141   REAL, INTENT(IN) :: msfty( ims:ime , jms:jme )  ! Map scale on t-grid
 142   REAL, INTENT(IN) :: msfux( ims:ime , jms:jme )  ! Map scale on u-grid
 143   REAL, INTENT(IN) :: msfuy( ims:ime , jms:jme )  ! Map scale on u-grid
 144   REAL, INTENT(IN) :: msfvx( ims:ime , jms:jme )  ! Map scale on v-grid
 145   REAL, INTENT(IN) :: msfvy( ims:ime , jms:jme )  ! Map scale on v-grid
 146
 147   REAL, INTENT(IN) :: c1h( kms:kme )   ! Hybrid coordinate, weighting function
 148   REAL, INTENT(IN) :: c2h( kms:kme )   ! Hybrid coordinate, weighting function
 149
 150   REAL, INTENT(IN) :: p_phy( ims:ime, kms:kme, jms:jme )
 151   REAL, INTENT(INOUT) :: t_tendf( ims:ime, kms:kme, jms:jme )
 152   REAL, INTENT(IN) :: t0
 153   REAL, INTENT(INOUT) :: savwt( nobs_ndg_vars, ims:ime, kms:kme, jms:jme )
 154   REAL, INTENT(INOUT) :: regime( ims:ime, jms:jme )
 155   REAL, INTENT(IN) :: pblh( ims:ime, jms:jme )
 156   REAL, INTENT(IN) :: z_at_w( ims:ime, kms:kme, jms:jme ) ! Model ht(m) asl, f-levs
 157   REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme )  ! Model ht(m) asl, h-levs
 158
 159 #if ( EM_CORE == 1 )
 160   TYPE(fdob_type), intent(inout)  :: fdob
 161 #endif
 162
 163   REAL,   INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
 164   REAL,   INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
 165   REAL,   INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
 166   REAL,   INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
 167   REAL,   INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme ) ! Base press. (Pa)
 168   REAL,   INTENT(IN) :: ptop
 169   REAL,   INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme )  ! Press. pert. (Pa)
 170   REAL,   INTENT(IN) :: phb( ims:ime, kms:kme, jms:jme ) ! Base geopotential
 171   REAL,   INTENT(IN) :: ph( ims:ime, kms:kme, jms:jme )  ! Perturbation geopotential
 172   REAL,   INTENT(IN) :: uratx( ims:ime, jms:jme )     ! On mass points
 173   REAL,   INTENT(IN) :: vratx( ims:ime, jms:jme )     !       "
 174   REAL,   INTENT(IN) :: tratx( ims:ime, jms:jme )     !       "
 175   REAL,   INTENT(INOUT) :: ru_tendf( ims:ime, kms:kme, jms:jme )
 176   REAL,   INTENT(INOUT) :: rv_tendf( ims:ime, kms:kme, jms:jme )
 177   REAL,   INTENT(INOUT) :: moist_tend( ims:ime, kms:kme, jms:jme )
 178
 179 ! Local variables
 180   logical            :: nudge_flag   ! Flag for doing nudging
 181   integer            :: KTAU         ! Forecast timestep
 182   real               :: dtmin        ! dt in minutes
 183   integer            :: i, j, k      ! Loop counters.
 184   integer            :: idom         ! Loop counter.
 185   integer            :: nsta         ! Number of observation stations
 186   integer            :: infr         ! Frequency for obs input and error calc
 187   integer            :: idarst       ! Flag for calling sub errob on restart
 188   real               :: dtr          ! Abs value of dtramp (for dynamic init)
 189   real               :: tconst       ! Reciprocal of dtr
 190   real    :: vih_uv(its:ite,jts:jte,2) ! Vert infl heights abv grd for LML obs (wind)
 191   real    :: vih_t (its:ite,jts:jte,2) ! Vert infl heights abv grd for LML obs (temp)
 192   real    :: vih_q (its:ite,jts:jte,2) ! Vert infl heights abv grd for LML obs (mois)
 193   integer :: vik_uv(its:ite,jts:jte,2) ! Vert infl k-levels for LML obs (wind)
 194   integer :: vik_t (its:ite,jts:jte,2) ! Vert infl k-levels for LML obs (temp)
 195   integer :: vik_q (its:ite,jts:jte,2) ! Vert infl k-levels for LML obs (mois)
 196   real    :: z_at_p( kms:kme )       ! Height at p levels
 197 #ifdef RAL
 198   real    :: HTIJ(ids:ide, jds:jde) = 0.  ! Terrain ht on global grid
 199 #endif
 200   character(len=200) :: msg  ! Argument to wrf_message
 201
 202 #if ( EM_CORE == 1 )
 203   nudge_flag = (nudge_opt  .eq. 1)
 204
 205   if (.not. nudge_flag) return
 206
 207 !----------------------------------------------------------------------
 208 ! ***************       BEGIN FDDA SETUP SECTION        ***************
 209
 210 ! Calculate forecast time.
 211   dtmin = dt/60.
 212   ktau  = itimestep - 1        !ktau corresponds to xtime
 213
 214 ! Set NSTA to zero on startup, or else retrieve value from last pass.
 215   IF(ktau.EQ.fdob%ktaur) THEN
 216      if (iprt_nudob) then
 217         write(msg,'(a,i2,a)') 'OBS NUDGING is requested on a total of ',   &
 218                               fdob%domain_tot,' domain(s).'
 219         call wrf_message(msg)
 220      endif
 221      nsta=0.
 222   ELSE
 223      nsta=fdob%nstat
 224   ENDIF
 225
 226   infr = ionf*(parent_grid_ratio**fdob%levidn(inest))
 227   nsta=fdob%nstat
 228   idarst = 0
 229   IF(restart .AND. ktau.EQ.fdob%ktaur) idarst=1
 230
 231   CALL wrf_debug(100,'in PSU FDDA scheme')
 232
 233 ! Make sure regime array is set over entire grid
 234 ! (ajb: Copied code from fddagd)
 235     IF( config_flags%sf_sfclay_physics /= sfclayscheme &
 236   .AND. config_flags%sf_sfclay_physics /= mynnsfcscheme &
 237   .AND. config_flags%sf_sfclay_physics /= pxsfcscheme &
 238   .AND. config_flags%sf_sfclay_physics /= sfclayrevscheme ) THEN
 239       DO j = jts, jte
 240       DO i = its, ite
 241            IF( pblh(i,j) > z_at_w(i,2,j)-ht(i,j) ) THEN
 242              regime(i,j) = 4.0
 243            ELSE
 244              regime(i,j) = 1.0
 245            ENDIF
 246       ENDDO
 247       ENDDO
 248     ENDIF
 249
 250 ! Compute VIF heights for each grid column (used for LML obs)
 251    if(fdob%sfc_scheme_vert.EQ.0) then
 252      if(nudge_wind.EQ.1 .AND. NSTA.GT.0)  then
 253        CALL compute_VIH( fdob%vif_uv, fdob%vif_max,                &
 254                          fdob%vif_fullmin, fdob%vif_rampmin,       &
 255                          regime, pblh,                             &
 256                          ht, z, vih_uv,                            &
 257                          ids,ide, jds,jde, kds,kde,                &
 258                          ims,ime, jms,jme, kms,kme,                &
 259                          its,ite, jts,jte, kts,kte )
 260      endif
 261      if(nudge_temp.EQ.1 .AND. NSTA.GT.0)  then
 262        CALL compute_VIH( fdob%vif_t, fdob%vif_max,                 &
 263                          fdob%vif_fullmin, fdob%vif_rampmin,       &
 264                          regime, pblh,                             &
 265                          ht, z, vih_t,                             &
 266                          ids,ide, jds,jde, kds,kde,                &
 267                          ims,ime, jms,jme, kms,kme,                &
 268                          its,ite, jts,jte, kts,kte )
 269      endif
 270      if(nudge_mois.EQ.1 .AND. NSTA.GT.0)  then
 271        CALL compute_VIH( fdob%vif_q, fdob%vif_max,                 &
 272                          fdob%vif_fullmin, fdob%vif_rampmin,       &
 273                          regime, pblh,                             &
 274                          ht, z, vih_q,                             &
 275                          ids,ide, jds,jde, kds,kde,                &
 276                          ims,ime, jms,jme, kms,kme,                &
 277                          its,ite, jts,jte, kts,kte )
 278      endif
 279    endif
 280 !********************* END AJB MOVE TO SLAB ***************************
 281
 282 ! COMPUTE ERROR BETWEEN OBSERVATIONS and MODEL
 283   IF( nsta.GT.0 ) THEN
 284     IF( MOD(ktau,infr).EQ.0 .OR. idarst.EQ.1) THEN
 285
 286         CALL errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rcp,       &
 287                    z,                                                &
 288                    uratx, vratx, tratx, kpbl,                        &
 289                    nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,    &
 290                    fdob%levidn, parid, fdob%nstat, fdob%nstaw,       &
 291                    nudge_wind, nudge_temp, nudge_mois, nudge_pstr,   &
 292                    fdob%timeob, fdob%rio, fdob%rjo, fdob%rko,        &
 293                    varobs, errf, ktau, xtime,                        &
 294                    parent_grid_ratio, npfi,                          &
 295                    obs_prt_max, obs_prt_freq, iprt_errob,            &
 296                    fdob%obsprt, fdob%stnidprt,                       &
 297                    fdob%latprt, fdob%lonprt,                         &
 298                    fdob%mlatprt, fdob%mlonprt,                       &
 299                    ids,ide, jds,jde, kds,kde,                        &
 300                    ims,ime, jms,jme, kms,kme,                        &
 301                    its,ite, jts,jte, kts,kte)
 302     ENDIF
 303   ENDIF
 304
 305   fdob%tfaci=1.0
 306   IF(idynin.EQ.1.AND.nudge_opt.EQ.1) THEN
 307     dtr=ABS(dtramp)
 308     tconst=1./dtr
 309 ! FDAEND(IN) IS THE TIME IN MINUTES TO END THE DYNAMIC INITIALIZATION CY
 310     IF(xtime.LT.fdaend-dtr)THEN
 311       fdob%tfaci=1.
 312     ELSEIF(xtime.GE.fdaend-dtr.AND.xtime.LE.fdaend) THEN
 313       fdob%tfaci=(fdaend-xtime)*tconst
 314     ELSE
 315       fdob%tfaci=0.0
 316     ENDIF
 317     IF(ktau.EQ.fdob%ktaur.OR.MOD(ktau,10).EQ.0) THEN
 318       IF (iprt_nudob)                                                  &
 319          PRINT*,' DYNINOBS: IN,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST',   &
 320          ',TFACI: ',INEST,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST,         &
 321          fdob%TFACI
 322     ENDIF
 323   ENDIF
 324
 325 #ifdef RAL
 326 ! MEIXU: collect terrain array HT into a global array HTIJ
 327   CALL loc2glob(HT, HTIJ, "2D", "REAL",                  &
 328                 ids,ide, jds,jde, kds,kde,               &
 329                 ims,ime, jms,jme, kms,kme )
 330 ! MEIXU end
 331 #endif
 332 !
 333 ! ***************        END FDDA SETUP SECTION         ***************
 334 !----------------------------------------------------------------------
 335
 336 !----------------------------------------------------------------------
 337 ! ***************         BEGIN NUDGING SECTION         ***************
 338
 339   DO J = jts, jte
 340 !
 341 ! IF NUDGING SURFACE WINDS IN THE BOUNDARY LAYER, IF IWINDS(INEST+2)=1
 342 ! USE A SIMILARITY CORRECTION BASED ON ROUGHNESS TO APPLY 10M
 343 ! WIND TO THE SURFACE LAYER (K=KL) AT 40M.  TO DO THIS WE MUST
 344 ! STORE ROUGHNESS AND REGIME FOR EACH J SLICE AFTER THE CALL TO
 345 ! HIRPBL FOR LATER USE IN BLNUDGD.
 346 !
 347 !--- OBS NUDGING FOR TEMP AND MOISTURE
 348 !
 349      NSTA=NSTAT
 350      IF(J .GT. 2 .and. J .LT. jde-1) THEN
 351        IF(nudge_temp.EQ.1 .AND. NSTA.GT.0)  &
 352        THEN
 353 !         write(6,*) 'calling nudob: IVAR=3, J = ',j
 354           CALL nudob(J, 3, t_tendf(ims,kms,j),                       &
 355                   inest, restart, ktau, fdob%ktaur, xtime,           &
 356                   mut(ims,j), c1h, c2h, msftx(ims,j), msfty(ims,j),  &
 357                   nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
 358                   npfi, ionf, rinxy, fdob%window,                    &
 359                   fdob%nudge_t_pbl,                                  &
 360                   fdob%sfcfact, fdob%sfcfacr,                        &
 361                   fdob%levidn,                                       &
 362                   parid, nstat,                                      &
 363                   fdob%rinfmn, fdob%rinfmx, fdob%pfree,              &
 364                   fdob%dcon, fdob%tfaci,                             &
 365                   fdob%sfc_scheme_horiz, fdob%sfc_scheme_vert,       &
 366                   fdob%max_sndng_gap,                                &
 367                   fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,          &
 368                   parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
 369                   fdob%rko, fdob%timeob, varobs, errf,               &
 370                   pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
 371                   nudge_wind, nudge_temp, nudge_mois,                &
 372                   coef_wind, coef_temp, coef_mois,                   &
 373                   savwt(1,ims,kms,j), kpbl(ims,j), 0,                &
 374                   vih_t(its,j,1), vih_t(its,j,2), ht(ims,j),         &
 375                   z(ims,kms,j),                                      &
 376                   iprt_nudob,                                        &
 377                   ids,ide, jds,jde, kds,kde,                         & ! domain dims
 378                   ims,ime, jms,jme, kms,kme,                         & ! memory dims
 379                   its,ite, jts,jte, kts,kte,                         & ! tile   dims
 380                   qvb, config_flags%obs_scl_neg_qv_innov ) ! Water vapor mixing ratio / scale negative qv innovations
 381 !         write(6,*) 'return from nudob: IVAR=3, J = ',j
 382        ENDIF
 383
 384        IF(nudge_mois.EQ.1 .AND. NSTA.GT.0)  &
 385        THEN
 386 !         write(6,*) 'calling nudob: IVAR=4, J = ',j
 387           CALL nudob(J, 4, moist_tend(ims,kms,j),                    &
 388                   inest, restart, ktau, fdob%ktaur, xtime,           &
 389                   mut(ims,j), c1h, c2h, msftx(ims,j), msfty(ims,j),  &
 390                   nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
 391                   npfi, ionf, rinxy, fdob%window,                    &
 392                   fdob%nudge_q_pbl,                                  &
 393                   fdob%sfcfact, fdob%sfcfacr,                        &
 394                   fdob%levidn,                                       &
 395                   parid, nstat,                                      &
 396                   fdob%rinfmn, fdob%rinfmx, fdob%pfree,              &
 397                   fdob%dcon, fdob%tfaci,                             &
 398                   fdob%sfc_scheme_horiz, fdob%sfc_scheme_vert,       &
 399                   fdob%max_sndng_gap,                                &
 400                   fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,          &
 401                   parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
 402                   fdob%rko, fdob%timeob, varobs, errf,               &
 403                   pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
 404                   nudge_wind, nudge_temp, nudge_mois,                &
 405                   coef_wind, coef_temp, coef_mois,                   &
 406                   savwt(1,ims,kms,j), kpbl(ims,j), 0,                &
 407                   vih_q(its,j,1), vih_q(its,j,2), ht(ims,j),         &
 408                   z(ims,kms,j),                                      &
 409                   iprt_nudob,                                        &
 410                   ids,ide, jds,jde, kds,kde,                         & ! domain dims
 411                   ims,ime, jms,jme, kms,kme,                         & ! memory dims
 412                   its,ite, jts,jte, kts,kte,                         & ! tile   dims
 413                   qvb, config_flags%obs_scl_neg_qv_innov ) ! Water vapor mixing ratio / scale negative qv innovations
 414 !         write(6,*) 'return from nudob: IVAR=4, J = ',j
 415        ENDIF
 416      ENDIF
 417
 418      IF(nudge_wind.EQ.1 .AND. NSTA.GT.0)    &
 419      THEN
 420 !         write(6,*) 'calling nudob: IVAR=1, J = ',j
 421         CALL nudob(J, 1, ru_tendf(ims,kms,j),                        &
 422                 inest, restart, ktau, fdob%ktaur, xtime,             &
 423                 muu(ims,j), c1h, c2h, msfux(ims,j), msfuy(ims,j),    &
 424                 nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
 425                 npfi, ionf, rinxy, fdob%window,                      &
 426                 fdob%nudge_uv_pbl,                                   &
 427                 fdob%sfcfact, fdob%sfcfacr,                          &
 428                 fdob%levidn,                                         &
 429                 parid, nstat,                                        &
 430                 fdob%rinfmn, fdob%rinfmx, fdob%pfree,                &
 431                 fdob%dcon, fdob%tfaci,                               &
 432                 fdob%sfc_scheme_horiz, fdob%sfc_scheme_vert,         &
 433                 fdob%max_sndng_gap,                                  &
 434                 fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,            &
 435                 parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
 436                 fdob%rko, fdob%timeob, varobs, errf,                 &
 437                 pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
 438                 nudge_wind, nudge_temp, nudge_mois,                  &
 439                 coef_wind, coef_temp, coef_mois,                     &
 440                 savwt(1,ims,kms,j), kpbl(ims,j), 0,                  &
 441                 vih_uv(its,j,1), vih_uv(its,j,2), ht(ims,j),         &
 442                 z(ims,kms,j),                                        &
 443                 iprt_nudob,                                          &
 444                 ids,ide, jds,jde, kds,kde,                           & ! domain dims
 445                 ims,ime, jms,jme, kms,kme,                           & ! memory dims
 446                 its,ite, jts,jte, kts,kte,                           & ! tile   dims
 447                 qvb, config_flags%obs_scl_neg_qv_innov ) ! Water vapor mixing ratio / scale negative qv innovations
 448 !       write(6,*) 'return from nudob: IVAR=1, J = ',j
 449
 450 !       write(6,*) 'calling nudob: IVAR=2, J = ',j
 451         CALL nudob(J, 2, rv_tendf(ims,kms,j),                        &
 452                 inest, restart, ktau, fdob%ktaur, xtime,             &
 453                 muv(ims,j), c1h, c2h, msfvx(ims,j), msfvy(ims,j),    &
 454                 nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
 455                 npfi, ionf, rinxy, fdob%window,                      &
 456                 fdob%nudge_uv_pbl,                                   &
 457                 fdob%sfcfact, fdob%sfcfacr,                          &
 458                 fdob%levidn,                                         &
 459                 parid, nstat,                                        &
 460                 fdob%rinfmn, fdob%rinfmx, fdob%pfree,                &
 461                 fdob%dcon, fdob%tfaci,                               &
 462                 fdob%sfc_scheme_horiz, fdob%sfc_scheme_vert,         &
 463                 fdob%max_sndng_gap,                                  &
 464                 fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,            &
 465                 parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
 466                 fdob%rko, fdob%timeob, varobs, errf,                 &
 467                 pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
 468                 nudge_wind, nudge_temp, nudge_mois,                  &
 469                 coef_wind, coef_temp, coef_mois,                     &
 470                 savwt(1,ims,kms,j), kpbl(ims,j), 0,                  &
 471                 vih_uv(its,j,1), vih_uv(its,j,2), ht(ims,j),         &
 472                 z(ims,kms,j),                                        &
 473                 iprt_nudob,                                          &
 474                 ids,ide, jds,jde, kds,kde,                           & ! domain dims
 475                 ims,ime, jms,jme, kms,kme,                           & ! memory dims
 476                 its,ite, jts,jte, kts,kte,                           & ! tile   dims
 477                 qvb, config_flags%obs_scl_neg_qv_innov ) ! Water vapor mixing ratio / scale negative qv innovations
 478 !       write(6,*) 'return from nudob: IVAR=2, J = ',j
 479      ENDIF
 480   ENDDO
 481 !
 482 ! --- END OF 4DDA
 483 !
 484   RETURN
 485 #endif
 486   END SUBROUTINE fddaobs_driver
 487
 488   SUBROUTINE compute_VIH(vif, hmax, fullmin, rampmin,       &
 489                          regime, pblh, terrh, z, vih,       &
 490                          ids,ide, jds,jde, kds,kde,         & ! domain dims
 491                          ims,ime, jms,jme, kms,kme,         &
 492                          its,ite, jts,jte, kts,kte)
 493
 494   USE module_fddaobs_rtfdda
 495
 496   IMPLICIT NONE
 497 !*******************************************************************************
 498 !*****     COMPUTE HEIGHTS FOR SURFACE OBS VERTICAL INFLUENCE FUNCTION     *****
 499 !*******************************************************************************
 500
 501   REAL,    INTENT(IN)  :: vif(6)                     ! Vert infl params for regimes
 502   REAL,    INTENT(IN)  :: hmax                       ! Max height to apply nudging
 503   REAL,    INTENT(IN)  :: fullmin                    ! Min height of full nudging
 504   REAL,    INTENT(IN)  :: rampmin                    ! Min height to ramp full-to-0
 505   REAL,    INTENT(IN)  :: regime(ims:ime,jms:jme)    ! Stability regime
 506   REAL,    INTENT(IN)  :: pblh(ims:ime,jms:jme)      ! PBL height (m)
 507   REAL,    INTENT(IN)  :: terrh(ims:ime,jms:jme)     ! Terrain ht (m)
 508   REAL,    INTENT(IN)  :: z(ims:ime,kms:kme,jms:jme) ! Ht (m) above sl (half levs)
 509   REAL,    INTENT(OUT) :: vih(its:ite,jts:jte,2)     ! Vt infl hts abv grd for LML obs
 510 ! INTEGER, INTENT(OUT) :: vik(its:ite,jts:jte,2)     ! Vert infl k-levels for LML obs
 511   INTEGER, INTENT(IN)  :: ids,ide, jds,jde, kds,kde  ! domain dims.
 512   INTEGER, INTENT(IN)  :: ims,ime, jms,jme, kms,kme  ! memory dims.
 513   INTEGER, INTENT(IN)  :: its,ite, jts,jte, kts,kte  ! tile   dims.
 514
 515 ! Local variables
 516   real  :: fullr(its:ite)    ! Height up to full vertical weighting
 517   real  :: rampr(its:ite)    ! Height of ramp-down to zero weighting
 518   character(len=200) :: msg  ! Argument to wrf_message
 519   integer:: i, j             ! Loop counters
 520
 521   integer k     ! ajb test
 522
 523 ! Do J-slabs
 524   do j = jts, jte
 525
 526 !   Set fullr and rampr values according to regimes
 527     do i = its, ite
 528
 529       if(regime(i,j).eq.1.0) then        ! REGIME 1
 530         fullr(i) = vif(1)
 531         rampr(i) = vif(2)
 532       elseif(regime(i,j).eq.2.0) then    ! REGIME 2
 533         fullr(i) = vif(3)
 534         rampr(i) = vif(4)
 535       elseif(regime(i,j).eq.3.0 .or. regime(i,j).eq.4.0) then    ! REGIME 4
 536         fullr(i) = vif(5)
 537         rampr(i) = vif(6)
 538       else
 539         write(msg,'(a,f5.1,2(a,i4))') 'Unknown regime type ', regime(i,j),    &
 540                                  ' at grid coordinate i = ',i,' j = ',j
 541         call wrf_message(msg)
 542         call wrf_error_fatal ( 'fddaobs_driver: compute_VIH STOP' )
 543
 544       endif
 545
 546     enddo
 547
 548 !   Get vert infl heights for LML obs, from fullr, rampr, and pblh
 549     CALL get_vif_hts_slab(fullr, rampr, pblh(ims,j),         &
 550                           hmax, fullmin, rampmin,            &
 551                           vih(its,j,1), vih(its,j,2),        &
 552                           ims,ime, its,ite)
 553   enddo
 554   END SUBROUTINE compute_VIH
 555
 556   SUBROUTINE get_vif_hts_slab(fullr, rampr, pblh, hmax, fullmin, rampmin, &
 557                               ht1, ht2, ims,ime, its,ite)
 558 ! Compute VIF heights
 559
 560   IMPLICIT NONE
 561
 562   REAL, INTENT(IN)    :: fullr(its:ite)   ! Height up to full vertical weighting
 563   REAL, INTENT(IN)    :: rampr(its:ite)   ! Height of ramp-down to zero weighting
 564   REAL, INTENT(IN)    :: pblh(ims:ime)    ! PBL height (m)
 565   REAL, INTENT(IN)    :: hmax             ! Max height to apply nudging
 566   REAL, INTENT(IN)    :: fullmin          ! Min height of full nudging
 567   REAL, INTENT(IN)    :: rampmin          ! Min height to ramp full-to-0
 568   REAL, INTENT(OUT)   :: ht1(its:ite)     ! Vert infl fcn height 1
 569   REAL, INTENT(OUT)   :: ht2(its:ite)     ! Vert infl fcn height 2
 570   INTEGER, INTENT(IN) :: ims,ime          ! Memory dims.
 571   INTEGER, INTENT(IN) :: its,ite          ! Tile   dims.
 572
 573 ! Local variables
 574   integer :: i
 575
 576   do i = its, ite
 577
 578 ! Determine lower height (below which the VIF=1 for full weighting)
 579     if(fullr(i).ge.0.0) then        ! fullr is height from ground
 580       ht1(i) = fullr(i)
 581     else                            ! fullr is relative to pbl (-5000 bias)
 582       ht1(i) = pblh(i) - (fullr(i)+5000.)
 583     endif
 584
 585 !   Height ht1 can be no smaller than fullmin
 586     ht1(i) = max(fullmin,ht1(i))
 587
 588 ! Determine upper height (to which the VIF ramps down to zero weighting)
 589 ! NOTE: Height of ramp-down (ht2-ht1) can be no smaller than rampmin
 590
 591     if(rampr(i).ge.0.0) then
 592 !     rampr is height from ht1
 593       ht2(i) = ht1(i) + max(rampmin,rampr(i))
 594     else
 595 !     rampr is relative to pbl (-5000 bias)
 596       ht2(i) = max( ht1(i)+rampmin, pblh(i)-(rampr(i)+5000.) )
 597     endif
 598
 599 ! Apply hmax
 600     ht1(i) = min(ht1(i), hmax-rampmin)
 601     ht2(i) = min(ht2(i), hmax)
 602   enddo
 603   END SUBROUTINE get_vif_hts_slab
 604
 605   SUBROUTINE get_vik_slab( h, hlevs, ht, vik, ims,ime, kms,kme, its,ite, kts,kte )
 606
 607 ! Compute VIK values from heights on j-slab
 608
 609   IMPLICIT NONE
 610
 611   REAL,    INTENT(IN) :: h(its:ite)          ! height (m) above ground, on j-slab
 612   REAL,    INTENT(IN) :: hlevs(ims:ime,kms:kme) ! hgt (m) abv grd at modl levs (slab)
 613   REAL,    INTENT(IN) :: ht(ims:ime)         ! terrain height (m) (slab)
 614   INTEGER, INTENT(OUT):: vik(its:ite)        ! vert infl k levels (slab)
 615   INTEGER, INTENT(IN) :: ims,ime, kms,kme    ! memory dims
 616   INTEGER, INTENT(IN) :: its,ite, kts,kte    ! tile   dims
 617
 618 ! Local variables
 619   integer :: i
 620   integer :: k
 621   real    :: ht_ag(kts:kte)
 622
 623   do i = its, ite
 624
 625 !   Get column of height-above-ground values for this i coord
 626     do k = kts,kte
 627        ht_ag(k) = hlevs(i,k) - ht(i)
 628     enddo
 629 !   Get k levels that correspond to height values
 630     vik(i) = ht_to_k( h(i), ht_ag, kts,kte )
 631   enddo
 632   END SUBROUTINE get_vik_slab
 633
 634   INTEGER FUNCTION ht_to_k( h, hlevs, kts,kte )
 635   IMPLICIT NONE
 636
 637   REAL,    INTENT(IN)  :: h                     ! height value (m)
 638   REAL,    INTENT(IN)  :: hlevs(kts:kte)        ! model height levels
 639   INTEGER, INTENT(IN)  :: kts,kte               ! tile dims
 640
 641 ! Local variables
 642   INTEGER :: k               ! loop counter
 643   INTEGER :: klo             ! lower k bound
 644
 645   KLEVS: do k = kts, kte
 646     klo = k-1
 647     if(h .le. hlevs(k)) then
 648       EXIT KLEVS
 649     endif
 650   enddo KLEVS
 651   klo = max0(1,klo)
 652   ht_to_k = min0(kte,klo)
 653   RETURN
 654   END FUNCTION ht_to_k
 655
 656 END MODULE module_fddaobs_driver