emcsfc_ice_blend.f90

Go to the documentation of this file.

 program emcsfc_ice_blend

 use grib_mod  ! grib 2 libraries

 implicit none

 type(gribfield)        :: ims
 type(gribfield)        :: mask
 type(gribfield)        :: mmab
 character(len=200)     :: infile, outfile
 integer, parameter     :: imax=4320
 integer, parameter     :: jmax=2160
 integer                :: i,j, istat, iunit
 integer                :: ii, iii, jj, jjj, count
 integer                :: lugi
 integer                :: jdisc
 integer                :: jgdtn
 integer                :: jpdtn
 integer                :: k
 integer                :: jids(200)
 integer                :: jgdt(200)
 integer                :: jpdt(200)
 integer, allocatable   :: mask_5min(:,:), mask_ims(:,:)

 logical*1, allocatable :: lbms_ims(:,:)
 logical                :: unpack

 real, allocatable      :: dummy(:,:)
 real, allocatable      :: ice_ims(:,:), ice_5min(:,:), ice_blend(:,:)

!--------------------------------------------------------------------------------
! Read the 5-minute mmab land mask file.  Required because the 5-minute data does
! not have a bitmap.
!--------------------------------------------------------------------------------

 call w3tagb('EMCSFC_ICE_BLEND',2014,75,0000,'EMC')

 call getenv("FORT17", infile)
 iunit=17
 print*,"- OPEN 5-MINUTE LAND-SEA MASK FILE: ", trim(infile)
 call baopenr (iunit, infile, istat)
 if (istat /= 0) then
   print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
   stop 1
 endif
 
 nullify(mask%idsect)
 nullify(mask%local)
 nullify(mask%list_opt)
 nullify(mask%igdtmpl)
 nullify(mask%ipdtmpl)
 nullify(mask%coord_list)
 nullify(mask%idrtmpl)
 nullify(mask%bmap)
 nullify(mask%fld)

 j       = 0      ! search at beginning of file
 lugi    = 0      ! no grib index file
 jdisc   = 2      ! search for discipline
 jpdtn   = 0      ! search for product definition template number
 jgdtn   = 0      ! search for grid definition template number; 0 - lat/lon grid
 jids    = -9999  ! array of values in identification section, set to wildcard
 jgdt    = -9999  ! array of values in grid definition template 3.m
 jpdt    = -9999  ! array of values in product definition template 4.n
 jpdt(1) = 0      ! search for parameter category 
 jpdt(2) = 0      ! search for parameter number 
 unpack  = .true. ! unpack data

 print*,"- DEGRIB DATA"
 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
             unpack, k, mask, istat)

 if (istat /= 0) then
   print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
   stop 2
 endif

 call baclose (iunit,istat)

!--------------------------------------------------------------------------------
!  The MMAB mask codes are:
!
!  0    - water
!  1.57 - land
!  1.95 - coast
!
!  Convert these to something easier to work with.
!--------------------------------------------------------------------------------

 allocate(dummy(imax,jmax))
 dummy=reshape(mask%fld , (/imax,jmax/) )

 if(associated(mask%idsect))     deallocate(mask%idsect)
 if(associated(mask%local))      deallocate(mask%local)
 if(associated(mask%list_opt))   deallocate(mask%list_opt)
 if(associated(mask%igdtmpl))    deallocate(mask%igdtmpl)
 if(associated(mask%ipdtmpl))    deallocate(mask%ipdtmpl)
 if(associated(mask%coord_list)) deallocate(mask%coord_list)
 if(associated(mask%idrtmpl))    deallocate(mask%idrtmpl)
 if(associated(mask%bmap))       deallocate(mask%bmap)
 if(associated(mask%fld))        deallocate(mask%fld)

 allocate(mask_5min(imax,jmax))

 do j = 1, jmax
 do i = 1, imax
   if (dummy(i,j) < 0.1) then
      mask_5min(i,j)=0   ! water
   elseif (dummy(i,j) > 1.94) then
      mask_5min(i,j)=1   ! coast
   else
      mask_5min(i,j)=2   ! land
   endif
 enddo
 enddo

 deallocate(dummy)

!--------------------------------------------------------------------------------
! Read ims data that has been interpolated to the 5-min grid.
!--------------------------------------------------------------------------------

 call getenv("FORT11", infile)
 iunit=11
 print*,"- OPEN IMS ICE DATA: ", trim(infile)
 call baopenr (iunit, infile, istat)
 if (istat /= 0) then
   print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
   stop 5
 endif
 
 nullify(ims%idsect)
 nullify(ims%local)
 nullify(ims%list_opt)
 nullify(ims%igdtmpl)
 nullify(ims%ipdtmpl)
 nullify(ims%coord_list)
 nullify(ims%idrtmpl)
 nullify(ims%bmap)
 nullify(ims%fld)

 j       = 0      ! search at beginning of file
 lugi    = 0      ! no grib index file
 jdisc   = 10     ! search for discipline, ocean products
 jpdtn   = 0      ! search for product definition template number
 jgdtn   = 0      ! search for grid definition template number; 0 - lat/lon grid
 jids    = -9999  ! array of values in identification section, set to wildcard
 jgdt    = -9999  ! array of values in grid definition template 3.m
 jpdt    = -9999  ! array of values in product definition template 4.n
 jpdt(1) = 2      ! search for parameter category, ice
 jpdt(2) = 0      ! search for parameter number, ice cover
 unpack  = .true. ! unpack data

 print*,"- DEGRIB DATA"
 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
             unpack, k, ims, istat)

 if (istat /= 0) then
   print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
   stop 7
 endif

 call baclose (iunit,istat)

 allocate (lbms_ims(imax,jmax))
 lbms_ims=reshape(ims%bmap , (/imax,jmax/) )
 allocate (ice_ims(imax,jmax))
 ice_ims=reshape(ims%fld , (/imax,jmax/) )

 print*,"- CREATE IMS LAND-SEA MASK FROM BITMAP."
 allocate(mask_ims(imax,jmax))
 mask_ims = 0  ! water
 do j = 1, jmax
 do i = 1, imax
   if (.not.lbms_ims(i,j)) mask_ims(i,j) = 2  ! land
 enddo
 enddo

 deallocate(lbms_ims)

 if(associated(ims%idsect))     deallocate(ims%idsect)
 if(associated(ims%local))      deallocate(ims%local)
 if(associated(ims%list_opt))   deallocate(ims%list_opt)
 if(associated(ims%igdtmpl))    deallocate(ims%igdtmpl)
 if(associated(ims%ipdtmpl))    deallocate(ims%ipdtmpl)
 if(associated(ims%coord_list)) deallocate(ims%coord_list)
 if(associated(ims%idrtmpl))    deallocate(ims%idrtmpl)
 if(associated(ims%bmap))       deallocate(ims%bmap)
 if(associated(ims%fld))        deallocate(ims%fld)

!--------------------------------------------------------------------------------
! Now read the MMAB 5-minute data.
!--------------------------------------------------------------------------------

 call getenv("FORT15", infile)
 iunit=15
 print*,"- OPEN 5-MINUTE ICE CONCENTRATION DATA: ", trim(infile)
 call baopenr (iunit, infile, istat)
 if (istat /= 0) then
   print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
   stop 8
 endif
 
 nullify(mmab%idsect)
 nullify(mmab%local)
 nullify(mmab%list_opt)
 nullify(mmab%igdtmpl)
 nullify(mmab%ipdtmpl)
 nullify(mmab%coord_list)
 nullify(mmab%idrtmpl)
 nullify(mmab%bmap)
 nullify(mmab%fld)

 j       = 0      ! search at beginning of file
 lugi    = 0      ! no grib index file
 jdisc   = 10     ! search for discipline, ocean products
 jpdtn   = 0      ! search for product definition template number
 jgdtn   = 0      ! search for grid definition template number; 0 - lat/lon grid
 jids    = -9999  ! array of values in identification section, set to wildcard
 jgdt    = -9999  ! array of values in grid definition template 3.m
 jpdt    = -9999  ! array of values in product definition template 4.n
 jpdt(1) = 2      ! search for parameter category, ice
 jpdt(2) = 0      ! search for parameter number, ice cover
 unpack  = .true. ! unpack data

 print*,"- DEGRIB DATA"
 call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
             unpack, k, mmab, istat)

 if (istat /= 0) then
   print*,'FATAL ERROR: BAD DEGRIB OF DATA. ISTAT: ', istat
   stop 9
 endif

 allocate (ice_5min(imax,jmax))
 ice_5min=reshape(mmab%fld , (/imax,jmax/) )

 call baclose (iunit,istat)
 
!--------------------------------------------------------------------------------
! Blend IMS and 5-minute data in northern hemisphere.
!--------------------------------------------------------------------------------

 print*,"- BLEND IMS AND 5-MINUTE DATA IN NH."
 allocate(ice_blend(imax,jmax))
 ice_blend=-9.   ! flag for land.  will be bitmapped out.
 do j = 1, (jmax/2)
 do i = 1, imax
   if (mask_ims(i,j) == 0) then   ! ims water point
     if (mask_5min(i,j) > 0) then ! 5-min land
       ice_blend(i,j)=ice_ims(i,j)   ! use ims value
     else  ! ims and 5min mask indicate water point
       if (ice_ims(i,j) > .5) then   ! ims indicates ice
         count = 0
         do jj = -1, 1
         do ii = -1, 1
           if (ii == 0 .and. jj == 0) cycle
           jjj = j + jj
           if (jjj < 1) cycle
           iii = ii + i
           if (iii < 1) iii = iii + imax
           if (iii > imax) iii = iii - imax
           if (mask_5min(iii,jjj) == 0) then ! 5-min water
             if (ice_5min(iii,jjj) >= 0.5) then
               count = count + 1
             endif
           endif
         enddo
         enddo
         if (count > 0 .and. ice_5min(i,j) == 0.0) then
           ice_blend(i,j) = ice_ims(i,j)
         else
           ice_blend(i,j) = max(ice_5min(i,j),0.15)
         endif
       else  ! ims indicates open water.
         ice_blend(i,j) = 0.
       endif
     endif
   end if
 enddo
 enddo

 deallocate(mask_ims, ice_ims)

!--------------------------------------------------------------------------------
! In the SH, the blend is simply the 5-minute data.  Only consider 'water'
! points.
!--------------------------------------------------------------------------------

 do j = (jmax/2)+1, jmax
 do i = 1, imax
   if (mask_5min(i,j) == 0) then  ! 'water'
     ice_blend(i,j) = ice_5min(i,j)
   endif
 enddo
 enddo

 deallocate(mask_5min, ice_5min)

!--------------------------------------------------------------------------------
! Output blended data to a grib 2 file.
!--------------------------------------------------------------------------------

 call getenv("FORT51", outfile)
 iunit=51
 print*,"- OUTPUT BLENDED ICE DATA TO ", trim(outfile)
 print*,"- OPEN FILE."
 call baopenw(iunit, outfile, istat)
 if (istat /= 0) then
   print*,'FATAL ERROR: BAD OPEN. ISTAT: ', istat
   stop 16
 endif

!--------------------------------------------------------------------------------
! Use grib header information from the mmab file (stored in the mmab data
! structure) with the following exceptions:
!
! 1) Increase precision of corner point lat/lons and dx/dy.
! 2) Use simple packing instead of jpeg compression.  Grads has problems
!    with the latter.
! 3) Use a bitmap to mask out land points instead of the mmab convention
!    of using '0' at land points.
!--------------------------------------------------------------------------------

 mmab%igdtmpl(12)=89958333
 mmab%igdtmpl(13)=41667
 mmab%igdtmpl(15)=-89958333
 mmab%igdtmpl(16)=359958333
 mmab%igdtmpl(17)=83333
 mmab%igdtmpl(18)=83333

 deallocate (mmab%idrtmpl)   ! use simple packing
 mmab%idrtnum=0  
 mmab%idrtlen=5
 allocate(mmab%idrtmpl(mmab%idrtlen))
 mmab%idrtmpl=0
 mmab%idrtmpl(3) = 2         ! use two decimal points.

 mmab%fld=reshape(ice_blend, (/imax*jmax/) )

 mmab%ibmap=0  ! use bitmap
 allocate (mmab%bmap(imax*jmax))
 mmab%bmap=.true.
 where (mmab%fld < -8.) mmab%bmap=.false.   ! land

 print*,"- GRIB DATA."
 call putgb2(iunit, mmab, istat)
 if (istat /= 0) then
   print*,'FATAL ERROR: BAD WRITE. ISTAT: ', istat
   stop 17
 endif

 call baclose(iunit, istat)

 deallocate(ice_blend)

 if(associated(mmab%idsect))     deallocate(mmab%idsect)
 if(associated(mmab%local))      deallocate(mmab%local)
 if(associated(mmab%list_opt))   deallocate(mmab%list_opt)
 if(associated(mmab%igdtmpl))    deallocate(mmab%igdtmpl)
 if(associated(mmab%ipdtmpl))    deallocate(mmab%ipdtmpl)
 if(associated(mmab%coord_list)) deallocate(mmab%coord_list)
 if(associated(mmab%idrtmpl))    deallocate(mmab%idrtmpl)
 if(associated(mmab%bmap))       deallocate(mmab%bmap)
 if(associated(mmab%fld))        deallocate(mmab%fld)

 print*,''
 print*,'****************************'
 print*,'**** NORMAL TERMINATION ****'
 print*,'****************************'

 call w3tage('EMCSFC_ICE_BLEND')

 stop

 end program emcsfc_ice_blend