Adding functionality to use ECMWF forecast or era5 grid instead of pre-processing grid + additional dustflag from Gareth

common/Makefile

@@ -18,6 +18,7 @@ COMMON_OBJ_F90 = $(OBJS)/chunk_utils.o \
                  $(OBJS)/common_constants.o \
                  $(OBJS)/constants_cloud_typing_pavolonis.o \
                  $(OBJS)/global_attributes.o \
+                 $(OBJS)/morphology.o \
                  $(OBJS)/interpol.o \
                  $(OBJS)/hdf.o \
                  $(OBJS)/orac_indexing.o \

common/constants_cloud_typing_pavolonis.F90

@@ -116,6 +116,8 @@ module constants_cloud_typing_pavolonis_m
    integer(sint), parameter :: OVERLAP_TYPE = 8
    integer(sint), parameter :: PROB_OPAQUE_ICE_TYPE = 9 ! missing ch3.7 due to low S/N
    integer(sint), parameter :: PROB_CLEAR_TYPE = 10 ! cold antarctic
+   integer(sint), parameter :: DUST_CLEAR_TYPE = 11
+   integer(sint), parameter :: DUST_SWITCHED_FROM_CLOUD_TYPE = 12
 
    !--- used for sunglint_mask, nise_mask
    integer(sint), parameter :: NO = 0

common/def_output_primary.F90

@@ -83,6 +83,8 @@
 ! 2018/01/19, GT: Removed QCFlag scale_factor and add_offset values, as these
 !    should only be used for packed floating point data.
 ! 2018/06/08, SP: Add satellite azimuth angle to output.
+! 2024/03/13, GT: Updated cldtype to include dust values and removed
+!    scale_factor and add_offset from the various cloud mask/type variables
 !
 ! Bugs:
 ! None known.
@@ -1823,7 +1825,10 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, &
                'opaque_ice ' // &
                'cirrus ' // &
                'overlap ' // &
-               'prob_opaque_ice'
+               'prob_opaque_ice ' // &
+               'N/A ' // &
+               'dust_clear ' // &
+               'dust_switched_from_cloud'
 
    call ncdf_def_var_byte_packed_byte( &
            ncid, &
@@ -1839,7 +1844,7 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, &
            valid_min     = output_data%cldtype_vmin, &
            valid_max     = output_data%cldtype_vmax, &
            units         = '1', &
-           flag_values   = '0b 1b 2b 3b 4b 5b 6b 7b 8b 9b', &
+           flag_values   = '0b 1b 2b 3b 4b 5b 6b 7b 8b 9b 10b 11b 12b', &
            flag_meanings = trim(adjustl(input_dummy)), &
            deflate_level = deflate_level, &
            shuffle       = shuffle_flag)

common/morphology.F90

@@ -0,0 +1,217 @@
+!-------------------------------------------------------------------------------
+! Name: morphology.F90
+!
+! Purpose:
+! Module which provides simple morphological image transformation functions.
+!
+! Description and Algorithm details:
+! Four standard mathematical morphology functions are provided:
+! erode:  "Shrinks" brighter (more positive) objects in an image, by an amount
+!         defined by some kernel matrix (completely removing objects smaller
+!         than the kernel)
+! dilate: "Grows" brighter objects by an amount defined by a kernel matrix.
+! open:   Apply erode and dilate sequentially to the same image, using the same
+!         kernel. This will have the effect of removing any bright features
+!         smaller than the kernel - i.e. can be used to remove "salt-noise" or
+!         spikes in a image.
+! close:  The opposite of open. This has the effect of removing any dark
+!         features smaller than the kernel - i.e. can be used to remove
+!         "pepper-noise" or holes in an image.
+! All four functions accept two input arguments and return an integer matrix
+! with the same dimension as the first input array:
+! A(:,:) : The image to be transformed. Must be a 2-dimensional array (matrix)
+!          of integers. Note that binary data will work (i.e. bolean
+!          transformations will be done) but the output will still be integer
+!          type.
+! B(:,:) : The transform kernel to apply. This should also be a 2-dimensional
+!          integer array, smaller than A, with an odd number of elements in
+!          both dimensions.
+!
+! Arguments:
+! None (see above)
+!
+! History:
+! 2024/03/08, GT: First version
+!
+! Bugs:
+! None known
+!-------------------------------------------------------------------------------
+
+module morphology_m
+
+   implicit none
+
+!!$  interface
+!!$     function morph_erode(A, B) result(Ae)
+!!$       ! Input variables
+!!$       integer, dimension(:,:), intent(in)        :: A
+!!$       integer, dimension(:,:), intent(in)        :: B
+!!$       ! Returned variable
+!!$       integer, dimension( size(A,1), size(A,2) ) :: Ae
+!!$     end function morph_erode
+!!$     function morph_dilate(A, B) result(Ad)
+!!$       ! Input variables
+!!$       integer, dimension(:,:), intent(in)        :: A
+!!$       integer, dimension(:,:), intent(in)        :: B
+!!$       ! Returned variable
+!!$       integer, dimension( size(A,1), size(A,2) ) :: Ad
+!!$     end function morph_dilate
+!!$     function morph_open(A, B) result(Ao)
+!!$       ! Input variables
+!!$       integer, dimension(:,:), intent(in)        :: A
+!!$       integer, dimension(:,:), intent(in)        :: B
+!!$       ! Returned variable
+!!$       integer, dimension( size(A,1), size(A,2) ) :: Ao
+!!$     end function morph_open
+!!$     function morph_close(A, B) result(Ac)
+!!$       ! Input variables
+!!$       integer, dimension(:,:), intent(in)        :: A
+!!$       integer, dimension(:,:), intent(in)        :: B
+!!$       ! Returned variable
+!!$       integer, dimension( size(A,1), size(A,2) ) :: Ac
+!!$     end function morph_close
+!!$  end interface
+
+contains
+
+   subroutine morph_check_range(x, dx, minx, maxx, xmdx, xpdx, o0, o1)
+
+      implicit none
+
+      ! Input variables
+      integer, intent(in)  :: x, dx, minx, maxx
+      ! Output variables
+      integer, intent(out) :: xmdx, xpdx, o0, o1
+
+      if (x-dx .lt. minx) then
+         xmdx = minx
+         o0 = 1 + minx - (x-dx)
+      else
+         xmdx = x-dx
+         o0 = 1
+      end if
+      if (x+dx .gt. maxx) then
+         xpdx = maxx
+         o1 = 2*dx + 1 - (x+dx) + maxx
+      else
+         xpdx = x+dx
+         o1 = 2*dx + 1
+      end if
+
+   end subroutine morph_check_range
+
+   function morph_erode(A, B) result(Ae)
+
+      implicit none
+
+      ! Input variables
+      integer, dimension(:,:), intent(in)        :: A
+      integer, dimension(:,:), intent(in)        :: B
+      ! Returned variable
+      integer, dimension( size(A,1), size(A,2) ) :: Ae
+      ! Local variables
+      ! Array dimensions
+      integer                                    :: nx, ny, mx, my
+      ! Indices in x-direction
+      integer                                    :: i, i0, i1, bi0, bi1
+      ! Indices in y-direction
+      integer                                    :: j, j0, j1, bj0, bj1
+      ! Temporary array combining input and kernel
+      integer, dimension( size(B,1), size(B,2) ) :: Ab
+
+      ! Check that B has odd dimensions
+      if ((MOD(size(B,1), 2) .ne. 1) .or. (MOD(size(B,2), 2) .ne. 1)) then
+         write(*,*) 'ERROR: morph_erode(): Kernel array, B, must have odd dimensions'
+         stop 1
+      end if
+
+      ! Define dimension variables
+      nx = size(A,1)
+      ny = size(A,2)
+      mx = size(B,1) / 2
+      my = size(B,2) / 2
+
+      ! Now step through the A array, applying the erode transform
+      do i = 1, nx
+         do j = 1, ny
+            call morph_check_range(i, mx, 1, nx, i0, i1, bi0, bi1)
+            call morph_check_range(j, my, 1, ny, j0, j1, bj0, bj1)
+            Ab(:,:) = 0
+            Ab(bi0:bi1, bj0:bj1) = A(i0:i1, j0:j1) * B(bi0:bi1, bj0:bj1)
+            Ae(i, j) = minval(Ab, mask = B .gt. 0)
+         end do
+      end do
+
+   end function morph_erode
+
+   function morph_dilate(A, B) result(Ad)
+
+      implicit none
+      ! Input variables
+      integer, dimension(:,:), intent(in)        :: A
+      integer, dimension(:,:), intent(in)        :: B
+      ! Returned variable
+      integer, dimension( size(A,1), size(A,2) ) :: Ad
+      ! Local variables
+      ! Array dimensions
+      integer                                    :: nx, ny, mx, my
+      ! Indices in x-direction
+      integer                                    :: i, i0, i1, bi0, bi1
+      ! Indices in y-direction
+      integer                                    :: j, j0, j1, bj0, bj1
+      ! Temporary array combining input and kernel
+      integer, dimension( size(B,1), size(B,2) ) :: Ab
+
+      ! Check that B has odd dimensions
+      if ((MOD(size(B,1), 2) .ne. 1) .or. (MOD(size(B,2), 2) .ne. 1)) then
+         write(*,*) 'ERROR: morph_dilate(): Kernel array, B, must have odd dimensions'
+         stop 1
+      end if
+
+      ! Define dimension variables
+      nx = size(A,1)
+      ny = size(A,2)
+      mx = size(B,1) / 2
+      my = size(B,2) / 2
+
+      ! Now step through the A array, applying the dilate transform
+      do i = 1, nx
+         do j = 1, ny
+            call morph_check_range(i, mx, 1, nx, i0, i1, bi0, bi1)
+            call morph_check_range(j, my, 1, ny, j0, j1, bj0, bj1)
+            Ab(:,:) = 0
+            Ab(bi0:bi1, bj0:bj1) = A(i0:i1, j0:j1) * B(bi0:bi1, bj0:bj1)
+            Ad(i, j) = maxval(Ab, mask = B .gt. 0)
+         end do
+      end do
+
+   end function morph_dilate
+
+   function morph_open(A, B) result(Ao)
+      implicit none
+      ! Input variables
+      integer, dimension(:,:), intent(in)        :: A
+      integer, dimension(:,:), intent(in)        :: B
+      ! Returned variable
+      integer, dimension( size(A,1), size(A,2) ) :: Ao
+
+      ! Simply call dilate on the erode of A, using the kernel B
+      Ao = morph_dilate( morph_erode(A, B), B)
+
+   end function morph_open
+
+   function morph_close(A, B) result(Ac)
+      implicit none
+      ! Input variables
+      integer, dimension(:,:), intent(in)        :: A
+      integer, dimension(:,:), intent(in)        :: B
+      ! Returned variable
+      integer, dimension( size(A,1), size(A,2) ) :: Ac
+
+      ! Simply call dilate on the erode of A, using the kernel B
+      Ac = morph_erode( morph_dilate(A, B), B)
+
+   end function morph_close
+
+
+end module morphology_m

common/orac_output.F90

@@ -405,7 +405,7 @@ module orac_output_m
       integer(byte) :: cldtype_scale  = 1
       integer(byte) :: cldtype_offset = 0
       integer(byte) :: cldtype_vmin   = 0
-      integer(byte) :: cldtype_vmax   = 9
+      integer(byte) :: cldtype_vmax   = 12
 
       integer(byte) :: cldmask_scale  = 1
       integer(byte) :: cldmask_offset = 0

config/lib.conda.inc

@@ -3,7 +3,7 @@
 
 # Set up libraries and includes
 LIBS = -L$(CONDA_PREFIX)/lib \
-       -lemosR64 -lemos -lfftw3 \
+       -lfftw3 \
        -lhdfeos -lGctp \
        -leccodes_f90 -leccodes \
        -lmfhdf -lmfhdf_fortran -lmfhdf_fcstub -lhdf \
@@ -38,6 +38,10 @@ CINC += -DINCLUDE_ATSR_SUPPORT
 LIBS += -lnr
 INC  += -DINCLUDE_NR
 
+# Uncomment if LIBEMOS is available for meteorological field interpolation.
+#LIBS += -lemosR64 -lemos
+#INC  += -DINCLUDE_EMOS
+
 # Uncomment if Fu_Liou support is desired for broadband fluxes.
 LIBS += -lEd3Fu_201212
 INC  += -DINCLUDE_FU_LIOU_SUPPORT
@@ -52,5 +56,5 @@ INC  += -DINCLUDE_SEVIRI_SUPPORT
 
 # Uncomment if SEVIRI_ML support is desired.
 # Please modify Python -lpython3.X to your Python version.
-LIBS += -lsevann -lpython3.10
-INC  += -DINCLUDE_SEVIRI_NEURALNET
+#LIBS += -lsevann -lpython3.10
+#INC  += -DINCLUDE_SEVIRI_NEURALNET

config/lib.inc

@@ -59,7 +59,7 @@ SEVIRI_UTIL_LIB = $(LIBBASE_FORTRAN)/seviri_util
 SEVIRI_UTIL_INCLUDE = $(LIBBASE_FORTRAN)/seviri_util
 
 # seviri_ml
-# run seviri_ml/get_py_config.sh to obtain your SEVIRI_ML_PYTHON_LIB 
+# run seviri_ml/get_py_config.sh to obtain your SEVIRI_ML_PYTHON_LIB
 # and SEVIRI_ML_PYTHON_INCLUDE
 SEVIRI_ML_LIB = $(LIBBASE_FORTRAN)/seviri_ml
 SEVIRI_ML_INCLUDE = $(LIBBASE_FORTRAN)/seviri_ml
@@ -69,8 +69,7 @@ SEVIRI_ML_PYTHON_INCLUDE = /path/to/your/python/include
 SZLIB = $(LIBBASE)/szip/lib
 
 # Set up libraries and includes
-LIBS = -L$(EMOSLIB) -lemosR64 -lemos -lfftw3 \
-       -L$(EOSLIB) -lhdfeos -lGctp \
+LIBS = -L$(EOSLIB) -lhdfeos -lGctp \
        -L$(GRIBLIB) -leccodes_f90 -leccodes \
        -L$(HDFLIB) -lmfhdf -ldf \
        -L$(NCDF_FORTRAN_LIB) -lnetcdff \
@@ -109,6 +108,10 @@ CINC = -I./
 #LIBS += -L$(NRLIB) -lnr
 #INC  += -I$(NRINCLUDE) -DINCLUDE_NR
 
+# Uncomment if LIBEMOS is available for meteorological field interpolation.
+#LIBS += -L$(EMOSLIB) -lemosR64 -lemos -lfftw3
+#INC  += -DINCLUDE_EMOS
+
 # Uncomment if Fu_Liou support is desired for broadband fluxes.
 #LIBS += -L$(FULIOULIB) -lEd3Fu_201212
 #INC  += -I$(FULIOUINCLUDE) -DINCLUDE_FU_LIOU_SUPPORT

derived_products/broadband_fluxes/dependencies.inc

@@ -13,7 +13,7 @@ $(OBJS)/two_rt_lw.o: $(OBJS)/kinds.o
 $(OBJS)/two_rt_lw_iter.o: $(OBJS)/kinds.o
 $(OBJS)/two_rt_lw_ocastrndm.o: $(OBJS)/kinds.o
 $(OBJS)/two_rt_lw_sel.o: $(OBJS)/kinds.o
-$(OBJS)/two_rt_sw_bs.o: $(OBJS)/kinds.o $(OBJS)/bandsolve.o
+$(OBJS)/two_rt_sw_bs.o: $(OBJS)/bandsolve.o $(OBJS)/kinds.o
 $(OBJS)/two_rt_sw.o: $(OBJS)/kinds.o
 $(OBJS)/two_rt_sw_ocastrndm.o: $(OBJS)/kinds.o
 $(OBJS)/comscp1.o: $(OBJS)/kinds.o

pre_processing/Makefile

@@ -51,6 +51,7 @@ PREPROC_OBJ_F90 = $(OBJS)/aatsr_corrections.o \
                   $(OBJS)/cloud_emis.o \
                   $(OBJS)/cloud_typing_pavolonis.o \
                   $(OBJS)/correct_for_ice_snow.o \
+                  $(OBJS)/correct_for_dust.o \
                   $(OBJS)/cox_munk.o \
                   $(OBJS)/cox_munk_constants.o \
                   $(OBJS)/define_preprop_grid.o \