Minor changes to support FPT, static analysis regression test

Registry/registry.var

@@ -451,6 +451,8 @@ rconfig   real      varbc_factor            namelist,wrfvar14  1  1.0      - "va
 rconfig   integer   varbc_nbgerr            namelist,wrfvar14  1  5000     - "varbc_nbgerr"            ""  ""
 rconfig   integer   varbc_nobsmin           namelist,wrfvar14  1  10       - "varbc_nobsmin"           ""  ""
 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   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"       ""  ""

var/build/da.make

@@ -69,6 +69,7 @@ WRFVAR_OBJS = \
    gsi_kinds.o \
    gsi_constants.o \
    gsi_thinning.o \
+   mod_clddet_geoir.o \
    da_wrfvar_io.o \
    da_airsr.o \
    da_wrfvar_top.o \

var/build/depend.txt

@@ -110,7 +110,7 @@ da_blas.o : da_blas.f90 dgemv.inc dgemm.inc xerbla.inc lsame.inc dswap.inc dtrmm
 da_bogus.o : da_bogus.f90 da_calculate_grady_bogus.inc da_get_innov_vector_bogus.inc da_check_max_iv_bogus.inc da_transform_xtoy_bogus_adj.inc da_transform_xtoy_bogus.inc da_print_stats_bogus.inc da_oi_stats_bogus.inc da_residual_bogus.inc da_jo_and_grady_bogus.inc da_ao_stats_bogus.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_par_util1.o da_par_util.o da_interpolation.o da_define_structures.o da_control.o module_domain.o 
 da_buoy.o : da_buoy.f90 da_calculate_grady_buoy.inc da_get_innov_vector_buoy.inc da_check_max_iv_buoy.inc da_transform_xtoy_buoy_adj.inc da_transform_xtoy_buoy.inc da_print_stats_buoy.inc da_oi_stats_buoy.inc da_residual_buoy.inc da_jo_and_grady_buoy.inc da_ao_stats_buoy.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_grid_definitions.o da_par_util.o da_par_util1.o da_interpolation.o da_define_structures.o da_control.o module_domain.o 
 da_control.o : da_control.f90 module_driver_constants.o 
-da_crtm.o : da_crtm.f90 da_det_crtm_climat.inc da_crtm_sensor_descriptor.inc da_crtm_init.inc da_crtm_ad.inc da_crtm_direct.inc da_crtm_k.inc da_crtm_tl.inc da_get_innov_vector_crtm.inc da_transform_xtoy_crtm_adj.inc da_transform_xtoy_crtm.inc da_tracing.o da_tools.o da_tools_serial.o da_reporting.o da_radiance1.o module_dm.o da_interpolation.o da_control.o module_radiance.o da_define_structures.o module_domain.o 
+da_crtm.o : da_crtm.f90 da_det_crtm_climat.inc da_crtm_sensor_descriptor.inc da_crtm_init.inc da_crtm_ad.inc da_crtm_direct.inc da_crtm_k.inc da_crtm_tl.inc da_get_innov_vector_crtm.inc da_transform_xtoy_crtm_adj.inc da_transform_xtoy_crtm.inc da_tracing.o da_tools.o da_tools_serial.o da_reporting.o da_radiance1.o module_dm.o da_physics.o da_interpolation.o da_control.o module_radiance.o da_define_structures.o module_domain.o 
 da_define_structures.o : da_define_structures.f90 da_gauss_noise.inc da_random_seed.inc da_initialize_cv.inc da_zero_vp_type.inc da_zero_y.inc da_zero_x.inc da_deallocate_y.inc da_deallocate_observations.inc da_deallocate_background_errors.inc da_allocate_y.inc da_allocate_observations.inc da_allocate_background_errors.inc da_wavelet.o da_reporting.o da_tools_serial.o da_tracing.o da_control.o module_domain.o da_allocate_y_rain.inc da_allocate_y_radar.inc da_allocate_observations_rain.inc da_allocate_obs_info.inc 
 da_dynamics.o : da_dynamics.f90 da_wz_base.inc da_uv_to_vorticity.inc da_w_adjustment_adj.inc da_w_adjustment_lin.inc da_uv_to_divergence_adj.inc da_uv_to_divergence.inc da_psichi_to_uv_adj.inc da_psichi_to_uv.inc da_hydrostaticp_to_rho_lin.inc da_hydrostaticp_to_rho_adj.inc da_balance_geoterm_lin.inc da_balance_geoterm_adj.inc da_balance_equation_lin.inc da_balance_equation_adj.inc da_balance_cycloterm_lin.inc da_balance_cycloterm_adj.inc da_balance_cycloterm.inc da_wpec_constraint.inc da_wpec_constraint_adj.inc da_wpec_constraint_cycloterm.inc da_wpec_constraint_geoterm.inc da_wpec_constraint_lin.inc da_tools.o da_tracing.o da_ffts.o da_reporting.o da_define_structures.o module_comm_dm.o module_dm.o module_domain.o da_control.o da_divergence_constraint.inc da_divergence_constraint_adj.inc 
 da_etkf.o : da_etkf.f90 da_solve_etkf.inc da_matmultiover.inc da_matmulti.inc da_innerprod.inc da_lapack.o da_gen_be.o da_control.o 
@@ -230,6 +230,7 @@ gen_spectra.o : gen_spectra.f90 da_tracing.o da_spectral.o da_control.o
 gsi_constants.o : gsi_constants.f90 init_constants_derived.inc gsi_kinds.o 
 gsi_kinds.o : gsi_kinds.f90 
 gsi_thinning.o : gsi_thinning.f90 gsi_constants.o gsi_kinds.o 
+mod_clddet_geoir.o : mod_clddet_geoir.f90 da_control.o
 kma2netcdf.o : kma2netcdf.f90 module_kma_wave2grid.o module_kma2netcdf_interface.o da_wrfvar_io.o module_dm.o module_wrf_error.o module_timing.o module_driver_constants.o module_machine.o da_control.o module_domain.o module_configure.o 
 module_ffts.o : module_ffts.f90 qpassm.inc fft661.inc fft551.inc 
 module_kma2netcdf_interface.o : module_kma2netcdf_interface.f90 module_tiles.o module_configure.o module_driver_constants.o module_timing.o module_domain.o module_kma_wave2grid.o 

var/da/da_define_structures/da_define_structures.f90

@@ -508,7 +508,16 @@ module da_define_structures
       real,    pointer     :: bgerr(:) 
       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(:,:,:)
+   end type clddet_geoir_type   
+   type superob_type
+     real, allocatable :: tb_obs(:,:)
+     type(clddet_geoir_type), allocatable :: cld_qc(:)
+   end type superob_type
    type cv_index_type
       integer              :: ts
       integer              :: nclouds
@@ -544,10 +553,12 @@ module da_define_structures
       integer, pointer     :: cloud_flag(:,:)
       integer, pointer     :: cloudflag(:)
       integer, pointer     :: rain_flag(:)
+      real, allocatable    :: cloud_frac(:)
       real,    pointer     :: satzen(:) 
       real,    pointer     :: satazi(:) 
       real,    pointer     :: solzen(:) 
       real,    pointer     :: solazi(:) 
+      real,    pointer     :: tropt(:)  !! Tropopause temperature, K.
       real,    pointer     :: t(:,:)
       real,    pointer     :: q(:,:)
       real,    pointer     :: mr(:,:)
@@ -623,11 +634,12 @@ module da_define_structures
       real,    pointer     :: ice_coverage(:)
       real,    pointer     :: snow_coverage(:)
       integer, pointer     :: crtm_climat(:) ! CRTM only
-
+      integer              :: superob_width = 1
       type (varbc_info_type)        :: varbc_info
       type (varbc_type),pointer     :: varbc(:)
       type (cv_index_type), pointer :: cv_index(:)
       type (infa_type)              :: info
+      type (superob_type), allocatable :: superob(:,:)
    end type instid_type
 
    type iv_type

var/da/da_monitor/da_rad_diags.f90

@@ -58,8 +58,8 @@ program da_rad_diags
    integer, dimension(:), allocatable     :: landsea_mask, soiltyp, vegtyp
    real*4,  dimension(:), allocatable     :: lat, lon, elv, elev
    real*4,  dimension(:), allocatable     :: ret_clw
-   real*4,  dimension(:), allocatable     :: satzen, satazi, t2m, mr2m, u10, v10, ps, ts
-   real*4,  dimension(:), allocatable     :: smois, tslb, snowh, vegfra, clwp
+   real*4,  dimension(:), allocatable     :: satzen, satazi, solzen, solazi, t2m, mr2m, u10, v10, ps, ts
+   real*4,  dimension(:), allocatable     :: smois, tslb, snowh, vegfra, clwp, cloud_frac
    integer, dimension(:,:), allocatable   :: tb_qc
    real*4,  dimension(:,:), allocatable   :: tb_obs, tb_bak, tb_inv, tb_oma, tb_err, ems, ems_jac
    real*4,  dimension(:,:), allocatable   :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water
@@ -230,6 +230,8 @@ program da_rad_diags
             allocate (          lon(1:total_npixel) )
             allocate (       satzen(1:total_npixel) )
             allocate (       satazi(1:total_npixel) )
+            allocate (       solzen(1:total_npixel) )
+            allocate (       solazi(1:total_npixel) )
             allocate (      ret_clw(1:total_npixel) ) !obs retrieved clw
             allocate (          t2m(1:total_npixel) )
             allocate (         mr2m(1:total_npixel) )
@@ -246,6 +248,7 @@ program da_rad_diags
             allocate (       vegfra(1:total_npixel) )
             allocate (         elev(1:total_npixel) )
             allocate (         clwp(1:total_npixel) ) !model/guess clwp
+            allocate (   cloud_frac(1:total_npixel) )
             allocate ( tb_obs(1:nchan,1:total_npixel) )
             allocate ( tb_bak(1:nchan,1:total_npixel) )
             allocate ( tb_inv(1:nchan,1:total_npixel) )
@@ -341,14 +344,15 @@ program da_rad_diags
 
          npixel_loop: do ipixel = ips, ipe
 
-            read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i6,i3,f6.0,4f8.2,f8.3)',iostat=ios)  &
+            read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i6,i3,f6.0,6f8.2,f8.3)',iostat=ios)  &
                n, datestr2(ipixel), scanpos(ipixel), landsea_mask(ipixel), elv(ipixel), &
-               lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel), ret_clw(ipixel)
+               lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel), solzen(ipixel), &
+               solazi(ipixel), ret_clw(ipixel)
 
-            read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2)',iostat=ios)  &
+            read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2,f15.5)',iostat=ios)  &
                t2m(ipixel), mr2m(ipixel), u10(ipixel), v10(ipixel), ps(ipixel), ts(ipixel), &
                smois(ipixel), tslb(ipixel), snowh(ipixel), isflg(ipixel), soiltyp(ipixel),  &
-               vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel)
+               vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel), cloud_frac(ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! OBS
             read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_obs(:,ipixel)
             read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf           ! BAK
@@ -580,6 +584,8 @@ program da_rad_diags
       ios = NF_DEF_VAR(ncid, 'lon',          NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'satzen',       NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'satazi',       NF_FLOAT, 1, ishape(1), varid)
+      ios = NF_DEF_VAR(ncid, 'solzen',       NF_FLOAT, 1, ishape(1), varid)
+      ios = NF_DEF_VAR(ncid, 'solazi',       NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 't2m',          NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'mr2m',         NF_FLOAT, 1, ishape(1), varid)
       ios = NF_DEF_VAR(ncid, 'u10',          NF_FLOAT, 1, ishape(1), varid)
@@ -598,6 +604,7 @@ program da_rad_diags
       if ( amsr2 ) then
          ios = NF_DEF_VAR(ncid, 'ret_clw',   NF_FLOAT, 1, ishape(1), varid)
       end if
+      ios = NF_DEF_VAR(ncid, 'cloud_frac',   NF_FLOAT, 1, ishape(1), varid)
 
       ios = NF_ENDDEF(ncid)
       !
@@ -768,6 +775,10 @@ program da_rad_diags
       ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), satzen)
       ios = NF_INQ_VARID (ncid, 'satazi', varid)
       ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), satazi)
+      ios = NF_INQ_VARID (ncid, 'solzen', varid)
+      ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), solzen)
+      ios = NF_INQ_VARID (ncid, 'solazi', varid)
+      ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), solazi)
       ios = NF_INQ_VARID (ncid, 't2m', varid)
       ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), t2m)
       ios = NF_INQ_VARID (ncid, 'mr2m', varid)
@@ -802,6 +813,8 @@ program da_rad_diags
          ios = NF_INQ_VARID (ncid, 'ret_clw', varid)
          ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), ret_clw)
       end if
+      ios = NF_INQ_VARID (ncid, 'cloud_frac', varid)
+      ios = NF_PUT_VARA_REAL(ncid,  varid, istart(2), icount(2), cloud_frac)
       !
       ios = NF_CLOSE(ncid)
 
@@ -818,6 +831,8 @@ program da_rad_diags
       deallocate ( lon )
       deallocate ( satzen )
       deallocate ( satazi )
+      deallocate ( solzen )
+      deallocate ( solazi )
       deallocate ( t2m )
       deallocate ( mr2m )
       deallocate ( u10 )
@@ -833,6 +848,7 @@ program da_rad_diags
       deallocate ( vegfra )
       deallocate ( elev )
       deallocate ( clwp )
+      deallocate ( cloud_frac )
       deallocate ( ret_clw )
       deallocate ( tb_obs )
       deallocate ( tb_bak )

var/da/da_physics/da_physics.f90

@@ -107,6 +107,6 @@ module da_physics
 #include "da_integrat_dz.inc"
 #include "da_uv_to_sd_lin.inc"
 #include "da_uv_to_sd_adj.inc"
-
+#include "da_trop_wmo.inc"
 end module da_physics
 

var/da/da_physics/da_trop_wmo.inc

@@ -0,0 +1,114 @@
+subroutine da_trop_wmo (t, z, p, nlev, tropt, tropp, tropk)
+
+   !----------------------------------------------------------------------------
+   ! * Computes tropopause T, P, and/or level based on code from Cameron Homeyer
+   !     and WMO definition
+   !
+   ! * WMO tropopause definition:
+   !     The boundary between the troposphere and the stratosphere, where an 
+   !     abrupt change in lapse rate usually occurs. It is defined as the lowest 
+   !     level at which the lapse rate decreases to 2 °C/km or less, provided 
+   !     that the average lapse rate between this level and all higher levels 
+   !     within 2 km does not exceed 2 °C/km.
+   !----------------------------------------------------------------------------
+
+   implicit none 
+
+   ! Assumed shape inputs for single column (size=nlev)
+   !  ordered from bottom to top of model
+   real, intent(in) :: t(:)             ! Temperature, K. (3D)
+
+!JJG: "z" is supposed to be height, not geopotential height. Does it matter?
+   real, intent(in) :: z(:)             ! Geopotential height above m.s.l., m.
+   real, intent(in) :: p(:)             ! Pressure, mb.
+   real, optional, intent(out) :: tropt ! Tropopause temperature, K.
+   real, optional, intent(out) :: tropp ! Tropopause pressure, mb.
+   integer, optional, intent(out) :: tropk
+   integer, intent(in) :: nlev
+
+   real    :: dtdz,laps  !LAPS
+   integer :: dtdztest(nlev), ztest(nlev)
+   integer :: i, j, k, kk, ktrop
+
+!   real, parameter :: tropz_min = 5000.0
+!   real, parameter :: tropz_max = 19000.0
+
+   if (.not.present(tropt) .and. &
+       .not.present(tropp) .and. &
+       .not.present(tropk)) return
+
+   if (present(tropt)) tropt = missing_r
+   if (present(tropp)) tropp = missing_r
+
+   !Loop over levels to find tropopause (single column)
+   ktrop = nlev-1
+   trop_loop: do  k = 1, nlev-1
+      if ( p(k) .le. 500.0 ) then
+         ! Compute lapse rate (-dT/dz)
+         dtdz = ( t(k+1) - t(k)   ) / &
+                ( z(k)   - z(k+1) )
+      else
+         ! Set lapse rate for p > 500 hPa
+         dtdz = 999.9
+      endif
+      !Check if local lapse rate <= 2 K/km
+      if (dtdz .le. 0.002) then
+         ! Initialize lapse rate and altitude test arrays
+         dtdztest = 0
+         ztest = 0
+
+         ! Compute average lapse rate across levels above current candidate 
+         do kk = k+1, nlev-1
+            dtdz = ( t(kk+1) - t(k)   ) / &
+                   ( z(k)    - z(kk+1) )
+
+            !If avg. lapse rate <= 2 K/km and z <= trop + 2 km, set pass flag
+            if ( ( dtdz .le. 0.002 ) .and. &
+                 ( (z(k) - z(kk)) .le. 2000. ) ) THEN
+               dtdztest(kk) = 1                        
+            endif
+
+            ! If z <= trop + 2 km, set pass flag
+            IF ( (z(k) - z(kk)) .le. 2000.0 ) THEN
+               ztest(kk) = 1                           
+            endif
+         enddo   !kk loop
+         laps=dtdz !LAPS
+         IF (SUM(dtdztest) .eq. SUM(ztest)) THEN
+            ! If qualified as tropopause, set altitude index and return value
+            ktrop = k  
+            exit trop_loop
+         ENDIF
+      ENDIF
+   end do trop_loop
+
+!   ! Filter ktrop using tpause height thresholds
+!   ztest = 0
+!   if ( z(ktrop) .gt. tropz_max ) then
+!      where ( z.le.tropz_max )
+!         ztest = 1
+!      end where
+!      do  k = nlev, 1, -1
+!         if (ztest(k) .eq. 1) then
+!            ktrop = k
+!            exit
+!         end if
+!      end do
+!   else if ( z(ktrop) .lt. tropz_min ) then
+!      where ( z.ge.tropz_min )
+!         ztest = 1
+!      end where
+!      do  k = 1, nlev
+!         if (ztest(k) .eq. 1) then
+!            ktrop = k
+!            exit
+!         end if
+!      end do
+!   end if
+
+   if (present(tropt)) tropt = t(ktrop)
+   if (present(tropp)) tropp = p(ktrop)
+   if (present(tropk)) tropk = ktrop
+
+end subroutine da_trop_wmo
+

var/da/da_radiance/da_allocate_rad_iv.inc

@@ -10,7 +10,7 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
 
    integer           ,  intent (in)    :: i, nchan
    type (iv_type)    ,  intent (inout) :: iv
-
+   integer :: n, ix, iy
    call da_trace_entry("da_allocate_rad_iv")
 
       allocate (iv%instid(i)%info%date_char(iv%instid(i)%num_rad))
@@ -104,7 +104,20 @@ subroutine da_allocate_rad_iv (i, nchan, iv)
       allocate (iv%instid(i)%satazi(iv%instid(i)%num_rad))
       allocate (iv%instid(i)%solzen(iv%instid(i)%num_rad))
       allocate (iv%instid(i)%solazi(iv%instid(i)%num_rad))
+      allocate (iv%instid(i)%tropt(iv%instid(i)%num_rad))	  
       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
+         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))
+         end do
+         end do
+      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_crtm.f90

@@ -32,7 +32,7 @@ module da_crtm
       biasprep, qc_rad,missing_r,rtminit_sensor,rtminit_nsensor, filename_len, &
       use_error_factor_rad,read_biascoef, analysis_date,time_window_max, &
       time_window_min,num_fgat_time,rtminit_platform, print_detail_rad, &
-      rtminit_satid, global,kms,kme,ims,ime,jms,jme,kts,kte,use_clddet, &
+      rtminit_satid, global,kms,kme,ims,ime,jms,jme,kts,kte,use_clddet, use_clddet_zz,&
       use_crtm_kmatrix, use_varbc, freeze_varbc, use_pseudo_rad, &
       use_antcorr, time_slots, use_satcv, use_simulated_rad, simulated_rad_io, &
       simulated_rad_ngrid, interp_option, use_mspps_emis, use_mspps_ts, calc_weightfunc, &
@@ -41,6 +41,7 @@ module da_crtm
       cloud_cv_options
    use da_interpolation, only : da_interp_lin_2d_partial,da_interp_lin_2d_adj_partial, &
       da_interp_2d_partial
+   use da_physics, only: da_trop_wmo
    use module_dm, only : wrf_dm_sum_real, wrf_dm_sum_reals
    use da_radiance1, only : da_biasprep,da_detsurtyp,da_biascorr, &
        da_biasprep,da_cld_eff_radius, da_mspps_emis, da_mspps_ts