- Timestamp:
- Nov 30, 2016, 1:28:41 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2665-2668,2670-2674,2677-2681,2683-2684,2686,2690-2719
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/cosp/modis_simulator.F90
r2435 r2720 2 2 ! Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences 3 3 ! All rights reserved. 4 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov.2013) $4 ! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $ 5 5 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $ 6 6 ! … … 79 79 real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90. 80 80 integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate 81 logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? 82 81 ! DJS2015: Remove unused parameter 82 ! logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? 83 ! DJS2015 END 83 84 ! 84 85 ! Precompute near-IR optical params vs size for retrieval scheme … … 125 126 nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + & 126 127 nominalPressureHistogramBoundaries(2, :) ) / 2. 128 ! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is 129 ! identical to what is done in COSPv.2.0.0 for histogram bin initialization. 130 integer :: j 131 integer,parameter :: & 132 numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram 133 numMODISReffIceBins = 6 ! Number of bins for tau/ReffICE joint-histogram 134 real,parameter,dimension(numMODISReffLiqBins+1) :: & 135 reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/) 136 real,parameter,dimension(numMODISReffIceBins+1) :: & 137 reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/) 138 real,parameter,dimension(2,numMODISReffIceBins) :: & 139 reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k), & 140 l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/), & 141 shape = (/2,numMODISReffIceBins/)) 142 real,parameter,dimension(2,numMODISReffLiqBins) :: & 143 reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k), & 144 l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/), & 145 shape = (/2,numMODISReffLiqBins/)) 146 real,parameter,dimension(numMODISReffIceBins) :: & 147 reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2. 148 real,parameter,dimension(numMODISReffLiqBins) :: & 149 reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2. 150 ! DJS2015 END 151 127 152 ! ------------------------------ 128 153 ! There are two ways to call the MODIS simulator: … … 384 409 385 410 end subroutine modis_L2_simulator_oneTau 411 412 ! ######################################################################################## 413 subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, & 414 Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & 415 Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & 416 Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & 417 Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,& 418 Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, & 419 Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & 420 Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq) 421 422 ! INPUTS 423 integer,intent(in) :: & 424 nPoints, & ! Number of horizontal gridpoints 425 nSubCols ! Number of subcolumns 426 integer,intent(in), dimension(:,:) :: & 427 !ds integer,intent(in), dimension(nPoints, nSubCols) :: & 428 phase 429 real,intent(in),dimension(:,:) :: & 430 !ds real,intent(in),dimension(nPoints, nSubCols) :: & 431 cloud_top_pressure, & 432 optical_thickness, & 433 particle_size 434 435 ! OUTPUTS 436 real,intent(inout),dimension(:) :: & ! 437 !ds real,intent(inout),dimension(nPoints) :: & ! 438 Cloud_Fraction_Total_Mean, & ! 439 Cloud_Fraction_Water_Mean, & ! 440 Cloud_Fraction_Ice_Mean, & ! 441 Cloud_Fraction_High_Mean, & ! 442 Cloud_Fraction_Mid_Mean, & ! 443 Cloud_Fraction_Low_Mean, & ! 444 Optical_Thickness_Total_Mean, & ! 445 Optical_Thickness_Water_Mean, & ! 446 Optical_Thickness_Ice_Mean, & ! 447 Optical_Thickness_Total_MeanLog10, & ! 448 Optical_Thickness_Water_MeanLog10, & ! 449 Optical_Thickness_Ice_MeanLog10, & ! 450 Cloud_Particle_Size_Water_Mean, & ! 451 Cloud_Particle_Size_Ice_Mean, & ! 452 Cloud_Top_Pressure_Total_Mean, & ! 453 Liquid_Water_Path_Mean, & ! 454 Ice_Water_Path_Mean ! 455 real,intent(inout),dimension(:,:,:) :: & 456 !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: & 457 Optical_Thickness_vs_Cloud_Top_Pressure 458 real,intent(inout),dimension(:,:,:) :: & 459 !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: & 460 Optical_Thickness_vs_ReffIce 461 real,intent(inout),dimension(:,:,:) :: & 462 !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: & 463 Optical_Thickness_vs_ReffLiq 464 465 ! LOCAL VARIABLES 466 real, parameter :: & 467 LWP_conversion = 2./3. * 1000. ! MKS units 468 integer :: i, j 469 logical, dimension(nPoints,nSubCols) :: & 470 cloudMask, & 471 waterCloudMask, & 472 iceCloudMask, & 473 validRetrievalMask 474 real,dimension(nPoints,nSubCols) :: & 475 tauWRK,ctpWRK,reffIceWRK,reffLiqWRK 476 477 ! ######################################################################################## 478 ! Include only those pixels with successful retrievals in the statistics 479 ! ######################################################################################## 480 validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0. 481 cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. & 482 validRetrievalMask(1:nPoints,1:nSubCols) 483 waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. & 484 validRetrievalMask(1:nPoints,1:nSubCols) 485 iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. & 486 validRetrievalMask(1:nPoints,1:nSubCols) 487 488 ! ######################################################################################## 489 ! Use these as pixel counts at first 490 ! ######################################################################################## 491 Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2)) 492 Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2)) 493 Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2)) 494 Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= & 495 highCloudPressureLimit, dim = 2)) 496 Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > & 497 lowCloudPressureLimit, dim = 2)) 498 Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)& 499 - Cloud_Fraction_Low_Mean(1:nPoints) 500 501 ! ######################################################################################## 502 ! Compute mean optical thickness. 503 ! ######################################################################################## 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) 532 533 ! ######################################################################################## 534 ! Normalize pixel counts to fraction. The first three cloud fractions have been set to -1 535 ! in cloud-free areas, so set those places to 0. 536 ! ######################################################################################## 537 Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols 538 Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols 539 Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols 540 541 ! ######################################################################################## 542 ! Set clear-scenes to undefined 543 ! ######################################################################################## 544 where (Cloud_Fraction_Total_Mean == 0) 545 Optical_Thickness_Total_Mean = R_UNDEF 546 Optical_Thickness_Total_MeanLog10 = R_UNDEF 547 Cloud_Top_Pressure_Total_Mean = R_UNDEF 548 endwhere 549 where (Cloud_Fraction_Water_Mean == 0) 550 Optical_Thickness_Water_Mean = R_UNDEF 551 Optical_Thickness_Water_MeanLog10 = R_UNDEF 552 Cloud_Particle_Size_Water_Mean = R_UNDEF 553 Liquid_Water_Path_Mean = R_UNDEF 554 endwhere 555 where (Cloud_Fraction_Ice_Mean == 0) 556 Optical_Thickness_Ice_Mean = R_UNDEF 557 Optical_Thickness_Ice_MeanLog10 = R_UNDEF 558 Cloud_Particle_Size_Ice_Mean = R_UNDEF 559 Ice_Water_Path_Mean = R_UNDEF 560 endwhere 561 where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF 562 where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF 563 where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF 564 565 ! ######################################################################################## 566 ! Joint histogram 567 ! ######################################################################################## 568 569 ! Loop over all points 570 tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols) 571 ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols) 572 reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask) 573 reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask) 574 do j=1,nPoints 575 576 ! Fill clear and optically thin subcolumns with fill 577 where(.not. cloudMask(j,1:nSubCols)) 578 tauWRK(j,1:nSubCols) = -999. 579 ctpWRK(j,1:nSubCols) = -999. 580 endwhere 581 ! Joint histogram of tau/CTP 582 call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,& 583 tauHistogramBoundaries,numTauHistogramBins,& 584 pressureHistogramBoundaries,numPressureHistogramBins,& 585 Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins)) 586 ! Joint histogram of tau/ReffICE 587 call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, & 588 tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds, & 589 numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins)) 590 ! Joint histogram of tau/ReffLIQ 591 call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, & 592 tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds, & 593 numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins)) 594 595 enddo 596 Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = & 597 Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols 598 Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = & 599 Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols 600 Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = & 601 Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols 602 603 end subroutine modis_column 604 ! ###################################################################################### 605 ! SUBROUTINE hist2D 606 ! ###################################################################################### 607 subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist) 608 implicit none 609 610 ! INPUTS 611 integer, intent(in) :: & 612 npts, & ! Number of data points to be sorted 613 nbin1, & ! Number of bins in histogram direction 1 614 nbin2 ! Number of bins in histogram direction 2 615 real,intent(in),dimension(npts) :: & 616 var1, & ! Variable 1 to be sorted into bins 617 var2 ! variable 2 to be sorted into bins 618 real,intent(in),dimension(nbin1+1) :: & 619 bin1 ! Histogram bin 1 boundaries 620 real,intent(in),dimension(nbin2+1) :: & 621 bin2 ! Histogram bin 2 boundaries 622 ! OUTPUTS 623 real,intent(out),dimension(nbin1,nbin2) :: & 624 jointHist 625 626 ! LOCAL VARIABLES 627 integer :: ij,ik 628 629 do ij=2,nbin1+1 630 do ik=2,nbin2+1 631 jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. & 632 var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik)) 633 enddo 634 enddo 635 end subroutine hist2D 636 386 637 !------------------------------------------------------------------------------------------------ 387 638 subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, & … … 666 917 ! If first retrieval works, can try 2nd iteration using greater re resolution 667 918 ! 668 if(use_two_re_iterations .and. retrieve_re > 0.) then 669 re_min = retrieve_re - delta_re 670 re_max = retrieve_re + delta_re 671 delta_re = (re_max - re_min)/real(num_trial_res-1) 672 673 trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) 674 g(:) = get_g_nir( phase, trial_re(:)) 675 w0(:) = get_ssa_nir(phase, trial_re(:)) 676 predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) 677 retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 678 end if 919 ! DJS2015: Remove unused piece of code 920 ! if(use_two_re_iterations .and. retrieve_re > 0.) then 921 ! re_min = retrieve_re - delta_re 922 ! re_max = retrieve_re + delta_re 923 ! delta_re = (re_max - re_min)/real(num_trial_res-1) 924 ! 925 ! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) 926 ! g(:) = get_g_nir( phase, trial_re(:)) 927 ! w0(:) = get_ssa_nir(phase, trial_re(:)) 928 ! predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) 929 ! retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 930 ! end if 931 ! DJS2015 END 679 932 else 680 933 retrieve_re = re_fill … … 739 992 real, intent(in) :: re 740 993 real :: get_g_nir 741 742 real, dimension(3), parameter :: ice_coefficients = (/ 0.7432, 4.5563e-3, -2.8697e-5 /), &743 small_water_coefficients = (/ 0.8027, -1.0496e-2, 1.7071e-3 /), &744 big_water_coefficients = (/ 0.7931, 5.3087e-3, -7.4995e-5 /)745 746 ! approx. fits from MODIS Collection 5 LUT scattering calculations747 if(phase == phaseIsLiquid) then 748 if(re < 8.) then749 get_g_nir = fit_to_quadratic(re, small_water_coefficients)750 if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)751 else752 get_g_nir = fit_to_quadratic(re,big_water_coefficients)753 if(re > re_water_max) get_g_nir = fit_to_quadratic(re_water_max, big_water_coefficients)754 end if994 995 real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), & 996 small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /) 997 real, dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /) 998 999 ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals 1000 if(phase == phaseIsLiquid) then 1001 if(re < 7.) then 1002 get_g_nir = fit_to_quadratic(re, small_water_coefficients) 1003 if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients) 1004 else 1005 get_g_nir = fit_to_cubic(re, big_water_coefficients) 1006 if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients) 1007 end if 755 1008 else 756 get_g_nir = fit_to_quadratic(re, ice_coefficients)1009 get_g_nir = fit_to_quadratic(re, ice_coefficients) 757 1010 if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients) 758 1011 if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients) … … 771 1024 ! Fits from Steve Platnick 772 1025 ! 1026 real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/) 1027 real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /) 773 1028 774 real, dimension(4), parameter :: ice_coefficients = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /) 775 real, dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /) 776 777 ! approx. fits from MODIS Collection 5 LUT scattering calculations 1029 ! approx. fits from MODIS Collection 6 LUT scattering calculations 778 1030 if(phase == phaseIsLiquid) then 779 1031 get_ssa_nir = fit_to_quadratic(re, water_coefficients)
Note: See TracChangeset
for help on using the changeset viewer.