diff --git a/bin.rams/RAMSIN.arctictest b/bin.rams/RAMSIN.arctictest
index a8720fc..b7c6eb5 100644
--- a/bin.rams/RAMSIN.arctictest
+++ b/bin.rams/RAMSIN.arctictest
@@ -214,10 +214,26 @@
 !-----------------------------------------------------------------------
    FRQLITE  = 0.,            ! Frequency for "lite" files
                              !   = 0 : no lite files
-   NLITE_VARS=0,
+NLITE_VARS=21,
 
    LITE_VARS= 'GLAT','GLON','TOPT','PATCH_AREA', !Keep these always
               'UP','VP','WP', !This would be useful for hi-frequency trajectories
+              'THETA','PI','RV','RTP','RCP','RRP','RDP',
+               'CCP','CRP','CDP',
+               'SWUP','SWDN','LWUP','LWDN',
+
+   ! Accuracy of the lite vars if using ZFP/lossy compression
+   ! Use 0 for lossless compression
+   ! Using nonzero values requires compiling with -DENABLE_ZFP_COMPRESSION
+   ! on the C compiler.
+
+   ACC_LT_VAR= 
+               0.0, 0.0, 0.1,0.01,              ! 'GLAT','GLON','TOPT','PATCH_AREA',
+               0.001, 0.001, 0.001,             !'UP','VP','WP',
+               0.001,0.001,0.00001,0.00001,0.000001,0.000001,0.000001,   ! 'THETA','PI','RV','RTP','RCP','RRP','RDP',
+               0.00001,0.00001,0.000001,0.000001,0.000001,                                 ! 'CCP','CRP','CDP',
+               0.001, 0.001, 0.001,0.001,      ! 'SWUP','SWDN','LWUP','LWDN',
+
 
 !-----------------------------------------------------------------------
    AVGTIM   = 0.,            ! Averaging time for analysis variables
diff --git a/include.mk b/include.mk
index c97acfe..b126d3b 100644
--- a/include.mk
+++ b/include.mk
@@ -29,8 +29,9 @@ RAMS_VERSION=6.3.04
 # Typically can use "parallel" for either, but some supercomputers require
 # use of the serial executable.
 #############################################################################
-HDF5_ROOT=
-#/share/apps/hdf5-1.10.1/intel
+HDF5_ROOT=#/share/apps/hdf5-1.10.1/intel
+#ZFP Compression requites ZFP
+HDZ_ZFP_ROOT=
 
 #############################################################################
 # Set root locations for parallel processing MPI software.
@@ -65,8 +66,8 @@ UTILS_INCS=-I$(MODEL)/include
 #HDF5_LIBS=-L$(HDF5_ROOT)/lib -lhdf5_hl -lhdf5
 #HDF5_INCS=-I$(HDF5_ROOT)/include
 #HDF5_DEFS=
-HDF5_LIBS= -lhdf5_hl -lhdf5
-HDF5_INCS=
+HDF5_LIBS= -lhdf5_hl -lhdf5 #
+HDF5_INCS=-I$(HDF5_ROOT)/include #-I$(H5Z_ZFP_ROOT)/include 
 HDF5_DEFS=
 
 #############################################################################
@@ -175,6 +176,7 @@ LIBS=-L/usr/lib/x86_64-linux-gnu -lrt -lpthread -lz -lsz
 #C_COMP=gcc
 C_COMP=/home/smsaleeb/software/mpich-3.3.2/bin/mpicc
 C_OPTS=-O3 -DUNDERSCORE -DLITTLE -std=gnu99 -DENABLE_PARALLEL_COMPRESSION -w
+#-DENABLE_ZFP_COMPRESSION
 #C_OPTS=-O3 -DUNDERSCORE -DLITTLE -std=gnu99 -DRAMS_DOUBLE_PREC \
 #  -DENABLE_PARALLEL_COMPRESSION -w
 
diff --git a/src/bc/rbnd_nonscalar.f90 b/src/bc/rbnd_nonscalar.f90
index 20f9afe..6f5b017 100644
--- a/src/bc/rbnd_nonscalar.f90
+++ b/src/bc/rbnd_nonscalar.f90
@@ -892,6 +892,9 @@ Subroutine rad_bcond (m1,m2,m3,radiate,ibcon)
  do j = 1,m3
   do k = 1,m1
    radiate%fthrd(k,1,j)  = radiate%fthrd(k,2,j)
+   ! GRL 2024-03-22 added lw and sw heating rates
+   radiate%fthrdlw(k,1,j)  = radiate%fthrdlw(k,2,j)
+   radiate%fthrdsw(k,1,j)  = radiate%fthrdsw(k,2,j)
    if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
      radiate%bext(k,1,j)   = radiate%bext(k,2,j)
    endif
@@ -912,6 +915,9 @@ Subroutine rad_bcond (m1,m2,m3,radiate,ibcon)
  do j = 1,m3
   do k = 1,m1
    radiate%fthrd(k,m2,j) = radiate%fthrd(k,m2-1,j)
+   ! GRL 2024-03-22 added lw and sw heating rates
+   radiate%fthrdlw(k,m2,j) = radiate%fthrdlw(k,m2-1,j)
+   radiate%fthrdsw(k,m2,j) = radiate%fthrdsw(k,m2-1,j)
    if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
      radiate%bext(k,m2,j)  = radiate%bext(k,m2-1,j)
    endif
@@ -932,6 +938,9 @@ Subroutine rad_bcond (m1,m2,m3,radiate,ibcon)
   do i = 1,m2
    do k = 1,m1
      radiate%fthrd(k,i,1)   = radiate%fthrd(k,i,2)
+     ! GRL 2024-03-22 added lw and sw heating rates
+     radiate%fthrdlw(k,i,1)   = radiate%fthrdlw(k,i,2)
+     radiate%fthrdsw(k,i,1)   = radiate%fthrdsw(k,i,2)
      if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
        radiate%bext(k,i,1)    = radiate%bext(k,i,2)
      endif
@@ -952,6 +961,9 @@ Subroutine rad_bcond (m1,m2,m3,radiate,ibcon)
   do i = 1,m2
    do k = 1,m1
      radiate%fthrd(k,i,m3)  = radiate%fthrd(k,i,m3-1)
+     ! GRL 2024-03-22 added lw and sw heating rates
+     radiate%fthrdlw(k,i,m3)  = radiate%fthrdlw(k,i,m3-1)
+     radiate%fthrdsw(k,i,m3)  = radiate%fthrdsw(k,i,m3-1)
      if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
        radiate%bext(k,i,m3)   = radiate%bext(k,i,m3-1)
      endif
@@ -1047,6 +1059,7 @@ Subroutine sfcrad_bcond (m2,m3,radiate,ibcon)
   radiate%rshort(1,j)  = radiate%rshort(2,j)
   radiate%rlong(1,j)   = radiate%rlong(2,j)
   radiate%rlongup(1,j) = radiate%rlongup(2,j)
+  radiate%rlontop(1,j) = radiate%rlontop(2,j)
   radiate%albedt(1,j)  = radiate%albedt(2,j)
   radiate%cosz(1,j)    = radiate%cosz(2,j)
   radiate%aodt(1,j)    = radiate%aodt(2,j)
@@ -1059,6 +1072,7 @@ Subroutine sfcrad_bcond (m2,m3,radiate,ibcon)
   radiate%rshort(m2,j)  = radiate%rshort(m2-1,j)
   radiate%rlong(m2,j)   = radiate%rlong(m2-1,j)
   radiate%rlongup(m2,j) = radiate%rlongup(m2-1,j)
+  radiate%rlontop(m2,j) = radiate%rlontop(m2-1,j)
   radiate%albedt(m2,j)  = radiate%albedt(m2-1,j)
   radiate%cosz(m2,j)    = radiate%cosz(m2-1,j)
   radiate%aodt(m2,j)    = radiate%aodt(m2-1,j)
@@ -1071,6 +1085,7 @@ Subroutine sfcrad_bcond (m2,m3,radiate,ibcon)
      radiate%rshort(i,1)  = radiate%rshort(i,2)
      radiate%rlong(i,1)   = radiate%rlong(i,2)
      radiate%rlongup(i,1) = radiate%rlongup(i,2)
+     radiate%rlontop(i,1) = radiate%rlontop(i,2)
      radiate%albedt(i,1)  = radiate%albedt(i,2)
      radiate%cosz(i,1)    = radiate%cosz(i,2)
      radiate%aodt(i,1)    = radiate%aodt(i,2)
@@ -1083,6 +1098,7 @@ Subroutine sfcrad_bcond (m2,m3,radiate,ibcon)
      radiate%rshort(i,m3)  = radiate%rshort(i,m3-1)
      radiate%rlong(i,m3)   = radiate%rlong(i,m3-1)
      radiate%rlongup(i,m3) = radiate%rlongup(i,m3-1)
+     radiate%rlontop(i,m3) = radiate%rlontop(i,m3-1)
      radiate%albedt(i,m3)  = radiate%albedt(i,m3-1)
      radiate%cosz(i,m3)    = radiate%cosz(i,m3-1)
      radiate%aodt(i,m3)    = radiate%aodt(i,m3-1)
diff --git a/src/bc/rbnd_trsetsns.f90 b/src/bc/rbnd_trsetsns.f90
index d619582..8e52b46 100644
--- a/src/bc/rbnd_trsetsns.f90
+++ b/src/bc/rbnd_trsetsns.f90
@@ -15,6 +15,9 @@ Subroutine trsets_ns ()
 
 if(ilwrtyp+iswrtyp > 0) then
  CALL trsetns (mzp,mxp,myp,radiate_g(ngrid)%fthrd,1)
+ !GRL 2024-03-22 added lw and sw heating rates
+ CALL trsetns (mzp,mxp,myp,radiate_g(ngrid)%fthrdlw,1)
+ CALL trsetns (mzp,mxp,myp,radiate_g(ngrid)%fthrdsw,1)
 endif
 
 if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
diff --git a/src/io/anal_write.f90 b/src/io/anal_write.f90
index 11d038f..08d0482 100644
--- a/src/io/anal_write.f90
+++ b/src/io/anal_write.f90
@@ -37,6 +37,9 @@ Subroutine anal_write (vtype)
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
 real, dimension(:,:,:,:), allocatable :: temp_var1, temp_var2
+! Holding variable for ZFP accuracy. This will be 0 for all but lite files, where 
+! it is user set. If it is 0, that means to not run ZFP/lossy compression. 
+real :: zfp_accuracy
 
 type (head_table), allocatable,save :: aw_table(:)
 
@@ -139,17 +142,21 @@ Subroutine anal_write (vtype)
       iwrite=0
       if(vtype == 'INST' .and. vtab_r(nv,ngr)%ianal == 1) then
          iwrite=1
+         zfp_accuracy = 0
          v_pointer => vtab_r(nv,ngr)%var_p
       elseif(vtype == 'LITE' .and. vtab_r(nv,ngr)%ilite == 1) then
          iwrite=1
+         zfp_accuracy = vtab_r(nv,ngr)%var_acc
          v_pointer => vtab_r(nv,ngr)%var_p
       elseif(vtype == 'MEAN' .and. vtab_r(nv,ngr)%imean == 1 .and. &
                                    vtab_r(nv,ngr)%ianal == 1) then
          iwrite=1
+          zfp_accuracy = 0
          v_pointer => vtab_r(nv,ngr)%var_m
       elseif(vtype == 'BOTH' .and. vtab_r(nv,ngr)%ilite == 1 .and. &
                                    vtab_r(nv,ngr)%imean == 1) then
          iwrite=1
+          zfp_accuracy = 0
          v_pointer => vtab_r(nv,ngr)%var_m
       endif
 
@@ -240,9 +247,9 @@ Subroutine anal_write (vtype)
 
          CALL shdf5_set_hs_select (vtab_r(nv,ngr)%idim_type,'W',ngr &
                 ,mem_select,file_select,file_chunks)
-
+         
          CALL shdf5_orec (h5_fid,iphdf5,varn &
-                ,mem_select,file_select,file_chunks,rvara=temp_var1)
+                ,mem_select,file_select,file_chunks,zfp_accuracy,rvara=temp_var1)
          !print*,'done all ',vtype,' ',trim(varn)
 
          deallocate(temp_var1)
@@ -368,7 +375,7 @@ Subroutine anal_write (vtype)
          CALL shdf5_set_hs_select (idtype,'W',ngr &
                  ,mem_select,file_select,file_chunks)
          CALL shdf5_orec (h5_fid,iphdf5,varn &
-                 ,mem_select,file_select,file_chunks,rvara=temp_var2)
+                 ,mem_select,file_select,file_chunks, zfp_accuracy,rvara=temp_var2)
          !print*,'done xtra ',vtype,' ',trim(varn)
 
       endif
diff --git a/src/io/io_params.f90 b/src/io/io_params.f90
index 46ff9c1..8c3a5e3 100644
--- a/src/io/io_params.f90
+++ b/src/io/io_params.f90
@@ -6,6 +6,8 @@ Module io_params
 implicit none
 
 character(len=32) :: lite_vars(maxlite)
+! this must be fixed at a 32 bit float because that is what the C library expects
+real(kind=4) :: lite_var_acc(maxlite)
 character(len=strl1) :: hfilin,afilepref
 
 integer :: ipast_sfc
diff --git a/src/io/rname.f90 b/src/io/rname.f90
index a081a55..84cefda 100644
--- a/src/io/rname.f90
+++ b/src/io/rname.f90
@@ -25,7 +25,7 @@ Subroutine nvfillm (group,vr,ii,ff,cc,nv)
 real :: ff
 integer :: ii,nv
 integer :: inrflg
-integer, parameter ::nvgrid=38,nvstrt=77,nvindat=152,nvsound=13
+integer, parameter ::nvgrid=38,nvstrt=78,nvindat=152,nvsound=13
 integer ::  igrids(nvgrid),istart(nvstrt),iindat(nvindat),isound(nvsound)
 character(len=16) :: grids(nvgrid),start(nvstrt),indat(nvindat),sound(nvsound)
 data igrids/nvgrid*0/,istart/nvstrt*0/,iindat/nvindat*0/,isound/nvsound*0/
@@ -48,7 +48,7 @@ Subroutine nvfillm (group,vr,ii,ff,cc,nv)
      ,'RODA_UPA0','RODA_HGT','RODA_ZFAC','ODA_SFC_TIL','ODA_SFC_TEL'     &
      ,'ODA_UPA_TIL','ODA_UPA_TEL','HFILIN','IPAST_SFC','ICLOBBER'        &
      ,'IOUTPUT','AFILEPREF','FRQSTATE','FRQST_KEEP','FRQLITE'            &
-     ,'NLITE_VARS','LITE_VARS','AVGTIM','FRQMEAN','FRQBOTH','TOPFILES'   &
+     ,'NLITE_VARS','LITE_VARS','ACC_LT_VAR','AVGTIM','FRQMEAN','FRQBOTH','TOPFILES'   &
      ,'SFCFILES','SSTFPFX','NDVIFPFX','ITOPTFLG','ISSTFLG','IVEGTFLG'    &
      ,'ISOILFLG','NDVIFLG','IUPDNDVI','IUPDSST','ITOPTFN'                &
      ,'ISSTFN','IVEGTFN','ISOILFN','NDVIFN','ITOPSFLG','TOPTENH'         &
@@ -208,6 +208,7 @@ END SUBROUTINE varseti
  IF(VR.EQ.'FRQLITE')     CALL varsetf (VR,FRQLITE,NV,1,FF,0.,1.E20)
  IF(VR.EQ.'NLITE_VARS')  CALL varseti (VR,NLITE_VARS,NV,1,II,0,MAXLITE)
  IF(VR.EQ.'LITE_VARS')   CALL varsetc (VR,LITE_VARS(NV),NV,MAXLITE,CC,1,32)
+ IF(VR.EQ.'ACC_LT_VAR')  CALL varsetf (VR,LITE_VAR_ACC(NV),NV,MAXLITE,FF,0.,1.e20)
  IF(VR.EQ.'AVGTIM')      CALL varsetf (VR,AVGTIM,NV,1,FF,-1.E20,1.E20)
  IF(VR.EQ.'FRQMEAN')     CALL varsetf (VR,FRQMEAN,NV,1,FF,0.,1.E20)
  IF(VR.EQ.'FRQBOTH')     CALL varsetf (VR,FRQBOTH,NV,1,FF,0.,1.E20)
diff --git a/src/isan/isan_io.f90 b/src/isan/isan_io.f90
index eaa0508..e4c322a 100644
--- a/src/isan/isan_io.f90
+++ b/src/isan/isan_io.f90
@@ -12,6 +12,8 @@ Subroutine isenio (h5_fid,iphdf5,inout,n1,n2,igrid)
 character(len=*) :: inout
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
+! this will always be 0 for ISAN
+real(kind=4):: zfp_accuracy
 
 if(inout == 'IN') THEN
 
@@ -110,19 +112,19 @@ Subroutine isenio (h5_fid,iphdf5,inout,n1,n2,igrid)
    ! scalar values
    CALL shdf5_set_hs_select (1,'W',igrid,mem_select,file_select,file_chunks)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_year',mem_select,file_select  &
-                   ,file_chunks,ivars=iyear)
+                   ,file_chunks, zfp_accuracy,ivars=iyear)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_month',mem_select,file_select &
-                   ,file_chunks,ivars=imonth)
+                   ,file_chunks, zfp_accuracy,ivars=imonth)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_date',mem_select,file_select  &
-                   ,file_chunks,ivars=idate)
+                   ,file_chunks, zfp_accuracy,ivars=idate)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_hour',mem_select,file_select  &
-                   ,file_chunks,ivars=ihour)
+                   ,file_chunks, zfp_accuracy,ivars=ihour)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_nx',mem_select,file_select    &
-                   ,file_chunks,ivars=n1)
+                   ,file_chunks, zfp_accuracy,ivars=n1)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_ny',mem_select,file_select    &
-                   ,file_chunks,ivars=n2)
+                   ,file_chunks, zfp_accuracy,ivars=n2)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_nisn',mem_select,file_select  &
-                   ,file_chunks,ivars=nisn)
+                   ,file_chunks, zfp_accuracy,ivars=nisn)
 
    ! Vector, length is nisn 
    ! To indicate to shdf5_set_h5_select that the variable is a vector of
@@ -130,73 +132,73 @@ Subroutine isenio (h5_fid,iphdf5,inout,n1,n2,igrid)
    CALL shdf5_set_hs_select (-nisn,'W',igrid,mem_select,file_select &
                             ,file_chunks)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_levth',mem_select,file_select &
-                   ,file_chunks,ivara=levth)
+                   ,file_chunks, zfp_accuracy,ivara=levth)
 
    npts=n1*n2*nisn
    ! ISAN 3D isentropic variables
    CALL shdf5_set_hs_select (8,'W',igrid,mem_select,file_select,file_chunks)
    CALL vmissw (pi_u,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_u',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (pi_v,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_v',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (pi_p,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_p',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (pi_s,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_s',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (pi_r,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'isen_r',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
 
    npts=n1*n2
    ! ISAN 2D variables
    CALL shdf5_set_hs_select (2,'W',igrid,mem_select,file_select,file_chunks)
    CALL vmissw (rs_u,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_u',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_v,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_v',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_p,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_p',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_t,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_t',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_r,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_r',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_s,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_s',mem_select,file_select    &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_top,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_topo',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_qual,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_qual',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    
    CALL vmissw (rs_soilmoist1,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_soilmoist1',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_soilmoist2,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_soilmoist2',mem_select,file_select &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_soiltemp1,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_soiltemp1',mem_select,file_select  &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_soiltemp2,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_soiltemp2',mem_select,file_select  &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_snowmass,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_snowmass',mem_select,file_select   &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
    CALL vmissw (rs_snowdepth,npts,pi_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sfc_snowdepth',mem_select,file_select  &
-                   ,file_chunks,rvara=pi_scra)
+                   ,file_chunks, zfp_accuracy,rvara=pi_scra)
 
    print 201,' *****  Isentropic file written *************'  &
         ,iyear,imonth,idate,ihour,n1,n2,nisn  &
@@ -227,6 +229,7 @@ Subroutine sigzio (h5_fid,iphdf5,inout,n1,n2,igrid)
 character(len=*) :: inout
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
+real(kind=4):: zfp_accuracy = 0
 
 if(inout == 'IN') then
 
@@ -286,7 +289,7 @@ Subroutine sigzio (h5_fid,iphdf5,inout,n1,n2,igrid)
    CALL shdf5_set_hs_select (1,'W',igrid,mem_select,file_select &
                             ,file_chunks)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_nsigz',mem_select,file_select &
-                   ,file_chunks,ivars=nsigz)
+                   ,file_chunks, zfp_accuracy,ivars=nsigz)
 
    ! Vector, length is nsigz.
    ! To indicate to shdf5_set_h5_select that the variable is a vector of
@@ -294,7 +297,7 @@ Subroutine sigzio (h5_fid,iphdf5,inout,n1,n2,igrid)
    CALL shdf5_set_hs_select (-nsigz,'R',igrid,mem_select,file_select &
                             ,file_chunks)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_sigz',mem_select,file_select &
-                   ,file_chunks,rvara=sigz)
+                   ,file_chunks, zfp_accuracy,rvara=sigz)
 
    npts=n1*n2*nsigz
    ! ISAN 3D sigma-z variables
@@ -302,19 +305,19 @@ Subroutine sigzio (h5_fid,iphdf5,inout,n1,n2,igrid)
                             ,file_chunks)
    CALL vmissw (ps_u,npts,ps_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_u',mem_select,file_select &
-                   ,file_chunks,rvara=ps_scra)
+                   ,file_chunks, zfp_accuracy,rvara=ps_scra)
    CALL vmissw (ps_v,npts,ps_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_v',mem_select,file_select &
-                   ,file_chunks,rvara=ps_scra)
+                   ,file_chunks, zfp_accuracy,rvara=ps_scra)
    CALL vmissw (ps_p,npts,ps_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_p',mem_select,file_select &
-                   ,file_chunks,rvara=ps_scra)
+                   ,file_chunks, zfp_accuracy,rvara=ps_scra)
    CALL vmissw (ps_t,npts,ps_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_t',mem_select,file_select &
-                   ,file_chunks,rvara=ps_scra)
+                   ,file_chunks, zfp_accuracy,rvara=ps_scra)
    CALL vmissw (ps_r,npts,ps_scra,1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'sigz_r',mem_select,file_select &
-                   ,file_chunks,rvara=ps_scra)
+                   ,file_chunks, zfp_accuracy,rvara=ps_scra)
 
    print 201,' *****  Sigma-z file written *************'  &
         ,iyear,imonth,idate,ihour,n1,n2,nsigz   &
diff --git a/src/isan/write_varf.f90 b/src/isan/write_varf.f90
index b72d8b3..14a4fd2 100644
--- a/src/isan/write_varf.f90
+++ b/src/isan/write_varf.f90
@@ -16,6 +16,7 @@ Subroutine write_varf ()
 integer :: iphdf5
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
+real(kind=4):: zfp_accuracy = 0
 
 if (nmachs .gt. 1) then
   iphdf5 = 1
@@ -36,74 +37,74 @@ Subroutine write_varf ()
    ! scalar values
    CALL shdf5_set_hs_select (1,'W',ng,mem_select,file_select,file_chunks)
    CALL shdf5_orec (h5_fid,iphdf5,'year'   ,mem_select,file_select &
-       ,file_chunks,ivars=iyear)     
+       ,file_chunks, zfp_accuracy,ivars=iyear)     
    CALL shdf5_orec (h5_fid,iphdf5,'month'  ,mem_select,file_select &
-       ,file_chunks,ivars=imonth)    
+       ,file_chunks, zfp_accuracy,ivars=imonth)    
    CALL shdf5_orec (h5_fid,iphdf5,'day'    ,mem_select,file_select &
-       ,file_chunks,ivars=idate)     
+       ,file_chunks, zfp_accuracy,ivars=idate)     
    CALL shdf5_orec (h5_fid,iphdf5,'hour'   ,mem_select,file_select &
-       ,file_chunks,ivars=ihour)     
+       ,file_chunks, zfp_accuracy,ivars=ihour)     
    CALL shdf5_orec (h5_fid,iphdf5,'nx'     ,mem_select,file_select &
-       ,file_chunks,ivars=nnxp(ng))  
+       ,file_chunks, zfp_accuracy,ivars=nnxp(ng))  
    CALL shdf5_orec (h5_fid,iphdf5,'ny'     ,mem_select,file_select &
-       ,file_chunks,ivars=nnyp(ng))  
+       ,file_chunks, zfp_accuracy,ivars=nnyp(ng))  
    CALL shdf5_orec (h5_fid,iphdf5,'nz'     ,mem_select,file_select &
-       ,file_chunks,ivars=nnzp(ng))  
+       ,file_chunks, zfp_accuracy,ivars=nnzp(ng))  
    CALL shdf5_orec (h5_fid,iphdf5,'polelat',mem_select,file_select &
-       ,file_chunks,rvars=polelat)        
+       ,file_chunks, zfp_accuracy,rvars=polelat)        
    CALL shdf5_orec (h5_fid,iphdf5,'polelon',mem_select,file_select &
-       ,file_chunks,rvars=polelon)       
+       ,file_chunks, zfp_accuracy,rvars=polelon)       
    CALL shdf5_orec (h5_fid,iphdf5,'dx'     ,mem_select,file_select &
-       ,file_chunks,rvars=deltaxn(ng))         
+       ,file_chunks, zfp_accuracy,rvars=deltaxn(ng))         
    CALL shdf5_orec (h5_fid,iphdf5,'dz'     ,mem_select,file_select &
-       ,file_chunks,rvars=deltazn(ng))     
+       ,file_chunks, zfp_accuracy,rvars=deltazn(ng))     
    CALL shdf5_orec (h5_fid,iphdf5,'dzrat'  ,mem_select,file_select &
-       ,file_chunks,rvars=dzrat)     
+       ,file_chunks, zfp_accuracy,rvars=dzrat)     
    CALL shdf5_orec (h5_fid,iphdf5,'dzmax'  ,mem_select,file_select &
-       ,file_chunks,rvars=dzmax)
+       ,file_chunks, zfp_accuracy,rvars=dzmax)
 
    ! Atmospheric 3D vars
    CALL shdf5_set_hs_select (3,'W',ng,mem_select,file_select,file_chunks)
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_u,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'UP',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_v,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'VP',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_p,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'PI',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_t,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'THETA',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_r,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'RV',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
 
    CALL rearrange (nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_cond,rr_scr3)
    CALL shdf5_orec (h5_fid,iphdf5,'COND',mem_select,file_select &
-       ,file_chunks,rvara=rr_scr3)
+       ,file_chunks, zfp_accuracy,rvara=rr_scr3)
    
    ! Atmospheric 2D vars
    CALL shdf5_set_hs_select (2,'W',ng,mem_select,file_select,file_chunks)
    CALL vmissw (is_grids(ng)%rr_soilmoist1(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SOILMOIST1',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
    CALL vmissw (is_grids(ng)%rr_soilmoist2(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SOILMOIST2',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
    CALL vmissw (is_grids(ng)%rr_soiltemp1(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SOILTEMP1',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
    CALL vmissw (is_grids(ng)%rr_soiltemp2(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SOILTEMP2',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
    CALL vmissw (is_grids(ng)%rr_snowmass(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SNOWMASS',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
    CALL vmissw (is_grids(ng)%rr_snowdepth(1,1),nxyp,rr_vt2da(1),1E30,-9999.)
    CALL shdf5_orec (h5_fid,iphdf5,'SNOWDEPTH',mem_select,file_select &
-       ,file_chunks,rvara=rr_vt2da)
+       ,file_chunks, zfp_accuracy,rvara=rr_vt2da)
 
    CALL shdf5_close (h5_fid)
 
diff --git a/src/lib/hdf5_f2c.c b/src/lib/hdf5_f2c.c
index 005102d..34b4d2c 100644
--- a/src/lib/hdf5_f2c.c
+++ b/src/lib/hdf5_f2c.c
@@ -6,6 +6,11 @@
 
 #include "hdf5.h"
 
+// for ZFP compression
+#ifdef ENABLE_ZFP_COMPRESSION
+#include "H5Zzfp_plugin.h"
+#endif
+
 // Size of temporary hsize_t arrays
 #define HDF5_MAX_DIMS 10
 
@@ -294,8 +299,8 @@ hid_t  fh5_set_hdf5_dtype (int type_code);
 // See the comments for fh5d_read() for info about the purpose of the hsize_t arrays.
 //
  void fh5d_write (int64_t *fileid,char *dname,int *h5type,int *iphdf5,int *m_ndims
-                ,int *m_select,int *f_ndims,int *f_select,int *f_csize,void *buf
-                ,int *hdferr)
+                ,int *m_select,int *f_ndims,int *f_select,int *f_csize, int *enable_truncation,
+                float *truncation_accuracy, void *buf, int *hdferr)
   {
   hid_t dsetid, fspcid, mspcid, propid;
   int i;
@@ -318,6 +323,12 @@ hid_t  fh5_set_hdf5_dtype (int type_code);
   hsize_t f_stride[HDF5_MAX_DIMS];
 
   hsize_t chunk_size[HDF5_MAX_DIMS]; // for the file data space
+  //bool enable_truncation = true;
+  // ZFP Parameter array
+  int cd_nelmts = 10;
+  unsigned int cd_vals[cd_nelmts];
+  //printf("ZFP Options: enable? %i accuracy? %f \n", *enable_truncation, *truncation_accuracy);
+  //float accuracy = 0.1;
   
   //******************* CONVERT INPUT ARGUMENTS *************************
 
@@ -361,11 +372,37 @@ hid_t  fh5_set_hdf5_dtype (int type_code);
     *hdferr = propid; 
     printf("Error in H5Pcreate, propid: %d", propid);
     return;}
-  
+
+  // Enable compression options. 
+  if(*enable_truncation == 1){
+    #ifdef ENABLE_ZFP_COMPRESSION
+    herr = H5Pset_chunk (propid,*f_ndims,chunk_size);
+      if (herr < 0) { 
+    *hdferr = herr; 
+    printf("Error in chunking, herr: %d", herr);
+    return;}
+
+    //printf("Truncating variable %s to accuracy %f\n", dname, *truncation_accuracy);
+    H5Pset_zfp_accuracy_cdata(*truncation_accuracy, cd_nelmts, cd_vals);
+    herr = H5Pset_filter(propid, H5Z_FILTER_ZFP, H5Z_FLAG_MANDATORY, cd_nelmts, cd_vals);
+
+      if (herr < 0) { 
+    *hdferr = herr; 
+    printf("Error in ZFP, herr: %d", herr);
+    return;}
+
+    #else
+    // we want truncation on, but the ZFP library isn't included.
+    *hdferr=-9999;
+    printf("Error: ZFP compression requested, but ZFP library not included.\n");
+    return;
+    #endif
+  }
+  else{
   herr = H5Pset_chunk (propid,*f_ndims,chunk_size);
   #if H5_VERSION_GE(1,10,3)
   #ifdef ENABLE_PARALLEL_COMPRESSION 
-  
+    
     herr = H5Pset_shuffle (propid);
     herr = H5Pset_deflate (propid,6);
   
@@ -383,6 +420,7 @@ hid_t  fh5_set_hdf5_dtype (int type_code);
     herr = H5Pset_deflate (propid,6);
     }
     #endif //H5_VERSION_GE if
+  }
   if (herr < 0) { 
     *hdferr = herr; 
     printf("Error in chunking, herr: %d", herr);
diff --git a/src/lib/hdf5_utils.f90 b/src/lib/hdf5_utils.f90
index 57837ab..a7a4a07 100644
--- a/src/lib/hdf5_utils.f90
+++ b/src/lib/hdf5_utils.f90
@@ -178,7 +178,7 @@ END SUBROUTINE shdf5_info
 
 !##############################################################################
 Subroutine shdf5_orec (fileid,iphdf5,dsetname,mem_select,file_select &
-                      ,file_chunks &
+                      ,file_chunks, zfp_accuracy &
                       ,ivara,rvara,cvara,dvara,lvara  &
                       ,ivars,rvars,cvars,dvars,lvars)
 
@@ -206,6 +206,22 @@ Subroutine shdf5_orec (fileid,iphdf5,dsetname,mem_select,file_select &
 character(len=2) :: ctype    ! Variable type: int, real, char
 integer :: hdferr ! Error flag
 
+! 1 to enable truncation, 0 to disable. 
+integer :: enable_zfp_truncation
+real*4 :: zfp_accuracy
+!print*, "Variable: ", dsetname, " ZFP accuracy: ", zfp_accuracy
+if(zfp_accuracy == 0.) then
+enable_zfp_truncation = 0
+elseif(zfp_accuracy>0.) then
+enable_zfp_truncation = 1
+else 
+print*, "Incorrect ZFP Accuracy parameter", zfp_accuracy
+stop 'zfp accuracy wrong'
+endif
+
+!print*, "Enable ZFP Truncation? ", enable_zfp_truncation
+!print*, "ZFP accuracy? ", zfp_accuracy
+
 ! Find which data type is input
     if(present(ivars)) then ; ctype='is'
 elseif(present(rvars)) then ; ctype='rs'
@@ -277,34 +293,44 @@ Subroutine shdf5_orec (fileid,iphdf5,dsetname,mem_select,file_select &
 ! Write the dataset.
 if (ctype == 'is') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,ivars,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,ivars,hdferr)
 elseif (ctype == 'rs') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,rvars,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,rvars,hdferr)
 elseif (ctype == 'cs') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,cvars,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,cvars,hdferr)
 elseif (ctype == 'ds') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,dvars,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,dvars,hdferr)
 elseif (ctype == 'ls') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,lvars,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,lvars,hdferr)
 elseif (ctype == 'ia') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,ivara,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,ivara,hdferr)
 elseif (ctype == 'ra') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,rvara,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, enable_zfp_truncation, &
+                   zfp_accuracy,rvara,hdferr)
 elseif (ctype == 'ca') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,cvara,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,cvara,hdferr)
 elseif (ctype == 'da') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,dvara,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, enable_zfp_truncation, &
+                   zfp_accuracy,dvara,hdferr)
 elseif (ctype == 'la') then
    CALL fh5d_write (fileid,trim(dsetname)//char(0),h5_type,iphdf5 &
-                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks,lvara,hdferr)
+                   ,m_ndims,m_sel,f_ndims,f_sel,f_chunks, 0, &
+                   zfp_accuracy,lvara,hdferr)
 endif
 
 if (hdferr /= 0) then
diff --git a/src/lib/parlib.c b/src/lib/parlib.c
index 013a37a..894e25e 100644
--- a/src/lib/parlib.c
+++ b/src/lib/parlib.c
@@ -9,6 +9,7 @@
   #define MPI_Request int
 #endif
 
+
 int flag_msgtag=0;
 MPI_Request mpi_msgtags[MAX_MSGTAG];
 
diff --git a/src/memory/mem_radiate.f90 b/src/memory/mem_radiate.f90
index bc084f2..07fa5be 100644
--- a/src/memory/mem_radiate.f90
+++ b/src/memory/mem_radiate.f90
@@ -7,11 +7,13 @@ Module mem_radiate
    
       ! Variables to be dimensioned by (nzp,nxp,nyp)
    real, allocatable, dimension(:,:,:) :: &
-                          fthrd,fthrdp,bext,swup,swdn,lwup,lwdn
+                          fthrd,fthrdp,bext,swup,swdn,lwup,lwdn,fthrdsw,fthrdlw
+                          !GRL 2024-03-22 added variables for lw and sw heating rates
+
                           
       ! Variables to be dimensioned by (nxp,nyp)
    real, allocatable, dimension(:,:) :: &
-                          rshort,rlong,rlongup,albedt,cosz,aodt
+                          rshort,rlong,rlongup,rlontop,albedt,cosz,aodt 
 
    End Type
    
@@ -38,9 +40,13 @@ Subroutine alloc_radiate (radiate,n1,n2,n3)
                          allocate (radiate%rshort(n2,n3))
                          allocate (radiate%rlong(n2,n3))
                          allocate (radiate%rlongup(n2,n3))
+                         allocate (radiate%rlontop(n2,n3))
                          allocate (radiate%albedt(n2,n3))
                          allocate (radiate%cosz(n2,n3))
                          allocate (radiate%aodt(n2,n3))
+                         !GRL 2024-03-22 added variables for lw and sw heating rates
+                         allocate (radiate%fthrdsw(n1,n2,n3))
+                         allocate (radiate%fthrdlw(n1,n2,n3))
       endif
       if(ilwrtyp >= 3 .or. iswrtyp >= 3) then
          allocate (radiate%bext(n1,n2,n3))
@@ -69,6 +75,7 @@ Subroutine dealloc_radiate (radiate)
    if (allocated(radiate%rshort))   deallocate (radiate%rshort)
    if (allocated(radiate%rlong))    deallocate (radiate%rlong)
    if (allocated(radiate%rlongup))  deallocate (radiate%rlongup)
+   if (allocated(radiate%rlontop))  deallocate (radiate%rlontop)
    if (allocated(radiate%albedt))   deallocate (radiate%albedt)
    if (allocated(radiate%cosz))     deallocate (radiate%cosz)
    if (allocated(radiate%aodt))     deallocate (radiate%aodt)
@@ -77,6 +84,9 @@ Subroutine dealloc_radiate (radiate)
    if (allocated(radiate%swdn))     deallocate (radiate%swdn)
    if (allocated(radiate%lwup))     deallocate (radiate%lwup)
    if (allocated(radiate%lwdn))     deallocate (radiate%lwdn)
+   !GRL 2024-03-22 added variables for lw and sw heating rates
+   if (allocated(radiate%fthrdsw))  deallocate (radiate%fthrdsw)
+   if (allocated(radiate%fthrdlw))  deallocate (radiate%fthrdlw)
 
 return
 END SUBROUTINE dealloc_radiate
@@ -120,7 +130,15 @@ Subroutine filltab_radiate (radiate,radiatem,imean,n1,n2,n3,ng)
       CALL vtables2 (radiate%lwdn(1,1,1),radiatem%lwdn(1,1,1)  &
                  ,ng, npts, imean,  &
                  'LWDN :3:anal:mpti')
-
+   !GRL 2024-03-22 added variables for lw and sw heating rates
+   if (allocated(radiate%fthrdlw))  &
+      CALL vtables2 (radiate%fthrdlw(1,1,1),radiatem%fthrdlw(1,1,1)  &
+                 ,ng, npts, imean,  &
+                 'FTHRDLW :3:anal:mpti')
+   if (allocated(radiate%fthrdsw))  &
+      CALL vtables2 (radiate%fthrdsw(1,1,1),radiatem%fthrdsw(1,1,1)  &
+                 ,ng, npts, imean,  &
+                 'FTHRDSW :3:anal:mpti')
 
    npts=n2*n3
    if (allocated(radiate%rshort))  &
@@ -135,6 +153,10 @@ Subroutine filltab_radiate (radiate,radiatem,imean,n1,n2,n3,ng)
       CALL vtables2 (radiate%rlongup(1,1),radiatem%rlongup(1,1)  &
                  ,ng, npts, imean,  &
                  'RLONGUP :2:anal:mpti')
+   if (allocated(radiate%rlontop))  &
+      CALL vtables2 (radiate%rlontop(1,1),radiatem%rlontop(1,1)  &
+                 ,ng, npts, imean,  &
+                 'RLONTOP :2:anal:mpti')
    if (allocated(radiate%albedt))  &
       CALL vtables2 (radiate%albedt(1,1),radiatem%albedt(1,1)  &
                  ,ng, npts, imean,  &
diff --git a/src/memory/var_tables.f90 b/src/memory/var_tables.f90
index f34202f..f308419 100644
--- a/src/memory/var_tables.f90
+++ b/src/memory/var_tables.f90
@@ -10,6 +10,10 @@ Module var_tables
 type var_tables_r
    
    real, pointer :: var_p,var_m
+   ! How accurate to save the variable if using ZFP compression
+   ! if this is 0, use lossless compression
+   ! This parameter is used ONLY for lite file output
+   real(kind=4) :: var_acc
    integer :: npts, idim_type
    integer :: ianal,imean,ilite,impti,impt1,irecycle_sfc
    character(len=32) :: name
diff --git a/src/memory/vtab_fill.f90 b/src/memory/vtab_fill.f90
index 3d7358c..ea4432f 100644
--- a/src/memory/vtab_fill.f90
+++ b/src/memory/vtab_fill.f90
@@ -83,6 +83,8 @@ Subroutine lite_varset ()
    
          if (vtab_r(nv,ng)%name == lite_vars(nvl) ) then
             vtab_r(nv,ng)%ilite = 1
+            vtab_r(nv,ng)%var_acc = lite_var_acc(nvl)
+            !print*, "Loading in vtable. Variable: ", vtab_r(nv,ng)%name, " ZFP Accuracy: ", vtab_r(nv,ng)%var_acc
             ifound=1
          endif
          
diff --git a/src/micro/mic_driv.f90 b/src/micro/mic_driv.f90
index 61c2596..e3495e6 100644
--- a/src/micro/mic_driv.f90
+++ b/src/micro/mic_driv.f90
@@ -98,9 +98,13 @@ Subroutine micro ()
 
 !Zero out the radiative heating rate "fthrd" if this this a radiation timestep
 !and if running Harrington Radiation as called below with microphysics loop.
+! GRL 2024-08-28 Note that we are not zeroing out fthrd, fthrdlw, fthrdsw for SW/LWtype=5 (RTE), but this should be okay because these are zeroed 
+! out and calculated different in radcalc5
 if (iswrtyp .eq. 3 .or. ilwrtyp .eq. 3 .or. iswrtyp .eq. 4 .or. ilwrtyp .eq. 4) then
   if (mod(time + .001,radfrq) .lt. dtlt .or. time .lt. .001) then
     CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrd(1,1,1))
+    CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrdlw(1,1,1))
+    CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrdsw(1,1,1))
   endif
 endif
 
@@ -283,10 +287,13 @@ Subroutine mcphys (m1,k1,k2,k3,i,j,ngr,maxnzp          &
          ,radiate%albedt (i,j)    &
          ,radiate%cosz   (i,j)    &
          ,radiate%rlongup(i,j)    &
+         ,radiate%rlontop(i,j)    &
          ,radiate%rshort (i,j)    &
          ,radiate%rlong  (i,j)    &
          ,radiate%aodt   (i,j)    &
          ,radiate%fthrd(1,i,j)    &
+         ,radiate%fthrdlw(1,i,j)    &
+         ,radiate%fthrdsw(1,i,j)    & !GRL 2024-03-22 added lw and sw heating rates
          ,radiate%bext (1,i,j)    &
          ,radiate%swup (1,i,j)    &
          ,radiate%swdn (1,i,j)    &
@@ -298,8 +305,8 @@ Subroutine mcphys (m1,k1,k2,k3,i,j,ngr,maxnzp          &
       CALL radcalc4 (m1,maxnzp,7,iswrtyp,ilwrtyp  &
          ,glat,rtgt,topt  &
          ,radiate%albedt  (i,j) ,radiate%cosz  (i,j)  &
-         ,radiate%rlongup (i,j) ,radiate%rshort(i,j)  &
-         ,radiate%rlong   (i,j)  &
+         ,radiate%rlongup (i,j), radiate%rlontop(i,j) &
+         ,radiate%rshort(i,j)  ,radiate%rlong   (i,j)  &
          ,zm,zt,rv(1),dn0(1),pi0(1),pp(1),radiate%fthrd(1,i,j),i,j,ngr &
          ,radiate%bext(1,i,j),radiate%swup(1,i,j),radiate%swdn(1,i,j) &
          ,radiate%lwup(1,i,j),radiate%lwdn(1,i,j))
@@ -307,9 +314,11 @@ Subroutine mcphys (m1,k1,k2,k3,i,j,ngr,maxnzp          &
       CALL radcalc5 (m1,maxnzp,iswrtyp,ilwrtyp  &
          ,glat,rtgt,topt  &
          ,radiate%albedt  (i,j) ,radiate%cosz  (i,j)  &
-         ,radiate%rlongup (i,j) ,radiate%rshort(i,j)  &
-         ,radiate%rlong   (i,j) ,radiate%aodt  (i,j)  &
-         ,zm,zt,rv(1),dn0(1),pi0(1),pp(1),radiate%fthrd(1,i,j),i,j,ngr &
+         ,radiate%rlongup (i,j) ,radiate%rlontop(i,j)    &
+         ,radiate%rshort(i,j) ,radiate%rlong   (i,j) &
+         ,radiate%aodt  (i,j)  &
+         ,zm,zt,rv(1),dn0(1),pi0(1),pp(1),radiate%fthrd(1,i,j) &
+         ,radiate%fthrdlw(1,i,j),radiate%fthrdsw(1,i,j),i,j,ngr & !GRL 2024-03-22 added lw and sw heating rates
          ,radiate%bext(1,i,j),radiate%swup(1,i,j),radiate%swdn(1,i,j) &
          ,radiate%lwup(1,i,j),radiate%lwdn(1,i,j))
    endif
diff --git a/src/mksfc/mksfc_ndvi.f90 b/src/mksfc/mksfc_ndvi.f90
index 25a9230..c9a22f3 100644
--- a/src/mksfc/mksfc_ndvi.f90
+++ b/src/mksfc/mksfc_ndvi.f90
@@ -200,6 +200,7 @@ Subroutine ndvi_write (ifm,ivt)
 integer :: iphdf5
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
+real(kind=4):: zfp_accuracy = 0.
 
 if (nmachs .gt. 1) then
   iphdf5 = 1
@@ -220,34 +221,34 @@ Subroutine ndvi_write (ifm,ivt)
 ! Scalar vars
 CALL shdf5_set_hs_select (1,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'year',mem_select,file_select    &
-                ,file_chunks,ivars=iyearvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=iyearvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'month',mem_select,file_select   &
-                ,file_chunks,ivars=imonthvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=imonthvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'day',mem_select,file_select     &
-                ,file_chunks,ivars=idatevn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=idatevn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'hour',mem_select,file_select    &
-                ,file_chunks,ivars=ihourvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=ihourvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'nx',mem_select,file_select      &
-                ,file_chunks,ivars=nnxp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnxp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'ny',mem_select,file_select      &
-                ,file_chunks,ivars=nnyp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnyp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'npatch',mem_select,file_select  &
-                ,file_chunks,ivars=npatch)
+                ,file_chunks, zfp_accuracy,ivars=npatch)
 CALL shdf5_orec (h5_fid,iphdf5,'dx',mem_select,file_select      &
-                ,file_chunks,rvars=deltaxn(ifm))
+                ,file_chunks, zfp_accuracy,rvars=deltaxn(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'polelat',mem_select,file_select &
-                ,file_chunks,rvars=polelat)
+                ,file_chunks, zfp_accuracy,rvars=polelat)
 CALL shdf5_orec (h5_fid,iphdf5,'polelon',mem_select,file_select &
-                ,file_chunks,rvars=polelon)
+                ,file_chunks, zfp_accuracy,rvars=polelon)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lat',mem_select,file_select  &
-                ,file_chunks,rvars=glatr)
+                ,file_chunks, zfp_accuracy,rvars=glatr)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lon',mem_select,file_select  &
-                ,file_chunks,rvars=glonr)
+                ,file_chunks, zfp_accuracy,rvars=glonr)
 
 ! Leaf3 2D vars
 CALL shdf5_set_hs_select (6,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'VEG_NDVIF',mem_select,file_select &
-                ,file_chunks,rvara=sfcfile_p(ifm)%veg_ndvif)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%veg_ndvif)
 CALL shdf5_close (h5_fid)
 
 return
diff --git a/src/mksfc/mksfc_sfc.f90 b/src/mksfc/mksfc_sfc.f90
index 029efe6..ba0fbcb 100644
--- a/src/mksfc/mksfc_sfc.f90
+++ b/src/mksfc/mksfc_sfc.f90
@@ -230,6 +230,7 @@ Subroutine sfc_write (ifm)
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
 
 real, dimension(:), allocatable :: r_scratch
+real(kind=4):: zfp_accuracy = 0.
 
 if (nmachs .gt. 1) then
   iphdf5 = 1
@@ -252,34 +253,34 @@ Subroutine sfc_write (ifm)
 ! Scalar variables
 CALL shdf5_set_hs_select (1,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'nx',mem_select,file_select       &
-                ,file_chunks,ivars=nnxp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnxp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'ny',mem_select,file_select       &
-                ,file_chunks,ivars=nnyp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnyp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'nzg',mem_select,file_select      &
-                ,file_chunks,ivars=nzg)
+                ,file_chunks, zfp_accuracy,ivars=nzg)
 CALL shdf5_orec (h5_fid,iphdf5,'npatch',mem_select,file_select   &
-                ,file_chunks,ivars=npatch)
+                ,file_chunks, zfp_accuracy,ivars=npatch)
 CALL shdf5_orec (h5_fid,iphdf5,'dx',mem_select,file_select       &
-                ,file_chunks,rvars=deltaxn(ifm))
+                ,file_chunks, zfp_accuracy,rvars=deltaxn(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'polelat',mem_select,file_select  &
-                ,file_chunks,rvars=polelat)
+                ,file_chunks, zfp_accuracy,rvars=polelat)
 CALL shdf5_orec (h5_fid,iphdf5,'polelon',mem_select,file_select  &
-                ,file_chunks,rvars=polelon)
+                ,file_chunks, zfp_accuracy,rvars=polelon)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lat',mem_select,file_select   &
-                ,file_chunks,rvars=glatr)
+                ,file_chunks, zfp_accuracy,rvars=glatr)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lon',mem_select,file_select   &
-                ,file_chunks,rvars=glonr)
+                ,file_chunks, zfp_accuracy,rvars=glonr)
 CALL shdf5_orec (h5_fid,iphdf5,'ivegtflg',mem_select,file_select &
-                ,file_chunks,ivars=ivegtflg(ifm))
+                ,file_chunks, zfp_accuracy,ivars=ivegtflg(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'isoilflg',mem_select,file_select &
-                ,file_chunks,ivars=isoilflg(ifm))
+                ,file_chunks, zfp_accuracy,ivars=isoilflg(ifm))
 
 ! Leaf3 2D vars
 CALL shdf5_set_hs_select (6,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'PATCH_AREA',mem_select,file_select &
-                ,file_chunks,rvara=sfcfile_p(ifm)%patch_area)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%patch_area)
 CALL shdf5_orec (h5_fid,iphdf5,'LEAF_CLASS',mem_select,file_select &
-                ,file_chunks,rvara=sfcfile_p(ifm)%leaf_class)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%leaf_class)
 
 ! Leaf3 Soil 3D vars
 allocate(r_scratch(mmxysp(ifm)))
@@ -287,7 +288,7 @@ Subroutine sfc_write (ifm)
                  ,r_scratch)
 CALL shdf5_set_hs_select (4,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'SOIL_TEXT',mem_select,file_select &
-                ,file_chunks,rvara=r_scratch)
+                ,file_chunks, zfp_accuracy,rvara=r_scratch)
 deallocate(r_scratch)
 
 CALL shdf5_close (h5_fid)
diff --git a/src/mksfc/mksfc_sst.f90 b/src/mksfc/mksfc_sst.f90
index 7e4a035..4fd7e8d 100644
--- a/src/mksfc/mksfc_sst.f90
+++ b/src/mksfc/mksfc_sst.f90
@@ -193,7 +193,7 @@ Subroutine sst_write (ifm,ivt)
 integer :: iphdf5
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
-
+real(kind=4)::zfp_accuracy =0.
 if (nmachs .gt. 1) then
   iphdf5 = 1
 else
@@ -213,32 +213,32 @@ Subroutine sst_write (ifm,ivt)
 ! Scalar vars
 CALL shdf5_set_hs_select (1,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'year',mem_select,file_select    &
-                ,file_chunks,ivars=iyearvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=iyearvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'month',mem_select,file_select   &
-                ,file_chunks,ivars=imonthvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=imonthvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'day',mem_select,file_select     &
-                ,file_chunks,ivars=idatevn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=idatevn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'hour',mem_select,file_select    &
-                ,file_chunks,ivars=ihourvn(ivt,ifm))
+                ,file_chunks, zfp_accuracy,ivars=ihourvn(ivt,ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'nx',mem_select,file_select      &
-                ,file_chunks,ivars=nnxp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnxp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'ny',mem_select,file_select      &
-                ,file_chunks,ivars=nnyp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnyp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'dx',mem_select,file_select      &
-                ,file_chunks,rvars=deltaxn(ifm))
+                ,file_chunks, zfp_accuracy,rvars=deltaxn(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'polelat',mem_select,file_select &
-                ,file_chunks,rvars=polelat)
+                ,file_chunks, zfp_accuracy,rvars=polelat)
 CALL shdf5_orec (h5_fid,iphdf5,'polelon',mem_select,file_select &
-                ,file_chunks,rvars=polelon)
+                ,file_chunks, zfp_accuracy,rvars=polelon)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lat',mem_select,file_select  &
-                ,file_chunks,rvars=glatr)
+                ,file_chunks, zfp_accuracy,rvars=glatr)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lon',mem_select,file_select  &
-                ,file_chunks,rvars=glonr)
+                ,file_chunks, zfp_accuracy,rvars=glonr)
 
 ! Atmos 2D vars
 CALL shdf5_set_hs_select (2,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'SEATF',mem_select,file_select &
-                ,file_chunks,rvara=sfcfile_p(ifm)%seatf)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%seatf)
 
 CALL shdf5_close (h5_fid)
 
diff --git a/src/mksfc/mksfc_top.f90 b/src/mksfc/mksfc_top.f90
index b714eaa..28eb14e 100644
--- a/src/mksfc/mksfc_top.f90
+++ b/src/mksfc/mksfc_top.f90
@@ -219,6 +219,7 @@ Subroutine top_write (ifm)
 integer :: iphdf5
 type (hdf5_select_type) :: mem_select,file_select
 integer, dimension(HDF5_MAX_DIMS) :: file_chunks
+real(kind=4):: zfp_accuracy =0.
 
 if (nmachs .gt. 1) then
   iphdf5 = 1
@@ -239,40 +240,40 @@ Subroutine top_write (ifm)
 ! Scalar vars
 CALL shdf5_set_hs_select (1,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'nx',mem_select,file_select       &
-                ,file_chunks,ivars=nnxp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnxp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'ny',mem_select,file_select       &
-                ,file_chunks,ivars=nnyp(ifm))
+                ,file_chunks, zfp_accuracy,ivars=nnyp(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'dx',mem_select,file_select       &
-                ,file_chunks,rvars=deltaxn(ifm))
+                ,file_chunks, zfp_accuracy,rvars=deltaxn(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'polelat',mem_select,file_select  &
-                ,file_chunks,rvars=polelat)
+                ,file_chunks, zfp_accuracy,rvars=polelat)
 CALL shdf5_orec (h5_fid,iphdf5,'polelon',mem_select,file_select  &
-                ,file_chunks,rvars=polelon)
+                ,file_chunks, zfp_accuracy,rvars=polelon)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lat',mem_select,file_select   &
-                ,file_chunks,rvars=glatr)
+                ,file_chunks, zfp_accuracy,rvars=glatr)
 CALL shdf5_orec (h5_fid,iphdf5,'sw_lon',mem_select,file_select   &
-                ,file_chunks,rvars=glonr)
+                ,file_chunks, zfp_accuracy,rvars=glonr)
 CALL shdf5_orec (h5_fid,iphdf5,'itoptflg',mem_select,file_select &
-                ,file_chunks,ivars=itoptflg(ifm))
+                ,file_chunks, zfp_accuracy,ivars=itoptflg(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'itopsflg',mem_select,file_select &
-                ,file_chunks,ivars=itopsflg(ifm))
+                ,file_chunks, zfp_accuracy,ivars=itopsflg(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'toptenh',mem_select,file_select  &
-                ,file_chunks,rvars=toptenh(ifm))
+                ,file_chunks, zfp_accuracy,rvars=toptenh(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'toptwvl',mem_select,file_select  &
-                ,file_chunks,rvars=toptwvl(ifm))
+                ,file_chunks, zfp_accuracy,rvars=toptwvl(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'iz0flg',mem_select,file_select   &
-                ,file_chunks,ivars=iz0flg(ifm))
+                ,file_chunks, zfp_accuracy,ivars=iz0flg(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'z0max',mem_select,file_select    &
-                ,file_chunks,rvars=z0max(ifm))
+                ,file_chunks, zfp_accuracy,rvars=z0max(ifm))
 CALL shdf5_orec (h5_fid,iphdf5,'z0fact',mem_select,file_select   &
-                ,file_chunks,rvars=z0fact)
+                ,file_chunks, zfp_accuracy,rvars=z0fact)
 
 ! Atmos 2D vars
 CALL shdf5_set_hs_select (2,'W',ifm,mem_select,file_select,file_chunks)
 CALL shdf5_orec (h5_fid,iphdf5,'TOPT',mem_select,file_select  &
-                ,file_chunks,rvara=sfcfile_p(ifm)%topt)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%topt)
 CALL shdf5_orec (h5_fid,iphdf5,'TOPZO',mem_select,file_select &
-                ,file_chunks,rvara=sfcfile_p(ifm)%topzo)
+                ,file_chunks, zfp_accuracy,rvara=sfcfile_p(ifm)%topzo)
 
 CALL shdf5_close (h5_fid)
 
diff --git a/src/mpi/mpass_init.f90 b/src/mpi/mpass_init.f90
index e4cddd2..6106f2f 100644
--- a/src/mpi/mpass_init.f90
+++ b/src/mpi/mpass_init.f90
@@ -161,6 +161,7 @@ Subroutine broadcast_config ()
        print*,'lite pack:',nm,trim(LITE_VARS(nm))
        CALL par_put_char (LITE_VARS(nm),32)
     enddo
+    CALL par_put_float (LITE_VAR_ACC, MAXLITE)
     CALL par_put_float (AVGTIM,1)
     CALL par_put_float (FRQMEAN,1)
     CALL par_put_float (FRQBOTH,1)
@@ -503,6 +504,7 @@ Subroutine broadcast_config ()
     do nm = 1, nlite_vars
        CALL par_get_char (LITE_VARS(nm),32)
     enddo
+    CALL par_get_float (LITE_VAR_ACC, MAXLITE)
     CALL par_get_float (AVGTIM,1)
     CALL par_get_float (FRQMEAN,1)
     CALL par_get_float (FRQBOTH,1)
diff --git a/src/radiate/bugsrad/driver_read.f90 b/src/radiate/bugsrad/driver_read.f90
index 21c5b0b..fd336e5 100644
--- a/src/radiate/bugsrad/driver_read.f90
+++ b/src/radiate/bugsrad/driver_read.f90
@@ -144,7 +144,7 @@ subroutine bugs_driver(nlm,amu0,alvdr,pl,tl,ql,qcwl,ncwl,qcil &
       alvdf = alvdr
       alndf = alvdr
 
-      umco2=400.
+      umco2=420.
       umch4=0. !Methane
       umn2o=0. !N2O
       fdsw=0.
diff --git a/src/radiate/rad_driv.f90 b/src/radiate/rad_driv.f90
index 3111eb3..54b754d 100644
--- a/src/radiate/rad_driv.f90
+++ b/src/radiate/rad_driv.f90
@@ -159,6 +159,9 @@ Subroutine radiate (mzp,mxp,myp,ia,iz,ja,jz)
            ! Zero out the radiative heating rate "fthrd" if this this a 
            ! radiation timestep and fthrd will be updated.
            CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrd(1,1,1))
+           CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrdlw(1,1,1))
+           CALL azero (mzp*mxp*myp,radiate_g(ngrid)%fthrdsw(1,1,1))
+
            ! Run the Harrington radiation for non-LEVEL=3 micro
            CALL radcomp3 (mzp,mxp,myp,ia,iz,ja,jz  &
             ,grid_g(ngrid)%glat       (1,1)    &
@@ -167,12 +170,15 @@ Subroutine radiate (mzp,mxp,myp,ia,iz,ja,jz)
             ,radiate_g(ngrid)%albedt  (1,1)    &
             ,radiate_g(ngrid)%cosz    (1,1)    &
             ,radiate_g(ngrid)%rlongup (1,1)    &
+            ,radiate_g(ngrid)%rlontop (1,1)    &
             ,radiate_g(ngrid)%rshort  (1,1)    &
             ,radiate_g(ngrid)%rlong   (1,1)    &
             ,radiate_g(ngrid)%aodt    (1,1)    &
             ,basic_g(ngrid)%rv        (1,1,1)  &
             ,basic_g(ngrid)%dn0       (1,1,1)  &
             ,radiate_g(ngrid)%fthrd   (1,1,1)  &
+            ,radiate_g(ngrid)%fthrdlw   (1,1,1)  &
+            ,radiate_g(ngrid)%fthrdsw   (1,1,1)  &
             ,basic_g(ngrid)%pi0       (1,1,1)  &
             ,basic_g(ngrid)%pp        (1,1,1)  &
             ,basic_g(ngrid)%theta     (1,1,1)  &
@@ -479,8 +485,8 @@ END SUBROUTINE radcomp
 
 !##############################################################################
 Subroutine radcomp3 (m1,m2,m3,ia,iz,ja,jz  &
-   ,glat,rtgt,topt,albedt,cosz,rlongup,rshort,rlong,aodt  &
-   ,rv,dn0,fthrd,pi0,pp,theta,rcp &
+   ,glat,rtgt,topt,albedt,cosz,rlongup,rlontop,rshort,rlong,aodt  &
+   ,rv,dn0,fthrd,fthrdlw,fthrdsw,pi0,pp,theta,rcp &
    ,bext,swup,swdn,lwup,lwdn &
    ,cn1np,cn1mp,cn2np,cn2mp,md1np,md1mp,md2np,md2mp &
    ,salt_film_np,salt_film_mp,salt_jet_np,salt_jet_mp &
@@ -499,8 +505,8 @@ Subroutine radcomp3 (m1,m2,m3,ia,iz,ja,jz  &
 integer :: m1,m2,m3,ia,iz,ja,jz,mcat,i,j,k,kk,k0
 
 real :: cfmasi,cparmi,glg,glgm,picpi
-real, dimension(m2,m3) :: glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong,aodt
-real, dimension(m1,m2,m3) :: dn0,rv,fthrd,pi0,pp,theta,rcp
+real, dimension(m2,m3) :: glat,rtgt,topt,cosz,albedt,rlongup,rlontop,rshort,rlong,aodt
+real, dimension(m1,m2,m3) :: dn0,rv,fthrd,fthrdlw,fthrdsw,pi0,pp,theta,rcp
 real, dimension(m1,m2,m3) :: bext,swup,swdn,lwup,lwdn
 real, dimension(m1,m2,m3) :: cn1np,cn1mp,cn2np,cn2mp,md1np,md1mp,md2np,md2mp &
   ,salt_film_np,salt_film_mp,salt_jet_np,salt_jet_mp,salt_spum_np,salt_spum_mp &
@@ -573,10 +579,13 @@ Subroutine radcomp3 (m1,m2,m3,ia,iz,ja,jz  &
             ,albedt(i,j)          &
             ,cosz(i,j)            &
             ,rlongup(i,j)         &
+            ,rlontop(i,j)         &
             ,rshort(i,j)          &
             ,rlong(i,j)           &
             ,aodt(i,j)            &
             ,fthrd(1,i,j)         &
+            ,fthrdlw(1,i,j)         &
+            ,fthrdsw(1,i,j)         &
             ,bext(1,i,j)          &
             ,swup(1,i,j)          &
             ,swdn(1,i,j)          &
@@ -680,8 +689,8 @@ END SUBROUTINE zen
 
 !##############################################################################
 Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
-   ,glat,rtgt,topt,rv,albedt,cosz,rlongup,rshort,rlong,aodt &
-   ,fthrd,bext,swup,swdn,lwup,lwdn &
+   ,glat,rtgt,topt,rv,albedt,cosz,rlongup,rlontop,rshort,rlong,aodt &
+   ,fthrd,fthrdlw,fthrdsw,bext,swup,swdn,lwup,lwdn &
    ,dn0 &
    )
 
@@ -736,6 +745,7 @@ Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
 !  albedt           : surface albedo
 !  cosz             : solar zenith angle
 !  rlongup          : upward longwave radiation at surface (W/m^2)
+!  rlontop          : upward longwave at top radiation level (W/m^2)
 !  rshort           : downward shortwave radiation at surface (W/m^2)
 !  rlong            : downward longwave radiation at surface (W/m^2)
 !  aodt             : total aerosol optical depth (band=3)
@@ -763,7 +773,7 @@ Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
 !  dl (nrad)        : air density of all radiation levels (kg/m^3)
 !  rl (nrad)        : vapor density of all radiation levels (kg/m^3)
 !  vp (nrad)        : vapor pressure (Pa)
-!  o3l (nrad)       : stores the calculated ozone profile (g/m^3)
+!  o3l (nrad)       : stores the calculated ozone profile (g/m^3) !GRL pretty sure this is also kg/m3 from values
 !  flxu (nrad)      : Total upwelling flux (W/m^2)
 !  flxd (nrad)      : Total downwelling flux (W/m^2)
 !  t (nrad)         : layer transmission func
@@ -864,8 +874,8 @@ Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
 !       ngas(3) =  O3
 
 real, save :: eps=1.e-15
-real :: glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong,aodt
-real :: zm(m1),zt(m1),dn0(m1),rv(m1),fthrd(m1)
+real :: glat,rtgt,topt,cosz,albedt,rlongup,rlontop,rshort,rlong,aodt
+real :: zm(m1),zt(m1),dn0(m1),rv(m1),fthrd(m1),fthrdsw(m1),fthrdlw(m1)
 real :: bext(m1),swup(m1),swdn(m1),lwup(m1),lwdn(m1)
 
 real, allocatable, save, dimension(:)     :: zml,ztl,dzl,pl,tl,dl,rl,o3l  &
@@ -964,9 +974,9 @@ Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
    do k = 2,m1-1
       exner(k) = (press(k)*p00i)**rocp
       !divide by exner to get potential temp heating rate
-      fthrd(k) = fthrd(k)  &
-         + (flxds(k) - flxds(k-1) + flxus(k-1) - flxus(k)) &
+      fthrdsw(k) = (flxds(k) - flxds(k-1) + flxus(k-1) - flxus(k)) &
             / (dl(k) * dzl(k) * cp * exner(k))
+      fthrd(k) = fthrd(k) + fthrdsw(k)
       swup(k) = flxus(k)
       swdn(k) = flxds(k)
     enddo
@@ -998,15 +1008,19 @@ Subroutine radcalc3 (m1,i,j,ngrid,maxnzp,mcat,iswrtyp,ilwrtyp,zm,zt &
    !Set rlong to surface level downward longwave flux.
    rlong = flxdl(1)
 
+   ! Save upwelling longwave radiation (OLR) at top radiation level
+   rlontop = flxul(nrad)
+
+
    !Make lowest level upward longwave flux equal to rlongup
    !produced from land surface models (LEAF,SiB).
    flxul(1) = rlongup
 
    do k = 2,m1-1
       !divide by exner to get potential temp heating rate
-      fthrd(k) = fthrd(k)  &
-         + (flxdl(k) - flxdl(k-1) + flxul(k-1) - flxul(k)) &
+      fthrdlw(k) = (flxdl(k) - flxdl(k-1) + flxul(k-1) - flxul(k)) &
             / (dl(k) * dzl(k) * cp * exner(k))
+      fthrd(k) = fthrd(k) +fthrdlw(k)
       lwup(k) = flxul(k)
       lwdn(k) = flxdl(k)
    enddo
@@ -1246,6 +1260,7 @@ END SUBROUTINE sum_opt
 Subroutine path_lengths (nrad,u,rl,dzl,dl,o3l,vp,pl,eps)
 
 ! Get the path lengths for the various gases...
+! GRL 2024-03-14: Converting from mass path length to pressure path length 
 
 implicit none
 
@@ -1255,8 +1270,16 @@ Subroutine path_lengths (nrad,u,rl,dzl,dl,o3l,vp,pl,eps)
 real :: rvk0,rvk1,dzl9,rmix,eps
 integer :: k
 
+! water vapor
 u(1,1) = .5 * (rl(2) + rl(1)) * 9.81 * dzl(1)
-u(1,2) = .5 * (dl(2) + dl(1)) * .45575e-3 * 9.81 * dzl(1)
+
+! carbon dioxide
+! GRL edited CO2 concentration to be 420ppm 2024-08-27
+! constant value before 9.81 (g) is CO2(parts per million)/(1,000,000) * (M_co2/M_d)
+! where M_co2/M_d is ratio of molar mass CO2 to dry air = 44.01/28.964= 1.519472
+u(1,2) = .5 * (dl(2) + dl(1)) * .6382e-3 * 9.81 * dzl(1)
+
+! ozone
 u(1,3) = o3l(1) * 9.81 * dzl(1)
 
 rvk0 = rl(1)
@@ -1266,8 +1289,9 @@ Subroutine path_lengths (nrad,u,rl,dzl,dl,o3l,vp,pl,eps)
    rmix = rvk1 / dl(k)
    vp(k) = pl(k) * rmix / (.622 + rmix)
    u(k,1) = (rvk1 - rvk0) / (log(rvk1 / rvk0) + eps) * dzl9
+   !GRL edited CO2 profile 2024-08-27
    u(k,2) = (dl(k) - dl(k-1)) / (log(dl(k) / dl(k-1)) + eps)  &
-       * dzl9 * 0.45575e-3
+         * dzl9 * 0.6283e-3
    u(k,3) = 0.5 * dzl9 * (o3l(k) + o3l(k-1))
    rvk0 = rvk1
 enddo
@@ -1408,7 +1432,7 @@ Subroutine cloud_prep_lev4 (m1,i,j,ng)
 END SUBROUTINE cloud_prep_lev4
 !##############################################################################
 Subroutine radcalc4 (m1,maxnzp,mcat,iswrtyp,ilwrtyp  &
-   ,glat,rtgt,topt,albedt,cosz,rlongup,rshort,rlong  &
+   ,glat,rtgt,topt,albedt,cosz,rlongup,rlontop,rshort,rlong  &
    ,zm,zt,rv,dn0,pi0,pp,fthrd,i,j,ngrid &
    ,bext,swup,swdn,lwup,lwdn)
 
@@ -1438,7 +1462,7 @@ Subroutine radcalc4 (m1,maxnzp,mcat,iswrtyp,ilwrtyp  &
 
 real, save :: eps=1.e-15
 real :: prsnz,prsnzp
-real :: glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong
+real :: glat,rtgt,topt,cosz,albedt,rlongup,rlontop,rshort,rlong
 real :: zm(m1),zt(m1),dn0(m1),rv(m1),pi0(m1),pp(m1),fthrd(m1)
 real :: bext(m1),swup(m1),swdn(m1),lwup(m1),lwdn(m1)
 
@@ -1515,6 +1539,8 @@ Subroutine radcalc4 (m1,maxnzp,mcat,iswrtyp,ilwrtyp  &
 
 rlong = flxdl(1)
 
+rlontop = flxul(nrad)
+
 !not integrated with BUGSrad yet. just diagnostic anyway
 bext(:)=0.
 
@@ -1522,8 +1548,8 @@ Subroutine radcalc4 (m1,maxnzp,mcat,iswrtyp,ilwrtyp  &
 END SUBROUTINE radcalc4
 !##############################################################################
 Subroutine radcalc5 (m1,maxnzp,iswrtyp,ilwrtyp  &
-   ,glat,rtgt,topt,albedt,cosz,rlongup,rshort,rlong,aodt  &
-   ,zm,zt,rv,dn0,pi0,pp,fthrd,i,j,ngrid &
+   ,glat,rtgt,topt,albedt,cosz,rlongup,rlontop,rshort,rlong,aodt  &
+   ,zm,zt,rv,dn0,pi0,pp,fthrd,fthrdlw,fthrdsw,i,j,ngrid &
    ,bext,swup,swdn,lwup,lwdn)
 
 ! RTE+RRTMGP radiation
@@ -1542,8 +1568,8 @@ Subroutine radcalc5 (m1,maxnzp,iswrtyp,ilwrtyp  &
 integer i,j,k,ii,printsound
 integer, save :: ncall = 0,nradmax
 
-real :: glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong,aodt
-real :: zm(m1),zt(m1),dn0(m1),rv(m1),pi0(m1),pp(m1),fthrd(m1)
+real :: glat,rtgt,topt,cosz,albedt,rlongup,rlontop,rshort,rlong,aodt
+real :: zm(m1),zt(m1),dn0(m1),rv(m1),pi0(m1),pp(m1),fthrd(m1),fthrdlw(m1),fthrdsw(m1)
 real :: bext(m1),swup(m1),swdn(m1),lwup(m1),lwdn(m1)
 
 real, allocatable, save, dimension(:) :: zml,ztl,dzl,pl,tl,dl,rl,o3l  &
@@ -1729,12 +1755,19 @@ Subroutine radcalc5 (m1,maxnzp,iswrtyp,ilwrtyp  &
    exner(k) = (pi0(k)+pp(k))/cp
    !divide by exner to get potential temp heating rate
    fthrd(k) = (fthsw(k)+fthlw(k))/exner(k)
+
+   !GRL 2024-03-22 added lw and sw heating rates
+   fthrdlw(k) = fthlw(k)/exner(k)
+   fthrdsw(k) = fthsw(k)/exner(k)
 enddo
 
 rshort = flxds(1)
 
 rlong = flxdl(1)
 
+rlontop = flxul(nrad)
+
+
 !not integrated with RTE+RRTMGP yet. just diagnostic anyway
 bext(:)=0.
 
diff --git a/src/radiate/rte-rrtmgp/interface/rams_rrtmgp_driver.F90 b/src/radiate/rte-rrtmgp/interface/rams_rrtmgp_driver.F90
index c1b3621..749e3ff 100644
--- a/src/radiate/rte-rrtmgp/interface/rams_rrtmgp_driver.F90
+++ b/src/radiate/rte-rrtmgp/interface/rams_rrtmgp_driver.F90
@@ -35,7 +35,7 @@ subroutine rte_rrtmgp_init()
   call stop_on_err(gas_concs%set_vmr('ch4', 0.0))
   call stop_on_err(gas_concs%set_vmr('n2o', 0.0))
   call stop_on_err(gas_concs%set_vmr('o2 ', 0.209))
-  call stop_on_err(gas_concs%set_vmr('co2', 400.e-6))
+  call stop_on_err(gas_concs%set_vmr('co2', 420.e-6))
 
   filename = trim(hucmfile)//'/../RTE-RRTMGP/rrtmgp-data-lw-g128-210809.nc'
   call load_and_init(k_dist_lw, filename, gas_concs)
diff --git a/src/turb/turb_diff.f90 b/src/turb/turb_diff.f90
index 13ae4dd..3cf04e8 100644
--- a/src/turb/turb_diff.f90
+++ b/src/turb/turb_diff.f90
@@ -525,8 +525,11 @@ Subroutine diffsclr (m1,m2,m3,ia,iz,ja,jz,jd  &
 
 !     compute 2 horizontal scalar gradients needed for dscp/dt
 
-if (ihorgrad .eq. 1 .and. scalar_tab(n,ngrid)%name /= 'THP') then
-   ! only for ihorgrad = 1 for all scalars except thp
+!GRL 2024-06-04 Made change for IHORGRAD=1 
+! Originally when IHORGRAD=1 different equations were called for THP and other scalars.
+! This seems inconsistnet and couldn't think of any reason this could be the case. 
+! In tests for RCEMIP, this helped to fix energy balance.
+if (ihorgrad .eq. 1) then
    ! vt3df=d(scp)/dx(?)
    CALL grad (m1,m2,m3,1,iz,ja,jz,scp_diffuse,vt3df,'XDIR','TPNT')
    ! vt3dg=d(scp)/dy(?)
@@ -555,17 +558,16 @@ Subroutine diffsclr (m1,m2,m3,ia,iz,ja,jz,jd  &
 
 !         horizontal flux divergence for scalars
 
-if (ihorgrad .eq. 1 .and. scalar_tab(n,ngrid)%name /= 'THP') then
-   ! This is called for ihorgrad=1 for all scalars except thp
+if (ihorgrad .eq. 1) then
+   ! This is called for ihorgrad=1 for all scalars 
    ! fluxdivx = scalar flux divergence in x direction [scp*kg/m3/s]
    ! fluxdivy = scalar flux divergence in y direction
    CALL divcart (m1,m2,m3,ia,iz,ja,jz,vt3df,fluxdivx,'XDIR','UPNT')
 
    CALL divcart (m1,m2,m3,ia,iz,ja,jz,vt3dg,fluxdivy,'YDIR','VPNT')
 
-elseif (ihorgrad .eq. 2 .or. scalar_tab(n,ngrid)%name == 'THP') then
-   ! This is called for ihorgrad=2 for all scalars, or just for 
-   !  thp if ihorgrad=1 - haven't gone through this subroutine
+elseif (ihorgrad .eq. 2) then
+   ! This is called for ihorgrad=2 for all scalars
    CALL truhor (m1,m2,m3,ia,iz,ja,jz  &
               ,scp_diffuse,fluxdivx,'xdir','dxu',grid_g(ngrid)%dxu(1,1)  &
               ,grid_g(ngrid)%topt(1,1),grid_g(ngrid)%rtgt(1,1)  &