source: branches/jerome/src_flexwrf_v3.1/readwind_nests_timeav.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 15.0 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23
24      subroutine readwind_nests_timeav(indj,n,uuhn,vvhn,wwhn)
25!**********************************************************************
26!                                                                     *
27!             TRAJECTORY MODEL SUBROUTINE READWIND                    *
28!                                                                     *
29!**********************************************************************
30!                                                                     *
31!  April 2012, J. Brioude: This routine handles the difference in time
32! for time-average fields.
33!**********************************************************************
34!                                                                     *
35! Note:  This is the FLEXPART_WRF version of subroutine readwind.     *
36!    The met fields are read from WRF netcdf output files.            *
37!    There are many differences from the FLEXPART version.            *
38!                                                                     *
39! DESCRIPTION:                                                        *
40!                                                                     *
41! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
42! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
43!                                                                     *
44! INPUT:                                                              *
45! indj               indicates number of the wind field to be read in *
46! n                  temporal index for meteorological fields (1 to 3)*
47!                                                                     *
48! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
49!                                                                     *
50! wfname             File name of data to be read in                  *
51! nx,ny,nuvz,nwz     expected field dimensions                        *
52! nlev_ec            number of "T-grid" vertical levels wwf model     *
53!                    (the unstaggered "bottom_top" dimension)         *
54! uu,vv,ww           wind fields                                      *
55! tt,qv              temperature and specific humidity                *
56! ps                 surface pressure                                 *
57!                                                                     *
58!**********************************************************************
59!
60
61!      include 'includepar'
62!      include 'includecom'
63  use par_mod
64  use com_mod
65
66! subr arguments
67      integer :: indj, n
68
69      real(kind=4) :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
70      real(kind=4) :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
71      real(kind=4) :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
72      real(kind=4) :: uuhn2(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
73      real(kind=4) :: vvhn2(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
74      real(kind=4) :: wwhn2(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
75
76! local variables
77      integer,parameter :: ndims_max=4
78
79      integer :: i, idiagaa, ierr, ifn, itime
80      integer :: iduma,indj2,l
81      integer :: j, jhhmmss, jyyyymmdd
82      integer :: k, kbgn
83      integer :: lendim(ndims_max), lendim_exp(ndims_max), &
84          lendim_max(ndims_max)
85      integer :: levdiff2,deltat,deltat2
86      integer :: ndims, ndims_exp
87      integer :: n_west_east, n_south_north, n_bottom_top
88      integer :: m_grid_id_dum, m_parent_grid_id_dum, &
89        m_parent_grid_ratio_dum,  &
90        i_parent_start_dum, j_parent_start_dum, &
91        map_proj_id_dum,  &
92        ext_scalar,pbl_physics,mp_physics_dum
93
94      real :: dx_met, dy_met
95      real :: duma, dumb, dumc, dumd, dume
96      real :: dumdz
97!      real(kind=4) :: dumarray_aa(nwzmax+1)
98!      real(kind=4) :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
99      real :: dumarray_aa(nwzmax+1)
100      real :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
101      real :: ewater_mb, esatwater_mb
102      real :: ew      ! this is an external function
103      real :: map_stdlon_dum, map_truelat1_dum, map_truelat2_dum
104      real :: pint
105      real :: toler
106
107!      real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
108      real :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
109      real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
110
111      real(kind=dp) :: jul,juldate
112
113      character(len=160) :: fnamenc, varname,fnamenc2
114
115      logical :: hflswitch
116
117!
118!   get grid info from the wrf netcdf file
119!   and check it for consistency against values from gridcheck
120!
121
122!
123! main loop -- process each nest
124!
125      do l=1,numbnests
126
127      m = numpath+2*(l-1)+1
128      fnamenc = path(m)(1:length(m)) // wfnamen(l,indj)
129
130!      fnamenc = path(2)(1:length(2))//wfname(indj)
131      idiagaa = 0
132
133      call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, &
134        n_west_east, n_south_north, n_bottom_top, &
135        dx_met, dy_met,  &
136        m_grid_id_dum, m_parent_grid_id_dum, m_parent_grid_ratio_dum, &
137        i_parent_start_dum, j_parent_start_dum, &
138        map_proj_id_dum, map_stdlon_dum,  &
139        map_truelat1_dum, map_truelat2_dum, &
140        ext_scalar,pbl_physics,mp_physics_dum )
141      if (ierr .ne. 0) then
142          write(*,9100) 'error getting gridinfor for met file', fnamenc
143          stop
144      end if
145
146      i_parent_start_dum = i_parent_start_dum-1
147      j_parent_start_dum = j_parent_start_dum-1
148
149
1509100  format( / '*** readwind -- ', a )
1519110  format( / '*** readwind -- ', a, 1x, i8 / &
152        'file = ', a )
1539120  format( / '*** readwind -- ', a, 2(1x,i8) / &
154        'file = ', a )
1559130  format( / '*** readwind -- ', a, 3(1x,i8) / &
156        'file = ', a )
1579115  format( / '*** readwind -- ', a / a, 1x, i8 / &
158        'file = ', a )
1599125  format( / '*** readwind -- ', a / a, 2(1x,i8) / &
160        'file = ', a )
1619135  format( / '*** readwind -- ', a / a, 3(1x,i8) / &
162        'file = ', a )
163
164
165! READ THE FIRST FILE
166       if (ldirect.eq.1) indj2=indj
167       if (ldirect.eq.-1) indj2=indj
168       deltat=wfdt(indj2)
169      fnamenc = path(m)(1:length(m)) // wfnamen(l,indj2)
170! locate the date/time in the file
171      itime = 0
1721100  itime = itime + 1
173      call read_ncwrfout_1datetime( ierr, fnamenc, &
174          itime, jyyyymmdd, jhhmmss )
175      if (ierr .eq. -1) then
176          write(*,9100) 'error reading time from met file', fnamenc
177          stop
178      else if (ierr .ne. 0) then
179          write(*,9125) 'unable to locate date/time in met file',  &
180              'indj, itime =', indj2, itime, fnamenc
181          stop
182      else
183          jul = juldate( jyyyymmdd, jhhmmss )
184          duma = (jul-bdate)*86400.
185          iduma = nint(duma)
186          if (iduma .ne. wftime(indj2)) goto 1100
187      end if
188      write(*,*)
189      write(*,*) 'readwind_nests processing wrfout file ='
190      write(*,*) fnamenc
191      write(*,*) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss
192
193      kbgn = 1 + add_sfc_level
194
195      if (wind_option.eq.1) varname = 'AVGFLX_RUM'
196      do i = 1, ndims_max
197          lendim_exp(i) = 0
198          lendim_max(i) = 1
199      end do
200      lendim_exp(1) = nxn(l)+1
201      lendim_max(1) = nxmaxn
202      lendim_exp(2) = nyn(l)
203      lendim_max(2) = nymaxn
204      lendim_exp(3) = nuvz-add_sfc_level
205      lendim_max(3) = nuvzmax
206      ndims_exp = 4
207      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
208          varname, uuhn(0,0,kbgn,l), &
209          itime, &
210          ndims, ndims_exp, ndims_max, &
211          lendim, lendim_exp, lendim_max )
212      if (ierr .ne. 0) then
213          write(*,9100) 'error doing ncread of U', fnamenc
214      if (wind_option.le.0) print*,'you asked snapshot winds'
215      if (wind_option.eq.1) print*,'you asked mean winds'
216        print*,'change wind_option'
217
218          stop
219      end if
220
221
222! v wind velocity
223!   the wrf output file contains (nuvz-add_sfc_level) levels
224!   read the data into k=kbgn,nuvz
225!   (interpolate it from "V-grid" to "T-grid" later)
226      if (wind_option.eq.1) varname = 'AVGFLX_RVM'
227      lendim_exp(1) = nxn(l)
228      lendim_max(1) = nxmaxn
229      lendim_exp(2) = nyn(l)+1
230      lendim_max(2) = nymaxn
231      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
232          varname, vvhn(0,0,kbgn,l), &
233          itime, &
234          ndims, ndims_exp, ndims_max, &
235          lendim, lendim_exp, lendim_max )
236      if (ierr .ne. 0) then
237          write(*,9100) 'error doing ncread of V', fnamenc
238      if (wind_option.eq.0) print*,'you asked snapshot winds'
239      if (wind_option.eq.1) print*,'you asked mean winds'
240        print*,'change wind_option'
241          stop
242      end if
243
244
245! w wind velocity
246!   this is on the "W-grid", and
247!   the wrf output file contains nwz levels, so no shifting needed
248      if (wind_option.eq.1) varname = 'AVGFLX_WWM'
249!      print*,'varname',varname
250      lendim_exp(1) = nxn(l)
251      lendim_max(1) = nxmaxn
252      lendim_exp(2) = nyn(l)
253      lendim_max(2) = nymaxn
254      lendim_exp(3) = nwz
255      lendim_max(3) = nwzmax
256      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
257          varname, wwhn(0,0,1,l), &
258          itime, &
259          ndims, ndims_exp, ndims_max, &
260          lendim, lendim_exp, lendim_max )
261      if (ierr .ne. 0) then
262          write(*,9100) 'error doing ncread of W', fnamenc
263      if (wind_option.eq.0) print*,'you asked snapshot winds'
264      if (wind_option.eq.1) print*,'you asked mean winds'
265        print*,'change wind_option'
266          stop
267      end if
268
269! READ THE SECOND FILE and average
270! deltat must be equal to deltat2, otherwise the time-average wind cannot be
271! fixed.
272! deltat is assumed to be the same than the WRF output.
273
274       if (ldirect.eq.1) indj2=indj+1
275       if (ldirect.eq.-1) indj2=indj-1
276       deltat2=wfdt(indj2)
277      fnamenc = path(m)(1:length(m)) // wfnamen(l,indj2)
278
279     if (deltat.eq.deltat2 .and. indj2.ge.1 .and.indj2.le.numbwf) then
280! locate the date/time in the file
281      itime = 0
2821101  itime = itime + 1
283      call read_ncwrfout_1datetime( ierr, fnamenc, &
284          itime, jyyyymmdd, jhhmmss )
285      if (ierr .eq. -1) then
286          write(*,9100) 'error reading time from met file', fnamenc
287          stop
288      else if (ierr .ne. 0) then
289          write(*,9125) 'unable to locate date/time in met file',  &
290              'indj, itime =', indj2, itime, fnamenc
291          stop
292      else
293          jul = juldate( jyyyymmdd, jhhmmss )
294          duma = (jul-bdate)*86400.
295          iduma = nint(duma)
296          if (iduma .ne. wftime(indj2)) goto 1101
297      end if
298      if (option_verbose.eq.1) then
299      write(*,*)
300      write(*,*) 'readwind processing wrfout file ='
301      write(*,*) fnamenc
302      write(*,*) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss
303      endif
304      kbgn = 1 + add_sfc_level
305
306      if (wind_option.eq.1) varname = 'AVGFLX_RUM'
307      do i = 1, ndims_max
308          lendim_exp(i) = 0
309          lendim_max(i) = 1
310      end do
311      lendim_exp(1) = nxn(l)+1
312      lendim_max(1) = nxmaxn
313      lendim_exp(2) = nyn(l)
314      lendim_max(2) = nymaxn
315      lendim_exp(3) = nuvz-add_sfc_level
316      lendim_max(3) = nuvzmax
317      ndims_exp = 4
318      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
319          varname, uuhn2(0,0,kbgn,l), &
320          itime, &
321          ndims, ndims_exp, ndims_max, &
322          lendim, lendim_exp, lendim_max )
323      if (ierr .ne. 0) then
324          write(*,9100) 'error doing ncread of U', fnamenc
325      if (wind_option.le.0) print*,'you asked snapshot winds'
326      if (wind_option.eq.1) print*,'you asked mean winds'
327        print*,'change wind_option'
328
329          stop
330      end if
331
332
333! v wind velocity
334!   the wrf output file contains (nuvz-add_sfc_level) levels
335!   read the data into k=kbgn,nuvz
336!   (interpolate it from "V-grid" to "T-grid" later)
337      if (wind_option.eq.1) varname = 'AVGFLX_RVM'
338      lendim_exp(1) = nxn(l)
339      lendim_max(1) = nxmaxn
340      lendim_exp(2) = nyn(l)+1
341      lendim_max(2) = nymaxn
342      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
343          varname, vvhn2(0,0,kbgn,l), &
344          itime, &
345          ndims, ndims_exp, ndims_max, &
346          lendim, lendim_exp, lendim_max )
347      if (ierr .ne. 0) then
348          write(*,9100) 'error doing ncread of V', fnamenc
349      if (wind_option.eq.0) print*,'you asked snapshot winds'
350      if (wind_option.eq.1) print*,'you asked mean winds'
351        print*,'change wind_option'
352          stop
353      end if
354
355
356! w wind velocity
357!   this is on the "W-grid", and
358!   the wrf output file contains nwz levels, so no shifting needed
359      if (wind_option.eq.1) varname = 'AVGFLX_WWM'
360!     print*,'varname',varname
361      lendim_exp(1) = nxn(l)
362      lendim_max(1) = nxmaxn
363      lendim_exp(2) = nyn(l)
364      lendim_max(2) = nymaxn
365      lendim_exp(3) = nwz
366      lendim_max(3) = nwzmax
367      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
368          varname, wwhn2(0,0,1,l), &
369          itime, &
370          ndims, ndims_exp, ndims_max, &
371          lendim, lendim_exp, lendim_max )
372      if (ierr .ne. 0) then
373          write(*,9100) 'error doing ncread of W', fnamenc
374      if (wind_option.eq.0) print*,'you asked snapshot winds'
375      if (wind_option.eq.1) print*,'you asked mean winds'
376        print*,'change wind_option'
377          stop
378      end if
379
380
381         do k=1,nuvzmax
382        do j=0,nymax-1
383       do i=0,nxmax-1
384           uuhn(i,j,k,l)=0.5*(uuhn(i,j,k,l)+uuhn2(i,j,k,l))
385           vvhn(i,j,k,l)=0.5*(vvhn(i,j,k,l)+vvhn2(i,j,k,l))
386         enddo
387        enddo
388       enddo
389         do k=1,nwzmax
390        do j=0,nymax-1
391       do i=0,nxmax-1
392           wwhn(i,j,k,l)=0.5*(wwhn(i,j,k,l)+wwhn2(i,j,k,l))
393         enddo
394        enddo
395       enddo
396
397
398      endif ! test on deltat
399
400      enddo ! loop over the nests
401
402      return   
403      end subroutine readwind_nests_timeav
404
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG