da_update_bc.f90
References to this file elsewhere.
1 program da_update_bc
2
3 !-----------------------------------------------------------------------
4 ! Purpose: update BC file from wrfvar output.
5 ! current version reads only wrf-netcdf file format
6 !-----------------------------------------------------------------------
7
8 use da_netcdf_interface, only : da_get_var_3d_real_cdf, &
9 da_put_var_3d_real_cdf, da_get_dims_cdf, da_put_var_2d_real_cdf, &
10 da_get_var_2d_real_cdf, da_get_var_2d_int_cdf, da_get_bdytimestr_cdf, &
11 da_get_times_cdf, da_get_bdyfrq, stderr, stdout
12
13 use da_module_couple_uv, only : da_couple_uv
14
15 implicit none
16
17 integer, parameter :: max_3d_variables = 20, &
18 max_2d_variables = 20
19
20 character(len=512) :: wrfvar_output_file, &
21 wrf_bdy_file, &
22 wrf_input
23
24 character(len=20) :: var_pref, var_name, vbt_name
25
26 character(len=20) :: var3d(max_3d_variables), &
27 var2d(max_2d_variables)
28
29 character(len=10), dimension(4) :: bdyname, tenname
30
31 integer :: ids, ide, jds, jde, kds, kde
32 integer :: num3d, num2d, ndims
33 integer :: time_level
34 integer :: i,j,k,l,m,n
35
36 integer, dimension(4) :: dims
37
38 real, allocatable, dimension(:,:,:) :: tend3d, scnd3d, frst3d, full3d
39
40 real, allocatable, dimension(:,:,:) :: u, v
41
42 real, allocatable, dimension(:, :) :: mu, mub, msfu, msfv, msfm, &
43 tend2d, scnd2d, frst2d, full2d
44
45 real, allocatable, dimension(:, :) :: tsk, tsk_wrfvar
46
47 integer, allocatable, dimension(:,:) :: ivgtyp
48
49 character(len=80), allocatable, dimension(:) :: times, &
50 thisbdytime, nextbdytime
51
52 integer :: east_end, north_end, io_status
53
54 logical :: cycling, debug, low_bdy_only
55
56 real :: bdyfrq
57
58 integer, parameter :: namelist_unit = 7, &
59 ori_unit = 11, &
60 new_unit = 12
61
62 namelist /control_param/ wrfvar_output_file, &
63 wrf_bdy_file, &
64 wrf_input, &
65 cycling, debug, low_bdy_only
66
67 wrfvar_output_file = 'wrfvar_output'
68 wrf_bdy_file = 'wrfbdy_d01'
69 wrf_input = 'wrfinput_d01'
70
71 cycling = .false.
72 debug = .false.
73 low_bdy_only = .false.
74 !---------------------------------------------------------------------
75 ! Read namelist
76 !---------------------------------------------------------------------
77 io_status = 0
78
79 open(unit = namelist_unit, file = 'parame.in', &
80 status = 'old' , access = 'sequential', &
81 form = 'formatted', action = 'read', &
82 iostat = io_status)
83
84 if (io_status /= 0) then
85 write(unit=stdout,fmt=*) 'Error to open namelist file: parame.in.'
86 write(unit=stdout,fmt=*) 'Will work for updating lateral boundary only.'
87 else
88 read(unit=namelist_unit, nml = control_param , iostat = io_status)
89
90 if (io_status /= 0) then
91 write(unit=stdout,fmt=*) 'Error to read control_param. Stopped.'
92 stop
93 end if
94
95 WRITE(unit=stdout, fmt='(2a)') &
96 'wrfvar_output_file = ', trim(wrfvar_output_file), &
97 'wrf_bdy_file = ', trim(wrf_bdy_file), &
98 'wrf_input = ', trim(wrf_input)
99
100 WRITE(unit=stdout, fmt='(a, L10)') &
101 'cycling = ', cycling
102
103 close(unit=namelist_unit)
104 end if
105
106 ! 3D need update
107 num3d=6
108 var3d(1)='U'
109 var3d(2)='V'
110 var3d(3)='W'
111 var3d(4)='T'
112 var3d(5)='PH'
113 var3d(6)='QVAPOR'
114
115 ! 2D need update
116 num2d=10
117 var2d(1)='MUB'
118 var2d(2)='MU'
119 var2d(3)='MAPFAC_U'
120 var2d(4)='MAPFAC_V'
121 var2d(5)='MAPFAC_M'
122 var2d(6)='TMN'
123 var2d(7)='SST'
124 var2d(8)='TSK'
125 var2d(9)='VEGFRA'
126 var2d(10)='ALBBCK'
127
128 ! First, the boundary times
129 call da_get_dims_cdf(wrf_bdy_file, 'Times', dims, ndims, debug)
130
131 if (debug) then
132 write(unit=stdout, fmt='(a,i2,2x,a,4i6)') &
133 'Times: ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
134 end if
135
136 time_level = dims(2)
137
138 if (time_level < 1) then
139 write(unit=stdout, fmt='(a,i2/a)') &
140 'time_level = ', time_level, &
141 'We need at least one time-level BDY.'
142 stop 'Wrong BDY file.'
143 end if
144
145 allocate(times(dims(2)))
146 allocate(thisbdytime(dims(2)))
147 allocate(nextbdytime(dims(2)))
148
149 call da_get_times_cdf(wrf_bdy_file, times, dims(2), dims(2), debug)
150
151 call da_get_bdytimestr_cdf(wrf_bdy_file, 'thisbdytime', thisbdytime, dims(2), debug)
152 call da_get_bdytimestr_cdf(wrf_bdy_file, 'nextbdytime', nextbdytime, dims(2), debug)
153
154 call da_get_bdyfrq(thisbdytime(1), nextbdytime(1), bdyfrq, debug)
155
156 if (debug) then
157 do n=1, dims(2)
158 write(unit=stdout, fmt='(3(a, i2, 2a,2x))') &
159 ' times(', n, ')=', trim(times(n)), &
160 'thisbdytime (', n, ')=', trim(thisbdytime(n)), &
161 'nextbdytime (', n, ')=', trim(nextbdytime(n))
162 end do
163 end if
164
165 east_end=0
166 north_end=0
167
168 ! For 2D variables
169 ! Get mu, mub, msfu, and msfv
170
171 do n=1,num2d
172 call da_get_dims_cdf( wrfvar_output_file, trim(var2d(n)), dims, &
173 ndims, debug)
174
175 select case(trim(var2d(n)))
176 case ('MU') ;
177 if (low_bdy_only) cycle
178
179 allocate(mu(dims(1), dims(2)))
180
181 call da_get_var_2d_real_cdf( wrfvar_output_file, &
182 trim(var2d(n)), mu, dims(1), dims(2), 1, debug)
183
184 east_end=dims(1)+1
185 north_end=dims(2)+1
186 case ('MUB') ;
187 if (low_bdy_only) cycle
188
189 allocate(mub(dims(1), dims(2)))
190
191 call da_get_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), mub, &
192 dims(1), dims(2), 1, debug)
193 case ('MAPFAC_U') ;
194 if (low_bdy_only) cycle
195
196 allocate(msfu(dims(1), dims(2)))
197
198 call da_get_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), msfu, &
199 dims(1), dims(2), 1, debug)
200 case ('MAPFAC_V') ;
201 if (low_bdy_only) cycle
202
203 allocate(msfv(dims(1), dims(2)))
204
205 call da_get_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), msfv, &
206 dims(1), dims(2), 1, debug)
207 case ('MAPFAC_M') ;
208 if (low_bdy_only) cycle
209
210 allocate(msfm(dims(1), dims(2)))
211
212 call da_get_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), msfm, &
213 dims(1), dims(2), 1, debug)
214 case ('TSK') ;
215 if (.not. cycling) cycle
216
217 allocate(tsk(dims(1), dims(2)))
218 allocate(tsk_wrfvar(dims(1), dims(2)))
219 allocate(ivgtyp(dims(1), dims(2)))
220
221 call da_get_var_2d_real_cdf( wrf_input, trim(var2d(n)), tsk, &
222 dims(1), dims(2), 1, debug)
223 call da_get_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), tsk_wrfvar, &
224 dims(1), dims(2), 1, debug)
225 call da_get_var_2d_int_cdf( wrfvar_output_file, 'IVGTYP', ivgtyp, &
226 dims(1), dims(2), 1, debug)
227
228 ! update TSK.
229 do j=1,dims(2)
230 do i=1,dims(1)
231 if (ivgtyp(i,j) /= 16) &
232 tsk(i,j)=tsk_wrfvar(i,j)
233 end do
234 end do
235
236 call da_put_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), tsk, &
237 dims(1), dims(2), 1, debug)
238 deallocate(tsk)
239 deallocate(ivgtyp)
240 deallocate(tsk_wrfvar)
241 case ('TMN', 'SST', 'VEGFRA', 'ALBBCK') ;
242 if (.not. cycling) cycle
243
244 allocate(full2d(dims(1), dims(2)))
245
246 call da_get_var_2d_real_cdf( wrf_input, trim(var2d(n)), full2d, &
247 dims(1), dims(2), 1, debug)
248
249 call da_put_var_2d_real_cdf( wrfvar_output_file, trim(var2d(n)), full2d, &
250 dims(1), dims(2), 1, debug)
251 deallocate(full2d)
252 case default ;
253 write(unit=stdout,fmt=*) 'It is impossible here. var2d(n)=', trim(var2d(n))
254 end select
255 end do
256
257 if (low_bdy_only) then
258 write(unit=stdout,fmt=*) 'only low boundary updated.'
259 stop
260 end if
261
262 if (east_end < 1 .or. north_end < 1) then
263 write(unit=stdout, fmt='(a)') 'Wrong data for Boundary.'
264 stop
265 end if
266
267 ! boundary variables
268 bdyname(1)='_BXS'
269 bdyname(2)='_BXE'
270 bdyname(3)='_BYS'
271 bdyname(4)='_BYE'
272
273 ! boundary tendancy variables
274 tenname(1)='_BTXS'
275 tenname(2)='_BTXE'
276 tenname(3)='_BTYS'
277 tenname(4)='_BTYE'
278
279 do m=1,4
280 var_name='MU' // trim(bdyname(m))
281 vbt_name='MU' // trim(tenname(m))
282
283 call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
284
285 allocate(frst2d(dims(1), dims(2)))
286 allocate(scnd2d(dims(1), dims(2)))
287 allocate(tend2d(dims(1), dims(2)))
288
289 ! Get variable at second time level
290 if (time_level > 1) then
291 call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), scnd2d, &
292 dims(1), dims(2), 2, debug)
293 else
294 call da_get_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
295 dims(1), dims(2), 1, debug)
296 call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
297 dims(1), dims(2), 1, debug)
298 end if
299
300 if (debug) then
301 write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
302 'No.', m, 'Variable: ', trim(vbt_name), &
303 'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
304
305 call da_get_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
306 dims(1), dims(2), 1, debug)
307
308 write(unit=ori_unit, fmt='(a, 10i12)') &
309 ' old ', (i, i=1,dims(2))
310 do j=1,dims(1)
311 write(unit=ori_unit, fmt='(i4, 1x, 5e20.7)') &
312 j, (tend2d(j,i), i=1,dims(2))
313 end do
314 end if
315
316 ! calculate variable at first time level
317 select case(m)
318 case (1) ; ! West boundary
319 do l=1,dims(2)
320 do j=1,dims(1)
321 if (time_level < 2) &
322 scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
323 frst2d(j,l)=mu(l,j)
324 end do
325 end do
326 case (2) ; ! East boundary
327 do l=1,dims(2)
328 do j=1,dims(1)
329 if (time_level < 2) &
330 scnd2d(j,l)=frst2d(j,l)+tend2d(j,l)*bdyfrq
331 frst2d(j,l)=mu(east_end-l,j)
332 end do
333 end do
334 case (3) ; ! South boundary
335 do l=1,dims(2)
336 do i=1,dims(1)
337 if (time_level < 2) &
338 scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
339 frst2d(i,l)=mu(i,l)
340 end do
341 end do
342 case (4) ; ! North boundary
343 do l=1,dims(2)
344 do i=1,dims(1)
345 if (time_level < 2) &
346 scnd2d(i,l)=frst2d(i,l)+tend2d(i,l)*bdyfrq
347 frst2d(i,l)=mu(i,north_end-l)
348 end do
349 end do
350 case default ;
351 write(unit=stdout,fmt=*) 'It is impossible here. mu, m=', m
352 end select
353
354 ! calculate new tendancy
355 do l=1,dims(2)
356 do i=1,dims(1)
357 tend2d(i,l)=(scnd2d(i,l)-frst2d(i,l))/bdyfrq
358 end do
359 end do
360
361 if (debug) then
362 write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
363 'No.', m, 'Variable: ', trim(vbt_name), &
364 'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
365
366 write(unit=new_unit, fmt='(a, 10i12)') &
367 ' new ', (i, i=1,dims(2))
368
369 do j=1,dims(1)
370 write(unit=new_unit, fmt='(i4, 1x, 5e20.7)') &
371 j, (tend2d(j,i), i=1,dims(2))
372 end do
373 end if
374
375 ! output new variable at first time level
376 call da_put_var_2d_real_cdf( wrf_bdy_file, trim(var_name), frst2d, &
377 dims(1), dims(2), 1, debug)
378 ! output new tendancy
379 call da_put_var_2d_real_cdf( wrf_bdy_file, trim(vbt_name), tend2d, &
380 dims(1), dims(2), 1, debug)
381
382 deallocate(frst2d)
383 deallocate(scnd2d)
384 deallocate(tend2d)
385 end do
386
387 !---------------------------------------------------------------------
388 ! For 3D variables
389
390 ! Get U
391 call da_get_dims_cdf( wrfvar_output_file, 'U', dims, ndims, debug)
392
393 ! call da_get_att_cdf( wrfvar_output_file, 'U', debug)
394
395 allocate(u(dims(1), dims(2), dims(3)))
396
397 ids=1
398 ide=dims(1)-1
399 jds=1
400 jde=dims(2)
401 kds=1
402 kde=dims(3)
403
404 call da_get_var_3d_real_cdf( wrfvar_output_file, 'U', u, &
405 dims(1), dims(2), dims(3), 1, debug)
406
407 ! do j=1,dims(2)
408 ! write(unit=stdout, fmt='(2(a,i5), a, f12.8)') &
409 ! 'u(', dims(1), ',', j, ',1)=', u(dims(1),j,1)
410 ! end do
411
412 ! Get V
413 call da_get_dims_cdf( wrfvar_output_file, 'V', dims, ndims, debug)
414
415 ! call da_get_att_cdf( wrfvar_output_file, 'V', debug)
416
417 allocate(v(dims(1), dims(2), dims(3)))
418
419 call da_get_var_3d_real_cdf( wrfvar_output_file, 'V', v, &
420 dims(1), dims(2), dims(3), 1, debug)
421
422 ! do i=1,dims(1)
423 ! write(unit=stdout, fmt='(2(a,i5), a, f12.8)') &
424 ! 'v(', i, ',', dims(2), ',1)=', v(i,dims(2),1)
425 ! end do
426
427 if (debug) then
428 write(unit=stdout, fmt='(a,e20.12,4x)') &
429 'Before couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
430 'Before couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
431 end if
432
433 !---------------------------------------------------------------------
434 ! Couple u, v.
435 call da_couple_uv ( u, v, mu, mub, msfu, msfv, ids, ide, jds, jde, kds, kde)
436
437 if (debug) then
438 write(unit=stdout, fmt='(a,e20.12,4x)') &
439 'After couple Sample u=', u(dims(1)/2,dims(2)/2,dims(3)/2), &
440 'After couple Sample v=', v(dims(1)/2,dims(2)/2,dims(3)/2)
441 end if
442
443 !---------------------------------------------------------------------
444 !For 3D variables
445
446 do n=1,num3d
447 write(unit=stdout, fmt='(a, i3, 2a)') 'Processing: var3d(', n, ')=', trim(var3d(n))
448
449 call da_get_dims_cdf( wrfvar_output_file, trim(var3d(n)), dims, ndims, debug)
450
451 allocate(full3d(dims(1), dims(2), dims(3)))
452
453 east_end=dims(1)+1
454 north_end=dims(2)+1
455
456 select case(trim(var3d(n)))
457 case ('U') ; ! U
458 ! var_pref='R' // trim(var3d(n))
459 var_pref=trim(var3d(n))
460 full3d(:,:,:)=u(:,:,:)
461 case ('V') ; ! V
462 ! var_pref='R' // trim(var3d(n))
463 var_pref=trim(var3d(n))
464 full3d(:,:,:)=v(:,:,:)
465 case ('W') ;
466 ! var_pref = 'R' // trim(var3d(n))
467 var_pref = trim(var3d(n))
468
469 call da_get_var_3d_real_cdf( wrfvar_output_file, trim(var3d(n)), &
470 full3d, dims(1), dims(2), dims(3), 1, debug)
471
472 if (debug) then
473 write(unit=stdout, fmt='(3a,e20.12,4x)') &
474 'Before couple Sample ', trim(var3d(n)), &
475 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
476 end if
477
478 do k=1,dims(3)
479 do j=1,dims(2)
480 do i=1,dims(1)
481 full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))/msfm(i,j)
482 end do
483 end do
484 end do
485
486 if (debug) then
487 write(unit=stdout, fmt='(3a,e20.12,4x)') &
488 'After couple Sample ', trim(var3d(n)), &
489 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
490 end if
491 case ('T', 'PH') ;
492 var_pref=trim(var3d(n))
493
494 call da_get_var_3d_real_cdf( wrfvar_output_file, trim(var3d(n)), &
495 full3d, dims(1), dims(2), dims(3), 1, debug)
496
497 if (debug) then
498 write(unit=stdout, fmt='(3a,e20.12,4x)') &
499 'Before couple Sample ', trim(var3d(n)), &
500 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
501 end if
502
503 do k=1,dims(3)
504 do j=1,dims(2)
505 do i=1,dims(1)
506 full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))
507 end do
508 end do
509 end do
510
511 if (debug) then
512 write(unit=stdout, fmt='(3a,e20.12,4x)') &
513 'After couple Sample ', trim(var3d(n)), &
514 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
515 end if
516 case ('QVAPOR', 'QCLOUD', 'QRAin', 'QICE', 'QSNOW', 'QGRAUP') ;
517 ! var_pref='R' // var3d(n)(1:2)
518 ! var_pref=var3d(n)(1:2)
519 var_pref=var3d(n)
520
521 call da_get_var_3d_real_cdf( wrfvar_output_file, trim(var3d(n)), &
522 full3d, dims(1), dims(2), dims(3), 1, debug)
523
524 if (debug) then
525 write(unit=stdout, fmt='(3a,e20.12,4x)') &
526 'Before couple Sample ', trim(var3d(n)), &
527 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
528 end if
529
530 do k=1,dims(3)
531 do j=1,dims(2)
532 do i=1,dims(1)
533 full3d(i,j,k)=full3d(i,j,k)*(mu(i,j)+mub(i,j))
534 end do
535 end do
536 end do
537
538 if (debug) then
539 write(unit=stdout, fmt='(3a,e20.12,4x)') &
540 'After couple Sample ', trim(var3d(n)), &
541 '=', full3d(dims(1)/2,dims(2)/2,dims(3)/2)
542 end if
543 case default ;
544 write(unit=stdout,fmt=*) 'It is impossible here. var3d(', n, ')=', trim(var3d(n))
545 end select
546
547 do m=1,4
548 var_name=trim(var_pref) // trim(bdyname(m))
549 vbt_name=trim(var_pref) // trim(tenname(m))
550
551 write(unit=stdout, fmt='(a, i3, 2a)') &
552 'Processing: bdyname(', m, ')=', trim(var_name)
553
554 call da_get_dims_cdf( wrf_bdy_file, trim(var_name), dims, ndims, debug)
555
556 allocate(frst3d(dims(1), dims(2), dims(3)))
557 allocate(scnd3d(dims(1), dims(2), dims(3)))
558 allocate(tend3d(dims(1), dims(2), dims(3)))
559
560 ! Get variable at second time level
561 if (time_level > 1) then
562 call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), scnd3d, &
563 dims(1), dims(2), dims(3), 2, debug)
564 else
565 call da_get_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
566 dims(1), dims(2), dims(3), 1, debug)
567 call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
568 dims(1), dims(2), dims(3), 1, debug)
569 end if
570
571 if (debug) then
572 write(unit=ori_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
573 'No.', m, 'Variable: ', trim(vbt_name), &
574 'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
575
576 call da_get_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
577 dims(1), dims(2), dims(3), 1, debug)
578
579 write(unit=ori_unit, fmt='(a, 10i12)') &
580 ' old ', (i, i=1,dims(3))
581 do j=1,dims(1)
582 write(unit=ori_unit, fmt='(i4, 1x, 5e20.7)') &
583 j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
584 end do
585 end if
586
587 select case(trim(bdyname(m)))
588 case ('_BXS') ; ! West boundary
589 do l=1,dims(3)
590 do k=1,dims(2)
591 do j=1,dims(1)
592 if (time_level < 2) &
593 scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
594 frst3d(j,k,l)=full3d(l,j,k)
595 end do
596 end do
597 end do
598 case ('_BXE') ; ! East boundary
599 do l=1,dims(3)
600 do k=1,dims(2)
601 do j=1,dims(1)
602 if (time_level < 2) &
603 scnd3d(j,k,l)=frst3d(j,k,l)+tend3d(j,k,l)*bdyfrq
604 frst3d(j,k,l)=full3d(east_end-l,j,k)
605 end do
606 end do
607 end do
608 case ('_BYS') ; ! South boundary
609 do l=1,dims(3)
610 do k=1,dims(2)
611 do i=1,dims(1)
612 if (time_level < 2) &
613 scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
614 frst3d(i,k,l)=full3d(i,l,k)
615 end do
616 end do
617 end do
618 case ('_BYE') ; ! North boundary
619 do l=1,dims(3)
620 do k=1,dims(2)
621 do i=1,dims(1)
622 if (time_level < 2) &
623 scnd3d(i,k,l)=frst3d(i,k,l)+tend3d(i,k,l)*bdyfrq
624 frst3d(i,k,l)=full3d(i,north_end-l,k)
625 end do
626 end do
627 end do
628 case default ;
629 write(unit=stdout,fmt=*) 'It is impossible here.'
630 write(unit=stdout,fmt=*) 'bdyname(', m, ')=', trim(bdyname(m))
631 stop
632 end select
633
634 write(unit=stdout, fmt='(a, i3, 2a)') &
635 'cal. tend: bdyname(', m, ')=', trim(vbt_name)
636
637 ! calculate new tendancy
638 do l=1,dims(3)
639 do k=1,dims(2)
640 do i=1,dims(1)
641 tend3d(i,k,l)=(scnd3d(i,k,l)-frst3d(i,k,l))/bdyfrq
642 end do
643 end do
644 end do
645
646 if (debug) then
647 write(unit=new_unit, fmt='(a,i2,2x,2a/a,i2,2x,a,4i6)') &
648 'No.', m, 'Variable: ', trim(vbt_name), &
649 'ndims=', ndims, 'dims=', (dims(i), i=1,ndims)
650
651 write(unit=new_unit, fmt='(a, 10i12)') &
652 ' new ', (i, i=1,dims(3))
653
654 do j=1,dims(1)
655 write(unit=new_unit, fmt='(i4, 1x, 5e20.7)') &
656 j, (tend3d(j,dims(2)/2,i), i=1,dims(3))
657 end do
658 end if
659
660 ! output new variable at first time level
661 call da_put_var_3d_real_cdf( wrf_bdy_file, trim(var_name), frst3d, &
662 dims(1), dims(2), dims(3), 1, debug)
663 call da_put_var_3d_real_cdf( wrf_bdy_file, trim(vbt_name), tend3d, &
664 dims(1), dims(2), dims(3), 1, debug)
665
666 deallocate(frst3d)
667 deallocate(scnd3d)
668 deallocate(tend3d)
669 end do
670
671 deallocate(full3d)
672 end do
673
674 deallocate(mu)
675 deallocate(mub)
676 deallocate(msfu)
677 deallocate(msfv)
678 deallocate(u)
679 deallocate(v)
680
681 if (io_status /= 0) then
682 write(unit=stdout,fmt=*) 'only lateral boundary updated.'
683 else
684 if (low_bdy_only) then
685 if (cycling) then
686 write(unit=stdout,fmt=*) 'Both low boudary and lateral boundary updated.'
687 else
688 write(unit=stdout,fmt=*) 'only low boudary updated.'
689 end if
690 else
691 write(unit=stdout,fmt=*) 'only lateral boundary updated.'
692 end if
693 end if
694
695 write (unit=stdout,fmt=*) "*** Update_bc completed successfully ***"
696
697 end program da_update_bc
698