Merge pull request #73 from imoradi/CRTMV3_active_sensor

Crtmv3 active sensor
JCSDA · Sep 22, 2023 · 0ee8a41 · 0ee8a41
2 parents 6053b42 + b42fe2e
commit 0ee8a41
Show file tree

Hide file tree

Showing 51 changed files with 4,243 additions and 3,836 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,2 +1,10 @@
 .DS_Store
 *.o
+*.tgz
+.ipynb*
+build
+fix/
+**/Build/
+build.bash
+rebuild.bash
+.gitignore
diff --git a/configuration/ifort.setup b/configuration/ifort.setup
@@ -5,9 +5,6 @@
 
 hostname=$(hostname)
 export FC="ifort"
-export NC4_DIR="/opt/netcdf4/4.6.2-intel-18.0.3"
-export HDF_DIR="/opt/hdf5/1.8.21-intel-18.0.3"
-#export HDF4_DIR="/opt/hdf4/4.2.14-intel-18.0.3"
 
 export OMP_NUM_THREADS=4  #set here for testing purposes, probably want to override this
 
@@ -25,11 +22,11 @@ HDF5DISC="/usr/local/other/hdf5/1.13.0/intel-19.1.0.16_impi-20.0.0.166/"
 
 if [[ "$hostname" == "s4-submit.ssec.wisc.edu" ]] ; then #S4
   export LIBS="\
-  -L $NC4_DIR/lib -lnetcdf -lnetcdff \
-  -L $HDF_DIR/lib -lhdf5 \
+  -L ${NC4_LIB} -lnetcdf -lnetcdff \
+  -L ${HDF5_LIB} -lhdf5 \
   -g -fpp -traceback -check bounds \
-  -I $NC4_DIR/include \
-  -I $HDF_DIR/include "
+  -I ${NC4_INC} \
+  -I ${HDF5_INC} "
 elif [[ "$hostname" == "discover"* ]] ; then # NCCS Discover
   export LIBS="\
   -L $NETCDF/lib -lnetcdf -lnetcdff \
@@ -40,4 +37,3 @@ elif [[ "$hostname" == "discover"* ]] ; then # NCCS Discover
 else
   export LIBS=""
 fi
-
diff --git a/src/AtmAbsorption/ODPS/ODPS_CoordinateMapping.f90 b/src/AtmAbsorption/ODPS/ODPS_CoordinateMapping.f90
@@ -43,6 +43,7 @@ MODULE ODPS_CoordinateMapping
   PUBLIC :: Interpolate_Profile_F1_TL
   PUBLIC :: Interpolate_Profile_F1_AD
   PUBLIC :: Compute_Interp_Index
+  PUBLIC :: Geopotential_Height
 
   ! Parameters used in the geopotential height calculation routines
   ! a factor used in the virtual temperature Tv calculation

diff --git a/src/AtmOptics/AOvar_Define.f90 b/src/AtmOptics/AOvar_Define.f90
@@ -9,6 +9,10 @@
 !       Written by:     Paul van Delst, 01-Jul-2013
 !                       [email protected]
 !                       
+!       Modified by:    Isaac Moradi,    [email protected]
+!                       14-Nov-2021
+!                       Added the backscat_coefficient
+!
 
 MODULE AOvar_Define
 
@@ -90,6 +94,7 @@ MODULE AOvar_Define
     REAL(fp), ALLOCATABLE :: optical_depth(:)
     REAL(fp), ALLOCATABLE :: bs(:)
     REAL(fp), ALLOCATABLE :: w(:)
+    REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:)
   END TYPE AOvar_type
 
 
@@ -134,6 +139,7 @@ ELEMENTAL SUBROUTINE AOvar_Create( &
     ALLOCATE( self%optical_depth(n_Layers), &
               self%bs(n_Layers), &
               self%w(n_Layers), &
+              self%Backscat_Coefficient(n_layers), &
               STAT = alloc_stat )
     IF ( alloc_stat /= 0 ) RETURN
 
@@ -164,6 +170,8 @@ SUBROUTINE AOvar_Inspect( self)
     WRITE(*,'(5(1x,es22.15,:))') self%bs
     WRITE(*,'(3x,"Single scatter albedo (w) :")')
     WRITE(*,'(5(1x,es22.15,:))') self%w
+    WRITE(*,'(3x,"Backscattering coefficient :")')
+    WRITE(*,'(5(1x,es22.15,:))') self%Backscat_Coefficient
   END SUBROUTINE AOvar_Inspect
 
 
@@ -471,6 +479,12 @@ FUNCTION AOvar_ReadFile( &
       msg = 'Error reading single scatter albedo - '//TRIM(io_msg)
       CALL Read_Cleanup(); RETURN
     END IF
+    ! ...Backscattering coefficient
+    READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) AOvar%Backscat_Coefficient
+    IF ( io_stat /= 0 ) THEN
+      msg = 'Error reading backscattering coefficient - '//TRIM(io_msg)
+      CALL Read_Cleanup(); RETURN
+    END IF
 
 
     ! Close the file
@@ -631,6 +645,12 @@ FUNCTION AOvar_WriteFile( &
       msg = 'Error writing single scatter albedo - '//TRIM(io_msg)
       CALL Write_Cleanup(); RETURN
     END IF
+    ! ...backscattering coefficient
+    WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) AOvar%Backscat_Coefficient
+    IF ( io_stat /= 0 ) THEN
+      msg = 'Error writing backscattering coefficient - '//TRIM(io_msg)
+      CALL Write_Cleanup(); RETURN
+    END IF
 
 
     ! Close the file
@@ -690,10 +710,12 @@ ELEMENTAL FUNCTION AOvar_Equal( x, y ) RESULT( is_equal )
     ! ...Dimensions
     IF ( (x%n_Layers /= y%n_Layers ) ) RETURN
     ! ...Data
-    IF (    (x%transmittance .EqualTo. y%transmittance ) .AND. &
-         ALL(x%optical_depth .EqualTo. y%optical_depth ) .AND. &
-         ALL(x%bs            .EqualTo. y%bs            ) .AND. &
-         ALL(x%w             .EqualTo. y%w             ) ) &
+    IF (    (x%transmittance         .EqualTo. y%transmittance )         .AND. &
+         ALL(x%optical_depth         .EqualTo. y%optical_depth )         .AND. &
+         ALL(x%bs                    .EqualTo. y%bs            )         .AND. &
+         ALL(x%w                     .EqualTo. y%w             )         .AND. &
+         ALL(x%Backscat_Coefficient .EqualTo. y%Backscat_Coefficient )       &
+         ) &
       is_equal = .TRUE.
   END FUNCTION AOvar_Equal
 

diff --git a/src/AtmOptics/CRTM_AtmOptics.f90 b/src/AtmOptics/CRTM_AtmOptics.f90
@@ -9,6 +9,12 @@
 !                       Paul van Delst, [email protected]
 !                       08-Jun-2005
 !
+!       Modified by:    Isaac Moradi,    [email protected]
+!                       14-Nov-2021
+!                       Added the backscat_coefficient as well as layer to
+!                       TOA transmittance with its AD and TL
+!
+
 
 MODULE CRTM_AtmOptics
 
@@ -209,7 +215,6 @@ SUBROUTINE CRTM_Compute_Transmittance( atmoptics, transmittance )
     transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k)))
   END SUBROUTINE CRTM_Compute_Transmittance
 
-
 !--------------------------------------------------------------------------------
 !:sdoc+:
 !
@@ -386,8 +391,9 @@ SUBROUTINE CRTM_AtmOptics_Combine( &
 
 
       ! Save the unmodified optical parameters
-      AOvar%Optical_Depth(k) = AtmOptics%Optical_Depth(k)
-      AOvar%bs(k)            = AtmOptics%Single_Scatter_Albedo(k)
+      AOvar%Optical_Depth(k)         = AtmOptics%Optical_Depth(k)
+      AOvar%bs(k)                    = AtmOptics%Single_Scatter_Albedo(k)
+      AOvar%Backscat_Coefficient(k) = AtmOptics%Backscat_Coefficient(k)
       ! ...Initialise scattering dependent terms
       AOvar%w(k) = ZERO
 
@@ -404,16 +410,15 @@ SUBROUTINE CRTM_AtmOptics_Combine( &
           ! ...Normalization requirement for energy conservation
           AtmOptics%Phase_Coefficient(0,i,k) = POINT_5
         END DO
-        AtmOptics%Delta_Truncation(k) = AtmOptics%Phase_Coefficient(AtmOptics%n_Legendre_Terms,1,k)
 
+        AtmOptics%Delta_Truncation(k) = AtmOptics%Phase_Coefficient(AtmOptics%n_Legendre_Terms,1,k)
 
         ! Redfine the total optical depth and single scattering
         ! albedo for the delta-function adjustment
         AtmOptics%Optical_Depth(k) = ( ONE - ( AtmOptics%Delta_Truncation(k) * AOvar%w(k) )) * &
                                      AtmOptics%Optical_Depth(k)
         AtmOptics%Single_Scatter_Albedo(k) = ( ONE - AtmOptics%Delta_Truncation(k) ) * AOvar%w(k)   / &
                                              ( ONE - ( AtmOptics%Delta_Truncation(k) * AOvar%w(k) ) )
-
       END IF Significant_Scattering
 
 

diff --git a/src/AtmOptics/CRTM_AtmOptics_Define.f90 b/src/AtmOptics/CRTM_AtmOptics_Define.f90
@@ -8,6 +8,10 @@
 !       Written by:     Paul van Delst, 11-Oct-2011
 !                       [email protected]
 !
+!       Modified by:    Isaac Moradi, 11-Nov-2021
+!                       [email protected]
+!                       Modifications to include back scattering information
+!
 
 MODULE CRTM_AtmOptics_Define
 
@@ -114,6 +118,7 @@ MODULE CRTM_AtmOptics_Define
     REAL(fp), ALLOCATABLE :: Optical_Depth(:)         ! K-Max
     REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:) ! K-Max
     REAL(fp), ALLOCATABLE :: Asymmetry_Factor(:)      ! K-Max
+    REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:)  ! K-Max
     REAL(fp), ALLOCATABLE :: Delta_Truncation(:)      ! K-Max
     REAL(fp), ALLOCATABLE :: Phase_Coefficient(:,:,:) ! 0:Ic-Max x Ip-Max x K-Max
   END TYPE CRTM_AtmOptics_type
@@ -302,6 +307,7 @@ PURE SUBROUTINE AtmOptics_Allocate(self,alloc_stat)
       ALLOCATE( self%Optical_Depth( n_Layers ), &
                 self%Single_Scatter_Albedo( n_Layers ), &
                 self%Asymmetry_Factor( n_Layers ), &
+                self%Backscat_Coefficient( n_Layers ), &
                 self%Delta_Truncation( n_Layers ), &
                 self%Phase_Coefficient( 0:n_Legendre_Terms, n_Phase_Elements, n_Layers ), &
                 STAT = alloc_stat )
@@ -347,6 +353,7 @@ ELEMENTAL SUBROUTINE CRTM_AtmOptics_Zero( self )
     self%Optical_Depth         = ZERO
     self%Single_Scatter_Albedo = ZERO
     self%Asymmetry_Factor      = ZERO
+    self%Backscat_Coefficient  = ZERO
     self%Delta_Truncation      = ZERO
     self%Phase_Coefficient     = ZERO
   END SUBROUTINE CRTM_AtmOptics_Zero
@@ -393,6 +400,8 @@ SUBROUTINE Scalar_Inspect(self)
     WRITE(*,'(5(1x,es22.15,:))') self%Single_Scatter_Albedo(1:self%n_Layers)
     WRITE(*,'(3x,"Asymmetry_Factor :")')
     WRITE(*,'(5(1x,es22.15,:))') self%Asymmetry_Factor(1:self%n_Layers)
+    WRITE(*,'(3x,"Backscat_Coefficient Backscattering coefficient :")')
+    WRITE(*,'(5(1x,es22.15,:))') self%Backscat_Coefficient(1:self%n_Layers)
     WRITE(*,'(3x,"Delta_Truncation :")')
     WRITE(*,'(5(1x,es22.15,:))') self%Delta_Truncation(1:self%n_Layers)
     WRITE(*,'(3x,"Phase_Coefficient Legendre polynomial coefficients :")')
@@ -628,6 +637,10 @@ ELEMENTAL FUNCTION CRTM_AtmOptics_Compare( &
                       x%Asymmetry_Factor(1:k), &
                       y%Asymmetry_Factor(1:k), &
                       n                      ))) .OR. &
+         (.NOT. ALL(Compares_Within_Tolerance( &
+                      x%Backscat_Coefficient(1:k), &
+                      y%Backscat_Coefficient(1:k), &
+                      n                      ))) .OR. &
          (.NOT. ALL(Compares_Within_Tolerance( &
                       x%Delta_Truncation(1:k), &
                       y%Delta_Truncation(1:k), &
@@ -1022,6 +1035,7 @@ FUNCTION CRTM_AtmOptics_ReadFile( &
         AtmOptics(n)%Optical_Depth        , &
         AtmOptics(n)%Single_Scatter_Albedo, &
         AtmOptics(n)%Asymmetry_Factor     , &
+        AtmOptics(n)%Backscat_Coefficient , &
         AtmOptics(n)%Delta_Truncation
       IF ( io_stat /= 0 ) THEN
         msg = 'Error reading profile data '//TRIM(count_msg)//' from '//&
@@ -1277,6 +1291,7 @@ FUNCTION CRTM_AtmOptics_WriteFile( &
         AtmOptics(n)%Optical_Depth(1:AtmOptics(n)%n_Layers)        , &
         AtmOptics(n)%Single_Scatter_Albedo(1:AtmOptics(n)%n_Layers), &
         AtmOptics(n)%Asymmetry_Factor(1:AtmOptics(n)%n_Layers)     , &
+        AtmOptics(n)%Backscat_Coefficient(1:AtmOptics(n)%n_Layers)     , &
         AtmOptics(n)%Delta_Truncation(1:AtmOptics(n)%n_Layers)
       IF ( io_stat /= 0 ) THEN
         msg = 'Error writing profile data '//TRIM(count_msg)//' to '//&
@@ -1400,6 +1415,7 @@ ELEMENTAL FUNCTION CRTM_AtmOptics_Equal( x, y ) RESULT( is_equal )
       IF ( .NOT. (ALL(x%Optical_Depth(1:k)         .EqualTo. y%Optical_Depth(1:k)        ) .AND. &
                   ALL(x%Single_Scatter_Albedo(1:k) .EqualTo. y%Single_Scatter_Albedo(1:k)) .AND. &
                   ALL(x%Asymmetry_Factor(1:k)      .EqualTo. y%Asymmetry_Factor(1:k)     ) .AND. &
+                  ALL(x%Backscat_Coefficient(1:k)  .EqualTo. y%Backscat_Coefficient(1:k) ) .AND. &
                   ALL(x%Delta_Truncation(1:k)      .EqualTo. y%Delta_Truncation(1:k)     ) .AND. &
                   ALL(x%Phase_Coefficient(0:ic, 1:ip, 1:k) .EqualTo. &
                       y%Phase_Coefficient(0:ic, 1:ip, 1:k))) ) RETURN
@@ -1471,6 +1487,7 @@ ELEMENTAL FUNCTION CRTM_AtmOptics_Add( ao1, ao2 ) RESULT( aosum )
     aosum%Optical_Depth(1:k)               = aosum%Optical_Depth(1:k)               + ao2%Optical_Depth(1:k)
     aosum%Single_Scatter_Albedo(1:k)       = aosum%Single_Scatter_Albedo(1:k)       + ao2%Single_Scatter_Albedo(1:k)
     aosum%Asymmetry_Factor(1:k)            = aosum%Asymmetry_Factor(1:k)            + ao2%Asymmetry_Factor(1:k)
+    aosum%Backscat_Coefficient(1:k)        = aosum%Backscat_Coefficient(1:k)        + ao2%Backscat_Coefficient(1:k)
     aosum%Delta_Truncation(1:k)            = aosum%Delta_Truncation(1:k)            + ao2%Delta_Truncation(1:k)
     aosum%Phase_Coefficient(0:ic,1:ip,1:k) = aosum%Phase_Coefficient(0:ic,1:ip,1:k) + ao2%Phase_Coefficient(0:ic,1:ip,1:k)
 
@@ -1536,6 +1553,7 @@ ELEMENTAL FUNCTION CRTM_AtmOptics_Subtract( ao1, ao2 ) RESULT( aodiff )
     aodiff%Optical_Depth(1:k)               = aodiff%Optical_Depth(1:k)               - ao2%Optical_Depth(1:k)
     aodiff%Single_Scatter_Albedo(1:k)       = aodiff%Single_Scatter_Albedo(1:k)       - ao2%Single_Scatter_Albedo(1:k)
     aodiff%Asymmetry_Factor(1:k)            = aodiff%Asymmetry_Factor(1:k)            - ao2%Asymmetry_Factor(1:k)
+    aodiff%Backscat_Coefficient(1:k)        = aodiff%Backscat_Coefficient(1:k)        - ao2%Backscat_Coefficient(1:k)
     aodiff%Delta_Truncation(1:k)            = aodiff%Delta_Truncation(1:k)            - ao2%Delta_Truncation(1:k)
     aodiff%Phase_Coefficient(0:ic,1:ip,1:k) = aodiff%Phase_Coefficient(0:ic,1:ip,1:k) - ao2%Phase_Coefficient(0:ic,1:ip,1:k)
 

diff --git a/src/AtmScatter/CRTM_AOD_Module.f90 b/src/AtmScatter/CRTM_AOD_Module.f90
@@ -225,7 +225,6 @@ FUNCTION CRTM_AOD( &
     ! ------------
     Profile_Loop: DO m = 1, n_Profiles
 
-
       ! Check the aerosol coeff. data for cases with aerosols
       IF( Atmosphere(m)%n_Aerosols > 0 .AND. .NOT. CRTM_AerosolCoeff_IsLoaded() )THEN
          Error_Status = FAILURE