GSI changes to read and assimilate IASI-NG

src/gsi/calc_fov_crosstrk.f90

@@ -1289,6 +1289,8 @@ subroutine get_sat_height(satid, height, valid)
        height=840._r_kind
     case('n20', 'n21', 'n22', 'n23')
        height=840._r_kind
+    case('metop-sg-a1', 'metop-sg-a2', 'metop-sg-a3') 
+       height=840._r_kind
     case default
        write(6,*) 'GET_SAT_HEIGHT: ERROR, unrecognized satellite id: ', trim(satid)
        valid=.false.

src/gsi/crtm_interface.f90

@@ -667,6 +667,21 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo
    error_status = crtm_channelinfo_subset(channelinfo(1), &
         channel_subset = nuchan(subset_start:subset_end))
 
+!TODO JAJ the correct one    else if (channelinfo(1)%sensor_id(1:7) == 'iasi-ng' .AND. isis(1:7) == 'iasi-ng') then
+else if (channelinfo(1)%sensor_id(1:3) == '999' .AND. isis(1:7) == 'iasi-ng') then  ! TODO To be removed.
+   sensorindex = 1
+   subset_start = 0
+   subset_end = 0
+   do k=1, jpch_rad
+     if (isis == nusis(k)) then
+       if (subset_start == 0) subset_start = k
+       subset_end = k
+     endif
+   end do
+
+   error_status = crtm_channelinfo_subset(channelinfo(1), &
+        channel_subset = nuchan(subset_start:subset_end))
+
 else if (channelinfo(1)%sensor_id(1:4) == 'iasi' .AND. isis(1:4) == 'iasi') then
    sensorindex = 1
    subset_start = 0

src/gsi/gsi_files.cmake

@@ -497,6 +497,7 @@ read_goesndr.f90
 read_gps.f90
 read_guess.F90
 read_iasi.f90
+read_iasing.f90
 read_l2bufr_mod.f90
 read_lag.f90
 read_lidar.f90

src/gsi/gsi_obOperTypeManager.F90

@@ -327,6 +327,7 @@ function dtype2index_(dtype) result(index_)
     case("hsb"    ); index_= iobOper_rad
         !
     case("iasi"   ); index_= iobOper_rad
+    case("iasi-ng"); index_= iobOper_rad
     case("cris"   ); index_= iobOper_rad
     case("cris-fsr"  ); index_= iobOper_rad
         !
@@ -357,6 +358,7 @@ function dtype2index_(dtype) result(index_)
         !
     case("avhrr_navy"); index_= iobOper_rad
     case("avhrr"  ); index_= iobOper_rad
+    case("metimage" ); index_= iobOper_rad
     case("viirs-m"  ); index_= iobOper_rad
 
   case("tcp"    ,"[tcpoper]"    ); index_= iobOper_tcp

src/gsi/gsimod.F90

@@ -94,7 +94,7 @@ module gsimod
       init_qcvars,vadfile,noiqc,c_varqc,gps_jacqc,qc_noirjaco3,qc_noirjaco3_pole,&
       buddycheck_t,buddydiag_save,njqc,vqc,nvqc,hub_norm,vadwnd_l2rw_qc, &
       pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check, &
-      cris_cads, iasi_cads, airs_cads
+      cris_cads, iasi_cads, iasing_cads, airs_cads
   use qcmod, only: troflg,lat_c,nrand
   use cads, only: M__Sensor,N__Num_Bands,N__GradChkInterval,N__Band_Size,N__Bands,N__Window_Width, &
       N__Window_Bounds,R__BT_Threshold,R__Grad_Threshold,R__Window_Grad_Threshold, L__Do_Quick_Exit, &
@@ -1064,9 +1064,10 @@ module gsimod
 !
 !     Flags to use the new IR cloud detection routine.  Flag must be set to true to use the new routine.  The default
 !     (no flag or .false.) will use the default.
-!     airs_cads: use the clod and aerosool detection software for the AIRS instrument
-!     cris_cads: use the cloud and aerosol detection software for CrIS instruments
-!     iasi_cads: use the cloud and aerosol detection software for IASI instruments
+!     airs_cads  : use the cloud and aerosol detection software for AIRS instrument
+!     cris_cads  : use the cloud and aerosol detection software for CrIS instruments
+!     iasi_cads  : use the cloud and aerosol detection software for IASI instruments
+!     iasing_cads: use the cloud and aerosol detection software for IASI-NG instruments
 !     
 
   namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,&
@@ -1078,7 +1079,7 @@ module gsimod
        q_doe_a_136,q_doe_a_137,q_doe_b_136,q_doe_b_137, &
        t_doe_a_136,t_doe_a_137,t_doe_b_136,t_doe_b_137, &
        uv_doe_a_236,uv_doe_a_237,uv_doe_a_213,uv_doe_b_236,uv_doe_b_237,uv_doe_b_213, &
-       vad_near_analtime,airs_cads,cris_cads,iasi_cads
+       vad_near_analtime,airs_cads,cris_cads,iasi_cads, iasing_cads
 
 ! OBS_INPUT (controls input data):
 !      dmesh(max(dthin))- thinning mesh for each group
@@ -2320,7 +2321,7 @@ subroutine gsimain_initialize
      endif
      do i=1,ndat
         write(6,401)dfile(i),dtype(i),dplat(i),dsis(i),dval(i),dthin(i),dsfcalc(i),time_window(i)
- 401    format(1x,a20,1x,a10,1x,a11,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2)
+ 401    format(1x,a20,1x,a10,1x,a12,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2)
      end do
      write(6,superob_radar)
      write(6,lag_data)

src/gsi/obsmod.F90

@@ -1186,7 +1186,7 @@ subroutine init_obsmod_vars(nhr_assim,mype)
        endif
 ! for cris, iasi, atms, regional analysis may want shorter time window
        if (index(dtype(ii),'cris') /= 0 .or. index(dtype(ii),'atms') /= 0 .or. &
-           index(dtype(ii),'iasi') /= 0 ) then
+           index(dtype(ii),'iasi') /= 0 .or. index(dtype(ii),'iasi-ng') /= 0) then
           if(time_window(ii)>time_window_rad) then
              time_window(ii) = time_window_rad
              if (mype==0) write(6,*) 'INIT_OBSMOD_VARS: reset time window for ',dtype(ii),&

src/gsi/qcmod.f90

@@ -118,6 +118,7 @@ module qcmod
 !   def airs_cads       - if true, use the cloud and aerosol detection routine for Aqua/AIRS instrument
 !   def cris_cads       - if true, use the cloud and aerosol detection routine for CrIS instruments
 !   def iasi_cads       - if true, use the cloud and aerosol detection routine for IASI instruments
+!   def iasing_cads     - if true, use the cloud and aerosol detection routine for IASI-NG instruments
 !
 ! following used for nonlinear qc:
 !
@@ -204,7 +205,7 @@ module qcmod
   public :: troflg
   public :: lat_c
   public :: nrand 
-  public :: airs_cads, cris_cads, iasi_cads
+  public :: airs_cads, cris_cads, iasi_cads, iasing_cads
 
   logical nlnqc_iter,njqc,vqc,nvqc,hub_norm
   logical noiqc
@@ -220,7 +221,7 @@ module qcmod
   logical vadwnd_l2rw_qc
   logical troflg
   logical cao_check
-  logical airs_cads, cris_cads, iasi_cads
+  logical airs_cads, cris_cads, iasi_cads, iasing_cads
 
   character(10):: vadfile
   integer(i_kind) npres_print
@@ -461,9 +462,10 @@ subroutine init_qcvars
     lat_c=21.0_r_kind
     nrand=13
 
-    airs_cads = .false.
-    cris_cads = .false.
-    iasi_cads = .false.
+    airs_cads   = .false.
+    cris_cads   = .false.
+    iasi_cads   = .false.
+    iasing_cads = .false.
 
     return
   end subroutine init_qcvars
@@ -598,8 +600,8 @@ subroutine setup_tzr_qc(obstype)
       tzchk = 0.85_r_kind
     elseif (  obstype == 'hirs2' .or. obstype == 'hirs3' .or. obstype == 'hirs4' .or. & 
               obstype == 'sndr' .or. obstype == 'sndrd1' .or. obstype == 'sndrd2'.or. &
-              obstype == 'sndrd3' .or. obstype == 'sndrd4' .or.  &
-              obstype == 'goes_img' .or. obstype == 'ahi' .or. obstype == 'airs' .or. obstype == 'iasi' .or. &
+              obstype == 'sndrd3' .or. obstype == 'sndrd4' .or.  obstype == 'goes_img' .or. &
+              obstype == 'ahi' .or. obstype == 'airs' .or. obstype == 'iasi' .or. obstype == 'iasi-ng' .or.&
               obstype == 'cris' .or. obstype == 'cris-fsr' .or. obstype == 'seviri'  .or. obstype == 'abi') then
       tzchk = 0.85_r_kind
     endif
@@ -2076,20 +2078,20 @@ subroutine qc_saphir(nchanl,sfchgt,luse,sea, &
 end subroutine qc_saphir
 
 subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs,                         &
-     cris,iasi,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, &
+     cris,iasi,iasing,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, &
      wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim,                   &
      id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,cluster_fraction,    &
      cluster_bt, chan_stdev, model_bt)
 !    id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,radmod) ! all-sky
 
 !$$$ subprogram documentation block
 !               .      .    .
-! subprogram:  qc_irsnd    QC for ir sounder data(hirs,goessndr,airs,iasi,cris)
+! subprogram:  qc_irsnd    QC for ir sounder data(hirs,goessndr,airs,iasi,iasing,cris)
 !
 !   prgmmr: derber           org: np23            date: 2010-08-20
 !
 ! abstract: set quality control criteria for ir sounder data (hirs, 
-!          goessndr, airs, iasi, cris)
+!          goessndr, airs, iasi, iasing, cris)
 !
 ! program history log:
 !     2010-08-10  derber transfered from setuprad
@@ -2110,6 +2112,8 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
 !     goessndr     - logical flag - if goessndr data - true
 !     cris         - logical flag - if cris data - true
 !     avhrr        - logical flag - if avhrr data - true
+!     iasi         - logical flag - if iasi data - true
+!     iasing       - logical flag - if iasing data - true
 !     zsges        - elevation of guess
 !     cenlat       - latitude of observation
 !     frac_sea     - fraction of grid box covered with water
@@ -2164,7 +2168,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
 
 ! Declare passed variables
 
-  logical,                            intent(in   ) :: sea,land,ice,snow,luse,goessndr,airs,cris,hirs,iasi
+  logical,                            intent(in   ) :: sea,land,ice,snow,luse,goessndr,airs,cris,hirs,iasi,iasing
   logical,                            intent(inout) :: zero_irjaco3_pole
   integer(i_kind),                    intent(in   ) :: nsig,nchanl,ndat,is
   integer(i_kind),dimension(nchanl),  intent(in   ) :: ich
@@ -2326,6 +2330,21 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
              tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, &
              cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp )
 
+  elseif ( iasing .and. iasing_cads ) then
+      I_Sensor_ID = 59
+      chan_array = nuchan(ich)                  ! channel numbers
+      tb_bc = tbc + tsim                        ! observation BT with bias correction
+      boundary_layer_pres = nint(0.8_r_kind*prsltmp(1))  !  boundary layer set to be 80% of surface pressure
+      tropopause_height = nint(trop5)
+      imager_chans = (/18,19/)                  ! imager channel numbers (from satinfo)
+      isurface_chan = 2539                      ! surface channel
+      ichan_10_micron = 2343                    ! ~10.7 micron channel for low level cloud test
+      ichan_12_micron = 1509                    ! ~12.0 micron channel for low level cloud test
+
+      call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, &
+             tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, &
+             cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp )
+
   elseif ( airs .and. airs_cads ) then
       I_Sensor_ID = 11
       chan_array = nuchan(ich)                  ! channel numbers
@@ -2348,7 +2367,8 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
 
 ! compute cloud stats 
 ! If using CADS
-  if ((cris .and. cris_cads) .or. (iasi .and. iasi_cads) .or. (airs .and. airs_cads)) then
+  if ((cris .and. cris_cads) .or. (iasi .and. iasi_cads) .or. (airs .and. airs_cads) .or. &
+       iasing .and. iasing_cads ) then
 
 !   Reject channels affected by clouds
     do i=1, nchanl
@@ -2416,7 +2436,6 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
 ! default compute cloud stats, emc_legacy_cloud_detect 
   else  
     if ( lcloud > 0 ) then
-
       do i=1,nchanl
 !       reject channels with iuse_rad(j)=-1 when they are peaking below the cloud
         j=ich(i)

src/gsi/radiance_mod.f90

@@ -412,7 +412,8 @@ subroutine radiance_obstype_init
           rtype(i) == 'avhrr'  .or. rtype(i) == 'amsre'    .or.  rtype(i) == 'ssmis'  .or. & 
           rtype(i) == 'ssmi'   .or. rtype(i) == 'atms'     .or.  rtype(i) == 'cris'   .or. & 
           rtype(i) == 'amsr2'  .or. rtype(i) == 'gmi'      .or.  rtype(i) == 'saphir' .or. &
-          rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi'    .or.  rtype(i) == 'viirs' ) then
+          rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi'    .or.  rtype(i) == 'viirs'  .or. &
+          rtype(i) == 'iasi-ng' )then
           drtype(i)='rads'
        end if
     end do

src/gsi/radinfo.f90

@@ -768,7 +768,8 @@ subroutine radinfo_read
        end select 
 
        if ( .not. diag_rad .and. iuse_rad(j) < 0 .and. iextra_det(j) < 0 .and. &
-          ( nusis(j)(1:4) == 'cris' .or. nusis(j)(1:4) == 'iasi' .or. nusis(j)(1:4) == 'airs')) cycle
+          ( nusis(j)(1:4) == 'cris' .or. nusis(j)(1:7) == 'iasi-ng' .or. nusis(j)(1:4) == 'iasi'  &
+            .or. nusis(j)(1:4) == 'airs')) cycle
 
        if(iuse_rad(j) == 4 .or. iuse_rad(j) == 2) air_rad(j)=zero
        if(iuse_rad(j) == 4 .or. iuse_rad(j) == 3) ang_rad(j)=zero
@@ -869,7 +870,7 @@ subroutine radinfo_read
 !            The second part of the if statement keeps from printing them.
              if ( .not. cfound ) then
                 if ((diag_rad .and. mype ==0) .or. &
-                   (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:4) /= 'iasi')) &
+                   (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:7) /= 'iasi-ng' .and. isis(1:4) /= 'iasi')) &
                    write(6,*) '***WARNING instrument/channel ',isis,ichan,'found in satbias_pc file but not found in satinfo'
              endif
 
@@ -1112,7 +1113,7 @@ subroutine radinfo_read
 !         The second part of the if statement keeps from printing them.
           if ( .not. cfound ) then
              if ((diag_rad .and. mype ==0) .or. &
-                (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:4) /= 'iasi')) &
+                (.not. diag_rad .and. isis(1:4)/='airs' .and. isis(1:4) /= 'cris' .and. isis(1:7) /= 'iasi-ng' .and. isis(1:4) /= 'iasi')) &
                 write(6,*) '***WARNING instrument/channel ',isis,ichan,'found in satbias_in file but not found in satinfo'
           endif