Skip to content

Commit

Permalink
Fixing CloudSat Precip flag surface elevation bug (#19)
Browse files Browse the repository at this point in the history 
* bug fix on OPAQ variables

* DEBUG comments erased

* CloudSat Precip flag surface elevation bug fixed

* Fixing CloudSat Precip flag surface elevation bug

* CloudSat Precip flag bug fix reviewed, magic number 480 replaced by 'zstep' as defined in 'cosp_init', but replacing the parameter cloudsat_preclvl=39 by a more general expression is out of the scope of this bug fix

* For land ponits, use 2 levels above the surface for retreival. This is consistent w/ the operational product and what Ken Kay did in DOI:10.1002/2017JD028213

* Change index adjustment for land.

* Added computation of path integrated attenutation for both land/ocean points.

* Added lower limit on reflectivty for no precipitation class.

* Add temporary diagnostics.

* Some small changes.

* Change to indexing.

* Reverted some changes

* Removed lower bound on unclassified precipitation over land, apparently reflectivities can go that low, something I did not know.

* Added comment
  • Loading branch information
@rodrigoguzman-lmd @dustinswales
rodrigoguzman-lmd authored and dustinswales committed Aug 2, 2019
1 parent 4915173 commit f282c66
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 47 deletions.
18 changes: 15 additions & 3 deletions driver/src/cosp2_test.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -547,7 +547,9 @@ subroutine subsample_and_optics(nPoints, nLevels, nColumns, nHydro, overlap,
! Local variables
type(rng_state),allocatable,dimension(:) :: rngs ! Seeds for random number generator
integer,dimension(:),allocatable :: seed
integer,dimension(:),allocatable :: cloudsat_preclvl_index
integer :: i,j,k
real(wp) :: zstep
real(wp),dimension(:,:), allocatable :: &
ls_p_rate, cv_p_rate, frac_ls, frac_cv, prec_ls, prec_cv,g_vol
real(wp),dimension(:,:,:), allocatable :: &
Expand Down Expand Up @@ -857,11 +859,21 @@ subroutine subsample_and_optics(nPoints, nLevels, nColumns, nHydro, overlap,
fracPrecipIce_statGrid(:,:,:) = 0._wp
call cosp_change_vertical_grid(Npoints, Ncolumns, Nlevels, cospstateIN%hgt_matrix(:,Nlevels:1:-1), &
cospstateIN%hgt_matrix_half(:,Nlevels:1:-1), fracPrecipIce(:,:,Nlevels:1:-1), Nlvgrid_local, &
vgrid_zl(Nlvgrid_local:1:-1), vgrid_zu(Nlvgrid_local:1:-1), fracPrecipIce_statGrid)
vgrid_zl(Nlvgrid_local:1:-1), vgrid_zu(Nlvgrid_local:1:-1), fracPrecipIce_statGrid(:,:,Nlvgrid_local:1:-1))

! Find proper layer above de surface elevation to compute precip flags in Cloudsat/Calipso statistical grid
allocate(cloudsat_preclvl_index(nPoints))
cloudsat_preclvl_index(:) = 0._wp
! Compute the zstep distance between two atmopsheric layers
zstep = vgrid_zl(1)-vgrid_zl(2)
! Computing altitude index for precip flags calculation (one layer above surfelev layer)
cloudsat_preclvl_index(:) = cloudsat_preclvl - floor( cospstateIN%surfelev(:)/zstep )

! For near-surface diagnostics, we only need the frozen fraction at one layer.
cospIN%fracPrecipIce(:,:) = fracPrecipIce_statGrid(:,:,cloudsat_preclvl)

do i=1,nPoints
cospIN%fracPrecipIce(i,:) = fracPrecipIce_statGrid(i,:,cloudsat_preclvl_index(i))
enddo

endif

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expand Down
4 changes: 3 additions & 1 deletion driver/test_cosp2Imp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,9 @@
# Compute relative difference.
diff = var/var_ref-1
diff[numpy.where(var_ref == 0.0)] = 0.0

# Relative difference = 1 if var_ref = 0 and var != var_ref
diff[(var_ref == 0.0) & (var != var_ref)] = 1.0

# Flag any point which exceeds error threshold.
error = numpy.argwhere(abs(diff) >= args.zeroThresh)

Expand Down
6 changes: 3 additions & 3 deletions src/cosp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ MODULE MOD_COSP
ntau,modis_histTau,tau_binBounds, &
modis_histTauEdges,tau_binEdges,nCloudsatPrecipClass,&
modis_histTauCenters,tau_binCenters, &
cloudsat_preclvl,grLidar532_histBsct,atlid_histBsct
cloudsat_preclvl,grLidar532_histBsct,atlid_histBsct
USE MOD_COSP_MODIS_INTERFACE, ONLY: cosp_modis_init, modis_IN
USE MOD_COSP_RTTOV_INTERFACE, ONLY: cosp_rttov_init, rttov_IN
USE MOD_COSP_MISR_INTERFACE, ONLY: cosp_misr_init, misr_IN
Expand Down Expand Up @@ -1241,8 +1241,8 @@ function COSP_SIMULATOR(cospIN,cospgridIN,cospOUT,start_idx,stop_idx,debug)
! Call simulator
call quickbeam_column(cloudsatIN%Npoints, cloudsatIN%Ncolumns, cloudsatIN%Nlevels,&
Nlvgrid, cloudsat_DBZE_BINS, 'cloudsat', cloudsatDBZe, cloudsatZe_non, &
cospgridIN%land(:), cospgridIN%at(:,cospIN%Nlevels), cospIN%fracPrecipIce, &
cospgridIN%hgt_matrix, cospgridIN%hgt_matrix_half, &
cospgridIN%land(:), cospgridIN%surfelev(:), cospgridIN%at(:,cospIN%Nlevels), &
cospIN%fracPrecipIce, cospgridIN%hgt_matrix, cospgridIN%hgt_matrix_half, &
cospOUT%cloudsat_cfad_ze(ij:ik,:,:), cospOUT%cloudsat_precip_cover, &
cospOUT%cloudsat_pia)
! Free up memory (if necessary)
Expand Down
105 changes: 65 additions & 40 deletions src/simulator/quickbeam/quickbeam.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ module quickbeam
pClass_Snow1, pClass_Snow2, pClass_Mixed1, pClass_Mixed2, &
pClass_Rain4, pClass_default, Zenonbinval, Zbinvallnd, &
N_HYDRO,nCloudsatPrecipClass,cloudsat_preclvl

USE MOD_COSP_STATS, ONLY: COSP_LIDAR_ONLY_CLOUD,hist1D,COSP_CHANGE_VERTICAL_GRID
implicit none

Expand Down Expand Up @@ -232,7 +233,7 @@ end subroutine quickbeam_subcolumn
! SUBROUTINE quickbeam_column
! ######################################################################################
subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform, &
Ze_tot, Ze_tot_non, land, t2m, fracPrecipIce, zlev, zlev_half, cfad_ze, &
Ze_tot, Ze_tot_non, land, surfelev, t2m, fracPrecipIce, zlev, zlev_half, cfad_ze, &
cloudsat_precip_cover, cloudsat_pia)
! Inputs
integer,intent(in) :: &
Expand All @@ -245,6 +246,7 @@ subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform
platform ! Name of platform (e.g. cloudsat)
real(wp),dimension(Npoints),intent(in) :: &
land, & ! Land/Sea mask. (1/0)
surfelev, & ! Surface Elevation (m)
t2m ! Near-surface temperature
real(wp),dimension(Npoints,Ncolumns),intent(in) :: &
fracPrecipIce ! Fraction of precipitation which is frozen. (1)
Expand All @@ -266,6 +268,7 @@ subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform

! Local variables
integer :: i,j
real(wp) :: zstep
real(wp),dimension(npoints,ncolumns,llm) :: ze_toti,ze_noni
logical :: lcloudsat = .false.

Expand Down Expand Up @@ -294,9 +297,11 @@ subroutine quickbeam_column(npoints, ncolumns, nlevels, llm, DBZE_BINS, platform
call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zlev(:,nlevels:1:-1),&
zlev_half(:,nlevels:1:-1),Ze_tot_non(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),&
vgrid_zu(llm:1:-1),Ze_noni(:,:,llm:1:-1),log_units=.true.)
! Not call routine to generate diagnostics.
! Compute the zstep distance between two atmopsheric layers
zstep = vgrid_zl(1)-vgrid_zl(2)
! Now call routine to generate diagnostics.
call cloudsat_precipOccurence(Npoints, Ncolumns, llm, N_HYDRO, Ze_toti, Ze_noni, &
land, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia)
land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia, zstep)
else
! Effective reflectivity histogram
do i=1,Npoints
Expand Down Expand Up @@ -345,7 +350,7 @@ end subroutine quickbeam_column
! parameter cloudsat_preclvl, defined in src/cosp_config.F90
! ######################################################################################
subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_non_out, &
land, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia)
land, surfelev, t2m, fracPrecipIce, cloudsat_precip_cover, cloudsat_pia, zstep)

! Inputs
integer,intent(in) :: &
Expand All @@ -355,23 +360,29 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
llm ! Number of levels
real(wp),dimension(Npoints),intent(in) :: &
land, & ! Land/Sea mask. (1/0)
surfelev, & ! Surface Elevation (m)
t2m ! Near-surface temperature
real(wp),dimension(Npoints,Ncolumns,llm),intent(in) :: &
Ze_out, & ! Effective reflectivity factor (dBZ)
Ze_non_out ! Effective reflectivity factor, w/o attenuation (dBZ)
real(wp),dimension(Npoints,Ncolumns),intent(in) :: &
fracPrecipIce ! Fraction of precipitation which is frozen. (1)
real(wp),intent(in) :: &
zstep ! Distance between two atmopsheric layers (m)

! Outputs
real(wp),dimension(Npoints,nCloudsatPrecipClass),intent(out) :: &
cloudsat_precip_cover ! Model precip rate in by CloudSat precip flag
real(wp),dimension(Npoints),intent(out) :: &
cloudsat_pia ! Cloudsat path integrated attenuation
cloudsat_pia ! Cloudsat path integrated attenuation

! Local variables
integer,dimension(Npoints,Ncolumns) :: &
cloudsat_pflag, & ! Subcolumn precipitation flag
cloudsat_precip_pia ! Subcolumn path integrated attenutation.
integer,dimension(Npoints) :: &
cloudsat_preclvl_index ! Altitude index for precip flags calculation
! in 40-level grid (one layer above surfelev)
integer :: pr,i,k,m,j
real(wp) :: Zmax

Expand All @@ -380,57 +391,62 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
cloudsat_precip_pia(:,:) = 0._wp
cloudsat_precip_cover(:,:) = 0._wp
cloudsat_pia(:) = 0._wp
cloudsat_preclvl_index(:) = 0._wp

! Computing altitude index for precip flags calculation
cloudsat_preclvl_index(:) = cloudsat_preclvl - floor( surfelev(:)/zstep )

! ######################################################################################
! SUBCOLUMN processing
! ######################################################################################
do i=1, Npoints
do pr=1,Ncolumns
! 1) Compute the PIA in all profiles containing hydrometeors
if ( (Ze_non_out(i,pr,cloudsat_preclvl).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then
if ( (Ze_non_out(i,pr,cloudsat_preclvl).lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then
cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl)
endif
endif

! 2) Compute precipitation flag
! Compute precipitation flag
! ################################################################################
! 2a) Oceanic points.
! 1) Oceanic points.
! ################################################################################
if (land(i) .eq. 0) then

! 1a) Compute the PIA in all profiles containing hydrometeors
if ( (Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl_index(i)).gt.-100) ) then
if ( (Ze_non_out(i,pr,cloudsat_preclvl_index(i)).lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl_index(i)).lt.100) ) then
cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl_index(i)) - Ze_out(i,pr,cloudsat_preclvl_index(i))
endif
endif

! Snow
if(fracPrecipIce(i,pr).gt.0.9) then
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(2)) then
cloudsat_pflag(i,pr) = pClass_Snow2 ! TSL: Snow certain
endif
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl_index(i)).le.Zenonbinval(2)) then
cloudsat_pflag(i,pr) = pClass_Snow1 ! TSL: Snow possible
endif
endif

! Mixed
if(fracPrecipIce(i,pr).gt.0.1.and.fracPrecipIce(i,pr).le.0.9) then
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(2)) then
cloudsat_pflag(i,pr) = pClass_Mixed2 ! TSL: Mixed certain
endif
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl_index(i)).le.Zenonbinval(2)) then
cloudsat_pflag(i,pr) = pClass_Mixed1 ! TSL: Mixed possible
endif
endif

! Rain
if(fracPrecipIce(i,pr).le.0.1) then
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(1)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(1)) then
cloudsat_pflag(i,pr) = pClass_Rain3 ! TSL: Rain certain
endif
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(3).and. &
Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(1)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(3).and. &
Ze_non_out(i,pr,cloudsat_preclvl_index(i)).le.Zenonbinval(1)) then
cloudsat_pflag(i,pr) = pClass_Rain2 ! TSL: Rain probable
endif
if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(3)) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).gt.Zenonbinval(4).and. &
Ze_non_out(i,pr,cloudsat_preclvl_index(i)).le.Zenonbinval(3)) then
cloudsat_pflag(i,pr) = pClass_Rain1 ! TSL: Rain possible
endif
if(cloudsat_precip_pia(i,pr).gt.40) then
Expand All @@ -439,37 +455,46 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
endif

! No precipitation
if(Ze_non_out(i,pr,cloudsat_preclvl).le.-15) then
if(Ze_non_out(i,pr,cloudsat_preclvl_index(i)).le.-15) then
cloudsat_pflag(i,pr) = pClass_noPrecip ! TSL: Not Raining
endif
endif ! Ocean points

! ################################################################################
! 2b) Land points.
! 2) Land points.
! *NOTE* For land points we go up a layer higher, so cloudsat_preclvl_index(i)-1
!
! ################################################################################
if (land(i) .eq. 1) then
if (land(i) .eq. 1) then
! 2a) Compute the PIA in all profiles containing hydrometeors
if ( (Ze_non_out(i,pr,cloudsat_preclvl_index(i)-1).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl_index(i)-1).gt.-100) ) then
if ( (Ze_non_out(i,pr,cloudsat_preclvl_index(i)-1).lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl_index(i)-1).lt.100) ) then
cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl_index(i)-1) - Ze_out(i,pr,cloudsat_preclvl_index(i)-1)
endif
endif

! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out);
Zmax=maxval(Ze_out(i,pr,:))

! Snow (T<273)
if(t2m(i) .lt. 273._wp) then
if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(5)) then
if(Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(5)) then
cloudsat_pflag(i,pr) = pClass_Snow2 ! JEK: Snow certain
endif
if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .and. &
Ze_out(i,pr,cloudsat_preclvl).le.Zbinvallnd(5)) then
if(Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(6) .and. &
Ze_out(i,pr,cloudsat_preclvl_index(i)-1).le.Zbinvallnd(5)) then
cloudsat_pflag(i,pr) = pClass_Snow1 ! JEK: Snow possible
endif
endif

! Mized phase (273<T<275)
if(t2m(i) .ge. 273._wp .and. t2m(i) .le. 275._wp) then
if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &
(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(4))) then
(Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(4))) then
cloudsat_pflag(i,pr) = pClass_Mixed2 ! JEK: Mixed certain
endif
if ((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .and. &
Ze_out(i,pr,cloudsat_preclvl) .le. Zbinvallnd(4)) .and. &
if ((Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(6) .and. &
Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .le. Zbinvallnd(4)) .and. &
(Zmax .gt. Zbinvallnd(5)) ) then
cloudsat_pflag(i,pr) = pClass_Mixed1 ! JEK: Mixed possible
endif
Expand All @@ -478,14 +503,14 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
! Rain (T>275)
if(t2m(i) .gt. 275) then
if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &
(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(2))) then
(Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(2))) then
cloudsat_pflag(i,pr) = pClass_Rain3 ! JEK: Rain certain
endif
if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. &
if((Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(6)) .and. &
(Zmax .gt. Zbinvallnd(3))) then
cloudsat_pflag(i,pr) = pClass_Rain2 ! JEK: Rain probable
endif
if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. &
if((Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .gt. Zbinvallnd(6)) .and. &
(Zmax.lt.Zbinvallnd(3))) then
cloudsat_pflag(i,pr) = pClass_Rain1 ! JEK: Rain possible
endif
Expand All @@ -495,7 +520,7 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
endif

! No precipitation
if(Ze_out(i,pr,cloudsat_preclvl).le.-15) then
if(Ze_out(i,pr,cloudsat_preclvl_index(i)-1) .le. -15) then
cloudsat_pflag(i,pr) = pClass_noPrecip ! JEK: Not Precipitating
endif
endif ! Land points
Expand All @@ -514,7 +539,7 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n
cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq. k-1)
endif
enddo

! Gridmean path integrated attenuation (pia)
cloudsat_pia(i)=sum(cloudsat_precip_pia(i,:))
enddo
Expand Down

0 comments on commit f282c66

Please sign in to comment.