?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / var / mri4dvar / da_vp_bilin.f90
blob3855e0f95f9ebd15432b55d51cc598de59ebae41
1 program da_vp_bilin
2
3 !----------------------------------------------------------------------
4 ! Purpose: Regridding from low to high resolution in control variable space
5 ! by using bilinear interpolation
6 !
7 ! where n is the grid number in x or y
8 ! ns is the refinement ratio between two resulotions
9 !
10 ! Method: follow da_bilin.f90
11 !
12 ! Compile:
13 !
14 ! pgf90 -o da_vp_bilin.exe da_vp_bilin.f90
15 !
16 ! liuz@ucar.edu , 2016-08, NCAR/MMM
17 !----------------------------------------------------------------------
18
19 !use netcdf
20
21 implicit none
22
23 !These variables' incremental will be regridded by default
24 character (len=6), dimension(1:19) :: vNam
25
26 integer :: ix, jy, kz, k, status
27 integer :: ixh, jyh, kzh
28 integer :: nLat, nLon, oLat, oLon
29 integer :: sLat, eLat, sLon, eLon
30 integer :: rLat, rLon
31
32 real, dimension(:,:,:), allocatable :: v1, v2, v3, v4, v5
33 real, dimension(:,:,:), allocatable :: v6, v7, v8, v9, v10, v11
34 real, dimension(:,:,:), allocatable :: v1h, v2h, v3h, v4h, v5h
35 real, dimension(:,:,:), allocatable :: v6h, v7h, v8h, v9h, v10h, v11h
36 real, dimension(:,:), allocatable :: iVar, oVar
37
38 character (len = 255) :: appname = ""
39 character (len = 255) :: arg = ""
40 character (len = 19) :: analysis_date
41 character (len = 255) :: input_file= "vp_output.global"
42 character (len = 255) :: output_file= "vp_output.global_hires"
43
44 integer, parameter :: vp_unit = 8
45 integer, parameter :: vp_hires_unit = 9
46 integer :: ratio ! resolution ratio
47 integer :: cloud_cv_options ! 2 or 3 with cloud cv variables
48 integer :: use_cv_w ! =1 for w control variable
49 integer :: io_status
50 integer iargc
51
52 LOGICAL :: file_exists
53
54 !These variables' incremental will be regridded by default
55
56 !call getarg(0, appname)
57 !n=index(appname, '/', BACK=.true.)
58 !appname = trim(appname(n+1:))
59
60 call getarg(1, arg)
61 call getarg(2, arg)
62 read(arg, '(i3)') ratio
63
64 call getarg(3, arg)
65 call getarg(4, arg)
66 read(arg, '(i3)') cloud_cv_options
67
68 call getarg(5, arg)
69 call getarg(6, arg)
70 read(arg, '(i3)') use_cv_w
71
72
73 write (*, *) 'ratio = ', ratio, 'cloud_cv_options = ', cloud_cv_options, &
74 'use_cv_w = ', use_cv_w
75
76
77 ! read vp file
78 !--------------------
79 inquire(FILE=trim(input_file), EXIST=file_exists)
80
81 if ( .not. file_exists ) then
82 Write(*,*) "\nError: "//trim(input_file)//" not exists\n"
83 call exit(-1)
84 else
85 Write(*,*) "Found: "//trim(input_file)
86 endif
87
88 open(unit=vp_unit,file=trim(input_file),iostat=io_status,form='UNFORMATTED',status='OLD')
89 if (io_status /= 0) then
90 write(*,*) "Error ",io_status," opening vp file "//trim(input_file)
91 call exit(-1)
92 end if
93 write(*,*) 'Reading vp from : '//trim(input_file)
94 !read(vp_unit) analysis_date
95 !print *, 'analysis_date = ', analysis_date
96 read(vp_unit) ix, jy, kz ! domain dimension (unstagered)
97 print *, "input file: ix, jy, kz = ", ix, jy, kz
98
99 allocate ( v1 (1:ix,1:jy,1:kz))
100 allocate ( v2 (1:ix,1:jy,1:kz))
101 allocate ( v3 (1:ix,1:jy,1:kz))
102 allocate ( v4 (1:ix,1:jy,1:kz))
103 allocate ( v5 (1:ix,1:jy,1:kz))
104
105 read(vp_unit) v1, v2, v3, v4, v5
106
107 if ( cloud_cv_options >= 2 ) then
108 allocate ( v6 (1:ix,1:jy,1:kz))
109 allocate ( v7 (1:ix,1:jy,1:kz))
110 allocate ( v8 (1:ix,1:jy,1:kz))
111 allocate ( v9 (1:ix,1:jy,1:kz))
112 allocate ( v10 (1:ix,1:jy,1:kz))
113 read(vp_unit) v6, v7, v8, v9, v10
114 end if
115
116 if ( use_cv_w == 1 ) then
117 allocate ( v11 (1:ix,1:jy,1:kz))
118 read(vp_unit) v11
119 end if
120
121 write(*,*) 'End Reading vp from : '//trim(input_file)
122 close(vp_unit)
123 !-----------------------------
124 ! end read vp file
125 !----------------------
126
127 nLon = ix + 2 ! 52
128 nLat = jy + 2 ! 52
129
130 rLon = ix * ratio ! 150
131 rLat = jy * ratio ! 150
132
133 oLon = ( nLon - 1 ) * ratio + 1 ! 154
134 oLat = ( nLat - 1 ) * ratio + 1
135
136 elat = (oLat - rLat) / 2 ! 2
137 slat = oLat - rLat - elat + 1 ! 3
138
139 elon = (oLon - rLon) / 2
140 slon = oLon - rLon - elon + 1
141
142 allocate(iVar(nLon, nLat), stat=status)
143 allocate(oVar(oLon, oLat), stat=status)
144
145
146 ixh = ix*ratio
147 jyh = jy*ratio
148
149 allocate ( v1h (1:ixh,1:jyh,1:kz))
150 allocate ( v2h (1:ixh,1:jyh,1:kz))
151 allocate ( v3h (1:ixh,1:jyh,1:kz))
152 allocate ( v4h (1:ixh,1:jyh,1:kz))
153 allocate ( v5h (1:ixh,1:jyh,1:kz))
154
155 if ( cloud_cv_options >= 2 ) then
156 allocate ( v6h (1:ixh,1:jyh,1:kz))
157 allocate ( v7h (1:ixh,1:jyh,1:kz))
158 allocate ( v8h (1:ixh,1:jyh,1:kz))
159 allocate ( v9h (1:ixh,1:jyh,1:kz))
160 allocate ( v10h (1:ixh,1:jyh,1:kz))
161 end if
162
163 if ( use_cv_w == 1 ) then
164 allocate ( v11h (1:ixh,1:jyh,1:kz))
165 end if
166
167 do k = 1, kz
168 iVar(2:nlon-1,2:nlat-1) = v1(:,:,k)
169 iVar(1,:) = iVar(2,:)
170 iVar(nlon,:) = iVar(nlon-1,:)
171 iVar(:,1) = iVar(:,2)
172 iVar(:,nlat) = iVar(:,nlat-1)
173 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
174 v1h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
175
176 iVar(2:nlon-1,2:nlat-1) = v2(:,:,k)
177 iVar(1,:) = iVar(2,:)
178 iVar(nlon,:) = iVar(nlon-1,:)
179 iVar(:,1) = iVar(:,2)
180 iVar(:,nlat) = iVar(:,nlat-1)
181 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
182 v2h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
183
184 iVar(2:nlon-1,2:nlat-1) = v3(:,:,k)
185 iVar(1,:) = iVar(2,:)
186 iVar(nlon,:) = iVar(nlon-1,:)
187 iVar(:,1) = iVar(:,2)
188 iVar(:,nlat) = iVar(:,nlat-1)
189 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
190 v3h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
191
192 iVar(2:nlon-1,2:nlat-1) = v4(:,:,k)
193 iVar(1,:) = iVar(2,:)
194 iVar(nlon,:) = iVar(nlon-1,:)
195 iVar(:,1) = iVar(:,2)
196 iVar(:,nlat) = iVar(:,nlat-1)
197 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
198 v4h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
199
200 iVar(2:nlon-1,2:nlat-1) = v5(:,:,k)
201 iVar(1,:) = iVar(2,:)
202 iVar(nlon,:) = iVar(nlon-1,:)
203 iVar(:,1) = iVar(:,2)
204 iVar(:,nlat) = iVar(:,nlat-1)
205 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
206 v5h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
207
208 if ( cloud_cv_options >= 2 ) then
209 iVar(2:nlon-1,2:nlat-1) = v6(:,:,k)
210 iVar(1,:) = iVar(2,:)
211 iVar(nlon,:) = iVar(nlon-1,:)
212 iVar(:,1) = iVar(:,2)
213 iVar(:,nlat) = iVar(:,nlat-1)
214 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
215 v6h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
216
217 iVar(2:nlon-1,2:nlat-1) = v7(:,:,k)
218 iVar(1,:) = iVar(2,:)
219 iVar(nlon,:) = iVar(nlon-1,:)
220 iVar(:,1) = iVar(:,2)
221 iVar(:,nlat) = iVar(:,nlat-1)
222 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
223 v7h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
224
225 iVar(2:nlon-1,2:nlat-1) = v8(:,:,k)
226 iVar(1,:) = iVar(2,:)
227 iVar(nlon,:) = iVar(nlon-1,:)
228 iVar(:,1) = iVar(:,2)
229 iVar(:,nlat) = iVar(:,nlat-1)
230 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
231 v8h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
232
233 iVar(2:nlon-1,2:nlat-1) = v9(:,:,k)
234 iVar(1,:) = iVar(2,:)
235 iVar(nlon,:) = iVar(nlon-1,:)
236 iVar(:,1) = iVar(:,2)
237 iVar(:,nlat) = iVar(:,nlat-1)
238 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
239 v9h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
240
241 iVar(2:nlon-1,2:nlat-1) = v10(:,:,k)
242 iVar(1,:) = iVar(2,:)
243 iVar(nlon,:) = iVar(nlon-1,:)
244 iVar(:,1) = iVar(:,2)
245 iVar(:,nlat) = iVar(:,nlat-1)
246 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
247 v10h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
248 end if
249
250 if ( use_cv_w == 1 ) then
251 iVar(2:nlon-1,2:nlat-1) = v11(:,:,k)
252 iVar(1,:) = iVar(2,:)
253 iVar(nlon,:) = iVar(nlon-1,:)
254 iVar(:,1) = iVar(:,2)
255 iVar(:,nlat) = iVar(:,nlat-1)
256 call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
257 v11h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
258 end if
259 enddo
260
261 open(unit=vp_hires_unit,file=trim(output_file),iostat=io_status,form='UNFORMATTED',status='UNKNOWN')
262 if (io_status /= 0) then
263 write(*,*) "Error ",io_status," opening vp file "//trim(output_file)
264 call exit(-1)
265 end if
266 write(*,*) 'Writting vp on hires. to : '//trim(output_file)
267
268 print *, 'output file: ixh, jyh, kz=', ixh, jyh, kz
269 write(vp_hires_unit) ixh, jyh, kz
270 write(vp_hires_unit) v1h,v2h,v3h,v4h,v5h
271 if ( cloud_cv_options >= 2 ) then
272 write(vp_hires_unit) v6h,v7h,v8h,v9h,v10h
273 end if
274 if ( use_cv_w == 1 ) then
275 write(vp_hires_unit) v11h
276 end if
277
278 deallocate(v1, stat=status)
279 deallocate(v2, stat=status)
280 deallocate(v3, stat=status)
281 deallocate(v4, stat=status)
282 deallocate(v5, stat=status)
283
284 deallocate(v1h, stat=status)
285 deallocate(v2h, stat=status)
286 deallocate(v3h, stat=status)
287 deallocate(v4h, stat=status)
288 deallocate(v5h, stat=status)
289
290 if ( cloud_cv_options >= 2 ) then
291 deallocate(v6, stat=status)
292 deallocate(v7, stat=status)
293 deallocate(v8, stat=status)
294 deallocate(v9, stat=status)
295 deallocate(v10, stat=status)
296
297 deallocate(v6h, stat=status)
298 deallocate(v7h, stat=status)
299 deallocate(v8h, stat=status)
300 deallocate(v9h, stat=status)
301 deallocate(v10h, stat=status)
302 end if
303
304 if ( use_cv_w == 1 ) then
305 deallocate(v11, stat=status)
306 deallocate(v11h, stat=status)
307 end if
308
309 Write(*,*) "Regridding increment completed successfully"
310
311 contains
312 subroutine show_usage()
313 Write(*,*) 'Usage :'//trim(appname)// &
314 '[-h] [-fg_lores filename] [-an_lores filename] [-fg_hires filename] [-ns n ] [-o outputfile]'
315 Write(*,*) " -fg_lores Optional, low resulotion first guess file, default - fg"
316 Write(*,*) " -an_lores Optional, low resulotion analysis file comes from wrfvar, default - wrfvar_output"
317 Write(*,*) " -fg_hires Optional, high resultion first guess file, default - wrfinput_hires"
318 Write(*,*) " -ns Optional, the refinement ratio between two resulotions, default - 3"
319 Write(*,*) " -o Optional, output high resulotion analysis file, default - wrfvar_output_hires"
320 Write(*,*) " -h Show this help"
321 end subroutine show_usage
322
323 !subroutine nf90_handle_err(status, errmsg)
324 ! integer, intent(in) :: status
325 ! character(len=*), intent(in) :: errmsg
326 !
327 ! if(status /= nf90_noerr) then
328 ! print *, trim(nf90_strerror(status))//" : "//trim(errmsg)
329 ! Stop
330 ! end if
331 ! end subroutine nf90_handle_err
332
333 subroutine bilin(old,xi,yi,ns,new,xo,yo)
334
335 ! assume: xo = (xi-1)*ns + 1, xi=50, xo=49*3+1=148
336 ! yo = (yi-1)*ns + 1
337
338 implicit none
339
340 integer, intent(in) :: xi,yi,xo,yo
341 real, dimension(xi,yi), intent(in) :: old
342 integer, intent(in) :: ns
343 real, dimension(xo,yo), intent(out):: new
344
345 real :: im(1:ns+1,2)
346 ! real :: imm(1:ns+3,2)
347 integer:: i,j,jm1,im1,ix1,ix2,iy1,iy2
348
349 forall(i=1:ns+1) im(i,2) = real(i-1)/ns
350 im(:,1) = 1 - im(:,2)
351
352 do j=2,yi
353 jm1 = j - 1
354 iy2 = jm1 * ns + 1
355 iy1 = iy2 - ns
356 do i=2,xi
357 im1 = i - 1
358 ix2 = im1 * ns + 1
359 ix1 = ix2 - ns
360 new(ix1:ix2,iy1:iy2) = matmul(im,matmul(old(im1:i,jm1:j),transpose(im)))
361 end do
362 end do
363
364
365 ! ns = ns + 2
366 ! forall(i=1:ns+1) imm(i,2) = real(i-1)/ns
367 ! imm(:,1) = 1 - imm(:,2)
368 !
369 ! j=yi
370 ! jm1 = j - 1
371 ! iy2 = jm1 * ns + 1
372 ! iy1 = iy2 - ns
373 !
374 ! i=xi
375 ! im1 = i - 1
376 ! ix2 = im1 * ns + 1
377 ! ix1 = ix2 - ns
378 ! new(ix1:ix2,iy1:iy2) = matmul(imm,matmul(old(im1:i,jm1:j),transpose(imm)))
379 ! end do
380 ! end do
381
382 end subroutine bilin
383
384 end program da_vp_bilin
Mirror of the official repository at https://github.com/wrf-model/WRF