Add GeoVaLS for radiance and ozone data types. (NOAA-EMC#4)

Co-authored-by: Emily Liu <[email protected]>

src/gsi/crtm_interface.f90

@@ -85,6 +85,7 @@ module crtm_interface
 public destroy_crtm         ! Subroutine destroys initialization for crtm
 public sensorindex
 public surface
+public atmosphere           !emily 
 public isatid               ! = 1  index of satellite id
 public itime                ! = 2  index of analysis relative obs time
 public ilon                 ! = 3  index of grid relative obs location (x)
@@ -125,6 +126,10 @@ module crtm_interface
 public idtw                 ! = 35/37 index of d(Tw)
 public idtc                 ! = 36/38 index of d(Tc)
 public itz_tr               ! = 37/39 index of d(Tz)/d(Tr)
+!>>emily
+public n_clouds_fwd_wk
+public n_absorbers
+!<<emily
 
 ! For TMI and GMI
 public iedge_log            ! = 32  ! index, if obs is to be obleted beause of locating near scan edges.
@@ -203,6 +208,7 @@ module crtm_interface
   logical        ,save :: lprecip_wk 
   logical        ,save :: mixed_use
   integer(i_kind), parameter :: min_n_absorbers = 2
+  integer(i_kind) :: n_absorbers  !emily
 
   integer(i_kind),save :: iedge_log
   integer(i_kind),save :: ilzen_ang2,ilazi_ang2,iscan_ang2,iszen_ang2,isazi_ang2
@@ -355,7 +361,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo
 ! ...all "additional absorber" variables
   integer(i_kind) :: j,icount
   integer(i_kind) :: ig
-  integer(i_kind) :: n_absorbers
+! integer(i_kind) :: n_absorbers  !orig
   logical quiet
   logical print_verbose
 

src/gsi/setupoz.f90

@@ -113,7 +113,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
   use state_vectors, only: svars3d, levels
 
   use constants, only : zero,half,one,two,tiny_r_kind
-  use constants, only : rozcon,cg_term,wgtlim,h300,r10
+  use constants, only : rozcon,cg_term,wgtlim,h300,r10,r1000,constoz 
 
   use m_obsdiagNode, only : obs_diag
   use m_obsdiagNode, only : obs_diags
@@ -201,6 +201,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
   real(r_kind),dimension(nlevs):: pobs,gross,tnoise
   real(r_kind),dimension(nreal+nlevs,nobs):: data
   real(r_kind),dimension(nsig+1)::prsitmp
+  real(r_kind),dimension(nsig)::ozgestmp   ! GeoVaLs for JEDI/UFO
   real(r_single),dimension(nlevs):: pob4,grs4,err4
   real(r_single),dimension(ireal,nobs):: diagbuf
   real(r_single),allocatable,dimension(:,:,:)::rdiagbuf
@@ -409,7 +410,13 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
            ozp_omi(nloz_omi) = prsitmp(1)
            call grdcrd1(ozp_omi(nloz_omi),prsitmp,nsig+1,-1)
         end if
-
+
+       ! GeoVaLs for JEDI/UFO
+       call tintrp2a1(ges_oz,ozgestmp,dlat,dlon,dtime,hrdifsig,&
+                      nsig,mype,nfldsig)
+       call tintrp2a1(ges_prsi, prsitmp,dlat,dlon,dtime,hrdifsig,&
+                      nsig+1,mype,nfldsig)
+
         if (obstype /= 'omieff' .and. obstype /= 'tomseff') then
            call intrp3oz1(ges_oz,ozges,dlat,dlon,ozp,dtime,&
                 nlevs,mype,doz_dz)
@@ -593,6 +600,9 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,&
                    call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind)
                    call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
                  endif
+                 ! GeoVaLs for JEDI/UFO
+                 call nc_diag_data2d("mole_fraction_of_ozone_in_air", sngl(ozgestmp*constoz)) !emily 
+                 call nc_diag_data2d("air_pressure_levels",sngl(prsitmp*r1000))
               endif
            endif
 

src/gsi/setuprad.f90

@@ -280,9 +280,9 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&
       finish_sst_retrieval,spline_cub
   use m_dtime, only: dtime_setup, dtime_check
   use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex,surface, &
-      itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang,iscan_pos,iszen_ang,isazi_ang, &
+      atmosphere,itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang,iscan_pos,iszen_ang,isazi_ang, &
       ifrac_sea,ifrac_lnd,ifrac_ice,ifrac_sno,itsavg, &
-      izz,idomsfc,isfcr,iff10,ilone,ilate, &
+      izz,idomsfc,isfcr,iff10,ilone,ilate, n_clouds_fwd_wk, n_absorbers, &
       isst_hires,isst_navy,idata_type,iclr_sky,itref,idtw,idtc,itz_tr
   use qcmod, only: qc_ssmi,qc_geocsr,qc_ssu,qc_avhrr,qc_goesimg,qc_msu,qc_irsnd,qc_amsua,qc_mhs,qc_atms
   use crtm_interface, only: ilzen_ang2,iscan_ang2,iszen_ang2,isazi_ang2
@@ -2487,6 +2487,8 @@ subroutine contents_netcdf_diag_(odiags,idv,iob)
   real(r_single),parameter::  missing = -9.99e9_r_single
   integer(i_kind),parameter:: imissing = -999999
   real(r_kind),dimension(:),allocatable :: predbias_angord
+  character(128) :: fieldname
+  integer(i_kind) :: iabsorb, icloud
 
   if (adp_anglebc) then
     allocate(predbias_angord(angord) )
@@ -2662,7 +2664,32 @@ subroutine contents_netcdf_diag_(odiags,idv,iob)
                        call nc_diag_data2d("BCPred_angord",   sngl(predbias_angord)                                )
                     endif
                  end if
-
+                 ! GeoVaLs for JEDI/UFO
+                 ! Get GeoVaLs for surface 
+                 call nc_diag_metadata("Vegetation_Type",    sngl(surface(1)%vegetation_type)                   )
+                 call nc_diag_metadata("Lai",                sngl(surface(1)%lai)                               )
+                 call nc_diag_metadata("Soil_Type",          sngl(surface(1)%soil_type)                         )
+
+                 call nc_diag_metadata("Sfc_Wind_Direction", sngl(surface(1)%wind_direction)                    )
+                 call nc_diag_metadata("Sfc_Height",         sngl(zsges    )                                    )   ! do we need this for geoval? I think we do not
+
+                 ! Get GeoVaLs for atmosphere
+                 call nc_diag_data2d("air_temperature",      sngl(atmosphere(1)%temperature)                    )   ! K 
+                 call nc_diag_data2d("air_pressure",         sngl(atmosphere(1)%pressure*r100)                  )
+                 call nc_diag_data2d("air_pressure_levels",  sngl(atmosphere(1)%level_pressure*r100)            )
+
+                 ! Get GeoVaLs for atmospheric absorbers
+                 do iabsorb = 1, n_absorbers
+                   write (fieldname, "(A,I0.2)") "atmosphere_absorber_", atmosphere(1)%absorber_id(iabsorb)
+                   call nc_diag_data2d(trim(fieldname),      sngl(atmosphere(1)%absorber(:,iabsorb))            )   ! check %absorber_units
+                 enddo
+                 ! Get GeoVaLs for hydrometeors 
+                 do icloud = 1, n_clouds_fwd_wk
+                   write (fieldname, "(A,I0.2)") "atmosphere_mass_content_of_cloud_", atmosphere(1)%Cloud(icloud)%Type
+                   call nc_diag_data2d(trim(fieldname),      sngl(atmosphere(1)%Cloud(icloud)%Water_Content)    )
+                   write (fieldname, "(A,I0.2)") "effective_radius_of_cloud_particle_", atmosphere(1)%Cloud(icloud)%Type
+                   call nc_diag_data2d(trim(fieldname),      sngl(atmosphere(1)%Cloud(icloud)%Effective_Radius) )
+                 enddo
               enddo
 !  if (adp_anglebc) then
   if (.true.) then