wrfv2_fire/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                nudge_opt, iprt_errob, iprt_nudob,                &
  40                fdasta, fdaend,                                   &
  41                nudge_wind, nudge_temp, nudge_mois,               &
  42                nudge_pstr,                                       &
  43                coef_wind, coef_temp, coef_mois,                  &
  44                coef_pstr, rinxy, rinsig,                         &
  45                npfi, ionf, nobs_prt, idynin, dtramp,             &
  46                xlatc_cg, xlonc_cg, true_lat1, true_lat2,         &
  47                map_proj, i_parent_start, j_parent_start,         &
  48                parent_grid_ratio, maxdom, itimestep,             &
  49                dt, gmt, julday,                                  &
  50 #if ( EM_CORE == 1 )
  51                fdob,                                             &
  52 #endif
  53                max_obs, nobs_ndg_vars,    &
  54                nobs_err_flds, nstat, varobs, errf, dx,           &
  55                KPBL, HT, mut, muu, muv,                          &
  56                msftx, msfty, msfux, msfuy, msfvx, msfvy, p_phy, t_tendf, t0,             &
  57                ub, vb, tb, qvb, pbase, ptop, pp,                 &
  58                uratx, vratx, tratx, ru_tendf, rv_tendf,          &
  59                moist_tend, savwt,                                &
  60                ids,ide, jds,jde, kds,kde,                        & ! domain dims
  61                ims,ime, jms,jme, kms,kme,                        & ! memory dims
  62                its,ite, jts,jte, kts,kte                         ) ! tile   dims
  63
  64 !-----------------------------------------------------------------------
  65   USE module_domain
  66   USE module_bc
  67   USE module_model_constants, ONLY : rovg, rcp
  68   USE module_fddaobs_rtfdda
  69
  70 ! This driver calls subroutines for fdda obs-nudging and
  71 ! returns computed tendencies
  72
  73 !
  74 !-----------------------------------------------------------------------
  75   IMPLICIT NONE
  76 !-----------------------------------------------------------------------
  77 ! taken from MM5 code - 03 Feb 2004.
  78 !-----------------------------------------------------------------------
  79
  80 !=======================================================================
  81 ! Definitions
  82 !-----------
  83 !-- KPBL          vertical layer index for PBL top
  84 !-- HT            terrain height (m)
  85 !-- p_phy         pressure (Pa)
  86 !-- t_tendf       temperature tendency
  87
  88   INTEGER, intent(in)  :: ids,ide, jds,jde, kds,kde  ! domain dims.
  89   INTEGER, intent(in)  :: ims,ime, jms,jme, kms,kme  ! memory dims.
  90   INTEGER, intent(in)  :: its,ite, jts,jte, kts,kte  ! tile   dims.
  91
  92   INTEGER, intent(in)  :: inest
  93   INTEGER, intent(in)  :: maxdom
  94   INTEGER, intent(in)  :: domid(maxdom)           ! Domain IDs
  95   INTEGER, intent(in)  :: parid(maxdom)           ! Parent domain IDs
  96   LOGICAL, intent(in)  :: restart
  97   INTEGER, intent(in)  :: itimestep
  98   INTEGER, intent(in)  :: nudge_opt
  99   LOGICAL, intent(in)  :: iprt_errob
 100   LOGICAL, intent(in)  :: iprt_nudob
 101   REAL, intent(in)     :: fdasta
 102   REAL, intent(in)     :: fdaend
 103   INTEGER, intent(in)  :: nudge_wind
 104   INTEGER, intent(in)  :: nudge_temp
 105   INTEGER, intent(in)  :: nudge_mois
 106   INTEGER, intent(in)  :: nudge_pstr
 107   REAL, intent(in) :: coef_wind
 108   REAL, intent(in) :: coef_temp
 109   REAL, intent(in) :: coef_mois
 110   REAL, intent(in) :: coef_pstr
 111   REAL, intent(inout)  :: rinxy
 112   REAL, intent(inout)  :: rinsig
 113   INTEGER, intent(in) :: npfi
 114   INTEGER, intent(in) :: ionf
 115   INTEGER, intent(in) :: nobs_prt      ! Number of current obs to print info.
 116   INTEGER, intent(in) :: idynin
 117   REAL, intent(inout) :: dtramp
 118   REAL, intent(in) :: xlatc_cg         ! center latitude  of coarse grid
 119   REAL, intent(in) :: xlonc_cg         ! center longitude of coarse grid
 120   REAL, intent(in) :: true_lat1
 121   REAL, intent(in) :: true_lat2
 122   INTEGER, intent(in) :: map_proj
 123   INTEGER, intent(in) :: i_parent_start(maxdom)
 124   INTEGER, intent(in) :: j_parent_start(maxdom)
 125   INTEGER, intent(in) :: parent_grid_ratio
 126   REAL, intent(in)     :: dt
 127   REAL, intent(in)     :: gmt
 128   INTEGER, intent(in)  :: julday
 129   INTEGER, intent(in)  :: max_obs         ! max number of observations
 130   INTEGER, intent(in)  :: nobs_ndg_vars
 131   INTEGER, intent(in)  :: nobs_err_flds
 132   INTEGER, intent(in)  :: nstat
 133   REAL, intent(inout)  :: varobs(nobs_ndg_vars, max_obs)
 134   REAL, intent(inout)  :: errf(nobs_err_flds, max_obs)
 135   REAL, intent(in)     :: dx           ! this-domain grid cell-size (m)
 136   INTEGER, INTENT(IN) :: kpbl( ims:ime, jms:jme )
 137   REAL, INTENT(IN) :: ht( ims:ime, jms:jme )
 138   REAL, INTENT(IN) :: mut( ims:ime , jms:jme )   ! Air mass on t-grid
 139   REAL, INTENT(IN) :: muu( ims:ime , jms:jme )   ! Air mass on u-grid
 140   REAL, INTENT(IN) :: muv( ims:ime , jms:jme )   ! Air mass on v-grid
 141   REAL, INTENT(IN) :: msftx( ims:ime , jms:jme )  ! Map scale on t-grid
 142   REAL, INTENT(IN) :: msfty( ims:ime , jms:jme )  ! Map scale on t-grid
 143   REAL, INTENT(IN) :: msfux( ims:ime , jms:jme )  ! Map scale on u-grid
 144   REAL, INTENT(IN) :: msfuy( ims:ime , jms:jme )  ! Map scale on u-grid
 145   REAL, INTENT(IN) :: msfvx( ims:ime , jms:jme )  ! Map scale on v-grid
 146   REAL, INTENT(IN) :: msfvy( ims:ime , jms:jme )  ! Map scale on v-grid
 147
 148   REAL, INTENT(IN) :: p_phy( ims:ime, kms:kme, jms:jme )
 149   REAL, INTENT(INOUT) :: t_tendf( ims:ime, kms:kme, jms:jme )
 150   REAL, INTENT(IN) :: t0
 151   REAL, INTENT(INOUT) :: savwt( nobs_ndg_vars, ims:ime, kms:kme, jms:jme )
 152
 153 #if ( EM_CORE == 1 )
 154   TYPE(fdob_type), intent(inout)  :: fdob
 155 #endif
 156
 157   REAL,   INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
 158   REAL,   INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
 159   REAL,   INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
 160   REAL,   INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
 161   REAL,   INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme ) ! Base press. (Pa)
 162   REAL,   INTENT(IN) :: ptop
 163   REAL,   INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme )  ! Press. pert. (Pa)
 164   REAL,   INTENT(IN) :: uratx( ims:ime, jms:jme )     ! On mass points
 165   REAL,   INTENT(IN) :: vratx( ims:ime, jms:jme )     !       "
 166   REAL,   INTENT(IN) :: tratx( ims:ime, jms:jme )     !       "
 167   REAL,   INTENT(INOUT) :: ru_tendf( ims:ime, kms:kme, jms:jme )
 168   REAL,   INTENT(INOUT) :: rv_tendf( ims:ime, kms:kme, jms:jme )
 169   REAL,   INTENT(INOUT) :: moist_tend( ims:ime, kms:kme, jms:jme )
 170
 171 ! Local variables
 172   logical            :: nudge_flag   ! Flag for doing nudging
 173   integer            :: KTAU         ! Forecast timestep
 174   real               :: xtime        ! Forecast time in minutes
 175   real               :: dtmin        ! dt in minutes
 176   integer            :: i, j, k      ! Loop counters.
 177   integer            :: idom         ! Loop counter.
 178   integer            :: nsta         ! Number of observation stations
 179   integer            :: infr         ! Frequency for obs input and error calc
 180   integer            :: idarst       ! Flag for calling sub errob on restart
 181   real               :: dtr          ! Abs value of dtramp (for dynamic init)
 182   real               :: tconst       ! Reciprocal of dtr
 183   integer :: KPBLJ(its:ite)          ! 1D temp array.
 184 #ifdef RAL
 185   real    :: HTIJ(ids:ide, jds:jde) = 0.  ! Terrain ht on global grid.
 186 #endif
 187
 188 #if ( EM_CORE == 1 )
 189   nudge_flag = (nudge_opt  .eq. 1)
 190
 191   if (.not. nudge_flag) return
 192
 193 !----------------------------------------------------------------------
 194 ! ***************       BEGIN FDDA SETUP SECTION        ***************
 195
 196 ! Calculate forecast time.
 197   dtmin = dt/60.
 198   xtime = dtmin*(itimestep-1)
 199
 200   ktau  = itimestep - 1        !ktau corresponds to xtime
 201
 202 ! DEFINE NSTA WHEN NOT NUDGING TO IND. OBS.
 203 ! print *,'in fddaobs_driver, xtime=',xtime
 204   IF(ktau.EQ.fdob%ktaur) THEN
 205      IF (iprt_nudob) PRINT *,3333,fdob%domain_tot
 206 !    print *,'ktau,ktaur,inest=',ktau,fdob%ktaur,inest
 207 3333 FORMAT(1X,'IN fddaobs_driver: I4DITOT = ',I2)
 208      nsta=0.
 209   ELSE
 210      nsta=fdob%nstat
 211   ENDIF
 212
 213   infr = ionf*(parent_grid_ratio**fdob%levidn(inest))
 214   nsta=fdob%nstat
 215   idarst = 0
 216   IF(restart .AND. ktau.EQ.fdob%ktaur) idarst=1
 217
 218 ! COMPUTE ERROR BETWEEN OBSERVATIONS and MODEL
 219   IF( nsta.GT.0 ) THEN
 220     IF( MOD(ktau,infr).EQ.0 .OR. idarst.EQ.1) THEN
 221
 222         CALL ERROB(inest, ub, vb, tb, t0, qvb, pbase, pp, rcp,       &
 223                    uratx, vratx, tratx,                              &
 224                    nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,    &
 225                    fdob%levidn, parid, fdob%nstat, fdob%nstaw,       &
 226                    nudge_wind, nudge_temp, nudge_mois, nudge_pstr,   &
 227                    fdob%rio, fdob%rjo, fdob%rko, varobs, errf,       &
 228                    i_parent_start, j_parent_start, ktau,             &
 229                    parent_grid_ratio, npfi, nobs_prt, iprt_errob,    &
 230                    ids,ide, jds,jde, kds,kde,                        &
 231                    ims,ime, jms,jme, kms,kme,                        &
 232                    its,ite, jts,jte, kts,kte)
 233     ENDIF
 234   ENDIF
 235
 236   fdob%tfaci=1.0
 237   IF(idynin.EQ.1.AND.nudge_opt.EQ.1) THEN
 238     dtr=ABS(dtramp)
 239     tconst=1./dtr
 240 ! FDAEND(IN) IS THE TIME IN MINUTES TO END THE DYNAMIC INITIALIZATION CY
 241     IF(xtime.LT.fdaend-dtr)THEN
 242       fdob%tfaci=1.
 243     ELSEIF(xtime.GE.fdaend-dtr.AND.xtime.LE.fdaend) THEN
 244       fdob%tfaci=(fdaend-xtime)*tconst
 245     ELSE
 246       fdob%tfaci=0.0
 247     ENDIF
 248     IF(ktau.EQ.fdob%ktaur.OR.MOD(ktau,10).EQ.0) THEN
 249       IF (iprt_nudob)                                                  &
 250          PRINT*,' DYNINOBS: IN,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST',   &
 251          ',TFACI: ',INEST,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST,         &
 252          fdob%TFACI
 253     ENDIF
 254   ENDIF
 255
 256 #ifdef RAL
 257 ! MEIXU: collect terrain array HT into a global array HTIJ
 258   CALL loc2glob(HT, HTIJ, "2D", "REAL",                  &
 259                 ids,ide, jds,jde, kds,kde,               &
 260                 ims,ime, jms,jme, kms,kme )
 261 ! MEIXU end
 262 #endif
 263 !
 264 ! ***************        END FDDA SETUP SECTION         ***************
 265 !----------------------------------------------------------------------
 266
 267 !----------------------------------------------------------------------
 268 ! ***************         BEGIN NUDGING SECTION         ***************
 269
 270   DO J = jts, jte
 271 !
 272 ! IF NUDGING SURFACE WINDS IN THE BOUNDARY LAYER, IF IWINDS(INEST+2)=1
 273 ! USE A SIMILARITY CORRECTION BASED ON ROUGHNESS TO APPLY 10M
 274 ! WIND TO THE SURFACE LAYER (K=KL) AT 40M.  TO DO THIS WE MUST
 275 ! STORE ROUGHNESS AND REGIME FOR EACH J SLICE AFTER THE CALL TO
 276 ! HIRPBL FOR LATER USE IN BLNUDGD.
 277 !
 278      DO I=its,ite
 279        KPBLJ(I)=KPBL(I,J)
 280      ENDDO
 281 !
 282 !--- OBS NUDGING FOR TEMP AND MOISTURE
 283 !
 284      NSTA=NSTAT
 285      IF(J .GT. 2 .and. J .LT.fdob%sn_end) THEN
 286        IF(nudge_temp.EQ.1 .AND. NSTA.GT.0)  &
 287        THEN
 288 !         write(6,*) 'calling nudob: IVAR=3, J = ',j
 289           CALL NUDOB(J, 3, t_tendf(ims,kms,j),                       &
 290                   inest, restart, ktau, fdob%ktaur, xtime,           &
 291                   mut(ims,j), msftx(ims,j), msfty(ims,j),            &
 292                   nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
 293                   npfi, ionf, rinxy, fdob%window,                    &
 294                   fdob%levidn,                                       &
 295                   parid, nstat, i_parent_start, j_parent_start,      &
 296                   fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,    &
 297                   parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
 298                   fdob%rko, fdob%timeob, varobs, errf,               &
 299                   pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
 300                   nudge_wind, nudge_temp, nudge_mois,                &
 301                   coef_wind, coef_temp, coef_mois,                   &
 302                   savwt(1,ims,kms,j), kpblj, 0,                      &
 303                   iprt_nudob,                                        &
 304                   ids,ide, jds,jde, kds,kde,                         & ! domain dims
 305                   ims,ime, jms,jme, kms,kme,                         & ! memory dims
 306                   its,ite, jts,jte, kts,kte         )                  ! tile   dims
 307 !         write(6,*) 'return from nudob: IVAR=3, J = ',j
 308        ENDIF
 309
 310        IF(nudge_mois.EQ.1 .AND. NSTA.GT.0)  &
 311        THEN
 312 !         write(6,*) 'calling nudob: IVAR=4, J = ',j
 313           CALL NUDOB(J, 4, moist_tend(ims,kms,j),                    &
 314                   inest, restart, ktau, fdob%ktaur, xtime,           &
 315                   mut(ims,j), msftx(ims,j), msfty(ims,j),            &
 316                   nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
 317                   npfi, ionf, rinxy, fdob%window,                    &
 318                   fdob%levidn,                                       &
 319                   parid, nstat, i_parent_start, j_parent_start,      &
 320                   fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,    &
 321                   parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
 322                   fdob%rko, fdob%timeob, varobs, errf,               &
 323                   pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
 324                   nudge_wind, nudge_temp, nudge_mois,                &
 325                   coef_wind, coef_temp, coef_mois,                   &
 326                   savwt(1,ims,kms,j), kpblj, 0,                      &
 327                   iprt_nudob,                                        &
 328                   ids,ide, jds,jde, kds,kde,                         & ! domain dims
 329                   ims,ime, jms,jme, kms,kme,                         & ! memory dims
 330                   its,ite, jts,jte, kts,kte         )                  ! tile   dims
 331 !         write(6,*) 'return from nudob: IVAR=4, J = ',j
 332        ENDIF
 333      ENDIF
 334
 335      IF(nudge_wind.EQ.1 .AND. NSTA.GT.0)    &
 336      THEN
 337 !       write(6,*) 'calling nudob: IVAR=1, J = ',j
 338         CALL NUDOB(J, 1, ru_tendf(ims,kms,j),                        &
 339                 inest, restart, ktau, fdob%ktaur, xtime,             &
 340                 muu(ims,j), msfux(ims,j), msfuy(ims,j),              &
 341                 nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
 342                 npfi, ionf, rinxy, fdob%window,                      &
 343                 fdob%levidn,                                         &
 344                 parid, nstat, i_parent_start, j_parent_start,        &
 345                 fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,      &
 346                 parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
 347                 fdob%rko, fdob%timeob, varobs, errf,                 &
 348                 pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
 349                 nudge_wind, nudge_temp, nudge_mois,                  &
 350                 coef_wind, coef_temp, coef_mois,                     &
 351                 savwt(1,ims,kms,j), kpblj, 0,                        &
 352                 iprt_nudob,                                          &
 353                 ids,ide, jds,jde, kds,kde,                           & ! domain dims
 354                 ims,ime, jms,jme, kms,kme,                           & ! memory dims
 355                 its,ite, jts,jte, kts,kte         )                    ! tile   dims
 356 !       write(6,*) 'return from nudob: IVAR=1, J = ',j
 357
 358 !       write(6,*) 'calling nudob: IVAR=2, J = ',j
 359         CALL NUDOB(J, 2, rv_tendf(ims,kms,j),                        &
 360                 inest, restart, ktau, fdob%ktaur, xtime,             &
 361                 muv(ims,j), msfvx(ims,j), msfvy(ims,j),              &
 362                 nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
 363                 npfi, ionf, rinxy, fdob%window,                      &
 364                 fdob%levidn,                                         &
 365                 parid, nstat, i_parent_start, j_parent_start,        &
 366                 fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,      &
 367                 parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
 368                 fdob%rko, fdob%timeob, varobs, errf,                 &
 369                 pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
 370                 nudge_wind, nudge_temp, nudge_mois,                  &
 371                 coef_wind, coef_temp, coef_mois,                     &
 372                 savwt(1,ims,kms,j), kpblj, 0,                        &
 373                 iprt_nudob,                                          &
 374                 ids,ide, jds,jde, kds,kde,                           & ! domain dims
 375                 ims,ime, jms,jme, kms,kme,                           & ! memory dims
 376                 its,ite, jts,jte, kts,kte         )                    ! tile   dims
 377 !       write(6,*) 'return from nudob: IVAR=2, J = ',j
 378      ENDIF
 379   ENDDO
 380 !
 381 ! --- END OF 4DDA
 382 !
 383   RETURN
 384 #endif
 385   END SUBROUTINE fddaobs_driver
 386
 387 END MODULE module_fddaobs_driver