Add the capability for assimilating GOES-ABI radiance data

Registry/registry.var

@@ -202,6 +202,7 @@ rconfig   logical   use_amsr2obs            namelist,wrfvar4  1  .false.  - "use
 rconfig   logical   use_ahiobs              namelist,wrfvar4  1  .false.  - "use_ahiobs"               ""  ""
 rconfig   logical   use_gmiobs              namelist,wrfvar4  1  .false.  - "use_gmiobs"               ""  ""
 rconfig   logical   use_goesimgobs          namelist,wrfvar4  1  .false.  - "use_goesimgobs"           ""  ""
+rconfig   logical   use_goesabiobs          namelist,wrfvar4  1  .false.  - "use_goesabiobs"           ""  ""
 rconfig   logical   use_kma1dvar            namelist,wrfvar4  1  .false.  - "use_kma1dvar"             ""  ""
 rconfig   logical   use_filtered_rad        namelist,wrfvar4  1  .false.  - "use_filtered_rad"         ""  ""
 rconfig   logical   use_obs_errfac          namelist,wrfvar4  1  .false.  - "use_obs_errfac"           ""  ""
@@ -468,6 +469,7 @@ rconfig   integer   varbc_nobsmin           namelist,wrfvar14  1  10       - "va
 rconfig   integer   use_clddet              namelist,wrfvar14  1  2        - "use_clddet"        "0: off, 1: mmr, 2: pf, 3: ecmwf"  ""
 rconfig   logical   use_clddet_zz           namelist,wrfvar14  1  .false.  - "use_clddet_zz"     "cloud detection scheme from Zhuge X. and Zou X. JAMC, 2016."  ""
 rconfig   integer   ahi_superob_halfwidth   namelist,wrfvar14  1  0        - "ahi_superob_halfwidth"   ""  ""
+rconfig   integer   abi_superob_halfwidth   namelist,wrfvar14  1  0        - "abi_superob_halfwidth"   ""  ""
 rconfig   logical   airs_warmest_fov        namelist,wrfvar14  1  .false.  - "airs_warmest_fov"        ""  ""
 rconfig   logical   use_satcv               namelist,wrfvar14  2  .false.  - "use_satcv"               ""  ""
 rconfig   logical   use_blacklist_rad       namelist,wrfvar14  1  .true.   - "use_blacklist_rad"       ""  ""
@@ -477,6 +479,7 @@ rconfig   character crtm_irwater_coef       namelist,wrfvar14  1  "Nalli.IRwater
 rconfig   character crtm_mwwater_coef       namelist,wrfvar14  1  "FASTEM5.MWwater.EmisCoeff.bin" - "crtm_mwwater_coef" "" ""
 rconfig   character crtm_irland_coef        namelist,wrfvar14  1  "USGS.IRland.EmisCoeff.bin"     - "crtm_irland_coef"  "" ""
 rconfig   character crtm_visland_coef       namelist,wrfvar14  1  "USGS.VISland.EmisCoeff.bin"    - "crtm_visland_coef" "" ""
+rconfig   logical   abi_use_symm_obs_err    namelist,wrfvar14  1  .false.  - "abi_use_symm_obs_err"         ""  ""
 rconfig   logical   ahi_use_symm_obs_err    namelist,wrfvar14  1  .false.  - "ahi_use_symm_obs_err"         ""  ""
 rconfig   logical   ahi_apply_clrsky_bias   namelist,wrfvar14  1  .false.  - "ahi_apply_clrsky_bias"         ""  ""
 rconfig   integer   num_pseudo              namelist,wrfvar15  1  0        - "num_pseudo"              ""  ""   

var/da/da_define_structures/da_define_structures.f90

@@ -574,10 +574,12 @@ module da_define_structures
       real,    pointer     :: vtox(:,:)
    end type varbc_type
    type clddet_geoir_type
-     real  :: RTCT, RFMFT, TEMPIR, terr_hgt
-     real  :: tb_stddev_10, tb_stddev_13,tb_stddev_14
-     real  :: CIRH2O
-     !real, allocatable :: CIRH2O(:,:,:)
+     real  :: RTCT, RFMFT, TEMPIR, terr_hgt ! for both ABI and AHI
+     real  :: tb_stddev_10, tb_stddev_13,tb_stddev_14 ! only for AHI
+     real  :: CIRH2O ! for both ABI and AHI
+     real, allocatable :: CIRH2O_abi(:,:,:) ! only for ABI
+     real, allocatable :: tb_stddev_3x3(:)  ! only for ABI
+     integer :: RFMFT_ij(2) ! only for ABI
    end type clddet_geoir_type   
    type superob_type
      real, allocatable :: tb_obs(:,:)
@@ -618,6 +620,8 @@ module da_define_structures
       integer, pointer     :: cloud_flag(:,:)
       integer, pointer     :: cloudflag(:)
       integer, pointer     :: rain_flag(:)
+      real,    pointer     :: cloud_mod(:,:) ! only for ABI
+      real,    pointer     :: cloud_obs(:,:) ! only for ABI
       real, allocatable    :: cloud_frac(:)
       real,    pointer     :: satzen(:) 
       real,    pointer     :: satazi(:) 
@@ -632,10 +636,10 @@ module da_define_structures
       real,    pointer     :: lod(:,:,:)       ! layer_optical_depth
       real,    pointer     :: trans(:,:,:)     ! layer transmittance
       real,    pointer     :: der_trans(:,:,:) ! d(transmittance)/dp
-      real,    pointer     :: kmin_t(:)	
-      real,    pointer     :: kmax_p(:)	  
-      real,    pointer     :: sensitivity_ratio(:,:,:)	  
-      real,    pointer     :: p_chan_level(:,:)	  
+      real,    pointer     :: kmin_t(:)
+      real,    pointer     :: kmax_p(:)
+      real,    pointer     :: sensitivity_ratio(:,:,:)
+      real,    pointer     :: p_chan_level(:,:)
       real,    pointer     :: qrn(:,:)
       real,    pointer     :: qcw(:,:)
       real,    pointer     :: qci(:,:)

var/da/da_monitor/da_rad_diags.f90

@@ -42,7 +42,7 @@ program da_rad_diags
    integer                                :: ncid, dimid, varid
    integer, dimension(3)                  :: ishape, istart, icount
 !
-   logical                                :: amsr2
+   logical                                :: amsr2, abi
    logical                                :: isfile, prf_found, jac_found
    integer, parameter                     :: datelen1 = 10
    integer, parameter                     :: datelen2 = 19
@@ -62,9 +62,9 @@ program da_rad_diags
    real*4,  dimension(:), allocatable     :: smois, tslb, snowh, vegfra, clwp, cloud_frac
    real*4,  dimension(:), allocatable     :: cip ! cloud-ice path
    integer, dimension(:), allocatable     :: cloudflag ! cloudflag from L2 AHI
-   integer, dimension(:,:), allocatable   :: tb_qc
+   integer, dimension(:,:), allocatable   :: tb_qc, cloud_flag
    real*4,  dimension(:,:), allocatable   :: tb_obs, tb_bak, tb_inv, tb_oma, tb_err, ems, ems_jac
-   real*4,  dimension(:,:), allocatable   :: tb_bak_clr ! clear-sky brightness temp
+   real*4,  dimension(:,:), allocatable   :: cloud_mod, cloud_obs, tb_bak_clr ! clear-sky brightness temp
    real*4,  dimension(:,:), allocatable   :: weightfunc_peak ! peak of weighting function
    real*4,  dimension(:,:), allocatable   :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water
    real*4,  dimension(:,:), allocatable   :: prf_ice, prf_rain, prf_snow, prf_grau, prf_hail
@@ -139,6 +139,7 @@ program da_rad_diags
       write(0,*) trim(instid(iinst))
 
       amsr2 = index(instid(iinst),'amsr2') > 0
+      abi   = index(instid(iinst),'abi') > 0
 
       nerr = 0
       total_npixel = 0
@@ -263,6 +264,12 @@ program da_rad_diags
             allocate ( tb_oma(1:nchan,1:total_npixel) )
             allocate ( tb_err(1:nchan,1:total_npixel) )
             allocate (  tb_qc(1:nchan,1:total_npixel) )
+            if ( abi ) then
+               allocate ( cloud_mod(1:nchan,1:total_npixel) )
+               allocate ( cloud_obs(1:nchan,1:total_npixel) )
+               allocate ( cloud_flag(1:nchan,1:total_npixel))
+               cloud_flag = 0
+            end if
             allocate (    ems(1:nchan,1:total_npixel) )
             if ( jac_found ) then
                allocate ( ems_jac(1:nchan,1:total_npixel) )
@@ -333,6 +340,11 @@ program da_rad_diags
             tb_inv = missing_r
             tb_oma = missing_r
             tb_err = missing_r
+            if ( abi ) then
+               cloud_mod = missing_r
+               cloud_obs = missing_r
+            end if
+
             ncname = 'diags_'//trim(instid(iinst))//"_"//datestr1(itime)//'.nc'
             ios = NF_CREATE(trim(ncname), NF_NETCDF4, ncid)  ! Change to output netcdf4 files
            !ios = NF_CREATE(trim(ncname), NF_CLOBBER, ncid)  ! NF_CLOBBER specifies the default behavior of
@@ -392,7 +404,15 @@ program da_rad_diags
             read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_err(:,ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! QC
             read(unit=iunit(iproc),fmt='(10i11)',iostat=ios  ) tb_qc(:,ipixel)
-            read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf
+            read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf 
+            if ( abi .and. buf(1:4) == "CMOD" ) then ! read cloud_mod, cloud_obs, cloud_flag for abi
+               read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) cloud_mod(:,ipixel)
+               read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf        ! CMOD
+               read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) cloud_obs(:,ipixel)
+               read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf        ! COBS
+               read(unit=iunit(iproc),fmt='(10i11)',iostat=ios  ) cloud_flag(:,ipixel)
+               read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf        ! cloud_flag
+            end if
             if ( buf(1:4) == "INFO" ) then
                backspace(iunit(iproc))
                cycle npixel_loop
@@ -523,6 +543,13 @@ program da_rad_diags
       end if
       ios = NF_DEF_VAR(ncid, 'tb_err', NF_FLOAT, 2, ishape(1:2), varid)
       ios = NF_DEF_VAR(ncid, 'tb_qc',  NF_INT,   2, ishape(1:2), varid)
+      if ( abi ) then
+         ios = NF_DEF_VAR(ncid, 'cloud_mod', NF_FLOAT, 2, ishape(1:2), varid)
+         ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
+         ios = NF_DEF_VAR(ncid, 'cloud_obs', NF_FLOAT, 2, ishape(1:2), varid)
+         ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
+         ios = NF_DEF_VAR(ncid, 'cloud_flag',  NF_INT,   2, ishape(1:2), varid)
+      end if
       !
       ! define 2-D array with dimensions nlev * total_npixel
       !
@@ -669,6 +696,14 @@ program da_rad_diags
       ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_err)
       ios = NF_INQ_VARID (ncid, 'tb_qc', varid)
       ios = NF_PUT_VARA_INT(ncid,  varid, istart(1:2), icount(1:2), tb_qc)
+      if ( abi ) then
+         ios = NF_INQ_VARID (ncid, 'cloud_mod', varid)
+         ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), cloud_mod)
+         ios = NF_INQ_VARID (ncid, 'cloud_obs', varid)
+         ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), cloud_obs)
+         ios = NF_INQ_VARID (ncid, 'cloud_flag', varid)
+         ios = NF_PUT_VARA_INT(ncid,  varid, istart(1:2), icount(1:2), cloud_flag)
+      end if
       !
       ! output 2-D array with dimensions nlev * total_npixel
       !
@@ -890,6 +925,11 @@ program da_rad_diags
       deallocate ( tb_bak_clr )
       deallocate ( weightfunc_peak )
       deallocate ( tb_inv )
+      if ( abi ) then
+         deallocate ( cloud_mod )
+         deallocate ( cloud_obs )
+         deallocate ( cloud_flag )
+      end if
       deallocate ( tb_oma )
       deallocate ( ems )
       if ( jac_found ) deallocate ( ems_jac )

var/da/da_radiance/da_allocate_rad_iv.inc

@@ -80,6 +80,10 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
       end if
       if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then
          allocate (iv%instid(i)%cloudflag(iv%instid(i)%num_rad))
+      end if
+      if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
+         allocate (iv%instid(i)%cloud_mod(nchan,iv%instid(i)%num_rad))
+         allocate (iv%instid(i)%cloud_obs(nchan,iv%instid(i)%num_rad))
       end if	  
       if ( index(iv%instid(i)%rttovid_string, 'gmi') > 0 ) then
          allocate (iv%instid(i)%clw(iv%instid(i)%num_rad))
@@ -112,16 +116,26 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
       allocate (iv%instid(i)%gamma_jacobian(nchan,iv%instid(i)%num_rad))
       allocate (iv%instid(i)%cloud_frac(iv%instid(i)%num_rad))
       if ( use_clddet_zz ) then
-         iv%instid(i)%superob_width = 2*ahi_superob_halfwidth+1
+         ! here we assume AHI and ABI (they cover different regions) are not used simultaneously
+         if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) &
+              iv%instid(i)%superob_width = 2*ahi_superob_halfwidth+1
+         if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) &
+              iv%instid(i)%superob_width = 2*abi_superob_halfwidth+1
+
          allocate (iv%instid(i)%superob(iv%instid(i)%superob_width, &
                                         iv%instid(i)%superob_width))
          do iy = 1, iv%instid(i)%superob_width
          do ix = 1, iv%instid(i)%superob_width
             allocate (iv%instid(i)%superob(ix,iy)%cld_qc(iv%instid(i)%num_rad))
             allocate (iv%instid(i)%superob(ix,iy)%tb_obs(nchan,iv%instid(i)%num_rad))
+            if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
+              do n = 1, iv%instid(i)%num_rad
+                 allocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3(nchan))
+              end do
+            end if
          end do
          end do
-      end if	  
+      end if
       if ( use_rttov_kmatrix .or. use_crtm_kmatrix ) then
          allocate(iv%instid(i)%ts_jacobian(nchan,iv%instid(i)%num_rad))
          allocate(iv%instid(i)%ps_jacobian(nchan,iv%instid(i)%num_rad))

var/da/da_radiance/da_deallocate_radiance.inc

@@ -38,6 +38,13 @@
 	    deallocate ( satinfo(i) % clearSkyBias)
          endif
 
+         ! Deallocate extra variables for ABI
+         if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
+	    deallocate (satinfo(i) % error_cld_y)
+            deallocate (satinfo(i) % error_cld_x)
+         endif
+
+
          if (use_error_factor_rad) then
             deallocate (satinfo(i) % error_factor)
          endif
@@ -115,6 +122,10 @@
             end if
             if ( index(iv%instid(i)%rttovid_string, 'ahi') > 0 ) then
                deallocate (iv%instid(i)%cloudflag)
+            end if
+            if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
+               deallocate (iv%instid(i)%cloud_mod)
+               deallocate (iv%instid(i)%cloud_obs)
             end if			
             if ( index(iv%instid(i)%rttovid_string,'gmi') > 0 ) then
                deallocate (iv%instid(i)%clw)
@@ -149,8 +160,16 @@
             if ( use_clddet_zz ) then
                do iy = 1, iv%instid(i)%superob_width
                do ix = 1, iv%instid(i)%superob_width
-                  deallocate (iv%instid(i)%superob(ix,iy)%cld_qc)
-                  deallocate (iv%instid(i)%superob(ix,iy)%tb_obs)
+                 if ( index(iv%instid(i)%rttovid_string, 'abi') > 0 ) then
+                  do n = 1,iv%instid(i)%num_rad
+                     if ( allocated (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3) ) &
+                        deallocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%tb_stddev_3x3)
+                     if ( allocated (iv%instid(i)%superob(ix,iy)%cld_qc(n)%CIRH2O_abi) ) &
+                        deallocate (iv%instid(i)%superob(ix,iy)%cld_qc(n)%CIRH2O_abi)
+                  end do
+                 end if
+                 deallocate (iv%instid(i)%superob(ix,iy)%cld_qc)
+                 deallocate (iv%instid(i)%superob(ix,iy)%tb_obs)
                end do
                end do
                deallocate (iv%instid(i)%superob)

var/da/da_radiance/da_get_innov_vector_crtm.inc

@@ -92,7 +92,7 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
    real, allocatable :: hessian(:,:)
    real*8, allocatable :: eignvec(:,:), eignval(:)
    real              :: rad_clr, rad_ovc_ilev, rad_ovc_jlev
-   
+
    integer           :: Band_Size(5), Bands(AIRS_Max_Channels,5) 
 !For Zhuge and Zou cloud detection
    real, allocatable :: geoht_full(:,:,:)
@@ -243,9 +243,10 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
       calc_tb_clr = .false.
       if ( crtm_cloud .and. &
            ( trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' .or. &
+             trim( crtm_sensor_name(rtminit_sensor(inst))) == 'abi' .or. &
              trim( crtm_sensor_name(rtminit_sensor(inst))) == 'ahi') ) then
          !Tb_clear_sky is only needed for symmetric obs error model
-         !symmetric obs error model only implemented for amsr2 for now
+         !symmetric obs error model only implemented for amsr2 & abi/ahi for now
          calc_tb_clr = .true.
       end if
 
@@ -443,7 +444,6 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
             call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
          end if
 
-
          call da_interp_2d_partial (grid%xb%u10,  iv%instid(inst)%info, 1, n, n, model_u10(n:n))
          call da_interp_2d_partial (grid%xb%v10,  iv%instid(inst)%info, 1, n, n, model_v10(n:n))
          call da_interp_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
@@ -476,6 +476,14 @@ subroutine da_get_innov_vector_crtm ( it, grid, ob, iv )
                cycle pixel_loop
             end if
          end do
+         !if ( all(ob%instid(inst)%tb(1:nchanl,n) < 0.) ) then
+         !   write(message(1),'(a,2i5.0,a)') ' Skipping the pixel at loc ', i, j,  &
+         !      ' where all observed BTs are < 0'
+         !   call da_warning(__FILE__,__LINE__,message(1:1))
+         !   iv%instid(inst)%tb_inv(:,n) = missing_r
+         !   iv%instid(inst)%info%proc_domain(:,n) = .false.
+         !   cycle pixel_loop
+         !end if
 
          ! convert cloud content unit from kg/kg to kg/m^2        
          if (crtm_cloud) then

var/da/da_radiance/da_get_innov_vector_rttov.inc

@@ -49,12 +49,30 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
    real, allocatable :: em_mspps(:)    ! emissivity caluclated using MSPPS algorithm
    real              :: ts_mspps       ! surface temperature calcualted using MSPPS algorithm
 
+   !For Zhuge and Zou cloud detection
+   real, allocatable :: geoht_full(:,:,:)
+   real              :: geoht_pixel(kts:min(kte,kme-1))
+   real              :: tt_pixel(kts:min(kte,kme-1))
+   real              :: pp_pixel(kts:min(kte,kme-1))
+
    if (trace_use) call da_trace_entry("da_get_innov_vector_rttov")
 
    !------------------------------------------------------
    ! [1.0] calculate the background bright temperature
    !-------------------------------------------------------
 
+   if ( use_clddet_zz ) then
+      allocate ( geoht_full(ims:ime,jms:jme,kms:kme-1) )
+      do k = kms, kme-1
+         do j = jms, jme
+            do i = ims, ime
+               geoht_full(i,j,k) = 0.5 * ( grid%ph_2(i,j,k)   + grid%phb(i,j,k) + &
+                                           grid%ph_2(i,j,k+1) + grid%phb(i,j,k+1) ) / gravity
+            end do
+         end do
+      end do
+   end if
+
    do inst = 1, iv%num_inst                 ! loop for sensor
       if ( iv%instid(inst)%num_rad < 1 ) cycle
       nlevels  = iv%instid(inst)%nlevels
@@ -99,7 +117,6 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
       call da_interp_lin_3d (grid%xb%t, iv%instid(inst)%info, iv%instid(inst)%t (:,n1:n2))
       call da_interp_lin_3d (grid%xb%q, iv%instid(inst)%info, iv%instid(inst)%mr(:,n1:n2))
 
-
       do n= n1,n2
          do k=1, nlevels
             if (iv%instid(inst)%info%zk(k,n) <= 0.0) then
@@ -132,6 +149,19 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
             iv%instid(inst)%surftype(n) = 0
          end if
 
+         if ( use_clddet_zz ) then
+            ! Find tropopause temperature for Zhuge and Zou Cloud Detection
+            do k = kts, min(kte,kme-1)
+               call da_interp_2d_partial ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
+               call da_interp_2d_partial ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
+               call da_interp_2d_partial ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
+
+!               call da_interp_lin_2d ( grid%xb%t(:,:,k), iv%instid(inst)%info, k, n, n, tt_pixel(k) )
+!               call da_interp_lin_2d ( grid%xb%p(:,:,k), iv%instid(inst)%info, k, n, n, pp_pixel(k) )
+!               call da_interp_lin_2d ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
+            end do
+            call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
+         end if
       end do
 
       call da_interp_lin_2d (grid%xb % u10,  iv%instid(inst)%info, 1, iv%instid(inst)%u10(n1:n2))
@@ -381,6 +411,8 @@ real,allocatable :: temp(:), temp2(:), temp3(:,:)
 
    end do        ! end loop for sensor
 
+   if ( use_clddet_zz ) deallocate ( geoht_full )
+
    if (trace_use) call da_trace_exit("da_get_innov_vector_rttov")
 #else
     call da_error(__FILE__,__LINE__, &