From f282c66a564909716b7043c2bf70af3ae6ad2dd0 Mon Sep 17 00:00:00 2001 From: Rodrigo Guzman Date: Fri, 2 Aug 2019 18:54:30 +0200 Subject: [PATCH] Fixing CloudSat Precip flag surface elevation bug (#19) * 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 --- driver/src/cosp2_test.f90 | 18 ++++- driver/test_cosp2Imp.py | 4 +- src/cosp.F90 | 6 +- src/simulator/quickbeam/quickbeam.F90 | 105 ++++++++++++++++---------- 4 files changed, 86 insertions(+), 47 deletions(-) diff --git a/driver/src/cosp2_test.f90 b/driver/src/cosp2_test.f90 index 68cf419d9d..239d003e67 100644 --- a/driver/src/cosp2_test.f90 +++ b/driver/src/cosp2_test.f90 @@ -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 :: & @@ -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 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/driver/test_cosp2Imp.py b/driver/test_cosp2Imp.py index 846c572134..1dcb8dce77 100644 --- a/driver/test_cosp2Imp.py +++ b/driver/test_cosp2Imp.py @@ -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) diff --git a/src/cosp.F90 b/src/cosp.F90 index 12045a6f06..92580f65ff 100644 --- a/src/cosp.F90 +++ b/src/cosp.F90 @@ -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 @@ -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) diff --git a/src/simulator/quickbeam/quickbeam.F90 b/src/simulator/quickbeam/quickbeam.F90 index cede9a49dd..54a4fd8152 100644 --- a/src/simulator/quickbeam/quickbeam.F90 +++ b/src/simulator/quickbeam/quickbeam.F90 @@ -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 @@ -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) :: & @@ -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) @@ -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. @@ -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 @@ -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) :: & @@ -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 @@ -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 @@ -439,25 +455,34 @@ 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 @@ -465,11 +490,11 @@ subroutine cloudsat_precipOccurence(Npoints, Ncolumns, llm, Nhydro, Ze_out, Ze_n ! Mized phase (273275) 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 @@ -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 @@ -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