?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / var / mri4dvar / da_bdy.f90
blob8567d04e2b5de8003319dab9696d457eb56813f9
1 program da_bdy
2
3 !----------------------------------------------------------------------
4 ! Purpose: Generates boundary file by using wrfinput
5 !
6 ! Input : fg -- first time level wrfinput generated by real
7 ! fg02 -- second time level wrfinput generated by real
8 ! wrfbdy_ref -- reference boundary file generated by real
9 !
10 ! Output : wrfbdy_out -- the output boundary file
11 !
12 ! Notes : 1. variable name and attributes, dimension name, bdy_width
13 ! come from wrfbdy.
14 ! 2. domain size and time come from fg
15 ! 3. boundary and tendency are calculated by using fg & fg02
16 ! 4. the output boundary file only contain the 1st time level
17 !
18 ! jliu@ucar.edu , 2011-12-15
19 !----------------------------------------------------------------------
20
21 use netcdf
22
23 implicit none
24
25 integer :: i, n, offset, bdyfrq, domainsize, fg_jd, fg02_jd
26
27 integer :: ncid, ncidfg, ncidfg02, ncidwrfbdy, ncidvarbdy, varid, varid_out, status
28 integer :: nDims, nVars, nGlobalAtts, numsAtts
29 integer :: dLen, attLen, xtype, unlimDimID
30 integer :: bdy_width, varbdy_dimID, wrfbdy_dimID, fg_dimID, vTimes_ID, MSF_ID
31 integer :: MU_fgID, MU_fg02ID, MUB_fgID, MUB_fg02ID, fg_varid, fg02_varid, tenid
32
33 integer, dimension(4) :: dsizes
34 integer, dimension(4), target :: start_u, start_v, start_mass
35 integer, dimension(4) :: cnt_4d, map_4d
36 integer, dimension(3) :: start_3d, cnt_3d, map_3d
37 integer, dimension(3), target :: start_msfu, start_msfv, cnt_msfu, cnt_msfv, map_msfu, map_msfv
38 integer, dimension(:), pointer :: start_msf, cnt_msf, map_msf, start_4d
39
40 integer :: south_north, south_north_stag
41 integer :: west_east, west_east_stag
42 integer :: bottom_top, bottom_top_stag
43
44 integer, dimension(nf90_max_var_dims) :: vDimIDs
45 integer, dimension(:), allocatable :: vdimsizes
46 integer, dimension(:,:,:,:), allocatable :: iVar
47
48 real, dimension(:,:,:,:), allocatable :: fVar_fg, fVar_fg02, Tend
49 real, dimension(:,:,:), allocatable , target :: MU_fg, MU_fg02, MUB_fg, MUB_fg02, MSF
50
51 real, dimension(:,:,:), pointer :: MU_fgptr, MU_fg02ptr, MUB_fgptr, MUB_fg02ptr, MSF_ptr
52
53 character (len = 19), dimension(:), allocatable :: times
54 character (len = 19) :: fg_time, fg02_time
55 character (len = 5) :: tenname
56 character (len = NF90_MAX_NAME) :: vNam, dNam, attNam
57 character (len = 9) :: MSF_NAME
58 character (len = 255) :: err_msg=""
59 character (len=8) :: i_char
60 character (len=255) :: arg = ""
61 character (len=255) :: appname =""
62 character (len=255) :: fg = "fg"
63 character (len=255) :: fg02 = "fg02"
64 character (len=255) :: wrfbdy = "wrfbdy_ref"
65 character (len=255) :: varbdy = "wrfbdy_out"
66
67 logical :: reverse, couple, stag
68
69 integer iargc
70
71 call getarg(0, appname)
72 n=index(appname, '/', BACK=.true.)
73 appname = trim(appname(n+1:))
74
75 DO i = 1, iargc(), 2
76 call getarg(i, arg)
77 select case ( trim(arg) )
78 case ("-fg")
79 call getarg(i+1, arg)
80 fg=trim(arg)
81 case ("-fg02")
82 call getarg(i+1, arg)
83 fg02=trim(arg)
84 case ("-bdy")
85 call getarg(i+1, arg)
86 wrfbdy=trim(arg)
87 case ("-o")
88 call getarg(i+1, arg)
89 varbdy=trim(arg)
90 case default
91 Write(*,*) "Usage : "//trim(appname)//" [-fg filename] [-fg02 filename] [-bdy filename] [-o outputfile] [-h]"
92 Write(*,*) " -fg Optional, 1st time levle first guess file, default - fg"
93 Write(*,*) " -fg02 Optional, 2nd time levle first guess file, default - fg02"
94 Write(*,*) " -bdy Optional, reference boundary file comes from real, default - wrfbdy_ref"
95 Write(*,*) " -o Optional, output boundary file, default - varbdy_out"
96 Write(*,*) " -h Show this usage"
97 call exit(0)
98 end select
99 END DO
100
101
102 status = nf90_open(fg, NF90_NOWRITE, ncidfg)
103 if ( status /= nf90_noerr ) then
104 err_msg="Failed to open "//trim(fg)
105 call nf90_handle_err(status, err_msg)
106 endif
107
108 status = nf90_open(fg02, NF90_NOWRITE, ncidfg02)
109 if ( status /= nf90_noerr ) then
110 err_msg="Failed to open "//trim(fg02)
111 call nf90_handle_err(status, err_msg)
112 endif
113
114 status = nf90_inq_varid(ncidfg, "Times", vTimes_ID )
115 if ( status /= nf90_noerr ) then
116 err_msg="Please make sure fg has a vaild Times variable"
117 call nf90_handle_err(status, err_msg)
118 endif
119
120 status = nf90_get_var(ncidfg, vTimes_ID, fg_time)
121 if ( status /= nf90_noerr ) then
122 err_msg="Please make sure fg has a vaild Time value"
123 call nf90_handle_err(status, err_msg)
124 endif
125
126 status = nf90_inq_varid(ncidfg02, "Times", vTimes_ID )
127 if ( status /= nf90_noerr ) then
128 err_msg="Please make sure fg02 has a vaild Times variable"
129 call nf90_handle_err(status, err_msg)
130 endif
131
132 status = nf90_get_var(ncidfg02, vTimes_ID, fg02_time)
133 if ( status /= nf90_noerr ) then
134 err_msg="Please make sure fg02 has a vaild Time value"
135 call nf90_handle_err(status, err_msg)
136 endif
137
138 status = nf90_open(wrfbdy, NF90_NOWRITE, ncidwrfbdy)
139 if ( status /= nf90_noerr ) then
140 err_msg="Failed to open "//trim(wrfbdy)
141 call nf90_handle_err(status, err_msg)
142 endif
143
144 status = nf90_create(varbdy, NF90_CLOBBER, ncidvarbdy)
145 if ( status /= nf90_noerr ) then
146 err_msg="Please make sure have write access"
147 call nf90_handle_err(status, err_msg)
148 endif
149
150 bdyfrq = datediff(fg_time, fg02_time)
151
152 select case ( bdyfrq )
153 case ( 0 )
154 bdyfrq = 1
155 case ( : -1 )
156 Write (*,*) "***WARNNING : time levle of fg is LATER then fg02's.***"
157 end select
158
159 write(i_char, '(i8)') bdyfrq
160
161 Write(*,*) " Input :"
162 Write(*,*) " fg "//fg_time
163 Write(*,*) " fg02 "//fg02_time
164 Write(*,*) " Reference bdy "//trim(wrfbdy)
165 Write(*,*) "Output : "
166 Write(*,*) " wrfbdy_out "//fg_time
167 Write(*,*) " bdyfrq ",adjustl(i_char)
168
169 status = nf90_inquire(ncidfg, nAttributes=nGlobalAtts)
170 do i=1, nGlobalAtts
171 status = nf90_inq_attname(ncidfg, NF90_GLOBAL, i, attNam)
172 status = nf90_copy_att(ncidfg, NF90_GLOBAL, attNam, ncidvarbdy, NF90_GLOBAL)
173 end do
174
175 status = nf90_inquire(ncidwrfbdy, nDims, nVars, nGlobalAtts, unlimDimID)
176 if ( status /= nf90_noerr ) then
177 err_msg="Please make sure have a valid wrf boundary file"
178 call nf90_handle_err(status, err_msg)
179 endif
180
181 allocate (vdimsizes(nDims), stat=status)
182
183 do i=1, nDims
184
185 status = nf90_inquire_dimension(ncidwrfbdy, i, name=dNam, len = dLen)
186
187 vdimsizes(i) = dLen
188 select case (trim(dNam))
189 case ("south_north")
190 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
191 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
192 vdimsizes(i) = dLen
193 south_north = vdimsizes(i)
194 case ("west_east")
195 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
196 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
197 vdimsizes(i) = dLen
198 west_east = vdimsizes(i)
199 case ("south_north_stag")
200 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
201 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
202 vdimsizes(i) = dLen
203 south_north_stag = vdimsizes(i)
204 case ("west_east_stag")
205 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
206 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
207 vdimsizes(i) = dLen
208 west_east_stag = vdimsizes(i)
209 case ("bottom_top")
210 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
211 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
212 vdimsizes(i) = dLen
213 bottom_top = vdimsizes(i)
214 case ("bottom_top_stag")
215 status = nf90_inq_dimid(ncidfg, dNam, fg_dimID)
216 status = nf90_inquire_dimension(ncidfg, fg_dimID, len=dLen)
217 vdimsizes(i) = dLen
218 bottom_top_stag = vdimsizes(i)
219 case ("Time")
220 vdimsizes(i) = 1
221 allocate(times(vdimsizes(i)), stat=status)
222 case ("bdy_width")
223 bdy_width = dLen
224 end select
225
226 if ( i == unlimDimID ) dLen = NF90_UNLIMITED
227
228 status = nf90_def_dim(ncidvarbdy, dNam, dLen, varbdy_dimID)
229
230 end do
231
232 status = nf90_inq_varid(ncidfg , "MU" , MU_fgID )
233 status = nf90_inq_varid(ncidfg , "MUB", MUB_fgID )
234 status = nf90_inq_varid(ncidfg02, "MU" , MU_fg02ID )
235 status = nf90_inq_varid(ncidfg02, "MUB", MUB_fg02ID)
236
237 status = nf90_inq_varid(ncidfg, "Times", vTimes_ID )
238
239 do varid=1, nVars
240
241 status = nf90_inquire_variable(ncidwrfbdy,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs,natts=numsAtts)
242 status = nf90_def_var(ncidvarbdy, trim(vNam), xtype, vDimIDs(1:nDims), varid_out)
243 if ( status /= nf90_noerr ) then
244 err_msg="Failed to define variable : "//trim(vNam)
245 call nf90_handle_err(status, err_msg)
246 endif
247
248 do i=1, numsAtts
249 status = nf90_inq_attname(ncidwrfbdy, varid, i, attNam)
250 status = nf90_copy_att(ncidwrfbdy, varid, trim(attNam), ncidvarbdy, varid_out)
251 if ( status /= nf90_noerr ) then
252 err_msg="Failed to copy att : "//trim(attNam)
253 call nf90_handle_err(status, err_msg)
254 endif
255 end do
256
257 end do
258
259 status = nf90_enddef(ncidvarbdy)
260
261 do varid=1, nVars
262
263 status = nf90_inquire_variable(ncidwrfbdy,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs)
264 if ( status /= nf90_noerr ) then
265 err_msg="Failed to inquire varialbe '"//trim(vNam)//"' for wrfbdy"
266 call nf90_handle_err(status, err_msg)
267 endif
268
269 dsizes = 1
270 do i = 1 , nDims
271 dsizes(i) = vdimsizes(vDimIDs(i))
272 end do
273
274 offset = index(vNam, '_', BACK=.True.)
275 if ( offset <= 0 ) offset = Len(Trim(vNam))
276
277 ! fg
278 ! U (west_east_stag, south_north, bottom_top, time)
279 ! V (west_east, south_north_stag, bottom_top, time)
280 ! T, QVAPOR (west_east, south_north, bottom_top, time)
281 ! PH (west_east, south_north, bottom_top_stag, time)
282 ! MU (west_east, south_north, time)
283 ! MAPFAC_U (west_east_stag, south_north, time)
284 ! MAPFAC_V (west_east, south_north_stag, time)
285 ! bdy
286 ! west & east
287 ! U (south_north, bottom_top, bdy_width, time)
288 ! V (south_north_stag, bottom_top, bdy_width, time)
289 ! T, QVAPOR (south_north, bottom_top, bdy_width, time)
290 ! PH (south_north, bottom_top_stag, bdy_width, time)
291 ! MU (south_north, bdy_width, time)
292 ! north & south
293 ! U (west_east_stag, bottom_top, bdy_width, time)
294 ! V (west_east, bottom_top, bdy_width, time)
295 ! T, QVAPOR (west_east, bottom_top, bdy_width, time)
296 ! PH (west_east, bottom_top_stag, bdy_width, time)
297 ! MU (west_east, bdy_width, time)
298
299 select case (Trim(vNam(offset:)))
300 case ("_BXS") ! West Boundary
301 start_u = (/1,1,1,1/)
302 start_v = (/1,1,1,1/)
303 start_mass = (/1,1,1,1/)
304 start_3d = (/1,1,1/)
305 start_msfu = (/1,1,1/)
306 start_msfv = (/1,1,1/)
307
308 cnt_4d = (/dsizes(3),dsizes(1),dsizes(2),1/)
309 cnt_3d = (/bdy_width,south_north,1/)
310 cnt_msfu = (/bdy_width,south_north,1/)
311 cnt_msfv = (/bdy_width,south_north_stag,1/)
312
313 map_4d = (/dsizes(1)*dsizes(2), 1, dsizes(1), dsizes(1)*dsizes(2)*dsizes(3)/)
314 map_3d = (/south_north, 1, bdy_width*south_north/)
315 map_msfu = (/south_north, 1, bdy_width*south_north/)
316 map_msfv = (/south_north_stag, 1, bdy_width*south_north_stag/)
317
318 reverse = .False.
319 tenname = "_BTXS"
320 case ("_BXE") ! East Boundary
321 start_u = (/west_east_stag - bdy_width + 1, 1, 1, 1/)
322 start_v = (/west_east - bdy_width + 1, 1, 1, 1/)
323 start_mass = (/west_east - bdy_width + 1, 1, 1, 1/)
324 start_3d = (/west_east - bdy_width + 1, 1, 1/)
325 start_msfu = (/west_east_stag - bdy_width + 1, 1, 1/)
326 start_msfv = (/west_east - bdy_width + 1, 1, 1/)
327
328 cnt_4d = (/dsizes(3),dsizes(1),dsizes(2),1/)
329 cnt_3d = (/bdy_width,south_north,1/)
330 cnt_msfu = (/bdy_width,south_north,1/)
331 cnt_msfv = (/bdy_width,south_north_stag,1/)
332
333 map_4d = (/dsizes(1)*dsizes(2), 1, dsizes(1), dsizes(1)*dsizes(2)*dsizes(3)/)
334 map_3d = (/south_north, 1, bdy_width*south_north/)
335 map_msfu = (/south_north, 1, bdy_width*south_north/)
336 map_msfv = (/south_north_stag, 1, bdy_width*south_north_stag/)
337
338 reverse = .True.
339 tenname = "_BTXE"
340 case ("_BYE") ! North Boundary
341 start_u = (/1, south_north - bdy_width + 1, 1, 1/)
342 start_v = (/1, south_north_stag - bdy_width + 1, 1, 1/)
343 start_mass = (/1, south_north - bdy_width + 1, 1, 1/)
344 start_3d = (/1, south_north - bdy_width + 1, 1/)
345 start_msfu = (/1, south_north - bdy_width + 1, 1/)
346 start_msfv = (/1, south_north_stag - bdy_width + 1, 1/)
347
348 cnt_4d = (/dsizes(1),dsizes(3),dsizes(2),1/)
349 cnt_3d = (/west_east, bdy_width,1/)
350 cnt_msfu = (/west_east_stag, bdy_width,1/)
351 cnt_msfv = (/west_east, bdy_width,1/)
352
353 map_4d = (/1, dsizes(1)*dsizes(2), dsizes(1), dsizes(3)*dsizes(1)*dsizes(2)/)
354 map_3d = (/1, west_east, west_east*bdy_width/)
355 map_msfu = (/1, west_east_stag, west_east_stag*bdy_width/)
356 map_msfv = (/1, west_east, west_east*bdy_width/)
357
358 reverse = .True.
359 tenname = "_BTYE"
360
361 case ("_BYS") ! South Boundary
362 start_u = (/1, 1, 1, 1/)
363 start_v = (/1, 1, 1, 1/)
364 start_mass = (/1, 1, 1, 1/)
365 start_3d = (/1, 1, 1/)
366 start_msfu = (/1, 1, 1/)
367 start_msfv = (/1, 1, 1/)
368
369 cnt_4d = (/dsizes(1),dsizes(3),dsizes(2),1/)
370 cnt_3d = (/west_east, bdy_width,1/)
371 cnt_msfu = (/west_east_stag, bdy_width,1/)
372 cnt_msfv = (/west_east, bdy_width,1/)
373
374 map_4d = (/1, dsizes(1)*dsizes(2), dsizes(1), dsizes(3)*dsizes(1)*dsizes(2)/)
375 map_3d = (/1, west_east, west_east*bdy_width/)
376 map_msfu = (/1, west_east_stag, west_east_stag*bdy_width/)
377 map_msfv = (/1, west_east, west_east*bdy_width/)
378
379 reverse = .False.
380 tenname = "_BTYS"
381
382 case ("_BTXS", "_BTXE","_BTYS","_BTYE")
383 cycle
384 end select
385
386 select case (nDims)
387 case (2)
388 if (vNam(1:offset) == "Times") then
389 ncid = ncidfg
390 else
391 n = index(vNam, "bdytime")
392 if ( n <= 0 ) cycle
393 select case (vNam(n-4:n-1))
394 case ("this")
395 ncid = ncidfg
396 case ("next")
397 ncid = ncidfg02
398 case default
399 cycle
400 end select
401 end if
402 status = nf90_get_var(ncid, vTimes_ID, times)
403 status = nf90_put_var(ncidvarbdy, varid, times)
404 case (3,4)
405
406 Write(*,*) "Processing for "//trim(vNam)
407
408 couple = .true.
409
410 allocate(MU_fg (dsizes(1),bdy_width,1), stat=status)
411 allocate(MU_fg02 (dsizes(1),bdy_width,1), stat=status)
412 allocate(MUB_fg (dsizes(1),bdy_width,1), stat=status)
413 allocate(MUB_fg02(dsizes(1),bdy_width,1), stat=status)
414 allocate(MSF (dsizes(1),bdy_width,1), stat=status)
415
416 allocate(Tend(dsizes(1), dsizes(2), dsizes(3), dsizes(4)), stat=status)
417
418 if ( dsizes(1) == west_east_stag .or. dsizes(1) == south_north_stag ) then
419 MU_fgptr => MU_fg (2:,:,:)
420 MU_fg02ptr => MU_fg02 (2:,:,:)
421 MUB_fgptr => MUB_fg (2:,:,:)
422 MUB_fg02ptr => MUB_fg02(2:,:,:)
423 stag = .True.
424 else
425 MU_fgptr => MU_fg
426 MU_fg02ptr => MU_fg02
427 MUB_fgptr => MUB_fg
428 MUB_fg02ptr => MUB_fg02
429 stag = .False.
430 end if
431
432 err_msg="Failed to get variable : "//trim(vNam)
433 status = nf90_get_var(ncidfg, MU_fgID, MU_fgptr, start=start_3d, count=cnt_3d, map=map_3d)
434 if ( status /= nf90_noerr ) call nf90_handle_err(status, err_msg)
435
436 status = nf90_get_var(ncidfg02, MU_fg02ID, MU_fg02ptr, start=start_3d,count=cnt_3d, map=map_3d)
437 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
438
439 status = nf90_get_var(ncidfg, MUB_fgID, MUB_fgptr, start=start_3d, count=cnt_3d,map=map_3d)
440 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
441
442 status = nf90_get_var(ncidfg02, MUB_fg02ID, MUB_fg02ptr, start=start_3d, count=cnt_3d, map=map_3d)
443 if(status /= nf90_noerr) call nf90_handle_err(status, err_msg)
444
445 err_msg="Failed to inquire tendency id for "//trim(vNam)//" for output file"
446 status = nf90_inq_varid(ncidvarbdy, vNam(1:offset-1)//tenname, tenid)
447 if(status /= nf90_noerr) call nf90_handle_err(status, err_msg)
448
449 if ( reverse ) then
450 MU_fg = MU_fg (:,bdy_width:1:-1,:)
451 MU_fg02 = MU_fg02 (:,bdy_width:1:-1,:)
452 MUB_fg = MUB_fg (:,bdy_width:1:-1,:)
453 MUB_fg02 = MUB_fg02(:,bdy_width:1:-1,:)
454 end if
455
456 select case (vNam(1:offset))
457 case ("U_", "V_")
458 if ( stag ) then
459 MU_fg (1,:,:) = MU_fg (2,:,:)
460 MU_fg02 (1,:,:) = MU_fg02 (2,:,:)
461 MUB_fg (1,:,:) = MUB_fg (2,:,:)
462 MUB_fg02(1,:,:) = MUB_fg02(2,:,:)
463
464 MU_fg (2:dsizes(1)-1,:,:) = (MU_fg (2:dsizes(1)-1,:,:) + MU_fg (3:dsizes(1),:,:))*0.5
465 MU_fg02 (2:dsizes(1)-1,:,:) = (MU_fg02 (2:dsizes(1)-1,:,:) + MU_fg02 (3:dsizes(1),:,:))*0.5
466 MUB_fg (2:dsizes(1)-1,:,:) = (MUB_fg (2:dsizes(1)-1,:,:) + MUB_fg (3:dsizes(1),:,:))*0.5
467 MUB_fg02(2:dsizes(1)-1,:,:) = (MUB_fg02(2:dsizes(1)-1,:,:) + MUB_fg02(3:dsizes(1),:,:))*0.5
468 else
469 MU_fg (:,2:bdy_width,:) = (MU_fg (:,1:bdy_width-1,:) + MU_fg (:,2:bdy_width,:))*0.5
470 MU_fg02 (:,2:bdy_width,:) = (MU_fg02 (:,1:bdy_width-1,:) + MU_fg02 (:,2:bdy_width,:))*0.5
471 MUB_fg (:,2:bdy_width,:) = (MUB_fg (:,1:bdy_width-1,:) + MUB_fg (:,2:bdy_width,:))*0.5
472 MUB_fg02(:,2:bdy_width,:) = (MUB_fg02(:,1:bdy_width-1,:) + MUB_fg02(:,2:bdy_width,:))*0.5
473 end if
474
475 if ( vNam(1:offset) == "U_" ) then
476 start_4d => start_u
477 start_msf => start_msfu
478 cnt_msf => cnt_msfu
479 map_msf => map_msfu
480 MSF_NAME = "MAPFAC_U"
481 else
482 start_4d => start_v
483 start_msf => start_msfv
484 cnt_msf => cnt_msfv
485 map_msf => map_msfv
486 MSF_NAME = "MAPFAC_V"
487 end if
488
489 status = nf90_inq_varid(ncidfg , MSF_NAME , MSF_ID )
490 err_msg="Failed to get varialbe MSF"
491 status = nf90_get_var(ncidfg, MSF_ID, MSF, start=start_msf, count=cnt_msf, map=map_msf)
492 if(status /= nf90_noerr) call nf90_handle_err(status, err_msg)
493
494 if ( reverse ) MSF = MSF(:,bdy_width:1:-1,:)
495
496 case ("T_","PH_","QVAPOR_")
497 MSF = 1.0
498 start_4d => start_mass
499 case ("MU_")
500 status = nf90_inq_varid(ncidvarbdy, "MU"//tenname, tenid)
501 Tend(:,:,:,1) = ( MU_fg02 - MU_fg ) / bdyfrq
502 status = nf90_put_var(ncidvarbdy, varid, MU_fg)
503 !status = nf90_put_var(ncidvarbdy, varid, MU_fg02)
504 err_msg="Failed to put variable "//trim(vNam)
505 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
506 status = nf90_put_var(ncidvarbdy, tenid, Tend(:,:,:,1))
507 err_msg="Failed to put tendency for "//trim(vNam)
508 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
509 couple = .false.
510
511 case default
512 Tend = 0.0
513 couple = .false.
514 select case (xtype)
515 case (nf90_float)
516 allocate(fVar_fg( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
517 fVar_fg = 0.0
518 status = nf90_put_var(ncidvarbdy, varid, fVar_fg)
519 err_msg="Failed to put variable "//trim(vNam)
520 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
521 status = nf90_put_var(ncidvarbdy, tenid, Tend)
522 err_msg="Failed to put tendency for "//trim(vNam)
523 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
524 deallocate (fVar_fg)
525 case (nf90_int)
526 allocate(iVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
527 iVar = 0
528 status = nf90_put_var(ncidvarbdy, varid, iVar)
529 err_msg="Failed to put variable "//trim(vNam)
530 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
531 status = nf90_put_var(ncidvarbdy, tenid, Tend)
532 err_msg="Failed to put tendency for "//trim(vNam)
533 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
534 deallocate (iVar)
535 end select ! end of xtype
536
537 end select ! end of vNam
538
539 if ( couple ) then
540
541 allocate( fVar_fg(dsizes(1), dsizes(2), dsizes(3), dsizes(4)), stat=status)
542 allocate(fVar_fg02(dsizes(1), dsizes(2), dsizes(3), dsizes(4)), stat=status)
543
544 err_msg="Failed to inquire variable id for "//vNam(1:offset-1)//" for fg"
545 status = nf90_inq_varid(ncidfg, vNam(1:offset-1), fg_varid)
546 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
547
548 err_msg="Failed to inquire variable id for "//vNam(1:offset-1)//" for fg02"
549 status = nf90_inq_varid(ncidfg02, vNam(1:offset-1), fg02_varid)
550 if(status /= nf90_noerr) call nf90_handle_err(status, err_msg)
551
552 err_msg="Failed to inquire tendency id for "//trim(vNam(1:offset-1))//" for output file"
553 status = nf90_inq_varid(ncidvarbdy, vNam(1:offset-1)//tenname, tenid)
554 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
555
556 err_msg="Failed to get variable "//vNam(1:offset-1)//" from fg"
557 status = nf90_get_var(ncidfg, fg_varid, fVar_fg, start=start_4d, count=cnt_4d, map=map_4d)
558 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
559
560 err_msg="Failed to get variable "//vNam(1:offset-1)//" from fg02"
561 status = nf90_get_var(ncidfg02, fg02_varid, fVar_fg02, start=start_4d, count=cnt_4d, map=map_4d)
562 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
563
564 MU_fg = MU_fg + MUB_fg
565 MU_fg02 = MU_fg02 + MUB_fg
566 !MU_fg02 = MU_fg02 + MUB_fg02
567
568 if ( reverse ) then
569 fVar_fg = fVar_fg (:,:,bdy_width:1:-1,:)
570 fVar_fg02 = fVar_fg02(:,:,bdy_width:1:-1,:)
571 end if
572
573 do i = 1, dsizes(2)
574 fVar_fg(:,i,:,:) = (fVar_fg (:,i,:,:) * MU_fg ) / MSF
575 fVar_fg02(:,i,:,:) = (fVar_fg02(:,i,:,:) * MU_fg02) / MSF
576 end do
577
578 Tend = ( fVar_fg02 - fVar_fg ) / bdyfrq
579
580 err_msg="Failed to put variable "//trim(vNam)
581 status = nf90_put_var(ncidvarbdy, varid, fVar_fg)
582 !status = nf90_put_var(ncidvarbdy, varid, fVar_fg02)
583 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
584
585 err_msg="Failed to put tendency for "//trim(vNam)
586 status = nf90_put_var(ncidvarbdy, tenid, Tend)
587 if(status /= nf90_noerr) call nf90_handle_err(status,err_msg)
588
589 deallocate (fVar_fg)
590 deallocate (fVar_fg02)
591
592 end if
593
594 NULLIFY (MU_fgptr)
595 NULLIFY (MU_fg02ptr)
596 NULLIFY (MUB_fgptr)
597 NULLIFY (MUB_fg02ptr)
598 NULLIFY (MSF_ptr)
599
600 deallocate (Tend)
601 deallocate (MU_fg)
602 deallocate (MU_fg02)
603 deallocate (MUB_fg)
604 deallocate (MUB_fg02)
605 deallocate (MSF)
606 case default
607 cycle
608 end select ! end of nDims
609
610 end do
611
612 deallocate (times)
613
614 status = nf90_close(ncidfg)
615 status = nf90_close(ncidfg02)
616 status = nf90_close(ncidwrfbdy)
617 status = nf90_close(ncidvarbdy)
618
619 Write(*,*) "Boundary file generated successfully"
620
621 contains
622
623 subroutine nf90_handle_err(status, err_msg)
624 integer, intent (in) :: status
625 character (len=*), intent(in) :: err_msg
626
627 if(status /= nf90_noerr) then
628 print *, trim(nf90_strerror(status))
629 print *, trim(err_msg)
630 call exit(-1)
631 end if
632 end subroutine nf90_handle_err
633
634 function jd(yyyy, mm, dd) result(ival)
635
636 integer, intent(in) :: yyyy
637 integer, intent(in) :: mm
638 integer, intent(in) :: dd
639 integer :: ival
640
641 ! DATE ROUTINE JD(YYYY, MM, DD) CONVERTS CALENDER DATE TO
642 ! JULIAN DATE. SEE CACM 1968 11(10):657, LETTER TO THE
643 ! EDITOR BY HENRY F. FLIEGEL AND THOMAS C. VAN FLANDERN.
644 ! EXAMPLE JD(1970, 1, 1) = 2440588
645
646 ival = dd - 32075 + 1461*(yyyy+4800+(mm-14)/12)/4 + &
647 367*(mm-2-((mm-14)/12)*12)/12 - 3*((yyyy+4900+(mm-14)/12)/100)/4
648
649 return
650 end function jd
651
652 function datediff(date_1, date_2) result(ival)
653
654 character(len=*), intent(in) :: date_1
655 character(len=*), intent(in) :: date_2
656 integer :: ival
657 integer :: jd1, jd2
658 integer :: yyyy,mm,dd
659 integer :: hh1,nn1,ss1
660 integer :: hh2,nn2,ss2
661
662
663 ! date string : yyyy-mm-dd_hh:mm:ss
664 ! calculate the difference between date_1 and date_2 in seconds
665
666 read(date_1(1:19), '(i4,5(1x,i2))') &
667 yyyy, mm, dd, hh1, nn1, ss1
668
669 jd1=jd(yyyy,mm,dd)
670
671 read(date_2(1:19), '(i4,5(1x,i2))') &
672 yyyy, mm, dd, hh2, nn2, ss2
673
674 jd2=jd(yyyy,mm,dd)
675
676 ival=(jd2-jd1)*86400 + ( hh2-hh1)*3600 + (nn2-nn1)*60 + (ss2-ss1)
677
678 return
679 end function datediff
680
681 end program da_bdy
Mirror of the official repository at https://github.com/wrf-model/WRF