Skip to content

Commit 08f4cce

Browse files
dustinswalesdustinswales
authored
Merge pull request #4 from CFMIP/modisLidarBugFix
Modis lidar bug fix
2 parents 7fc6532 + 5b9356b commit 08f4cce

File tree

3 files changed

+119
-98
lines changed

3 files changed

+119
-98
lines changed

‎MISR_simulator/MISR_simulator.f‎

Lines changed: 53 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -283,60 +283,63 @@ SUBROUTINE MISR_simulator(
283283

284284
!
285285
! Modify MISR CTH for satellite spatial / pattern matcher effects
286+
!
287+
! Code in this region added by roj 5/2006 to account
288+
! for spatial effect of the MISR pattern matcher.
289+
! Basically, if a column is found between two neighbors
290+
! at the same CTH, and that column has no hieght or
291+
! a lower CTH, THEN misr will tend to but place the
292+
! odd column at the same height as it neighbors.
293+
!
294+
! This setup assumes the columns represent a about a 1 to 4 km scale
295+
! it will need to be modified significantly, otherwise
296+
297+
! Commented out. Should only be used with subcolumns are not randomly ordered. Roj 4/2018
286298
!
287-
! Code in this region added by roj 5/2006 to account
288-
! for spatial effect of the MISR pattern matcher.
289-
! Basically, if a column is found between two neighbors
290-
! at the same CTH, and that column has no hieght or
291-
! a lower CTH, THEN misr will tend to but place the
292-
! odd column at the same height as it neighbors.
299+
! if(ncol.eq.1) then
300+
!
301+
! adjust based on neightboring points ... i.e. only 2D grid was input
302+
! do j=2,npoints-1
303+
!
304+
! if(box_MISR_ztop(j-1,1).gt.0 .and.
305+
! & box_MISR_ztop(j+1,1).gt.0 ) then
293306
!
294-
! This setup assumes the columns represent a about a 1 to 4 km scale
295-
! it will need to be modified significantly, otherwise
296-
if(ncol.eq.1) then
297-
298-
! adjust based on neightboring points ... i.e. only 2D grid was input
299-
do j=2,npoints-1
300-
301-
if(box_MISR_ztop(j-1,1).gt.0 .and.
302-
& box_MISR_ztop(j+1,1).gt.0 ) then
303-
304-
if( abs( box_MISR_ztop(j-1,1) -
305-
& box_MISR_ztop(j+1,1) ) .lt. 500
306-
& .and.
307-
& box_MISR_ztop(j,1) .lt.
308-
& box_MISR_ztop(j+1,1) ) then
309-
310-
box_MISR_ztop(j,1) =
311-
& box_MISR_ztop(j+1,1)
312-
endif
313-
314-
endif
315-
enddo
316-
else
317-
307+
! if( abs( box_MISR_ztop(j-1,1) -
308+
! & box_MISR_ztop(j+1,1) ) .lt. 500
309+
! & .and.
310+
! & box_MISR_ztop(j,1) .lt.
311+
! & box_MISR_ztop(j+1,1) ) then
312+
!
313+
! box_MISR_ztop(j,1) =
314+
! & box_MISR_ztop(j+1,1)
315+
! endif
316+
!
317+
! endif
318+
! enddo
319+
! else
320+
!
318321
! adjust based on neighboring subcolumns ....
319-
do ibox=2,ncol-1
320-
321-
if(box_MISR_ztop(1,ibox-1).gt.0 .and.
322-
& box_MISR_ztop(1,ibox+1).gt.0 ) then
323-
324-
if( abs( box_MISR_ztop(1,ibox-1) -
325-
& box_MISR_ztop(1,ibox+1) ) .lt. 500
326-
& .and.
327-
& box_MISR_ztop(1,ibox) .lt.
328-
& box_MISR_ztop(1,ibox+1) ) then
329-
330-
box_MISR_ztop(1,ibox) =
331-
& box_MISR_ztop(1,ibox+1)
332-
endif
333-
334-
endif
335-
enddo
336-
337-
endif
322+
! do ibox=2,ncol-1
323+
!
324+
! if(box_MISR_ztop(1,ibox-1).gt.0 .and.
325+
! & box_MISR_ztop(1,ibox+1).gt.0 ) then
326+
!
327+
! if( abs( box_MISR_ztop(1,ibox-1) -
328+
! & box_MISR_ztop(1,ibox+1) ) .lt. 500
329+
! & .and.
330+
! & box_MISR_ztop(1,ibox) .lt.
331+
! & box_MISR_ztop(1,ibox+1) ) then
332+
!
333+
! box_MISR_ztop(1,ibox) =
334+
! & box_MISR_ztop(1,ibox+1)
335+
! endif
336+
!
337+
! endif
338+
! enddo
339+
!
340+
! endif
338341

339-
!
342+
!
340343
! DETERMINE CLOUD TYPE FREQUENCIES
341344
!
342345
! Now that ztop and tau have been determined,

‎MODIS_simulator/modis_simulator.F90‎

Lines changed: 42 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -501,34 +501,47 @@ subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thic
501501
! ########################################################################################
502502
! Compute mean optical thickness.
503503
! ########################################################################################
504-
Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
505-
Cloud_Fraction_Total_Mean(1:nPoints)
506-
Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
507-
Cloud_Fraction_Water_Mean(1:nPoints)
508-
Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
509-
Cloud_Fraction_Ice_Mean(1:nPoints)
510-
511-
! ########################################################################################
512-
! We take the absolute value of optical thickness here to satisfy compilers that complains
513-
! when we evaluate the logarithm of a negative number, even though it's not included in
514-
! the sum.
515-
! ########################################################################################
516-
Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
517-
dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
518-
Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
519-
dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
520-
Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
521-
dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
522-
Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
523-
Cloud_Fraction_Water_Mean(1:nPoints)
524-
Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
525-
Cloud_Fraction_Ice_Mean(1:nPoints)
526-
Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
527-
max(1, count(cloudMask, dim = 2))
528-
Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
529-
mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
530-
Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
531-
mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
504+
where(Cloud_Fraction_Total_Mean(1:nPoints) > 0)
505+
Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
506+
Cloud_Fraction_Total_Mean(1:nPoints)
507+
Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
508+
dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
509+
elsewhere
510+
Optical_Thickness_Total_Mean = 0.
511+
Optical_Thickness_Total_MeanLog10 = 0.
512+
endwhere
513+
where(Cloud_Fraction_Water_Mean(1:nPoints) > 0)
514+
Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
515+
Cloud_Fraction_Water_Mean(1:nPoints)
516+
Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
517+
mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
518+
Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
519+
dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
520+
Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
521+
Cloud_Fraction_Water_Mean(1:nPoints)
522+
elsewhere
523+
Optical_Thickness_Water_Mean = 0.
524+
Optical_Thickness_Water_MeanLog10 = 0.
525+
Cloud_Particle_Size_Water_Mean = 0.
526+
Liquid_Water_Path_Mean = 0.
527+
endwhere
528+
where(Cloud_Fraction_Ice_Mean(1:nPoints) > 0)
529+
Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
530+
Cloud_Fraction_Ice_Mean(1:nPoints)
531+
Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
532+
mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
533+
Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
534+
dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
535+
Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
536+
Cloud_Fraction_Ice_Mean(1:nPoints)
537+
elsewhere
538+
Optical_Thickness_Ice_Mean = 0.
539+
Optical_Thickness_Ice_MeanLog10 = 0.
540+
Cloud_Particle_Size_Ice_Mean = 0.
541+
Ice_Water_Path_Mean = 0.
542+
endwhere
543+
Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
544+
max(1, count(cloudMask, dim = 2))
532545

533546
! ########################################################################################
534547
! Normalize pixel counts to fraction.
@@ -560,10 +573,7 @@ subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thic
560573
Cloud_Particle_Size_Ice_Mean = R_UNDEF
561574
Ice_Water_Path_Mean = R_UNDEF
562575
endwhere
563-
where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF
564-
where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF
565-
where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF
566-
576+
567577
! ########################################################################################
568578
! Joint histogram
569579
! ########################################################################################

‎actsim/lidar_simulator.F90‎

Lines changed: 24 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -199,7 +199,8 @@ SUBROUTINE lidar_simulator(npoints,nlev,npart,nrefl &
199199
REAL pnorm_ice(npoints,nlev) ! lidar backscattered signal power for ice
200200
REAL pnorm_perp_ice(npoints,nlev) ! perpendicular lidar backscattered signal power for ice
201201
REAL pnorm_perp_liq(npoints,nlev) ! perpendicular lidar backscattered signal power for liq
202-
202+
real epsreal
203+
203204
! Output variable
204205
REAL pnorm_perp_tot (npoints,nlev) ! perpendicular lidar backscattered signal power
205206

@@ -498,20 +499,24 @@ SUBROUTINE lidar_simulator(npoints,nlev,npart,nrefl &
498499
ENDDO
499500

500501
! Computation of beta_perp_ice/liq using the lidar equation
502+
epsreal = epsilon(1.)
501503
! Ice only
502504
! Upper layer
503505
beta_perp_ice(:,nlev) = pnorm_perp_ice(:,nlev) * (2.*tautot_ice(:,nlev)) &
504506
& / (1.-exp(-2.0*tautot_ice(:,nlev)))
505507

506508
DO k= nlev-1, 1, -1
507-
tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
508-
WHERE (tautot_lay_ice(:).GT.0.)
509-
beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
510-
& / (1.-exp(-2.0*tautot_lay_ice(:)))
511-
512-
ELSEWHERE
513-
beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/EXP(-2.0*tautot_ice(:,k+1))
514-
END WHERE
509+
tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
510+
WHERE ( EXP(-2.0*tautot_ice(:,k+1)) .gt. epsreal )
511+
WHERE (tautot_lay_ice(:).GT.0.)
512+
beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
513+
& / (1.-exp(-2.0*tautot_lay_ice(:)))
514+
ELSEWHERE
515+
beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/EXP(-2.0*tautot_ice(:,k+1))
516+
END WHERE
517+
elsewhere
518+
beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/epsreal
519+
endwhere
515520
ENDDO
516521

517522
! Liquid only
@@ -521,13 +526,16 @@ SUBROUTINE lidar_simulator(npoints,nlev,npart,nrefl &
521526

522527
DO k= nlev-1, 1, -1
523528
tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1)
524-
WHERE (tautot_lay_liq(:).GT.0.)
525-
beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ EXP(-2.0*tautot_liq(:,k+1)) * (2.*tautot_lay_liq(:)) &
526-
& / (1.-exp(-2.0*tautot_lay_liq(:)))
527-
528-
ELSEWHERE
529-
beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/EXP(-2.0*tautot_liq(:,k+1))
530-
END WHERE
529+
WHERE ( EXP(-2.0*tautot_liq(:,k+1)) .gt. epsreal )
530+
WHERE (tautot_lay_liq(:).GT.0.)
531+
beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ EXP(-2.0*tautot_liq(:,k+1)) * (2.*tautot_lay_liq(:)) &
532+
& / (1.-exp(-2.0*tautot_lay_liq(:)))
533+
ELSEWHERE
534+
beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/EXP(-2.0*tautot_liq(:,k+1))
535+
END WHERE
536+
elsewhere
537+
beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/epsreal
538+
endwhere
531539
ENDDO
532540

533541

0 commit comments

Comments
 (0)