Add ensemble component for RRFS-SD DA

 Changes to be committed:
	modified:   src/gsi/control2state.f90
        modified:   src/gsi/cplr_get_fv3_regional_ensperts.f90
	modified:   src/gsi/ensctl2model.f90
	modified:   src/gsi/ensctl2model_ad.f90
	modified:   src/gsi/ensctl2state.f90
	modified:   src/gsi/ensctl2state_ad.f90
	modified:   src/gsi/gsi_rfv3io_mod.f90
	modified:   src/gsi/gsimod.F90
	modified:   src/gsi/obsmod.F90
	modified:   src/gsi/wrf_vars_mod.f90

src/gsi/control2state.f90

@@ -35,6 +35,7 @@ module control2state_mod
 use gridmod, only: lat2,lon2,nsig,nlat,nlon            
 use chemmod, only: laeroana_fv3cmaq, naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,icvt_cmaq_fv3
 use mpeu_util, only: getindex
+use wrf_vars_mod, only : sdp_exist
 
 implicit none
 
@@ -283,6 +284,15 @@ subroutine control2state(xhat,sval,bval)
       enddo
    end if
 
+! Add smoke 
+   if (sdp_exist) then
+   call gsi_bundlegetpointer (sval(jj),'smoke',sv_rank3,istatus)
+   call gsi_bundlegetvar     (wbundle, 'smoke',sv_rank3,istatus)
+   call gsi_bundlegetpointer (sval(jj),'dust' ,sv_rank3,istatus)
+   call gsi_bundlegetvar     (wbundle, 'dust' ,sv_rank3,istatus)
+   call gsi_bundlegetpointer (sval(jj),'coarsepm',sv_rank3,istatus)
+   call gsi_bundlegetvar     (wbundle, 'coarsepm',sv_rank3,istatus)
+   end if
    call gsi_bundlegetpointer (sval(jj),'ps'  ,sv_ps,  istatus)
    call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps,  istatus )
 
@@ -732,6 +742,16 @@ subroutine control2state_ad(rval,bval,grad)
       enddo
    end if
 
+! Add smoke 
+   if (sdp_exist) then
+   call gsi_bundlegetpointer (rval(jj),'smoke',rv_rank3,istatus)
+   call gsi_bundleputvar     (wbundle, 'smoke',rv_rank3,istatus)
+   call gsi_bundlegetpointer (rval(jj),'dust' ,rv_rank3,istatus)
+   call gsi_bundleputvar     (wbundle, 'dust' ,rv_rank3,istatus)
+   call gsi_bundlegetpointer (rval(jj),'coarsepm',rv_rank3,istatus)
+   call gsi_bundleputvar     (wbundle, 'coarsepm',rv_rank3,istatus)
+   end if
+
 !  Calculate sensible temperature
    if(do_tv_to_tsen) call tv_to_tsen_ad(cv_t,rv_q,rv_tsen)
 

src/gsi/ensctl2model.f90

@@ -44,6 +44,7 @@ subroutine ensctl2model(xhat,mval,eval)
 use mod_strong, only: tlnmc_option
 use timermod, only: timer_ini,timer_fnl
 use hybrid_ensemble_parameters,only: naensgrp
+use wrf_vars_mod, only : sdp_exist
 implicit none
 
 ! Declare passed variables
@@ -212,7 +213,28 @@ subroutine ensctl2model(xhat,mval,eval)
           call gsi_bundlegetvar     (wbundle_c, clouds(ic),sv_rank3,istatus)
       endif
    enddo
-
+! add smoke 
+   if (sdp_exist) then
+     print*,"SMOKE_ensctl2model.f90"
+     call gsi_bundlegetpointer (eval(jj),'smoke',sv_rank3,istatus)
+     if(istatus/=0) then
+       write(6,*) trim(myname), ': trouble_get_pointer: smoke '
+     else
+       call gsi_bundlegetvar     (wbundle_c, 'smoke',sv_rank3,istatus)
+     end if
+     call gsi_bundlegetpointer (eval(jj),'dust',sv_rank3,istatus)
+     if(istatus/=0) then
+       write(6,*) trim(myname), ': trouble_get_pointer: dust '
+     else
+       call gsi_bundlegetvar     (wbundle_c, 'dust',sv_rank3,istatus)
+     end if
+     call gsi_bundlegetpointer (eval(jj),'coarsepm',sv_rank3,istatus)
+     if(istatus/=0) then
+       write(6,*) trim(myname), ': trouble_get_pointer: coarsepm '
+     else
+       call gsi_bundlegetvar     (wbundle_c, 'coarsepm',sv_rank3,istatus)
+     end if
+  end if
 !  Add contribution from static B, if necessary
    call self_add(eval(jj),mval)
 

src/gsi/ensctl2model_ad.f90

@@ -43,6 +43,7 @@ subroutine ensctl2model_ad(eval,mval,grad)
 use mod_strong, only: tlnmc_option
 use timermod, only: timer_ini,timer_fnl
 use hybrid_ensemble_parameters,only: naensgrp
+use wrf_vars_mod, only : sdp_exist
 implicit none
 
 ! Declare passed variables
@@ -190,6 +191,29 @@ subroutine ensctl2model_ad(eval,mval,grad)
           call gsi_bundleputvar     (wbundle_c, clouds(ic),rv_rank3,istatus)
       endif
    enddo
+
+! add smoke 
+  if (sdp_exist) then
+  print*,"SMOKE_ensctl2model_ad.f90"
+  call gsi_bundlegetpointer (eval(jj), 'smoke',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: smoke '
+  else
+    call gsi_bundleputvar     (wbundle_c,'smoke',rv_rank3,istatus)
+  endif
+  call gsi_bundlegetpointer (eval(jj), 'dust',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: dust '
+  else
+    call gsi_bundleputvar     (wbundle_c,'dust',rv_rank3,istatus)
+  endif
+  call gsi_bundlegetpointer (eval(jj), 'coarsepm',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: coarsepm '
+  else
+    call gsi_bundleputvar     (wbundle_c,'coarsepm',rv_rank3,istatus)
+  endif
+  end if
 !  Convert RHS calculations for u,v to st/vp
    if (do_getuv) then
       if(uv_hyb_ens) then

src/gsi/ensctl2state.f90

@@ -47,6 +47,8 @@ subroutine ensctl2state(xhat,mval,eval)
 use cwhydromod, only: cw2hydro_tl_hwrf
 use timermod, only: timer_ini,timer_fnl
 use gridmod, only: nems_nmmb_regional
+use wrf_vars_mod, only : sdp_exist
+
 implicit none
 
 ! Declare passed variables
@@ -236,6 +238,28 @@ subroutine ensctl2state(xhat,mval,eval)
       enddo
    endif
 
+! add smoke 
+  if (sdp_exist) then
+  print*,"SMOKE_ensctl2state.f90"
+  call gsi_bundlegetpointer (eval(jj),'smoke',sv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: smoke '
+  else
+   call gsi_bundlegetvar     (wbundle_c, 'smoke',sv_rank3,istatus)
+  end if
+  call gsi_bundlegetpointer (eval(jj),'dust',sv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: dust '
+  else
+   call gsi_bundlegetvar     (wbundle_c, 'dust',sv_rank3,istatus)
+  end if
+  call gsi_bundlegetpointer (eval(jj),'coarsepm',sv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: coarsepm '
+  else
+   call gsi_bundlegetvar     (wbundle_c, 'coarsepm',sv_rank3,istatus)
+  end if
+  end if
 !$omp section
 
 !  Get pointers to required state variables

src/gsi/ensctl2state_ad.f90

@@ -46,6 +46,8 @@ subroutine ensctl2state_ad(eval,mval,grad)
 use cwhydromod, only: cw2hydro_ad_hwrf
 use timermod, only: timer_ini,timer_fnl
 use gridmod, only: nems_nmmb_regional
+use wrf_vars_mod, only : sdp_exist
+
 implicit none
 
 ! Declare passed variables
@@ -245,6 +247,28 @@ subroutine ensctl2state_ad(eval,mval,grad)
       enddo
    endif
 
+! add smoke 
+  if (sdp_exist) then
+  print*,"SMOKE_ensctl2state_ad.f90"
+  call gsi_bundlegetpointer (eval(jj), 'smoke',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: smoke '
+  else
+    call gsi_bundleputvar     (wbundle_c,'smoke',rv_rank3,istatus)
+  endif
+  call gsi_bundlegetpointer (eval(jj), 'dust',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: dust '
+  else
+    call gsi_bundleputvar     (wbundle_c,'dust',rv_rank3,istatus)
+  endif
+  call gsi_bundlegetpointer (eval(jj), 'coarsepm',rv_rank3,istatus)
+  if(istatus/=0) then
+    write(6,*) trim(myname), ': trouble_get_pointer: coarsepm '
+  else
+    call gsi_bundleputvar     (wbundle_c,'coarsepm',rv_rank3,istatus)
+  endif
+  end if
 !  Calculate sensible temperature
    if(do_q_copy) then
       call gsi_bundleputvar (wbundle_c, 'q', rv_q, istatus )

src/gsi/gsi_rfv3io_mod.f90

@@ -1014,22 +1014,22 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
        endif
       end do 
       if (iuv /= 2.or. ndynvario3d<=0.or.ntracerio3d<=0 ) then
-        write(6,*)"the set up for met variable is not as expected, abort"
+        write(6,*)"The set up for met variable is not as expected, abort"
         call stop2(222)
       endif
       if ( if_model_dbz .and. if_model_fed ) then
          if( nphyvario3d<=1 ) then
-            write(6,*)"the set up for met variable (dbz and fed in phyvar) is not as expected,abort"
+            write(6,*)"The set up for met variable (dbz and fed in phyvar) is not as expected,abort"
             call stop2(223)
          end if
       elseif ( if_model_dbz ) then 
          if( nphyvario3d<=0 ) then
-            write(6,*)"the set up for met variable (dbz in phyvar) is not as expected, abort"
+            write(6,*)"The set up for met variable (dbz in phyvar) is not as expected, abort"
             call stop2(223)
          end if
      elseif ( if_model_fed ) then
          if( nphyvario3d<=0 ) then
-            write(6,*)"the set up for met variable (fed in phyvar) is not as expected, abort"
+            write(6,*)"The set up for met variable (fed in phyvar) is not as expected, abort"
             call stop2(223)
          end if
       endif
@@ -2862,7 +2862,7 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin)
 end subroutine gsi_fv3ncdf_readuv_v1
 
 subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
-           delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed,iope)
+           delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed,smoke,dust,coarsepm,iope)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    gsi_fv3ncdf_read_ens_parallel_over_ens    
@@ -2897,14 +2897,15 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
     use mod_fv3_lola, only: fv3_h_to_ll
     use gsi_bundlemod, only: gsi_bundle
     use general_sub2grid_mod, only: sub2grid_info,general_grid2sub
+    use gsi_io, only: verbose
 
     implicit none
     character(*),intent(in):: filenamein
     type (type_fv3regfilenameg),intent(in) ::fv3filenamegin
     integer(i_kind)   ,intent(in   ) :: iope
     real(r_kind),allocatable,dimension(:,:):: uu2d, uu2d_tmp
     real(r_kind),dimension(nlat,nlon,nsig):: hwork
-    real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed
+    real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed,smoke,dust,coarsepm
     character(len=max_varname_length) :: varname
     character(len=max_varname_length) :: name
     character(len=max_filename_length), allocatable,dimension(:) :: varname_files
@@ -2921,43 +2922,78 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
     integer(i_kind),allocatable :: gfile_loc_layout(:)
     character(len=180)  :: filename_layout
 
+    logical print_verbose
+
+
+    print_verbose = .false.
+    if(verbose)print_verbose=.true.
+
     mm1=mype+1
     nloncase=nlon
     nlatcase=nlat
     nxcase=nx
     nycase=ny
     kbgn=1
     kend=nsig
-
+    !if (print_verbose) then
+      write(6,*)"G1_READ_ENS_Present_delp ",present(delp)
+      write(6,*)"G1_READ_ENS_Present_tsen ",present(tsen)
+      write(6,*)"G1_READ_ENS_Present_w    ",present(w)
+      write(6,*)"G2_READ_ENS_Present_q    ",present(q)
+      write(6,*)"G2_READ_ENS_Present_oz   ",present(oz)
+      write(6,*)"G2_READ_ENS_Present_ql   ",present(ql)
+      write(6,*)"G2_READ_ENS_Present_qr   ",present(qr)
+      write(6,*)"G2_READ_ENS_Present_qs   ",present(qs)
+      write(6,*)"G2_READ_ENS_Present_qi   ",present(qi)
+      write(6,*)"G2_READ_ENS_Present_qg   ",present(qg)
+      write(6,*)"G3_READ_ENS_Present_dbz  ",present(dbz)
+      write(6,*)"G4_READ_ENS_Present_fed  ",present(fed)
+      write(6,*)"G5_READ_ENS_Present_smoke",present(smoke)
+      write(6,*)"G5_READ_ENS_Present_dust ",present(dust)
+      write(6,*)"G5_READ_ENS_Present_copm ",present(coarsepm)
+    !end if
     if( mype == iope )then
        allocate(uu2d(nxcase,nycase))
        if( present(delp).or.present(tsen).or.present(w) )then  ! dynvars
+          print*,"READ_ENS_Size_T_p_W ", present(delp),present(tsen),present(w)
           if( present(w) )then
              allocate(varname_files(3))
              varname_files = (/'T   ','delp','W   '/)
           else
              allocate(varname_files(2))
              varname_files = (/'T   ','delp'/)
           end if
+          print*,"READ_ENS_Size_T_p_W ",size(varname_files),varname_files
        end if
        if( present(q).or.present(ql).or.present(qr) )then ! tracers
+       print*,"READ_ENS_Size_q_o3_hydro ", present(q),present(ql),present(qr)
           if(present(qr))then
              allocate(varname_files(7))
              varname_files = (/'sphum  ','o3mr   ','liq_wat','ice_wat','rainwat','snowwat','graupel'/)
           else
              allocate(varname_files(2))
              varname_files = (/'sphum',' o3mr'/)
           end if
+          print*,"READ_ENS_Size_q_o3_hydro ",size(varname_files),varname_files
        end if
+       print*,"READ_ENS_Size_dbz ", present(dbz)
        if( present(dbz) )then            ! phyvars: dbz
           allocate(varname_files(1))
           varname_files = (/'ref_f3d'/)
+          print*,"READ_ENS_Size_dbz ",size(varname_files),varname_files
        end if
+       print*,"READ_ENS_Size_fed ", present(fed)
        if( present(fed) )then            ! phyvars: fed 
           allocate(varname_files(1))
           varname_files = (/'flash_extent_density'/)
+          print*,"READ_ENS_Size_fed ",size(varname_files),varname_files
+       end if
+       print*,"READ_ENS_Size_sdp ", present(smoke),present(dust),present(coarsepm)
+       if( present(smoke).or.present(dust).or.present(coarsepm) )then            ! tracers
+          allocate(varname_files(3))
+          varname_files = (/'smoke','dust','coarsepm'/)
+          print*,"READ_ENS_Size_sdp ",size(varname_files),varname_files
        end if
-
 
        if(fv3_io_layout_y > 1) then
           allocate(gfile_loc_layout(0:fv3_io_layout_y-1))
@@ -2978,7 +3014,14 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
              call stop2(333)
           endif
        endif
+       print*,"READ_ENS_size_varname_files ",size(varname_files),varname_files 
+       if (size(varname_files) == 0) then 
+          !if (print_verbose)
+           write(6,*)"Please check your subroutine call, No variale is setup to read !!!! "
+          return 
+       end if
        do ivar = 1, size(varname_files)
+          print*,"READ_ENS: ivar= ",ivar,trim(varname_files(ivar))
           do ilevtot=kbgn,kend
              ilev=ilevtot
              nz=nsig
@@ -3069,7 +3112,15 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
           if( present(fed) )then            ! phyvars: fed 
             fed = hwork
           end if
-
+          if( present(smoke) )then   
+            smoke = hwork
+          end if
+          if( present(dust) )then                 
+            dust = hwork
+          end if
+          if( present(coarsepm) )then                 
+            coarsepm = hwork
+          end if
        end do
 
        if(fv3_io_layout_y > 1) then
@@ -3080,7 +3131,14 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, &
         else
           iret=nf90_close(gfile_loc)
        endif
-
+       !if(ALLOCATED(uu2d))then 
+       !  deallocate (uu2d)
+       !  print*,"ASSOCIATED_uu2d"
+       !end if
+       !if(ALLOCATED(varname_files))then
+       !  deallocate(varname_files)
+       ! print*,"ASSOCIATED_varname_files"
+       !end if
        deallocate (uu2d,varname_files)
     end if
 
@@ -5152,6 +5210,7 @@ subroutine gsi_copy_bundle(bundi,bundo)
     call gsi_bundleinquire(bundo,'shortnames::3d',target_name_vars3d,istatus)
     call gsi_bundleinquire(bundo,'shortnames::2d',target_name_vars2d,istatus)
     do ivar=1,src_nc3d
+    print*,"COPY_BUND_ivar= ",ivar,trim(src_name_vars3d(ivar))
       varname=trim(src_name_vars3d(ivar))
       do jvar=1,target_nc3d
         if(index(target_name_vars3d(jvar),varname) > 0)  then