diff --git a/common/Makefile b/common/Makefile index e5546f09..e5c6328b 100644 --- a/common/Makefile +++ b/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 \ diff --git a/common/constants_cloud_typing_pavolonis.F90 b/common/constants_cloud_typing_pavolonis.F90 index 195706ec..84503765 100644 --- a/common/constants_cloud_typing_pavolonis.F90 +++ b/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 diff --git a/common/def_output_primary.F90 b/common/def_output_primary.F90 index f9a45688..1478099d 100644 --- a/common/def_output_primary.F90 +++ b/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) diff --git a/common/morphology.F90 b/common/morphology.F90 new file mode 100644 index 00000000..13358c89 --- /dev/null +++ b/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 diff --git a/common/orac_output.F90 b/common/orac_output.F90 index 99c43785..4f05a152 100644 --- a/common/orac_output.F90 +++ b/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 diff --git a/config/lib.conda.inc b/config/lib.conda.inc index 3bc566c7..1fc86293 100644 --- a/config/lib.conda.inc +++ b/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 diff --git a/config/lib.inc b/config/lib.inc index 881ad7d0..0d15ee18 100644 --- a/config/lib.inc +++ b/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 diff --git a/derived_products/broadband_fluxes/dependencies.inc b/derived_products/broadband_fluxes/dependencies.inc index aa6bb22a..302d61a4 100644 --- a/derived_products/broadband_fluxes/dependencies.inc +++ b/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 diff --git a/pre_processing/Makefile b/pre_processing/Makefile index c16ce189..0d209a6e 100644 --- a/pre_processing/Makefile +++ b/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 \ diff --git a/pre_processing/allocate_preproc_structures.F90 b/pre_processing/allocate_preproc_structures.F90 index 49444c41..32465581 100644 --- a/pre_processing/allocate_preproc_structures.F90 +++ b/pre_processing/allocate_preproc_structures.F90 @@ -40,7 +40,7 @@ ! preproc_lwrtm%phi_lay and preproc_lwrtm%phi_lev ! 2013/10/16, CP: removed kdim_pre and removed ecmwf_dims form subroutine call, ! hanged array definition of players and plev -! 2013/11/08, GM: Removed double allocations (allocated twice in a row) of +! 2013/11/08, GM: Removed double allocations (associated twice in a row) of ! preproc_prtm%temperature through preproc_prtm%ozone. ! 2014/02/10, AP: Simplifying variable names ! 2014/05/01, GM: Add some allocations that were being done outside. @@ -53,6 +53,7 @@ ! 2017/03/29, SP: Add ability to calculate tropospheric cloud emissivity (ExtWork) ! 2017/11/15, SP: Add feature to give access to sensor azimuth angle ! 2018/07/18, DE: Add tropoopause temperature +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. @@ -67,67 +68,60 @@ subroutine allocate_preproc_prtm(preproc_dims, preproc_prtm) type(preproc_dims_t), intent(in) :: preproc_dims type(preproc_prtm_t), intent(out) :: preproc_prtm - integer :: sx, ex, sy, ey - - sx=preproc_dims%min_lon - ex=preproc_dims%max_lon - sy=preproc_dims%min_lat - ey=preproc_dims%max_lat - - allocate(preproc_prtm%pressure(sx:ex,sy:ey,preproc_dims%kdim)) + allocate(preproc_prtm%pressure(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim)) preproc_prtm%pressure=sreal_fill_value - allocate(preproc_prtm%temperature(sx:ex,sy:ey,preproc_dims%kdim)) + allocate(preproc_prtm%temperature(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim)) preproc_prtm%temperature=sreal_fill_value - allocate(preproc_prtm%spec_hum(sx:ex,sy:ey,preproc_dims%kdim)) + allocate(preproc_prtm%spec_hum(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim)) preproc_prtm%spec_hum=sreal_fill_value - allocate(preproc_prtm%ozone(sx:ex,sy:ey,preproc_dims%kdim)) + allocate(preproc_prtm%ozone(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim)) preproc_prtm%ozone=sreal_fill_value - allocate(preproc_prtm%phi_lay(sx:ex,sy:ey,preproc_dims%kdim)) + allocate(preproc_prtm%phi_lay(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim)) preproc_prtm%phi_lay=sreal_fill_value - allocate(preproc_prtm%phi_lev(sx:ex,sy:ey,preproc_dims%kdim+1)) + allocate(preproc_prtm%phi_lev(preproc_dims%xdim,preproc_dims%ydim,preproc_dims%kdim+1)) preproc_prtm%phi_lev=sreal_fill_value - allocate(preproc_prtm%geopot(sx:ex,sy:ey)) + allocate(preproc_prtm%geopot(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%geopot=sreal_fill_value - allocate(preproc_prtm%lnsp(sx:ex,sy:ey)) + allocate(preproc_prtm%lnsp(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%lnsp=sreal_fill_value - allocate(preproc_prtm%u10(sx:ex,sy:ey)) + allocate(preproc_prtm%u10(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%u10=sreal_fill_value - allocate(preproc_prtm%v10(sx:ex,sy:ey)) + allocate(preproc_prtm%v10(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%v10=sreal_fill_value - allocate(preproc_prtm%land_sea_mask(sx:ex,sy:ey)) + allocate(preproc_prtm%land_sea_mask(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%land_sea_mask=sreal_fill_value - allocate(preproc_prtm%temp2(sx:ex,sy:ey)) + allocate(preproc_prtm%temp2(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%temp2=sreal_fill_value - allocate(preproc_prtm%skin_temp(sx:ex,sy:ey)) + allocate(preproc_prtm%skin_temp(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%skin_temp=sreal_fill_value - allocate(preproc_prtm%snow_albedo(sx:ex,sy:ey)) + allocate(preproc_prtm%snow_albedo(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%snow_albedo=sreal_fill_value - allocate(preproc_prtm%snow_depth(sx:ex,sy:ey)) + allocate(preproc_prtm%snow_depth(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%snow_depth=sreal_fill_value - allocate(preproc_prtm%sst(sx:ex,sy:ey)) + allocate(preproc_prtm%sst(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%sst=sreal_fill_value - allocate(preproc_prtm%sea_ice_cover(sx:ex,sy:ey)) - preproc_prtm%sea_ice_cover=sreal_fill_value + allocate(preproc_prtm%sea_ice_cover(preproc_dims%xdim,preproc_dims%ydim)) + preproc_prtm%sea_ice_cover=sreal_fill_value - allocate(preproc_prtm%totcolwv(sx:ex,sy:ey)) + allocate(preproc_prtm%totcolwv(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%totcolwv=sreal_fill_value #ifdef INCLUDE_SATWX - allocate(preproc_prtm%trop_p(sx:ex,sy:ey)) + allocate(preproc_prtm%trop_p(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%trop_p=sreal_fill_value - allocate(preproc_prtm%trop_t(sx:ex,sy:ey)) + allocate(preproc_prtm%trop_t(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%trop_t=sreal_fill_value - allocate(preproc_prtm%cape(sx:ex,sy:ey)) + allocate(preproc_prtm%cape(preproc_dims%xdim,preproc_dims%ydim)) preproc_prtm%cape=sreal_fill_value #endif @@ -153,50 +147,46 @@ subroutine allocate_preproc_structures(imager_angles,preproc_dims, & type(preproc_cld_t), intent(out) :: preproc_cld type(channel_info_t), intent(inout) :: channel_info - integer :: nchan_sw, nchan_lw, sx, ex, sy, ey + integer :: nchan_sw, nchan_lw nchan_sw=channel_info%nchannels_sw nchan_lw=channel_info%nchannels_lw - sx=preproc_dims%min_lon - ex=preproc_dims%max_lon - sy=preproc_dims%min_lat - ey=preproc_dims%max_lat ! preproc_dims - allocate(preproc_dims%counter_sw(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_dims%counter_sw(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_dims%counter_sw=0 - allocate(preproc_dims%counter_lw(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_dims%counter_lw(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_dims%counter_lw=0 ! preproc_geoloc - allocate(preproc_geoloc%longitude(sx:ex)) + allocate(preproc_geoloc%longitude(preproc_dims%xdim)) preproc_geoloc%longitude=sreal_fill_value - allocate(preproc_geoloc%latitude(sy:ey)) + allocate(preproc_geoloc%latitude(preproc_dims%ydim)) preproc_geoloc%latitude=sreal_fill_value ! preproc_geo (init to 0 as used for summation in build_preproc_fields) - allocate(preproc_geo%solza(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_geo%solza(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_geo%solza=0.0 - allocate(preproc_geo%solazi(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_geo%solazi(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_geo%solazi=0.0 - allocate(preproc_geo%satza(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_geo%satza(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_geo%satza=0.0 - allocate(preproc_geo%satazi(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_geo%satazi(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_geo%satazi=0.0 - allocate(preproc_geo%relazi(sx:ex,sy:ey,imager_angles%nviews)) + allocate(preproc_geo%relazi(preproc_dims%xdim,preproc_dims%ydim,imager_angles%nviews)) preproc_geo%relazi=0.0 ! preproc_prtm call allocate_preproc_prtm(preproc_dims, preproc_prtm) ! preproc_surf - allocate(preproc_surf%emissivity(sx:ex,sy:ey,nchan_lw)) + allocate(preproc_surf%emissivity(preproc_dims%xdim,preproc_dims%ydim,nchan_lw)) preproc_surf%emissivity=sreal_fill_value ! preproc_cld - allocate(preproc_cld%clear_bt(sx:ex,sy:ey,nchan_lw)) + allocate(preproc_cld%clear_bt(preproc_dims%xdim,preproc_dims%ydim,nchan_lw)) preproc_cld%clear_bt=0.0 - allocate(preproc_cld%cloud_bt(sx:ex,sy:ey,nchan_lw)) + allocate(preproc_cld%cloud_bt(preproc_dims%xdim,preproc_dims%ydim,nchan_lw)) preproc_cld%cloud_bt=0.0 end subroutine allocate_preproc_structures diff --git a/pre_processing/build_preproc_fields.F90 b/pre_processing/build_preproc_fields.F90 index 43ac220a..9625ff3f 100644 --- a/pre_processing/build_preproc_fields.F90 +++ b/pre_processing/build_preproc_fields.F90 @@ -33,17 +33,20 @@ ! 2015/21/01, OS: bug fix in setting lon_i/lat_i min/max limits ! 2015/01/30, AP: Remove uscan and vscan as unnecessary. ! 2017/11/15, SP: Add feature to give access to sensor azimuth angle +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions and +! use native grid of ECMWF ! ! Bugs: ! None known. !------------------------------------------------------------------------------- subroutine build_preproc_fields(preproc_dims, preproc_geoloc, preproc_geo, & - imager_geolocation, imager_angles) + imager_geolocation, imager_angles, ecmwf, use_ecmwf_preproc_grid) use imager_structures_m use preproc_constants_m use preproc_structures_m + use ecmwf_m implicit none @@ -52,27 +55,34 @@ subroutine build_preproc_fields(preproc_dims, preproc_geoloc, preproc_geo, & type(preproc_geo_t), intent(inout) :: preproc_geo type(imager_geolocation_t), intent(inout) :: imager_geolocation type(imager_angles_t), intent(inout) :: imager_angles - + type(ecmwf_t), intent(inout) :: ecmwf + logical, intent(inout) :: use_ecmwf_preproc_grid integer(kind=lint) :: i,j,k,lon_i,lat_j real(sreal) :: fac - ! build the arrays for the regular grid - ! create grid resolution lat - fac = 1. / preproc_dims%dellat - preproc_geoloc%latitude(preproc_dims%min_lat) = & - (preproc_dims%min_lat-0.5)*fac - real(preproc_dims%lat_offset,sreal) - do i = preproc_dims%min_lat+1, preproc_dims%max_lat - preproc_geoloc%latitude(i) = preproc_geoloc%latitude(i-1) + fac - end do - ! create grid resolution lon - fac = 1. / preproc_dims%dellon - preproc_geoloc%longitude(preproc_dims%min_lon) = & - (preproc_dims%min_lon-0.5)*fac - real(preproc_dims%lon_offset,sreal) - do i = preproc_dims%min_lon+1, preproc_dims%max_lon - preproc_geoloc%longitude(i) = preproc_geoloc%longitude(i-1) + fac - end do + ! build the arrays for the regular grid + if (use_ecmwf_preproc_grid) then + preproc_geoloc%latitude = ecmwf%lat(preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + preproc_geoloc%longitude = ecmwf%lon(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind) + else + ! create grid resolution lat + fac = 1. / preproc_dims%dellat + preproc_geoloc%latitude(1) = & + (preproc_dims%min_lat-0.5)*fac - real(preproc_dims%lat_offset,sreal) + do i = 2, preproc_dims%ydim + preproc_geoloc%latitude(i) = preproc_geoloc%latitude(i-1) + fac + end do + + ! create grid resolution lon + fac = 1. / preproc_dims%dellon + preproc_geoloc%longitude(1) = & + (preproc_dims%min_lon-0.5)*fac - real(preproc_dims%lon_offset,sreal) + do i = 2, preproc_dims%xdim + preproc_geoloc%longitude(i) = preproc_geoloc%longitude(i-1) + fac + end do + end if ! imager resolution is always higher than preprocessing resolution ! =>average imager properties to this coarser resolution grid. @@ -89,16 +99,13 @@ subroutine build_preproc_fields(preproc_dims, preproc_geoloc, preproc_geo, & imager_geolocation%longitude(i,j) .eq. sreal_fill_value) cycle ! find grid cell coordinates into which L1b pixel falls - lon_i = floor((imager_geolocation%longitude(i,j) + & - preproc_dims%lon_offset)*preproc_dims%dellon, kind=lint) + 1 - lat_j = floor((imager_geolocation%latitude(i,j) + & - preproc_dims%lat_offset)*preproc_dims%dellat, kind=lint) + 1 - - if (lon_i .lt. preproc_dims%min_lon) lon_i = preproc_dims%min_lon - if (lat_j .lt. preproc_dims%min_lat) lat_j = preproc_dims%min_lat - if (lon_i .gt. preproc_dims%max_lon) lon_i = preproc_dims%max_lon - if (lat_j .gt. preproc_dims%max_lat) lat_j = preproc_dims%max_lat + lat_j = minloc(abs(preproc_geoloc%latitude - imager_geolocation%latitude(i,j)),1) + lon_i = minloc(abs(preproc_geoloc%longitude - imager_geolocation%longitude(i,j)),1) + if (imager_geolocation%longitude(i,j) .lt. minval(preproc_geoloc%longitude)) lon_i = 1 + if (imager_geolocation%latitude(i,j) .lt. minval(preproc_geoloc%latitude)) lat_j = 1 + if (imager_geolocation%longitude(i,j) .gt. maxval(preproc_geoloc%longitude)) lon_i = preproc_dims%xdim + if (imager_geolocation%latitude(i,j) .gt. maxval(preproc_geoloc%latitude)) lat_j = preproc_dims%ydim do k = 1, imager_angles%nviews if (imager_angles%satzen(i,j,k) .ne. sreal_fill_value) then preproc_geo%satza(lon_i,lat_j,k) = & @@ -132,8 +139,8 @@ subroutine build_preproc_fields(preproc_dims, preproc_geoloc, preproc_geo, & end do ! loop over preprocessor data i.e reduced resolution - do j = preproc_dims%min_lat, preproc_dims%max_lat - do i = preproc_dims%min_lon, preproc_dims%max_lon + do j = 1, preproc_dims%ydim + do i = 1, preproc_dims%xdim do k = 1, imager_angles%nviews if (preproc_dims%counter_lw(i,j,k) .gt. 0) then ! if this is a good preprocessing pixel, calculate the average diff --git a/pre_processing/cloud_emis.F90 b/pre_processing/cloud_emis.F90 index cc61cd9d..64bf35cf 100644 --- a/pre_processing/cloud_emis.F90 +++ b/pre_processing/cloud_emis.F90 @@ -14,6 +14,7 @@ ! 2018/07/18, DE: Add tropopause temperature ! 2018/11/05, SP: Add CAPE ! 2019/08/14, SP: Add Fengyun-4A support. +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. @@ -38,7 +39,7 @@ subroutine get_trop_tp(preproc_prtm, preproc_dims) real, parameter :: max_tropopause = 450.0 ! Lowest p allowed for trop integer, parameter :: depth = 2 ! # layers added to inversions - integer :: nx, ny, nz ! Number of vertical levels, pixels + integer :: nz ! Number of vertical levels, pixels integer :: x, y ! Looping variables over preproc integer :: k, l ! Indexing variables integer :: step ! Direction of search @@ -47,13 +48,11 @@ subroutine get_trop_tp(preproc_prtm, preproc_dims) real, dimension(preproc_dims%kdim) :: h ! Height profile - nx = preproc_dims%xdim - ny = preproc_dims%ydim nz = preproc_dims%kdim - do x = preproc_dims%min_lon, preproc_dims%max_lon - do y = preproc_dims%min_lat, preproc_dims%max_lat + do x = 1, preproc_dims%xdim + do y = 1, preproc_dims%ydim k = nz t = preproc_prtm%temperature(x,y,:) @@ -103,7 +102,7 @@ subroutine get_trop_tp(preproc_prtm, preproc_dims) end subroutine get_trop_tp subroutine get_cloud_emis(channel_info, imager_measurements, imager_geolocation, & - preproc_dims, preproc_geoloc, preproc_cld, preproc_prtm, imager_cloud, ecmwf, & + preproc_dims, preproc_geoloc, preproc_cld, preproc_prtm, imager_cloud, & sensor, verbose) use channel_structures_m @@ -127,7 +126,6 @@ subroutine get_cloud_emis(channel_info, imager_measurements, imager_geolocation, type(preproc_cld_t), intent(in) :: preproc_cld type(preproc_prtm_t), intent(in) :: preproc_prtm type(imager_cloud_t), intent(out) :: imager_cloud - type(ecmwf_t), intent(in) :: ecmwf character(len=*), intent(in) :: sensor logical, intent(in) :: verbose @@ -188,12 +186,12 @@ subroutine get_cloud_emis(channel_info, imager_measurements, imager_geolocation, allocate(interp(1)) ! Needed for grid interpolation, adapted from ../src/ReadPRTM_nc.F90 - NLat = abs(preproc_dims%min_lat - preproc_dims%max_lat) - NLon = abs(preproc_dims%min_lon - preproc_dims%max_lon) - Lat0 = real(preproc_geoloc%latitude(preproc_dims%min_lat), kind=8) - LatN = real(preproc_geoloc%latitude(preproc_dims%max_lat), kind=8) - Lon0 = real(preproc_geoloc%longitude(preproc_dims%min_lon), kind=8) - LonN = real(preproc_geoloc%longitude(preproc_dims%max_lon), kind=8) + NLat = preproc_dims%ydim + NLon = preproc_dims%xdim + Lat0 = real(preproc_geoloc%latitude(1), kind=8) + LatN = real(preproc_geoloc%latitude(preproc_dims%ydim), kind=8) + Lon0 = real(preproc_geoloc%longitude(1), kind=8) + LonN = real(preproc_geoloc%longitude(preproc_dims%xdim), kind=8) ! Grid spacing and inverse delta_Lat = (LatN - Lat0) / (NLat-1) diff --git a/pre_processing/compute_geopot_coordinate.F90 b/pre_processing/compute_geopot_coordinate.F90 index 36e397ef..59865ee0 100644 --- a/pre_processing/compute_geopot_coordinate.F90 +++ b/pre_processing/compute_geopot_coordinate.F90 @@ -21,6 +21,7 @@ ! 2014/05/08, AP: Updated to new ecmwf structure. ! 2019/05/23, GT: Added check for valid data in each grid cell before doing ! calculation +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. @@ -41,57 +42,57 @@ subroutine compute_geopot_coordinate(preproc_prtm, preproc_dims, ecmwf) real(kind=sreal) :: virt_temp,p,pp1,logpp,r_ratio,alpha,sp real(kind=sreal) :: sum_term,add_term - r_ratio=r_water_vap/(r_dry_air-1.0_sreal) + r_ratio = r_water_vap / (r_dry_air - 1.0_sreal) ! compute the summation terms of the sum in (2.21) and necessary terms in ! (2.22) & (2.23) from TOA down (index ik represents cell centers and cell ! upper boundaries (wrt height)) - do ij=preproc_dims%min_lat,preproc_dims%max_lat - do ii=preproc_dims%min_lon,preproc_dims%max_lon + do ij = 1, preproc_dims%ydim + do ii = 1, preproc_dims%xdim ! Check to see we have data for this particular pixel if (preproc_prtm%lnsp(ii,ij) .ne. sreal_fill_value) then ! this is the lowest level=surface, it also has the surface pressure - preproc_prtm%phi_lev(ii,ij,ecmwf%kdim+1)=preproc_prtm%geopot(ii,ij) - sp=exp(preproc_prtm%lnsp(ii,ij)) + preproc_prtm%phi_lev(ii,ij,ecmwf%kdim+1) = preproc_prtm%geopot(ii,ij) + sp = exp(preproc_prtm%lnsp(ii,ij)) !pressure at cell lower boundary - pp1=ecmwf%avec(ecmwf%kdim+1)+ecmwf%bvec(ecmwf%kdim+1)*sp + pp1 = ecmwf%avec(ecmwf%kdim+1) + ecmwf%bvec(ecmwf%kdim+1) * sp ! sum from surface up according to (2.21) - do ik=ecmwf%kdim,1,-1 + do ik = ecmwf%kdim, 1, -1 !pressure at cell upper boundary - p=ecmwf%avec(ik)+ecmwf%bvec(ik)*sp - preproc_prtm%pressure(ii,ij,ik)=0.5*(p+pp1) + p = ecmwf%avec(ik) + ecmwf%bvec(ik) * sp + preproc_prtm%pressure(ii,ij,ik) = 0.5 * (p + pp1) !logpp is logarithmic pressure difference, defined on cell centers if (p .gt. dither) then - logpp=log(pp1/p) + logpp = log(pp1/p) else !TOA has zero pressure, therefore: - logpp=log(pp1) + logpp = log(pp1) end if !virtual temperature at cell centers - virt_temp=preproc_prtm%temperature(ii,ij,ik)*(1.0_sreal + & - r_ratio*preproc_prtm%spec_hum(ii,ij,ik)) - sum_term=r_dry_air*virt_temp*logpp + virt_temp = preproc_prtm%temperature(ii,ij,ik) * (1.0_sreal + & + r_ratio * preproc_prtm%spec_hum(ii,ij,ik)) + sum_term = r_dry_air * virt_temp * logpp !alpha term used later to put gph on cell centers, s.b. if (ik .eq. 1) then - alpha=log(2.0_sreal) + alpha = log(2.0_sreal) else - alpha=1.0_sreal-p/(pp1-p)*logpp + alpha = 1.0_sreal - p / (pp1 - p) * logpp end if !add_term dito to alpha term, s.b. - add_term=alpha*r_dry_air*virt_temp + add_term = alpha*r_dry_air*virt_temp ! perform sum - preproc_prtm%phi_lev(ii,ij,ik)=preproc_prtm%phi_lev(ii,ij,ik+1) + & + preproc_prtm%phi_lev(ii,ij,ik) = preproc_prtm%phi_lev(ii,ij,ik+1) + & sum_term - preproc_prtm%phi_lay(ii,ij,ik)=preproc_prtm%phi_lev(ii,ij,ik+1) + & + preproc_prtm%phi_lay(ii,ij,ik) = preproc_prtm%phi_lev(ii,ij,ik+1) + & add_term - pp1=p + pp1 = p end do end if end do diff --git a/pre_processing/correct_for_dust.F90 b/pre_processing/correct_for_dust.F90 new file mode 100644 index 00000000..4c0df3a8 --- /dev/null +++ b/pre_processing/correct_for_dust.F90 @@ -0,0 +1,105 @@ +!------------------------------------------------------------------------------- +! Name: correct_for_dust.F90 +! +! Purpose: +! This is a very simple correction for detecting areas masked as cloud +! which are actually dust. Might also be effective for volcanic ash. +! It selects pixels which have been flagged as cloud, but with a high +! uncertainty, and then uses the 11-12 micron BTD to flag possible +! dust. Then, a morphological opening transform is applied to the +! resulting dust mask to remove random false detections. +! +! History: +! 2024/03/11, GT: Initial version +! +! Bugs: +! None known. +!------------------------------------------------------------------------------- + +subroutine correct_for_dust(channel_info, imager_measurements, imager_angles, & + imager_geolocation, imager_pavolonis, verbose) + + use common_constants_m + use channel_structures_m + use imager_structures_m + use morphology_m + + ! Input (and output) variables + type(channel_info_t), intent(in) :: channel_info + type(imager_measurements_t), intent(in) :: imager_measurements + type(imager_angles_t), intent(in) :: imager_angles + type(imager_geolocation_t), intent(in) :: imager_geolocation + type(imager_pavolonis_t), intent(inout) :: imager_pavolonis + logical, intent(in) :: verbose + ! Local variables + integer(kind=lint) :: i, j + integer(kind=lint) :: bt11=-1, bt12=-1 + integer, dimension( imager_geolocation%startx:imager_geolocation%endx, & + 1:imager_geolocation%ny ) :: dustmask + integer, dimension(5,5) :: kernel=reshape( (/ 0, 1, 1, 1, 0, & + 1, 1, 1, 1, 1, & + 1, 1, 1, 1, 1, & + 1, 1, 1, 1, 1, & + 0, 1, 1, 1, 0 /), & + shape(kernel) ) + + ! First we need to determine if we have the required channels: + do i = 1, channel_info%nchannels_total + if (ANINT(channel_info%channel_wl_abs(i)) .eq. 11.0) bt11 = i + if (ANINT(channel_info%channel_wl_abs(i)) .eq. 12.0) bt12 = i + if (bt11 .gt. 0 .and. bt12 .gt. 0) exit + end do + + ! If we do have the requisite brightness temperatures, then do our test + ! The tests are: + ! * Is the result from the NN cloud-mask uncertain? + ! * Do we have valid BTs? + ! * BT difference test + ! * Limit on solar-zenith. Deals with false positives at high latitudes + ! in LEO instruments (need to check if this is suitable for Geo). + if (bt11 .gt. 0 .and. bt12 .gt. 0) then + where(imager_pavolonis%cldmask_uncertainty(:,:,1) .gt. 10 .and. & + imager_measurements%data(:,:,bt11) .gt. 0 .and. & + (imager_measurements%data(:,:,bt11) - & + imager_measurements%data(:,:,bt12)) .lt. 0.3 .and. & + (imager_angles%solzen(:,:,1) .lt. 50 .or. & + abs(imager_geolocation%latitude) .lt. 40)) + dustmask = 1 + elsewhere + dustmask = 0 + end where + ! Now do perform an opening transform on our new mask, which should + ! remove scattered false positives + if (verbose) write(*,*) count(dustmask .eq. 1),' pixels flagged as dust before open' + dustmask = morph_open(dustmask, kernel) + + if (verbose) write(*,*) count(dustmask .eq. 1),' pixels flagged as dust' + ! Finally, correct the cldmask with the detected dust pixels and add + ! two new values to the Pavalonis cloud-type mask, one indicating where + ! we think a clear pixel is dust, and where previously flagged cloud + ! has been switched to dust + do i = imager_geolocation%startx, imager_geolocation%endx + do j = 1, size(imager_pavolonis%cldmask,2) + if (dustmask(i,j) .eq. 1) then + if (maxval(imager_pavolonis%cldmask(i,j,:)) .eq. 1) then + imager_pavolonis%cldtype(i,j,1) = DUST_SWITCHED_FROM_CLOUD_TYPE + else + imager_pavolonis%cldtype(i,j,1) = DUST_CLEAR_TYPE + end if + imager_pavolonis%cldmask(i,j,:) = CLEAR_TYPE + end if + end do + end do + if (verbose) write(*,*) & + count(imager_pavolonis%cldtype(:,:,1) .eq. DUST_CLEAR_TYPE), & + ' clear pixels set as dust' + if (verbose) write(*,*) & + count(imager_pavolonis%cldtype(:,:,1) .eq. DUST_SWITCHED_FROM_CLOUD_TYPE), & + ' cloud pixels set as dust' + else + ! We don't have the required BTs, so issue a warning and return without + ! doing anything + write(*,*) 'WARNING: correct_for_dust(): Need 11 and 12 micron BTs for dust correction. No correction performed.' + end if + +end subroutine correct_for_dust diff --git a/pre_processing/correct_for_ice_snow.F90 b/pre_processing/correct_for_ice_snow.F90 index 610726ba..46e1d90f 100644 --- a/pre_processing/correct_for_ice_snow.F90 +++ b/pre_processing/correct_for_ice_snow.F90 @@ -15,6 +15,7 @@ ! struct (defined in imager_structures) ! preproc_dims struct in Preprocessing dimensions, including sw and ! lw channel counts +! preproc_geoloc struct in Summary of preprocessing lat/lon ! surface struct both Surface properties structure ! cyear string in Year, as a 4 character string. ! cmonth string in Month of year, as a 2 character string. @@ -92,6 +93,7 @@ ! 2016/02/18, OS: ECMWF snow/ice mask now corrected by USGS land/sea mask ! 2016/02/23, OS: previous commit on ECMWF snow/ice mask was incomplete ! 2018/10/01, SP: Introduce a more comprehensive snow albedo dataset +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. @@ -586,7 +588,7 @@ end subroutine apply_ice_correction !------------------------------------------------------------------------------ subroutine correct_for_ice_snow_nwp(nwp_path, imager_geolocation, & - channel_info, imager_flags, preproc_dims, preproc_prtm, surface, & + channel_info, imager_flags, preproc_dims, preproc_prtm, preproc_geoloc, surface, & include_full_brdf, source_atts, verbose) use channel_structures_m @@ -606,6 +608,7 @@ subroutine correct_for_ice_snow_nwp(nwp_path, imager_geolocation, & type(imager_flags_t), intent(in) :: imager_flags type(preproc_dims_t), intent(in) :: preproc_dims type(preproc_prtm_t), intent(in) :: preproc_prtm + type(preproc_geoloc_t), intent(in) :: preproc_geoloc logical, intent(in) :: include_full_brdf type(source_attributes_t), intent(inout) :: source_atts logical, intent(in) :: verbose @@ -614,7 +617,7 @@ subroutine correct_for_ice_snow_nwp(nwp_path, imager_geolocation, & ! Local variables logical :: flag - integer(kind=4) :: i, j,lon_i, lat_j + integer(kind=4) :: i, j, lon_i, lat_j real(kind=sreal), dimension(:), allocatable :: tmp_albedo, tmp_snow, tmp_ice real(kind=sreal), dimension(4) :: snow_albedo, ice_albedo real(kind=sreal) :: snow_threshold, ice_threshold @@ -649,15 +652,12 @@ subroutine correct_for_ice_snow_nwp(nwp_path, imager_geolocation, & imager_geolocation%longitude(i,j) .eq. sreal_fill_value) cycle ! find grid cell coordinates into which L1b pixel falls - lon_i = floor((imager_geolocation%longitude(i,j) + & - preproc_dims%lon_offset)*preproc_dims%dellon, kind=lint) + 1 - lat_j = floor((imager_geolocation%latitude(i,j) + & - preproc_dims%lat_offset)*preproc_dims%dellat, kind=lint) + 1 - - if (lon_i .lt. preproc_dims%min_lon) lon_i = preproc_dims%min_lon - if (lat_j .lt. preproc_dims%min_lat) lat_j = preproc_dims%min_lat - if (lon_i .gt. preproc_dims%max_lon) lon_i = preproc_dims%max_lon - if (lat_j .gt. preproc_dims%max_lat) lat_j = preproc_dims%max_lat + lat_j = minloc(abs(preproc_geoloc%latitude - imager_geolocation%latitude(i,j)),1) + lon_i = minloc(abs(preproc_geoloc%longitude - imager_geolocation%longitude(i,j)),1) + if (imager_geolocation%longitude(i,j) .lt. minval(preproc_geoloc%longitude)) lon_i = 1 + if (imager_geolocation%latitude(i,j) .lt. minval(preproc_geoloc%latitude)) lat_j = 1 + if (imager_geolocation%longitude(i,j) .gt. maxval(preproc_geoloc%longitude)) lon_i = preproc_dims%xdim + if (imager_geolocation%latitude(i,j) .gt. maxval(preproc_geoloc%latitude)) lat_j = preproc_dims%ydim tmp_albedo = surface%albedo(i,j,:) diff --git a/pre_processing/define_preprop_grid.F90 b/pre_processing/define_preprop_grid.F90 index f9c9a9d5..84861460 100644 --- a/pre_processing/define_preprop_grid.F90 +++ b/pre_processing/define_preprop_grid.F90 @@ -56,13 +56,13 @@ subroutine define_preprop_grid(imager_geolocation,preproc_dims,verbose) 1:imager_geolocation%ny)) ! determine which preproc grid points each pixel falls in (start at 1) - lat = (imager_geolocation%latitude + preproc_dims%lat_offset)* & + lat = (imager_geolocation%latitude + preproc_dims%lat_offset) * & preproc_dims%dellat + 1. - lon = (imager_geolocation%longitude + preproc_dims%lon_offset)* & + lon = (imager_geolocation%longitude + preproc_dims%lon_offset) * & preproc_dims%dellon + 1. ! only consider valid lat/lon values - mask = imager_geolocation%latitude.ne.sreal_fill_value .and. & + mask = imager_geolocation%latitude.ne.sreal_fill_value .and. & imager_geolocation%longitude.ne.sreal_fill_value ! take one more pixel than is required for interpolation @@ -85,9 +85,9 @@ subroutine define_preprop_grid(imager_geolocation,preproc_dims,verbose) preproc_dims%ydim = preproc_dims%max_lat - preproc_dims%min_lat + 1 if (verbose) then - write(*,*) 'preproc_dims: ',preproc_dims%xdim, & - preproc_dims%ydim,preproc_dims%kdim - write(*,*) 'dellon, dellat: ',preproc_dims%dellon,preproc_dims%dellat + print*, 'preproc_dims: ', preproc_dims%xdim, & + preproc_dims%ydim, preproc_dims%kdim + print*, 'dellon, dellat: ',preproc_dims%dellon, preproc_dims%dellat print*, 'preproc_dims%min_lat: ', preproc_dims%min_lat print*, 'preproc_dims%max_lat: ', preproc_dims%max_lat print*, 'preproc_dims%min_lon: ', preproc_dims%min_lon @@ -99,3 +99,77 @@ subroutine define_preprop_grid(imager_geolocation,preproc_dims,verbose) deallocate(mask) end subroutine define_preprop_grid + +!------------------------------------------------------------------------------- +! Name: define_preprop_grid_ecmwf.F90 +! +! Purpose: +! Use the high resolution ecmwf grid for preprocessing without interpolation. +! +! Description and Algorithm details: +! 1) Determine minimum and maximum lat/lon of image from the ECMWF lat/lon grid. +! +! Arguments: +! Name Type In/Out/Both Description +! ------------------------------------------------------------------------------ +! imager_geolocation both Summary of pixel positions +! struct +! preproc_dims struct Out Structure summarising dimensions of preprocessing. +! ecmwf struct Both Summary of contents of ECMWF files +! verbose logic In F: minimise information printed to screen; T: don't +! +! History: +! 2024/07/01, DH: Initial version +! +! Bugs: +! None known. +!------------------------------------------------------------------------------- + +subroutine define_preproc_grid_ecmwf(imager_geolocation, preproc_dims, ecmwf, verbose) + + use preproc_constants_m + use preproc_structures_m + use ecmwf_m + use imager_structures_m + + implicit none + + type(ecmwf_t), intent(in) :: ecmwf + type(preproc_dims_t), intent(inout) :: preproc_dims + type(imager_geolocation_t), intent(in) :: imager_geolocation + logical, intent(in) :: verbose + + integer :: min_lon_ind, max_lon_ind, min_lat_ind, max_lat_ind + + ! get the ecmwf index of min and max lat and lon + min_lon_ind = minloc(abs(ecmwf%lon - minval(merge(imager_geolocation%longitude, HUGE(0.0), imager_geolocation%longitude /= -999))),1) + max_lon_ind = minloc(abs(ecmwf%lon - maxval(imager_geolocation%longitude)),1) + min_lat_ind = minloc(abs(ecmwf%lat - minval(merge(imager_geolocation%latitude, HUGE(0.0), imager_geolocation%latitude /= -999))),1) + max_lat_ind = minloc(abs(ecmwf%lat - maxval(imager_geolocation%latitude)),1) + + ! use ecmwf grid for preproc grid + preproc_dims%kdim = ecmwf%kdim + preproc_dims%xdim = max_lon_ind - min_lon_ind + 1 + preproc_dims%ydim = max_lat_ind - min_lat_ind + 1 + preproc_dims%min_lat = ecmwf%lat(min_lat_ind) + preproc_dims%min_lat_ind = min_lat_ind + preproc_dims%max_lat = ecmwf%lat(max_lat_ind) + preproc_dims%max_lat_ind = max_lat_ind + preproc_dims%min_lon = ecmwf%lon(min_lon_ind) + preproc_dims%min_lon_ind = min_lon_ind + preproc_dims%max_lon = ecmwf%lon(max_lon_ind) + preproc_dims%max_lon_ind = max_lon_ind + preproc_dims%dellat = (preproc_dims%ydim - 1) / (preproc_dims%max_lat - preproc_dims%min_lat) + preproc_dims%dellon = (preproc_dims%xdim - 1) / (preproc_dims%max_lon - preproc_dims%min_lon) + + if (verbose) then + print*, 'preproc_dims: ', preproc_dims%xdim, & + preproc_dims%ydim, preproc_dims%kdim + print*, 'dellon, dellat: ', preproc_dims%dellon, preproc_dims%dellat + print*, 'preproc_dims%min_lat: ', preproc_dims%min_lat + print*, 'preproc_dims%max_lat: ', preproc_dims%max_lat + print*, 'preproc_dims%min_lon: ', preproc_dims%min_lon + print*, 'preproc_dims%max_lon: ', preproc_dims%max_lon + end if + +end subroutine define_preproc_grid_ecmwf diff --git a/pre_processing/dependencies.inc b/pre_processing/dependencies.inc index 70dacb34..a66fa342 100644 --- a/pre_processing/dependencies.inc +++ b/pre_processing/dependencies.inc @@ -1,6 +1,6 @@ $(OBJS)/aatsr_corrections.o: $(OBJS)/calender.o $(OBJS)/preproc_constants.o $(OBJS)/bright_m.o: fundamental_constants.inc -$(OBJS)/build_preproc_fields.o: $(OBJS)/imager_structures.o \ +$(OBJS)/build_preproc_fields.o: $(OBJS)/ecmwf.o $(OBJS)/imager_structures.o \ $(OBJS)/preproc_constants.o $(OBJS)/preproc_structures.o $(OBJS)/calender.o: $(OBJS)/preproc_constants.o $(OBJS)/channel_structures.o: $(OBJS)/preproc_constants.o \ @@ -16,10 +16,12 @@ $(OBJS)/correct_for_ice_snow.o: $(OBJS)/channel_structures.o \ $(OBJS)/imager_structures.o $(OBJS)/nsidc_nise.o \ $(OBJS)/preproc_constants.o $(OBJS)/preproc_structures.o \ $(OBJS)/surface_structures.o +$(OBJS)/correct_for_dust.o: $(OBJS)/channel_structures.o \ + $(OBJS)/imager_structures.o $(OBJS)/cox_munk.o: $(OBJS)/cox_munk_constants.o $(OBJS)/gauss_leg_quad.o \ $(OBJS)/ocean_colour.o $(OBJS)/preproc_constants.o $(OBJS)/cox_munk_constants.o: $(OBJS)/preproc_constants.o -$(OBJS)/define_preprop_grid.o: $(OBJS)/imager_structures.o \ +$(OBJS)/define_preprop_grid.o: $(OBJS)/ecmwf.o $(OBJS)/imager_structures.o \ $(OBJS)/preproc_constants.o $(OBJS)/preproc_structures.o $(OBJS)/ecmwf.o: $(OBJS)/preproc_constants.o $(OBJS)/preproc_structures.o \ compute_geopot_coordinate.F90 deallocate_ecmwf_structures.F90 \ diff --git a/pre_processing/get_surface_emissivity.F90 b/pre_processing/get_surface_emissivity.F90 index 4156bd81..99a233c7 100644 --- a/pre_processing/get_surface_emissivity.F90 +++ b/pre_processing/get_surface_emissivity.F90 @@ -22,6 +22,7 @@ ! channel counts ! preproc_dims struct in Preprocessing dimensions, including sw and lw ! channel counts +! preproc_geoloc struct in Summary of preprocessing lat/lon ! assume_full_path logic in T: inputs are filenames; F: folder names ! verbose logic in T: print status information; F: don't ! surface struct both Surface properties structure @@ -70,14 +71,15 @@ ! 2015/01/13, AP: Alter channel indexing to allow channels in arbitrary order. ! 2015/10/19, GM: Turn back on reading of unused emissivity fields which are now ! optionally required. +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. !------------------------------------------------------------------------------- subroutine get_surface_emissivity(cyear, cdoy, cimss_emis_path, imager_flags, & - imager_geolocation, channel_info, preproc_dims, & - assume_full_path, verbose, surface, preproc_surf, source_atts) + imager_geolocation, channel_info, preproc_dims, preproc_geoloc, & + assume_full_path, verbose, surface, preproc_surf, source_atts) use channel_structures_m use cimss_emissivity_m @@ -98,6 +100,7 @@ subroutine get_surface_emissivity(cyear, cdoy, cimss_emis_path, imager_flags, & type(imager_geolocation_t), intent(in) :: imager_geolocation type(channel_info_t), intent(in) :: channel_info type(preproc_dims_t), intent(in) :: preproc_dims + type(preproc_geoloc_t), intent(in) :: preproc_geoloc logical, intent(in) :: assume_full_path logical, intent(in) :: verbose type(surface_t), intent(inout) :: surface @@ -175,10 +178,10 @@ subroutine get_surface_emissivity(cyear, cdoy, cimss_emis_path, imager_flags, & ! Read the data itself if (read_cimss_emissivity(cimss_emis_path_file, emis, & - channel_info%channel_wl_abs(ch_total_index), verbose) .ne. 0) then - write(*,*) 'ERROR: read_cimss_emissivity(), problem reading CIMSS ' // & - 'emissivity file: ', cimss_emis_path_file - stop error_stop_code + channel_info%channel_wl_abs(ch_total_index), verbose) .ne. 0) then + write(*,*) 'ERROR: read_cimss_emissivity(), problem reading CIMSS ' // & + 'emissivity file: ', cimss_emis_path_file + stop error_stop_code end if ! This emissivity data has very few missing values, but there are some. Set @@ -203,29 +206,26 @@ subroutine get_surface_emissivity(cyear, cdoy, cimss_emis_path, imager_flags, & imager_geolocation%latitude(i,j), interp) do k = 1, n_chans call interp_field(transemis(:,:,k), & - surface%emissivity(i,j,ch_lw_index(k)), interp) + surface%emissivity(i,j,ch_lw_index(k)), interp) end do end if end do end do ! calculate the mean emissivity in each preproc grid - allocate(counter(preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat)) - allocate(summat(preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat, n_chans)) + allocate(counter(1:preproc_dims%xdim, 1:preproc_dims%ydim)) + allocate(summat(1:preproc_dims%xdim, 1:preproc_dims%ydim, n_chans)) counter = 0 summat = 0. do j = 1, emis%nlat - lat = floor((emis%lat0+(j-1)*emis%lat_del+preproc_dims%lat_offset)* & - preproc_dims%dellat)+1 - if (lat.ge.preproc_dims%min_lat .and. lat.le.preproc_dims%max_lat) then + if ((emis%lat0+(j-1)*emis%lat_del).ge.minval(preproc_geoloc%latitude) .and. & + (emis%lat0+(j-1)*emis%lat_del).le.maxval(preproc_geoloc%latitude)) then + lat = minloc(abs(preproc_geoloc%latitude - (emis%lat0+(j-1)*emis%lat_del)),1) do i = 1, emis%nlon - lon = floor((emis%lon0+(i-1)*emis%lon_del+preproc_dims%lon_offset)* & - preproc_dims%dellon)+1 - if (lon.ge.preproc_dims%min_lon .and. & - lon.le.preproc_dims%max_lon) then + if ((emis%lon0+(i-1)*emis%lon_del).ge.minval(preproc_geoloc%longitude) .and. & + (emis%lon0+(i-1)*emis%lon_del).le.maxval(preproc_geoloc%longitude)) then + lon = minloc(abs(preproc_geoloc%longitude - (emis%lon0+(i-1)*emis%lon_del)),1) summat(lon,lat,:) = summat(lon,lat,:)+transemis(i,j,:) counter(lon,lat) = counter(lon,lat)+1 end if @@ -233,11 +233,11 @@ subroutine get_surface_emissivity(cyear, cdoy, cimss_emis_path, imager_flags, & end if end do - do j = preproc_dims%min_lat, preproc_dims%max_lat - do i = preproc_dims%min_lon, preproc_dims%max_lon + do j = 1, preproc_dims%ydim + do i = 1, preproc_dims%xdim if (counter(i,j) .gt. 0) then preproc_surf%emissivity(i,j,ch_lw_index) = summat(i,j,:) / & - real(counter(i,j)) + real(counter(i,j)) end if end do end do @@ -259,8 +259,8 @@ end subroutine get_surface_emissivity !------------------------------------------------------------------------------- subroutine get_camel_emissivity(cyear, cmonth, camel_emis_path, imager_flags, & - imager_geolocation, channel_info, preproc_dims, & - assume_full_path, verbose, surface, preproc_surf, source_atts) + imager_geolocation, channel_info, preproc_dims, preproc_geoloc, & + assume_full_path, verbose, surface, preproc_surf, source_atts) use channel_structures_m use camel_emissivity_m @@ -281,6 +281,7 @@ subroutine get_camel_emissivity(cyear, cmonth, camel_emis_path, imager_flags, & type(imager_geolocation_t), intent(in) :: imager_geolocation type(channel_info_t), intent(in) :: channel_info type(preproc_dims_t), intent(in) :: preproc_dims + type(preproc_geoloc_t), intent(in) :: preproc_geoloc logical, intent(in) :: assume_full_path logical, intent(in) :: verbose type(surface_t), intent(inout) :: surface @@ -358,10 +359,10 @@ subroutine get_camel_emissivity(cyear, cmonth, camel_emis_path, imager_flags, & ! Read the data itself if (read_camel_emissivity(camel_emis_path_file, emis, & - channel_info%channel_wl_abs(ch_total_index), verbose) .ne. 0) then - write(*,*) 'ERROR: read_camel_emissivity(), problem reading camel ' // & - 'emissivity file: ', camel_emis_path_file - stop error_stop_code + channel_info%channel_wl_abs(ch_total_index), verbose) .ne. 0) then + write(*,*) 'ERROR: read_camel_emissivity(), problem reading camel ' // & + 'emissivity file: ', camel_emis_path_file + stop error_stop_code end if ! This emissivity data has very few missing values, but there are some. Set @@ -382,29 +383,26 @@ subroutine get_camel_emissivity(cyear, cmonth, camel_emis_path, imager_flags, & imager_geolocation%latitude(i,j), interp) do k = 1, n_chans call interp_field(emis%emissivity(:,:,k), & - surface%emissivity(i,j,ch_lw_index(k)), interp) + surface%emissivity(i,j,ch_lw_index(k)), interp) end do end if end do end do ! calculate the mean emissivity in each preproc grid - allocate(counter(preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat)) - allocate(summat(preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat, n_chans)) + allocate(counter(1:preproc_dims%xdim, 1:preproc_dims%ydim)) + allocate(summat(1:preproc_dims%xdim, 1:preproc_dims%ydim, n_chans)) counter = 0 summat = 0. do j = 1, emis%nlat - lat = floor((emis%lat0+(j-1)*emis%lat_del+preproc_dims%lat_offset)* & - preproc_dims%dellat)+1 - if (lat.ge.preproc_dims%min_lat .and. lat.le.preproc_dims%max_lat) then + if ((emis%lat0+(j-1)*emis%lat_del).ge.minval(preproc_geoloc%latitude) .and. & + (emis%lat0+(j-1)*emis%lat_del).le.maxval(preproc_geoloc%latitude)) then + lat = minloc(abs(preproc_geoloc%latitude - (emis%lat0+(j-1)*emis%lat_del)),1) do i = 1, emis%nlon - lon = floor((emis%lon0+(i-1)*emis%lon_del+preproc_dims%lon_offset)* & - preproc_dims%dellon)+1 - if (lon.ge.preproc_dims%min_lon .and. & - lon.le.preproc_dims%max_lon) then + if ((emis%lon0+(i-1)*emis%lon_del).ge.minval(preproc_geoloc%longitude) .and. & + (emis%lon0+(i-1)*emis%lon_del).le.maxval(preproc_geoloc%longitude)) then + lon = minloc(abs(preproc_geoloc%longitude - (emis%lon0+(i-1)*emis%lon_del)),1) summat(lon,lat,:) = summat(lon,lat,:)+emis%emissivity(i,j,:) counter(lon,lat) = counter(lon,lat)+1 end if @@ -412,11 +410,11 @@ subroutine get_camel_emissivity(cyear, cmonth, camel_emis_path, imager_flags, & end if end do - do j = preproc_dims%min_lat, preproc_dims%max_lat - do i = preproc_dims%min_lon, preproc_dims%max_lon + do j = 1, preproc_dims%ydim + do i = 1, preproc_dims%xdim if (counter(i,j) .gt. 0) then preproc_surf%emissivity(i,j,ch_lw_index) = summat(i,j,:) / & - real(counter(i,j)) + real(counter(i,j)) end if end do end do diff --git a/pre_processing/netcdf_output_create_file.F90 b/pre_processing/netcdf_output_create_file.F90 index 73e4aa2b..c23f8670 100644 --- a/pre_processing/netcdf_output_create_file.F90 +++ b/pre_processing/netcdf_output_create_file.F90 @@ -99,6 +99,7 @@ ! 2018/04/29, SP: Add cloud emissivity support for ECMWF profiles (ExtWork) ! 2018/07/18, DE: Add tropoopause temperature ! 2018/11/05, SP: Add CAPE +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! None known. @@ -148,8 +149,8 @@ subroutine netcdf_create_rtm(global_atts, source_atts, cyear, cmonth, cday, chou integer(lint) :: nlon, nlat, kdim - nlon = preproc_dims%max_lon-preproc_dims%min_lon+1 - nlat = preproc_dims%max_lat-preproc_dims%min_lat+1 + nlon = preproc_dims%xdim + nlat = preproc_dims%ydim ! Set number of vertical levels/layers here, as GFS is different to ECMWF diff --git a/pre_processing/orac_preproc.F90 b/pre_processing/orac_preproc.F90 index e8b34505..eca4e1ba 100644 --- a/pre_processing/orac_preproc.F90 +++ b/pre_processing/orac_preproc.F90 @@ -450,9 +450,10 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi integer :: index_space real :: ecmwf_time_int_fac + integer :: date, ind -! Temporary variables for the aerosol_cci dust mask hack - real(kind=sreal), allocatable :: tot_cldmask_uncertainty(:,:) + ! Temporary variables for the aerosol_cci dust mask hack + !real(kind=sreal), allocatable :: tot_cldmask_uncertainty(:,:) ! integer, dimension(8) :: values @@ -505,6 +506,8 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi preproc_opts%use_seviri_ann_ctp_fg = .false. preproc_opts%use_seviri_ann_mlay = .false. preproc_opts%mcd43_max_qaflag = 5 + preproc_opts%do_dust_correction = .true. + preproc_opts%use_ecmwf_preproc_grid = .false. ! When true, the offset between the nadir and oblique views is read from ! the track_offset global attribute. Otherwise, the two longitude fields are @@ -877,16 +880,15 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi write(*,*) 'nwp_path_file3: ', trim(preproc_opts%nwp_fnames%nwp_path_file3(1)) end if end if - ! NOAA GFS has limited (pressure) levels and no HR, so set these. if (nwp_flag .eq. 0) preproc_opts%nwp_nlevels = 31 ! read surface wind fields and ECMWF dimensions if (preproc_opts%ecmwf_time_int_method .ne. 2) then - call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_opts%nwp_nlevels, verbose) + call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_opts%nwp_nlevels, date, ind, verbose) else - call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf1, preproc_opts%nwp_nlevels, verbose) - call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 2, ecmwf2, preproc_opts%nwp_nlevels, verbose) + call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf1, preproc_opts%nwp_nlevels, date, ind, verbose) + call read_ecmwf_wind(nwp_flag, preproc_opts%nwp_fnames, 2, ecmwf2, preproc_opts%nwp_nlevels, date, ind, verbose) call dup_ecmwf_allocation(ecmwf1, ecmwf) @@ -899,39 +901,52 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi ! define preprocessing grid from user grid spacing and satellite limits if (verbose) write(*,*) 'Define preprocessing grid' - preproc_dims%kdim = ecmwf%kdim - call define_preprop_grid(imager_geolocation, preproc_dims, verbose) + if (preproc_opts%use_ecmwf_preproc_grid) then + call define_preproc_grid_ecmwf(imager_geolocation, preproc_dims, ecmwf, verbose) + else + preproc_dims%kdim = ecmwf%kdim + call define_preprop_grid(imager_geolocation, preproc_dims, verbose) + + end if ! allocate preprocessing structures if (verbose) write(*,*) 'Allocate preprocessing structures' call allocate_preproc_structures(imager_angles, preproc_dims, & preproc_geoloc, preproc_geo, preproc_prtm, preproc_surf, preproc_cld, & channel_info) - ! now read the actual data and interpolate it to the preprocessing grid if (verbose) write(*,*) 'Build preprocessing grid' call build_preproc_fields(preproc_dims, preproc_geoloc, preproc_geo, & - imager_geolocation, imager_angles) + imager_geolocation, imager_angles, ecmwf, preproc_opts%use_ecmwf_preproc_grid) ! read ecmwf era interim file if (verbose) write(*,*) 'Read and interpolate NWP / Reanalysis data.' - if (preproc_opts%ecmwf_time_int_method .ne. 2) then - call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, preproc_geoloc, & - preproc_prtm, verbose) + + if (preproc_opts%use_ecmwf_preproc_grid.and.nwp_flag.gt.0.and.nwp_flag.lt.4) then + if (verbose) write(*,*) 'Using ECMWF as preproc grid' + call ecmwf_for_preproc_structures(preproc_opts, ecmwf, preproc_geoloc, & + preproc_prtm, preproc_dims, verbose, ecmwf_time_int_fac, date, ind, & + nwp_flag) else - call allocate_preproc_prtm(preproc_dims, preproc_prtm1) - call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, preproc_geoloc, & - preproc_prtm1, verbose) + if (preproc_opts%ecmwf_time_int_method .ne. 2) then + call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, & + preproc_geoloc, preproc_prtm, preproc_opts, date, ind, verbose) + else + call allocate_preproc_prtm(preproc_dims, preproc_prtm1) + call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, & + preproc_geoloc, preproc_prtm1, preproc_opts, date, ind, verbose) + + call allocate_preproc_prtm(preproc_dims, preproc_prtm2) + call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 2, ecmwf, preproc_dims, & + preproc_geoloc, preproc_prtm2, preproc_opts, date, ind, verbose) - call allocate_preproc_prtm(preproc_dims, preproc_prtm2) - call read_ecmwf(nwp_flag, preproc_opts%nwp_fnames, 2, ecmwf, preproc_dims, preproc_geoloc, & - preproc_prtm2, verbose) + call linearly_combine_prtms(1.-ecmwf_time_int_fac, ecmwf_time_int_fac, & + preproc_prtm1, preproc_prtm2, preproc_prtm) - call linearly_combine_prtms(1.-ecmwf_time_int_fac, ecmwf_time_int_fac, & - preproc_prtm1, preproc_prtm2, preproc_prtm) + call deallocate_preproc_prtm(preproc_prtm1) + call deallocate_preproc_prtm(preproc_prtm2) + end if - call deallocate_preproc_prtm(preproc_prtm1) - call deallocate_preproc_prtm(preproc_prtm2) end if if (verbose) write(*,*) 'Compute geopotential vertical coords' @@ -952,11 +967,11 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi if (verbose) write(*,*) 'Get surface emissivity' if (.not. preproc_opts%use_camel_emis) then call get_surface_emissivity(granule%cyear, granule%cdoy, cimss_emiss_path, & - imager_flags, imager_geolocation, channel_info, preproc_dims, & + imager_flags, imager_geolocation, channel_info, preproc_dims, preproc_geoloc, & assume_full_paths, verbose, surface, preproc_surf, source_atts) else call get_camel_emissivity(granule%cyear, granule%cmonth, cimss_emiss_path, & - imager_flags, imager_geolocation, channel_info, preproc_dims, & + imager_flags, imager_geolocation, channel_info, preproc_dims, preproc_geoloc, & assume_full_paths, verbose, surface, preproc_surf, source_atts) end if @@ -984,7 +999,7 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi else call correct_for_ice_snow_nwp(preproc_opts%nwp_fnames%nwp_path_file(1), & imager_geolocation, channel_info, imager_flags, preproc_dims, & - preproc_prtm, surface, include_full_brdf, source_atts, & + preproc_prtm, preproc_geoloc, surface, include_full_brdf, source_atts, & verbose) end if end if @@ -1012,58 +1027,10 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi end if end if - if (imager_angles%nviews .gt. 1) then - ! A temporary hack for Aerosol_cci: - ! Due to the cloud masking being very effective at detecting dust, - ! we'll try and re-introduce it - if (trim(adjustl(granule%sensor)) .eq. 'AATSR' .or. & - trim(adjustl(granule%sensor)) .eq. 'ATSR2' .or. & - trim(adjustl(granule%sensor)) .eq. 'SLSTR') then - if (1 .eq. 1 .and. & - minval(imager_geolocation%latitude) .lt. 40.0 .and. & - maxval(imager_geolocation%latitude) .gt. 0.0 .and. & - minval(imager_geolocation%longitude) .lt. 75.0 .and. & - maxval(imager_geolocation%longitude) .gt. -40.0) then - if (verbose) write(*,*) 'Aerosol_cci dust correction hack is underway' - allocate(tot_cldmask_uncertainty( & - imager_geolocation%startx:imager_geolocation%endx, & - 1:imager_geolocation%ny) ) - ! product a smoothed version of the cldmask uncertainty - if (verbose) then - write(*,*) minval(imager_pavolonis%cldmask_uncertainty(:,:,1)), & - maxval(imager_pavolonis%cldmask_uncertainty(:,:,1)) - write(*,*) minval(imager_pavolonis%cldmask_uncertainty(:,:,2)), & - maxval(imager_pavolonis%cldmask_uncertainty(:,:,2)) - end if - - tot_cldmask_uncertainty(:,:) = & - imager_pavolonis%cldmask_uncertainty(:,:,1) + & - imager_pavolonis%cldmask_uncertainty(:,:,2) - - if (verbose) write(*,*) 'Total cldmask uncertainty: min-max', minval(tot_cldmask_uncertainty), maxval(tot_cldmask_uncertainty) - ! Now use this smoothed mask, and the pavolonis cloud type - ! to "de-mask" possibly dust-filled pixels - ! Note that we leave the cldtype alone, so we can still tell - ! that the pixels were originally flagged as cloud - if (verbose) write(*,*) 'Total clouds before correction: ', & - count(imager_pavolonis%cldmask(:,:,1) .gt. 0), & - count(imager_pavolonis%cldmask(:,:,2) .gt. 0) - where(tot_cldmask_uncertainty .gt. 70 .and. & - (imager_pavolonis%cldtype(:,:,1) .eq. 3 .or. & - imager_pavolonis%cldtype(:,:,2) .eq. 3) .and. & - imager_geolocation%latitude .gt. 0.0 .and. & - imager_geolocation%latitude .lt. 40.0 .and. & - imager_geolocation%longitude .gt. -40.0 .and. & - imager_geolocation%longitude .lt. 75.0) - imager_pavolonis%cldmask(:,:,1) = 0 - imager_pavolonis%cldmask(:,:,2) = 0 - end where - if (verbose) write(*,*) 'Total clouds after correction: ', & - count(imager_pavolonis%cldmask(:,:,1) .gt. 0), & - count(imager_pavolonis%cldmask(:,:,2) .gt. 0) - deallocate(tot_cldmask_uncertainty) - end if - end if + if (preproc_opts%do_dust_correction) then + if (verbose) write(*,*) 'Apply dust-detection correction to cloud mask' + call correct_for_dust(channel_info, imager_measurements, imager_angles, & + imager_geolocation, imager_pavolonis, verbose) end if ! create output netcdf files. @@ -1088,7 +1055,7 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi if (preproc_opts%do_cloud_emis) then call get_cloud_emis(channel_info, imager_measurements, & imager_geolocation, preproc_dims, preproc_geoloc, & - preproc_cld, preproc_prtm, imager_cloud, ecmwf, & + preproc_cld, preproc_prtm, imager_cloud, & granule%sensor, verbose) end if #endif @@ -1105,7 +1072,7 @@ subroutine orac_preproc(mytask, ntasks, lower_bound, upper_bound, driver_path_fi #ifdef INCLUDE_SATWX call get_cloud_emis(channel_info, imager_measurements, & imager_geolocation, preproc_dims, preproc_geoloc, & - preproc_cld, preproc_prtm, imager_cloud, ecmwf, & + preproc_cld, preproc_prtm, imager_cloud, & granule%sensor, verbose) call do_cb_detect(channel_info, imager_measurements, & imager_geolocation, imager_cloud, imager_pavolonis, & diff --git a/pre_processing/preproc_structures.F90 b/pre_processing/preproc_structures.F90 index 00dc3671..1162789d 100644 --- a/pre_processing/preproc_structures.F90 +++ b/pre_processing/preproc_structures.F90 @@ -30,6 +30,7 @@ ! 2018/11/05, SP: Add CAPE ! 2021/03/09, AP: Consolidate path arguments into preproc_paths_t structure ! 2021/03/10, AP: Consolidate paths/dates into setup_args_t structure +! 2024/07/01, DH: Add option for using native ecmwf grid for preprocessing ! ! Bugs: ! None known. @@ -45,6 +46,7 @@ module preproc_structures_m integer(kind=lint) :: xdim, ydim, kdim integer(kind=lint) :: nchan_sw, nchan_lw integer(kind=lint) :: min_lat, max_lat, min_lon, max_lon + integer(kind=lint) :: min_lat_ind, max_lat_ind, min_lon_ind, max_lon_ind real(kind=sreal) :: dellon, dellat real(kind=sreal) :: lat_offset=90.0, lon_offset=180.0 @@ -105,6 +107,7 @@ module preproc_structures_m logical :: do_cloud_emis logical :: do_cloud_type logical :: do_ironly + logical :: do_dust_correction integer :: nwp_nlevels integer :: ecmwf_time_int_method integer, pointer :: channel_ids(:) @@ -123,6 +126,7 @@ module preproc_structures_m logical :: use_seviri_ann_cma_cph logical :: use_seviri_ann_ctp_fg logical :: use_seviri_ann_mlay + logical :: use_ecmwf_preproc_grid integer :: mcd43_max_qaflag character(len=path_length) :: ext_lsm_path diff --git a/pre_processing/read_ecmwf.F90 b/pre_processing/read_ecmwf.F90 index 1e9bad0b..2c6661d3 100644 --- a/pre_processing/read_ecmwf.F90 +++ b/pre_processing/read_ecmwf.F90 @@ -35,7 +35,7 @@ !------------------------------------------------------------------------------- subroutine read_ecmwf_wind(nwp_flag, nwp_fnames, idx, ecmwf, & - nwp_nlevels, verbose) + nwp_nlevels, date, ind, verbose) use preproc_structures_m @@ -47,6 +47,7 @@ subroutine read_ecmwf_wind(nwp_flag, nwp_fnames, idx, ecmwf, & type(ecmwf_t), intent(inout) :: ecmwf integer, intent(in) :: nwp_nlevels logical, intent(in) :: verbose + integer, intent(out) :: date, ind ! Set the number of levels in the input file, defaults to 61 select case(nwp_nlevels) @@ -93,9 +94,9 @@ subroutine read_ecmwf_wind(nwp_flag, nwp_fnames, idx, ecmwf, & ! It is possible for this field to have fill. It is set to -9e+33 for DWD ! and 9999.0 for Ox/RAL. Here we set it to ORAC's value. where (ecmwf%sea_ice_cover .lt. 0.0 .or. ecmwf%sea_ice_cover .gt. 1.0) & - ecmwf%sea_ice_cover = sreal_fill_value + ecmwf%sea_ice_cover = sreal_fill_value - call rearrange_ecmwf(ecmwf) + call rearrange_ecmwf(ecmwf, date, ind) end subroutine read_ecmwf_wind @@ -126,7 +127,7 @@ end subroutine read_ecmwf_wind !------------------------------------------------------------------------------- subroutine read_ecmwf(nwp_flag,nwp_fnames, idx, ecmwf, preproc_dims, & - preproc_geoloc, preproc_prtm, verbose) + preproc_geoloc, preproc_prtm, preproc_opts, date, ind, verbose) use preproc_structures_m @@ -134,12 +135,14 @@ subroutine read_ecmwf(nwp_flag,nwp_fnames, idx, ecmwf, preproc_dims, & integer, intent(in) :: nwp_flag type(preproc_nwp_fnames_t), intent(inout) :: nwp_fnames - integer, intent(in) :: idx - type(ecmwf_t), intent(in) :: ecmwf + integer, intent(in) :: idx + type(ecmwf_t), intent(inout) :: ecmwf type(preproc_dims_t), intent(in) :: preproc_dims type(preproc_geoloc_t), intent(in) :: preproc_geoloc type(preproc_prtm_t), intent(inout) :: preproc_prtm + type(preproc_opts_t), intent(inout) :: preproc_opts logical, intent(in) :: verbose + integer, intent(in) :: date, ind select case (nwp_flag) case(0) @@ -153,7 +156,7 @@ subroutine read_ecmwf(nwp_flag,nwp_fnames, idx, ecmwf, preproc_dims, & case(2) if (verbose) write(*,*) 'Reading JASMIN ERA5 path: ', trim(nwp_fnames%nwp_path_file(idx)) call read_era5_jasmin_nc(nwp_fnames, idx, ecmwf, preproc_dims, preproc_geoloc, & - preproc_prtm, verbose, nwp_flag) + preproc_prtm, preproc_opts, verbose, nwp_flag, date, ind) case(3) if (verbose) write(*,*) 'Reading ecmwf path: ', trim(nwp_fnames%nwp_path_file(idx)) call read_ecmwf_nc(nwp_fnames%nwp_path_file(idx), ecmwf, preproc_dims, preproc_geoloc, & @@ -172,5 +175,77 @@ subroutine read_ecmwf(nwp_flag,nwp_fnames, idx, ecmwf, preproc_dims, & preproc_prtm, verbose) end select - end subroutine read_ecmwf + +!------------------------------------------------------------------------------- +! Name: read_ecmwf +! +! Purpose: +! +! Description and Algorithm details: +! +! Arguments: +! Name Type In/Out/Both Description + +! Return value: +! Name Type Description +! ------------------------------------------------------------------------------ +! +! History: +! 2024/07/01, DH: Variation on the read_ecmwf subroutine when using the ecmwf +! grid as the preprocessing grid, as no interpolation is needed. +! +! Bugs: +! - only works for nwp_flag=1-3, not implemented for ERA-Interim and NOAA GFS +!------------------------------------------------------------------------------- + +subroutine ecmwf_for_preproc_structures(preproc_opts, ecmwf, preproc_geoloc, & + preproc_prtm, preproc_dims, verbose, ecmwf_time_int_fac, date, ind, nwp_flag) + use orac_ncdf_m + use preproc_constants_m + use preproc_structures_m + + implicit none + + + type(preproc_opts_t), intent(inout) :: preproc_opts + type(ecmwf_t), intent(inout) :: ecmwf + type(preproc_geoloc_t), intent(inout) :: preproc_geoloc + type(preproc_prtm_t), intent(inout) :: preproc_prtm + type(preproc_dims_t), intent(inout) :: preproc_dims + real, intent(in) :: ecmwf_time_int_fac + integer, intent(in) :: date, ind + logical, intent(in) :: verbose + integer, intent(in) :: nwp_flag + + type(preproc_prtm_t) :: preproc_prtm1 + type(preproc_prtm_t) :: preproc_prtm2 + + select case (nwp_flag) + case(1) + if (verbose) write(*,*) 'Reading ECMWF path: ', trim(preproc_opts%nwp_fnames%nwp_path_file(1)) + call ecmwf_nc_for_preproc_structures(preproc_opts, ecmwf, preproc_geoloc, & + preproc_prtm, preproc_dims, verbose, ecmwf_time_int_fac, date, ind) + case(2) + if (verbose) write(*,*) 'Reading JASMIN ERA5' + if (preproc_opts%ecmwf_time_int_method .ne. 2) then + call read_era5_jasmin_nc(preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm, preproc_opts, verbose, nwp_flag, date, ind) + else + call allocate_preproc_prtm(preproc_dims, preproc_prtm1) + call read_era5_jasmin_nc(preproc_opts%nwp_fnames, 1, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm1, preproc_opts, verbose, nwp_flag, date, ind) + call allocate_preproc_prtm(preproc_dims, preproc_prtm2) + call read_era5_jasmin_nc(preproc_opts%nwp_fnames, 2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm2, preproc_opts, verbose, nwp_flag, date, ind) + call linearly_combine_prtms(1.-ecmwf_time_int_fac, ecmwf_time_int_fac, & + preproc_prtm1, preproc_prtm2, preproc_prtm) + + call deallocate_preproc_prtm(preproc_prtm1) + call deallocate_preproc_prtm(preproc_prtm2) + + end if + case(3) + if (verbose) write(*,*) 'Reading ecmwf path: ', trim(preproc_opts%nwp_fnames%nwp_path_file(1)) + call ecmwf_nc_for_preproc_structures(preproc_opts, ecmwf, preproc_geoloc, & + preproc_prtm, preproc_dims, verbose, ecmwf_time_int_fac, date, ind) + end select + +end subroutine ecmwf_for_preproc_structures diff --git a/pre_processing/read_ecmwf_grib.F90 b/pre_processing/read_ecmwf_grib.F90 index 29ec6a48..3e3a486c 100644 --- a/pre_processing/read_ecmwf_grib.F90 +++ b/pre_processing/read_ecmwf_grib.F90 @@ -60,6 +60,7 @@ ! 2016/12/07, GT: Replaced call to INTF with INTF2, as INTF requires the GRIBEX ! subroutine from EMOSLIB when interpolating from or to a GRIB field, which ! is no-longer supported (or it appears functional) from v4.1.1 of the EMOSLIB +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! - If you're having problems with INTF, set the environment variable JDCNDBG=1 @@ -81,8 +82,10 @@ subroutine read_ecmwf_grib(ecmwf_file, preproc_dims, preproc_geoloc, & type(preproc_prtm_t), intent(inout) :: preproc_prtm logical, intent(in) :: verbose +#ifdef INCLUDE_EMOS integer(lint), parameter :: BUFFER = 3000000 integer(lint), external :: INTIN, INTOUT, INTF2 +#endif integer(lint) :: fu, stat, nbytes ! integer(lint) :: in_words, out_words integer(lint) :: out_bytes, out_words @@ -97,6 +100,7 @@ subroutine read_ecmwf_grib(ecmwf_file, preproc_dims, preproc_geoloc, & real(sreal), dimension(:), allocatable :: pl, val real(sreal), dimension(:,:), pointer :: array +#ifdef INCLUDE_EMOS ! open the ECMWF file call PBOPEN(fu, ecmwf_file, 'r', stat) if (stat .ne. 0) call h_e_e('grib', 'Failed to read file.') @@ -113,10 +117,10 @@ subroutine read_ecmwf_grib(ecmwf_file, preproc_dims, preproc_geoloc, & grid(2) = 0.5 / preproc_dims%dellat if (INTOUT('grid', iblank, grid, charv) .ne. 0) & call h_e_e('grib', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) if (INTOUT('area', iblank, area, charv) .ne. 0) & call h_e_e('grib', 'INTOUT area failed.') @@ -178,64 +182,64 @@ subroutine read_ecmwf_grib(ecmwf_file, preproc_dims, preproc_geoloc, & select case (param) case(130) array => preproc_prtm%temperature( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,level) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,level) case(133) array => preproc_prtm%spec_hum( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,level) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,level) case(203) array => preproc_prtm%ozone( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,level) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,level) case(129) array => preproc_prtm%geopot( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(152) array => preproc_prtm%lnsp( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(31) - array => preproc_prtm%sea_ice_cover( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%sea_ice_cover( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(32) - array => preproc_prtm%snow_albedo( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%snow_albedo( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(34) - array => preproc_prtm%sst( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%sst( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(137) - array => preproc_prtm%totcolwv( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%totcolwv( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(141) - array => preproc_prtm%snow_depth( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%snow_depth( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(165) - array => preproc_prtm%u10( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%u10( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(166) - array => preproc_prtm%v10( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%v10( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(167) - array => preproc_prtm%temp2( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%temp2( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(172) - array => preproc_prtm%land_sea_mask( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%land_sea_mask( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(235) - array => preproc_prtm%skin_temp( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + array => preproc_prtm%skin_temp( & + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case default cycle end select @@ -263,5 +267,9 @@ subroutine read_ecmwf_grib(ecmwf_file, preproc_dims, preproc_geoloc, & ! close ECMWF file call PBCLOSE(fu, stat) if (stat .ne. 0) call h_e_e('grib', 'Failed to close file.') +#else + write(*,*) 'ERROR: read_ecmwf_grib(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine read_ecmwf_grib diff --git a/pre_processing/read_ecmwf_nc.F90 b/pre_processing/read_ecmwf_nc.F90 index ca0e8aa8..2e1376c2 100644 --- a/pre_processing/read_ecmwf_nc.F90 +++ b/pre_processing/read_ecmwf_nc.F90 @@ -17,11 +17,11 @@ ! Arguments: ! Name Type In/Out/Both Description ! ------------------------------------------------------------------------------ -! nwp_path string In NetCDF ECMWF file to be opened. -! ecmwf struct both Structure summarising contents of ECMWF files. +! nwp_path string In NetCDF ECMWF file to be opened. +! ecmwf struct Both Structure summarising contents of ECMWF files. ! preproc_dims struct In Dimensions of the preprocessing grid. ! preproc_prtm struct Both Pressure-level information for RTTOV. -! verbose logic in T: Print min/max of each field; F: Don't. +! verbose logic In T: Print min/max of each field; F: Don't. ! ! History: ! 2012/08/06, ??: Initial version ecmwf code @@ -92,8 +92,10 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & logical, intent(in) :: verbose integer, intent(in) :: nwp_flag +#ifdef INCLUDE_EMOS integer(lint), external :: INTIN, INTOUT, INTF integer(lint), parameter :: BUFFER = 4000000 +#endif integer(lint), dimension(1) :: intv, old_grib, new_grib real(dreal) :: grid(2), area(4) @@ -101,7 +103,7 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & character(len=20), dimension(1) :: charv real(sreal), pointer :: array2d(:,:), array3d(:,:,:) - integer(4) :: n, ni, nj, i, j, k, ivar + integer(lint) :: n, ni, nj, i, j, k, ivar integer(4) :: fid, nvar integer(4) :: old_len, new_len character(len=20) :: name @@ -112,13 +114,13 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & #ifdef WRAPPER real(sreal) :: ecmwf_lon(ecmwf%xdim) real(sreal) :: ecmwf_lat(ecmwf%ydim) - integer(lint) :: pointer_x(preproc_dims%min_lon:preproc_dims%max_lon) - integer(lint) :: pointer_y(preproc_dims%min_lat:preproc_dims%max_lat) + integer(lint) :: pointer_x(preproc_dims%xdim) + integer(lint) :: pointer_y(preproc_dims%ydim) real(sreal) :: diff_lon(ecmwf%xdim), diff_lat(ecmwf%ydim) #endif - +#ifdef INCLUDE_EMOS n = ecmwf%xdim*ecmwf%ydim ! input details of new grid (see note in read_ecmwf_grib) @@ -138,10 +140,10 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & grid(2) = 0.5 / preproc_dims%dellat if (INTOUT('grid', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) if (INTOUT('area', intv, area, charv) .ne. 0) & call h_e_e('nc', 'INTOUT area failed.') ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 @@ -167,12 +169,12 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & where(ecmwf_lon .gt. 180.) ecmwf_lon = ecmwf_lon-360. - do i = preproc_dims%min_lon, preproc_dims%max_lon + do i = 1, preproc_dims%xdim diff_lon = abs(ecmwf_lon - preproc_geoloc%longitude(i)) pointer_x(i) = minloc(diff_lon, 1) end do - do j = preproc_dims%min_lat, preproc_dims%max_lat + do j = 1, preproc_dims%ydim diff_lat = abs(ecmwf_lat - preproc_geoloc%latitude(j)) pointer_y(j) = minloc(diff_lat, 1) end do @@ -259,19 +261,17 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & if (INTF(old_grib, old_len, old_data, new_grib, new_len, new_data).ne.0)& call h_e_e('nc', 'INTF failed.') if (new_len .ne. ni*nj) print*, '3D Interpolation grid wrong.' - ! copy data into preprocessing grid - do j = 1, nj, 2 - do i = 1, ni, 2 - array3d(preproc_dims%min_lon+i/2, & - preproc_dims%min_lat+(nj-j)/2, k) = & + do j = 1, nj,2 + do i = 1, ni,2 + array3d(1+i/2,1+(nj-j)/2, k) = & real(new_data(i+(j-1)*ni), kind=4) end do end do #else ! copy data into preprocessing grid - do i = preproc_dims%min_lon, preproc_dims%max_lon - do j = preproc_dims%min_lat, preproc_dims%max_lat + do i = 1, preproc_dims%xdim + do j = 1, preproc_dims%ydim array3d(i,j,k) = real(dummy3d(pointer_x(i), pointer_y(j), k, 1)) end do end do @@ -294,15 +294,14 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & ! copy data into preprocessing grid do j = 1, nj, 2 do i = 1, ni, 2 - array2d(preproc_dims%min_lon+i/2, & - preproc_dims%min_lat+(nj-j)/2) = & + array2d(1+i/2,1+(nj-j)/2) = & real(new_data(i+(j-1)*ni), kind=4) end do end do #else ! copy data into preprocessing grid - do i = preproc_dims%min_lon, preproc_dims%max_lon - do j = preproc_dims%min_lat, preproc_dims%max_lat + do i = 1, preproc_dims%xdim + do j = 1, preproc_dims%ydim array2d(i,j) = real(dummy2d(pointer_x(i), pointer_y(j), 1, 1)) end do end do @@ -316,5 +315,240 @@ subroutine read_ecmwf_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & deallocate(new_data) call ncdf_close(fid, 'read_ecmwf_nc()') +#else + write(*,*) 'ERROR: read_ecmwf_nc(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine read_ecmwf_nc + +!------------------------------------------------------------------------------- +! Name: read_ecmwf_nc.F90 +! +! Purpose: +! Reads NetCDF format ECMWF ERA interim data. It is interpolated onto the +! preprocessing grid using the EMOS package +! +! Description and Algorithm details: +! 1) Open file. +! 2) Loop over variables: +! a) Identify variable with desired data field. +! b) Read and copy data into preprocessor structure. +! 3) Close file. +! +! Arguments: +! Name Type In/Out/Both Description +! ------------------------------------------------------------------------------ +! nwp_path string In NetCDF ECMWF file to be opened. +! ecmwf struct Both Structure summarising contents of ECMWF files. +! preproc_dims struct In Dimensions of the preprocessing grid. +! preproc_prtm struct Both Pressure-level information for RTTOV. +! verbose logic In T: Print min/max of each field; F: Don't. +! +! History: +! 2024/07/01, DH: Initial version ecmwf code when using ecmwf grid for +! preprocessing +! +! Bugs: +! - you need to be careful with parameter naming as the variable names are not +! consistent across files for example the variable name could be lnsp or LNSP +!------------------------------------------------------------------------------- + + +subroutine ecmwf_nc_for_preproc_structures(preproc_opts, ecmwf, preproc_geoloc, preproc_prtm, preproc_dims, verbose, ecmwf_time_int_fac, date, ind) + use orac_ncdf_m + use preproc_constants_m + use preproc_structures_m + !use ecmwf_m + + implicit none + + type(preproc_opts_t), intent(inout) :: preproc_opts + type(ecmwf_t), intent(inout) :: ecmwf + type(preproc_geoloc_t), intent(inout) :: preproc_geoloc + type(preproc_prtm_t), intent(inout) :: preproc_prtm + type(preproc_dims_t), intent(inout) :: preproc_dims + real, intent(in) :: ecmwf_time_int_fac + type(preproc_prtm_t) :: preproc_prtm1 + type(preproc_prtm_t) :: preproc_prtm2 + integer, intent(in) :: date, ind + logical, intent(in) :: verbose + + character(len=20) :: name + integer(4) :: ivar + integer(4) :: fid,fid0, fid1,fid2, nvar + + + ! open file + call ncdf_open(fid, preproc_opts%nwp_fnames%nwp_path_file(1), 'read_ecmwf_nc()') + if (nf90_inquire(fid, nVariables=nvar) .ne. 0) & + call h_e_e('nc', 'NF INQ failed.') + + if (preproc_opts%ecmwf_time_int_method .eq. 2) then + call allocate_preproc_prtm(preproc_dims, preproc_prtm1) + call allocate_preproc_prtm(preproc_dims, preproc_prtm2) + end if + + ! loop over variables + do ivar = 1, nvar + + if(ivar .ge. 2) then + if (preproc_opts%ecmwf_time_int_method .ne. 2) then + call ecmwf_read_var(fid0, preproc_opts%nwp_fnames, 1, preproc_prtm, preproc_dims, ecmwf, ivar, date, ind, name) + else + call ecmwf_read_var(fid0, preproc_opts%nwp_fnames, 1, preproc_prtm, preproc_dims, ecmwf, ivar, date, ind, name) + + call ecmwf_read_var(fid1, preproc_opts%nwp_fnames, 1, preproc_prtm1, preproc_dims, ecmwf, ivar, date, ind, name) + + call ecmwf_read_var(fid2, preproc_opts%nwp_fnames, 2, preproc_prtm2, preproc_dims, ecmwf, ivar, date, ind, name) + + end if + end if + end do + if (preproc_opts%ecmwf_time_int_method .eq. 2) then + call linearly_combine_prtms(1.-ecmwf_time_int_fac, ecmwf_time_int_fac, & + preproc_prtm1, preproc_prtm2, preproc_prtm) + end if + call deallocate_preproc_prtm(preproc_prtm1) + call deallocate_preproc_prtm(preproc_prtm2) + call ncdf_close(fid, 'read_ecmwf_nc()') +end subroutine ecmwf_nc_for_preproc_structures + + +subroutine ecmwf_read_var(fid, nwp_fnames, idx, preproc_prtm, preproc_dims, ecmwf, ivar, date, ind, name) + use orac_ncdf_m + use preproc_constants_m + use preproc_structures_m + !use ecmwf_m + + implicit none + + integer, intent(in) :: idx + type(preproc_prtm_t), intent(inout) :: preproc_prtm + type(preproc_dims_t), intent(inout) :: preproc_dims + type(ecmwf_t), intent(inout) :: ecmwf + character(len=20), intent(out) :: name + type(preproc_nwp_fnames_t), intent(inout) :: nwp_fnames + integer, intent(in) :: date, ind + integer(4), intent(in) :: ivar + integer(4) :: fid + integer(4) :: nvar + real(sreal) :: dummy2d(ecmwf%xdim,ecmwf%ydim) + real(sreal) :: dummy3d(ecmwf%xdim,ecmwf%ydim,ecmwf%kdim) +#ifdef WRAPPER + real(sreal) :: ecmwf_lon(ecmwf%xdim) + real(sreal) :: ecmwf_lat(ecmwf%ydim) +#endif + + ! open file + call ncdf_open(fid, nwp_fnames%nwp_path_file(idx), 'read_ecmwf_nc()') + if (nf90_inquire(fid, nVariables=nvar) .ne. 0) & + call h_e_e('nc', 'NF INQ failed.') + if (nf90_inquire_variable(fid, ivar, name) .ne. 0) & + call h_e_e('nc', 'NF VAR INQUIRE failed.') + + ! determine if field exists + select case (name) + case('Z', 'z') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%geopot = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('Q', 'q') + call ncdf_read_array(fid, name, dummy3d) + call rearrange_ecmwf_var3d(ecmwf, dummy3d, date, ind) + preproc_prtm%spec_hum = dummy3d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + case('T') + call ncdf_read_array(fid, name, dummy3d) + call rearrange_ecmwf_var3d(ecmwf, dummy3d, date, ind) + preproc_prtm%temperature = dummy3d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + case('O3', 'o3') + call ncdf_read_array(fid, name, dummy3d) + call rearrange_ecmwf_var3d(ecmwf, dummy3d, date, ind) + preproc_prtm%ozone = dummy3d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + case('LNSP', 'lnsp') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%lnsp = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('CI', 'ci') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%sea_ice_cover = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('ASN', 'asn') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%snow_albedo = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('TCWV', 'tcwv') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%totcolwv = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('SD', 'sd') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%snow_depth = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('U10', 'u10', 'U10M') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%u10 = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('V10', 'v10', 'V10M') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%v10 = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('T2', 't2', 'T2M') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%temp2 = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('SKT', 'skt') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%skin_temp = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('SSTK', 'sstk') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%sst = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + case('AL', 'al', 'LSM') + call ncdf_read_array(fid, name, dummy2d) + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%land_sea_mask = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) +#ifdef INCLUDE_SATWX + case('cape', 'CAPE') + call ncdf_read_array(fid, name, dummy3d) + call rearrange_ecmwf_var3d(ecmwf, dummy3d, date, ind) + preproc_prtm%cape = dummy2d( & + preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) +#endif + end select + + call ncdf_close(fid, 'read_ecmwf_nc()') + +end subroutine ecmwf_read_var diff --git a/pre_processing/read_era5_jasmin.F90 b/pre_processing/read_era5_jasmin.F90 index 2e2bf517..ed247843 100644 --- a/pre_processing/read_era5_jasmin.F90 +++ b/pre_processing/read_era5_jasmin.F90 @@ -34,288 +34,392 @@ end subroutine read_era5_jasmin_wind_nc subroutine read_era5_jasmin_nc(nwp_fnames, idx, ecmwf, preproc_dims, preproc_geoloc, & - preproc_prtm, verbose, nwp_flag) - - use orac_ncdf_m - use preproc_constants_m - use preproc_structures_m - - implicit none - - type(preproc_nwp_fnames_t), intent(in) :: nwp_fnames - integer, intent(in) :: idx - type(ecmwf_t), intent(in) :: ecmwf - type(preproc_dims_t), intent(in) :: preproc_dims - type(preproc_geoloc_t), intent(in) :: preproc_geoloc - type(preproc_prtm_t), intent(inout) :: preproc_prtm - logical, intent(in) :: verbose - integer, intent(in) :: nwp_flag - - integer(lint), external :: INTIN, INTOUT, INTF - integer(lint), parameter :: BUFFER = 2000000 - - real(sreal) :: dummy2d(ecmwf%xdim,ecmwf%ydim) - real(sreal) :: dummy3d_2(ecmwf%xdim,ecmwf%ydim,ecmwf%kdim) - - ! open file - ! Do temperature on model levels - call load_era5_netcdf_3d(nwp_fnames%t_f(idx), 't', dummy3d_2) - call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%temperature) - if (verbose) print*, 'T) Min: ', minval(preproc_prtm%temperature), ', Max: ', maxval(preproc_prtm%temperature) - - ! Do specific humidity on model levels - call load_era5_netcdf_3d(nwp_fnames%q_f(idx), 'q', dummy3d_2) - call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%spec_hum) - if (verbose) print*, 'Q) Min: ', minval(preproc_prtm%spec_hum), ', Max: ', maxval(preproc_prtm%spec_hum) - - ! Do ozone on model levels - call load_era5_netcdf_3d(nwp_fnames%o3_f(idx), 'o3', dummy3d_2) - call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%ozone) - if (verbose) print*, 'O3) Min: ', minval(preproc_prtm%ozone), ', Max: ', maxval(preproc_prtm%ozone) - - ! Do logarithm of surface pressure - call load_era5_netcdf_2d(nwp_fnames%lnsp_f(idx), 'lnsp', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%lnsp) - if (verbose) print*, 'LNSP) Min: ', minval(preproc_prtm%lnsp), ', Max: ', maxval(preproc_prtm%lnsp) - - ! Do geopotential - call load_era5_netcdf_2d(nwp_fnames%z_f(idx), 'z', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%geopot) - if (verbose) print*, 'GEOP) Min: ', minval(preproc_prtm%geopot), ', Max: ', maxval(preproc_prtm%geopot) - - ! Do sea ice fraction - call load_era5_netcdf_2d(nwp_fnames%ci_f(idx), 'siconc', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%sea_ice_cover) - if (verbose) print*, 'CI) Min: ', minval(preproc_prtm%sea_ice_cover), ', Max: ', maxval(preproc_prtm%sea_ice_cover) - - ! Do snow albedo - call load_era5_netcdf_2d(nwp_fnames%asn_f(idx), 'asn', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%snow_albedo) - if (verbose) print*, 'ASN) Min: ', minval(preproc_prtm%snow_albedo), ', Max: ', maxval(preproc_prtm%snow_albedo) - - ! Do total column water vapour - call load_era5_netcdf_2d(nwp_fnames%tcwv_f(idx), 'tcwv', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%totcolwv) - if (verbose) print*, 'TCWV) Min: ', minval(preproc_prtm%totcolwv), ', Max: ', maxval(preproc_prtm%totcolwv) - - ! Do snow depth - call load_era5_netcdf_2d(nwp_fnames%sd_f(idx), 'sd', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%snow_depth) - if (verbose) print*, 'SD) Min: ', minval(preproc_prtm%snow_depth), ', Max: ', maxval(preproc_prtm%snow_depth) - - ! Do U-component of 10m wind - call load_era5_netcdf_2d(nwp_fnames%u10_f(idx), 'u10', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%u10) - if (verbose) print*, 'U10) Min: ', minval(preproc_prtm%u10), ', Max: ', maxval(preproc_prtm%u10) - - ! Do V-component of 10m wind - call load_era5_netcdf_2d(nwp_fnames%v10_f(idx), 'v10', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%v10) - if (verbose) print*, 'V10) Min: ', minval(preproc_prtm%v10), ', Max: ', maxval(preproc_prtm%v10) - - ! Do 2m temperature - call load_era5_netcdf_2d(nwp_fnames%t2_f(idx), 't2m', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%temp2) - if (verbose) print*, 'T2) Min: ', minval(preproc_prtm%temp2), ', Max: ', maxval(preproc_prtm%temp2) - - ! Do skin temperature - call load_era5_netcdf_2d(nwp_fnames%skt_f(idx), 'skt', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%skin_temp) - if (verbose) print*, 'SKT) Min: ', minval(preproc_prtm%skin_temp), ', Max: ', maxval(preproc_prtm%skin_temp) - - ! Do sea surface temperature - call load_era5_netcdf_2d(nwp_fnames%sstk_f(idx), 'sst', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%sst) - if (verbose) print*, 'SST) Min: ', minval(preproc_prtm%sst), ', Max: ', maxval(preproc_prtm%sst) + preproc_prtm, preproc_opts, verbose, nwp_flag, date, ind) + + use orac_ncdf_m + use preproc_constants_m + use preproc_structures_m + + implicit none + + type(preproc_opts_t), intent(inout) :: preproc_opts + type(preproc_nwp_fnames_t), intent(in) :: nwp_fnames + integer, intent(in) :: idx + type(ecmwf_t), intent(inout) :: ecmwf + type(preproc_dims_t), intent(in) :: preproc_dims + type(preproc_geoloc_t), intent(in) :: preproc_geoloc + type(preproc_prtm_t), intent(inout) :: preproc_prtm + logical, intent(in) :: verbose + integer, intent(in) :: nwp_flag + integer, intent(in) :: date, ind + + real(sreal) :: dummy2d(ecmwf%xdim,ecmwf%ydim) + real(sreal) :: dummy3d_2(ecmwf%xdim,ecmwf%ydim,ecmwf%kdim) + + ! open file + ! Do temperature on model levels + call load_era5_netcdf_3d(nwp_fnames%t_f(idx), 't', dummy3d_2) + if (verbose) print*, trim(adjustl(nwp_fnames%t_f(idx))) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%temperature) + else + call rearrange_ecmwf_var3d(ecmwf, dummy3d_2, date, ind) + preproc_prtm%temperature = dummy3d_2(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + end if + if (verbose) print*, 'T) Min: ', minval(preproc_prtm%temperature), ', Max: ', maxval(preproc_prtm%temperature) + + ! Do specific humidity on model levels + call load_era5_netcdf_3d(nwp_fnames%q_f(idx), 'q', dummy3d_2) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%spec_hum) + else + call rearrange_ecmwf_var3d(ecmwf, dummy3d_2, date, ind) + preproc_prtm%spec_hum = dummy3d_2(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + end if + if (verbose) print*, 'Q) Min: ', minval(preproc_prtm%spec_hum), ', Max: ', maxval(preproc_prtm%spec_hum) + + ! Do ozone on model levels + call load_era5_netcdf_3d(nwp_fnames%o3_f(idx), 'o3', dummy3d_2) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_3d_var(dummy3d_2, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%ozone) + else + call rearrange_ecmwf_var3d(ecmwf, dummy3d_2, date, ind) + preproc_prtm%ozone = dummy3d_2(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind, :) + end if + if (verbose) print*, 'O3) Min: ', minval(preproc_prtm%ozone), ', Max: ', maxval(preproc_prtm%ozone) + + ! Do logarithm of surface pressure + call load_era5_netcdf_2d(nwp_fnames%lnsp_f(idx), 'lnsp', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%lnsp) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%lnsp = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'LNSP) Min: ', minval(preproc_prtm%lnsp), ', Max: ', maxval(preproc_prtm%lnsp) + + ! Do geopotential + call load_era5_netcdf_2d(nwp_fnames%z_f(idx), 'z', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%geopot) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%geopot = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'GEOP) Min: ', minval(preproc_prtm%geopot), ', Max: ', maxval(preproc_prtm%geopot) + + ! Do sea ice fraction + call load_era5_netcdf_2d(nwp_fnames%ci_f(idx), 'siconc', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%sea_ice_cover) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%sea_ice_cover = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'CI) Min: ', minval(preproc_prtm%sea_ice_cover), ', Max: ', maxval(preproc_prtm%sea_ice_cover) + + ! Do snow albedo + call load_era5_netcdf_2d(nwp_fnames%asn_f(idx), 'asn', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%snow_albedo) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%snow_albedo = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'ASN) Min: ', minval(preproc_prtm%snow_albedo), ', Max: ', maxval(preproc_prtm%snow_albedo) + + ! Do total column water vapour + call load_era5_netcdf_2d(nwp_fnames%tcwv_f(idx), 'tcwv', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%totcolwv) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%totcolwv = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'TCWV) Min: ', minval(preproc_prtm%totcolwv), ', Max: ', maxval(preproc_prtm%totcolwv) + + ! Do snow depth + call load_era5_netcdf_2d(nwp_fnames%sd_f(idx), 'sd', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%snow_depth) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%snow_depth = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'SD) Min: ', minval(preproc_prtm%snow_depth), ', Max: ', maxval(preproc_prtm%snow_depth) + + ! Do U-component of 10m wind + call load_era5_netcdf_2d(nwp_fnames%u10_f(idx), 'u10', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%u10) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%u10 = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'U10) Min: ', minval(preproc_prtm%u10), ', Max: ', maxval(preproc_prtm%u10) + + ! Do V-component of 10m wind + call load_era5_netcdf_2d(nwp_fnames%v10_f(idx), 'v10', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%v10) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%v10 = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'V10) Min: ', minval(preproc_prtm%v10), ', Max: ', maxval(preproc_prtm%v10) + + ! Do 2m temperature + call load_era5_netcdf_2d(nwp_fnames%t2_f(idx), 't2m', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%temp2) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%temp2 = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'T2) Min: ', minval(preproc_prtm%temp2), ', Max: ', maxval(preproc_prtm%temp2) + + ! Do skin temperature + call load_era5_netcdf_2d(nwp_fnames%skt_f(idx), 'skt', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%skin_temp) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%skin_temp = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'SKT) Min: ', minval(preproc_prtm%skin_temp), ', Max: ', maxval(preproc_prtm%skin_temp) + + ! Do sea surface temperature + call load_era5_netcdf_2d(nwp_fnames%sstk_f(idx), 'sst', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%sst) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%sst = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'SST) Min: ', minval(preproc_prtm%sst), ', Max: ', maxval(preproc_prtm%sst) #ifdef INCLUDE_SATWX - ! Do convective available potential energy - call load_era5_netcdf_2d(nwp_fnames%cape_f(idx), 'cape', dummy2d) - call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%cape) - if (verbose) print*, 'CAPE) Min: ', minval(preproc_prtm%cape), ', Max: ', maxval(preproc_prtm%cape) + ! Do convective available potential energy + call load_era5_netcdf_2d(nwp_fnames%cape_f(idx), 'cape', dummy2d) + if (.not. preproc_opts%use_ecmwf_preproc_grid) then + call preproc_2d_var(dummy2d, ecmwf, preproc_dims, preproc_geoloc, preproc_prtm%cape) + else + call rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + preproc_prtm%cape = dummy2d(preproc_dims%min_lon_ind:preproc_dims%max_lon_ind, & + preproc_dims%min_lat_ind:preproc_dims%max_lat_ind) + end if + if (verbose) print*, 'CAPE) Min: ', minval(preproc_prtm%cape), ', Max: ', maxval(preproc_prtm%cape) #endif end subroutine read_era5_jasmin_nc subroutine load_era5_netcdf_3d(fname, varname, data_arr) - use orac_ncdf_m - implicit none + use orac_ncdf_m + implicit none - character(len=*), intent(in) :: fname - character(len=*), intent(in) :: varname - real(sreal), dimension(:,:,:), intent(inout) :: data_arr - integer(4) :: fid + character(len=*), intent(in) :: fname + character(len=*), intent(in) :: varname + real(sreal), dimension(:,:,:), intent(inout) :: data_arr + integer(4) :: fid - call ncdf_open(fid, trim(adjustl(fname)), 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') - call ncdf_read_array(fid, trim(adjustl(varname)), data_arr) - call ncdf_close(fid, 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') + call ncdf_open(fid, trim(adjustl(fname)), 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') + call ncdf_read_array(fid, trim(adjustl(varname)), data_arr) + call ncdf_close(fid, 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') end subroutine load_era5_netcdf_3d subroutine load_era5_netcdf_2d(fname, varname, data_arr) - use orac_ncdf_m - implicit none + use orac_ncdf_m + implicit none - character(len=*), intent(in) :: fname - character(len=*), intent(in) :: varname - real(sreal), dimension(:,:), intent(inout) :: data_arr - integer(4) :: fid + character(len=*), intent(in) :: fname + character(len=*), intent(in) :: varname + real(sreal), dimension(:,:), intent(inout) :: data_arr + integer(4) :: fid - call ncdf_open(fid, trim(adjustl(fname)), 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') - call ncdf_read_array(fid, trim(adjustl(varname)), data_arr) - call ncdf_close(fid, 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') + call ncdf_open(fid, trim(adjustl(fname)), 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') + call ncdf_read_array(fid, trim(adjustl(varname)), data_arr) + call ncdf_close(fid, 'read_era5_jasmin_nc_'//trim(adjustl(varname))//'()') end subroutine load_era5_netcdf_2d subroutine preproc_3d_var(dummy_in, ecmwf, preproc_dims, preproc_geoloc, out_arr3d) - use preproc_constants_m - use preproc_structures_m - - real(sreal), dimension(:,:,:), intent(in) :: dummy_in - type(ecmwf_t), intent(in) :: ecmwf - type(preproc_dims_t), intent(in) :: preproc_dims - type(preproc_geoloc_t), intent(in) :: preproc_geoloc - real(sreal), pointer, intent(inout) :: out_arr3d(:,:, :) - - integer(lint), dimension(1) :: intv, old_grib, new_grib - integer(lint), external :: INTIN, INTOUT, INTF - integer(lint), parameter :: BUFFER = 2000000 - character(len=20), dimension(1) :: charv - real(dreal) :: grid(2), area(4) - real(dreal), allocatable, dimension(:) :: old_data, new_data + use preproc_constants_m + use preproc_structures_m + + real(sreal), dimension(:,:,:), intent(in) :: dummy_in + type(ecmwf_t), intent(in) :: ecmwf + type(preproc_dims_t), intent(in) :: preproc_dims + type(preproc_geoloc_t), intent(in) :: preproc_geoloc + real(sreal), pointer, intent(inout) :: out_arr3d(:,:, :) + + integer(lint), dimension(1) :: intv, old_grib, new_grib +#ifdef INCLUDE_EMOS + integer(lint), external :: INTIN, INTOUT, INTF + integer(lint), parameter :: BUFFER = 2000000 +#endif + character(len=20), dimension(1) :: charv + real(dreal) :: grid(2), area(4) + real(dreal), allocatable, dimension(:) :: old_data, new_data - integer(4) :: n, ni, nj, i, j, k - integer(4) :: old_len, new_len + integer(4) :: n, ni, nj, i, j, k + integer(4) :: old_len, new_len - n = ecmwf%xdim*ecmwf%ydim +#ifdef INCLUDE_EMOS + n = ecmwf%xdim*ecmwf%ydim - ! input details of new grid (see note in read_ecmwf_grib) - charv(1) = 'yes' - grid(1) = sreal_fill_value - if (INTIN('missingvalue', intv, grid, charv) .ne. 0) & + ! input details of new grid (see note in read_ecmwf_grib) + charv(1) = 'yes' + grid(1) = sreal_fill_value + if (INTIN('missingvalue', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN missingvalue failed.') - charv(1) = 'unpacked' - if (INTIN('form', intv, grid, charv) .ne. 0) & + charv(1) = 'unpacked' + if (INTIN('form', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN form failed.') - if (INTOUT('form', intv, grid, charv) .ne. 0) & + if (INTOUT('form', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT form failed.') - intv(1) = ecmwf%ydim/2 - if (INTIN('regular', intv, grid, charv) .ne. 0) & + intv(1) = ecmwf%ydim/2 + if (INTIN('regular', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN reg failed.') - grid(1) = 0.5 / preproc_dims%dellon - grid(2) = 0.5 / preproc_dims%dellat - if (INTOUT('grid', intv, grid, charv) .ne. 0) & + grid(1) = 0.5 / preproc_dims%dellon + grid(2) = 0.5 / preproc_dims%dellat + if (INTOUT('grid', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) - if (INTOUT('area', intv, area, charv) .ne. 0) & + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) + if (INTOUT('area', intv, area, charv) .ne. 0) & call h_e_e('nc', 'INTOUT area failed.') - ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 - nj = ceiling((area(1)+90.)/grid(2)) - floor((area(3)+90.)/grid(2)) + 1 - - allocate(old_data(BUFFER)) - allocate(new_data(BUFFER)) + ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 + nj = ceiling((area(1)+90.)/grid(2)) - floor((area(3)+90.)/grid(2)) + 1 - do k = 1, ecmwf%kdim - old_len = n - old_data(1:n) = reshape(real(dummy_in(:,:,k), kind=8), [n]) + allocate(old_data(BUFFER)) + allocate(new_data(BUFFER)) - new_len = BUFFER + do k = 1, ecmwf%kdim + old_len = n + old_data(1:n) = reshape(real(dummy_in(:,:,k), kind=8), [n]) - if (INTF(old_grib, old_len, old_data, new_grib, new_len, new_data).ne.0)& - call h_e_e('nc', 'INTF failed.') + new_len = BUFFER - if (new_len .ne. ni*nj) print*, '3D Interpolation grid wrong, ' + if (INTF(old_grib, old_len, old_data, new_grib, new_len, new_data).ne.0)& + call h_e_e('nc', 'INTF failed.') - ! copy data into preprocessing grid - do j = 1, nj, 2 - do i = 1, ni, 2 - out_arr3d(preproc_dims%min_lon+i/2, preproc_dims%min_lat+(nj-j)/2, k) = real(new_data(i+(j-1)*ni), kind=4) - end do - end do - end do + if (new_len .ne. ni*nj) print*, '3D Interpolation grid wrong, ' - deallocate(old_data) - deallocate(new_data) + ! copy data into preprocessing grid + do j = 1, nj, 2 + do i = 1, ni, 2 + out_arr3d(1+i/2,1+(nj-j)/2, k) = real(new_data(i+(j-1)*ni), kind=4) + end do + end do + end do + + deallocate(old_data) + deallocate(new_data) +#else + write(*,*) 'ERROR: preproc_3d_var(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine preproc_3d_var subroutine preproc_2d_var(dummy_in, ecmwf, preproc_dims, preproc_geoloc, out_arr2d) - use preproc_constants_m - use preproc_structures_m - - real(sreal), dimension(:,:), intent(in) :: dummy_in - type(ecmwf_t), intent(in) :: ecmwf - type(preproc_dims_t), intent(in) :: preproc_dims - type(preproc_geoloc_t), intent(in) :: preproc_geoloc - real(sreal), pointer, intent(inout) :: out_arr2d(:,:) - - integer(lint), dimension(1) :: intv, old_grib, new_grib - integer(lint), external :: INTIN, INTOUT, INTF - integer(lint), parameter :: BUFFER = 2000000 - character(len=20), dimension(1) :: charv - real(dreal) :: grid(2), area(4) - real(dreal), allocatable, dimension(:) :: old_data, new_data + use preproc_constants_m + use preproc_structures_m + + real(sreal), dimension(:,:), intent(in) :: dummy_in + type(ecmwf_t), intent(in) :: ecmwf + type(preproc_dims_t), intent(in) :: preproc_dims + type(preproc_geoloc_t), intent(in) :: preproc_geoloc + real(sreal), pointer, intent(inout) :: out_arr2d(:,:) + + integer(lint), dimension(1) :: intv, old_grib, new_grib +#ifdef INCLUDE_EMOS + integer(lint), external :: INTIN, INTOUT, INTF + integer(lint), parameter :: BUFFER = 2000000 +#endif + character(len=20), dimension(1) :: charv + real(dreal) :: grid(2), area(4) + real(dreal), allocatable, dimension(:) :: old_data, new_data - integer(4) :: n, ni, nj, i, j - integer(4) :: old_len, new_len + integer(4) :: n, ni, nj, i, j + integer(4) :: old_len, new_len - n = ecmwf%xdim*ecmwf%ydim +#ifdef INCLUDE_EMOS + n = ecmwf%xdim*ecmwf%ydim - ! input details of new grid (see note in read_ecmwf_grib) - charv(1) = 'yes' - grid(1) = sreal_fill_value - if (INTIN('missingvalue', intv, grid, charv) .ne. 0) & + ! input details of new grid (see note in read_ecmwf_grib) + charv(1) = 'yes' + grid(1) = sreal_fill_value + if (INTIN('missingvalue', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN missingvalue failed.') - charv(1) = 'unpacked' - if (INTIN('form', intv, grid, charv) .ne. 0) & + charv(1) = 'unpacked' + if (INTIN('form', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN form failed.') - if (INTOUT('form', intv, grid, charv) .ne. 0) & + if (INTOUT('form', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT form failed.') - intv(1) = ecmwf%ydim/2 - if (INTIN('regular', intv, grid, charv) .ne. 0) & + intv(1) = ecmwf%ydim/2 + if (INTIN('regular', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTIN reg failed.') - grid(1) = 0.5 / preproc_dims%dellon - grid(2) = 0.5 / preproc_dims%dellat - if (INTOUT('grid', intv, grid, charv) .ne. 0) & + grid(1) = 0.5 / preproc_dims%dellon + grid(2) = 0.5 / preproc_dims%dellat + if (INTOUT('grid', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) - if (INTOUT('area', intv, area, charv) .ne. 0) & + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) + if (INTOUT('area', intv, area, charv) .ne. 0) & call h_e_e('nc', 'INTOUT area failed.') - ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 - nj = ceiling((area(1)+90.)/grid(2)) - floor((area(3)+90.)/grid(2)) + 1 + ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 + nj = ceiling((area(1)+90.)/grid(2)) - floor((area(3)+90.)/grid(2)) + 1 - allocate(old_data(BUFFER)) - allocate(new_data(BUFFER)) + allocate(old_data(BUFFER)) + allocate(new_data(BUFFER)) - old_len = n - old_data(1:n) = reshape(real(dummy_in(:,:), kind=8), [n]) + old_len = n + old_data(1:n) = reshape(real(dummy_in(:,:), kind=8), [n]) - new_len = BUFFER + new_len = BUFFER - if (INTF(old_grib, old_len, old_data, new_grib, new_len, new_data).ne.0)& + if (INTF(old_grib, old_len, old_data, new_grib, new_len, new_data).ne.0)& call h_e_e('nc', 'INTF failed.') - if (new_len .ne. ni*nj) print*, '2D Interpolation grid wrong, ' + if (new_len .ne. ni*nj) print*, '2D Interpolation grid wrong, ' - ! copy data into preprocessing grid - do j = 1, nj, 2 - do i = 1, ni, 2 - out_arr2d(preproc_dims%min_lon+i/2, preproc_dims%min_lat+(nj-j)/2) = real(new_data(i+(j-1)*ni), kind=4) - end do - end do + ! copy data into preprocessing grid + do j = 1, nj, 2 + do i = 1, ni, 2 + out_arr2d(1+i/2,1+(nj-j)/2) = real(new_data(i+(j-1)*ni), kind=4) + end do + end do - deallocate(old_data) - deallocate(new_data) + deallocate(old_data) + deallocate(new_data) +#else + write(*,*) 'ERROR: preproc_2d_var(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine preproc_2d_var diff --git a/pre_processing/read_gfs_grib.F90 b/pre_processing/read_gfs_grib.F90 index 8d85d3a0..2e4f7660 100644 --- a/pre_processing/read_gfs_grib.F90 +++ b/pre_processing/read_gfs_grib.F90 @@ -51,6 +51,7 @@ ! over high altitude land regions (Tibet, f.ex) (ExtWork) ! 2017/03/30, SP: Add ability to calculate tropospheric cloud emissivity (ExtWork) ! 2017/06/21, OS: line continuation symbol set to & +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! - If you're having problems with INTF, set the environment variable JDCNDBG=1 @@ -72,8 +73,10 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & type(preproc_prtm_t), intent(inout) :: preproc_prtm logical, intent(in) :: verbose +#ifdef INCLUDE_EMOS integer(lint), parameter :: BUFFER = 3000000 integer(lint), external :: INTIN,INTOUT,INTF2 +#endif integer(lint) :: fu,stat,nbytes integer(lint) :: out_bytes, out_words integer(lint), allocatable, dimension(:) :: in_data,out_data @@ -92,6 +95,7 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & integer(lint),dimension(31) :: gfs_levlist +#ifdef INCLUDE_EMOS gfs_levlist = (/1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,350, & 400,450,500,550,600,650,700,750,800,850,900,925,950,975,1000/) @@ -122,10 +126,10 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & grid(2) = 0.5 / preproc_dims%dellat if (INTOUT('grid',iblank,grid,charv) .ne. 0) & call h_e_e('grib', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) if (INTOUT('area',iblank,area,charv) .ne. 0) & call h_e_e('grib', 'INTOUT area failed.') @@ -165,8 +169,6 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & call h_e_e('grib', 'INTF2 failed.') out_words = out_bytes/lint -! print*,shape(out_data) -! print*,out_words,out_bytes,lint ! stop ! load grib data into grib_api call grib_new_from_message(gid,out_data(1:out_bytes),stat) @@ -203,8 +205,8 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & if (any(level .eq. gfs_levlist) .and. & trim(ltype) .eq. 'isobaricInhPa') then array => preproc_prtm%temperature( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,tlev) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,tlev) preproc_prtm%pressure(:,:,tlev)=level tlev=tlev+1 else @@ -215,58 +217,58 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & trim(ltype) .ne. 'isobaricInhPa') cycle ! Relative humidity array => preproc_prtm%spec_hum( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,qlev) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,qlev) qlev=qlev+1 case(156) if (all(level .ne. gfs_levlist) .or. & trim(ltype) .ne. 'isobaricInhPa') cycle ! Geopotential array => preproc_prtm%phi_lev( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,glev) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,glev) glev=glev+1 case(260131) ! Ozone if (all(level .ne. gfs_levlist) .or. & trim(ltype) .ne. 'isobaricInhPa') cycle array => preproc_prtm%ozone( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat,olev) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim,olev) olev=olev+1 case(134) array => preproc_prtm%lnsp( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(31) array => preproc_prtm%sea_ice_cover( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(3066) array => preproc_prtm%snow_depth( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(165) array => preproc_prtm%u10( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(166) array => preproc_prtm%v10( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(167) array => preproc_prtm%temp2( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(172) array => preproc_prtm%land_sea_mask( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case(54) if (trim(ltype) .ne. 'tropopause') cycle array => preproc_prtm%trop_p( & - preproc_dims%min_lon:preproc_dims%max_lon, & - preproc_dims%min_lat:preproc_dims%max_lat) + 1:preproc_dims%xdim, & + 1:preproc_dims%ydim) case default cycle end select @@ -314,6 +316,10 @@ subroutine read_gfs_grib(ecmwf_file,preproc_dims,preproc_geoloc, & ! Refactor all the GFS levels so that below-surface contributions are removed. call sort_gfs_levels(preproc_prtm,verbose) +#else + write(*,*) 'ERROR: read_gfs_grib(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine read_gfs_grib @@ -321,7 +327,7 @@ end subroutine read_gfs_grib ! This function transforms the GFS fixed pressure levels into surface-relative ! levels that are more similar to those from ECMWF. Needed to prevent below- ! surface contributions to the transmission and radiances. -subroutine sort_gfs_levels(preproc_prtm,verbose) +subroutine sort_gfs_levels(preproc_prtm, verbose) use preproc_constants_m use preproc_structures_m @@ -331,14 +337,14 @@ subroutine sort_gfs_levels(preproc_prtm,verbose) type(preproc_prtm_t), intent(inout) :: preproc_prtm logical, intent(in) :: verbose - integer :: sh(3),lb(3),ub(3),i_0,i_1,j_0,j_1,nl - integer :: i,j,l,stoplev + integer :: sh(3), lb(3), ub(3), i_0, i_1, j_0, j_1, nl + integer :: i, j, l, stoplev - real(dreal) :: surfp,interp - real,allocatable :: p(:),t(:),q(:),o(:),pl(:) + real(dreal) :: surfp, interp + real,allocatable :: p(:), t(:), q(:), o(:), pl(:) logical :: stopper - if (verbose)write(*,*)">>>>>>Sort_gfs_levels>>>>>>" + if (verbose) write(*,*) ">>>>>>Sort_gfs_levels>>>>>>" ! Get the array bounds sh = shape(preproc_prtm%pressure) diff --git a/pre_processing/read_gfs_nc.F90 b/pre_processing/read_gfs_nc.F90 index 8ac26a15..14b09ee4 100644 --- a/pre_processing/read_gfs_nc.F90 +++ b/pre_processing/read_gfs_nc.F90 @@ -25,6 +25,7 @@ ! ! History: ! 2017/07/20, SP: Initial version, cloned from read_ecmwf_nc.F90 +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! - you need to be careful with parameter naming as the variable names are not @@ -48,8 +49,10 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & logical, intent(in) :: verbose integer, intent(in) :: nwp_flag +#ifdef INCLUDE_EMOS integer(lint), external :: INTIN, INTOUT, INTF integer(lint), parameter :: BUFFER = 3000000 +#endif integer(lint), dimension(1) :: intv, old_grib, new_grib real(dreal) :: grid(2), area(4) @@ -67,6 +70,7 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & integer(lint), dimension(31) :: gfs_levlist +#ifdef INCLUDE_EMOS gfs_levlist = (/1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,350, & 400,450,500,550,600,650,700,750,800,850,900,925,950,975,1000/) @@ -93,10 +97,10 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & grid(2) = 0.5 / preproc_dims%dellat if (INTOUT('grid', intv, grid, charv) .ne. 0) & call h_e_e('nc', 'INTOUT grid failed.') - area(1) = preproc_geoloc%latitude(preproc_dims%max_lat) + 0.01*grid(2) - area(2) = preproc_geoloc%longitude(preproc_dims%min_lon) + 0.01*grid(1) - area(3) = preproc_geoloc%latitude(preproc_dims%min_lat) + 0.01*grid(2) - area(4) = preproc_geoloc%longitude(preproc_dims%max_lon) + 0.01*grid(1) + area(1) = preproc_geoloc%latitude(preproc_dims%ydim) + 0.01*grid(2) + area(2) = preproc_geoloc%longitude(1) + 0.01*grid(1) + area(3) = preproc_geoloc%latitude(1) + 0.01*grid(2) + area(4) = preproc_geoloc%longitude(preproc_dims%xdim) + 0.01*grid(1) if (INTOUT('area', intv, area, charv) .ne. 0) & call h_e_e('nc', 'INTOUT area failed.') ni = ceiling((area(4)+180.)/grid(1)) - floor((area(2)+180.)/grid(1)) + 1 @@ -186,8 +190,7 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & ! copy data into preprocessing grid do j = 1, nj, 2 do i = 1, ni, 2 - array3d(preproc_dims%min_lon+i/2, & - preproc_dims%min_lat+(nj-j)/2, k) = & + array3d(1+i/2,1+(nj-j)/2, k) = & real(new_data(i+(j-1)*ni), kind=4) end do end do @@ -208,8 +211,7 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & ! copy data into preprocessing grid do j = 1, nj, 2 do i = 1, ni, 2 - array2d(preproc_dims%min_lon+i/2, & - preproc_dims%min_lat+(nj-j)/2) = & + array2d(1+i/2,1+(nj-j)/2) = & real(new_data(i+(j-1)*ni), kind=4) end do end do @@ -238,5 +240,9 @@ subroutine read_gfs_nc(nwp_path, ecmwf, preproc_dims, preproc_geoloc, & call sort_gfs_levels(preproc_prtm, verbose) call ncdf_close(fid, 'read_gfs_nc()') +#else + write(*,*) 'ERROR: read_gfs_nc(): LIBEMOS is required for ' // & + 'use_ecmwf_preproc_grid = .false.' +#endif end subroutine read_gfs_nc diff --git a/pre_processing/rearrange_ecmwf.F90 b/pre_processing/rearrange_ecmwf.F90 index 3880d106..e56c7535 100644 --- a/pre_processing/rearrange_ecmwf.F90 +++ b/pre_processing/rearrange_ecmwf.F90 @@ -26,18 +26,21 @@ ! 2014/11/04, OS: Added skin temperature. ! 2015/11/17, OS: Added rearrangement of high resolution ERA-Interim data. ! 2018/07/26, AP: Switch to dynamic allocation to reduce stack requirements. +! 2024/07/01, DH: Keep dateline information for ecmwf grid, needed when using +! ecmwf grid as preproc grid ! ! Bugs: ! None known. !------------------------------------------------------------------------------- -subroutine rearrange_ecmwf(ecmwf) +subroutine rearrange_ecmwf(ecmwf, date, ind) implicit none type(ecmwf_t), intent(inout) :: ecmwf + integer, intent(out) :: date, ind - integer :: date, ind, i + integer :: i real(kind=sreal), allocatable, dimension(:,:) :: u, v real(kind=sreal), allocatable, dimension(:,:) :: skint, snow_depth real(kind=sreal), allocatable, dimension(:,:) :: sea_ice_cover @@ -101,3 +104,41 @@ subroutine rearrange_ecmwf(ecmwf) deallocate(lat) end subroutine rearrange_ecmwf + +subroutine rearrange_ecmwf_var2d(ecmwf, dummy2d, date, ind) + + implicit none + + type(ecmwf_t), intent(inout) :: ecmwf + real(sreal), intent(inout) :: dummy2d(ecmwf%xdim,ecmwf%ydim) + integer, intent(in) :: date, ind + real(sreal) :: dummy2d_new(ecmwf%xdim,ecmwf%ydim) + integer :: i + + dummy2d_new(1:ind,:) = dummy2d(date:,:) + dummy2d_new(ind+1:,:)= dummy2d(1:date-1,:) + + do i = 1, ecmwf%ydim + dummy2d(:,ecmwf%ydim+1-i) = dummy2d_new(:,i) + end do + +end subroutine rearrange_ecmwf_var2d + +subroutine rearrange_ecmwf_var3d(ecmwf, dummy3d, date, ind) + + implicit none + + type(ecmwf_t), intent(inout) :: ecmwf + real(sreal), intent(inout) :: dummy3d(ecmwf%xdim,ecmwf%ydim, ecmwf%kdim) + integer, intent(in) :: date, ind + real(sreal) :: dummy3d_new(ecmwf%xdim,ecmwf%ydim, ecmwf%kdim) + integer :: i + + dummy3d_new(1:ind,:, :) = dummy3d(date:,:, :) + dummy3d_new(ind+1:,:, :)= dummy3d(1:date-1,:, :) + + do i = 1, ecmwf%ydim + dummy3d(:,ecmwf%ydim+1-i, :) = dummy3d_new(:,i, :) + end do + +end subroutine rearrange_ecmwf_var3d diff --git a/pre_processing/rttov_driver.F90 b/pre_processing/rttov_driver.F90 index 34a8fab2..4abc3cce 100644 --- a/pre_processing/rttov_driver.F90 +++ b/pre_processing/rttov_driver.F90 @@ -159,6 +159,7 @@ ! RTTOV documentation). Removed declaration on nevals ! as it is not used anywhere else (and never had a ! value assigned to it)! +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! - BRDF not yet implemented here, so RTTOV internal calculation used. @@ -559,7 +560,6 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & write(*,*) 'ERROR: rttov_alloc_prof(), errorstatus = ', stat stop error_stop_code end if - profiles%id = 'standard' ! Compute the appropriate CO2 value for this scene @@ -572,10 +572,9 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & ! Copy preprocessing grid data into RTTOV profile structure ! Create a lowest layer from the surface properties count = 0 - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 - ! Check to see if the ECMWF data read in includes ozone ! profiles: forecast data does not if (maxval(preproc_prtm%ozone(idim,jdim,:)) == 0.0) then @@ -593,7 +592,6 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & profiles(count)%t(:nlayers) = preproc_prtm%temperature(idim,jdim,:) profiles(count)%q(:nlayers) = preproc_prtm%spec_hum(idim,jdim,:) profiles(count)%o3(:nlayers) = preproc_prtm%ozone(idim,jdim,:) - ! Add CO2 in kg/kg for each level if (pre_opts%do_co2) profiles(count)%co2(:) = co2_val @@ -640,8 +638,8 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & ! Write profiles structure to PRTM file (array operations needed to ! recast structure in form ncdf_write_array recognises) - i_ = idim - preproc_dims%min_lon + 1 - j_ = jdim - preproc_dims%min_lat + 1 + i_ = idim + j_ = jdim call ncdf_write_array(netcdf_info%ncid_prtm, 'lon_rtm', & netcdf_info%vid_lon_pw, & (/profiles(count)%longitude/), & @@ -688,8 +686,8 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & if (verbose) write(*,*) ' - Calculating for viewing geometry number', cview count = 0 - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 profiles(count)%zenangle = preproc_geo%satza(idim,jdim,cview) @@ -822,11 +820,9 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & #else if (verbose) write(*,*) 'Run RTTOV' #endif - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon - + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 - ! Process points that contain information and satisfy the zenith ! angle restrictions of the coefficient file if ((i_coef == 1 .and. & @@ -835,7 +831,6 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & (i_coef == 2 .and. & preproc_dims%counter_sw(idim,jdim,cview) > 0 .and. & profiles(count)%zenangle <= zenmaxv9)) then - ! Fetch emissivity from atlas call rttov_get_emis(stat, opts, chanprof, & profiles(count:count), coefs, emis_atlas, emis_data) @@ -845,7 +840,6 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & stop error_stop_code end if emissivity%emis_in = emis_data - ! Fetch emissivity from the MODIS CIMSS emissivity product if (i_coef == 1 .and. pre_opts%use_modis_emis_in_rttov) then where (preproc_surf%emissivity(idim,jdim,:) /= & @@ -854,9 +848,7 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & preproc_surf%emissivity(idim,jdim,chan_pos) end where end if - calcemis = emissivity%emis_in <= dither - if (preproc_dims%counter_lw(idim,jdim,cview) .gt. 0) then ! Call RTTOV for this profile #ifdef INCLUDE_RTTOV_OPENMP @@ -875,7 +867,6 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & write(*,*) 'ERROR: rttov_direct(), errorstatus = ', stat stop error_stop_code end if - ! Remove the Rayleigh component from the RTTOV tranmittances only if RTTOV ! ran with Rayleigh scattering switched on. if (i_coef == 2) then @@ -895,16 +886,16 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & if (i_coef == 1) then do i_ = 1, nchan call write_ir_rttov(netcdf_info, & - idim-preproc_dims%min_lon+1, & - jdim-preproc_dims%min_lat+1, & + idim, & + jdim, & profiles(count)%nlevels, emissivity, transmission, & radiance, radiance2, write_rttov, chan_pos(i_), i_) end do else do i_ = 1, nchan call write_solar_rttov(netcdf_info, coefs, & - idim-preproc_dims%min_lon+1, & - jdim-preproc_dims%min_lat+1, & + idim, & + jdim, & profiles(count)%nlevels, profiles(count)%zenangle, & transmission, write_rttov, chan_pos(i_), i_) end do @@ -921,8 +912,8 @@ subroutine rttov_driver(coef_path, emiss_path, granule, preproc_dims, & #else if (verbose) write(*,*) 'Run RTTOV for cloud' #endif - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 profiles(count)%cfraction = 1. profiles(count)%ctp = preproc_prtm%trop_p(idim,jdim) diff --git a/pre_processing/rttov_driver_gfs.F90 b/pre_processing/rttov_driver_gfs.F90 index fa280db6..cb2139d2 100644 --- a/pre_processing/rttov_driver_gfs.F90 +++ b/pre_processing/rttov_driver_gfs.F90 @@ -61,6 +61,7 @@ ! with the allocation statements for these arrays and ! RTTOV documentation. (nevals was not used anywhere ! else and wasn't even initialised to any value!) +! 2024/07/01, DH: Change indexing to use preproc_dims for all dimensions ! ! Bugs: ! - BRDF not yet implemented here, so RTTOV internal calculation used. @@ -474,8 +475,8 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & ! Copy preprocessing grid data into RTTOV profile structure ! Create a lowest layer from the surface properties count = 0 - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 ! set gas units to 1, specifying gas input in kg/kg @@ -530,8 +531,8 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & ! Write profiles structure to PRTM file (array operations needed to ! recast structure in form ncdf_write_array recognises) - i_ = idim - preproc_dims%min_lon + 1 - j_ = jdim - preproc_dims%min_lat + 1 + i_ = idim + j_ = jdim call ncdf_write_array(netcdf_info%ncid_prtm, 'lon_rtm', & netcdf_info%vid_lon_pw, & (/profiles(count)%longitude/), & @@ -577,8 +578,8 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & if (verbose) write(*,*) ' - Calculating for viewing geometry number', cview count = 0 - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 profiles(count)%zenangle = preproc_geo%satza(idim,jdim,cview) profiles(count)%azangle = preproc_geo%satazi(idim,jdim,cview) @@ -705,8 +706,8 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & if ((verbose) .and. i_coef .eq. 1) write(*,*) 'Run RTTOV Longwave' if ((verbose) .and. i_coef .eq. 2) write(*,*) 'Run RTTOV Shortwave' #endif - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 ! Process points that contain information and satisfy the zenith @@ -770,16 +771,16 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & if (i_coef == 1) then do i_ = 1, nchan call write_ir_rttov(netcdf_info, & - idim-preproc_dims%min_lon+1, & - jdim-preproc_dims%min_lat+1, & + idim, & + jdim, & profiles(count)%nlevels, emissivity, transmission, & radiance, radiance2, write_rttov, chan_pos(i_), i_) end do else do i_ = 1, nchan call write_solar_rttov(netcdf_info, coefs, & - idim-preproc_dims%min_lon+1, & - jdim-preproc_dims%min_lat+1, & + idim, & + jdim, & profiles(count)%nlevels, profiles(count)%zenangle, & transmission, write_rttov, chan_pos(i_), i_) end do @@ -797,8 +798,8 @@ subroutine rttov_driver_gfs(coef_path, emiss_path, granule, preproc_dims, & #else if (verbose) write(*,*) 'Run RTTOV for cloud' #endif - do jdim = preproc_dims%min_lat, preproc_dims%max_lat - do idim = preproc_dims%min_lon, preproc_dims%max_lon + do jdim = 1, preproc_dims%ydim + do idim = 1, preproc_dims%xdim count = count + 1 profiles(count)%cfraction = 1. profiles(count)%ctp = preproc_prtm%trop_p(idim,jdim) diff --git a/pre_processing/utils_for_main.F90 b/pre_processing/utils_for_main.F90 index a044deee..438007ad 100644 --- a/pre_processing/utils_for_main.F90 +++ b/pre_processing/utils_for_main.F90 @@ -33,6 +33,7 @@ ! USE_GSICS enables this to be disabled. ! 2018/08/30, SP: Allow variable CO2 in RTTOV, linear scaling from 2006 value ! 2021/12/14, DP: Added SEVIRI external ANN option +! 2024/07/01, DH: Added option to use native ECMWF grid for preprocessing ! ! Bugs: ! None known. @@ -181,6 +182,9 @@ subroutine parse_optional(label, value, preproc_opts) case('USE_SEVIRI_ANN_MLAY') if (parse_string(value, preproc_opts%use_seviri_ann_mlay) /= 0) & call handle_parse_error(label) + case('USE_ECMWF_PREPROC_GRID') + if (parse_string(value, preproc_opts%use_ecmwf_preproc_grid) /= 0) & + call handle_parse_error(label) case default write(*,*) 'ERROR: Unknown option: ', trim(label) stop error_stop_code diff --git a/src/orac_constants.F90 b/src/orac_constants.F90 index ede9c881..45b4d5cd 100644 --- a/src/orac_constants.F90 +++ b/src/orac_constants.F90 @@ -109,7 +109,7 @@ module ORAC_constants_m integer, parameter :: MaxNumThermal = 16 ! Max no. of thermal channels integer, parameter :: MaxCloudType = 5 ! Max no. of cloud types to be integer, parameter :: MaxCRProps = 14 ! Max no. of properties in SAD_LUT arrays - integer, parameter :: MaxTypes = 11 ! Number of possible cloud/aerosol types + integer, parameter :: MaxTypes = 13 ! Number of possible cloud/aerosol types ! Tolerance values diff --git a/src/read_driver.F90 b/src/read_driver.F90 index 392fc177..704d4724 100644 --- a/src/read_driver.F90 +++ b/src/read_driver.F90 @@ -670,7 +670,7 @@ subroutine Read_Driver(Ctrl, global_atts, source_atts) if (Ctrl%verbose) write(*,*) 'WARNING: Read_Driver(): '// & 'Ctrl%use_ann_phase=true with Ctrl%Approach=AppCld2L '// & 'processes all pixels.' - Ctrl%NTypes_to_process = 11 + Ctrl%NTypes_to_process = 13 Ctrl%Types_to_process(1) = CLEAR_TYPE Ctrl%Types_to_process(2) = SWITCHED_TO_WATER_TYPE Ctrl%Types_to_process(3) = FOG_TYPE @@ -682,6 +682,8 @@ subroutine Read_Driver(Ctrl, global_atts, source_atts) Ctrl%Types_to_process(9) = OVERLAP_TYPE Ctrl%Types_to_process(10) = PROB_OPAQUE_ICE_TYPE Ctrl%Types_to_process(11) = PROB_CLEAR_TYPE + Ctrl%Types_to_process(12) = DUST_CLEAR_TYPE + Ctrl%Types_to_process(13) = DUST_SWITCHED_FROM_CLOUD_TYPE else if (.not. Ctrl%use_ann_phase .and. Ctrl%Approach == AppCld2L) then Ctrl%NTypes_to_process = 1 diff --git a/tools/pyorac/definitions.py b/tools/pyorac/definitions.py index 73e950d1..1eb8bd62 100644 --- a/tools/pyorac/definitions.py +++ b/tools/pyorac/definitions.py @@ -7,7 +7,7 @@ # Names of the possible Pavolonis cloud classes ALL_TYPES = ('CLEAR', 'SWITCHED_TO_WATER', 'FOG', 'WATER', 'SUPERCOOLED', 'SWITCHED_TO_ICE', 'OPAQUE_ICE', 'CIRRUS', 'OVERLAP', - 'PROB_OPAQUE_ICE', 'PROB_CLEAR') + 'PROB_OPAQUE_ICE', 'PROB_CLEAR', 'DUST_CLEAR', 'DUST_SWITCHED_FROM_CLOUD') # Colours used when printing to screen COLOURING = {