emcsfc_ice_blend 1.14.0
Loading...
Searching...
No Matches
emcsfc_ice_blend.f90
Go to the documentation of this file.
1
6
93 program emcsfc_ice_blend
94
95 use grib_mod ! grib 2 libraries
96
97 implicit none
98
99 type(gribfield) :: ims
100 type(gribfield) :: mask
101 type(gribfield) :: mmab
102 character(len=200) :: infile, outfile
103 integer, parameter :: imax=4320
104 integer, parameter :: jmax=2160
105 integer :: i,j, istat, iunit
106 integer :: ii, iii, jj, jjj, count
107 integer :: lugi
108 integer :: jdisc
109 integer :: jgdtn
110 integer :: jpdtn
111 integer :: k
112 integer :: jids(200)
113 integer :: jgdt(200)
114 integer :: jpdt(200)
115 integer, allocatable :: mask_5min(:,:), mask_ims(:,:)
116
117 logical*1, allocatable :: lbms_ims(:,:)
118 logical :: unpack
119
120 real, allocatable :: dummy(:,:)
121 real, allocatable :: ice_ims(:,:), ice_5min(:,:), ice_blend(:,:)
122
123!--------------------------------------------------------------------------------
124! Read the 5-minute mmab land mask file. Required because the 5-minute data does
125! not have a bitmap.
126!--------------------------------------------------------------------------------
127
128 call w3tagb('EMCSFC_ICE_BLEND',2014,75,0000,'EMC')
129
130 call getenv("FORT17", infile)
131 iunit=17
132 print*,"- OPEN 5-MINUTE LAND-SEA MASK FILE: ", trim(infile)
133 call baopenr (iunit, infile, istat)
134 if (istat /= 0) then
135 print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
136 stop 1
137 endif
138
139 nullify(mask%idsect)
140 nullify(mask%local)
141 nullify(mask%list_opt)
142 nullify(mask%igdtmpl)
143 nullify(mask%ipdtmpl)
144 nullify(mask%coord_list)
145 nullify(mask%idrtmpl)
146 nullify(mask%bmap)
147 nullify(mask%fld)
148
149 j = 0 ! search at beginning of file
150 lugi = 0 ! no grib index file
151 jdisc = 2 ! search for discipline
152 jpdtn = 0 ! search for product definition template number
153 jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid
154 jids = -9999 ! array of values in identification section, set to wildcard
155 jgdt = -9999 ! array of values in grid definition template 3.m
156 jpdt = -9999 ! array of values in product definition template 4.n
157 jpdt(1) = 0 ! search for parameter category
158 jpdt(2) = 0 ! search for parameter number
159 unpack = .true. ! unpack data
160
161 print*,"- DEGRIB DATA"
162 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
163 unpack, k, mask, istat)
164
165 if (istat /= 0) then
166 print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
167 stop 2
168 endif
169
170 call baclose (iunit,istat)
171
172!--------------------------------------------------------------------------------
173! The MMAB mask codes are:
174!
175! 0 - water
176! 1.57 - land
177! 1.95 - coast
178!
179! Convert these to something easier to work with.
180!--------------------------------------------------------------------------------
181
182 allocate(dummy(imax,jmax))
183 dummy=reshape(mask%fld , (/imax,jmax/) )
184
185 if(associated(mask%idsect)) deallocate(mask%idsect)
186 if(associated(mask%local)) deallocate(mask%local)
187 if(associated(mask%list_opt)) deallocate(mask%list_opt)
188 if(associated(mask%igdtmpl)) deallocate(mask%igdtmpl)
189 if(associated(mask%ipdtmpl)) deallocate(mask%ipdtmpl)
190 if(associated(mask%coord_list)) deallocate(mask%coord_list)
191 if(associated(mask%idrtmpl)) deallocate(mask%idrtmpl)
192 if(associated(mask%bmap)) deallocate(mask%bmap)
193 if(associated(mask%fld)) deallocate(mask%fld)
194
195 allocate(mask_5min(imax,jmax))
196
197 do j = 1, jmax
198 do i = 1, imax
199 if (dummy(i,j) < 0.1) then
200 mask_5min(i,j)=0 ! water
201 elseif (dummy(i,j) > 1.94) then
202 mask_5min(i,j)=1 ! coast
203 else
204 mask_5min(i,j)=2 ! land
205 endif
206 enddo
207 enddo
208
209 deallocate(dummy)
210
211!--------------------------------------------------------------------------------
212! Read ims data that has been interpolated to the 5-min grid.
213!--------------------------------------------------------------------------------
214
215 call getenv("FORT11", infile)
216 iunit=11
217 print*,"- OPEN IMS ICE DATA: ", trim(infile)
218 call baopenr (iunit, infile, istat)
219 if (istat /= 0) then
220 print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
221 stop 5
222 endif
223
224 nullify(ims%idsect)
225 nullify(ims%local)
226 nullify(ims%list_opt)
227 nullify(ims%igdtmpl)
228 nullify(ims%ipdtmpl)
229 nullify(ims%coord_list)
230 nullify(ims%idrtmpl)
231 nullify(ims%bmap)
232 nullify(ims%fld)
233
234 j = 0 ! search at beginning of file
235 lugi = 0 ! no grib index file
236 jdisc = 10 ! search for discipline, ocean products
237 jpdtn = 0 ! search for product definition template number
238 jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid
239 jids = -9999 ! array of values in identification section, set to wildcard
240 jgdt = -9999 ! array of values in grid definition template 3.m
241 jpdt = -9999 ! array of values in product definition template 4.n
242 jpdt(1) = 2 ! search for parameter category, ice
243 jpdt(2) = 0 ! search for parameter number, ice cover
244 unpack = .true. ! unpack data
245
246 print*,"- DEGRIB DATA"
247 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
248 unpack, k, ims, istat)
249
250 if (istat /= 0) then
251 print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
252 stop 7
253 endif
254
255 call baclose (iunit,istat)
256
257 allocate (lbms_ims(imax,jmax))
258 lbms_ims=reshape(ims%bmap , (/imax,jmax/) )
259 allocate (ice_ims(imax,jmax))
260 ice_ims=reshape(ims%fld , (/imax,jmax/) )
261
262 print*,"- CREATE IMS LAND-SEA MASK FROM BITMAP."
263 allocate(mask_ims(imax,jmax))
264 mask_ims = 0 ! water
265 do j = 1, jmax
266 do i = 1, imax
267 if (.not.lbms_ims(i,j)) mask_ims(i,j) = 2 ! land
268 enddo
269 enddo
270
271 deallocate(lbms_ims)
272
273 if(associated(ims%idsect)) deallocate(ims%idsect)
274 if(associated(ims%local)) deallocate(ims%local)
275 if(associated(ims%list_opt)) deallocate(ims%list_opt)
276 if(associated(ims%igdtmpl)) deallocate(ims%igdtmpl)
277 if(associated(ims%ipdtmpl)) deallocate(ims%ipdtmpl)
278 if(associated(ims%coord_list)) deallocate(ims%coord_list)
279 if(associated(ims%idrtmpl)) deallocate(ims%idrtmpl)
280 if(associated(ims%bmap)) deallocate(ims%bmap)
281 if(associated(ims%fld)) deallocate(ims%fld)
282
283!--------------------------------------------------------------------------------
284! Now read the MMAB 5-minute data.
285!--------------------------------------------------------------------------------
286
287 call getenv("FORT15", infile)
288 iunit=15
289 print*,"- OPEN 5-MINUTE ICE CONCENTRATION DATA: ", trim(infile)
290 call baopenr (iunit, infile, istat)
291 if (istat /= 0) then
292 print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
293 stop 8
294 endif
295
296 nullify(mmab%idsect)
297 nullify(mmab%local)
298 nullify(mmab%list_opt)
299 nullify(mmab%igdtmpl)
300 nullify(mmab%ipdtmpl)
301 nullify(mmab%coord_list)
302 nullify(mmab%idrtmpl)
303 nullify(mmab%bmap)
304 nullify(mmab%fld)
305
306 j = 0 ! search at beginning of file
307 lugi = 0 ! no grib index file
308 jdisc = 10 ! search for discipline, ocean products
309 jpdtn = 0 ! search for product definition template number
310 jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid
311 jids = -9999 ! array of values in identification section, set to wildcard
312 jgdt = -9999 ! array of values in grid definition template 3.m
313 jpdt = -9999 ! array of values in product definition template 4.n
314 jpdt(1) = 2 ! search for parameter category, ice
315 jpdt(2) = 0 ! search for parameter number, ice cover
316 unpack = .true. ! unpack data
317
318 print*,"- DEGRIB DATA"
319 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
320 unpack, k, mmab, istat)
321
322 if (istat /= 0) then
323 print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
324 stop 9
325 endif
326
327 allocate (ice_5min(imax,jmax))
328 ice_5min=reshape(mmab%fld , (/imax,jmax/) )
329
330 call baclose (iunit,istat)
331
332!--------------------------------------------------------------------------------
333! Blend IMS and 5-minute data in northern hemisphere.
334!--------------------------------------------------------------------------------
335
336 print*,"- BLEND IMS AND 5-MINUTE DATA IN NH."
337 allocate(ice_blend(imax,jmax))
338 ice_blend=-9. ! flag for land. will be bitmapped out.
339 do j = 1, (jmax/2)
340 do i = 1, imax
341 if (mask_ims(i,j) == 0) then ! ims water point
342 if (mask_5min(i,j) > 0) then ! 5-min land
343 ice_blend(i,j)=ice_ims(i,j) ! use ims value
344 else ! ims and 5min mask indicate water point
345 if (ice_ims(i,j) > .5) then ! ims indicates ice
346 count = 0
347 do jj = -1, 1
348 do ii = -1, 1
349 if (ii == 0 .and. jj == 0) cycle
350 jjj = j + jj
351 if (jjj < 1) cycle
352 iii = ii + i
353 if (iii < 1) iii = iii + imax
354 if (iii > imax) iii = iii - imax
355 if (mask_5min(iii,jjj) == 0) then ! 5-min water
356 if (ice_5min(iii,jjj) >= 0.5) then
357 count = count + 1
358 endif
359 endif
360 enddo
361 enddo
362 if (count > 0 .and. ice_5min(i,j) == 0.0) then
363 ice_blend(i,j) = ice_ims(i,j)
364 else
365 ice_blend(i,j) = max(ice_5min(i,j),0.15)
366 endif
367 else ! ims indicates open water.
368 ice_blend(i,j) = 0.
369 endif
370 endif
371 end if
372 enddo
373 enddo
374
375 deallocate(mask_ims, ice_ims)
376
377!--------------------------------------------------------------------------------
378! In the SH, the blend is simply the 5-minute data. Only consider 'water'
379! points.
380!--------------------------------------------------------------------------------
381
382 do j = (jmax/2)+1, jmax
383 do i = 1, imax
384 if (mask_5min(i,j) == 0) then ! 'water'
385 ice_blend(i,j) = ice_5min(i,j)
386 endif
387 enddo
388 enddo
389
390 deallocate(mask_5min, ice_5min)
391
392!--------------------------------------------------------------------------------
393! Output blended data to a grib 2 file.
394!--------------------------------------------------------------------------------
395
396 call getenv("FORT51", outfile)
397 iunit=51
398 print*,"- OUTPUT BLENDED ICE DATA TO ", trim(outfile)
399 print*,"- OPEN FILE."
400 call baopenw(iunit, outfile, istat)
401 if (istat /= 0) then
402 print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
403 stop 16
404 endif
405
406!--------------------------------------------------------------------------------
407! Use grib header information from the mmab file (stored in the mmab data
408! structure) with the following exceptions:
409!
410! 1) Increase precision of corner point lat/lons and dx/dy.
411! 2) Use simple packing instead of jpeg compression. Grads has problems
412! with the latter.
413! 3) Use a bitmap to mask out land points instead of the mmab convention
414! of using '0' at land points.
415!--------------------------------------------------------------------------------
416
417 mmab%igdtmpl(12)=89958333
418 mmab%igdtmpl(13)=41667
419 mmab%igdtmpl(15)=-89958333
420 mmab%igdtmpl(16)=359958333
421 mmab%igdtmpl(17)=83333
422 mmab%igdtmpl(18)=83333
423
424 deallocate (mmab%idrtmpl) ! use simple packing
425 mmab%idrtnum=0
426 mmab%idrtlen=5
427 allocate(mmab%idrtmpl(mmab%idrtlen))
428 mmab%idrtmpl=0
429 mmab%idrtmpl(3) = 2 ! use two decimal points.
430
431 mmab%fld=reshape(ice_blend, (/imax*jmax/) )
432
433 mmab%ibmap=0 ! use bitmap
434 allocate (mmab%bmap(imax*jmax))
435 mmab%bmap=.true.
436 where (mmab%fld < -8.) mmab%bmap=.false. ! land
437
438 print*,"- GRIB DATA."
439 call putgb2(iunit, mmab, istat)
440 if (istat /= 0) then
441 print*,'FATAL ERROR: BAD WRITE. ISTAT: ', istat
442 stop 17
443 endif
444
445 call baclose(iunit, istat)
446
447 deallocate(ice_blend)
448
449 if(associated(mmab%idsect)) deallocate(mmab%idsect)
450 if(associated(mmab%local)) deallocate(mmab%local)
451 if(associated(mmab%list_opt)) deallocate(mmab%list_opt)
452 if(associated(mmab%igdtmpl)) deallocate(mmab%igdtmpl)
453 if(associated(mmab%ipdtmpl)) deallocate(mmab%ipdtmpl)
454 if(associated(mmab%coord_list)) deallocate(mmab%coord_list)
455 if(associated(mmab%idrtmpl)) deallocate(mmab%idrtmpl)
456 if(associated(mmab%bmap)) deallocate(mmab%bmap)
457 if(associated(mmab%fld)) deallocate(mmab%fld)
458
459 print*,''
460 print*,'****************************'
461 print*,'**** NORMAL TERMINATION ****'
462 print*,'****************************'
463
464 call w3tage('EMCSFC_ICE_BLEND')
465
466 stop
467
468 end program emcsfc_ice_blend
emcsfc_ice_blend
program emcsfc_ice_blend
Create a global 5-minute blended ice concentration dataset for use by GDAS/GFS.
Definition emcsfc_ice_blend.f90:93