Ticket #14: read_ncwrfout.f90

File read_ncwrfout.f90, 31.4 KB (added by jbrioude, 11 years ago)

read_ncwrf version 3

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!**********************************************************************
25! FLEXPART SOURCE FILE READ_NCWRFOUT - CONTAINS                       *
26!                                                                     *
27!    SUBROUTINE READ_NCWRFOUT_GRIDINFO                                *
28!    SUBROUTINE READ_NCWRFOUT_1REALFIELD                              *
29!    SUBROUTINE READ_NCWRFOUT_1DATETIME                               *
30!                                                                     *
31!**********************************************************************
32!                                                                     *
33!             AUTHOR:      R. Easter & J. Fast, PNNL                  *
34!             DATE:        2005-autumn-??                             *
35!             LAST UPDATE: same                                       *
36!                                                                     *
37!**********************************************************************
38!                                                                     *
39! DESCRIPTION:                                                        *
40!                                                                     *
41! These routines read the netcdf wrf output files.                    *
42!                                                                     *
43!       13 JUNE 2007, add more arguments-- ext_scalar,pbl_physcis
44!**********************************************************************
45
46
47!-----------------------------------------------------------------------
48        subroutine read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, &
49          n_west_east, n_south_north, n_bottom_top,  &
50          dx_met, dy_met,  &
51          m_grid_id, m_parent_grid_id, m_parent_grid_ratio, &
52          i_parent_start, j_parent_start, &
53          map_proj_id, map_stdlon, map_truelat1, map_truelat2, &
54          ext_scalar,pbl_physics,microphysics)
55!
56!   reads grid definition information from a netcdf wrf output file
57!
58!   arguments
59!       ierr - output - if non-zero, an error occurred
60!               while opening or reading from the file
61!       idiagaa - input - if positive, testing diagnostics are printed
62!       fnamenc - input - path+filename of the wrf output file
63!       n_west_east - output - east_west dimension of the "T-grid"
64!       n_south_north - output - south_north dimension of the "T-grid"
65!       n_bottom_top - output - bottom_top dimension of the "T-grid"
66!       dx_met, dy_met - output - horizontal grid spacing (m)
67!       m_grid_id - output - grid id number
68!       m_parent_grid_id - output - grid id number of parent grid
69!       m_parent_grid_ratio - output - ratio of parent grid dxy to current grid dxy
70!       i_parent_start, j_parent_start - output - location of lower left corner
71!               of current grid relative to the parent grid.
72!       (if there is no parent grid, then the above 4 "...parent..." variables
73!               area set to -987.)
74!       map_proj_id - WRF map projection id (2=polar stereographic)
75!       map_stdlon - map projection standard longitude (deg)
76!       map_truelat1, truelat2 - map projection true latitudes (deg)
77!
78!       ext_scalar  -- dimension of ex_scalar
79!       pbl_physics  -- type of PBL scheme
80!       microphysics   -- micorphysice scheme used
81
82        include 'netcdf.inc'
83!        use netcdf
84!       implicit none
85
86!   arguments
87        integer :: ierr, idiagaa, &
88        n_west_east, n_south_north, n_bottom_top,  &
89        m_grid_id, m_parent_grid_id, m_parent_grid_ratio, &
90        i_parent_start, j_parent_start, map_proj_id, &
91        ext_scalar,pbl_physics,microphysics
92
93        real :: dx_met, dy_met, map_stdlon, map_truelat1, map_truelat2
94
95!character*(*) fnamenc
96        character(len=160) :: fnamenc
97
98!   local variables
99        integer,parameter :: maxdim=20
100        integer,parameter :: ibadaa= -987
101        integer,parameter ::  xbadaa= -987
102
103        integer :: i, iatt, idimid_unlim, idum, iret, ivtype
104        integer :: l, lenatt, lendim(maxdim)
105        integer :: natts_tot, ncid, ndims_tot, nvars_tot
106        integer :: n_west_east_stag, n_south_north_stag, n_bottom_top_stag
107
108        real :: duma
109        real, allocatable, dimension(:) :: duma_alloc
110
111        character(len=80) :: dimname(maxdim)
112        character(len=80) :: attname
113        character(len=1000) :: dumch1000
114
115!   externals
116!       integer nf_close
117!       integer nf_inq
118!       integer nf_inq_dim
119!       integer nf_open
120
121
122
123!   initialize with "missing values"
124        n_west_east = ibadaa
125        n_south_north = ibadaa
126        n_bottom_top = ibadaa
127        dx_met = xbadaa
128        dy_met = xbadaa
129        m_grid_id = ibadaa
130        m_parent_grid_id = ibadaa
131        m_parent_grid_ratio = ibadaa
132        i_parent_start = ibadaa
133        j_parent_start = ibadaa
134        ext_scalar = ibadaa
135        pbl_physics = ibadaa
136        microphysics = ibadaa
137!
138!   open the netcdf file
139!
140        ncid = 10
141!C          write(*,*)'xxx inside read_ncwrfout.f fnamenc=',fnamenc
142!C          write(*,*)'mp_physics=',microphysics
143
144!       print*,'filename ',fnamenc
145        iret = nf_open( fnamenc, NF_NOWRITE, ncid )
146        if (iret .ne. nf_noerr) then
147            write(*,9100) 'error doing open 123',fnamenc
148!         print*, NF_NOWRITE, ncid,iret
149            ierr = -1
150            return
151        end if
152
1539100  format( / '*** read_ncwrfout_gridinfo -- ', a / &
154        'file = ', a )
1559110  format( / '*** read_ncwrfout_gridinfo -- ', a, 1x, i8 / &
156        'file = ', a )
1579120  format( / '*** read_ncwrfout_gridinfo -- ', a, 2(1x,i8) / &
158        'file = ', a )
159
16090030  format( a, 2i6, 2(2x,a) )
161 
162!
163! get information on dimensions
164!
165        iret = nf_inq( ncid,  &
166                ndims_tot, nvars_tot, natts_tot, idimid_unlim )
167        if (iret .ne. nf_noerr) then
168            write(*,9100) 'error inquiring dimensions', fnamenc
169            ierr = -2
170            return
171        end if
172
173        n_west_east_stag = ibadaa
174        n_south_north_stag = ibadaa
175        n_bottom_top_stag = ibadaa
176
177        do i = 1, min(ndims_tot,maxdim)
178            iret = nf_inq_dim( ncid, i, dimname(i), lendim(i) )
179            if (iret .ne. nf_noerr) then
180                write(*,9110) 'error inquiring dimensions for dim#', i, fnamenc
181                ierr = -2
182                return
183            end if
184        end do
185
186        do i = 1, min(ndims_tot,maxdim)
187            if (dimname(i) .eq. 'west_east')    &
188                               n_west_east = lendim(i)
189            if (dimname(i) .eq. 'south_north') &
190                               n_south_north = lendim(i)
191            if (dimname(i) .eq. 'bottom_top')  &
192                               n_bottom_top = lendim(i)
193            if (dimname(i) .eq. 'west_east_stag')   &
194                               n_west_east_stag = lendim(i)
195            if (dimname(i) .eq. 'south_north_stag') &
196                               n_south_north_stag = lendim(i)
197            if (dimname(i) .eq. 'bottom_top_stag')  &
198                               n_bottom_top_stag = lendim(i)
199            if (dimname(i) .eq. 'ext_scalar') &
200                               ext_scalar = lendim(i)
201
202        end do
203
204        if (idiagaa .gt. 0) then
205        write(*,9100) 'diagnostics', fnamenc
206        do i = 1, min(ndims_tot,maxdim)
207            write(*,90030) 'dim #, len, name =',  &
208                i, lendim(i), dimname(i)
209        end do
210        end if
211
212        if ((n_west_east .le. 0) .or. &
213            (n_west_east+1 .ne. n_west_east_stag)) then
214            write(*,9120) 'bad n_west_east, n_west_east_stag = ',  &
215                n_west_east, n_west_east_stag, fnamenc
216            ierr = -3
217            return
218        end if
219
220        if ((n_south_north .le. 0) .or. &
221            (n_south_north+1 .ne. n_south_north_stag)) then
222            write(*,9120) 'bad n_south_north, n_south_north_stag = ',  &
223                n_south_north, n_south_north_stag, fnamenc
224            ierr = -3
225            return
226        end if
227
228        if ((n_bottom_top .le. 0) .or. &
229            (n_bottom_top+1 .ne. n_bottom_top_stag)) then
230            write(*,9120) 'bad n_bottom_top, n_bottom_top_stag = ',  &
231                n_bottom_top, n_bottom_top_stag, fnamenc
232            ierr = -3
233            return
234        end if
235
236!
237!   get information on global attributes
238!
239
240!   first just do diagnostics
241        if (idiagaa .gt. 0) then
242            write(*,*)
243            write(*,*) 'attribute #, name, type, value'
244        end if
245        do iatt = 1, natts_tot
246            iret = nf_inq_attname( ncid, nf_global, iatt, attname)
247            if (iret .ne. nf_noerr) goto 3600
248            iret = nf_inq_att( ncid, nf_global, attname, ivtype, lenatt )
249            if (iret .ne. nf_noerr) goto 3600
250            if (ivtype .eq. 2) then
251                iret = nf_get_att_text( ncid, nf_global, attname, dumch1000 )
252                if (iret .ne. nf_noerr) goto 3600
253                i = max(1,min(1000,lenatt))
254                if (idiagaa .gt. 0) write(*,91010) &
255                        iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
256            else if (ivtype .eq. 4) then
257                iret = nf_get_att_int( ncid, nf_global, attname, idum )
258                if (iret .ne. nf_noerr) goto 3600
259                if (idiagaa .gt. 0) write(*,91020) &
260                        iatt, attname(1:40), ivtype, lenatt, idum
261            else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
262                iret = nf_get_att_real( ncid, nf_global, attname, duma )
263                if (iret .ne. nf_noerr) goto 3600
264                if (idiagaa .gt. 0) write(*,91030) &
265                        iatt, attname(1:40), ivtype, lenatt, duma
266            else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
267                allocate( duma_alloc(lenatt) )
268                iret = nf_get_att_real( ncid, nf_global, attname, duma_alloc )
269                if (iret .ne. nf_noerr) goto 3600
270                if (idiagaa .gt. 0) then
271                    write(*,91010) iatt, attname(1:40), ivtype, lenatt
272                    write(*,91040) (duma_alloc(i), i=1,lenatt)
273                end if
274                deallocate( duma_alloc )
275            else
276                if (idiagaa .gt. 0) write(*,'(i4,1x,a,2(1x,i6))')  &
277                        iatt, attname(1:40), ivtype, lenatt
278                goto 3400
279            endif
280
281            if (attname .eq. 'GRID_ID') then
282                m_grid_id = idum
283            else if (attname .eq. 'PARENT_ID') then
284                m_parent_grid_id = idum
285            else if (attname .eq. 'PARENT_GRID_RATIO') then
286                m_parent_grid_ratio = idum
287            else if (attname .eq. 'I_PARENT_START') then
288                i_parent_start = idum
289            else if (attname .eq. 'J_PARENT_START') then
290                j_parent_start = idum
291            else if (attname .eq. 'DX') then
292                dx_met = duma
293            else if (attname .eq. 'DY') then
294                dy_met = duma
295            else if (attname .eq. 'MAP_PROJ') then
296                map_proj_id = idum
297            else if (attname .eq. 'STAND_LON') then
298                map_stdlon = duma
299            else if (attname .eq. 'TRUELAT1') then
300                map_truelat1 = duma
301            else if (attname .eq. 'TRUELAT2') then
302                map_truelat2 = duma
303            else if (attname .eq. 'BL_PBL_PHYSICS') then
304                pbl_physics  = idum
305            else if (attname .eq. 'MP_PHYSICS') then
306                microphysics  = idum
307            end if
308        enddo
3093400    continue
31091010   format( i4, 1x, a, 2(1x,i6), 1x, a )
31191020   format( i4, 1x, a, 2(1x,i6), 1x, i10 )
31291030   format( i4, 1x, a, 2(1x,i6), 1x, 1pe12.4 )
31391040   format(( 12x, 5(1pe12.4) ))
314
315        goto 3900
316
3173600    write(*,9110) 'error inquiring attribute', iatt, fnamenc
318        stop
319
3203900    continue
321
322!C        write(*,*)'mp_physics=',microphysics,pbl_physics
323
324!
325!   close and return
326!
327        iret = nf_close( ncid )
328        ierr = 0
329
330        return
331        end subroutine read_ncwrfout_gridinfo
332
333
334
335!-----------------------------------------------------------------------
336        subroutine read_ncwrfout_1datetime( ierr, fnamenc, &
337          itime, jyyyymmdd, jhhmmss )
338!
339!   a wrf output file may contain data at multiple time.  This routine returns
340!       the date & time of the "itime" data group in the file.
341!
342!   arguments
343!       ierr - output - if non-zero, an error occurred
344!               while opening or reading from the file,
345!               or itime < 0, or itime > number of times in the file.
346!       fnamenc - input - path+filename of the wrf output file
347!       itime - input - specifies which data/time to return. 
348!               1 for first, 2 for second, ...
349!       jyyyymmdd - output - date as 8 decimal digits (yyyymmdd).
350!               yyyy=year, mm=month, dd=day of month.
351!       jhhmmss - output - time of day as 6 decimal digits (hhmmss).
352!               hh=hour, mm=minute, ss=second
353!       if (jyyyymmdd=jhhmmss=-1, then ierr is non-zero, and vice-versa)
354!
355
356!        use netcdf
357        include 'netcdf.inc'
358!implicit none
359
360
361!   arguments
362        integer :: ierr, itime, jyyyymmdd, jhhmmss
363!       character*(*) fnamenc
364        character(len=160) :: fnamenc
365
366!   local variables
367        integer,parameter :: ndims_maxbb=4 ! max number of dimensions for a variable
368
369        integer :: i, id_var, iret, itype_var
370        integer :: iduma, idumb, idumc
371        integer :: id_dim(ndims_maxbb)
372        integer :: istart(ndims_maxbb), icount(ndims_maxbb)
373        integer :: lendim(ndims_maxbb)
374        integer :: natts_tot, ncid, ndims
375
376        character(len=32) timetext
377        character(len=80) varname, varnamenc
378
379!   externals
380!       integer nf_close
381!       integer nf_inq
382!       integer nf_inq_dim
383!       integer nf_open
384
385
386
387        jyyyymmdd = -1
388        jhhmmss = -1
389
390!
391!   open the netcdf file
392!
393        ncid = 10
394        iret = nf_open( fnamenc, NF_NOWRITE, ncid )
395        if (iret .ne. nf_noerr) then
396            write(*,9100) 'error doing open 370', fnamenc
397            ierr = -1
398            goto 8100
399        end if
400
4019100    format( / '*** read_ncwrfout_1datetime -- ', a / &
402        'file = ', a )
4039110  format( / '*** read_ncwrfout_1datetime -- ', a, 1x, i8 / &
404        'file = ', a )
4059120  format( / '*** read_ncwrfout_1datetime -- ', a, 2(1x,i8) / &
406        'file = ', a )
4079130  format( / '*** read_ncwrfout_1datetime -- ', a, 3(1x,i8) / &
408        'file = ', a )
4099115  format( / '*** read_ncwrfout_1datetime -- ', a / a, 1x, i8 / &
410        'file = ', a )
4119125  format( / '*** read_ncwrfout_1datetime -- ', a / a, 2(1x,i8) / &
412        'file = ', a )
4139135  format( / '*** read_ncwrfout_1datetime -- ', a / a, 3(1x,i8) / &
414        'file = ', a )
415
41690030  format( a, 2i6, 2(2x,a) )
417 
418!
419! get information on the variable
420!
421  varname = 'Times'
422        iret = nf_inq_varid( ncid, varname, id_var )
423        if (iret .ne. nf_noerr) then
424            write(*,9100) 'error inquiring var id for ' // varname, fnamenc
425            ierr = -1
426            goto 8100
427        end if
428
429        iret = nf_inq_var( ncid, id_var,  &
430                varnamenc, itype_var, ndims, id_dim, natts_tot )
431        if (iret .ne. nf_noerr) then
432            write(*,9100) 'error inquiring var info for ' // varname, fnamenc
433            ierr = -1
434            goto 8100
435        end if
436
437!   check variable type
438        if (itype_var .ne. nf_char) then
439            write(*,9110) 'var type wrong for ' // varname,  &
440                itype_var, fnamenc
441            ierr = -1
442            goto 8100
443        end if
444
445
446!   check number of dimensions
447        if (ndims .ne. 2) then
448            write(*,9115) 'var ndims is wrong for ' // varname,  &
449                'ndims =', ndims, fnamenc
450            ierr = -1
451            goto 8100
452        end if
453
454!   get sizes of dimensions
455!   dimension 1 = # of characters in date/time string
456!   dimension 2 = # of times in the file
457        do i = 1, ndims
458            iret = nf_inq_dimlen( ncid, id_dim(i), lendim(i) )
459            if (iret .ne. nf_noerr) then
460                write(*,9115) 'error inquiring var dim len for ' // varname, &
461                        'idim =', i, fnamenc
462                ierr = -1
463                goto 8100
464            end if
465        end do
466
467        if (itime .lt. 1) then
468            ierr = -11
469            goto 8100
470        else if (itime .gt. lendim(2)) then
471            ierr = -12
472            goto 8100
473        end if
474
475!   get the data and extract the data & time
476        do i = 1, ndims_maxbb
477            istart(i) = 1
478            icount(i) = 1
479        end do
480        istart(1) = 1
481        icount(1) = lendim(1)
482        istart(2) = itime
483        icount(2) = 1
484        iret = nf_get_vara_text( ncid, id_var, istart, icount, timetext )
485        if (iret .ne. nf_noerr) then
486            write(*,9100) 'error reading var data for ' // varname,  &
487                fnamenc
488            ierr = -1
489            goto 8100
490        end if
491
492        read( timetext, '(i4,1x,i2,1x,i2)', iostat=iret )  &
493                iduma, idumb, idumc
494        if (iret .ne. 0) then
495            write(*,9125) &
496                'error reading from timetext = "' // timetext // '"', &
497                'itime, lendim(1) =', itime, lendim(1), fnamenc
498            ierr = -1
499            goto 8100
500        end if
501        jyyyymmdd = iduma*10000 + idumb*100 + idumc
502
503        read( timetext, '(11x,i2,1x,i2,1x,i2)', iostat=iret )  &
504                iduma, idumb, idumc
505        if (iret .ne. 0) then
506            write(*,9125) &
507                'error reading from timetext = "' // timetext // '"', &
508                'itime, lendim(1) =', itime, lendim(1), fnamenc
509            ierr = -1
510            goto 8100
511        end if
512        jhhmmss = iduma*10000 + idumb*100 + idumc
513
514!
515!   success - close and return
516!
517        iret = nf_close( ncid )
518        ierr = 0
519        return
520
521!
522!   failure - close and return
523!
5248100    iret = nf_close( ncid )
525        return
526
527        end subroutine read_ncwrfout_1datetime
528
529
530
531!-----------------------------------------------------------------------
532        subroutine read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
533          varname, vardata, &
534          itime, &
535          ndims, ndims_exp, ndims_max, &
536          lendim, lendim_exp, lendim_max )
537!
538!   reads of real (single precision) field at one time from a netcdf wrf output file
539!
540!   arguments
541!       ierr - output - if non-zero, an error occurred
542!               while opening or reading from the file
543!                -1 = error opening file
544!                -2 = requested variable is not in the file
545!                -3 = error while inquiring about the variable
546!                -4 = variable type is other than real (single precision)
547!               ... = check below, in the code, for explanation of other ierr values.
548!       idiagaa - input - if positive, testing diagnostics are printed
549!       fnamenc - input - path+filename of the wrf output file
550!       varname - input - field name
551!       vardata - output - the data for the field
552!       itime - input - specifies which time to read. 
553!               (1 for first time in the file, 2 for second, ...)
554!       ndims - output - number of (netcdf) dimensions of the field. 
555!               This includes the time dimension.
556!       ndims_exp - input - expected number of dimensions of the field. 
557!               An error occurs if ndims .ne. ndims_exp.
558!       ndims_max - input - The dimension/size of the lendim_... arrays.
559!       lendim - output - The size of each of the "ndims" dimensions.
560!       lendim_exp - input - The expected size of each of the "ndims" dimensions.
561!               If lendim_exp .gt. 0, then an error occurs if lendim .ne. lendim_exp.
562!       lendim_max - input - The maximum size of each dimension.  These are
563!               The dimensions of the vardata array.
564!
565!        use netcdf
566        include 'netcdf.inc'
567!implicit none
568
569
570!   arguments
571        integer :: ierr, idiagaa, itime, &
572          ndims, ndims_exp, ndims_max, &
573          lendim(ndims_max), lendim_exp(ndims_max),   &
574          lendim_max(ndims_max)
575
576        real :: vardata( lendim_max(1), lendim_max(2), lendim_max(3) )
577
578!       character*(*) fnamenc, varname
579        character(len=160) :: fnamenc, varname
580
581!   local variables
582        integer,parameter :: ndims_maxbb=4 ! max number of dimensions for a variable
583        integer,parameter :: ibadaa=-987
584        integer,parameter :: xbadaa=-987
585        integer,parameter :: ijktestmax=500
586
587        integer :: i, iatt, id_var, iret, itype_var
588        integer :: itot, ijktot, ii
589        integer :: id_dim(ndims_maxbb)
590        integer :: istart(ndims_maxbb), icount(ndims_maxbb)
591        integer :: j, jtot, jj
592        integer :: k, ktot, kk
593        integer :: l, lenatt
594        integer :: lendim_use(ndims_maxbb)
595        integer :: m
596        integer :: natts_tot, ncid, ndiffa, ndiffb
597
598        real :: duma, dumb
599        real :: testavg(ijktestmax,3,2), &
600            testmin(ijktestmax,3,2), testmax(ijktestmax,3,2)
601
602        character(len=80) varnamenc
603        character(len=80) dimname(ndims_maxbb)
604
605!   externals
606!       integer nf_close
607!       integer nf_inq
608!       integer nf_inq_dim
609!       integer nf_open
610
611!        print*,'dim array',lendim_max(1), lendim_max(2), lendim_max(3)
612
613!
614!   open the netcdf file
615!
616       print*,'file name',fnamenc
617        ncid = 10
618        iret = nf_open( fnamenc, NF_NOWRITE, ncid )
619
620        if (iret .ne. nf_noerr) then
621            write(*,9100) 'error doing open 592', fnamenc
622            ierr = -1
623            return
624        end if
625
6269100  format( / '*** read_ncwrfout_1realfield -- ', a / &
627        'file = ', a )
6289110  format( / '*** read_ncwrfout_1realfield -- ', a, 1x, i8 / &
629        'file = ', a )
6309120  format( / '*** read_ncwrfout_1realfield -- ', a, 2(1x,i8) / &
631        'file = ', a )
6329130  format( / '*** read_ncwrfout_1realfield -- ', a, 3(1x,i8) / &
633        'file = ', a )
6349115  format( / '*** read_ncwrfout_1realfield -- ', a / a, 1x, i8 / &
635       'file = ', a )
6369125  format( / '*** read_ncwrfout_1realfield -- ', a / a, 2(1x,i8) / &
637        'file = ', a )
6389135  format( / '*** read_ncwrfout_1realfield -- ', a / a, 3(1x,i8) / &
639        'file = ', a )
640
64190030  format( a, 2i6, 2(2x,a) )
642
643        print*,'variable begin',varname 
644!
645! get information on the variable
646!
647        iret = nf_inq_varid( ncid, varname, id_var )
648        if (iret .ne. nf_noerr) then
649!           write(*,9100) 'error inquiring var id for ' // varname, fnamenc
650            ierr = -2
651            goto 8100
652        end if
653
654        iret = nf_inq_var( ncid, id_var,  &
655                varnamenc, itype_var, ndims, id_dim, natts_tot )
656        if (iret .ne. nf_noerr) then
657            write(*,9100) 'error inquiring var info for ' // varname, fnamenc
658            ierr = -3
659            goto 8100
660        end if
661
662!   check variable type
663        if (itype_var .ne. nf_real) then
664            write(*,9110) 'var type wrong for ' // varname,  &
665                itype_var, fnamenc
666            ierr = -4
667            goto 8100
668        end if
669
670
671!   check number of dimensions
672        if (ndims_exp .le. 0) then
673            write(*,9115)  &
674                'bad ndims_exp for ' // varname,  &
675                'ndims_exp =', ndims_exp, fnamenc
676            ierr = -11
677            goto 8100
678        end if
679        if (ndims .ne. ndims_exp) then
680            write(*,9125) 'var ndims mismatch for ' // varname,  &
681                'ndims_exp, ndims =', ndims_exp, ndims, fnamenc
682            ierr = -12
683            goto 8100
684        end if
685        if (ndims .gt. ndims_max) then
686            write(*,9125) 'var ndims > ndims_max for ' // varname,  &
687                'ndims, ndims_max =', ndims, ndims_max, fnamenc
688            ierr = -13
689            goto 8100
690        end if
691        if (ndims .gt. ndims_maxbb) then
692            write(*,9125) 'var ndims > ndims_maxbb for ' // varname,  &
693                'ndims, ndims_maxbb =', ndims, ndims_maxbb, fnamenc
694            ierr = -14
695            goto 8100
696        end if
697
698!   check size of each dimension
699        do i = 1, ndims_exp
700            iret = nf_inq_dimlen( ncid, id_dim(i), lendim(i) )
701            if (iret .ne. nf_noerr) then
702                write(*,9110) 'error inquiring var dim len for ' // varname, &
703                        i, fnamenc
704                ierr = -15
705                goto 8100
706            end if
707            if ((i .lt. ndims_exp) .and. (lendim_exp(i) .gt. 0) .and.  &
708                (lendim(i) .ne. lendim_exp(i))) then
709!           print*,i,ndims_exp,lendim_exp(i),lendim(i),lendim_exp(i)
710                write(*,9130) 'var lendim mismatch for ' // varname,  &
711                    i, lendim_exp(i), lendim(i), fnamenc
712                ierr = -16
713                goto 8100
714            end if
715            if ((i .lt. ndims_exp) .and.  &
716                (lendim(i) .gt. lendim_max(i))) then
717                write(*,9130) 'var lendim too big for ' // varname,  &
718                    i, lendim_max(i), lendim(i), fnamenc
719                ierr = -17
720                goto 8100
721            end if
722            if ((i .eq. ndims_exp) .and. (lendim(i) .lt. itime)) then 
723                write(*,9130) 'var itime < ntimes for ' // varname,  &
724                    i, itime, lendim(i), fnamenc
725                ierr = -18
726                goto 8100
727            end if
728        end do
729
730!   do diagnostics on the dimensions
731        if (idiagaa .gt. 0) then
732          write(*,'(/a)')  &
733                'read_ncwrfout_1realfield - dim info for var = ' //  &
734                varname(1:20)
735          do i = 1, ndims
736            iret = nf_inq_dim( ncid, id_dim(i), dimname(i), lendim(i) )
737            if (iret .ne. nf_noerr) then
738                write(*,9115) 'error inquiring var dim info for ' // varname, &
739                        'idim =', i, fnamenc
740                ierr = -19
741                goto 8100
742            end if
743            write(*,'(a,3i5,2x,a)') '     i,id,len,name =',  &
744                i, id_dim(i), lendim(i), dimname(i)(1:32)
745          end do
746        end if
747
748!
749!   get the data
750!
751        do i = 1, ndims_maxbb
752            istart(i) = 1
753            icount(i) = 1
754        end do
755        do i = 1, ndims_exp - 1
756            istart(i) = 1
757            icount(i) = lendim(i)
758        end do
759!   in wrfout files, the last dimension should always be time
760        istart(ndims_exp) = itime
761        icount(ndims_exp) = 1
762
763        iret = nf_get_vara_real( ncid, id_var, istart, icount, vardata )
764        if (iret .ne. nf_noerr) then
765            write(*,9120) 'error reading var data for ' // varname,  &
766                fnamenc
767            ierr = -21
768            goto 8100
769        end if
770
771!
772!   reorder the data
773!
774        lendim_use(1) = lendim(1)
775        lendim_use(2) = 1
776        if (ndims_exp .ge. 3) lendim_use(2) = lendim(2)
777        lendim_use(3) = 1
778        if (ndims_exp .ge. 4) lendim_use(3) = lendim(3)
779
780        print*,'variable',varname
781!        print*,'value in',lendim_use, lendim_max
782        call reorder_ncwrfout_1realfield( ierr, idiagaa,  &
783          varname, vardata, vardata, &
784          lendim_use, lendim_max )
785
786        if (ierr .ne. 0) then
787            write(*,9120) 'error re-ordering var data for ' // varname,  &
788                fnamenc
789            ierr = -22
790            goto 8100
791        end if
792
793!
794!   success - close and return
795!
796        iret = nf_close( ncid )
797        ierr = 0
798        return
799
800!
801!   error - close and return
802!
8038100    iret = nf_close( ncid )
804        return
805
806        end subroutine read_ncwrfout_1realfield
807
808
809
810!-----------------------------------------------------------------------
811        subroutine reorder_ncwrfout_1realfield( ierr, idiagaa,  &
812          varname, vardata_in, vardata_out, &
813          lendim_use, lendim_max )
814!
815!   reorders a real (single precision) field
816!       the nf_get_vara_real loads the data for a field into
817!           a contiguous block of memory starting at vardata(1,1,1)
818!       it does not know if perhaps lendim() < lendim_max()
819!       this routine corrects for that, so that the data are
820!           loaded into non-contiguous blocks when lendim() < lendim_max()
821!
822!   arguments
823!       ierr - output - if non-zero, an error occurred while re-ordering
824!       idiagaa - input - if positive, testing diagnostics are printed
825!       varname - input - field name
826!
827!       vardata_in  - input  - the data for the field
828!       vardata_out - output - the data for the field
829!               In the calling program, vardata_in & vardata_out are usually
830!               the same array.  This routine "pretends" that they are
831!               different, and specifies their dimensions differently,
832!               to facilitate the reordering.
833!
834!       lendim_use - input - The actual size of the spatial dimensions of the field.
835!       lendim_max - input - The actual spatial dimensions of the vardata array.
836!               Most wrf fields are spatially either 1d (z), 2d (xy), or 3d (xyz).
837!               For a 1d spatial field (e.g., z only), set
838!                   lendim_use(1) = nz,  lendim_max(1) = nz_max
839!                   lendim_use(2) = 1,   lendim_max(2) = 1
840!                   lendim_use(3) = 1,   lendim_max(3) = 1
841!               For a 2d spatial field (e.g., xy only), set
842!                   lendim_use(1) = nx,  lendim_max(1) = nx_max
843!                   lendim_use(2) = ny   lendim_max(2) = ny_max
844!                   lendim_use(3) = 1,   lendim_max(3) = 1
845!               For a 3d spatial field (xyz), set
846!                   lendim_use(1) = nx,  lendim_max(1) = nx_max
847!                   lendim_use(2) = ny   lendim_max(2) = ny_max
848!                   lendim_use(3) = nz   lendim_max(3) = nz_max
849!
850
851!        use netcdf
852        include 'netcdf.inc'
853!implicit none
854
855
856!   arguments
857        integer :: ierr, idiagaa,  &
858          lendim_use(3), lendim_max(3)
859
860        real :: vardata_in(  lendim_use(1), lendim_use(2), lendim_use(3) )
861        real :: vardata_out( lendim_max(1), lendim_max(2), lendim_max(3) )
862
863        character*(*) varname
864
865!   local variables
866        integer,parameter :: ijktestmax=500
867        integer,parameter :: check_reordering=1
868
869        integer :: i, j, k, m, n
870        integer :: itestend, jtestend, ktestend
871        integer :: ijk, ijktestend
872        integer :: ndiffa, ndiffb
873
874        real :: duma, dumb
875        real :: testavg(ijktestmax,3,2),  &
876            testmin(ijktestmax,3,2), testmax(ijktestmax,3,2)
877
878!
879!   the testavg/min/max are avg, min, and max values for
880!       a i (or j or k) fixed and j,k (or i,k or i,j) varying
881!   they are computed before and after the data has been reordered
882!       then compared at the end
883!   an error occurs if they do not match
884!
885!       print*,'size in out',lendim_use(1), lendim_use(2), lendim_use(3)
886!       print*,'size in out',lendim_max(1), lendim_max(2), lendim_max(3)
887        if (check_reordering .gt. 0) then
888
889        do n = 1, 2
890        do m = 1, 3
891        do i = 1, ijktestmax
892            testavg(i,m,n) = 0.0
893            testmin(i,m,n) = +1.0e37
894            testmax(i,m,n) = -1.0e37
895        end do
896        end do
897        end do
898!          print*,varname
899        ktestend = min( ijktestmax, lendim_use(3) )
900        jtestend = min( ijktestmax, lendim_use(2) )
901        itestend = min( ijktestmax, lendim_use(1) )
902
903!   pass 1 -- compute the test---(:,:,1) from vardata_in
904        do k = 1, ktestend
905        do j = 1, jtestend
906        do i = 1, itestend
907            duma = vardata_in(i,j,k)
908            testavg(i,1,1) =      testavg(i,1,1) + duma
909            testmin(i,1,1) = min( testmin(i,1,1),  duma )
910            testmax(i,1,1) = max( testmax(i,1,1),  duma )
911            testavg(j,2,1) =      testavg(j,2,1) + duma
912            testmin(j,2,1) = min( testmin(j,2,1),  duma )
913            testmax(j,2,1) = max( testmax(j,2,1),  duma )
914            testavg(k,3,1) =      testavg(k,3,1) + duma
915            testmin(k,3,1) = min( testmin(k,3,1),  duma )
916            testmax(k,3,1) = max( testmax(k,3,1),  duma )
917        end do
918        end do
919        end do
920
921        end if    ! if (check_reordering .gt. 0) then
922
923!   pass 2 -- shift the data values
924!         print*,'max',lendim_max
925        do k = lendim_use(3), 1, -1
926        do j = lendim_use(2), 1, -1
927        do i = lendim_use(1), 1, -1
928!         print*,i,j,k
929            vardata_out(i,j,k) = vardata_in(i,j,k)
930        end do
931        end do
932        end do
933
934!   pass 3 -- compute the test---(:,:,2) from vardata_out
935        if (check_reordering .gt. 0) then
936
937        do k = 1, ktestend
938        do j = 1, jtestend
939        do i = 1, itestend
940            duma = vardata_out(i,j,k)
941            testavg(i,1,2) =      testavg(i,1,2) + duma
942            testmin(i,1,2) = min( testmin(i,1,2),  duma )
943            testmax(i,1,2) = max( testmax(i,1,2),  duma )
944            testavg(j,2,2) =      testavg(j,2,2) + duma
945            testmin(j,2,2) = min( testmin(j,2,2),  duma )
946            testmax(j,2,2) = max( testmax(j,2,2),  duma )
947            testavg(k,3,2) =      testavg(k,3,2) + duma
948            testmin(k,3,2) = min( testmin(k,3,2),  duma )
949            testmax(k,3,2) = max( testmax(k,3,2),  duma )
950        end do
951        end do
952        end do
953
954!   now compare the test---(:,:,1) & test---(:,:,2)
955        ndiffb = 0
956        do m = 1, 3
957            if (m .eq. 1) then
958                ijktestend = itestend
959                duma = 1.0/(jtestend*ktestend)
960                if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
961                        'i, testavg(i,1), testmin(i,1), testmax(i,1)'
962            else if (m .eq. 2) then
963                ijktestend = jtestend
964                duma = 1.0/(itestend*ktestend)
965                if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
966                        'j, testavg(j,2), testmin(j,2), testmax(j,2)'
967            else
968                ijktestend = ktestend
969                duma = 1.0/(itestend*jtestend)
970                if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
971                        'k, testavg(k,3), testmin(k,3), testmax(k,3)'
972            end if
973
974            ndiffa = 0
975            do ijk = 1, ijktestend
976                i = ijk
977                dumb = max( abs(testavg(i,m,1)), abs(testavg(i,m,2)) )*2.0e-7
978                if (abs(testavg(i,m,1)-testavg(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
979                dumb = max( abs(testmin(i,m,1)), abs(testmin(i,m,2)) )*2.0e-7
980                if (abs(testmin(i,m,1)-testmin(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
981                dumb = max( abs(testmax(i,m,1)), abs(testmax(i,m,2)) )*2.0e-7
982                if (abs(testmax(i,m,1)-testmax(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
983            end do
984
985            if (ndiffa .le. 0) then
986                if (idiagaa .gt. 0) write(*,*) '     *** no differences'
987            else
988              do ijk = 1, ijktestend
989                i = ijk
990                if (idiagaa .gt. 0) write(*,'(i3,1p,3(2x,2e11.3))') i, &
991                    testavg(i,m,1)*duma,  &
992                    (testavg(i,m,1)-testavg(i,m,2))*duma,  &
993                    testmin(i,m,1), (testmin(i,m,1)-testmin(i,m,2)), &
994                    testmax(i,m,1), (testmax(i,m,1)-testmax(i,m,2))
995              end do
996            end if
997
998            ndiffb = ndiffb + ndiffa
999        end do
1000
1001        if (ndiffb .gt. 0) then
1002            ierr = -12
1003            goto 8100
1004        end if
1005
1006        end if    ! if (check_reordering .gt. 0) then
1007
1008!
1009!   success
1010!
1011        ierr = 0
1012        return
1013
1014!
1015!   error
1016!
10178100    return
1018
1019        end subroutine reorder_ncwrfout_1realfield
1020
1021
1022
hosted by ZAMG