Merge pull request #216 from GEOS-ESM/feature/vbuchard/update_ana_aod…

…_code

vb:update ana_aod.F code to handle both AOD and LAOD obs at the same …

src/Applications/GAAS_App/ana.rc

@@ -51,6 +51,8 @@ do_you_want_to_skip_PSAS: no # yes or no
 
 alpha_for_bias_estimation: 0.010
 
+eps_for_log_transform_aod_on_input:  0.01
+
 range_of_wavenumbers_to_analyze_in_nm: 470 870
 
 # superob_IMxJM: 576 361 # superob dimension (if <0, same as input grid)

src/Applications/GAAS_App/ana_aod.F

@@ -67,9 +67,10 @@ program ana_aod
 !     ---------------
       integer :: im, jm, km, nch
       integer :: im_so, jm_so    ! dimension of superob space
-      integer :: rc, nobs, nobs_good, nODS, i
+      integer :: rc, nobs, nobs_good, nODS, i, naod, nlaod
       integer :: nymd=0, nhms=0, nymd_b, nhms_b
       integer :: myKTobs
+      integer :: kt_list(2)
 
       real*8  :: alpha         ! bias update coefficient
       real*8  :: chAOD(2)      ! Range of AOD waenumbers (nm) to analyze
@@ -255,34 +256,42 @@ program ana_aod
       if ( rc /= 0 ) then
          call die ( myname, 'could not read ODS files')
       end if
-      myKTobs=ktAOD
+      ! can read AOD and log AOD obs
+      kt_list = (/ktAOD,ktLogAOD/)
       call ods_select( ods_, ods_%data%nobs, nobs, rc, 
-     &                 qcexcl=0, kt=myKTobs, odss = ods ) ! keep what we need
+     &                 qcexcl=0, kt_list=kt_list, odss = ods ) ! keep what we need
       if ( rc .eq. 0 ) then
-         nobs =  ods%data%nobs
-         if ( nobs==0 ) then
-            call ods_clean(ods,rc)
-            ! try reading log(AOD)
-            print *, myname//': could not find obs with kt = ', myKTobs
-            myKTobs=ktLogAOD
-            print *, myname//': trying kt = ', myKTobs, ' ... '
-            call ods_select( ods_, ods_%data%nobs, nobs, rc, 
-     &                       qcexcl=0, kt=myKTobs, odss = ods ) ! keep what we need
+         nobs  = ods%data%nobs
+         naod  = count(ods%data%kt==ktAOD)
+         nlaod = count(ods%data%kt==ktLogAOD)
+         ! if Log AOD comes in as input, convert obs into AOD
+         if (eps<=0.0) then 
+             call die ( myname, 'eps not defined for conversion of log(AOD) obs into AOD')
+         else
+           if (nlaod >0) then 
+               where (ods%data%kt == ktLogAOD) 
+                      ods%data%obs = max(eps,exp(ods%data%obs)-eps)
+                      ods%data%kt  = ktAOD            
+               endwhere
+           endif
          endif
-         call ods_clean(ods_,rc) ! get rid of what we don't need
+         myKTobs = ktAOD  ! all obs are now supposed to be AOD 
+         call ods_clean(ods_,rc) ! get rid of what we do not need
          print *, myname//': read ODS files'
          print *, myname//': nobs = ', nobs
+         print *, myname//': naod = ', naod, ' and nlaod = ', nlaod
          print *
       else
          call die ( myname, 'could not perform ODS select')
       end if
 
+
 !     Prints out input ODS summary
 !     ----------------------------
       call fix_odsmeta(ods)
       call ODS_Tally ( 6, ods, nobs, rc )
 
-!     Call observer for this date/time
+!     Call observer for this date/time (the convert option becomes obsolete now)
 !     --------------------------------
                                                                                call zeit_ci ( 'Observer' )
       call Observer ( nymd, nhms, y_f, nobs, ods, nobs_good, convertOb2AOD=(myKTobs==ktLogAOD) )
@@ -782,6 +791,17 @@ subroutine Init_ ( expid, aerFile, pre_ods, mODS,
           end if
       end if
 
+!     Try reading eps, ok if not present but code will only handle
+!     AOD in the input file; code will fail is log AOD found and eps
+!     is not defined.
+      call i90_label ( 'eps_for_log_transform_aod_on_input:', ier )
+      if (ier==0) then
+         eps = i90_gfloat ( ier )
+      else
+         print *, trim(myname), ': WARNING, eps not defined early enough'
+         print *, trim(myname), ': WARNING, code can only handle AOD in input(ODS)'
+      endif
+
 !     Bias coefficient
 !     ----------------
       if ( do_SBC ) then

src/Applications/GAAS_App/m_obs.F90

@@ -67,13 +67,14 @@ MODULE  m_obs
 !  14dec2000  da Silva  Increased mBoxes from 32 to 128.
 !  05jun2002  Dee       Convert heights to layer-mean virtual temperatures
 !  09feb2010  da Silva  Adapted from old GEOS-4 observer for AOD.
+!  03oct2022  Todling   initilize eps in preamble.
 !
 !EOP
 !-------------------------------------------------------------------------
 
     logical :: do_sqc = .true.      ! whether or not to statistical qc
     logical :: log_transf = .false. ! whether to log-transform AOD
-    real    :: eps                  ! eps for log-transform
+    real    :: eps = -obs_missing   ! eps for log-transform
     logical :: do_dupelim = .true.  ! duplicate elimination?
 
 !   Resource file