vort.F90
References to this file elsewhere.
1 ! on linux, compile wrf then compile as:
2 ! pgf90 -Mfree -I ../../main -I ../../inc -I /usr/local/netcdf-pgi/include vort.F90 libwrfio_nf.a /usr/local/netcdf-pgi/lib/libnetcdf.a ../../main/libwrflib.a
3 ! on AIX, compile wrf then compile as:
4 ! /lib/cpp -C -P vort.F90 > vort.f
5 ! mpxlf -qfree=f90 -I ../../share -I ../../main -I ../../inc -I /usr/local/netcdf/include vort.f libwrfio_nf.a /usr/local/netcdf/lib/libnetcdf.a ../../main/libwrflib.a
6
7 module read_util_module
8
9 #ifdef crayx1
10 #define iargc ipxfargc
11 #endif
12
13 contains
14
15 #ifdef crayx1
16 subroutine getarg(i, harg)
17 implicit none
18 character(len=*) :: harg
19 integer :: ierr, ilen, i
20
21 call pxfgetarg(i, harg, ilen, ierr)
22 return
23 end subroutine getarg
24 #endif
25
26 subroutine arguments(v2file, lmore)
27 implicit none
28 character(len=*) :: v2file
29 character(len=120) :: harg
30 logical :: lmore
31
32 integer :: ierr, i, numarg
33 integer, external :: iargc
34
35 numarg = iargc()
36
37 i = 1
38 lmore = .false.
39
40 do while ( i < numarg)
41 call getarg(i, harg)
42 print*, 'harg = ', trim(harg)
43
44 if (harg == "-v") then
45 i = i + 1
46 lmore = .true.
47 elseif (harg == "-h") then
48 call help
49 endif
50
51 enddo
52
53 call getarg(i,harg)
54 v2file = harg
55 end subroutine arguments
56
57 subroutine help
58 implicit none
59 character(len=120) :: cmd
60 call getarg(0, cmd)
61
62 write(*,'(/,"Usage: ", A, " [-v] v2file ")') trim(cmd)
63 write(*,'(8x, "-v : Print extra info")')
64 write(*,'(8x, "v3file : MM5v3 file name to read.")')
65 write(*,'(8x, "-h : print this help message and exit.",/)')
66 stop
67 end subroutine help
68 end module read_util_module
69
70
71
72 program readv3
73 use wrf_data
74 use read_util_module
75 use module_compute_geop
76
77
78 implicit none
79 #include "wrf_status_codes.h"
80 #include <netcdf.inc>
81 character(len=120) :: flnm
82 character(len=120) :: flnm2
83 character(len=120) :: arg3
84 character(len=19) :: DateStr
85 character(len=19) :: DateStr2
86 character(len=31) :: VarName
87 character(len=31) :: VarName2
88 integer dh1, dh2
89
90 integer :: flag, flag2
91 integer :: iunit, iunit2
92
93 integer :: i,j,k
94 integer :: levlim
95 integer :: cross
96 integer :: ndim, ndim2
97 integer :: WrfType, WrfType2
98 real :: time, time2
99 real*8 :: a, b
100 real*8 :: sum1, sum2, diff1, diff2, serr, perr, rms
101 integer, dimension(4) :: start_index, end_index, start_index2, end_index2, end_index_u, end_index_uz
102 integer , Dimension(3) :: MemS,MemE,PatS,PatE
103 character (len= 4) :: staggering, staggering2
104 character (len= 3) :: ordering, ordering2, ord
105 character (len=24) :: start_date, start_date2
106 character (len=24) :: current_date, current_date2
107 character (len=31) :: name, name2, tmpname
108 character (len=25) :: units, units2
109 character (len=46) :: description, description2
110
111 real, allocatable, dimension(:,:,:) :: ph, phb, p, pb
112 real, allocatable, dimension(:,:) :: height
113
114 integer :: ids, ide, jds, jde, kds, kde, &
115 ims, ime, jms, jme, kms, kme, &
116 its, ite, jts, jte, kts, kte
117 integer outcount
118
119
120 character (len=80), dimension(3) :: dimnames
121 character (len=80) :: SysDepInfo
122
123 integer :: l, n
124 integer :: ikdiffs, ifdiffs
125
126 real, allocatable, dimension(:,:,:,:) :: data,data2
127
128 integer :: ierr, ierr2, ier, ier2, Status, Status_next_time, Status_next_time2, Status_next_var, Status_next_var_2
129
130 logical :: newtime = .TRUE.
131 logical :: justplot, efound
132
133 integer, external :: iargc
134 logical, external :: iveceq
135
136 levlim = -1
137
138 call ext_ncd_ioinit(SysDepInfo,Status)
139 call set_wrf_debug_level ( 1 )
140
141
142 Justplot = .true.
143
144 ! get arguments
145 ! if ( iargc() .ge. 2 ) then
146 call getarg(1,flnm)
147 ! call getarg(2,flnm2)
148 ierr = 0
149 call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
150 if ( Status /= 0 ) then
151 print*,'error opening ',flnm, ' Status = ', Status ; stop
152 endif
153 ! call ext_ncd_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
154 ! if ( Status /= 0 ) go to 923
155 ! goto 924
156 !923 continue
157 !
158 !! bounce here if second name is not openable -- this would mean that
159 !! it is a field name instead.
160 !
161 ! print*,'could not open ',flnm2
162 ! name = flnm2
163 ! Justplot = .true.
164 !924 continue
165 ! if ( iargc() .eq. 3 ) then
166 ! call getarg(3,arg3)
167 ! read(arg3,*)levlim
168 ! print*,'LEVLIM = ',LEVLIM
169 ! endif
170 ! else
171 ! print*,'Usage: command file1 file2'
172 ! stop
173 ! endif
174
175 !print*,'Just plot ',Justplot
176
177 start_index = 1
178 end_index = 0
179
180 CALL ext_ncd_get_dom_ti_integer(dh1,'WEST-EAST_GRID_DIMENSION',end_index(1),1,OutCount,Status)
181 CALL ext_ncd_get_dom_ti_integer(dh1,'BOTTOM-TOP_GRID_DIMENSION',end_index(2),1,OutCount,Status)
182 CALL ext_ncd_get_dom_ti_integer(dh1,'SOUTH-NORTH_GRID_DIMENSION',end_index(3),1,OutCount,Status)
183 ord = 'XZY'
184 staggering = ' '
185
186
187
188 allocate(ph(end_index(1),end_index(2),end_index(3)))
189 allocate(phb(end_index(1),end_index(2),end_index(3)))
190 allocate(p(end_index(1),end_index(2),end_index(3)))
191 allocate(pb(end_index(1),end_index(2),end_index(3)))
192 allocate(height(end_index(1),end_index(3)))
193
194 ids=start_index(1); ide=end_index(1); jds=start_index(3); jde=end_index(3); kds=start_index(2); kde=end_index(2)
195 ims=start_index(1); ime=end_index(1); jms=start_index(3); jme=end_index(3); kms=start_index(2); kme=end_index(2)
196 its=start_index(1); ite=end_index(1)-1; jts=start_index(3); jte=end_index(3)-1; kts=start_index(2); kte=end_index(2)-1
197
198 end_index_u = end_index - 1
199 end_index_uz = end_index - 1
200 end_index_uz(2) = end_index_uz(2) + 1
201
202
203
204 if ( Justplot ) then
205 print*, 'flnm = ', trim(flnm)
206
207 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
208
209 DO WHILE ( Status_next_time .eq. 0 )
210 write(*,*)'Next Time ',TRIM(Datestr)
211
212 staggering = 'Z'
213 name = 'PH'
214 call ext_ncd_read_field(dh1,DateStr,TRIM(name),ph,WRF_REAL,0,0,0,ord, &
215 staggering, dimnames , &
216 start_index,end_index_uz, & !dom
217 start_index,end_index, & !mem
218 start_index,end_index_uz, & !pat
219 ierr)
220 name = 'PHB'
221 call ext_ncd_read_field(dh1,DateStr,TRIM(name),phb,WRF_REAL,0,0,0,ord, &
222 staggering, dimnames , &
223 start_index,end_index_uz, & !dom
224 start_index,end_index, & !mem
225 start_index,end_index_uz, & !pat
226 ierr)
227 staggering = ' '
228 name = 'P'
229 call ext_ncd_read_field(dh1,DateStr,TRIM(name),p,WRF_REAL,0,0,0,ord, &
230 staggering, dimnames , &
231 start_index,end_index_u, & !dom
232 start_index,end_index, & !mem
233 start_index,end_index_u, & !pat
234 ierr)
235 name = 'PB'
236 call ext_ncd_read_field(dh1,DateStr,TRIM(name),pb,WRF_REAL,0,0,0,ord, &
237 staggering, dimnames , &
238 start_index,end_index_u, & !dom
239 start_index,end_index, & !mem
240 start_index,end_index_u, & !pat
241 ierr)
242
243 CALL compute_500mb_height ( ph, phb, p, pb, &
244 height, &
245 ids, ide, jds, jde, kds, kde, &
246 ims, ime, jms, jme, kms, kme, &
247 its, ite, jts, jte, kts, kte )
248
249 write(88,*)end_index_u(1),end_index_u(3),' height ',trim(Datestr)
250 do j = 1, end_index_u(3)
251 do i = 1, end_index_u(1)
252 write(88,*) height(i,j)
253 enddo
254 enddo
255
256 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
257 enddo
258 endif
259
260 end program readv3
261
262 ! stub for routine called by module_wrf_error (used by netcdf implementation of IO api)
263 SUBROUTINE wrf_abort
264 STOP
265 END SUBROUTINE wrf_abort