source: trunk/WRF.COMMON/WRFV2/phys/module_fddaobs_driver.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 17.4 KB
Line 
1!WRF:MODEL_LAYER:PHYSICS
2MODULE 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 
35CONTAINS
36
37!-----------------------------------------------------------------------
38SUBROUTINE 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, twindo,                 &
45               npfi, ionf, 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               msft, msfu, msfv, 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  REAL, intent(inout)  :: twindo
114  INTEGER, intent(in) :: npfi
115  INTEGER, intent(in) :: ionf
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) :: msft( ims:ime , jms:jme )  ! Map scale on t-grid
142  REAL, INTENT(IN) :: msfu( ims:ime , jms:jme )  ! Map scale on u-grid
143  REAL, INTENT(IN) :: msfv( ims:ime , jms:jme )  ! Map scale on v-grid
144
145  REAL, INTENT(IN) :: p_phy( ims:ime, kms:kme, jms:jme )
146  REAL, INTENT(INOUT) :: t_tendf( ims:ime, kms:kme, jms:jme )
147  REAL, INTENT(IN) :: t0
148  REAL, INTENT(INOUT) :: savwt( nobs_ndg_vars, ims:ime, kms:kme, jms:jme )
149
150#if ( EM_CORE == 1 )
151  TYPE(fdob_type), intent(inout)  :: fdob
152#endif
153
154  REAL,   INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
155  REAL,   INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
156  REAL,   INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
157  REAL,   INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
158  REAL,   INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme ) ! Base press. (Pa)
159  REAL,   INTENT(IN) :: ptop
160  REAL,   INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme )  ! Press. pert. (Pa)
161  REAL,   INTENT(IN) :: uratx( ims:ime, jms:jme )     ! On mass points
162  REAL,   INTENT(IN) :: vratx( ims:ime, jms:jme )     !       "
163  REAL,   INTENT(IN) :: tratx( ims:ime, jms:jme )     !       "
164  REAL,   INTENT(INOUT) :: ru_tendf( ims:ime, kms:kme, jms:jme )
165  REAL,   INTENT(INOUT) :: rv_tendf( ims:ime, kms:kme, jms:jme )
166  REAL,   INTENT(INOUT) :: moist_tend( ims:ime, kms:kme, jms:jme )
167
168! Local variables
169  logical            :: nudge_flag   ! Flag for doing nudging
170  integer            :: KTAU         ! Forecast timestep
171  real               :: xtime        ! Forecast time in minutes
172  real               :: dtmin        ! dt in minutes
173  integer            :: i, j, k      ! Loop counters.
174  integer            :: idom         ! Loop counter.
175  integer            :: nsta         ! Number of observation stations
176  integer            :: infr         ! Frequency for obs input and error calc
177  integer            :: idarst       ! Flag for calling sub errob on restart
178  real               :: dtr          ! Abs value of dtramp (for dynamic init)
179  real               :: tconst       ! Reciprocal of dtr
180  integer :: KPBLJ(its:ite)          ! 1D temp array.
181#ifdef RAL
182  real    :: HTIJ(ids:ide, jds:jde) = 0.  ! Terrain ht on global grid.
183#endif
184
185#if ( EM_CORE == 1 )
186  nudge_flag = (nudge_opt  .eq. 1)
187
188  if (.not. nudge_flag) return
189
190!----------------------------------------------------------------------
191! ***************       BEGIN FDDA SETUP SECTION        ***************
192
193! Calculate forecast time.
194  dtmin = dt/60.     
195  xtime = dtmin*(itimestep-1)
196
197  ktau  = itimestep - 1        !ktau corresponds to xtime
198
199! DEFINE NSTA WHEN NOT NUDGING TO IND. OBS.
200! print *,'in fddaobs_driver, xtime=',xtime
201  IF(ktau.EQ.fdob%ktaur) THEN
202     IF (iprt_nudob) PRINT *,3333,fdob%domain_tot
203!    print *,'ktau,ktaur,inest=',ktau,fdob%ktaur,inest
2043333 FORMAT(1X,'IN fddaobs_driver: I4DITOT = ',I2)
205     nsta=0.
206  ELSE
207     nsta=fdob%nstat
208  ENDIF
209 
210  infr = ionf*(parent_grid_ratio**fdob%levidn(inest))
211  nsta=fdob%nstat
212  idarst = 0
213  IF(restart .AND. ktau.EQ.fdob%ktaur) idarst=1
214
215! COMPUTE ERROR BETWEEN OBSERVATIONS and MODEL
216  IF( nsta.GT.0 ) THEN
217    IF( MOD(ktau,infr).EQ.0 .OR. idarst.EQ.1) THEN
218
219        CALL ERROB(inest, ub, vb, tb, t0, qvb, pbase, pp, rcp,       &
220                   uratx, vratx, tratx,                              &
221                   nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,    &
222                   fdob%levidn, parid, fdob%nstat,                   &
223                   nudge_wind, nudge_temp, nudge_mois, nudge_pstr,   &
224                   fdob%rio, fdob%rjo, fdob%rko, varobs, errf,       &
225                   i_parent_start, j_parent_start, ktau,             &
226                   parent_grid_ratio, npfi, iprt_errob,              &
227                   ids,ide, jds,jde, kds,kde,                        &
228                   ims,ime, jms,jme, kms,kme,                        &
229                   its,ite, jts,jte, kts,kte)
230    ENDIF
231  ENDIF
232
233  fdob%tfaci=1.0
234  IF(idynin.EQ.1.AND.nudge_opt.EQ.1) THEN
235    dtr=ABS(dtramp)
236    tconst=1./dtr
237! FDAEND(IN) IS THE TIME IN MINUTES TO END THE DYNAMIC INITIALIZATION CY
238    IF(xtime.LT.fdaend-dtr)THEN
239      fdob%tfaci=1.
240    ELSEIF(xtime.GE.fdaend-dtr.AND.xtime.LE.fdaend) THEN
241      fdob%tfaci=(fdaend-xtime)*tconst
242    ELSE
243      fdob%tfaci=0.0
244    ENDIF
245    IF(ktau.EQ.fdob%ktaur.OR.MOD(ktau,10).EQ.0) THEN
246      IF (iprt_nudob)                                                  &
247         PRINT*,' DYNINOBS: IN,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST',   &
248         ',TFACI: ',INEST,KTAU,XTIME,FDAEND,DTRAMP,DTR,TCONST,         &
249         fdob%TFACI
250    ENDIF
251  ENDIF
252
253#ifdef RAL
254! MEIXU: collect terrain array HT into a global array HTIJ
255  CALL loc2glob(HT, HTIJ, "2D", "REAL",                  &
256                ids,ide, jds,jde, kds,kde,               &
257                ims,ime, jms,jme, kms,kme )
258! MEIXU end
259#endif
260!
261! ***************        END FDDA SETUP SECTION         ***************
262!----------------------------------------------------------------------
263
264!----------------------------------------------------------------------
265! ***************         BEGIN NUDGING SECTION         ***************
266
267  DO J = jts, jte
268!
269! IF NUDGING SURFACE WINDS IN THE BOUNDARY LAYER, IF IWINDS(INEST+2)=1
270! USE A SIMILARITY CORRECTION BASED ON ROUGHNESS TO APPLY 10M
271! WIND TO THE SURFACE LAYER (K=KL) AT 40M.  TO DO THIS WE MUST
272! STORE ROUGHNESS AND REGIME FOR EACH J SLICE AFTER THE CALL TO
273! HIRPBL FOR LATER USE IN BLNUDGD.
274!
275     DO I=its,ite
276       KPBLJ(I)=KPBL(I,J)
277     ENDDO
278!
279!--- OBS NUDGING FOR TEMP AND MOISTURE
280!
281     NSTA=NSTAT
282     IF(J .GT. 2 .and. J .LT.fdob%sn_end) THEN
283       IF(nudge_temp.EQ.1 .AND. NSTA.GT.0)  &
284       THEN
285!         write(6,*) 'calling nudob: IVAR=3, J = ',j
286          CALL NUDOB(J, 3, t_tendf(ims,kms,j),                       &
287                  inest, restart, ktau, fdob%ktaur, xtime,           &
288                  mut(ims,j), msft(ims,j),                           &
289                  nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
290                  npfi, ionf, rinxy, twindo,                         &
291                  fdob%levidn,                                       &
292                  parid, nstat, i_parent_start, j_parent_start,      &
293                  fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,    &
294                  parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
295                  fdob%rko, fdob%timeob, varobs, errf,               &
296                  pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
297                  nudge_wind, nudge_temp, nudge_mois,                &
298                  coef_wind, coef_temp, coef_mois,                   &
299                  savwt(1,ims,kms,j), kpblj, 0,                      &
300                  iprt_nudob,                                        &
301                  ids,ide, jds,jde, kds,kde,                         & ! domain dims
302                  ims,ime, jms,jme, kms,kme,                         & ! memory dims
303                  its,ite, jts,jte, kts,kte         )                  ! tile   dims
304!         write(6,*) 'return from nudob: IVAR=3, J = ',j
305       ENDIF
306
307       IF(nudge_mois.EQ.1 .AND. NSTA.GT.0)  &
308       THEN
309!         write(6,*) 'calling nudob: IVAR=4, J = ',j
310          CALL NUDOB(J, 4, moist_tend(ims,kms,j),                    &
311                  inest, restart, ktau, fdob%ktaur, xtime,           &
312                  mut(ims,j), msft(ims,j),                           &
313                  nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,     &
314                  npfi, ionf, rinxy, twindo,                         &
315                  fdob%levidn,                                       &
316                  parid, nstat, i_parent_start, j_parent_start,      &
317                  fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,    &
318                  parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,  &
319                  fdob%rko, fdob%timeob, varobs, errf,               &
320                  pbase(ims,kms,j), ptop, pp(ims,kms,j),             &
321                  nudge_wind, nudge_temp, nudge_mois,                &
322                  coef_wind, coef_temp, coef_mois,                   &
323                  savwt(1,ims,kms,j), kpblj, 0,                      &
324                  iprt_nudob,                                        &
325                  ids,ide, jds,jde, kds,kde,                         & ! domain dims
326                  ims,ime, jms,jme, kms,kme,                         & ! memory dims
327                  its,ite, jts,jte, kts,kte         )                  ! tile   dims
328!         write(6,*) 'return from nudob: IVAR=4, J = ',j
329       ENDIF
330     ENDIF
331
332     IF(nudge_wind.EQ.1 .AND. NSTA.GT.0)    &
333     THEN
334!       write(6,*) 'calling nudob: IVAR=1, J = ',j
335        CALL NUDOB(J, 1, ru_tendf(ims,kms,j),                        &
336                inest, restart, ktau, fdob%ktaur, xtime,             &
337                muu(ims,j), msfu(ims,j),                             &
338                nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
339                npfi, ionf, rinxy, twindo,                           &
340                fdob%levidn,                                         &
341                parid, nstat, i_parent_start, j_parent_start,        &
342                fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,      &
343                parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
344                fdob%rko, fdob%timeob, varobs, errf,                 &
345                pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
346                nudge_wind, nudge_temp, nudge_mois,                  &
347                coef_wind, coef_temp, coef_mois,                     &
348                savwt(1,ims,kms,j), kpblj, 0,                        &
349                iprt_nudob,                                          &
350                ids,ide, jds,jde, kds,kde,                           & ! domain dims
351                ims,ime, jms,jme, kms,kme,                           & ! memory dims
352                its,ite, jts,jte, kts,kte         )                    ! tile   dims
353!       write(6,*) 'return from nudob: IVAR=1, J = ',j
354
355!       write(6,*) 'calling nudob: IVAR=2, J = ',j
356        CALL NUDOB(J, 2, rv_tendf(ims,kms,j),                        &
357                inest, restart, ktau, fdob%ktaur, xtime,             &
358                muv(ims,j), msfv(ims,j),                             &
359                nobs_ndg_vars, nobs_err_flds, max_obs, maxdom,       &
360                npfi, ionf, rinxy, twindo,                           &
361                fdob%levidn,                                         &
362                parid, nstat, i_parent_start, j_parent_start,        &
363                fdob, fdob%lev_in_ob, fdob%plfo, fdob%nlevs_ob,      &
364                parent_grid_ratio, dx, dtmin, fdob%rio, fdob%rjo,    &
365                fdob%rko, fdob%timeob, varobs, errf,                 &
366                pbase(ims,kms,j), ptop, pp(ims,kms,j),               &
367                nudge_wind, nudge_temp, nudge_mois,                  &
368                coef_wind, coef_temp, coef_mois,                     &
369                savwt(1,ims,kms,j), kpblj, 0,                        &
370                iprt_nudob,                                          &
371                ids,ide, jds,jde, kds,kde,                           & ! domain dims
372                ims,ime, jms,jme, kms,kme,                           & ! memory dims
373                its,ite, jts,jte, kts,kte         )                    ! tile   dims
374!       write(6,*) 'return from nudob: IVAR=2, J = ',j
375     ENDIF
376  ENDDO
377!
378! --- END OF 4DDA
379!
380  RETURN
381#endif
382  END SUBROUTINE fddaobs_driver
383
384END MODULE module_fddaobs_driver
Note: See TracBrowser for help on using the repository browser.