emcsfc_ice_blend  1.5.0
 All Files Functions
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