!******************************************************************************* !******************************************************************************* ! THIS FILE CONTAINS THE FOLLOWING SUBROUTINES AND FUNCTIONS IN THE ! ORDER LISTED: ! ! RADIANT ! RAD ! ! BASIC_PROPERTIES ! BUILD_LAYER ! COMBINE_GS ! COMBINE_GT_GR ! COMBINE_LAYERS ! COMBINE_LAYERS_2 ! DATA_INSPECTOR ! DPLKAVG ! GET_GLOBAL_SOURCES5 ! GET_GT_GR4 ! INTERMEDIATE_RESULTS ! PLANCK ! SPECTRAL_PERTURBATION ! SURF_REF1_3 ! SURF_REF2_2 ! !******************************************************************************* !******************************************************************************* MODULE RADIANT2S IMPLICIT NONE !PRIVATE DATA PRIVATE :: & MU,W,YPLEGP,YPLEGM DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE :: & MU,W DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: & YPLEGP,YPLEGM CONTAINS !******************************************************************************* !******************************************************************************* SUBROUTINE RADIANT(N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,SRC_FLAG,FSUN,SZA,& ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,PLANCK_FLAG,WVN,WVNLO,WVNHI,& PHI,QUAD,DG,TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,& OMEGA_CLD,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& PF_CLD,CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO,& SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS,& AZIMUTHAL_RAD_ONLY,FOURIER_TOL,FLUXES,WARNING_FLAG,& P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO,ITM,IBMTOT,IBPTOT,ITPTOT,& FIT_TOT,FOT_TOT,FOB_TOT,FIB_TOT,MU) !INPUT : N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,SRC_FLAG,FSUN,SZA, ! ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,PLANCK_FLAG,WVN, ! WVNLO,WVNHI,PHI,QUAD,DG, ! TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD, ! G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,PF_CLD, ! CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO, ! SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS, ! FLUXES,WARNING_FLAG,P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO !OUTPUT: ITM,IBMTOT,IBPTOT,ITPTOT,FIT_TOT,FOT_TOT,FOB_TOT,FIB_TOT !THIS PROGRAM SERVES AS A PREPROCESSOR FOR SUBROUTINE RAD !PROGRAMMER: MATT CHRISTI !DATE: 10/01/03 !DATA DICTIONARY**************************************************************** ! ! ALBEDO = SURFACE ALBEDO WHEN SURF=1 ! B_IO = SUBROUTINE BUILD_LAYER I/O FLAG ! (0=OFF, 1=ON) ! B_T = VECTOR OF SPECTRAL RADIANCES OR RADIANCES (DEPENDING ON VALUE ! OF PLANCK_FLAG) USING THE PLANCK FUNCTION OF EMISSION ! BSURF = SURFACE RADIANCE (OR SURFACE SPECTRAL RADIANCE) DERIVED FROM ! THE PLANCK FUNCTION OF EMISSION ! CHI_AERO = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_CLD = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_SURF = VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED) ! DELTA_M_AERO = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE AEROSOL PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DELTA_M_CLD = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE CLOUD PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS, 1=DOUBLE GAUSS) ! EPS = AMOUNT BY WHICH OPTICAL PROPERTIES ARE PERTURBED BY FOR THE ! COMPUTATION OF JACOBIAN ELEMENTS BY FINITE DIFFERENCE WHEN ! THE PERTURBATION SCHEME IS IN USE ! FIB_TOT = TOTAL FLUX INTO THE BOTTOM BOUNDARY ! FIT_TOT = TOTAL FLUX INTO THE TOP BOUNDARY ! FLUX_IO = FLUX TEST I/O FLAG FOR SUBROUTINE RAD ! (0=OFF, 1=ON) ! FLUXES = RADIANT FLUX OUTPUT I/O FLAG ! (0=OFF, 1=ON) ! FOB_TOT = TOTAL FLUX OUT OF THE BOTTOM BOUNDARY ! FOT_TOT = TOTAL FLUX OUT OF THE TOP BOUNDARY ! FSUN = FLUX AT THE TOP OF THE ATMOSPHERE AT A GIVEN WAVELENGTH ! G_AERO = VECTOR OF AEROSOL EFFECTIVE ASYMMETRY PARAMETERS ! G1_AERO = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR FORWARD SCATTERING ! G2_AERO = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING ! G_CLD = VECTOR OF CLOUD EFFECTIVE ASYMMETRY PARAMETERS ! G1_CLD = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR FORWARD SCATTERING ! G2_CLD = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING ! INV_IO = SUBROUTINE INVERT I/O FLAG (FOR MATRIX INVERSE) ! (0=OFF, 1=ON) ! IBMTOT = VECTOR OF TOTAL DOWNWELLING RADIANCES AT THE BOTTOM OF THE ! ATMOSPHERE FROM ALL FOURIER COMPONENTS COMPUTED ! IBPTOT = VECTOR OF TOTAL UPWELLING RADIANCE AT THE BOTTOM OF THE ! ATMOSPHERE FROM ALL FOURIER COMPONENTS COMPUTED ! INPUT_IO = I/O FLAG TO PERMIT DISPLAY OF USER INPUT ! (0=OFF, 1=ON) ! ITM = VECTOR OF TOTAL DOWNWELLING DIFFUSE RADIANCES FROM THE TOP ! BOUNDARY ! ITMS = VECTOR OF DOWNWELLING DIFFUSE SOLAR RADIANCES FROM THE TOP ! BOUNDARY ! ITMT = VECTOR OF DOWNWELLING DIFFUSE THERMAL RADIANCES FROM THE TOP ! BOUNDARY ! ITPTOT = VECTOR OF TOTAL UPWELLING RADIANCES AT THE TOP BOUNDARY ! FROM ALL FOURIER COMPONENTS COMPUTED ! KFLAG = FLAG TO INDICATE WHETHER OR NOT ADDITIONAL COMPUTATIONS WILL ! BE PERFORMED INSIDE RADIANT TO PREPARE FOR SUBSEQUENT CALLS ! TO RADIANT WHERE SAVED MATRICES WILL BE USED ! (0=NO, 1=YES) ! NOTE: KFLAG USED IN CONJUNCTION WITH WVN_FLAG ! L_IO = SUBROUTINE LOCAL I/O FLAG ! (0=OFF, 1=ON) ! LAYCHG = LAYER NUMBER OF CHANGED LAYER WHEN IN LAYER-SAVING MODE ! MU0 = COSINE OF THE SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! OMEGA_AERO = VECTOR OF SINGLE SCATTER ALBEDO FOR AEROSOL ! OMEGA_CLD = VECTOR OF SINGLE SCATTER ALBEDO FOR CLOUD ! OUTPUT_IO = I/O FLAG TO PERMIT DISPLAY OF OUTPUT RADIANCES TO SCREEN ! (0=OFF, 1=ON) ! P_IO = SUBROUTINE RADIANT I/O FLAG ! (0=OFF, 1=ON) ! PERT_FLAG = FLAG TO INDICATE WHETHER OR NOT THE PERTURBATION SCHEME WILL ! BE USED TO BUILD LAYERS ON SUBSEQUENT CALLS TO RADIANT DURING ! LAYER SAVING MODE ! (0=NO, 1=YES) ! PF_AERO = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_AERO) ! PF_CLD = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_CLD) ! PHI = AZIMUTH ANGLE FOR WHICH THE RADIANCE VECTOR IS COMPUTED ! (in degrees) ! NOTE: REFERENCE AZIMUTH ASSUMED TO BE 0 ! PLANCK_FLAG = FLAG TO CHOOSE WHETHER SPECTRAL RADIANCES OR RADIANCES WILL ! BE USED FOR THERMAL SOURCES ! (0=SPECTRAL RADIANCES (in mW/(m^2 cm^-1 ster)) ! 1=RADIANCES (in mW/(m^2 ster)) ) ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED ! (0=GAUSS, 1=LOBATTO) ! R_IO = SUBROUTINE RAD I/O FLAG ! (0=OFF, 1=ON) ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES ! (1=SOLAR ONLY, 2=SOLAR & THERMAL, 3=THERMAL ONLY) ! SURF = FLAG INDICATING TYPE OF REFLECTING SURFACE ! (0=NONE, 1=LAMBERTIAN, 2=DETERMINED BY CHI_SURF) ! SZA = THE SOLAR ZENITH ANGLE (in degrees from zenith) ! TAU_AERO = VECTOR OF OPTICAL DEPTHS FOR AEROSOL ! TAU_CLD = VECTOR OF OPTICAL DEPTHS FOR CLOUD ! TAU_GAS = VECTOR OF OPTICAL DEPTHS FOR ATMOSPHERIC GAS ! TAU_RAY = VECTOR OF OPTICAL DEPTHS FOR RAYLEIGH SCATTERING ! TEMISS = EMISSIVITY OF TOP BOUNDARY FOR THERMAL SOURCES ! TLEV = VECTOR (I.E. VERTICAL PROFILE) OF LEVEL-BASED ATMOSPHERIC ! TEMPERATURES (in K) ! TSURF = TEMPERATURE OF THE SURFACE (in K) ! TTOP = TEMPERATURE OF THE TOP BOUNDARY (in K) ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE ! (0=OFF, 1=ON) ! WVN = WAVENUMBER FOR WHICH SPECTRAL RADIANCE IS COMPUTED FOR PLANCK ! FUNCTION (in cm^-1) ! WVNLO = LOW WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS ! COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1) ! WVNHI = HIGH WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS ! COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1) ! WVN_FLAG = FLAG TO INDICATE WHETHER OR NOT SAVED MATRICES WILL BE ! USED TO COMPUTE RADIANCES ON CURRENT CALL TO SUBROUTINE ! RADIANT ! ('NEW'=NO, 'OLD'=YES) ! NOTE: WVN_FLAG USED IN CONJUNCTION WITH KFLAG ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG ! (0=OFF, 1=ON) ! !******************************************************************************* !INTRINSIC SUBPROGRAMS USED BY RADIANT****************************************** ! NONE !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY RADIANT******************************************* ! DATA_INSPECTOR,DPLKAVG,GETQUAD2,PLANCK,RAD !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & B_IO,DG,FLUX_IO,FLUXES,INV_IO,KFLAG,L_IO,N,& NUMLAY,P_IO,PLANCK_FLAG,QUAD,R_IO,SRC_FLAG,SURF,VIEW_FLAG,& WARNING_FLAG,PERT_FLAG,LAYCHG INTEGER, DIMENSION(NUMLAY) :: & DELTA_M_AERO,DELTA_M_CLD,PF_AERO,PF_CLD DOUBLE PRECISION :: & ALBEDO,FSUN,FOURIER_TOL,PHI,SZA,TEMISS,TSURF,& TTOP,VIEW_ANGLE,WVN,WVNLO,WVNHI,EPS DOUBLE PRECISION, DIMENSION(N) :: & ITMS,ITMT DOUBLE PRECISION, DIMENSION(NUMLAY) :: & G_AERO,G_CLD,G1_AERO,G1_CLD,G2_AERO,G2_CLD,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD DOUBLE PRECISION, DIMENSION(NUMLAY+1) :: & TLEV DOUBLE PRECISION, DIMENSION(0:(4*N)) :: & CHI_SURF DOUBLE PRECISION, DIMENSION(0:(4*N),NUMLAY) :: & CHI_AERO,CHI_CLD CHARACTER (LEN=3) :: & WVN_FLAG LOGICAL :: & AZIMUTHAL_RAD_ONLY !OUTPUT VARIABLES DOUBLE PRECISION :: & FIT_TOT,FOB_TOT,FIB_TOT,FOT_TOT DOUBLE PRECISION, DIMENSION(N) :: & ITM,IBMTOT,IBPTOT,ITPTOT !INTERNAL VARIABLES INTEGER :: & I,NUMTEMP DOUBLE PRECISION :: & BSURF,BTOP,MU0 DOUBLE PRECISION, DIMENSION(1) :: & DUMMY_IN,DUMMY_OUT DOUBLE PRECISION, DIMENSION(N) :: & MU,W DOUBLE PRECISION, DIMENSION(NUMLAY+1) :: & B_T INTEGER, SAVE :: & DG_SAVE,N_SAVE,NUMLAY_SAVE,& QUAD_SAVE,VIEW_FLAG_SAVE DOUBLE PRECISION, SAVE :: & MU0_PERT=-1.0D0,SZA_SAVE,VIEW_ANGLE_SAVE LOGICAL, SAVE :: & AZIMUTHAL_RAD_ONLY_SAVE,& FIRST_RUN = .TRUE.,& QUAD_4VAR_RESET = .FALSE.,& SZA_RESET = .FALSE.,& QUAD_RESET = .FALSE.,& DIMENSIONS_RESET = .FALSE. DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 DOUBLE PRECISION F_CLD !START PROGRAM IF (P_IO == 1) THEN PRINT* PRINT*,'ENTERING RADIANT' END IF !INSPECT INPUT DATA FOR SUITABILITY CALL DATA_INSPECTOR(SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV,& TSURF,PLANCK_FLAG,WVN,WVNLO,WVNHI,VIEW_FLAG,VIEW_ANGLE,PHI,& QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO,WVN_FLAG,KFLAG,& NUMLAY,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,& PF_CLD,CHI_CLD,PF_AERO,CHI_AERO,ALBEDO,SURF,CHI_SURF,& FOURIER_TOL,FLUXES,WARNING_FLAG,P_IO,R_IO,FLUX_IO,L_IO,& B_IO,INV_IO) !COMPUTE COSINE OF SZA MU0 = DCOS(SZA/180.0D0*PI) !SAVE CERTAIN VARIABLES ON FIRST CALL TO RADIANT AND CHECK THEM ON LATER !CALLS IF (FIRST_RUN) THEN FIRST_RUN = .FALSE. IF ( WVN_FLAG /= 'NEW' ) THEN PRINT* PRINT*,'ERROR: WVN_FLAG SHOULD = ''NEW'' THE FIRST TIME ' // & 'RADIANT IS CALLED' STOP END IF !SAVE CERTAIN VARIABLES TO CHECK ON FUTURE CALLS SZA_SAVE = SZA QUAD_SAVE = QUAD DG_SAVE = DG VIEW_FLAG_SAVE = VIEW_FLAG VIEW_ANGLE_SAVE = VIEW_ANGLE N_SAVE = N NUMLAY_SAVE = NUMLAY AZIMUTHAL_RAD_ONLY_SAVE = AZIMUTHAL_RAD_ONLY !SINGULARITY BUSTER FOR MU0 CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) ! PRINT*,MU, 'MU' DO I=1,N IF (DABS(MU0-MU(I)) < 1.0D-11) THEN !IF MU0 TOO CLOSE TO ONE OF THE MU(I), THEN PERTURB IT IF((MU0 - 1.0D-10) < 0.0D0) THEN MU0 = MU0 + 1.0D-10 ELSE MU0 = MU0 - 1.0D-10 END IF MU0_PERT = MU0 PRINT* PRINT*,'SZA PERTURBED TO AVOID SINGULARITY. NOW, SZA = ', & (DACOS(MU0_PERT)/PI)*180.0D0 END IF END DO ELSE !SZA_RESET CHECK IF ( (SZA < (SZA_SAVE-1.0D-12)) .OR. & (SZA > (SZA_SAVE+1.0D-12)) ) THEN PRINT* PRINT*,'SZA HAS CHANGED' SZA_SAVE = SZA SZA_RESET = .TRUE. END IF !QUAD_RESET CHECK IF ( (QUAD /= QUAD_SAVE) .OR. (DG /= DG_SAVE) .OR. & (N /= N_SAVE) .OR. (VIEW_FLAG /= VIEW_FLAG_SAVE) .OR. & (VIEW_ANGLE < (VIEW_ANGLE_SAVE-1.0D-12)) .OR. & (VIEW_ANGLE > (VIEW_ANGLE_SAVE+1.0D-12)) ) THEN PRINT* IF ( (QUAD /= QUAD_SAVE) .OR. (DG /= DG_SAVE) .OR. & (N /= N_SAVE) .OR. (VIEW_FLAG /= VIEW_FLAG_SAVE) ) THEN QUAD_4VAR_RESET = .TRUE. END IF IF (QUAD /= QUAD_SAVE) THEN PRINT*,'QUAD HAS CHANGED' QUAD_SAVE = QUAD END IF IF (DG /= DG_SAVE) THEN PRINT*,'DG HAS CHANGED' DG_SAVE = DG END IF IF (N /= N_SAVE) THEN PRINT*,'N HAS CHANGED' N_SAVE = N END IF IF (VIEW_FLAG /= VIEW_FLAG_SAVE) THEN PRINT*,'VIEW_FLAG HAS CHANGED' VIEW_FLAG_SAVE = VIEW_FLAG END IF IF ( (VIEW_ANGLE < (VIEW_ANGLE_SAVE-1.0D-12)) .OR. & (VIEW_ANGLE > (VIEW_ANGLE_SAVE+1.0D-12)) ) THEN PRINT*,'VIEW_ANGLE HAS CHANGED' VIEW_ANGLE_SAVE = VIEW_ANGLE END IF QUAD_RESET = .TRUE. END IF IF ( SZA_RESET .OR. QUAD_RESET ) THEN !SZA OR QUADRATURE PARAMETERS CHANGED BY THE USER. RECHECK !COMPATIBILITY OF CURRENT MU0 WITH CURRENT MU(I) MU0_PERT = -1.0D0 !SINGULARITY BUSTER FOR MU0 CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) PRINT*,MU, 'MU' DO I=1,N IF (DABS(MU0-MU(I)) < 1.0D-11) THEN !IF MU0 TOO CLOSE TO ONE OF THE MU(I), THEN PERTURB IT IF((MU0 - 1.0D-10) < 0.0D0) THEN MU0 = MU0 + 1.0D-10 ELSE MU0 = MU0 - 1.0D-10 END IF MU0_PERT = MU0 PRINT* PRINT*,'SZA PERTURBED TO AVOID SINGULARITY. NOW, SZA = ', & (DACOS(MU0_PERT)/PI)*180.0D0 END IF END DO ELSE IF (MU0_PERT > 0.0D0) THEN !CURRENT SZA PERTURBED ON A PRIOR CALL TO AVOID SINGULARITY. !RESET MU0 AGAIN TO PREVIOUSLY DETERMINED APPROPRIATE VALUE. MU0 = MU0_PERT END IF !DIMENSIONS_RESET CHECK IF ( QUAD_4VAR_RESET .OR. (NUMLAY /= NUMLAY_SAVE) .OR. & (AZIMUTHAL_RAD_ONLY .NEQV. AZIMUTHAL_RAD_ONLY_SAVE) ) THEN PRINT* !(IF BLOCKS FOR CHANGES IN QUAD, DG, N, AND VIEW_FLAG ! INSIDE QUAD_RESET CHECK) IF (NUMLAY /= NUMLAY_SAVE) THEN PRINT*,'NUMLAY HAS CHANGED' NUMLAY_SAVE = NUMLAY END IF IF (AZIMUTHAL_RAD_ONLY .NEQV. AZIMUTHAL_RAD_ONLY_SAVE) THEN PRINT*,'AZIMUTHAL_RAD_ONLY HAS CHANGED' AZIMUTHAL_RAD_ONLY_SAVE = AZIMUTHAL_RAD_ONLY END IF DIMENSIONS_RESET = .TRUE. END IF !print*,'DIMENSIONS_RESET=',DIMENSIONS_RESET !print*,'SZA_RESET=',SZA_RESET !print*,'QUAD_RESET=',QUAD_RESET IF ( (WVN_FLAG == 'OLD') .AND. & (DIMENSIONS_RESET .OR. SZA_RESET .OR. QUAD_RESET) ) THEN PRINT* PRINT*,'ERROR: WVN_FLAG = ''OLD'' AND LAYER-SAVING MODE ' // & 'ASSUMED; THEREFORE, N, NUMLAY, SZA, QUAD, ' // & 'DG, VIEW_FLAG, VIEW_ANGLE, AND AZIMUTHAL_RAD_ONLY ' // & 'MUST NOT BE CHANGED.' STOP END IF END IF !SINGULARITY BUSTER FOR CLOUD AND AEROSOL SINGLE SCATTER ALBEDO !(POSSIBLE RAYLEIGH SCATTERING PROBLEM TAKEN CARE OF IN SUBROUTINE BASIC !PROPERTIES) DO I=1,NUMLAY IF ((TAU_CLD(I) > 1.0D-12).AND.(TAU_AERO(I) > 1.0D-12).AND. & (TAU_GAS(I) <= 1.0D-12)) THEN !BOTH CLOUD AND AEROSOL PRESENT, BUT GAS NOT PRESENT IF ((OMEGA_CLD(I) + OMEGA_AERO(I)) > & 1.9999999999D0) THEN OMEGA_AERO(I) = OMEGA_AERO(I)*0.9999999999D0 PRINT* PRINT*,'OMEGA_AERO PERTURBED TO AVOID SINGULARITY.' PRINT*,'NOW, OMEGA_AERO = ',OMEGA_AERO(I),' FOR LAYER ',I END IF ELSE IF ((TAU_CLD(I) > 1.0D-12).AND.(TAU_AERO(I) <= 1.0D-12).AND. & (TAU_GAS(I) <= 1.0D-12)) THEN !CLOUD PRESENT, BUT AEROSOL AND GAS NOT PRESENT IF (OMEGA_CLD(I) > 0.9999999999D0) THEN OMEGA_CLD(I) = OMEGA_CLD(I)*0.9999999999D0 PRINT* PRINT*,'OMEGA_CLD PERTURBED TO AVOID SINGULARITY.' PRINT*,'NOW, OMEGA_CLD = ',OMEGA_CLD(I),' FOR LAYER ',I END IF ELSE IF ((TAU_CLD(I) <= 1.0D-12).AND.(TAU_AERO(I) > 1.0D-12).AND. & (TAU_GAS(I) <= 1.0D-12)) THEN !AEROSOL PRESENT, BUT CLOUD AND GAS NOT PRESENT IF (OMEGA_AERO(I) > 0.9999999999D0) THEN OMEGA_AERO(I) = OMEGA_AERO(I)*0.9999999999D0 PRINT* PRINT*,'OMEGA_AERO PERTURBED TO AVOID SINGULARITY.' PRINT*,'NOW, OMEGA_AERO = ',OMEGA_AERO(I),' FOR LAYER ',I END IF END IF END DO !COMPUTE QUANTITIES FOR ATMOSPHERIC AND SURFACE THERMAL EMISSION IF IN USE IF ((SRC_FLAG == 2).OR.(SRC_FLAG == 3)) THEN NUMTEMP = NUMLAY + 1 IF (PLANCK_FLAG == 1) THEN !COMPUTE VECTOR OF PLANCK SPECTRAL RADIANCES !AT DISCRETE WAVENUMBER IF (TTOP > 0.0D0) THEN DUMMY_IN(1) = TTOP CALL PLANCK(WVN,1,DUMMY_IN,DUMMY_OUT) BTOP = DUMMY_OUT(1) ITMT = TEMISS*BTOP END IF CALL PLANCK(WVN,NUMTEMP,TLEV,B_T) DUMMY_IN(1) = TSURF CALL PLANCK(WVN,1,DUMMY_IN,DUMMY_OUT) BSURF = DUMMY_OUT(1) ELSE IF (PLANCK_FLAG == 2) THEN !COMPUTE VECTOR OF PLANCK RADIANCES !OVER WAVENUMBER INTERVAL IF (TTOP > 0.0D0) THEN BTOP = DPLKAVG(WVNLO,WVNHI,TTOP) ITMT = TEMISS*BTOP END IF DO I=1,NUMTEMP B_T(I) = DPLKAVG(WVNLO,WVNHI,TLEV(I)) END DO BSURF = DPLKAVG(WVNLO,WVNHI,TSURF) END IF ELSE ITMT = 0.0D0 BSURF = 0.0D0 END IF ITM = ITMS + ITMT IF ((P_IO == 1).AND.((SRC_FLAG == 2).OR.(SRC_FLAG == 3))) THEN PRINT* WRITE(*,45) 'LEVEL','B_T(LEVEL)' DO I=1,NUMLAY+1 WRITE(*,10) I,B_T(I) END DO PRINT* PRINT*,'BSURF = ',BSURF 10 FORMAT(1X,I3,6X,E15.9E2) 45 FORMAT(2X,A5,6X,A10) END IF !COMPUTE RADIANCES ! IF (WVN_FLAG == 'NEW') THEN ! PRINT*,'SZA = ',SZA ! PRINT*,'SZA_SAVE = ',SZA_SAVE ! PRINT*,'QUAD = ',QUAD ! PRINT*,'QUAD_SAVE = ',QUAD_SAVE ! PRINT*,'DG = ',DG ! PRINT*,'DG_SAVE = ',DG_SAVE ! PRINT*,'VIEW_FLAG = ',VIEW_FLAG ! PRINT*,'VIEW_FLAG_SAVE = ',VIEW_FLAG_SAVE ! PRINT*,'VIEW_ANGLE = ',VIEW_ANGLE ! PRINT*,'VIEW_ANGLE_SAVE = ',VIEW_ANGLE_SAVE ! PRINT*,'N = ',N ! PRINT*,'N_SAVE = ',N_SAVE ! PRINT*,'NUMLAY = ',NUMLAY ! PRINT*,'NUMLAY_SAVE = ',NUMLAY_SAVE ! PRINT*,'AZIMUTHAL_RAD_ONLY = ',AZIMUTHAL_RAD_ONLY ! PRINT*,'AZIMUTHAL_RAD_ONLY_SAVE = ',AZIMUTHAL_RAD_ONLY_SAVE ! END IF !if (wvn_flag == 'NEW') then ! print*,'mu0=',mu0 !end if CALL RAD(FSUN,MU0,PHI,QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO,& NUMLAY,ALBEDO,SURF,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,& B_T,BSURF,FLUXES,WARNING_FLAG,R_IO,FLUX_IO,L_IO,B_IO,& INV_IO,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS,SRC_FLAG,VIEW_FLAG,& VIEW_ANGLE,PF_CLD,PF_AERO,CHI_CLD,CHI_AERO,CHI_SURF,& AZIMUTHAL_RAD_ONLY,& FOURIER_TOL,DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,& ITM,IBMTOT,IBPTOT,ITPTOT,FIT_TOT,FOT_TOT,FOB_TOT,FIB_TOT, F_CLD) QUAD_4VAR_RESET = .FALSE. SZA_RESET = .FALSE. QUAD_RESET = .FALSE. DIMENSIONS_RESET = .FALSE. IF (P_IO == 1) THEN PRINT* PRINT*,'LEAVING RADIANT' END IF END SUBROUTINE RADIANT !******************************************************************************* !******************************************************************************* SUBROUTINE RAD(FSUN,MU0,PHI,QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO,& NUMLAY,ALBEDO,SURF,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,& B_T,BSURF,FLUXES,WARNING_FLAG,R_IO,FLUX_IO,L_IO,B_IO,& INV_IO,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS,SRC_FLAG,VIEW_FLAG,& VIEW_ANGLE,PF_CLD,PF_AERO,CHI_CLD,CHI_AERO,CHI_SURF,& AZIMUTHAL_RAD_ONLY,& FOURIER_TOL,DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,& ITM,IBMTOT,IBPTOT,ITPTOT,FIT_TOT,FOT_TOT,FOB_TOT,FIB_TOT,F_CLD) !INPUT : FSUN,MU0,PHI,QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO, ! NUMLAY,ALBEDO,SURF,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO, ! TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD, ! B_T,BSURF,FLUXES,WARNING_FLAG,R_IO,FLUX_IO,L_IO,B_IO, ! INV_IO,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS,SRC_FLAG,VIEW_FLAG, ! VIEW_ANGLE,PF_CLD,PF_AERO,CHI_CLD,CHI_AERO,CHI_SURF,ITM !OUTPUT: IBMTOT,IBPTOT,ITPTOT,FIT_TOT,FOT_TOT,FOB_TOT,FIB_TOT !THIS PROGRAM BUILDS AN MULTI-LAYER PLANE PARALLEL ATMOSPHERE FROM !THE TOP DOWN. THIS VERSION HAS BEEN DESIGNED WITH CERTAIN SAVING FEATURES !TO SPEED UP THE COMPUTATION OF THE JACOBIAN WHILE PERFORMING A RETRIEVAL. !IT IS SET FOR 2N-STREAM (N UP AND N DOWN). !IT IS IN Z-SPACE. !PROGRAMMER: MATT CHRISTI !DATE: 9/27/04 !DATA DICTIONARY**************************************************************** ! ! ALBEDO = SURFACE ALBEDO ! B_IO = SUBROUTINE BUILD_LAYER I/O FLAG (0=OFF, 1=ON) ! B_T = VECTOR OF SPECTRAL RADIANCES OR RADIANCES (DEPENDING ON VALUE ! OF PLANCK_FLAG) USING THE PLANCK FUNCTION OF EMISSION ! BC_SRCS = SOURCE VECTOR ASSOCIATED WITH DIRECT SOLAR BEAM AND SURFACE ! EMISSION SOURCES ! BLC_FLAG = FLAG TO INDICATE IF THE BOTTOM LAYER HAS BEEN REBUILT ! (I.E. "BOTTOM LAYER CHANGED") IN ORDER TO COMPUTE THE ! RADIANCES FOR THE PERTURBED ATMOSPHERIC SCENE ! BSURF = SURFACE RADIANCE (OR SURFACE SPECTRAL RADIANCE) DERIVED FROM ! THE PLANCK FUNCTION OF EMISSION ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! CHI_AERO = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_CLD = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_SURF = VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED) ! DELTA_M_AERO = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE AEROSOL PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DELTA_M_CLD = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE CLOUD PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! E = IDENTITY MATRIX ! EMISS = VECTOR OF DIRECTIONAL EMISSIVITIES FOR SURFACE EMISSION ! EPS = AMOUNT BY WHICH OPTICAL PROPERTIES ARE PERTURBED BY FOR THE ! COMPUTATION OF JACOBIAN ELEMENTS BY FINITE DIFFERENCE WHEN ! THE PERTURBATION SCHEME IS IN USE ! FBOT = FLUX AT THE BOTTOM OF THE ATMOSPHERE AT A GIVEN WAVELENGTH ! FIB_TOT = TOTAL FLUX INTO THE BOTTOM BOUNDARY ! FIT_TOT = TOTAL FLUX INTO THE TOP BOUNDARY ! FLUX_IO = FLUX TEST PROGRAM SEGMENT I/O FLAG (0=OFF, 1=ON) ! FLUXES = RADIANT FLUX OUTPUT I/O FLAG (0=OFF, 1=ON) ! FOB_TOT = TOTAL FLUX OUT OF THE BOTTOM BOUNDARY ! FOT_TOT = TOTAL FLUX OUT OF THE TOP BOUNDARY ! FTOP = FLUX AT THE TOP OF A LAYER AT A GIVEN WAVELENGTH ! FSUN = FLUX AT THE TOP OF THE ATMOSPHERE AT A GIVEN WAVELENGTH ! G*_AERO = MATRIX OF AEROSOL ASYMMETRY PARAMETERS ! G*_CLD = MATRIX OF CLOUD ASYMMETRY PARAMETERS ! GRM = GLOBAL REFLECTION MATRIX MINUS (DOWNWELLING) ! GRP = GLOBAL REFLECTION MATRIX PLUS (UPWELLING) ! GSMS = GLOBAL SOURCE VECTOR MINUS SOLAR (DOWNWELLING) ! GSPS = GLOBAL SOURCE VECTOR PLUS SOLAR (UPWELLING) ! GSMT = GLOBAL SOURCE VECTOR MINUS THERMAL (DOWNWELLING) ! GSPT = GLOBAL SOURCE VECTOR PLUS THERMAL (UPWELLING) ! GTM = GLOBAL TRANSMISSION MATRIX MINUS (DOWNWELLING) ! GTP = GLOBAL TRANSMISSION MATRIX PLUS (UPWELLING) ! G**_TEMP = TEMPORARY STORAGE MATRICES ! IBM = DOWNWELLING RADIANCE AT THE BOTTOM OF THE ATMOSPHERE FOR A ! PARTICULAR FOURIER COMPONENT ! IBMTOT = TOTAL DOWNWELLING RADIANCE AT THE BOTTOM OF THE ATMOSPHERE ! FROM ALL FOURIER COMPONENTS COMPUTED ! IBP = UPWELLING RADIANCE AT THE BOTTOM OF THE ATMOSPHERE FOR A ! PARTICULAR FOURIER COMPONENT ! IBPTOT = TOTAL UPWELLING RADIANCE AT THE BOTTOM OF THE ATMOSPHERE ! FROM ALL FOURIER COMPONENTS COMPUTED ! INV_IO = SUBROUTINE INVERT I/O FLAG (0=OFF, 1=ON) (FOR MATRIX INVERSE) ! ITM = DOWNWELLING RADIANCE AT THE TOP OF THE ATMOSPHERE FOR A ! PARTICULAR FOURIER COMPONENT ! ITP = UPWELLING RADIANCE AT THE TOP OF THE ATMOSPHERE FOR A ! PARTICULAR FOURIER COMPONENT ! ITPTOT = TOTAL UPWELLING RADIANCE AT THE TOP OF THE ATMOSPHERE ! FROM ALL FOURIER COMPONENTS COMPUTED ! KFLAG = FLAG INDICATING WHETHER OR NOT LAYER SAVING FUNCTIONS ARE ! UTILIZED (E.G. WHEN COMPUTING ELEMENTS OF A JACOBIAN BY ! FINITE DIFFERENCE) (0=NO,1=YES) ! L_IO = SUBROUTINE LOCAL I/O FLAG (0=OFF, 1=ON) ! LAYCHG = LAYER NUMBER OF CHANGED LAYER WHEN IN LAYER-SAVING MODE ! LAYER = LAYER IN WHICH GLOBAL R, T, AND S MARTICES ARE BEING ! CALCULATED ! MAT* = WORK MATRICES ! MU = VECTOR HOLDING QUADRATURE ROOTS (COSINES OF QUADRATURE ! ANGLES) ! MU0 = COSINE OF THE SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTIONS (DEPENDENT ON NUMBER OF STREAMS AND TYPE OF ! QUADRATURE EMPLOYED) ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! OMEGA_AERO = VECTOR OF SINGLE SCATTER ALBEDO FOR AEROSOL ! OMEGA_CLD = VECTOR OF SINGLE SCATTER ALBEDO FOR CLOUD ! PERT_FLAG = FLAG TO INDICATE WHETHER OR NOT THE PERTURBATION SCHEME WILL ! BE USED TO BUILD LAYERS ON SUBSEQUENT CALLS TO RADIANT DURING ! LAYER SAVING MODE ! (0=NO, 1=YES) ! PF_AERO = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_AERO) ! PF_CLD = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_CLD) ! PHI = THE AZIMUTH ANGLE FOR WHICH THE RADIANCE VECTOR IS COMPUTED ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! SUBROUTINE LOCAL (0=GAUSS,1=LOBATTO) ! R_IO = RADIANCE SUBROUTINE I/O FLAG (0=OFF, 1=ON) ! RG = SURFACE REFLECTION MATRIX ("R GROUND") ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES (1=SOLAR ONLY,2=BOTH,3=THERMAL ONLY) ! SURF = FLAG INDICATING TYPE OF REFLECTING SURFACE ! TAU = VECTOR HOLDING LAYER OPTICAL DEPTHS ! TAU_AERO = VECTOR OF OPTICAL DEPTHS FOR AEROSOL ! TAU_CLD = VECTOR OF OPTICAL DEPTHS FOR CLOUD ! TAU_GAS = VECTOR OF OPTICAL DEPTHS FOR ATMOSPHERIC GAS ! TAU_RAY = VECTOR OF OPTICAL DEPTHS FOR RAYLEIGH SCATTERING ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! VEC* = WORK VECTORS ! W = VECTOR HOLDING QUADRATURE WEIGHTS ! WVN_FLAG = FLAG TO INDICATE WHETHER OR NOT SAVED MATRICES WILL BE ! USED TO COMPUTE RADIANCES ON CURRENT CALL TO THIS ! SUBROUTINE ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG (0=OFF, 1=ON) ! !******************************************************************************* !INTRINSIC SUBPROGRAMS USED BY RAD********************************************** ! MATMUL !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY RAD*********************************************** ! BUILD_LAYER,COMBINE_GS,COMBINE_GT_GR,COMBINE_LAYERS, ! COMBINE_LAYERS_2,GETQUAD2,MATIDENT,MM_IG1G2,SURF_REF1_3, ! SURF_REF2_2 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,NUMLAY,SURF,FLUXES,WARNING_FLAG,R_IO,FLUX_IO,& L_IO,B_IO,INV_IO,KFLAG,SRC_FLAG,VIEW_FLAG,PERT_FLAG,LAYCHG INTEGER, DIMENSION(NUMLAY) :: & DELTA_M_AERO,DELTA_M_CLD,PF_AERO,PF_CLD DOUBLE PRECISION :: & ALBEDO,BSURF,EPS,FSUN,FOURIER_TOL,MU0,PHI,VIEW_ANGLE DOUBLE PRECISION, DIMENSION(N) :: & ITM DOUBLE PRECISION, DIMENSION(NUMLAY) :: & G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD DOUBLE PRECISION, DIMENSION(NUMLAY+1) :: & B_T DOUBLE PRECISION, DIMENSION(0:(4*N)) :: & CHI_SURF DOUBLE PRECISION, DIMENSION(0:(4*N),NUMLAY) :: & CHI_AERO,CHI_CLD CHARACTER (LEN=3) :: & WVN_FLAG LOGICAL :: & AZIMUTHAL_RAD_ONLY,DIMENSIONS_RESET,SZA_RESET,QUAD_RESET !OUTPUT VARIABLES DOUBLE PRECISION :: & FIB_TOT,FOB_TOT,FIT_TOT,FOT_TOT DOUBLE PRECISION, DIMENSION(N) :: & IBMTOT,IBPTOT,ITPTOT !INTERNAL VARIABLES INTEGER :: & NS,NUMDEG,DEGREE,I,J,K,L,LAY,LAYER,LDA,LDVL,LDVR,LWORK,& STATION_FLAG,BLC_FLAG,NEXP,ATMOS_CHANGED DOUBLE PRECISION :: & SUM,CONSTANT3,SI,FLUXTEST,FLUXTEMP,FTOP_S,& TEST_FTOP,TEST_FBOT,BC_TAU,BC_FTOP DOUBLE PRECISION, DIMENSION(NUMLAY) :: & TAU,FTOP DOUBLE PRECISION, DIMENSION(N) :: & IBM,IBP,ITP,LAST_IBMTOT,LAST_IBPTOT,LAST_ITPTOT,& GSPS1,GSMS1,GSPS2,GSMS2,GSPS_TEMP,GSMS_TEMP,& GSPT1,GSMT1,GSPT2,GSMT2,GSPT_TEMP,GSMT_TEMP,& GSP1,GSM1,& BC_SRCS,MU,W,ROOT,WEIGHT,VEC1,GAMMA,Y,EMISS DOUBLE PRECISION, DIMENSION(N,N) :: & GT2,GR2,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2,& GTP_TEMP,GTM_TEMP,GRP_TEMP,GRM_TEMP,RG,MAT1 LOGICAL :: & ALL_RADIANCES_CONVERGED DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE :: & FTOP_SAVED,G1_CLD_SAVED,G2_CLD_SAVED,G_CLD_SAVED,& G1_AERO_SAVED,G2_AERO_SAVED,G_AERO_SAVED,& OMEGA_CLD_SAVED,OMEGA_AERO_SAVED,& TAU_CLD_SAVED,TAU_AERO_SAVED,& TAU_RAY_SAVED,TAU_GAS_SAVED,TAU_SAVED DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE, SAVE :: & E DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: & GSPS,GSMS,GSPS_UB,GSMS_UB,GSPS_LB,GSMS_LB,& GSPT,GSMT,GSPT_UB,GSMT_UB,GSPT_LB,GSMT_LB DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: & GT,GR,GTP_UB,GTM_UB,GRP_UB,GRM_UB,& GTP_LB,GTM_LB,GRP_LB,GRM_LB DOUBLE PRECISION, DIMENSION(:,:,:,:,:), ALLOCATABLE, SAVE :: & MAT3_SAVE LOGICAL, SAVE :: & FIRST_RUN = .TRUE. DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 DOUBLE PRECISION F_CLD !FOR FLUX TEST SECTION INTEGER :: & KPTS,Q_IO DOUBLE PRECISION :: & VIEW_MU,FIT_DIR,FIT_DIF,FOB_DIR,FOB_DIF,FTEST !START PROGRAM NS=2*N LDA = N LDVL = N LDVR = N LWORK = N*(N+6) IF (QUAD == 0) THEN IF (DG == 1) THEN NEXP = NS ELSE NEXP = 2*NS END IF ELSE IF (QUAD == 1) THEN NEXP = 2*NS - 2 END IF IF (VIEW_FLAG == 1) THEN IF ((QUAD == 0).AND.(DG == 1)) THEN NEXP = NEXP - 2 ELSE NEXP = NEXP - 4 END IF END IF IF (AZIMUTHAL_RAD_ONLY) THEN NUMDEG = 0 ELSE NUMDEG = NEXP-1 END IF IF (R_IO == 1) THEN PRINT* PRINT*,'ENTERING RAD' END IF !if (wvn_flag == 'NEW') then !print* !print*,'DIMENSIONS_RESET=',DIMENSIONS_RESET !print*,'SZA_RESET=',SZA_RESET !print*,'QUAD_RESET=',QUAD_RESET !print*,'NEXP=',NEXP !print*,'NUMDEG=',NUMDEG !end if !INITIALIZE SOME MATRICES AND VECTORS IBMTOT = 0.0D0 IBPTOT = 0.0D0 ITPTOT = 0.0D0 LAST_IBMTOT = 0.0D0 LAST_IBPTOT = 0.0D0 LAST_ITPTOT = 0.0D0 !print*,'FIRST_RUN=',FIRST_RUN !print*,'DIMENSIONS_RESET=',DIMENSIONS_RESET IF (FIRST_RUN .OR. DIMENSIONS_RESET) THEN IF (FIRST_RUN) THEN FIRST_RUN = .FALSE. ELSE !print*,'at 4' !print*,'DEALLOCATING INSIDE RAD' DEALLOCATE ( E,FTOP_SAVED,& G1_CLD_SAVED,G2_CLD_SAVED,& G_CLD_SAVED,G1_AERO_SAVED,& G2_AERO_SAVED,G_AERO_SAVED,& OMEGA_CLD_SAVED,OMEGA_AERO_SAVED,& TAU_CLD_SAVED,TAU_AERO_SAVED,& TAU_RAY_SAVED,TAU_GAS_SAVED,& TAU_SAVED,& GSPS,GSMS,& GSPS_UB,GSMS_UB,& GSPS_LB,GSMS_LB,& GSPT,GSMT,& GSPT_UB,GSMT_UB,& GSPT_LB,GSMT_LB,& GT,& GTP_UB,GTM_UB,& GTP_LB,GTM_LB,& GR,& GRP_UB,GRM_UB,& GRP_LB,GRM_LB,& MAT3_SAVE ) END IF !print*,'at 5' !print*,'ALLOCATING INSIDE RAD' ALLOCATE ( E(N,N),FTOP_SAVED(NUMLAY),& G1_CLD_SAVED(NUMLAY),G2_CLD_SAVED(NUMLAY),& G_CLD_SAVED(NUMLAY),G1_AERO_SAVED(NUMLAY),& G2_AERO_SAVED(NUMLAY),G_AERO_SAVED(NUMLAY),& OMEGA_CLD_SAVED(NUMLAY),OMEGA_AERO_SAVED(NUMLAY),& TAU_CLD_SAVED(NUMLAY),TAU_AERO_SAVED(NUMLAY),& TAU_RAY_SAVED(NUMLAY),TAU_GAS_SAVED(NUMLAY),& TAU_SAVED(NUMLAY),& GSPS(N,NUMLAY,0:NUMDEG),GSMS(N,NUMLAY,0:NUMDEG),& GSPS_UB(N,NUMLAY,0:NUMDEG),GSMS_UB(N,NUMLAY,0:NUMDEG),& GSPS_LB(N,NUMLAY,0:NUMDEG),GSMS_LB(N,NUMLAY,0:NUMDEG),& GSPT(N,NUMLAY,0:NUMDEG),GSMT(N,NUMLAY,0:NUMDEG),& GSPT_UB(N,NUMLAY,0:NUMDEG),GSMT_UB(N,NUMLAY,0:NUMDEG),& GSPT_LB(N,NUMLAY,0:NUMDEG),GSMT_LB(N,NUMLAY,0:NUMDEG),& GT(N,N,NUMLAY,0:NUMDEG),& GTP_UB(N,N,NUMLAY,0:NUMDEG),GTM_UB(N,N,NUMLAY,0:NUMDEG),& GTP_LB(N,N,NUMLAY,0:NUMDEG),GTM_LB(N,N,NUMLAY,0:NUMDEG),& GR(N,N,NUMLAY,0:NUMDEG),& GRP_UB(N,N,NUMLAY,0:NUMDEG),GRM_UB(N,N,NUMLAY,0:NUMDEG),& GRP_LB(N,N,NUMLAY,0:NUMDEG),GRM_LB(N,N,NUMLAY,0:NUMDEG),& MAT3_SAVE(N,N,NUMLAY,0:NUMDEG,2) ) !if (DIMENSIONS_RESET) then !print*,'inside allocate if block' !print*,allocated(GT),allocated(GR),& ! allocated(GTP_LB),allocated(GTM_LB),& ! allocated(GRP_LB),allocated(GRM_LB),& ! allocated(MAT3_SAVE) !stop !end if !DEFINE IDENTITY MATRIX CALL MATIDENT(N,E) END IF !if (quad == 1) stop !START OF FOURIER COMPONENT LOOP DEGREE = 0 DO !print*,'degree=',degree !START OF ATMOSPHERE-BUILDING BLOCK IF (WVN_FLAG == 'NEW') THEN !BUILD ATMOSPHERE FOR UNPERTURBED SCENE DO LAYER=1,NUMLAY IF (LAYER == 1) THEN FTOP(LAYER) = FSUN ELSE FTOP(LAYER) = FTOP(LAYER-1)*DEXP(-1.0D0*TAU(LAYER-1)/MU0) END IF IF (FTOP(LAYER) < 1.0D-250) FTOP(LAYER) = 0.0D0 FTOP_SAVED(LAYER) = FTOP(LAYER) G1_CLD_SAVED(LAYER) = G1_CLD(LAYER) G2_CLD_SAVED(LAYER) = G2_CLD(LAYER) G_CLD_SAVED(LAYER) = G_CLD(LAYER) G1_AERO_SAVED(LAYER) = G1_AERO(LAYER) G2_AERO_SAVED(LAYER) = G2_AERO(LAYER) G_AERO_SAVED(LAYER) = G_AERO(LAYER) OMEGA_CLD_SAVED(LAYER) = OMEGA_CLD(LAYER) TAU_CLD_SAVED(LAYER) = TAU_CLD(LAYER) OMEGA_AERO_SAVED(LAYER) = OMEGA_AERO(LAYER) TAU_AERO_SAVED(LAYER) = TAU_AERO(LAYER) TAU_RAY_SAVED(LAYER) = TAU_RAY(LAYER) TAU_GAS_SAVED(LAYER) = TAU_GAS(LAYER) CALL BUILD_LAYER(G1_CLD(LAYER),G2_CLD(LAYER),G_CLD(LAYER),& G1_AERO(LAYER),G2_AERO(LAYER),G_AERO(LAYER),& OMEGA_CLD(LAYER),TAU_CLD(LAYER),OMEGA_AERO(LAYER),& TAU_AERO(LAYER),TAU_RAY(LAYER),TAU_GAS(LAYER),& MU0,FTOP(LAYER),QUAD,DG,N,NEXP,DEGREE,NUMDEG,LAYER,NUMLAY,& B_T(LAYER:LAYER+1),DELTA_M_CLD(LAYER),DELTA_M_AERO(LAYER),& LDA,LDVL,LDVR,LWORK,WVN_FLAG,KFLAG,PERT_FLAG,EPS,& WARNING_FLAG,L_IO,B_IO,INV_IO,& SRC_FLAG,VIEW_FLAG,VIEW_ANGLE,& PF_CLD(LAYER),PF_AERO(LAYER),CHI_CLD(1,LAYER),& CHI_AERO(0,LAYER),DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,& TAU(LAYER),& GT(1,1,LAYER,DEGREE),GR(1,1,LAYER,DEGREE),& GSPS(1,LAYER,DEGREE),GSMS(1,LAYER,DEGREE),& GSPT(1,LAYER,DEGREE),GSMT(1,LAYER,DEGREE),F_CLD) !print*,'after building' TAU_SAVED(LAYER) = TAU(LAYER) !print*,'at 1a' IF (WVN_FLAG == 'DUD') THEN !IF ((WVN_FLAG == 'NEW') .AND. (DEGREE == 7)) THEN PRINT* PRINT*,'LAYER = ',LAYER PRINT* PRINT*,'DEGREE = ',DEGREE PRINT* PRINT*,'GT(:,:,LAYER,DEGREE) = ' PRINT*,GT(:,:,LAYER,DEGREE) PRINT* PRINT*,'GR(:,:,LAYER,DEGREE) = ' PRINT*,GR(:,:,LAYER,DEGREE) PRINT* PRINT*,'GSPS(:,LAYER,DEGREE) = ' PRINT*,GSPS(:,LAYER,DEGREE) PRINT* PRINT*,'GSMS(:,LAYER,DEGREE) = ' PRINT*,GSMS(:,LAYER,DEGREE) PRINT* PRINT*,'GSPT(:,LAYER,DEGREE) = ' PRINT*,GSPT(:,LAYER,DEGREE) PRINT* PRINT*,'GSMT(:,LAYER,DEGREE) = ' PRINT*,GSMT(:,LAYER,DEGREE) !STOP END IF IF (LAYER == 1) THEN !SET 1ST SET OF UPPER BLOCK MATRICES AND VECTORS TO 1ST SET !OF MATRICES AND VECTORS !print*,'at 1b' GTP_UB(:,:,LAYER,DEGREE) = GT(:,:,LAYER,DEGREE) GTM_UB(:,:,LAYER,DEGREE) = GT(:,:,LAYER,DEGREE) GRP_UB(:,:,LAYER,DEGREE) = GR(:,:,LAYER,DEGREE) GRM_UB(:,:,LAYER,DEGREE) = GR(:,:,LAYER,DEGREE) IF (SRC_FLAG == 1) THEN GSPS_UB(:,LAYER,DEGREE) = GSPS(:,LAYER,DEGREE) GSMS_UB(:,LAYER,DEGREE) = GSMS(:,LAYER,DEGREE) ELSE IF (SRC_FLAG == 3) THEN GSPT_UB(:,LAYER,DEGREE) = GSPT(:,LAYER,DEGREE) GSMT_UB(:,LAYER,DEGREE) = GSMT(:,LAYER,DEGREE) ELSE IF (SRC_FLAG == 2) THEN GSPS_UB(:,LAYER,DEGREE) = GSPS(:,LAYER,DEGREE) GSMS_UB(:,LAYER,DEGREE) = GSMS(:,LAYER,DEGREE) GSPT_UB(:,LAYER,DEGREE) = GSPT(:,LAYER,DEGREE) GSMT_UB(:,LAYER,DEGREE) = GSMT(:,LAYER,DEGREE) END IF ELSE !COMBINE UPPER BLOCKS CONSISTING OF LAYERS 1 THRU "LAYER-1" TO !MATRICES & VECTORS FROM THE CURRENT LAYER TO BUILD UPPER BLOCKS !CONSISTING OF LAYERS 1 THRU "LAYER" !print*,'at 1c' IF (SRC_FLAG == 1) THEN CALL COMBINE_LAYERS(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPS_UB(1,LAYER-1,DEGREE) ,GSMS_UB(1,LAYER-1,DEGREE),& GT(1,1,LAYER,DEGREE) ,GT(1,1,LAYER,DEGREE),& GR(1,1,LAYER,DEGREE) ,GR(1,1,LAYER,DEGREE),& GSPS(1,LAYER,DEGREE) ,GSMS(1,LAYER,DEGREE),& GTP_UB(1,1,LAYER,DEGREE) ,GTM_UB(1,1,LAYER,DEGREE),& GRP_UB(1,1,LAYER,DEGREE) ,GRM_UB(1,1,LAYER,DEGREE),& GSPS_UB(1,LAYER,DEGREE) ,GSMS_UB(1,LAYER,DEGREE)) ELSE IF (SRC_FLAG == 3) THEN CALL COMBINE_LAYERS(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPT_UB(1,LAYER-1,DEGREE) ,GSMT_UB(1,LAYER-1,DEGREE),& GT(1,1,LAYER,DEGREE) ,GT(1,1,LAYER,DEGREE),& GR(1,1,LAYER,DEGREE) ,GR(1,1,LAYER,DEGREE),& GSPT(1,LAYER,DEGREE) ,GSMT(1,LAYER,DEGREE),& GTP_UB(1,1,LAYER,DEGREE) ,GTM_UB(1,1,LAYER,DEGREE),& GRP_UB(1,1,LAYER,DEGREE) ,GRM_UB(1,1,LAYER,DEGREE),& GSPT_UB(1,LAYER,DEGREE) ,GSMT_UB(1,LAYER,DEGREE)) ELSE IF (SRC_FLAG == 2) THEN CALL COMBINE_LAYERS_2(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPS_UB(1,LAYER-1,DEGREE) ,GSMS_UB(1,LAYER-1,DEGREE),& GSPT_UB(1,LAYER-1,DEGREE) ,GSMT_UB(1,LAYER-1,DEGREE),& GT(1,1,LAYER,DEGREE) ,GT(1,1,LAYER,DEGREE),& GR(1,1,LAYER,DEGREE) ,GR(1,1,LAYER,DEGREE),& GSPS(1,LAYER,DEGREE) ,GSMS(1,LAYER,DEGREE),& GSPT(1,LAYER,DEGREE) ,GSMT(1,LAYER,DEGREE),& GTP_UB(1,1,LAYER,DEGREE) ,GTM_UB(1,1,LAYER,DEGREE),& GRP_UB(1,1,LAYER,DEGREE) ,GRM_UB(1,1,LAYER,DEGREE),& GSPS_UB(1,LAYER,DEGREE) ,GSMS_UB(1,LAYER,DEGREE),& GSPT_UB(1,LAYER,DEGREE) ,GSMT_UB(1,LAYER,DEGREE)) END IF END IF END DO !print*,'at 1d' GTP1 = GTP_UB(:,:,NUMLAY,DEGREE) GTM1 = GTM_UB(:,:,NUMLAY,DEGREE) GRP1 = GRP_UB(:,:,NUMLAY,DEGREE) GRM1 = GRM_UB(:,:,NUMLAY,DEGREE) IF (SRC_FLAG == 1) THEN GSPS1 = GSPS_UB(:,NUMLAY,DEGREE) GSMS1 = GSMS_UB(:,NUMLAY,DEGREE) ELSE IF (SRC_FLAG == 3) THEN GSPT1 = GSPT_UB(:,NUMLAY,DEGREE) GSMT1 = GSMT_UB(:,NUMLAY,DEGREE) ELSE IF (SRC_FLAG == 2) THEN GSPS1 = GSPS_UB(:,NUMLAY,DEGREE) GSMS1 = GSMS_UB(:,NUMLAY,DEGREE) GSPT1 = GSPT_UB(:,NUMLAY,DEGREE) GSMT1 = GSMT_UB(:,NUMLAY,DEGREE) END IF !print*,'at 1e' !print*,allocated(GT),allocated(GR),& ! allocated(GTP_LB),allocated(GTM_LB),& ! allocated(GRP_LB),allocated(GRM_LB),& ! allocated(MAT3_SAVE) IF (KFLAG == 1) THEN !BUILD LOWER BLOCKS FOR JACOBIAN CALCULATIONS DO LAYER=NUMLAY,1,-1 IF (LAYER == NUMLAY) THEN !SET LAST SET OF LOWER BLOCK MATRICES TO SET !OF MATRICES FROM LAST LAYER GTP_LB(:,:,LAYER,DEGREE) = GT(:,:,LAYER,DEGREE) GTM_LB(:,:,LAYER,DEGREE) = GT(:,:,LAYER,DEGREE) GRP_LB(:,:,LAYER,DEGREE) = GR(:,:,LAYER,DEGREE) GRM_LB(:,:,LAYER,DEGREE) = GR(:,:,LAYER,DEGREE) MAT3_SAVE(:,:,LAYER,DEGREE,:) = 0.0D0 ELSE !COMBINE MATRICES FROM CURRENT LAYER TO THOSE OF !PREVIOUS LOWER BLOCKS CONSISTING OF LAYERS "LAYER+1" THRU !"NUMLAY" TO BUILD LOWER BLOCK CONSISTING OF LAYERS "LAYER" !THRU "NUMLAY" CALL COMBINE_GT_GR(N,E,& GT(1,1,LAYER,DEGREE) ,GT(1,1,LAYER,DEGREE),& GR(1,1,LAYER,DEGREE) ,GR(1,1,LAYER,DEGREE),& GTP_LB(1,1,LAYER+1,DEGREE),GTM_LB(1,1,LAYER+1,DEGREE),& GRP_LB(1,1,LAYER+1,DEGREE),GRM_LB(1,1,LAYER+1,DEGREE),& GTP_LB(1,1,LAYER,DEGREE) ,GTM_LB(1,1,LAYER,DEGREE),& GRP_LB(1,1,LAYER,DEGREE) ,GRM_LB(1,1,LAYER,DEGREE),& MAT3_SAVE(1,1,LAYER,DEGREE,1),MAT3_SAVE(1,1,LAYER,DEGREE,2)) END IF END DO END IF !print*,'at 1f' ELSE IF (WVN_FLAG == 'OLD') THEN !PROVIDE ATMOSPHERE FOR PERTURBED SCENE ATMOS_CHANGED = 0 BLC_FLAG = 0 LAYER = 0 DO !SEARCH FOR A POSSIBLE CHANGED LAYER BY COMPARING THE FOLLOWING !11 OPTICAL PROPERTIES FOR EACH LAYER FROM THE PERTURBED SCENE !WITH THOSE FROM THE UNPERTURBED SCENE LAYER = LAYER + 1 IF( (G1_CLD(LAYER) /= G1_CLD_SAVED(LAYER)) .OR. & (G2_CLD(LAYER) /= G2_CLD_SAVED(LAYER)) .OR. & (G_CLD(LAYER) /= G_CLD_SAVED(LAYER)) .OR. & (G1_AERO(LAYER) /= G1_AERO_SAVED(LAYER)) .OR. & (G2_AERO(LAYER) /= G2_AERO_SAVED(LAYER)) .OR. & (G_AERO(LAYER) /= G_AERO_SAVED(LAYER)) .OR. & (OMEGA_CLD(LAYER) /= OMEGA_CLD_SAVED(LAYER)) .OR. & (TAU_CLD(LAYER) /= TAU_CLD_SAVED(LAYER)) .OR. & (OMEGA_AERO(LAYER) /= OMEGA_AERO_SAVED(LAYER)) .OR. & (TAU_AERO(LAYER) /= TAU_AERO_SAVED(LAYER)) .OR. & (TAU_RAY(LAYER) /= TAU_RAY_SAVED(LAYER)) .OR. & (TAU_GAS(LAYER) /= TAU_GAS_SAVED(LAYER)) ) THEN ! IF((DABS(G1_CLD(LAYER) - G1_CLD_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(G2_CLD(LAYER) - G2_CLD_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(G_CLD(LAYER) - G_CLD_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(G1_AERO(LAYER) - G1_AERO_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(G2_AERO(LAYER) - G2_AERO_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(G_AERO(LAYER) - G_AERO_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(OMEGA_CLD(LAYER) - OMEGA_CLD_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(TAU_CLD(LAYER) - TAU_CLD_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(OMEGA_AERO(LAYER) - OMEGA_AERO_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(TAU_AERO(LAYER) - TAU_AERO_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(TAU_RAY(LAYER) - TAU_RAY_SAVED(LAYER)) > 1.0D-12) ! .OR.& ! (DABS(TAU_GAS(LAYER) - TAU_GAS_SAVED(LAYER)) > 1.0D-12) ) THEN ATMOS_CHANGED = 1 EXIT ELSE IF(LAYER == NUMLAY) THEN EXIT END IF END DO ! IF (LAYCHG > 0) THEN ! !A LAYER HAS CHANGED ! ATMOS_CHANGED = 1 ! LAYER = LAYCHG ! END IF IF (ATMOS_CHANGED == 1) THEN !IF ATMOSPHERE HAS CHANGED, RE-BUILD IT BY ONLY RE-BUILDING !THE LAYER THAT HAS BEEN AFFECTED AND THEN COMBINING IT !WITH THOSE THAT HAVE BEEN UNAFFECTED AND ARE CURRENTLY SAVED IN !MEMORY FROM THE BUILDING OF THE UNPERTURBED SCENE !IF BOTTOM LAYER CHANGED, SET BLC_FLAG IF (LAYER == NUMLAY) BLC_FLAG = 1 !REBUILD THE AFFECTED LAYER FTOP(LAYER) = FTOP_SAVED(LAYER) CALL BUILD_LAYER(G1_CLD(LAYER),G2_CLD(LAYER),G_CLD(LAYER),& G1_AERO(LAYER),G2_AERO(LAYER),G_AERO(LAYER),& OMEGA_CLD(LAYER),TAU_CLD(LAYER),OMEGA_AERO(LAYER),& TAU_AERO(LAYER),TAU_RAY(LAYER),TAU_GAS(LAYER),& MU0,FTOP(LAYER),QUAD,DG,N,NEXP,DEGREE,NUMDEG,LAYER,NUMLAY,& B_T(LAYER:LAYER+1),DELTA_M_CLD(LAYER),DELTA_M_AERO(LAYER),& LDA,LDVL,LDVR,LWORK,WVN_FLAG,KFLAG,PERT_FLAG,EPS,& WARNING_FLAG,L_IO,B_IO,INV_IO,& SRC_FLAG,VIEW_FLAG,VIEW_ANGLE,& PF_CLD(LAYER),PF_AERO(LAYER),CHI_CLD(1,LAYER),& CHI_AERO(0,LAYER),DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,& TAU(LAYER),& GTP1,GRP1,GSPS1,GSMS1,GSPT1,GSMT1,F_CLD) !print*,'after re-building' GTM1 = GTP1 GRM1 = GRP1 !print*, gtm1, 'gtm1' !print*, grm1, 'grm1' IF (LAYER /= NUMLAY) THEN !IF LAYER JUST REBUILT IS NOT LAYER "NUMLAY" (I.E. THE BOTTOM ONE), !THEN COMBINE MATRICES AND VECTORS FROM CURRENT LAYER TO THOSE OF !THE LOWER BLOCK CONSISTING OF LAYERS "LAYER+1" THRU !"NUMLAY" TO BUILD LOWER BLOCK CONSISTING OF LAYERS "LAYER" !THRU "NUMLAY" IF ((SRC_FLAG == 1).OR.(SRC_FLAG == 2)) THEN !IF USING SOLAR SOURCES, THEN !COMPUTE NEW GLOBAL SOLAR SOURCE VECTORS FOR BLOCK OF LAYERS !"LAYER+1" THRU "NUMLAY" DUE TO CHANGE IN SOLAR FLUX EXITING !THE PERTURBED LAYER "LAYER" JUST BUILT ... DO LAY=LAYER+1,NUMLAY IF (LAY == LAYER+1) THEN FTOP(LAY) = FTOP(LAY-1)*DEXP(-1.0D0*TAU(LAY-1)/MU0) ELSE FTOP(LAY) = FTOP(LAY-1)*DEXP(-1.0D0*TAU_SAVED(LAY-1)/MU0) END IF IF (FTOP(LAY) < 1.0D-250) FTOP(LAY) = 0.0D0 END DO DO LAY=NUMLAY,LAYER+1,-1 !RECOMPUTE THE SOLAR SOURCES FOR THE LAYERS BENEATH IF (LAY == NUMLAY) THEN IF (FTOP_SAVED(LAY) >= 1.0D-250) THEN GSPS_LB(:,LAY,DEGREE) = (FTOP(LAY)/FTOP_SAVED(LAY))* & GSPS(:,LAY,DEGREE) GSMS_LB(:,LAY,DEGREE) = (FTOP(LAY)/FTOP_SAVED(LAY))* & GSMS(:,LAY,DEGREE) ELSE GSPS_LB(:,LAY,DEGREE) = 0.0D0 GSMS_LB(:,LAY,DEGREE) = 0.0D0 END IF ELSE IF (FTOP_SAVED(LAY) >= 1.0D-250) THEN GSPS2 = (FTOP(LAY)/FTOP_SAVED(LAY))* & GSPS(:,LAY,DEGREE) GSMS2 = (FTOP(LAY)/FTOP_SAVED(LAY))* & GSMS(:,LAY,DEGREE) ELSE GSPS2 = 0.0D0 GSMS2 = 0.0D0 END IF CALL COMBINE_GS(N,& GR(1,1,LAY,DEGREE),& GSPS2,GSMS2,& GRM_LB(1,1,LAY+1,DEGREE),& GSPS_LB(1,LAY+1,DEGREE),GSMS_LB(1,LAY+1,DEGREE),& MAT3_SAVE(1,1,LAY,DEGREE,1),& MAT3_SAVE(1,1,LAY,DEGREE,2),& GSPS_LB(1,LAY,DEGREE),GSMS_LB(1,LAY,DEGREE)) END IF END DO END IF IF ((SRC_FLAG == 2).OR.(SRC_FLAG == 3)) THEN !IF USING THERMAL SOURCES, THEN !COMPUTE GLOBAL THERMAL SOURCE VECTORS FOR BLOCK OF LAYERS !"LAYER+1" THRU "NUMLAY" DO LAY=NUMLAY,LAYER+1,-1 !COMPUTE THE THERMAL SOURCES IF (LAY == NUMLAY) THEN GSPT_LB(:,LAY,DEGREE) = GSPT(:,LAY,DEGREE) GSMT_LB(:,LAY,DEGREE) = GSMT(:,LAY,DEGREE) ELSE GSPT2 = GSPT(:,LAY,DEGREE) GSMT2 = GSMT(:,LAY,DEGREE) CALL COMBINE_GS(N,& GR(1,1,LAY,DEGREE),& GSPT2,GSMT2,& GRM_LB(1,1,LAY+1,DEGREE),& GSPT_LB(1,LAY+1,DEGREE),GSMT_LB(1,LAY+1,DEGREE),& MAT3_SAVE(1,1,LAY,DEGREE,1),MAT3_SAVE(1,1,LAY,DEGREE,2),& GSPT_LB(1,LAY,DEGREE),GSMT_LB(1,LAY,DEGREE)) END IF END DO END IF !COMBINE THE PERTURBED LAYER WITH BLOCK OF LAYERS "LAYER+1" !THRU "NUMLAY" IF (SRC_FLAG == 1) THEN CALL COMBINE_LAYERS(N,E,& GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1,& GTP_LB(1,1,LAYER+1,DEGREE),GTM_LB(1,1,LAYER+1,DEGREE),& GRP_LB(1,1,LAYER+1,DEGREE),GRM_LB(1,1,LAYER+1,DEGREE),& GSPS_LB(1,LAYER+1,DEGREE) ,GSMS_LB(1,LAYER+1,DEGREE),& GTP_TEMP,GTM_TEMP,GRP_TEMP,GRM_TEMP,GSPS_TEMP,GSMS_TEMP) ELSE IF (SRC_FLAG == 3) THEN CALL COMBINE_LAYERS(N,E,& GTP1,GTM1,GRP1,GRM1,GSPT1,GSMT1,& GTP_LB(1,1,LAYER+1,DEGREE),GTM_LB(1,1,LAYER+1,DEGREE),& GRP_LB(1,1,LAYER+1,DEGREE),GRM_LB(1,1,LAYER+1,DEGREE),& GSPT_LB(1,LAYER+1,DEGREE) ,GSMT_LB(1,LAYER+1,DEGREE),& GTP_TEMP,GTM_TEMP,GRP_TEMP,GRM_TEMP,GSPT_TEMP,GSMT_TEMP) ELSE IF (SRC_FLAG == 2) THEN CALL COMBINE_LAYERS_2(N,E,& GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1,GSPT1,GSMT1,& GTP_LB(1,1,LAYER+1,DEGREE),GTM_LB(1,1,LAYER+1,DEGREE),& GRP_LB(1,1,LAYER+1,DEGREE),GRM_LB(1,1,LAYER+1,DEGREE),& GSPS_LB(1,LAYER+1,DEGREE) ,GSMS_LB(1,LAYER+1,DEGREE),& GSPT_LB(1,LAYER+1,DEGREE) ,GSMT_LB(1,LAYER+1,DEGREE),& GTP_TEMP,GTM_TEMP,GRP_TEMP,GRM_TEMP,& GSPS_TEMP,GSMS_TEMP,GSPT_TEMP,GSMT_TEMP) END IF GTP1 = GTP_TEMP GTM1 = GTM_TEMP GRP1 = GRP_TEMP GRM1 = GRM_TEMP IF (SRC_FLAG == 1) THEN GSPS1 = GSPS_TEMP GSMS1 = GSMS_TEMP ELSE IF (SRC_FLAG == 3) THEN GSPT1 = GSPT_TEMP GSMT1 = GSMT_TEMP ELSE IF (SRC_FLAG == 2) THEN GSPS1 = GSPS_TEMP GSMS1 = GSMS_TEMP GSPT1 = GSPT_TEMP GSMT1 = GSMT_TEMP END IF END IF IF (LAYER /= 1) THEN !COMBINE BLOCK CONTAINING LAYERS 1 THRU "LAYER-1" WITH BLOCK OF !LAYERS "LAYER" THRU "NUMLAY" GTP2 = GTP1 GTM2 = GTM1 GRP2 = GRP1 GRM2 = GRM1 IF (SRC_FLAG == 1) THEN GSPS2 = GSPS1 GSMS2 = GSMS1 ELSE IF (SRC_FLAG == 3) THEN GSPT2 = GSPT1 GSMT2 = GSMT1 ELSE IF (SRC_FLAG == 2) THEN GSPS2 = GSPS1 GSMS2 = GSMS1 GSPT2 = GSPT1 GSMT2 = GSMT1 END IF IF (SRC_FLAG == 1) THEN CALL COMBINE_LAYERS(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPS_UB(1,LAYER-1,DEGREE),GSMS_UB(1,LAYER-1,DEGREE),& GTP2,GTM2,GRP2,GRM2,GSPS2,GSMS2,& GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1) ELSE IF (SRC_FLAG == 3) THEN CALL COMBINE_LAYERS(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPT_UB(1,LAYER-1,DEGREE),GSMT_UB(1,LAYER-1,DEGREE),& GTP2,GTM2,GRP2,GRM2,GSPT2,GSMT2,& GTP1,GTM1,GRP1,GRM1,GSPT1,GSMT1) ELSE IF (SRC_FLAG == 2) THEN CALL COMBINE_LAYERS_2(N,E,& GTP_UB(1,1,LAYER-1,DEGREE),GTM_UB(1,1,LAYER-1,DEGREE),& GRP_UB(1,1,LAYER-1,DEGREE),GRM_UB(1,1,LAYER-1,DEGREE),& GSPS_UB(1,LAYER-1,DEGREE),GSMS_UB(1,LAYER-1,DEGREE),& GSPT_UB(1,LAYER-1,DEGREE),GSMT_UB(1,LAYER-1,DEGREE),& GTP2,GTM2,GRP2,GRM2,GSPS2,GSMS2,GSPT2,GSMT2,& GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1,GSPT1,GSMT1) END IF END IF ELSE !USE THE ATMOSPHERE SAVED IN MEMORY FROM THE BUILDING OF THE !UNPERTURBED SCENE GTP1 = GTP_UB(:,:,NUMLAY,DEGREE) GTM1 = GTM_UB(:,:,NUMLAY,DEGREE) GRP1 = GRP_UB(:,:,NUMLAY,DEGREE) GRM1 = GRM_UB(:,:,NUMLAY,DEGREE) IF (SRC_FLAG == 1) THEN GSPS1 = GSPS_UB(:,NUMLAY,DEGREE) GSMS1 = GSMS_UB(:,NUMLAY,DEGREE) ELSE IF (SRC_FLAG == 3) THEN GSPT1 = GSPT_UB(:,NUMLAY,DEGREE) GSMT1 = GSMT_UB(:,NUMLAY,DEGREE) ELSE IF (SRC_FLAG == 2) THEN GSPS1 = GSPS_UB(:,NUMLAY,DEGREE) GSMS1 = GSMS_UB(:,NUMLAY,DEGREE) GSPT1 = GSPT_UB(:,NUMLAY,DEGREE) GSMT1 = GSMT_UB(:,NUMLAY,DEGREE) END IF !END OF PERTURBED ATMOSPHERE BLOCK END IF !END OF ATMOSPHERE-BUILDING BLOCK END IF !ADD SOLAR AND THERMAL GLOBAL SOURCES FOR TOTAL GLOBAL SOURCES IF NECESSARY IF (SRC_FLAG == 1) THEN GSP1 = GSPS1 GSM1 = GSMS1 ELSE IF (SRC_FLAG == 3) THEN GSP1 = GSPT1 GSM1 = GSMT1 ELSE IF (SRC_FLAG == 2) THEN GSP1 = GSPS1 + GSPT1 GSM1 = GSMS1 + GSMT1 END IF IF (WVN_FLAG == 'DUD') THEN !IF (WVN_FLAG == 'NEW') THEN !.AND.(LAYER == 11)) THEN PRINT* PRINT*,'DEGREE=',DEGREE PRINT* PRINT*,'GTP1 = ' PRINT*,GTP1 PRINT* PRINT*,'GTM1 = ' PRINT*,GTM1 PRINT* PRINT*,'GRP1 = ' PRINT*,GRP1 PRINT* PRINT*,'GRM1 = ' PRINT*,GRM1 PRINT* PRINT*,'GSP1 = ' PRINT*,GSP1 PRINT* PRINT*,'GSM1 = ' PRINT*,GSM1 !STOP END IF !PREPARE SURFACE BOUNDARY CONDITION BY... !...CONSTRUCTING SURFACE REFLECTION MATRIX SELECT CASE (SURF) CASE(0) !NO REFLECTING SURFACE RG = 0.0D0 CONSTANT3 = 0.0D0 CASE(1) !ISOTROPIC SURFACE (LAMBERTIAN) IF (DEGREE == 0) THEN CALL SURF_REF1_3(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,SI,RG) RG = (ALBEDO/SI)*RG CONSTANT3 = 0.5D0/SI ELSE RG = 0.0D0 CONSTANT3 = 0.0D0 END IF CASE(2) !SURFACE READ IN FROM FILE CALL SURF_REF2_2(MU0,QUAD,DG,N,NEXP,DEGREE,NUMDEG,VIEW_FLAG,& VIEW_ANGLE,WARNING_FLAG,CHI_SURF(0),& DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,SI,RG,GAMMA) IF (DEGREE == 0) THEN RG = (1.0D0/SI)*RG Y = 1.0D0 EMISS = Y - MATMUL(RG,Y) ELSE RG = (0.5D0/SI)*RG EMISS = 0.0D0 END IF CONSTANT3 = 0.5D0/SI END SELECT IF (R_IO == 1) PRINT*,'CONSTANT3 = ',CONSTANT3 !...AND CONSTRUCTING SURFACE SOURCE VECTOR IF (SURF /= 0) THEN !REFLECTING SURFACE PRESENT !INSURE FTOP & TAU FOR BOTTOM LAYER ARE APPROPRIATE FOR THE !SITUATION IF (WVN_FLAG == 'NEW') THEN !NEW SCENE BC_FTOP = FTOP(NUMLAY) BC_TAU = TAU(NUMLAY) ELSE IF ((WVN_FLAG == 'OLD').AND.(ATMOS_CHANGED == 1) & .AND.(BLC_FLAG == 0)) THEN !SAVED SCENE WITH CHANGE IN AN ATMOS. UPPER LAYER BC_FTOP = FTOP(NUMLAY) BC_TAU = TAU_SAVED(NUMLAY) ELSE IF ((WVN_FLAG == 'OLD').AND.(ATMOS_CHANGED == 1) & .AND.(BLC_FLAG == 1)) THEN !SAVED SCENE WITH CHANGE IN ATMOS. BOTTOM LAYER BC_FTOP = FTOP_SAVED(NUMLAY) BC_TAU = TAU(NUMLAY) ELSE IF ((WVN_FLAG == 'OLD').AND.(ATMOS_CHANGED == 0)) THEN !SAVED SCENE WITH NO CHANGE IN ATMOS. BC_FTOP = FTOP_SAVED(NUMLAY) BC_TAU = TAU_SAVED(NUMLAY) ELSE PRINT*,'EXPERIENCING UNEXPECTED CASE FOR BC_SRCS' STOP END IF IF (SURF == 1) THEN !ISOTROPIC SURFACE (LAMBERTIAN) IF (SRC_FLAG == 1) THEN BC_SRCS = CONSTANT3*(ALBEDO*MU0*BC_FTOP/PI)* & DEXP(-1.0D0*BC_TAU/MU0) ELSE IF (SRC_FLAG == 3) THEN BC_SRCS = (1.0D0 - ALBEDO)*BSURF ELSE IF (SRC_FLAG == 2) THEN BC_SRCS = CONSTANT3*(ALBEDO*MU0*BC_FTOP/PI)* & DEXP(-1.0D0*BC_TAU/MU0) + (1.0D0 - ALBEDO)*BSURF END IF ELSE IF (SURF == 2) THEN !SURFACE READ IN FROM FILE IF (SRC_FLAG == 1) THEN DO I=1,N BC_SRCS(I) = CONSTANT3*(GAMMA(I)*MU0*BC_FTOP/PI)* & DEXP(-1.0D0*BC_TAU/MU0) END DO ELSE IF (SRC_FLAG == 3) THEN DO I=1,N BC_SRCS(I) = EMISS(I)*BSURF END DO ELSE IF (SRC_FLAG == 2) THEN DO I=1,N BC_SRCS(I) = CONSTANT3*(GAMMA(I)*MU0*BC_FTOP/PI)* & DEXP(-1.0D0*BC_TAU/MU0) + EMISS(I)*BSURF END DO END IF END IF ELSE !NO REFLECTING SURFACE PRESENT BC_SRCS = 0.0D0 END IF IF (R_IO == 1) THEN PRINT* PRINT*,'BC_SRCS = ', BC_SRCS END IF !COMPUTE RADIANCES FOR THE CURRENT Mth ORDER OF SCATTER COMPONENT !(IBM, IBP, AND ITP). !PRINT* !PRINT*,'AT STOP 1' !PRINT* !PRINT*,'E = ',E !PRINT* !PRINT*,'GRP1 =',GRP1 !PRINT* !PRINT*,'RG =',RG !PRINT* !PRINT*,'GTM1 =',GTM1 !PRINT* !PRINT*,'ITM =',ITM !PRINT* !PRINT*,'BC_SRCS =',BC_SRCS !PRINT* !PRINT*,'GSM1 =',GSM1 ! print*, gtp1, 'gtp1' ! print*, ibp,'ibp' ! print*, grm1, 'grm1' ! print*, itm, 'itm' MAT1 = E - MATMUL(GRP1,RG) VEC1 = MATMUL(GTM1,ITM) + MATMUL(GRP1,BC_SRCS) + GSM1 CALL MM_IG1G2(N,1,MAT1,VEC1,IBM,INV_IO) IBP = MATMUL(RG,IBM) + BC_SRCS ITP = MATMUL(GTP1,IBP) + MATMUL(GRM1,ITM) + GSP1 !NOW, ACCUMULATE THEM IN VECTORS FOR THE TOTAL RADIANCES AT THE BOTTOM AND TOP !OF THE ATMOSPHERE (IBMTOT, IBPTOT, AND ITPTOT). IBMTOT = IBMTOT + IBM*DCOS(DBLE(DEGREE)*PHI/180.0D0*PI) IBPTOT = IBPTOT + IBP*DCOS(DBLE(DEGREE)*PHI/180.0D0*PI) ITPTOT = ITPTOT + ITP*DCOS(DBLE(DEGREE)*PHI/180.0D0*PI) !DISPLAY RESULTS IF DESIRED IF (R_IO == 1) THEN PRINT* PRINT*,'THE DIFFUSE I- RADIANCE VECTOR AT THE BOA ' // & 'FOR THE M = ',DEGREE,' COMPONENT LOOKS LIKE:' WRITE(*,*) (IBM(I),I=1,N) PRINT* PRINT*,'THE DIFFUSE I+ RADIANCE VECTOR AT THE BOA ' // & 'FOR THE M = ',DEGREE,' COMPONENT LOOKS LIKE:' WRITE(*,*) (IBP(I),I=1,N) PRINT* PRINT*,'THE DIFFUSE I+ RADIANCE VECTOR AT THE TOA ' // & 'FOR THE M = ',DEGREE,' COMPONENT LOOKS LIKE:' WRITE(*,*) (ITP(I),I=1,N) END IF !FLUX SECTION IF (DEGREE == 0) THEN IF (FLUX_IO == 1) THEN PRINT* PRINT*,'*********************' PRINT*,' FLUX TESTS ' PRINT*,'*********************' END IF !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN QUADRATURE OR LOBATTO !QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) !$$$$$$$ !DIRECT FLUX INTO THE TOP OF THE LAYER FIT_DIR = MU0*FSUN IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE DIRECT BEAM FLUX INTO THE TOP OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FIT_DIR END IF !DIFFUSE FLUX INTO THE TOP OF THE LAYER FIT_DIF = 0.0D0 DO I=1,N FIT_DIF = FIT_DIF + W(I)*MU(I)*ITM(I) END DO FIT_DIF = 2.0D0*PI*FIT_DIF IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE DIFFUSE FLUX INTO THE TOP OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FIT_DIF END IF !TOTAL FLUX INTO THE TOP OF THE LAYER FIT_TOT = FIT_DIR + FIT_DIF IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE TOTAL FLUX INTO THE TOP OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FIT_TOT PRINT* PRINT* END IF !$$$$$$$ !DIRECT FLUX OUT OF THE BOTTOM OF THE LAYER FOB_DIR = MU0*BC_FTOP* DEXP(-1.0D0*BC_TAU/MU0) IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE DIRECT BEAM FLUX OUT OF THE BOTTOM OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FOB_DIR END IF !DIFFUSE FLUX OUT OF THE BOTTOM OF THE LAYER FOB_DIF = 0.0D0 DO I=1,N FOB_DIF = FOB_DIF + W(I)*MU(I)*IBM(I) END DO FOB_DIF = 2.0D0*PI*FOB_DIF IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE DIFFUSE FLUX OUT OF THE BOTTOM OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FOB_DIF END IF !TOTAL FLUX OUT OF THE BOTTOM OF THE LAYER FOB_TOT = FOB_DIR + FOB_DIF IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE TOTAL FLUX OUT OF THE BOTTOM OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FOB_TOT PRINT* PRINT* END IF !$$$$$$$ !TOTAL (DIFFUSE) FLUX INTO THE BOTTOM OF THE LAYER FIB_TOT = 0.0D0 DO I=1,N FIB_TOT = FIB_TOT + W(I)*MU(I)*IBP(I) END DO FIB_TOT = 2.0D0*PI*FIB_TOT IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE TOTAL (DIFFUSE) FLUX INTO THE BOTTOM OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FIB_TOT END IF !TOTAL (DIFFUSE) FLUX OUT OF THE TOP OF THE LAYER FOT_TOT = 0.0D0 DO I=1,N FOT_TOT = FOT_TOT + W(I)*MU(I)*ITP(I) END DO FOT_TOT = 2.0D0*PI*FOT_TOT IF (FLUX_IO == 1) THEN PRINT* PRINT*,'THE TOTAL (DIFFUSE) FLUX OUT OF THE TOP OF THE LAYER ' // & 'LOOKS LIKE:' PRINT*,FOT_TOT PRINT* PRINT* END IF IF (FLUX_IO == 1) THEN !$$$$$$$ !FLUX TEST #1: TOA FLUX CONSERVATION (SHOULD BE 0 FOR CONSERVATIVELY !SCATTERING ATMOSPHERE WITH CONSERVATIVELY REFLECTING SURFACE AND !POSITIVE OTHERWISE) FTEST = FIT_TOT - FOT_TOT PRINT* PRINT*,'THE TOA FLUX CONSERVATION TEST LOOKS LIKE:' PRINT*, FTEST !FLUX TEST #2: SURFACE FLUX CONSERVATION (SHOULD BE 0 FOR !CONSERVATIVELY REFLECTING SURFACE AND NEGATIVE FOR !NON-CONSERVATIVELY REFLECTING SURFACE) FTEST = FIB_TOT - FOB_TOT PRINT* PRINT*,'THE SURFACE FLUX CONSERVATION TEST LOOKS LIKE:' PRINT*, FTEST !FLUX TEST #3: TOTAL FLUX CONSERVATION (SHOULD BE 0 FOR !CONSERVATIVELY SCATTERING ATMOSPHERE AND POSITIVE OTHERWISE) FTEST = FIT_TOT - FOT_TOT + FIB_TOT - FOB_TOT PRINT* PRINT*,'THE TOTAL FLUX CONSERVATION TEST LOOKS LIKE:' PRINT*, FTEST PRINT* PRINT*,'*********************' PRINT*,' END OF FLUX TESTS ' PRINT*,'*********************' PRINT* END IF !END OF FLUX SECTION END IF !CHECK FOR FOURIER SERIES CONVERGENCE ALL_RADIANCES_CONVERGED = .TRUE. DO I=1,N !PRINT* !PRINT*, (DABS(IBMTOT(I)-LAST_IBMTOT(I))/IBMTOT(I)) !PRINT*, (DABS(IBPTOT(I)-LAST_IBPTOT(I))/IBPTOT(I)) !PRINT*, (DABS(ITPTOT(I)-LAST_ITPTOT(I))/ITPTOT(I)) !PRINT*,( (DABS(IBMTOT(I)-LAST_IBMTOT(I))/IBMTOT(I)) > FOURIER_TOL) !PRINT*,( (DABS(IBPTOT(I)-LAST_IBPTOT(I))/IBPTOT(I)) > FOURIER_TOL) !PRINT*,( (DABS(ITPTOT(I)-LAST_ITPTOT(I))/ITPTOT(I)) > FOURIER_TOL) !PAUSE IF ( ( (DABS(IBMTOT(I)-LAST_IBMTOT(I))/IBMTOT(I)) & > FOURIER_TOL ).OR. & ( (DABS(IBPTOT(I)-LAST_IBPTOT(I))/IBPTOT(I)) & > FOURIER_TOL ).OR. & ( (DABS(ITPTOT(I)-LAST_ITPTOT(I))/ITPTOT(I)) & > FOURIER_TOL ) ) & ALL_RADIANCES_CONVERGED = .FALSE. END DO !print* !print*,'DEGREE = ',DEGREE !print*,'NUMDEG = ',NUMDEG !print*,'ALL_RADIANCES_CONVERGED = ',ALL_RADIANCES_CONVERGED !print*,'IBMTOT(N) = ',IBMTOT(N),' LAST_IBMTOT(N) = ',LAST_IBMTOT(N) !print*,'IBPTOT(N) = ',IBPTOT(N),' LAST_IBPTOT(N) = ',LAST_IBPTOT(N) !print*,'ITPTOT(N) = ',ITPTOT(N),' LAST_ITPTOT(N) = ',LAST_ITPTOT(N) LAST_IBMTOT = IBMTOT LAST_IBPTOT = IBPTOT LAST_ITPTOT = ITPTOT IF (ALL_RADIANCES_CONVERGED .OR. (DEGREE == NUMDEG)) EXIT DEGREE = DEGREE + 1 !END OF FOURIER COMPONENT LOOP END DO IF (R_IO == 1) THEN PRINT* PRINT*,'THE TOTAL DIFFUSE I- RADIANCE VECTOR AT THE BOA LOOKS LIKE:' WRITE(*,98) (IBMTOT(I),I=1,N) PRINT* PRINT*,'THE TOTAL DIFFUSE I+ RADIANCE VECTOR AT THE BOA LOOKS LIKE:' WRITE(*,98) (IBPTOT(I),I=1,N) PRINT* PRINT*,'THE TOTAL DIFFUSE I+ RADIANCE VECTOR AT THE TOA LOOKS LIKE:' WRITE(*,98) (ITPTOT(I),I=1,N) END IF 95 FORMAT(8(1X,E23.16E2)) 98 FORMAT(1X,E23.16E2) !SURFACE STATUS MESSAGE IF (R_IO == 1) THEN IF (SURF == 0) THEN PRINT* PRINT*,'CURRENT STATUS: TESTING WITH NO SURFACE' ELSE PRINT* PRINT*,'CURRENT STATUS: TESTING WITH SURFACE', SURF END IF END IF !END PROGRAM IF (R_IO == 1) THEN PRINT* PRINT*,'LEAVING RAD' END IF END SUBROUTINE RAD !******************************************************************************* !******************************************************************************* SUBROUTINE BASIC_PROPERTIES (NEXP,TAU_RAY_IN,TAU_GAS_IN,& OMEGA_CLD_IN,TAU_CLD_IN,DELTA_M_CLD,PF_CLD,G1_CLD,G2_CLD,G_CLD,& CHI_CLD,OMEGA_AERO_IN,TAU_AERO_IN,G1_AERO,G2_AERO,G_AERO,& DELTA_M_AERO,PF_AERO,CHI_AERO,WARNING_FLAG,SRC_FLAG,L_IO,& OMEGA,TAU,X,F_CLD) !INPUT : NEXP,TAU_RAY,TAU_GAS,OMEGA_CLD,TAU_CLD,DELTA_M_CLD, ! PF_CLD,G1_CLD,G2_CLD,G_CLD,CHI_CLD,OMEGA_AERO,TAU_AERO, ! G1_AERO,G2_AERO,G_AERO,DELTA_M_AERO,PF_AERO,CHI_AERO, ! WARNING_FLAG,SRC_FLAG,L_IO !OUTPUT: OMEGA,TAU,X !THIS PROGRAM COMPUTES THE BASIC OPTICAL PROPERTIES REQUIRED FOR THE LAYER !PROGRAMMER: MATT CHRISTI !DATE: 4/16/04 !DATA DICTIONARY******************************************************** ! ! CHI_AERO = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_CLD = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED) ! DELTA_M_AERO = FLAG TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED FOR AEROSOL PHASE FUNCTION IN SUBROUTINE LOCAL ! (0=NO,1=YES) ! DELTA_M_CLD = FLAG TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED FOR CLOUD PHASE FUNCTION IN SUBROUTINE LOCAL ! (0=NO,1=YES) ! G*_AERO = AEROSOL ASYMMETRY PARAMETERS ! G*_CLD = CLOUD ASYMMETRY PARAMETERS ! L_IO = I/O FLAG (0=OFF, 1=ON) ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTION ! OMEGA = SINGLE SCATTER ALBEDO ! OMEGA_* = SINGLE SCATTER ALBEDO FOR CONSTITUENT "*" ! PF_AERO = FLAG TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE ! FUNCTION EXPANSION ARE OBTAINED ! PF_CLD = FLAG TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE ! FUNCTION EXPANSION ARE OBTAINED ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES (1=SOLAR ONLY,2=BOTH,3=THERMAL ONLY) ! TAU = TOTAL OPTICAL DEPTH OF THE LAYER ! TAU_* = OPTICAL DEPTH OF THE LAYER FOR CONSTITUENT "*" ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG (0=OFF, 1=ON) ! X = VECTOR HOLDING COEFFICIENTS OF COMPOSITE PHASE FUNCTION ! X_* = VECTOR HOLDING COEFFICIENTS OF PHASE FUNCTION FOR ! CONSTITUENT "*" !*********************************************************************** !INTRINSIC SUBPROGRAMS USED BY BASIC_PROPERTIES************************* ! NONE !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY BASIC_PROPERTIES************************** ! NONE !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: & NEXP,DELTA_M_CLD,DELTA_M_AERO,L_IO,PF_CLD,PF_AERO,& WARNING_FLAG,SRC_FLAG DOUBLE PRECISION :: & G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,TAU_RAY_IN,& TAU_GAS_IN,TAU_AERO_IN,TAU_CLD_IN,OMEGA_AERO_IN,OMEGA_CLD_IN DOUBLE PRECISION, DIMENSION(0:NEXP) :: & CHI_AERO,CHI_CLD !OUTPUT VARIABLES DOUBLE PRECISION :: & OMEGA,TAU DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & X !INTERNAL VARIABLES INTEGER :: & I,L,DUMMY DOUBLE PRECISION :: & F_CLD,F_AERO,B_CLD,B_AERO,TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,& OMEGA_RAY,OMEGA_AERO,OMEGA_CLD CHARACTER (LEN=30) :: & PHASE_INFILE DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & X_RAY,X_CLD,X_AERO DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 !START PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'ENTERING BASIC_PROPERTIES' END IF !DEFINE COMPONENT OPTICAL DEPTHS AND SINGLE SCATTER ALBEDOS TAU_RAY = TAU_RAY_IN TAU_CLD = TAU_CLD_IN TAU_AERO = TAU_AERO_IN TAU_GAS = TAU_GAS_IN IF (L_IO == 1) THEN PRINT* PRINT*,'TAU_RAY = ',TAU_RAY PRINT*,'TAU_CLD = ',TAU_CLD PRINT*,'TAU_AERO = ',TAU_AERO PRINT*,'TAU_GAS = ',TAU_GAS END IF IF (TAU_RAY > 1.0D-12) THEN OMEGA_RAY = 1.0D0 ELSE OMEGA_RAY = 0.0D0 END IF IF (TAU_CLD > 1.0D-12) THEN OMEGA_CLD = OMEGA_CLD_IN ELSE OMEGA_CLD = 0.0D0 END IF IF (TAU_AERO > 1.0D-12) THEN OMEGA_AERO = OMEGA_AERO_IN ELSE OMEGA_AERO = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'OMEGA_RAY = ',OMEGA_RAY PRINT*,'OMEGA_CLD = ',OMEGA_CLD PRINT*,'OMEGA_AERO = ',OMEGA_AERO END IF !CREATE COEFFICIENTS OF RAYLEIGH PHASE FUNCTION FOR CLEAR ATMOSPHERE X_RAY = 0.0D0 X_RAY(0) = 1.0D0 X_RAY(2) = 0.5D0 IF (L_IO == 1) THEN PRINT* PRINT*,'THE X_RAY(L) ARE:' WRITE(*,20) (X_RAY(L),L=0,NEXP-1) END IF !CREATE COEFFICIENTS OF CLOUD PHASE FUNCTION IF NECESSARY IF (TAU_CLD > 0.0D0) THEN IF (PF_CLD == 1) THEN !CREATE COEFFICIENTS OF NORMAL OR DOUBLE HENYEY-GREENSTEIN PHASE !FUNCTION FOR CLOUD PHASE FUNCTION X_CLD = 0.0D0 X_CLD(0) = 1.0D0 !COMPUTE COMPOSITE COEFFICIENTS IF (G2_CLD < -1.0D-12) THEN !FOR DHG !print* !print*,'doing DHG' B_CLD = (G_CLD - G2_CLD)/(G1_CLD - G2_CLD) DO L=1,NEXP-1 X_CLD(L) = B_CLD*(G1_CLD**L) + (1.0D0-B_CLD)*(G2_CLD**L) END DO ELSE !FOR HG !print* !print*,'doing HG' DO L=1,NEXP-1 X_CLD(L) = G1_CLD**L END DO END IF IF (DELTA_M_CLD == 1) THEN IF (G2_CLD < -1.0D-12) THEN !FOR DHG F_CLD = B_CLD*(G1_CLD**NEXP) + (1.0D0-B_CLD)*(G2_CLD**NEXP) ELSE !FOR HG F_CLD = G1_CLD**NEXP END IF TAU_CLD = (1.0D0 - OMEGA_CLD*F_CLD)*TAU_CLD OMEGA_CLD = ((1.0D0 - F_CLD)/(1.0D0 - OMEGA_CLD*F_CLD))* & OMEGA_CLD DO L=1,NEXP-1 X_CLD(L) = (X_CLD(L) - F_CLD)/(1.0D0 - F_CLD) END DO END IF DO L=1,NEXP-1 X_CLD(L) = (2.0D0*DBLE(L)+1.0D0)*X_CLD(L) END DO ELSE IF (PF_CLD == 2) THEN omega_cld=0.9999 X_CLD = 0.0D0 !x_cld(0)= 1.00003 !x_cld(1)= 1.59844 !x_cld(2)= 2.05035 !x_cld(3)= 2.15446 !x_cld(4)= 2.17113 !x_cld(5)= 2.14180 !x_cld(6)= 2.00745 !x_cld(7)= 1.80988 !x_cld(8)= 1.63911 !x_cld(9)= 1.49284 !x_cld(10)= 1.31392 !x_cld(11)= 1.10192 !x_cld(12)= 0.904483 !x_cld(13)= 0.743844 !x_cld(14)= 0.617412 !x_cld(15)= 0.529541 !x_cld(16)= 0.476313 !x_cld(17)= 0.428447 !x_cld(18)= 0.365618 !x_cld(19)= 0.299827 !x_cld(20)= 0.247628 !x_cld(21)= 0.205038 !x_cld(22)= 0.164951 !x_cld(23)= 0.129507 !x_cld(24)= 0.0980822 !x_cld(25)= 0.0647192 !x_cld(26)= 0.0349279 !x_cld(27)= 0.0211352 !x_cld(28)= 0.0228290 !x_cld(29)= 0.0263821 !x_cld(30)= 0.0248289 !x_cld(31)= 0.019168 !DEFINE COEFFICIENTS FOR CLOUD PHASE FUNCTION USING CHI_CLD ! X_CLD = CHI_CLD(0:NEXP-1) TAU_CLD = (1.0D0 - OMEGA_CLD*F_CLD)*TAU_CLD OMEGA_CLD = ((1.0D0 - F_CLD)/(1.0D0 - OMEGA_CLD*F_CLD))* & OMEGA_CLD END IF END IF IF(L_IO == 1) THEN PRINT* PRINT*,'THE X_CLD(L) ARE:' WRITE(*,20) (X_CLD(L),L=0,NEXP-1) END IF !CREATE COEFFICIENTS OF AEROSOL PHASE FUNCTION IF NECESSARY IF (TAU_AERO > 0.0D0) THEN IF (PF_AERO == 1) THEN !CREATE COEFFICIENTS OF NORMAL OR DOUBLE HENYEY-GREENSTEIN PHASE !FUNCTION FOR CLOUD PHASE FUNCTION X_AERO = 0.0D0 X_AERO(0) = 1.0D0 !COMPUTE COMPOSITE COEFFICIENTS IF (G2_AERO < -1.0D-12) THEN !FOR DHG B_AERO = (G_AERO - G2_AERO)/(G1_AERO - G2_AERO) DO L=1,NEXP-1 X_AERO(L) = B_AERO*(G1_AERO**L) + (1.0D0-B_AERO)*(G2_AERO**L) END DO ELSE !FOR HG DO L=1,NEXP-1 X_AERO(L) = G1_AERO**L END DO END IF IF (DELTA_M_AERO == 1) THEN IF (G2_AERO < -1.0D-12) THEN !FOR DHG F_AERO = B_AERO*(G1_AERO**NEXP) + (1.0D0-B_AERO)*(G2_AERO**NEXP) ELSE !FOR HG F_AERO = G1_AERO**NEXP END IF TAU_AERO = (1.0D0 - OMEGA_AERO*F_AERO)*TAU_AERO OMEGA_AERO = ((1.0D0 - F_AERO)/(1.0D0 - OMEGA_AERO*F_AERO))* & OMEGA_AERO DO L=1,NEXP-1 X_AERO(L) = (X_AERO(L) - F_AERO)/(1.0D0 - F_AERO) END DO END IF DO L=1,NEXP-1 X_AERO(L) = (2.0D0*DBLE(L)+1.0D0)*X_AERO(L) END DO ELSE IF (PF_AERO == 2) THEN X_AERO = 0.0D0 !DEFINE COEFFICIENTS FOR AEROSOL PHASE FUNCTION USING CHI_AERO X_AERO = CHI_AERO(0:NEXP-1) IF (DELTA_M_AERO == 1) THEN F_AERO = CHI_AERO(NEXP) F_AERO = F_AERO/((2.0D0*DBLE(NEXP))+1.0D0) TAU_AERO = (1.0D0 - OMEGA_AERO*F_AERO)*TAU_AERO OMEGA_AERO = ((1.0D0 - F_AERO)/(1.0D0 - OMEGA_AERO*F_AERO))* & OMEGA_AERO DO L=1,NEXP-1 X_AERO(L) = X_AERO(L)/((2.0D0*DBLE(L))+1.0D0) X_AERO(L) = ((2.0D0*DBLE(L))+1.0D0)* & (X_AERO(L) - F_AERO)/(1.0D0 - F_AERO) END DO END IF END IF ELSE X_AERO = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE X_AERO(L) ARE:' WRITE(*,20) (X_AERO(L),L=0,NEXP-1) END IF !COMPUTE COMPOSITE OPTICAL DEPTH AND SINGLE SCATTER ALBEDO TAU = TAU_RAY + TAU_CLD + TAU_AERO + TAU_GAS IF (L_IO == 1) THEN PRINT* PRINT*,'TAU = ',TAU END IF IF ((TAU_CLD > 1.0D-12) .OR. (TAU_AERO > 1.0D-12) .OR. & (TAU_GAS > 1.0D-12)) THEN OMEGA = (OMEGA_CLD*TAU_CLD + OMEGA_AERO*TAU_AERO + & OMEGA_RAY*TAU_RAY)/TAU ELSE IF(((SRC_FLAG == 1).OR.(SRC_FLAG == 2)) .AND. & (TAU_RAY > 1.0D-12)) THEN OMEGA_RAY = 0.9999999999D0 OMEGA = 0.9999999999D0 OMEGA_RAY = 0.999999D0 OMEGA = 0.999999D0 ELSE OMEGA_RAY = 0.0D0 OMEGA = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'TAU_CLD OR TAU_AERO CONDITION MET FOR OMEGA' END IF END IF IF (L_IO == 1) THEN PRINT* PRINT*,'OMEGA = ',OMEGA END IF !CREATE COEFFICIENTS OF COMPOSITE PHASE FUNCTION IF (OMEGA > 0.0D0) THEN X = ((OMEGA_RAY*TAU_RAY)*X_RAY + (OMEGA_CLD*TAU_CLD)*X_CLD + & (OMEGA_AERO*TAU_AERO)*X_AERO) / (OMEGA*TAU) ELSE X = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE X(L) ARE:' WRITE(*,20) (X(L),L=0,NEXP-1) END IF 20 FORMAT(1X,F16.8) !END PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'LEAVING BASIC_PROPERTIES' END IF END SUBROUTINE BASIC_PROPERTIES !********************************************************************** !********************************************************************** SUBROUTINE BUILD_LAYER(G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,& G_AERO,OMEGA_CLD,TAU_CLD,OMEGA_AERO,TAU_AERO,TAU_RAY,& TAU_GAS,MU0,FTOP,QUAD,DG,N,NEXP,DEGREE,NUMDEG,LAYER,NUMLAY,& B_T,DELTA_M_CLD,DELTA_M_AERO,LDA,LDVL,LDVR,LWORK,& WVN_FLAG,KFLAG,PERT_FLAG,EPS,WARNING_FLAG,& L_IO,B_IO,INV_IO,SRC_FLAG,VIEW_FLAG,VIEW_ANGLE,PF_CLD,PF_AERO,& CHI_CLD,CHI_AERO,DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,& TAU,GT,GR,GSPS,GSMS,GSPT,GSMT, F_CLD) !INPUT : G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,OMEGA_CLD,TAU_CLD, ! OMEGA_AERO,TAU_AERO,TAU_RAY,TAU_GAS,MU0,FTOP,QUAD,DG,N,NEXP,DEGREE, ! NUMDEG,LAYER,NUMLAY,B_T,DELTA_M_CLD,DELTA_M_AERO, ! LDA,LDVL,LDVR,LWORK,WVN_FLAG,KFLAG,PERT_FLAG,EPS,WARNING_FLAG, ! L_IO,B_IO,INV_IO,SRC_FLAG,VIEW_FLAG,VIEW_ANGLE, ! PF_CLD,PF_AERO,CHI_CLD,CHI_AERO !OUTPUT: TAU,GT,GR,GSPS,GSMS,GSPT,GSMT !THIS PROGRAM BUILDS ONE LAYER OF ATMOSPHERE FOR THE PARAMETERS INPUT. !THIS IS DONE THROUGH CONSTRUCTING THE GLOBAL TRANSMISSION AND REFLECTION !MATRICES AND SOURCE VECTORS FOR THE LAYER. !IT IS IN Z-SPACE. !PROGRAMMER: MATT CHRISTI !DATE: 4/16/04 !DATA DICTIONARY******************************************************** ! ! B = MATRIX PRODUCT OF (T-R)*(T+R) ! B_IO = SUBROUTINE BUILD_LAYER I/O FLAG (0=OFF, 1=ON) ! B_T = VECTOR OF SPECTRAL RADIANCES OR RADIANCES (DEPENDING ON VALUE ! OF PLANCK_FLAG) USING THE PLANCK FUNCTION OF EMISSION ! CHI_AERO = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_CLD = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED) ! D = MATRIX OF EIGENVECTORS OF B ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! DELTA_M_AERO = FLAG TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED FOR AEROSOL PHASE FUNCTION IN SUBROUTINE LOCAL ! (0=NO,1=YES) ! DELTA_M_CLD = FLAG TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED FOR CLOUD PHASE FUNCTION IN SUBROUTINE LOCAL ! (0=NO,1=YES) ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! E = IDENTITY MATRIX ! EPS = AMOUNT BY WHICH OPTICAL PROPERTIES ARE PERTURBED BY FOR THE ! COMPUTATION OF JACOBIAN ELEMENTS BY FINITE DIFFERENCE WHEN ! THE PERTURBATION SCHEME IS IN USE ! FTOP = FLUX AT THE TOP OF THE LAYER AT A GIVEN WAVELENGTH ! G*_AERO = AEROSOL ASYMMETRY PARAMETERS ! G*_CLD = CLOUD ASYMMETRY PARAMETERS ! GR = GLOBAL REFLECTION MATRIX ! GSPS = GLOBAL SOURCE VECTOR PLUS SOLAR (UPWELLING) ! GSMS = GLOBAL SOURCE VECTOR MINUS SOLAR (DOWNWELLING) ! GSPT = GLOBAL SOURCE VECTOR PLUS THERMAL (UPWELLING) ! GSMT = GLOBAL SOURCE VECTOR MINUS THERMAL (DOWNWELLING) ! GT = GLOBAL TRANSMISSION MATRIX ! INV_IO = SUBROUTINE INVERT I/O FLAG (0=OFF, 1=ON) (FOR MATRIX INVERSE) ! L_IO = SUBROUTINE LOCAL I/O FLAG (0=OFF, 1=ON) ! LAMBP = VECTOR OF POSITIVE EIGENVALUES OF MATRIX A ! LAYER = LAYER WHOSE OPTICAL PROPERTIES ARE BEING COMPUTED ! MU0 = COSINE OF SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTION ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! OMEGA_* = SINGLE SCATTER ALBEDO FOR CONSTITUENT "*" ! PERT_FLAG = FLAG TO INDICATE WHETHER OR NOT THE PERTURBATION SCHEME WILL ! BE USED TO BUILD LAYERS ON SUBSEQUENT CALLS TO RADIANT DURING ! LAYER SAVING MODE ! (0=NO, 1=YES) ! PF_AERO = FLAG TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE ! FUNCTION EXPANSION ARE OBTAINED ! PF_CLD = FLAG TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE ! FUNCTION EXPANSION ARE OBTAINED ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! SUBROUTINE LOCAL (0=GAUSS,1=LOBATTO) ! S1M = S2M*DECAYING EXPONENTIAL ! S1P = S2P*DECAYING EXPONENTIAL ! S2M = MATRIX-VECTOR PRODUCT OF IMUA_I*SIGMAM ! S2P = MATRIX-VECTOR PRODUCT OF IMUA_I*SIGMAP ! SIGMAM = LOCAL DOWNWELLING SOLAR SOURCE ! SIGMAP = LOCAL UPWELLING SOLAR SOURCE ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES (1=SOLAR ONLY,2=BOTH,3=THERMAL ONLY) ! TAU = TOTAL OPTICAL DEPTH OF THE LAYER ! TAU_* = OPTICAL DEPTH FOR CONSTITUENT "*" ! TMR = MATRIX HOLDING DIFFERENCE OF T AND R MATRICES (T-R) ! TPR = MATRIX HOLDING SUM OF T AND R MATRICES (T+R) ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG (0=OFF, 1=ON) ! OMEGA = SINGLE SCATTER ALBEDO ! WVN_FLAG = FLAG TO INDICATE WHETHER OR NOT SAVED MATRICES WILL BE ! USED TO COMPUTE RADIANCES ON CURRENT CALL TO SUBROUTINE ! RADIANT ! ('NEW'=NO, 'OLD'=YES) ! NOTE: WVN_FLAG USED IN CONJUNCTION WITH KFLAG ! X = VECTOR HOLDING COEFFICIENTS OF COMPOSITE PHASE FUNCTION ! !******************************************************************************* !INTRINSIC SUBPROGRAMS USED BY BUILD_LAYER************************************** ! MATMUL !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY BUILD_LAYER*************************************** ! BASIC_PROPERTIES,GET_GLOBAL_SOURCES5,GET_GT_GR4 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,NEXP,DEGREE,NUMDEG,LAYER,DELTA_M_CLD,DELTA_M_AERO,& L_IO,B_IO,INV_IO,SRC_FLAG,VIEW_FLAG,PF_CLD,PF_AERO,WARNING_FLAG,& NUMLAY DOUBLE PRECISION :: & G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,& MU0,FTOP,VIEW_ANGLE DOUBLE PRECISION, DIMENSION(0:1) :: & B_T DOUBLE PRECISION, DIMENSION(0:NEXP) :: & CHI_AERO,CHI_CLD LOGICAL :: & DIMENSIONS_RESET,SZA_RESET,QUAD_RESET !OUTPUT VARIABLES DOUBLE PRECISION :: & TAU DOUBLE PRECISION, DIMENSION(N) :: & GSPS,GSMS,GSPT,GSMT DOUBLE PRECISION, DIMENSION(N,N) :: & GT,GR !INTERNAL VARIABLES INTEGER :: & DELTA_0_M DOUBLE PRECISION :: & OMEGA DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & X DOUBLE PRECISION, DIMENSION(N) :: & RECMU,PSOLP,PSOLM DOUBLE PRECISION, DIMENSION(N,N) :: & T,R DOUBLE PRECISION F_CLD !FOR SUBROUTINE SGEEVX (FOR COMPUTING EIGENVALUES AND EIGENVECTORS): !INPUT VARIABLES INTEGER :: & LDA,LDVL,LDVR,LWORK !FOR PERTURBATION SCHEME: !INPUT VARIABLES INTEGER :: & KFLAG,PERT_FLAG DOUBLE PRECISION :: & EPS CHARACTER(LEN=3) :: & WVN_FLAG !START PROGRAM IF (B_IO == 1) THEN PRINT* WRITE(*,5) LAYER 5 FORMAT(1X,'ENTERING BUILD_LAYER: BUILDING LAYER #',I3) END IF !DISPLAY INPUT DATA IF (B_IO == 1) THEN PRINT* PRINT*,'SRC_FLAG = ',SRC_FLAG PRINT*,'VIEW_FLAG = ',VIEW_FLAG PRINT*,'VIEW_ANGLE = ',VIEW_ANGLE PRINT*,'MU0 = ', MU0 PRINT*,'FTOP = ', FTOP PRINT*,'QUAD = ', QUAD PRINT*,'DG = ',DG PRINT*,'N = ', N PRINT*,'NEXP = ', NEXP PRINT*,'DEGREE = ', DEGREE PRINT*,'NUMDEG = ', NUMDEG PRINT*,'LAYER = ', LAYER PRINT*,'B_T = ',B_T PRINT* PRINT*,'OMEGA_CLD = ', OMEGA_CLD PRINT*,'TAU_CLD = ', TAU_CLD PRINT*,'G1_CLD = ', G1_CLD PRINT*,'G2_CLD = ', G2_CLD PRINT*,'G_CLD = ', G_CLD PRINT*,'PF_CLD = ',PF_CLD PRINT*,'CHI_CLD = ',CHI_CLD PRINT*,'DELTA_M_CLD = ', DELTA_M_CLD PRINT* PRINT*,'OMEGA_AERO = ', OMEGA_AERO PRINT*,'TAU_AERO = ', TAU_AERO PRINT*,'G1_AERO = ', G1_AERO PRINT*,'G2_AERO = ', G2_AERO PRINT*,'G_AERO = ', G_AERO PRINT*,'PF_AERO = ',PF_AERO PRINT*,'CHI_AERO = ',CHI_AERO PRINT*,'DELTA_M_AERO = ', DELTA_M_AERO PRINT*,'TAU_RAY = ', TAU_RAY PRINT*,'TAU_GAS = ', TAU_GAS END IF !DEFINE SOME PARAMETERS IF (DEGREE == 0) THEN DELTA_0_M = 1.0D0 ELSE DELTA_0_M = 0.0D0 END IF !OBTAIN BASIC OPTICAL PROPERTIES OF THE LAYER: COMPOSITE !SINGLE SCATTER ALBEDO, OPTICAL DEPTH, & PHASE FUNCTION COEFFICIENTS CALL BASIC_PROPERTIES (NEXP,TAU_RAY,TAU_GAS,& OMEGA_CLD,TAU_CLD,DELTA_M_CLD,PF_CLD,G1_CLD,G2_CLD,G_CLD,& CHI_CLD,OMEGA_AERO,TAU_AERO,G1_AERO,G2_AERO,G_AERO,& DELTA_M_AERO,PF_AERO,CHI_AERO,WARNING_FLAG,SRC_FLAG,L_IO,& OMEGA,TAU,X,F_CLD) !COMPUTE GLOBAL TRANSMISSION & REFLECTION MATRICES CALL GET_GT_GR4 (N,LDA,LDVL,LDVR,LWORK,VIEW_FLAG,VIEW_ANGLE,& QUAD,DG,NUMDEG,DEGREE,NEXP,DELTA_0_M,MU0,LAYER,NUMLAY,& WVN_FLAG,KFLAG,PERT_FLAG,EPS,OMEGA,TAU,X,L_IO,INV_IO,& DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,RECMU,PSOLP,PSOLM,T,R,GT,GR) !COMPUTE GLOBAL SOURCE VECTORS CALL GET_GLOBAL_SOURCES5 (SRC_FLAG,N,DEGREE,MU0,OMEGA,& TAU,FTOP,B_T,RECMU,PSOLP,PSOLM,T,R,GT,GR,L_IO,INV_IO,& GSPS,GSMS,GSPT,GSMT) !END PROGRAM IF (B_IO == 1) THEN PRINT* PRINT*,'LEAVING BUILD_LAYER' END IF END SUBROUTINE BUILD_LAYER !******************************************************************************* !******************************************************************************* SUBROUTINE COMBINE_GS(N,GRP1,GSP1,GSM1,GRM2,GSP2,GSM2,& MAT3_SAVE1,MAT3_SAVE2,GSPOUT,GSMOUT) !INPUT : N,GRP1,GSP1,GSM1,GRM2,GSP2,GSM2,MAT3_SAVE1,MAT3_SAVE2 !OUTPUT: GSPOUT,GSMOUT !THIS PROGRAM COMBINES THE GLOBAL SOURCE VECTORS FROM TWO SEPARATE !LAYERS FOR A COMBINED LAYER !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY COMBINE_GS*************************************** ! MATMUL !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N DOUBLE PRECISION, DIMENSION(N) :: & GSP1,GSM1,GSP2,GSM2 DOUBLE PRECISION, DIMENSION(N,N) :: & GRP1,GRM2,MAT3_SAVE1,MAT3_SAVE2 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & GSPOUT,GSMOUT !INTERNAL VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & MAT5,MAT6,MAT7 !START PROGRAM !COMPUTE GSP FOR THE COMBINED LAYER MAT5 = MATMUL(GRM2,GSM1) MAT6 = MAT5 + GSP2 MAT7 = MATMUL(MAT3_SAVE1,MAT6) GSPOUT = MAT7 + GSP1 ! GSPOUT = MATMUL(MAT3_SAVE1,MATMUL(GRM2,GSM1) + GSP2) + GSP1 !COMPUTE GSM FOR THE COMBINED LAYER MAT5 = MATMUL(GRP1,GSP2) MAT6 = MAT5 + GSM1 MAT7 = MATMUL(MAT3_SAVE2,MAT6) GSMOUT = MAT7 + GSM2 ! GSMOUT = MATMUL(MAT3_SAVE2,MATMUL(GRP1,GSP2) + GSM1) + GSM2 END SUBROUTINE COMBINE_GS !******************************************************************************* !******************************************************************************* SUBROUTINE COMBINE_GT_GR(N,E,& GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2,& GTPOUT,GTMOUT,GRPOUT,GRMOUT,MAT3_SAVE1,MAT3_SAVE2) !INPUT : N,E,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2 !OUTPUT: GTPOUT,GTMOUT,GRPOUT,GRMOUT,MAT3_SAVE1,MAT3_SAVE2 !THIS PROGRAM COMBINES ONLY THE GLOBAL TRANSMISSION & REFLECTION MATRICES !FROM TWO SEPARATE LAYERS FOR A COMBINED LAYER !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY COMBINE_GT_GR************************************ ! MATMUL !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY COMBINE_GT_GR************************************* ! MATIDENT,MM_IG1G2 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N DOUBLE PRECISION, DIMENSION(N,N) :: & E,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N,N) :: & GTPOUT,GTMOUT,GRPOUT,GRMOUT,MAT3_SAVE1,MAT3_SAVE2 !INTERNAL VARIABLES DOUBLE PRECISION, DIMENSION(N,N) :: & MAT1,MAT2,MAT3,MAT4,MAT_TEMP INTEGER :: & INV_IO = 0 !START PROGRAM !COMPUTE GTP AND GRM FOR THE COMBINED LAYER MAT1 = MATMUL(GRM2,GRP1) MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTP1) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3_SAVE1 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3_SAVE1,MATMUL(GRM2,GTM1)) GTPOUT = MATMUL(MAT3_SAVE1,GTP2) GRMOUT = MAT4 + GRM1 ! CALL MM_IG1G2(N,N,TRANSPOSE(E - MATMUL(GRM2,GRP1)),& ! TRANSPOSE(GTP1),MAT3,INV_IO) ! MAT3_SAVE1 = TRANSPOSE(MAT3) ! GTPOUT = MATMUL(MAT3_SAVE1,GTP2) ! GRMOUT = MATMUL(MAT3_SAVE1,MATMUL(GRM2,GTM1)) + GRM1 !COMPUTE GTM AND GRP FOR THE COMBINED LAYER MAT1 = MATMUL(GRP1,GRM2) MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTM2) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3_SAVE2 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3_SAVE2,MATMUL(GRP1,GTP2)) GTMOUT = MATMUL(MAT3_SAVE2,GTM1) GRPOUT = MAT4 + GRP2 ! CALL MM_IG1G2(N,N,TRANSPOSE(E - MATMUL(GRP1,GRM2)),& ! TRANSPOSE(GTM2),MAT3,INV_IO) ! MAT3_SAVE2 = TRANSPOSE(MAT3) ! GTMOUT = MATMUL(MAT3_SAVE2,GTM1) ! GRPOUT = MATMUL(MAT3_SAVE2,MATMUL(GRP1,GTP2)) + GRP2 END SUBROUTINE COMBINE_GT_GR !******************************************************************************* !******************************************************************************* SUBROUTINE COMBINE_LAYERS(N,E,& GTP1,GTM1,GRP1,GRM1,GSP1,GSM1,& GTP2,GTM2,GRP2,GRM2,GSP2,GSM2,& GTPOUT,GTMOUT,GRPOUT,GRMOUT,GSPOUT,GSMOUT) !INPUT : N,E,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2,GSP1,GSM1,GSP2,GSM2 !OUTPUT: GTPOUT,GTMOUT,GRPOUT,GRMOUT,GSPOUT,GSMOUT !THIS PROGRAM COMBINES THE GLOBAL TRANSMISSION & REFLECTION MATRICES AND !GLOBAL SOURCE VECTORS FROM TWO SEPARATE LAYERS FOR A COMBINED LAYER !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY COMBINE_LAYERS*********************************** ! MATMUL !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY COMBINE_LAYERS************************************ ! MATIDENT,MM_IG1G2 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N DOUBLE PRECISION, DIMENSION(N) :: & GSP1,GSM1,GSP2,GSM2 DOUBLE PRECISION, DIMENSION(N,N) :: & E,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & GSPOUT,GSMOUT DOUBLE PRECISION, DIMENSION(N,N) :: & GTPOUT,GTMOUT,GRPOUT,GRMOUT !INTERNAL VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & MAT5,MAT6,MAT7 DOUBLE PRECISION, DIMENSION(N,N) :: & MAT1,MAT2,MAT3,MAT4,MAT_TEMP INTEGER :: & INV_IO = 0 !START PROGRAM !COMPUTE GTP, GRM, AND GSP FOR THE COMBINED LAYER MAT1 = MATMUL(GRM2,GRP1) MAT5 = MATMUL(GRM2,GSM1) MAT6 = MAT5 + GSP2 MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTP1) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3,MATMUL(GRM2,GTM1)) MAT7 = MATMUL(MAT3,MAT6) GTPOUT = MATMUL(MAT3,GTP2) GRMOUT = MAT4 + GRM1 GSPOUT = MAT7 + GSP1 ! CALL MM_IG1G2(N,N,TRANSPOSE(E - MATMUL(GRM2,GRP1)),& ! TRANSPOSE(GTP1),MAT3,INV_IO) ! MAT3 = TRANSPOSE(MAT3) ! ! GTPOUT = MATMUL(MAT3,GTP2) ! GRMOUT = MATMUL(MAT3,MATMUL(GRM2,GTM1)) + GRM1 ! GSPOUT = MATMUL(MAT3,MATMUL(GRM2,GSM1) + GSP2) + GSP1 !COMPUTE GTM, GRP, AND GSM FOR THE COMBINED LAYER MAT1 = MATMUL(GRP1,GRM2) MAT5 = MATMUL(GRP1,GSP2) MAT6 = MAT5 + GSM1 MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTM2) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3,MATMUL(GRP1,GTP2)) MAT7 = MATMUL(MAT3,MAT6) GTMOUT = MATMUL(MAT3,GTM1) GRPOUT = MAT4 + GRP2 GSMOUT = MAT7 + GSM2 ! CALL MM_IG1G2(N,N,TRANSPOSE(E - MATMUL(GRP1,GRM2)),& ! TRANSPOSE(GTM2),MAT3,INV_IO) ! MAT3 = TRANSPOSE(MAT3) ! ! GTMOUT = MATMUL(MAT3,GTM1) ! GRPOUT = MATMUL(MAT3,MATMUL(GRP1,GTP2)) + GRP2 ! GSMOUT = MATMUL(MAT3,MATMUL(GRP1,GSP2) + GSM1) + GSM2 END SUBROUTINE COMBINE_LAYERS !******************************************************************************* !******************************************************************************* SUBROUTINE COMBINE_LAYERS_2(N,E,& GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1,GSPT1,GSMT1,& GTP2,GTM2,GRP2,GRM2,GSPS2,GSMS2,GSPT2,GSMT2,& GTPOUT,GTMOUT,GRPOUT,GRMOUT,GSPSOUT,GSMSOUT,& GSPTOUT,GSMTOUT) !INPUT : N,E,GTP1,GTM1,GRP1,GRM1,GSPS1,GSMS1,GSPT1,GSMT1, ! GTP2,GTM2,GRP2,GRM2,GSPS2,GSMS2,GSPT2,GSMT2, !OUTPUT: GTPOUT,GTMOUT,GRPOUT,GRMOUT,GSPSOUT,GSMSOUT,GSPTOUT,GSMTOUT !THIS PROGRAM COMBINES THE GLOBAL TRANSMISSION & REFLECTION MATRICES AND !GLOBAL SOURCE VECTORS FROM TWO SEPARATE LAYERS FOR A COMBINED LAYER !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY COMBINE_LAYERS_2********************************* ! MATMUL !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY COMBINE_LAYERS_2********************************** ! MATIDENT,MM_IG1G2 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N DOUBLE PRECISION, DIMENSION(N) :: & GSPS1,GSMS1,GSPT1,GSMT1,GSPS2,GSMS2,GSPT2,GSMT2 DOUBLE PRECISION, DIMENSION(N,N) :: & E,GTP1,GTM1,GRP1,GRM1,GTP2,GTM2,GRP2,GRM2 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & GSPSOUT,GSMSOUT,GSPTOUT,GSMTOUT DOUBLE PRECISION, DIMENSION(N,N) :: & GTPOUT,GTMOUT,GRPOUT,GRMOUT !INTERNAL VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & MAT5,MAT6,MAT7 DOUBLE PRECISION, DIMENSION(N,N) :: & MAT1,MAT2,MAT3,MAT4,MAT_TEMP INTEGER :: & INV_IO !START PROGRAM !COMPUTE GTP, GRM, GSPS, AND GSPT FOR THE COMBINED LAYER MAT1 = MATMUL(GRM2,GRP1) MAT5 = MATMUL(GRM2,GSMS1) MAT6 = MAT5 + GSPS2 MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTP1) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3,MATMUL(GRM2,GTM1)) MAT7 = MATMUL(MAT3,MAT6) GTPOUT = MATMUL(MAT3,GTP2) GRMOUT = MAT4 + GRM1 GSPSOUT = MAT7 + GSPS1 MAT5 = MATMUL(GRM2,GSMT1) MAT6 = MAT5 + GSPT2 MAT7 = MATMUL(MAT3,MAT6) GSPTOUT = MAT7 + GSPT1 !COMPUTE GTM, GRP, GSMS, AND GSMT FOR THE COMBINED LAYER MAT1 = MATMUL(GRP1,GRM2) MAT5 = MATMUL(GRP1,GSPS2) MAT6 = MAT5 + GSMS1 MAT2 = E - MAT1 MAT_TEMP = TRANSPOSE(GTM2) MAT2 = TRANSPOSE(MAT2) CALL MM_IG1G2(N,N,MAT2,MAT_TEMP,MAT3,INV_IO) MAT3 = TRANSPOSE(MAT3) MAT4 = MATMUL(MAT3,MATMUL(GRP1,GTP2)) MAT7 = MATMUL(MAT3,MAT6) GTMOUT = MATMUL(MAT3,GTM1) GRPOUT = MAT4 + GRP2 GSMSOUT = MAT7 + GSMS2 MAT5 = MATMUL(GRP1,GSPT2) MAT6 = MAT5 + GSMT1 MAT7 = MATMUL(MAT3,MAT6) GSMTOUT = MAT7 + GSMT2 END SUBROUTINE COMBINE_LAYERS_2 !******************************************************************************* !******************************************************************************* SUBROUTINE DATA_INSPECTOR(SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV,& TSURF,PLANCK_FLAG,WVN,WVNLO,WVNHI,VIEW_FLAG,VIEW_ANGLE,PHI,& QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO,WVN_FLAG,KFLAG,& NUMLAY,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,& PF_CLD,CHI_CLD,PF_AERO,CHI_AERO,ALBEDO,SURF,CHI_SURF,& FOURIER_TOL,FLUXES,WARNING_FLAG,P_IO,R_IO,FLUX_IO,L_IO,& B_IO,INV_IO) !INPUT : SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV, ! TSURF,PLANCK_FLAG,WVN,WVNLO,WVNHI,VIEW_FLAG,VIEW_ANGLE,PHI, ! QUAD,DG,N,DELTA_M_CLD,DELTA_M_AERO,WVN_FLAG,KFLAG, ! NUMLAY,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO, ! TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD, ! PF_CLD,CHI_CLD,PF_AERO,CHI_AERO,ALBEDO,SURF,CHI_SURF, ! FLUXES,WARNING_FLAG,P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO !OUTPUT: IF ALL DATA IS SUITABLE, THEN INSPECTOR DOES NOTHING; ! HOWEVER, IF ANY BAD DATA IS FOUND, IT STOPS PROGRAM EXECUTION. !THIS PROGRAM INSPECTS THE DATA TO BE USED BY SUBROUTINE RAD. !PROGRAMMER: MATT CHRISTI !DATA DICTIONARY**************************************************************** ! ! ALBEDO = SURFACE ALBEDO WHEN SURF=1 ! B_IO = SUBROUTINE BUILD_LAYER I/O FLAG ! (0=OFF, 1=ON) ! CHI_AERO = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_CLD = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED) ! CHI_SURF = VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED ! DELTA_M_AERO = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE AEROSOL PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DELTA_M_CLD = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS ! USED IN SUBROUTINE LOCAL FOR POSSIBLE CLOUD PHASE FUNCTION ! IN EACH LAYER ! (0=NO, 1=YES) ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS, 1=DOUBLE GAUSS) ! FLUX_IO = FLUX TEST I/O FLAG FOR SUBROUTINE RAD ! (0=OFF, 1=ON) ! FLUXES = RADIANT FLUX OUTPUT I/O FLAG ! (0=OFF, 1=ON) ! FSUN = FLUX AT THE TOP OF THE ATMOSPHERE AT A GIVEN WAVELENGTH ! G_AERO = VECTOR OF AEROSOL EFFECTIVE ASYMMETRY PARAMETERS ! G1_AERO = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR FORWARD SCATTERING ! G2_AERO = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING ! G_CLD = VECTOR OF CLOUD EFFECTIVE ASYMMETRY PARAMETERS ! G1_CLD = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR FORWARD SCATTERING ! G2_CLD = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING ! INV_IO = SUBROUTINE INVERT I/O FLAG (FOR MATRIX INVERSE) ! (0=OFF, 1=ON) ! ITMS = VECTOR OF DOWNWELLING DIFFUSE SOLAR RADIANCES FROM THE TOP ! BOUNDARY ! ITMT = VECTOR OF DOWNWELLING DIFFUSE THERMAL RADIANCES FROM THE TOP ! BOUNDARY ! KFLAG = FLAG TO INDICATE WHETHER OR NOT ADDITIONAL COMPUTATIONS WILL ! BE PERFORMED INSIDE RADIANT TO PREPARE FOR SUBSEQUENT CALLS ! TO RADIANT WHERE SAVED MATRICES WILL BE USED ! (0=NO, 1=YES) ! NOTE: KFLAG USED IN CONJUNCTION WITH WVN_FLAG ! L_IO = SUBROUTINE LOCAL I/O FLAG ! (0=OFF, 1=ON) ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! OMEGA_AERO = VECTOR OF SINGLE SCATTER ALBEDO FOR AEROSOL ! OMEGA_CLD = VECTOR OF SINGLE SCATTER ALBEDO FOR CLOUD ! P_IO = SUBROUTINE RADIANT I/O FLAG ! (0=OFF, 1=ON) ! PF_AERO = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_AERO) ! PF_CLD = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE ! FUNCTION EXPANSION ARE OBTAINED ! (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT, ! 2=DETERMINED BY CHI_CLD) ! PHI = AZIMUTH ANGLE FOR WHICH THE RADIANCE VECTOR IS COMPUTED ! (in degrees) ! NOTE: REFERENCE AZIMUTH ASSUMED TO BE 0 ! PLANCK_FLAG = FLAG TO CHOOSE WHETHER SPECTRAL RADIANCES OR RADIANCES WILL ! BE USED FOR THERMAL SOURCES ! (0=SPECTRAL RADIANCES (in mW/(m^2 cm^-1 ster)) ! 1=RADIANCES (in mW/(m^2 ster)) ) ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED ! (0=GAUSS, 1=LOBATTO) ! R_IO = SUBROUTINE RAD I/O FLAG ! (0=OFF, 1=ON) ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES ! (1=SOLAR ONLY, 2=SOLAR & THERMAL, 3=THERMAL ONLY) ! SURF = FLAG INDICATING TYPE OF REFLECTING SURFACE ! (0=NONE, 1=LAMBERTIAN, 2=DETERMINED BY CHI_SURF) ! SZA = THE SOLAR ZENITH ANGLE (in degrees from zenith) ! TAU_AERO = VECTOR OF OPTICAL DEPTHS FOR AEROSOL ! TAU_CLD = VECTOR OF OPTICAL DEPTHS FOR CLOUD ! TAU_GAS = VECTOR OF OPTICAL DEPTHS FOR ATMOSPHERIC GAS ! TAU_RAY = VECTOR OF OPTICAL DEPTHS FOR RAYLEIGH SCATTERING ! TEMISS = EMISSIVITY OF TOP BOUNDARY FOR THERMAL SOURCES ! TLEV = VECTOR (I.E. VERTICAL PROFILE) OF LEVEL-BASED ATMOSPHERIC ! TEMPERATURES (in K) ! TSURF = TEMPERATURE OF THE SURFACE (in K) ! TTOP = TEMPERATURE OF THE TOP BOUNDARY (in K) ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE ! (0=OFF, 1=ON) ! WVN = WAVENUMBER FOR WHICH SPECTRAL RADIANCE IS COMPUTED FOR PLANCK ! FUNCTION (in cm^-1) ! WVNLO = LOW WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS ! COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1) ! WVNHI = HIGH WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS ! COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1) ! WVN_FLAG = FLAG TO INDICATE WHETHER OR NOT SAVED MATRICES WILL BE ! USED TO COMPUTE RADIANCES ON CURRENT CALL TO SUBROUTINE ! RADIANT ! ('NEW'=NO, 'OLD'=YES) ! NOTE: WVN_FLAG USED IN CONJUNCTION WITH KFLAG ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG ! (0=OFF, 1=ON) ! !******************************************************************************* !INTRINSIC SUBPROGRAMS USED BY DATA_INSPECTOR*********************************** ! NONE !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY DATA_INSPECTOR************************************ ! NONE !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & B_IO,DG,FLUX_IO,FLUXES,INV_IO,KFLAG,L_IO,N,& NUMLAY,P_IO,PLANCK_FLAG,QUAD,R_IO,SRC_FLAG,SURF,VIEW_FLAG,& WARNING_FLAG INTEGER, DIMENSION(NUMLAY) :: & DELTA_M_AERO,DELTA_M_CLD,PF_AERO,PF_CLD DOUBLE PRECISION :: & ALBEDO,FSUN,FOURIER_TOL,PHI,SZA,TEMISS,TSURF,& TTOP,VIEW_ANGLE,WVN,WVNLO,WVNHI DOUBLE PRECISION, DIMENSION(N) :: & ITMS,ITMT DOUBLE PRECISION, DIMENSION(NUMLAY) :: & G_AERO,G_CLD,G1_AERO,G1_CLD,G2_AERO,G2_CLD,& TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD DOUBLE PRECISION, DIMENSION(NUMLAY+1) :: & TLEV DOUBLE PRECISION, DIMENSION(0:(4*N)) :: & CHI_SURF DOUBLE PRECISION, DIMENSION(0:(4*N),NUMLAY) :: & CHI_AERO,CHI_CLD CHARACTER(LEN=3) :: & WVN_FLAG !INTERNAL VARIABLES INTEGER :: & I,ERROR_FLAG CHARACTER(LEN=300), DIMENSION(56) :: & PRINT_ERROR LOGICAL, DIMENSION(56) :: & CONDITION !START PROGRAM ! PRINT* ! PRINT*,'ENTERING DATA_INSPECTOR' !TEST SUITABILITY OF INPUT DATA CONDITION(1) = (SRC_FLAG == 1) .OR. (SRC_FLAG == 2) .OR. & (SRC_FLAG == 3) PRINT_ERROR(1) = 'ERROR: SRC_FLAG /= 1, 2, OR 3' CONDITION(2) = FSUN >= 0.0D0 PRINT_ERROR(2) = 'ERROR: FSUN < 0' CONDITION(3) = (SZA >= 0.0D0) .AND. (SZA < 89.99D0) PRINT_ERROR(3) = 'ERROR: SZA < 0 OR SZA >= 89.99' CONDITION(4) = .TRUE. DO I=1,N IF (ITMS(I) < 0.0D0) CONDITION(4) = .FALSE. END DO PRINT_ERROR(4) = 'ERROR: ONE OR MORE OF THE DIFFUSE SOLAR ' // & 'RADIANCES IN ITMS IS < 0' CONDITION(5) = .TRUE. DO I=1,N IF (ITMT(I) < 0.0D0) CONDITION(5) = .FALSE. END DO PRINT_ERROR(5) = 'ERROR: ONE OR MORE OF THE DIFFUSE THERMAL ' // & 'RADIANCES IN ITMT IS < 0' CONDITION(6) = TTOP >= 0.0D0 PRINT_ERROR(6) = 'ERROR: TTOP < 0' CONDITION(7) = (TEMISS >= 0.0D0) .AND. (TEMISS <= 1.0D0) PRINT_ERROR(7) = 'ERROR: TEMISS < 0 OR TEMISS > 1' CONDITION(8) = .TRUE. DO I=1,NUMLAY+1 IF (TLEV(I) < 0.0D0) CONDITION(8) = .FALSE. END DO PRINT_ERROR(8) = 'ERROR: ONE OR MORE COMPONENTS IN TLEV IS < 0' CONDITION(9) = TSURF >= 0.0D0 PRINT_ERROR(9) = 'ERROR: TSURF < 0' CONDITION(10) = (PLANCK_FLAG == 1).OR.(PLANCK_FLAG == 2) PRINT_ERROR(10) = 'ERROR: PLANCK_FLAG /= 1 OR 2' CONDITION(11) = WVN >= 0.0D0 PRINT_ERROR(11) = 'ERROR: WVN < 0' CONDITION(12) = (WVNLO >= 0.0D0) .AND. (WVNLO < WVNHI) PRINT_ERROR(12) = 'ERROR: WVNLO < 0 OR WVNLO >= WVNHI' CONDITION(13) = WVNHI >= 0.0D0 PRINT_ERROR(13) = 'ERROR: WVNHI < 0' CONDITION(14) = (VIEW_FLAG == 0).OR.(VIEW_FLAG == 1) PRINT_ERROR(14) = 'ERROR: VIEW_FLAG /= 0 OR 1' CONDITION(15) = (VIEW_ANGLE >= 0.0D0) .AND. (VIEW_ANGLE < 90.0D0) PRINT_ERROR(15) = 'ERROR: VIEW_ANGLE < 0 OR VIEW_ANGLE >= 90' CONDITION(16) = (PHI >= -180.0D0) .AND. (PHI <= 360.0D0) PRINT_ERROR(16) = 'ERROR: PHI < -180 OR PHI > 360' CONDITION(17) = (QUAD == 0).OR.(QUAD == 1) PRINT_ERROR(17) = 'ERROR: QUAD /= 0 OR 1' CONDITION(18) = (DG == 0).OR.(DG == 1) PRINT_ERROR(18) = 'ERROR: DG /= 0 OR 1' CONDITION(19) = N >= 3 PRINT_ERROR(19) = 'ERROR: N < 3' CONDITION(20) = FOURIER_TOL >= 0.0D0 PRINT_ERROR(20) = 'ERROR: FOURIER_TOL < 0' CONDITION(21) = .TRUE. DO I=1,NUMLAY IF (.NOT.(DELTA_M_CLD(I) == 0) .AND. .NOT.(DELTA_M_CLD(I) == 1)) & CONDITION(21) = .FALSE. END DO PRINT_ERROR(21) = 'ERROR: ONE OR MORE COMPONENTS IN DELTA_M_CLD ' // & 'IS /= 0 OR 1' CONDITION(22) = .TRUE. DO I=1,NUMLAY IF (.NOT.(DELTA_M_AERO(I) == 0) .AND. .NOT.(DELTA_M_AERO(I) == 1)) & CONDITION(22) = .FALSE. END DO PRINT_ERROR(22) = 'ERROR: ONE OR MORE COMPONENTS IN DELTA_M_AERO ' // & 'IS /= 0 OR 1' CONDITION(23) = (WVN_FLAG == 'NEW').OR.(WVN_FLAG == 'OLD') PRINT_ERROR(23) = 'ERROR: WVN_FLAG /= NEW OR OLD' CONDITION(24) = (KFLAG == 0).OR.(KFLAG == 1) PRINT_ERROR(24) = 'ERROR: KFLAG /= 0 OR 1' CONDITION(25) = (NUMLAY >= 1).AND.(NUMLAY <= 100) PRINT_ERROR(25) = 'ERROR: NUMLAY < 1 OR NUMLAY > 100' !NOT USED NOW CONDITION(26) = .TRUE. PRINT_ERROR(26) = '' CONDITION(27) = .TRUE. DO I=1,NUMLAY IF ((G2_CLD(I) < -0.9999999999D0) .OR. (G2_CLD(I) > 0.0D0)) & CONDITION(27) = .FALSE. END DO PRINT_ERROR(27) = 'ERROR: ONE OR MORE COMPONENTS IN G2_CLD ARE ' // & '<= -1 OR > 0' CONDITION(28) = .TRUE. DO I=1,NUMLAY IF (G2_CLD(I) < -1.0D-12) THEN !DOUBLE HENYEY-GREENSTEIN ASSUMED ACTIVE IF ((G1_CLD(I) < 0.0D0) .OR. (G1_CLD(I) > 0.9999999999D0)) & CONDITION(28) = .FALSE. ELSE !NORMAL HENYEY-GREENSTEIN ASSUME ACTIVE IF ((G1_CLD(I) < -0.9999999999D0) .OR. & (G1_CLD(I) > 0.9999999999D0)) & CONDITION(28) = .FALSE. END IF END DO PRINT_ERROR(28) = 'ERROR: ONE OR MORE COMPONENTS IN G1_CLD ARE ' // & 'OUT OF RANGE (EACH COMPONENT OF G1_CLD ' // & 'SHOULD BE BETWEEN 0 AND 0.9999999999 ' // & '(INCLUSIVE) FOR DOUBLE HENYEY-GREENSTEIN OR ' // & 'BETWEEN -0.9999999999 AND 0.9999999999 ' // & '(INCLUSIVE) FOR NORMAL HENYEY-GREENSTEIN' CONDITION(29) = .TRUE. DO I=1,NUMLAY IF ((G_CLD(I) < -0.9999999999D0) .OR. (G_CLD(I) > 0.9999999999D0)) & CONDITION(29) = .FALSE. END DO PRINT_ERROR(29) = 'ERROR: ONE OR MORE COMPONENTS IN G_CLD ARE ' // & '<= -1 OR >= 1' CONDITION(30) = .TRUE. DO I=1,NUMLAY IF ((G2_AERO(I) < -0.9999999999D0) .OR. (G2_AERO(I) > 0.0D0)) & CONDITION(30) = .FALSE. END DO PRINT_ERROR(30) = 'ERROR: ONE OR MORE COMPONENTS IN G2_AERO ARE ' // & '<= -1 OR > 0' CONDITION(31) = .TRUE. DO I=1,NUMLAY IF (G2_AERO(I) < -1.0D-12) THEN !DOUBLE HENYEY-GREENSTEIN ASSUMED ACTIVE IF ((G1_AERO(I) < 0.0D0) .OR. (G1_AERO(I) > 0.9999999999D0)) & CONDITION(31) = .FALSE. ELSE !NORMAL HENYEY-GREENSTEIN ASSUME ACTIVE IF ((G1_AERO(I) < -0.9999999999D0) .OR. & (G1_AERO(I) > 0.9999999999D0)) & CONDITION(31) = .FALSE. END IF END DO PRINT_ERROR(31) = 'ERROR: ONE OR MORE COMPONENTS IN G1_AERO ARE ' // & 'OUT OF RANGE (EACH COMPONENT OF G1_AERO ' // & 'SHOULD BE BETWEEN 0 AND 0.9999999999 ' // & '(INCLUSIVE) FOR DOUBLE HENYEY-GREENSTEIN OR ' // & 'BETWEEN -0.9999999999 AND 0.9999999999 ' // & '(INCLUSIVE) FOR NORMAL HENYEY-GREENSTEIN' CONDITION(32) = .TRUE. DO I=1,NUMLAY IF ((G_AERO(I) < -0.9999999999D0) .OR. (G_AERO(I) > 0.9999999999D0)) & CONDITION(32) = .FALSE. END DO PRINT_ERROR(32) = 'ERROR: ONE OR MORE COMPONENTS IN G_AERO ARE ' // & '<= -1 OR >= 1' CONDITION(33) = .TRUE. DO I=1,NUMLAY IF ((OMEGA_CLD(I) < 0.0D0).OR.(OMEGA_CLD(I) > 1.0D0)) & CONDITION(33) = .FALSE. END DO PRINT_ERROR(33) = 'ERROR: ONE OR MORE COMPONENTS IN OMEGA_CLD IS ' // & 'EITHER < 0 OR > 1' CONDITION(34) = .TRUE. DO I=1,NUMLAY IF (TAU_CLD(I) < 0.0D0) CONDITION(34) = .FALSE. END DO PRINT_ERROR(34) = 'ERROR: ONE OR MORE COMPONENTS IN TAU_CLD IS < 0' !NOT USED NOW CONDITION(35) = .TRUE. PRINT_ERROR(35) = '' CONDITION(36) = .TRUE. DO I=1,NUMLAY IF ((OMEGA_AERO(I) < 0.0D0).OR.(OMEGA_AERO(I) > 1.0D0)) & CONDITION(36) = .FALSE. END DO PRINT_ERROR(36) = 'ERROR: ONE OR MORE COMPONENTS IN OMEGA_AERO IS ' // & 'EITHER < 0 OR > 1' CONDITION(37) = .TRUE. DO I=1,NUMLAY IF (TAU_AERO(I) < 0.0D0) CONDITION(37) = .FALSE. END DO PRINT_ERROR(37) = 'ERROR: ONE OR MORE COMPONENTS IN TAU_AERO IS < 0' !NOT USED NOW CONDITION(38) = .TRUE. PRINT_ERROR(38) = '' CONDITION(39) = .TRUE. DO I=1,NUMLAY IF (TAU_RAY(I) < 0.0D0) CONDITION(39) = .FALSE. END DO PRINT_ERROR(39) = 'ERROR: ONE OR MORE COMPONENTS IN TAU_RAY IS < 0' CONDITION(40) = .TRUE. DO I=1,NUMLAY IF (TAU_GAS(I) < 0.0D0) CONDITION(40) = .FALSE. END DO PRINT_ERROR(40) = 'ERROR: ONE OR MORE COMPONENTS IN TAU_GAS IS < 0' CONDITION(41) = .TRUE. DO I=1,NUMLAY IF (.NOT.(PF_CLD(I) == 1) .AND. .NOT.(PF_CLD(I) == 2)) & CONDITION(41) = .FALSE. END DO PRINT_ERROR(41) = 'ERROR: ONE OR MORE COMPONENTS IN PF_CLD IS /= 1 OR 2' CONDITION(42) = .TRUE. DO I=1,NUMLAY IF (PF_CLD(I) == 2) THEN IF (.NOT.(CHI_CLD(0,I) >= 0.9999999999D0) .AND. & .NOT.(CHI_CLD(0,I) <= 1.0000000001D0)) & CONDITION(42) = .FALSE. END IF END DO PRINT_ERROR(42) = 'ERROR: ONE OR MORE CLOUD PHASE FUNCTIONS ' // & 'HAVE A FIRST COMPONENT THAT IS /= 1' CONDITION(43) = .TRUE. DO I=1,NUMLAY IF (.NOT.(PF_AERO(I) == 1) .AND. .NOT.(PF_AERO(I) == 2)) & CONDITION(43) = .FALSE. END DO PRINT_ERROR(43) = 'ERROR: ONE OR MORE COMPONENTS IN PF_AERO IS ' // & '= 1 OR 2' CONDITION(44) = .TRUE. DO I=1,NUMLAY IF (PF_AERO(I) == 2) THEN IF (.NOT.(CHI_AERO(0,I) >= 0.9999999999D0) .AND. & .NOT.(CHI_AERO(0,I) <= 1.0000000001D0)) & CONDITION(44) = .FALSE. END IF END DO PRINT_ERROR(44) = 'ERROR: ONE OR MORE AEROSOL PHASE FUNCTIONS ' // & 'HAVE A FIRST COMPONENT THAT IS /= 1' CONDITION(45) = (ALBEDO >= 0.0D0) .AND. (ALBEDO <= 1.0D0) PRINT_ERROR(45) = 'ERROR: ALBEDO < 0 OR ALBEDO > 1' CONDITION(46) = (SURF == 0).OR.(SURF == 1).OR.(SURF == 2) PRINT_ERROR(46) = 'ERROR: SURF /= 0, 1, OR 2' CONDITION(47) = .TRUE. IF (SURF == 2) & CONDITION(47) = (CHI_SURF(0) >= 0.0D0) .AND.& (CHI_SURF(0) <= 1.0000000001D0) PRINT_ERROR(47) = 'ERROR: SURFACE REFLECTION FUNCTION ' // & 'HAS A FIRST COMPONENT THAT IS < 0 OR > 1' CONDITION(48) = (FLUXES == 0).OR.(FLUXES == 1) PRINT_ERROR(48) = 'ERROR: FLUXES /= 0 OR 1' CONDITION(49) = (WARNING_FLAG == 0).OR.(WARNING_FLAG == 1) PRINT_ERROR(49) = 'ERROR: WARNING_FLAG /= 0 OR 1' CONDITION(50) = (P_IO == 0).OR.(P_IO == 1) PRINT_ERROR(50) = 'ERROR: P_IO /= 0 OR 1' CONDITION(51) = (R_IO == 0).OR.(R_IO == 1) PRINT_ERROR(51) = 'ERROR: R_IO /= 0 OR 1' CONDITION(52) = (FLUX_IO == 0).OR.(FLUX_IO == 1) PRINT_ERROR(52) = 'ERROR: FLUX_IO /= 0 OR 1' CONDITION(53) = (L_IO == 0).OR.(L_IO == 1) PRINT_ERROR(53) = 'ERROR: L_IO /= 0 OR 1' CONDITION(54) = (B_IO == 0).OR.(B_IO == 1) PRINT_ERROR(54) = 'ERROR: B_IO /= 0 OR 1' CONDITION(55) = (INV_IO == 0).OR.(INV_IO == 1) PRINT_ERROR(55) = 'ERROR: INV_IO /= 0 OR 1' CONDITION(56) = .TRUE. IF ((QUAD == 1).AND.(DG == 1)) CONDITION(56) = .FALSE. PRINT_ERROR(56) = 'ERROR: DG SHOULD BE 0 WHEN QUAD IS 1' ERROR_FLAG = 0 DO I=1,56 IF (.NOT. CONDITION(I)) THEN PRINT* PRINT*,'I = ',I,'CONDITION(I) = ',CONDITION(I) PRINT*, PRINT_ERROR(I) ERROR_FLAG = 1 ENDIF END DO IF (ERROR_FLAG == 1) THEN PRINT*,'ERROR FLAG WAS SET' STOP END IF !END PROGRAM ! PRINT* ! PRINT*,'LEAVING DATA_INSPECTOR' END SUBROUTINE DATA_INSPECTOR !******************************************************************************* !******************************************************************************* DOUBLE PRECISION FUNCTION DPLKAVG(WNUMLO,WNUMHI,T) ! Computes Planck function integrated between two wavenumbers ! ! INPUT : WNUMLO : Lower wavenumber (inv cm) of spectral interval ! ! WNUMHI : Upper wavenumber ! ! T : Temperature (K) ! ! NOTE: These stated units lack "sr" (MJC) ! (actually does | ! Watts/(sq m * sr)) | ! | ! OUTPUT : DPLKAVG : Integrated Planck function ( Watts/sq m )<- ! = Integral (WNUMLO to WNUMHI) of ! 2h c**2 nu**3 / ( EXP(hc nu/kT) - 1) ! (where h=Plancks constant, c=speed of ! light, nu=wavenumber, T=temperature, ! and k = Boltzmann constant) ! ! Reference : Specifications of the Physical World: New Value ! of the Fundamental Constants, Dimensions/N.B.S., ! Jan. 1974 ! ! Method : For WNUMLO close to WNUMHI, a Simpson-rule quadrature ! is done to avoid ill-conditioning; otherwise ! ! (1) For WNUMLO or WNUMHI small, ! integral(0 to WNUMLO/HI) is calculated by expanding ! the integrand in a power series and integrating ! term by term; ! ! (2) Otherwise, integral(WNUMLO/HI to INFINITY) is ! calculated by expanding the denominator of the ! integrand in powers of the exponential and ! integrating term by term. ! ! Accuracy : At least 6 significant digits, assuming the ! physical constants are infinitely accurate ! ! ERRORS WHICH ARE NOT TRAPPED: ! ! * power or exponential series may underflow, giving no ! significant digits. This may or may not be of concern, ! depending on the application. ! ! * Simpson-rule special case is skipped when denominator of ! integrand will cause overflow. In that case the normal ! procedure is used, which may be inaccurate if the ! wavenumber limits (WNUMLO, WNUMHI) are close together. ! ! LOCAL VARIABLES ! ! A1,2,... : Power series coefficients ! C2 : h * c / k, in units cm*K (h = Plancks constant, ! c = speed of light, k = Boltzmann constant) ! D(I) : Exponential series expansion of integral of ! Planck function from WNUMLO (i=1) or WNUMHI ! (i=2) to infinity ! EPSIL : SMALLEST NUMBER SUCH THAT 1+EPSIL .GT. 1 on ! computer ! EX : EXP( - V(I) ) ! EXM : EX**M ! MMAX : No. of terms to take in exponential series ! MV : Multiples of V(I) ! P(I) : Power series expansion of integral of ! Planck function from zero to WNUMLO (I=1) or ! WNUMHI (I=2) ! PI : 3.14159... ! SIGMA : Stefan-Boltzmann constant (W/m**2/K**4) ! SIGDPI : SIGMA / PI ! SMALLV : Number of times the power series is used (0,1,2) ! V(I) : C2 * (WNUMLO(I=1) or WNUMHI(I=2)) / temperature ! VCUT : Power-series cutoff point ! VCP : Exponential series cutoff points ! VMAX : Largest allowable argument of EXP function ! ! Called by- DISORT ! Calls- D1MACH, ERRMSG ! ! ! NOTE: ORIGINAL DISORT SUBROUTINE CONVERTED TO A MORE F90 FORMAT ! AND DOUBLE PRECISION. PHYSICAL CONSTANTS ALSO UPDATED. ! BY : M.J. CHRISTI ! DATE: 4/3/03 ! ! ---------------------------------------------------------------------- ! .. Parameters .. DOUBLE PRECISION A1,A2,A3,A4,A5,A6 PARAMETER (A1 = 1.0D0/3.0D0, A2 = -1.0D0/8.0D0, A3 = 1.0D0/60.0D0,& A4 = -1.0D0/5040.0D0, A5 = 1.0D0/272160.0D0,& A6 = -1. / 13305600. ) ! .. Scalar Arguments .. DOUBLE PRECISION T,WNUMHI,WNUMLO ! .. Local Scalars .. INTEGER I,K,M,MMAX,N,SMALLV DOUBLE PRECISION C2,CONC,DEL,EPSIL,EX,EXM,HH,MV,OLDVAL,PI,& SIGDPI,SIGMA,VAL,VAL0,VCUT,VMAX,VSQ,X ! .. Local Arrays .. DOUBLE PRECISION D(2),P(2),V(2),VCP(7) ! .. External Functions .. DOUBLE PRECISION D1MACH EXTERNAL D1MACH ! .. External Subroutines .. EXTERNAL ERRMSG ! .. Intrinsic Functions .. INTRINSIC DABS,DASIN,DEXP,DLOG,DMOD ! .. Statement Functions .. DOUBLE PRECISION PLKF SAVE PI,CONC,VMAX,EPSIL,SIGDPI DATA C2/1.4387752D0/, SIGMA/5.6704D-8/, VCUT/1.5D0/,& VCP/10.25D0,5.7D0,3.9D0,2.9D0,2.3D0,1.9D0,0.0D0/ DATA PI/0.0D0/ ! .. Statement Function definitions .. PLKF(X) = X**3/(DEXP(X) - 1) !**** Begin Executable **** !**** Preliminary Checks **** IF(PI == 0.0D0) THEN PI = 2.0D0*DASIN(1.0D0) VMAX = DLOG(D1MACH(2)) EPSIL = D1MACH(4) SIGDPI = SIGMA/PI CONC = 15.0D0/(PI**4) END IF IF((T < 0.0D0) .OR. (WNUMHI <= WNUMLO) .OR. (WNUMLO < 0.0D0)) & CALL ERRMSG('DPLKAVG--temperature or wavenums. wrong',.TRUE.) !**** Special Cases **** IF(T < 1.0D-4) THEN ! ** Temperature near 0 K DPLKAVG = 0.0D0 RETURN END IF V(1) = C2*WNUMLO/T V(2) = C2*WNUMHI/T IF((V(1) > EPSIL) .AND. (V(2) < VMAX) .AND. & (((WNUMHI - WNUMLO)/WNUMHI) < 1.0D-2)) THEN ! ** Wavenumbers are very close. Get integral ! ** by iterating Simpson rule to convergence. HH = V(2) - V(1) OLDVAL = 0.0D0 VAL0 = PLKF(V(1)) + PLKF(V(2)) DO 20 N=1,10 DEL = HH/(2.0D0*DBLE(N)) VAL = VAL0 DO 10 K=1,2*N-1 VAL = VAL + 2.0D0*(1.0D0 + & DMOD(DBLE(K),2.0D0))*PLKF(V(1) + DBLE(K)*DEL) 10 CONTINUE VAL = DEL/3.0D0*VAL IF(DABS((VAL - OLDVAL)/VAL) <= 1.0D-8) GO TO 30 OLDVAL = VAL 20 CONTINUE CALL ERRMSG( 'DPLKAVG--Simpson rule didnt converge',.FALSE.) 30 CONTINUE DPLKAVG = SIGDPI*(T**4)*CONC*VAL RETURN END IF !**** General case **** SMALLV = 0 DO 60 I=1,2 IF(V(I) < VCUT) THEN ! ** Use power series SMALLV = SMALLV + 1 VSQ = V(I)**2 P(I) = CONC*VSQ*V(I)*(A1 + V(I)*(A2 + V(I)*(A3 + & VSQ*(A4 + VSQ*(A5 + VSQ*A6))))) ELSE ! ** Use exponential series MMAX = 0 ! ** Find upper limit of series 40 CONTINUE MMAX = MMAX + 1 IF(V(I) < VCP(MMAX)) GO TO 40 EX = DEXP(-V(I)) EXM = 1.0D0 D(I) = 0.0D0 DO 50 M=1,MMAX MV = M*V(I) EXM = EX*EXM D(I) = D(I) + EXM*(6.0D0 + MV*(6.0D0 + MV*(3.0D0 + MV)))/M**4 50 CONTINUE D(I) = CONC*D(I) END IF 60 CONTINUE !**** Handle ill-conditioning **** IF(SMALLV == 2) THEN ! ** WNUMLO and WNUMHI both small DPLKAVG = P(2) - P(1) ELSE IF(SMALLV == 1) THEN ! ** WNUMLO small, WNUMHI large DPLKAVG = 1.0D0 - P(1) - D(2) ELSE ! ** WNUMLO and WNUMHI both large DPLKAVG = D(1) - D(2) END IF DPLKAVG = SIGDPI * T**4 * DPLKAVG IF(DPLKAVG == 0.0D0) & CALL ERRMSG('DPLKAVG--returns zero; possible underflow', & .FALSE.) END FUNCTION DPLKAVG !******************************************************************************* !******************************************************************************* SUBROUTINE GET_GLOBAL_SOURCES5 (SRC_FLAG,N,DEGREE,MU0,OMEGA,& TAU,FTOP,B_T,RECMU,PSOLP,PSOLM,T,R,GT,GR,L_IO,INV_IO,& GSPS,GSMS,GSPT,GSMT) !INPUT : SRC_FLAG,N,DEGREE,MU0,OMEGA,TAU,FTOP,B_T,RECMU,PSOLP,PSOLM,T,R, ! GT,GR,L_IO,INV_IO !OUTPUT: GSPS,GSMS,GSPT,GSMT !THIS PROGRAM OBTAINS THE GLOBAL SOURCE VECTORS FOR THE LAYER BEING BUILT !PROGRAMMER: MATT CHRISTI !DATE: 8/9/04 !DATA DICTIONARY******************************************************** ! ! A = THE A MATRIX ! B_T = VECTOR OF SPECTRAL RADIANCES OR RADIANCES (DEPENDING ON VALUE ! OF PLANCK_FLAG) USING THE PLANCK FUNCTION OF EMISSION ! BC = MATRIX OF CONSTANTS DERIVED FROM USING THE PLANCK FUNCTION ! OF EMISSION ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! FTOP = FLUX AT THE TOP OF THE LAYER AT A GIVEN WAVENUMBER ! GR = GLOBAL REFLECTION MATRIX ! GSPS = GLOBAL SOURCE VECTOR PLUS SOLAR (UPWELLING) ! GSMS = GLOBAL SOURCE VECTOR MINUS SOLAR (DOWNWELLING) ! GSPT = GLOBAL SOURCE VECTOR PLUS THERMAL (UPWELLING) ! GSMT = GLOBAL SOURCE VECTOR MINUS THERMAL (DOWNWELLING) ! GT = GLOBAL TRANSMISSION MATRIX ! IMUA_I = INVERSE OF MATRIX (I/MU0 - A) ! INV_IO = SUBROUTINE INVERT I/O FLAG (FOR MATRIX INVERSE) ! (0=OFF, 1=ON) ! L_IO = SUBROUTINE LOCAL I/O FLAG (0=OFF, 1=ON) ! MPM = MATRIX PRODUCT OF MINV*PSOLM ! MPP = MATRIX PRODUCT OF MINV*PSOLP ! MU0 = COSINE OF SOLAR ZENITH ANGLE ! OMEGA = SINGLE SCATTER ALBEDO ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTION ! NS = TOTAL NUMBER OF STREAMS ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! PSOLM = SOLAR SOURCE VECTOR CONSISTING OF COMPONENTS FROM +MU ! DIRECTIONS ! PSOLP = SOLAR SOURCE VECTOR CONSISTING OF COMPONENTS FROM -MU ! DIRECTIONS ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! SUBROUTINE LOCAL (0=GAUSS,1=LOBATTO) ! R = LOCAL REFLECTION MATRIX ! RECMU = VECTOR HOLDING RECIPROCALS OF THE QUADRATURE ROOTS ! S1M = S2M*DECAYING EXPONENTIAL ! S1P = S2P*DECAYING EXPONENTIAL ! S2M = MATRIX-VECTOR PRODUCT OF IMUA_I*SIGMAM ! S2P = MATRIX-VECTOR PRODUCT OF IMUA_I*SIGMAP ! S3M = MATRIX-VECTOR PRODUCT OF (-1)*(THERM_VEC1 + THERM_VEC2) ! S3P = MATRIX-VECTOR PRODUCT OF THERM_VEC1 + THERM_VEC2 ! S4M = VECTOR SUM OF S3M - (TAU*B_T(1))*THERM_VEC1 ! S4P = VECTOR SUM OF S3P - (TAU*B_T(1))*THERM_VEC1 ! SIGMAM = LOCAL DOWNWELLING SOLAR SOURCE ! SIGMAP = LOCAL UPWELLING SOLAR SOURCE ! SRC_FLAG = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING ! RADIANCES (1=SOLAR ONLY,2=BOTH,3=THERMAL ONLY) ! T = LOCAL TRANSMISSION MATRIX ! TAU = TOTAL OPTICAL DEPTH OF THE LAYER ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! X = VECTOR HOLDING COEFFICIENTS OF COMPOSITE PHASE FUNCTION ! !*********************************************************************** !INTRINSIC SUBPROGRAMS USED BY GET_GLOBAL_SOURCES5********************** ! DEXP,MATMUL !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY GET_GLOBAL_SOURCES5*********************** ! MATIDENT,MM_IG1G2,MV_DV !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N,DEGREE,L_IO,SRC_FLAG,INV_IO DOUBLE PRECISION :: & MU0,FTOP,TAU,OMEGA DOUBLE PRECISION, DIMENSION(0:1) :: & B_T DOUBLE PRECISION, DIMENSION(N) :: & PSOLM,PSOLP,RECMU DOUBLE PRECISION, DIMENSION(N,N) :: & GT,GR,T,R !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & GSPS,GSMS,GSPT,GSMT !INTERNAL VARIABLES INTEGER :: & I,J,NS DOUBLE PRECISION :: & CONSTANT2 DOUBLE PRECISION, DIMENSION(N) :: & MPM,MPP,SIGMAP,SIGMAM,& S1P,S1M,S2P,S2M,S3P,S3M,S4P,S4M DOUBLE PRECISION, DIMENSION(2*N) :: & SIGMA,S1,S2,THERM_VEC,THERM_VEC1,THERM_VEC2,S3,S4 DOUBLE PRECISION, DIMENSION(0:1) :: & BC DOUBLE PRECISION, DIMENSION(2*N,2*N) :: & A,IMU,IMUA,IMUA_I DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 !START PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'ENTERING GET_GLOBAL_SOURCES5' END IF NS=2*N !CONSTRUCT "A" MATRIX TO COMPUTE SOURCES DO I=1,N DO J=1,N A(I,J) = T(I,J) A(I,N+J) = -1*R(I,J) A(N+I,J) = R(I,J) A(N+I,N+J) = -1*T(I,J) END DO END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE A MATRIX LOOKS LIKE:' WRITE(*,50) ((A(I,J),J=1,NS),I=1,NS) 50 FORMAT(16(1X,F16.8)) END IF !COMPUTE LOCAL SOLAR SOURCES IF NECESSARY IF ((SRC_FLAG == 1).OR.(SRC_FLAG == 2)) THEN !CONSTRUCT SIGMA+ AND SIGMA- TO COMPUTE SOLAR SOURCES CALL MV_DV(N,RECMU,PSOLM,MPM) CALL MV_DV(N,RECMU,PSOLP,MPP) CONSTANT2 = (OMEGA*FTOP)/(4.0D0*PI) DO I=1,N SIGMAP(I) = CONSTANT2*MPM(I) SIGMA(I) = SIGMAP(I) SIGMAM(I) = -1.0D0*CONSTANT2*MPP(I) SIGMA(N+I) = SIGMAM(I) END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE SIGMA+ SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (SIGMAP(I),I=1,N) PRINT* PRINT*,'THE SIGMA- SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (SIGMAM(I),I=1,N) END IF !CONSTRUCT I/MU_NOT MATRIX CALL MATIDENT(NS,IMU) DO I=1,N IMU(I,I) = 1.0D0/MU0 IMU(N+I,N+I) = 1.0D0/MU0 END DO IMUA = IMU - A !CONSTRUCT SOLAR SOURCES CALL MM_IG1G2(NS,1,IMUA,SIGMA,S2,INV_IO) DO I=1,NS S1(I) = S2(I)*DEXP(-1.0D0*TAU/MU0) END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE S2 SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S2(I),I=1,NS) PRINT* PRINT*,'THE S1 SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S1(I),I=1,NS) END IF DO I=1,N S1P(I) = S1(I) S1M(I) = S1(N+I) S2P(I) = S2(I) S2M(I) = S2(N+I) END DO ELSE S1P = 0.0D0 S1M = 0.0D0 S2P = 0.0D0 S2M = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE S1P SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S1P(I),I=1,N) PRINT* PRINT*,'THE S1M SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S1M(I),I=1,N) PRINT* PRINT*,'THE S2P SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S2P(I),I=1,N) PRINT* PRINT*,'THE S2M SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S2M(I),I=1,N) END IF !COMPUTE LOCAL THERMAL SOURCES IF NECESSARY IF (((SRC_FLAG == 2).OR.(SRC_FLAG == 3)).AND.(DEGREE == 0)) THEN !CONSTRUCT THERMAL SOURCES DO I=1,N THERM_VEC(I) = (1.0D0-OMEGA)*RECMU(I) THERM_VEC(N+I) = -1.0D0*THERM_VEC(I) END DO CALL MM_IG1G2(NS,1,A,THERM_VEC,THERM_VEC1,INV_IO) CALL MM_IG1G2(NS,1,A,THERM_VEC1,THERM_VEC2,INV_IO) BC(0) = B_T(1) BC(1) = (B_T(0) - B_T(1))/TAU IF (L_IO == 1) THEN PRINT* WRITE(*,45) 'B_T(0) = ',B_T(0),'B_T(1) = ',B_T(1) PRINT* WRITE(*,45) ' BC(0) = ',BC(0),' BC(1) = ',BC(1) 45 FORMAT(2(2X,A9,2X,E15.9E2)) END IF S3 = BC(0)*THERM_VEC1 + BC(1)*THERM_VEC2 S4 = S3 + (BC(1)*TAU)*THERM_VEC1 DO I=1,N S3P(I) = S3(I) S3M(I) = S3(N+I) S4P(I) = S4(I) S4M(I) = S4(N+I) END DO ELSE S3P = 0.0D0 S3M = 0.0D0 S4P = 0.0D0 S4M = 0.0D0 END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE S3P SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S3P(I),I=1,N) PRINT* PRINT*,'THE S3M SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S3M(I),I=1,N) PRINT* PRINT*,'THE S4P SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S4P(I),I=1,N) PRINT* PRINT*,'THE S4M SOURCE VECTOR LOOKS LIKE:' WRITE(*,20) (S4M(I),I=1,N) END IF !COMPUTE GLOBAL SOURCES !COMPUTE GLOBAL SOLAR SOURCE VECTORS IF ((SRC_FLAG == 1).OR.(SRC_FLAG == 2)) THEN GSPS = S2P - MATMUL(GT,S1P) - MATMUL(GR,S2M) GSMS = S1M - MATMUL(GR,S1P) - MATMUL(GT,S2M) ELSE GSPS = 0.0D0 GSMS = 0.0D0 END IF ! print*, GR,'gr' !COMPUTE GLOBAL THERMAL SOURCE VECTORS IF (((SRC_FLAG == 2).OR.(SRC_FLAG == 3)).AND.(DEGREE == 0)) THEN GSPT = MATMUL(GT,S3P) + MATMUL(GR,S4M) - S4P GSMT = MATMUL(GR,S3P) + MATMUL(GT,S4M) - S3M ELSE GSPT = 0.0D0 GSMT = 0.0D0 END IF !DISPLAY GLOBAL SOURCE VECTORS IF (L_IO == 1) THEN PRINT* PRINT*,'THE GSPS VECTOR LOOKS LIKE:' WRITE(*,20) (GSPS(I),I=1,N) PRINT* PRINT*,'THE GSMS VECTOR LOOKS LIKE:' WRITE(*,20) (GSMS(I),I=1,N) PRINT* PRINT*,'THE GSPT VECTOR LOOKS LIKE:' WRITE(*,20) (GSPT(I),I=1,N) PRINT* PRINT*,'THE GSMT VECTOR LOOKS LIKE:' WRITE(*,20) (GSMT(I),I=1,N) END IF 20 FORMAT(1X,F16.8) !END PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'LEAVING GET_GLOBAL_SOURCES5' END IF END SUBROUTINE GET_GLOBAL_SOURCES5 !********************************************************************** !********************************************************************** SUBROUTINE GET_GT_GR4 (N,LDA,LDVL,LDVR,LWORK,VIEW_FLAG,VIEW_ANGLE,& QUAD,DG,NUMDEG,DEGREE,NEXP,DELTA_0_M,MU0,LAYER,NUMLAY,& WVN_FLAG,KFLAG,PERT_FLAG,EPS,OMEGA,TAU,X,L_IO,INV_IO,& DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,RECMU,PSOLP,PSOLM,T,R,GT,GR) !INPUT : N,LDA,LDVL,LDVR,LWORK,VIEW_FLAG,VIEW_ANGLE,QUAD,DG,NUMDEG,DEGREE, ! NEXP,DELTA_0_M,MU0,LAYER,NUMLAY,WVN_FLAG,KFLAG,PERT_FLAG,EPS, ! OMEGA,TAU,X,L_IO,INV_IO !OUTPUT: RECMU,PSOLP,PSOLM,T,R,GT,GR !THIS PROGRAM OBTAINS THE GLOBAL TRANSMISSION AND REFLECTION MATRICES !FOR THE LAYER BEING BUILT !PROGRAMMER: MATT CHRISTI !DATE: 5/5/04 !DATA DICTIONARY******************************************************** ! ! B = MATRIX PRODUCT OF (T-R)*(T+R) ! B_IO = SUBROUTINE BUILD_LAYER2 I/O FLAG (0=OFF, 1=ON) ! D = MATRIX OF EIGENVECTORS OF B ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! E = IDENTITY MATRIX ! EPS = AMOUNT BY WHICH OPTICAL PROPERTIES ARE PERTURBED BY FOR THE ! COMPUTATION OF JACOBIAN ELEMENTS BY FINITE DIFFERENCE WHEN ! THE PERTURBATION SCHEME IS IN USE ! GR = GLOBAL REFLECTION MATRIX ! GT = GLOBAL TRANSMISSION MATRIX ! INV_IO = SUBROUTINE MM_IG1G2 I/O FLAG (0=OFF, 1=ON) ! (FOR MATRIX MULTIPLICATION A^(-1)*B) ! L_IO = SUBROUTINE LOCAL I/O FLAG (0=OFF, 1=ON) ! LAMBP = VECTOR OF POSITIVE EIGENVALUES OF MATRIX A ! LAMBM = VECTOR OF NEGATIVE EIGENVALUES OF MATRIX A ! LAYER = LAYER WHOSE OPTICAL PROPERTIES ARE BEING COMPUTED ! MU0 = COSINE OF SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTION ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! PERT_FLAG = FLAG TO INDICATE WHETHER OR NOT THE PERTURBATION SCHEME WILL ! BE USED TO BUILD LAYERS ON SUBSEQUENT CALLS TO RADIANT DURING ! LAYER SAVING MODE (0=NO, 1=YES) ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! COMPUTING THE VALUES OF MU AND W (0=GAUSS,1=LOBATTO) ! R = LOCAL REFLECTION MATRIX ! T = LOCAL TRANSMISSION MATRIX ! TAU = TOTAL OPTICAL DEPTH OF THE LAYER ! TMR = MATRIX HOLDING DIFFERENCE OF T AND R MATRICES (T-R) ! TPR = MATRIX HOLDING SUM OF T AND R MATRICES (T+R) ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! OMEGA = SINGLE SCATTER ALBEDO ! X = VECTOR HOLDING COEFFICIENTS OF COMPOSITE PHASE FUNCTION ! !*********************************************************************** !INTRINSIC SUBPROGRAMS USED BY GET_GT_GR4******************************* ! DEXP,DSQRT,MATMUL,TRANSPOSE !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY GET_GT_GR4******************************** ! INTERMEDIATE_RESULTS,MATDIAG,MATIDENT,MM_IG1G2,SGEEVX !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,NEXP,DEGREE,LAYER,NUMDEG,L_IO,INV_IO,VIEW_FLAG,& DELTA_0_M,NUMLAY DOUBLE PRECISION :: & MU0,OMEGA,TAU,VIEW_ANGLE DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & X LOGICAL :: & DIMENSIONS_RESET,SZA_RESET,QUAD_RESET !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & RECMU,PSOLP,PSOLM DOUBLE PRECISION, DIMENSION(N,N) :: & T,R,GT,GR !INTERNAL VARIABLES INTEGER :: & I,J,K DOUBLE PRECISION, DIMENSION(N) :: & LAMBP,LAMBM,TEMP3,TEMP4 DOUBLE PRECISION, DIMENSION(N,N) :: & B,D,E,TMR,TPR,GRTEMP,GTTEMP,EXPO,XP,XM,UM,UP,& UMI,UMIUP,UPUMI,MAT1,MAT2,MAT3,MAT4 !FOR SUBROUTINE SGEEVX (FOR COMPUTING EIGENVALUES AND EIGENVECTORS): !INPUT VARIABLES INTEGER :: & LDA,LDVL,LDVR,LWORK !INTERNAL VARIABLES INTEGER :: & IWORK(2*N-2) INTEGER :: & IHI,ILO,INFO DOUBLE PRECISION :: & ABNRM,WR(N),WI(N),VL(LDVL,N) DOUBLE PRECISION :: & VR(LDVR,N),SCALE(N),RCONDE(N),RCONDV(N) DOUBLE PRECISION :: & WORK(LWORK) CHARACTER(LEN=1) :: & BALANC,JOBVL,JOBVR,SENSE !FOR SUBROUTINE ASYMTX (FOR COMPUTING EIGENVALUES AND EIGENVECTORS): !INTERNAL VARIABLES INTEGER :: & IER !FOR PERTURBATION SCHEME: !INPUT VARIABLES INTEGER :: & KFLAG,PERT_FLAG DOUBLE PRECISION :: & EPS CHARACTER(LEN=3) :: & WVN_FLAG !INTERNAL VARIABLES DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE :: & W1,W2 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE, SAVE :: & U1,U2 DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: & W DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: & BI,RI,TI,U,V LOGICAL, SAVE :: & FIRST_RUN = .TRUE. !START PROGRAM IF(L_IO == 1) THEN PRINT* PRINT*,'ENTERING GET_GT_GR4' END IF !COMPUTE IDENTITY MATRIX CALL MATIDENT(N,E) !COMPUTE SOME REQUIRED INTERMEDIATE RESULTS CALL INTERMEDIATE_RESULTS (N,VIEW_FLAG,VIEW_ANGLE,QUAD,& DG,NUMDEG,DEGREE,NEXP,DELTA_0_M,MU0,LAYER,OMEGA,TAU,X,L_IO,& DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,RECMU,PSOLP,PSOLM,T,R) !CONSTRUCT "B" MATRIX TPR = T + R TMR = T - R B = MATMUL(TMR,TPR) IF (L_IO == 1) THEN PRINT* PRINT*,'THE B MATRIX LOOKS LIKE:' WRITE(*,30) ((B(I,J),J=1,N),I=1,N) END IF !COMPUTE EIGENVALUE AND EIGENVECTOR INFO FOR GLOBAL TRANSMISSION AND !REFLECTION MATRICES. TO DO THIS, ... IF ( (WVN_FLAG == 'NEW') .OR. & (WVN_FLAG == 'OLD' .AND. PERT_FLAG == 0) ) THEN !... USE LAPACK SUBROUTINE IF (L_IO == 1) THEN PRINT* PRINT*,'IN LAPACK BLOCK' END IF !ALLOCATE MEMORY FOR PERTURBATION SCHEME IF NECESSARY ! IF ( (KFLAG == 1) .AND. (PERT_FLAG == 1) .AND. & ! (FIRST_RUN .OR. DIMENSIONS_RESET) ) THEN IF ( (KFLAG == 1) .AND. (PERT_FLAG == 1) .AND. & (FIRST_RUN .OR. DIMENSIONS_RESET) .AND. (LAYER == 1) ) THEN !print*,'allocating memory' IF (FIRST_RUN) THEN FIRST_RUN = .FALSE. ELSE !print*,'at 2' DEALLOCATE ( W1,W2,U1,U2,W,BI,TI,RI,U,V ) END IF ALLOCATE ( W1(N),W2(N),U1(N,N),U2(N,N),& W(N,NUMLAY,0:NUMDEG),BI(N,N,NUMLAY,0:NUMDEG),& TI(N,N,NUMLAY,0:NUMDEG),RI(N,N,NUMLAY,0:NUMDEG),& U(N,N,NUMLAY,0:NUMDEG),V(N,N,NUMLAY,0:NUMDEG) ) ENDIF IF ( (KFLAG == 1) .AND. (PERT_FLAG == 1) ) THEN !SET UP TO COMPUTE RIGHT & LEFT EIGENVECTORS BALANC = 'N' JOBVL = 'V' JOBVR = 'V' SENSE = 'V' !print*,'saving bi,ti,ri' !SAVE UNPERTURBED B MATRIX INFO FOR LATER CALLS BI(:,:,LAYER,DEGREE) = B TI(:,:,LAYER,DEGREE) = T RI(:,:,LAYER,DEGREE) = R ELSE !SET UP TO COMPUTE RIGHT EIGENVECTORS ONLY BALANC = 'N' JOBVL = 'N' JOBVR = 'V' SENSE = 'V' END IF !COMPUTE EIGENVALUES AND EIGENVECTORS OF MATRIX B !print*,'at sgeevx' ! CALL SGEEVX(BALANC, JOBVL, JOBVR, SENSE, N, B, & ! LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, & ! SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, & ! IWORK, INFO) CALL ASYMTX(B,N,N,N,IER,WR,VR) IF (L_IO == 1) THEN PRINT* DO J=1,N ! WRITE(*,*)' EIGENVALUE=', WR(J), WI(J) ! DO I=1,N ! WRITE(6,*)VR(I,J) ! END DO ! WRITE(*,*)' ' ! WRITE(*,*)' ' WRITE(*,*)' EIGENVALUE=', WR(J) END DO END IF !COMPUTE EIGENVALUES OF MATRIX A AND SET UP MATRIX !OF EIGENVECTORS IN MATRIX D DO J=1,N LAMBP(J) = DSQRT(WR(J)) LAMBM(N+1-J) = -1.0D0*LAMBP(J) DO I=1,N D(I,J) = VR(I,J) END DO ! IF (WR(J) < 0.0D0) THEN ! LAMBP(J) = DSQRT( DABS(WR(J)) ) ! ELSE ! LAMBP(J) = DSQRT(WR(J)) ! END IF END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE LAMBP VECTOR LOOKS LIKE:' WRITE(*,45) (LAMBP(I),I=1,N) END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE D MATRIX LOOKS LIKE:' WRITE(*,30) ((D(I,J),J=1,N),I=1,N) END IF !SAVE EIGEN INFO FOR UNPERTURBED B MATRIX FOR LATER CALLS IF ( (KFLAG == 1) .AND. (PERT_FLAG == 1) ) THEN W(:,LAYER,DEGREE) = WR !EIGENVALUES U(:,:,LAYER,DEGREE) = VR !RIGHT EIGENVECTORS V(:,:,LAYER,DEGREE) = VL !LEFT EIGENVECTORS !print* !print*, 'out of SGEEVX:' !print* !print*, LAYER !print* !print*, DEGREE !print* !print*,'w,u,v:' !print*, W(:,LAYER,DEGREE) !print* !print*, U(:,:,LAYER,DEGREE) !print* !print*, V(:,:,LAYER,DEGREE) END IF ELSE IF(WVN_FLAG == 'OLD' .AND. PERT_FLAG == 1) THEN !... USE PERTURBATION SCHEME SUBROUTINE IF (L_IO == 1) THEN PRINT* PRINT*,'IN PERTURBATION BLOCK' END IF !print* !print*, 'into SPECTRAL_PERTURBATION:' !print* !print*, LAYER !print* !print*, DEGREE !print* !print*, N !print* !print*, BI(:,:,LAYER,DEGREE) !print* !print*, TI(:,:,LAYER,DEGREE) !print* !print*, RI(:,:,LAYER,DEGREE) !print* !print*, B !print* !print*, T !print* !print*, R !print* !print*,'w,u,v:' !print*, W(:,LAYER,DEGREE) !print* !print*, U(:,:,LAYER,DEGREE) !print* !print*, V(:,:,LAYER,DEGREE) !print* !print*, EPS CALL SPECTRAL_PERTURBATION(N,& BI(1,1,LAYER,DEGREE),TI(1,1,LAYER,DEGREE),RI(1,1,LAYER,DEGREE),& B,T,R,& W(1,LAYER,DEGREE),U(1,1,LAYER,DEGREE),V(1,1,LAYER,DEGREE),& EPS,W1,U1,W2,U2) !print* !print*,'w = ',w(:,LAYER,DEGREE) !print* !print*,'w1 = ',w1 !print* !print*,'w2 = ',w2 !print* !print*,'eps = ',eps DO I=1,N ! !FIRST ORDER CORRECTION TO EIGENVALUES ! LAMBP(I) = DSQRT(W(I,LAYER,DEGREE) + EPS*W1(I)) !SECOND ORDER CORRECTION TO EIGENVALUES LAMBP(I) = DSQRT(W(I,LAYER,DEGREE) + EPS*(W1(I) + EPS*W2(I))) ENDDO IF (L_IO == 1) THEN PRINT* PRINT*,'THE LAMBP VECTOR LOOKS LIKE:' WRITE(*,45) (LAMBP(I),I=1,N) END IF ! !FIRST ORDER CORRECTION TO EIGENVECTORS ! D = U(:,:,LAYER,DEGREE) + EPS*U1 !SECOND ORDER CORRECTION TO EIGENVECTORS D = U(:,:,LAYER,DEGREE) + EPS*(U1 + EPS*U2) IF (L_IO == 1) THEN PRINT* PRINT*,'THE D MATRIX LOOKS LIKE:' WRITE(*,30) ((D(I,J),J=1,N),I=1,N) END IF END IF !COMPUTE EXPONENTIAL MATRIX DO J=1,N TEMP4(J) = DEXP(-1.0D0*LAMBP(J)*TAU) END DO CALL MATDIAG(N,TEMP4,EXPO) !COMPUTE U+ AND U- MATRICES UP = 0.0D0 UM = UP DO K=1,N !CONSTRUCT MATRICES 0.5*(E+(T+R)/LAMBP) AND 0.5*(E-(T+R)/LAMBP) !FOR EACH LAMBP DO J=1,N DO I=1,N XP(I,J) = 0.5D0*(E(I,J) + TPR(I,J)/LAMBP(K)) XM(I,J) = 0.5D0*(E(I,J) - TPR(I,J)/LAMBP(K)) END DO END DO !MULTIPLY THE RESULTING MATRICES BY EIGENVECTOR D(K) TO OBTAIN !CORRESPONDING COLUMN OF MATRICES U+ AND U- DO I=1,N DO J=1,N UP(I,K) = UP(I,K) + XP(I,J)*D(J,K) UM(I,K) = UM(I,K) + XM(I,J)*D(J,K) END DO END DO END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE UP MATRIX LOOKS LIKE:' WRITE(*,30) ((UP(I,J),J=1,N),I=1,N) PRINT* PRINT*,'THE UM MATRIX LOOKS LIKE:' WRITE(*,30) ((UM(I,J),J=1,N),I=1,N) END IF !COMPUTE VARIOUS MATRICES NEEDED FOR GLOBAL TRANSMISSION AND !REFLECTION MATRICES CALL MM_IG1G2(N,N,UM,UP,UMIUP,INV_IO) MAT1 = MATMUL(UMIUP,EXPO) MAT2 = TRANSPOSE(MATMUL(UM,MATMUL(EXPO,MAT1)) - UP) MAT3 = TRANSPOSE(MATMUL(UM - MATMUL(UP,UMIUP),EXPO)) MAT4 = TRANSPOSE(MATMUL(UM,E - MATMUL(MAT1,MAT1))) !COMPUTE GLOBAL TRANSMISSION AND REFLECTION MATRICES CALL MM_IG1G2(N,N,MAT4,MAT3,GT,INV_IO) GT = TRANSPOSE(GT) CALL MM_IG1G2(N,N,MAT4,MAT2,GR,INV_IO) GR = TRANSPOSE(GR) IF (L_IO == 1) THEN PRINT* PRINT*,'THE GT MATRIX LOOKS LIKE:' WRITE(*,30) ((GT(I,J),J=1,N),I=1,N) END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE GR MATRIX LOOKS LIKE:' WRITE(*,30) ((GR(I,J),J=1,N),I=1,N) END IF !COMPUTE ((U+)*(U_)INVERSE) TO CHECK GLOBAL R AT ASYMPTOTIC LIMIT AND !DISPLAY IF (L_IO == 1) THEN UPUMI = -1.0D0*MATMUL(UP,UMI) PRINT* PRINT*,'R TEST: THE UPUMI MATRIX LOOKS LIKE:' WRITE(*,30) ((UPUMI(I,J),J=1,N),I=1,N) END IF 30 FORMAT(8(1X,F16.8)) 45 FORMAT( 1X,E12.5) !END PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'LEAVING GET_GT_GR4' END IF END SUBROUTINE GET_GT_GR4 !********************************************************************** !********************************************************************** SUBROUTINE INTERMEDIATE_RESULTS (N,VIEW_FLAG,VIEW_ANGLE,QUAD,& DG,NUMDEG,DEGREE,NEXP,DELTA_0_M,MU0,LAYER,OMEGA,TAU,X,L_IO,& DIMENSIONS_RESET,SZA_RESET,QUAD_RESET,RECMU,PSOLP,PSOLM,T,R) !INPUT : N,VIEW_FLAG,VIEW_ANGLE,QUAD,DG,NUMDEG,DEGREE,NEXP,DELTA_0_M, ! MU0,OMEGA,TAU,X,L_IO !OUTPUT: RECMU,PSOLP,PSOLM,T,R !THIS PROGRAM COMPUTES PHASE FUNCTION AND LOCAL TRANSMISSION AND !REFLECTION QUANTITIES REQUIRED TO COMPUTE GLOBAL TRANSMISSION, !REFLECTION, AND SOURCES FOR A LAYER FOR A NEW SCENE. !PROGRAMMER: MATT CHRISTI !DATE: 4/16/04 !DATA DICTIONARY******************************************************** ! ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE ! USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! L_IO = SUBROUTINE LOCAL I/O FLAG (0=OFF, 1=ON) ! MINV = DIAGONAL MATRIX CONSISTING OF RECIPROCALS OF MU VALUES ! MPMC = MATRIX PRODUCT OF MINV*PM*C ! MPPC = MATRIX PRODUCT OF MINV*PP*C ! MU = VECTOR INITIALLY HOLDING QUADRATURE ROOTS (COSINES OF ! QUADRATURE ANGLES) AND LATER THEIR RECIPROCALS ! MU0 = COSINE OF SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR EQUIVALENTLY DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF PHASE ! FUNCTION ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! PM = PHASE MATRICES CONSISTING OF COMPONENTS FROM -MU DIRECTIONS ! PP = PHASE MATRICES CONSISTING OF COMPONENTS FROM +MU DIRECTIONS ! PSOLM = SOLAR SOURCE VECTOR CONSISTING OF COMPONENTS FROM +MU ! DIRECTIONS ! PSOLP = SOLAR SOURCE VECTOR CONSISTING OF COMPONENTS FROM -MU ! DIRECTIONS ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! COMPUTING THE VALUES OF MU AND W (0=GAUSS,1=LOBATTO) ! R = LOCAL REFLECTION MATRIX ! RECMU = VECTOR HOLDING RECIPROCALS OF THE QUADRATURE ROOTS ! T = LOCAL TRANSMISSION MATRIX ! TAU = TOTAL OPTICAL DEPTH OF THE LAYER ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! W = VECTOR HOLDING QUADRATURE WEIGHTS ! OMEGA = SINGLE SCATTER ALBEDO ! X = VECTOR HOLDING COEFFICIENTS OF COMPOSITE PHASE FUNCTION ! YPLEGM = MATRIX OF LEGENDRE POLYNOMIALS COMPUTED WITH -MU VALUES ! YPLEGP = MATRIX OF LEGENDRE POLYNOMIALS COMPUTED WITH +MU VALUES ! !*********************************************************************** !INTRINSIC SUBPROGRAMS USED BY INTERMEDIATE_RESULTS********************* ! NONE !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY INTERMEDIATE_RESULTS********************** ! GETQUAD2,MATDIAG,MM_D1GD2,PLEG !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,NEXP,DEGREE,NUMDEG,L_IO,LAYER,VIEW_FLAG,DELTA_0_M DOUBLE PRECISION :: & MU0,OMEGA,TAU,VIEW_ANGLE DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & X LOGICAL :: & DIMENSIONS_RESET,SZA_RESET,QUAD_RESET !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & RECMU,PSOLP,PSOLM DOUBLE PRECISION, DIMENSION(N,N) :: & T,R !INTERNAL VARIABLES INTEGER :: & I,J,L,DEG DOUBLE PRECISION :: & SUM,CONSTANT1,LEG_PROD DOUBLE PRECISION, DIMENSION(N,N) :: & MPMC,MPPC,MINV,PP,PM ! DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE :: & ! MU,W ! DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: & ! YPLEGP,YPLEGM LOGICAL, SAVE :: & FIRST_RUN = .TRUE. DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 !START PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'ENTERING INTERMEDIATE_RESULTS' END IF !print*,'layer=',layer !print*,allocated(mu),allocated(w),allocated(yplegp),allocated(yplegm) !START OF SET UP IF BLOCK ! IF ( FIRST_RUN .OR. DIMENSIONS_RESET .OR. & ! SZA_RESET .OR. QUAD_RESET ) THEN IF ( (FIRST_RUN .OR. DIMENSIONS_RESET .OR. & SZA_RESET .OR. QUAD_RESET) .AND. (LAYER == 1) ) THEN IF (FIRST_RUN) THEN FIRST_RUN = .FALSE. ALLOCATE ( MU(N),W(N),YPLEGP(0:NEXP-1,N+2,0:NUMDEG),& YPLEGM(0:NEXP-1,N+2,0:NUMDEG) ) !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN OR LOBATTO QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) ELSE IF (DIMENSIONS_RESET .OR. QUAD_RESET) THEN !print*,'at 3' !print*,'n=',n !print*,'nexp=',nexp !print*,'numdeg=',numdeg !print*,'should be allocated' !print*,allocated(mu),allocated(w),allocated(yplegp),allocated(yplegm) DEALLOCATE ( MU,W,YPLEGP,YPLEGM ) !print*,'should NOT be allocated' !print*,allocated(mu),allocated(w),allocated(yplegp),allocated(yplegm) ALLOCATE ( MU(N),W(N),YPLEGP(0:NEXP-1,N+2,0:NUMDEG),& YPLEGM(0:NEXP-1,N+2,0:NUMDEG) ) !print*,'should be allocated' !print*,allocated(mu),allocated(w),allocated(yplegp),allocated(yplegm) !print*,'at 3a' !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN OR LOBATTO QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE MUs ARE:' WRITE(*,20) (MU(I),I=1,N) PRINT* PRINT*,'THE WEIGHTS ARE:' WRITE(*,20) (W(I),I=1,N) END IF !CREATE LEGENDRE POLYNOMIALS FOR ALL DEGREES TO BE USED !IN CURRENT RUN DO DEG=0,NUMDEG !print*,'deg=',deg !print*,'at 3b' !CREATE LEGENDRE POLYNOMIALS FOR +MUs FOR DEGREE "DEG" CALL PLEG2(MU0,MU,DEG,NEXP,N,YPLEGP(0,1,DEG)) !print*,'at 3c' IF (L_IO == 1) THEN PRINT* PRINT*,'THE YPLEGP(L,MU,DEG) ARE FOR DEGREE = :', DEG WRITE(*,40) ((YPLEGP(L,J,DEG),J=1,N+2),L=0,NEXP-1) 40 FORMAT(7(1X,F17.13)) END IF !CREATE LEGENDRE POLYNOMIALS FOR -MUs FOR DEGREE "DEG" !(TAKING ADVANTAGE OF A LEGENDRE POLYNOMIAL RELATIONSHIP) DO J=1,N+1 DO L=0,NEXP-1 YPLEGM(L,J,DEG) = ((-1.0D0)**(L+DEG))*YPLEGP(L,J,DEG) END DO END DO END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE YPLEGP(L,MU,DEGREE) ARE FOR DEGREE = :', DEGREE WRITE(*,40) ((YPLEGP(L,J,DEGREE),J=1,N+2),L=0,NEXP-1) END IF !END OF SET UP IF BLOCK ENDIF !DO J=1,10 !PRINT* !END DO !PRINT*,'OUTSIDE PLEG2' !PRINT*,'THE YPLEGP(L,MU,DEGREE) ARE FOR DEGREE = :', DEGREE !WRITE(*,40) ((YPLEGP(L,J,DEGREE),J=1,N+2),L=0,NEXP-1) !print* !print*,'layer=',layer !print*,'should be allocated' !print*,allocated(mu),allocated(w),allocated(yplegp),allocated(yplegm) !if ( (.not. allocated(mu)) .or. (.not. allocated(w)) .or. & ! (.not. allocated(yplegp)) .or. (.not. allocated(yplegm)) ) then ! stop !end if !CREATE RECIPROCAL OF EACH MU AND 'M INVERSE' MATRIX DO I=1,N RECMU(I) = 1.0D0 / MU(I) END DO IF (L_IO == 1) THEN PRINT* PRINT*,'THE RECIPROCALS OF THE MUs ARE:' WRITE(*,20) (RECMU(I),I=1,N) END IF CALL MATDIAG(N,RECMU,MINV) !CONSTRUCT COMPOSITE PHASE FUNCTION ... !... FOR HOMOGENEOUS PORTION (PP AND PM) PP = 0.0D0 PM = 0.0D0 IF (OMEGA*TAU > 1.0D-25) THEN DO I=1,N DO J=1,N DO L=0,NEXP-1 LEG_PROD = YPLEGP(L,I,DEGREE)*YPLEGP(L,J,DEGREE) PP(I,J) = PP(I,J) + X(L)*LEG_PROD LEG_PROD = YPLEGM(L,I,DEGREE)*YPLEGP(L,J,DEGREE) PM(I,J) = PM(I,J) + X(L)*LEG_PROD END DO END DO END DO PP = (2.0D0-DELTA_0_M)*PP PM = (2.0D0-DELTA_0_M)*PM END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE P+ MATRIX LOOKS LIKE:' WRITE(*,30) ((PP(I,J),J=1,N),I=1,N) PRINT* PRINT*,'THE P- MATRIX LOOKS LIKE:' WRITE(*,30) ((PM(I,J),J=1,N),I=1,N) END IF !... FOR PARTICULAR PORTION (PSOLP AND PSOLM) PSOLP = 0.0D0 PSOLM = 0.0D0 IF (OMEGA*TAU > 1.0D-25) THEN DO I=1,N DO L=0,NEXP-1 LEG_PROD = YPLEGM(L,I,DEGREE)*YPLEGM(L,N+1,DEGREE) PSOLP(I) = PSOLP(I) + X(L)*LEG_PROD LEG_PROD = YPLEGP(L,I,DEGREE)*YPLEGM(L,N+1,DEGREE) PSOLM(I) = PSOLM(I) + X(L)*LEG_PROD END DO END DO PSOLP = (2.0D0-DELTA_0_M)*PSOLP PSOLM = (2.0D0-DELTA_0_M)*PSOLM END IF IF (L_IO == 1) THEN PRINT* PRINT*,'THE P+ SOLAR VECTOR LOOKS LIKE:' WRITE(*,20) (PSOLP(I),I=1,N) PRINT* PRINT*,'THE P- SOLAR VECTOR LOOKS LIKE:' WRITE(*,20) (PSOLM(I),I=1,N) END IF !PHASE FUNCTION NORM TEST FOR M=0 (SHOULD BE 1 FOR EACH QUADRATURE ANGLE) IF ((DEGREE == 0).AND.(L_IO == 1)) THEN PRINT* DO I=1,N SUM = 0.0D0 DO J=1,N SUM = SUM + 0.5D0*W(J)*(PP(I,J)+PM(I,J)) END DO WRITE(*,130) I,SUM 130 FORMAT(1X,'THE SUM OF NORM TEST ',I3,' IS: ',F18.12) END DO END IF !SOLAR PHASE FUNCTION NORM TEST FOR M=0 (SHOULD BE 1) IF ((DEGREE == 0).AND.(L_IO == 1)) THEN PRINT* SUM = 0.0D0 DO J=1,N SUM = SUM + 0.5D0*W(J)*(PSOLP(J)+PSOLM(J)) END DO WRITE(*,140) SUM 140 FORMAT(1X,'THE SUM OF SOLAR NORM TEST IS: ',F18.12) END IF !CONSTRUCT LOCAL REFLECTANCE AND TRANSMITTANCE MATRICES (r AND t) CALL MM_D1GD2(N,RECMU,PM,W,MPMC) CALL MM_D1GD2(N,RECMU,PP,W,MPPC) CONSTANT1 = (1.0D0+DELTA_0_M)*0.25D0*OMEGA DO I=1,N DO J=1,N R(I,J) = -1.0D0*CONSTANT1*MPMC(I,J) T(I,J) = -1.0D0*MINV(I,J) + CONSTANT1*MPPC(I,J) END DO END DO IF(L_IO == 1) THEN PRINT* PRINT*,'THE r MATRIX LOOKS LIKE:' WRITE(*,150) ((R(I,J),J=1,N),I=1,N) PRINT* PRINT*,'THE t MATRIX LOOKS LIKE:' WRITE(*,150) ((T(I,J),J=1,N),I=1,N) 150 FORMAT(8(1X,F14.8)) END IF 20 FORMAT(1X,F16.8) 30 FORMAT(8(1X,F16.8)) !END PROGRAM IF (L_IO == 1) THEN PRINT* PRINT*,'LEAVING INTERMEDIATE_RESULTS' END IF END SUBROUTINE INTERMEDIATE_RESULTS !********************************************************************** !********************************************************************** SUBROUTINE PLANCK(WVN,NUMLVL,TEMP,B_T) !INPUT : WVN,NUMLAY,TEMP !OUTPUT: B_T !THIS PROGRAM COMPUTES THERMALLY-EMITTED SPECTRAL RADIANCES USING THE !PLANCK FUNCTION OF EMISSION !PROGRAMMER: MATT CHRISTI !DATA DICTIONARY**************************************************************** ! ! B_T = VECTOR OF SPECTRAL RADIANCES USING THE PLANCK FUNCTION ! OF EMISSION (in mW/(m^2 cm^-1 ster)) ! C = SPEED OF LIGHT (in m/s) ! H = PLANCK'S CONSTANT (in J*s) ! K = BOLTZMANN'S CONSTANT (in J/K) ! NUMLAY = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE ! TEMP = VECTOR (I.E. VERTICAL PROFILE) OF ATMOSPHERIC TEMPERATURES ! (in K) ! WVN = WAVENUMBER AT WHICH PLANCK FUNCTION IS COMPUTED (in cm^-1) ! !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER NUMLVL DOUBLE PRECISION WVN DOUBLE PRECISION, DIMENSION(NUMLVL) :: & TEMP !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(NUMLVL) :: & B_T !INTERNAL VARIABLES INTEGER I DOUBLE PRECISION C,H,K,WVNUM C = 2.99792458D8 H = 6.62606876D-34 K = 1.3806503D-23 !CONVERT WAVENUMBER UNITS FROM cm^-1 TO m^-1 WVNUM = 100.0D0*WVN !COMPUTING PLANCK FUNCTION DO I=1,NUMLVL B_T(I) = (2.0D0*H*C**2*WVNUM**3)/ & (DEXP((H*C*WVNUM)/(K*TEMP(I))) - 1.0D0) END DO !CONVERT SPECTRAL RADIANCE UNITS FROM W/(m^2 m^-1 ster) !TO mW/(m^2 cm^-1 ster) B_T = 1.0D5*B_T END SUBROUTINE PLANCK !******************************************************************************* !******************************************************************************* SUBROUTINE SPECTRAL_PERTURBATION(N,BI,TI,RI,BF,TF,RF,W0,U0,V0,EPS,& W1,U1,W2,U2) !INPUT : N,BI,TI,RI,BF,TF,RF,W0,U0,V0,EPS !OUTPUT: W1,U1,W2,U2 !THIS PROGRAM GENERATES THE CORRECTIONS APPLIED TO THE EIGEN SPECTRUM OF !AN INITIAL MATRIX "BI" OF THE FORM BI = (TI-RI)(TI+RI) IN ORDER TO !OBTAIN THE EIGEN SPECTRUM OF THE PERTURBED MATRIX "BF" !PROGRAMMERS: PHIL GABRIEL & MATT CHRISTI !DATE: 5/4/04 !DATA DICTIONARY******************************************************** ! ! BI = INITIAL STATE OF MATRIX WHICH IS SUBJECTED TO PERTURBATION ! BF = FINAL STATE OF MATRIX AFTER SUBJECTED TO PERTURBATION ! EPS = MAGNITUDE BY WHICH OPTICAL PROPERTY ASSOCIATED WITH BI IS ! PERTURBED TO GENERATE MATRIX BF (1.0D-10 <= EPS <= 1.0D-6 ! IS RECOMMENDED) ! N = DIMENSION OF MATRICES IN THE SUBROUTINE ! RF = FINAL STATE OF LOCAL REFLECTION MATRIX ! RI = INITIAL STATE OF LOCAL REFLECTION MATRIX ! TF = FINAL STATE OF LOCAL TRANSMISSION MATRIX ! TI = INITIAL STATE OF LOCAL TRANSMISSION MATRIX ! U0 = RIGHT EIGENVECTORS ASSOCIATED WITH MATRIX BI ! U1 = 1ST ORDER CORRECTION APPLIED TO RIGHT EIGENVECTORS OF ! MATRIX BI TO OBTAIN RIGHT EIGENVECTORS OF MATRIX BF ! U2 = 2ND ORDER CORRECTION APPLIED TO RIGHT EIGENVECTORS OF ! MATRIX BI TO OBTAIN RIGHT EIGENVECTORS OF MATRIX BF ! V0 = LEFT EIGENVECTORS ASSOCIATED WITH MATRIX BI ! W0 = EIGENVALUES ASSOCIATED WITH MATRIX BI ! W1 = 1ST ORDER CORRECTION APPLIED TO EIGENVALUES OF MATRIX BI ! TO OBTAIN EIGENVALUES OF MATRIX BF ! W2 = 2ND ORDER CORRECTION APPLIED TO EIGENVALUES OF MATRIX BI ! TO OBTAIN EIGENVALUES OF MATRIX BF ! !*********************************************************************** !INTRINSIC SUBPROGRAMS USED BY SPECTRAL_PERTURBATION******************** ! DOT_PRODUCT,MATMUL,TRANSPOSE !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY SPECTRAL_PERTURBATION********************* ! NONE !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N DOUBLE PRECISION :: & EPS DOUBLE PRECISION, DIMENSION(N) :: & W0 DOUBLE PRECISION, DIMENSION(N,N) :: & BI,TI,RI,BF,TF,RF,U0,V0 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: & W1,W2 DOUBLE PRECISION, DIMENSION(N,N) :: & U1,U2 !INTERNAL VARIABLES INTEGER :: & I,J,K DOUBLE PRECISION, DIMENSION(N) :: & DP DOUBLE PRECISION, DIMENSION(N,N) :: & A1U1,A2U0,B1,B2,C,PD,R1,TERM1,TERM2,T1,T1MR1,T1PR1,V0DPI,& WKJ,W1U0 LOGICAL :: & SECOND_ORDER = .TRUE. !START PROGRAM ! PRINT* ! PRINT*,'ENTERING SPECTRAL_PERTURBATION' !DO PRELIMINARY COMPUTATIONS T1 = (TF - TI)/EPS R1 = (RF - RI)/EPS ! T1MR1 = T1 - R1 ! T1PR1 = T1 + R1 B2 = MATMUL(T1 - R1,T1 + R1) ! B1 = (BF - BI - EPS*EPS*B2)/EPS PD = MATMUL( (BF-BI-EPS*EPS*B2)/EPS , U0) !COMPUTE 1ST ORDER CORRECTION TO EIGENVALUES & EIGENVECTORS DO J=1,N ! DP(J) = DOT_PRODUCT(U0(:,J),V0(:,J)) ! V0DPI(:,J) = V0(:,J)/DP(J) V0DPI(:,J) = V0(:,J)/DOT_PRODUCT(U0(:,J),V0(:,J)) W1(J) = DOT_PRODUCT(PD(:,J),V0DPI(:,J)) W1U0(:,J) = W1(J)*U0(:,J) WKJ(J,J) = 0.0D0 END DO DO K=1,N DO J=K+1,N WKJ(K,J) = 1.0D0/(W0(K) - W0(J)) WKJ(J,K) = -WKJ(K,J) END DO END DO ! C = TRANSPOSE(WKJ)*MATMUL(TRANSPOSE(W1U0 - PD),V0DPI) ! U1 = MATMUL(U0,TRANSPOSE(C)) C = WKJ*TRANSPOSE( MATMUL(TRANSPOSE(W1U0 - PD),V0DPI) ) U1 = MATMUL(U0,C) !COMPUTE 2ND ORDER CORRECTION TO EIGENVALUES & EIGENVECTORS IF (SECOND_ORDER) THEN A1U1 = MATMUL(PD,C) A2U0 = MATMUL(B2,U0) DO J=1,N TERM2(:,J) = W2(J)*U0(:,J) TERM1(:,J) = A1U1(:,J) - W1(J)*U1(:,J) + A2U0(:,J) W2(J) = DOT_PRODUCT(TERM1(:,J),V0DPI(:,J)) END DO C = WKJ*TRANSPOSE( MATMUL(TRANSPOSE(TERM2 - TERM1),V0DPI) ) U2 = MATMUL(U0,C) END IF ! PRINT* ! PRINT*,'LEAVING SPECTRAL_PERTURBATION' END SUBROUTINE SPECTRAL_PERTURBATION !*************************************************************************** !*************************************************************************** SUBROUTINE SURF_REF1_3(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,SI,RG) !INPUT : QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE !OUTPUT: SI,RG !SURFACE TYPE: ISOTROPIC (LAMBERTIAN) !THIS PROGRAM RETURNS AN NxN MATRIX INDICATIVE OF THE REFLECTION !CHARACTERISTICS OF THE SURFACE FOR A MODEL ATMOSPHERE/SURFACE SYSTEM !PROGRAMMER: MATT CHRISTI !DATA DICTIONARY**************************************************************** ! ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL ! BE USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! COMPUTING ENTRIES OF RG (0=GAUSS,1=LOBATTO) ! R = ARRAY HOLDING QUADRATURE ROOTS ! RG = SURFACE REFLECTION MATRIX ("R GROUND") ! SI = NUMERICAL INTEGRATION OF MUs USING QUADRATURE ! ("SUM, ISOTROPIC") ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE ! (0=OFF, 1=ON) ! W = ARRAY HOLDING QUADRATURE WEIGHTS ! !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY SURF_REF1_3*************************************** ! GETQUAD2 !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,VIEW_FLAG DOUBLE PRECISION :: & VIEW_ANGLE !OUTPUT VARIABLES DOUBLE PRECISION :: & SI,RG(N,N) !INTERNAL VARIABLES INTEGER :: & I,J DOUBLE PRECISION :: & SUM DOUBLE PRECISION, DIMENSION(N) :: & MU,W INTEGER, SAVE :: & SURF_IO = 0 DOUBLE PRECISION, PARAMETER :: & PI=3.1415926535897932D0 !START PROGRAM IF (SURF_IO == 1) THEN PRINT* PRINT*,'ENTERING SURF_REF1_3' END IF !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN QUADRATURE OR LOBATTO !QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) IF (SURF_IO == 1) THEN PRINT* PRINT*,'THE MU ARE:' WRITE(*,*) (MU(I),I=1,N) PRINT* PRINT*,'THE WEIGHTS ARE:' WRITE(*,*) (W(I),I=1,N) END IF !COMPUTE SURFACE REFLECTION MATRIX RG DO I=1,N DO J=1,N RG(I,J) = MU(J)*W(J) END DO END DO !GAUSSIAN WEIGHT TEST AND COMPUTATION OF SI SUM = 0.0D0 SI = 0.0D0 DO J=1,N SUM = SUM + W(J) SI = SI + W(J)*MU(J) END DO IF (SURF_IO == 1) THEN PRINT* WRITE(*,35) SUM 35 FORMAT(1X,'THE SUM OF THE WEIGHTS = ',F14.12) PRINT* WRITE(*,40) SI 40 FORMAT(1X,'SI = ',F14.12) END IF !DISPLAY SURFACE REFLECTION MATRIX RG IF (SURF_IO == 1) THEN PRINT* PRINT*,'THE RG MATRIX LOOKS LIKE:' WRITE(*,50) ((RG(I,J),J=1,N),I=1,N) 50 FORMAT(8(1X,F9.5)) END IF !END PROGRAM IF (SURF_IO == 1) THEN PRINT* PRINT*,'LEAVING SURF_REF1_3' END IF END SUBROUTINE SURF_REF1_3 !******************************************************************************* !******************************************************************************* SUBROUTINE SURF_REF2_2(MU0,QUAD,DG,N,NEXP,DEGREE,NUMDEG,VIEW_FLAG,& VIEW_ANGLE,WARNING_FLAG,CHI_SURF,DIMENSIONS_RESET,SZA_RESET,& QUAD_RESET,SI,RG,GAMMA) !INPUT : MU0,QUAD,DG,N,NEXP,DEGREE,NUMDEG,VIEW_FLAG,VIEW_ANGLE,WARNING_FLAG, ! CHI_SURF !OUTPUT: SI,RG,GAMMA !SURFACE TYPE: SURFACE READ IN FROM FILE !THIS PROGRAM RETURNS AN NxN MATRIX INDICATIVE OF THE REFLECTION !CHARACTERISTICS OF THE SURFACE FOR A MODEL ATMOSPHERE/SURFACE SYSTEM !PROGRAMMER: MATT CHRISTI !DATA DICTIONARY**************************************************************** ! ! CHI_SURF = VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL ! EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED) ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL ! BE USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! DEGREE = FOURIER COMPONENT TO BE COMPUTED ! GAMMA = VECTOR OF BIDIRECTIONAL REFECTIVITIES FOR SOLAR BEAM ! MU0 = COSINE OF SOLAR ZENITH ANGLE ! N = NUMBER OF UPWARD (OR DOWNWARD) STREAMS ! NEXP = NUMBER OF TERMS IN LEGENDRE POLYNOMIAL EXPANSION OF SURFACE ! BIDIRECTIONAL REFLECTIVITY ! NUMDEG = TOTAL NUMBER OF FOURIER COMPONENTS TO BE COMPUTED - 1 ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! COMPUTING ENTRIES OF RG (0=GAUSS,1=LOBATTO) ! R = ARRAY HOLDING QUADRATURE ROOTS ! RG = SURFACE REFLECTION MATRIX ("R GROUND") ! SI = NUMERICAL INTEGRATION OF MUs USING QUADRATURE ! ("SUM, ISOTROPIC") ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE (0=OFF, 1=ON) ! W = ARRAY HOLDING QUADRATURE WEIGHTS ! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG (0=OFF, 1=ON) ! !******************************************************************************* !EXTERNAL SUBPROGRAMS USED BY SURF_REF2_2*************************************** ! GETQUAD2,PLEG !******************************************************************************* IMPLICIT NONE !INPUT VARIABLES INTEGER :: & QUAD,DG,N,NEXP,DEGREE,NUMDEG,VIEW_FLAG,WARNING_FLAG DOUBLE PRECISION :: & MU0,VIEW_ANGLE DOUBLE PRECISION, DIMENSION(0:NEXP-1) :: & CHI_SURF LOGICAL :: & DIMENSIONS_RESET,SZA_RESET,QUAD_RESET !OUTPUT VARIABLES DOUBLE PRECISION :: & SI DOUBLE PRECISION, DIMENSION(N) :: & GAMMA DOUBLE PRECISION, DIMENSION(N,N) :: & RG !INTERNAL VARIABLES INTEGER :: & I,J,L,DEG,NUM,NUMEXP DOUBLE PRECISION :: & SUM,LEG_PROD,DELTA_0_M ! DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE :: & ! MU,W ! DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: & ! YPLEGP,YPLEGM INTEGER, SAVE :: & SURF_IO = 0 LOGICAL, SAVE :: & FIRST_RUN = .TRUE. DOUBLE PRECISION, PARAMETER :: & PI = 3.1415926535897932D0 !START PROGRAM IF (SURF_IO == 1) THEN PRINT* PRINT*,'ENTERING SURF_REF2_2' END IF !THE FIRST TIME SUBROUTINE IS CALLED, SET UP AND SAVE VECTORS AND !MATRICES THAT WILL BE USED REPEATEDLY IF (SURF_IO == 1) THEN IF (FIRST_RUN .OR. DIMENSIONS_RESET .OR. & SZA_RESET .OR. QUAD_RESET) THEN IF (FIRST_RUN) THEN FIRST_RUN = .FALSE. ALLOCATE ( MU(N),W(N),YPLEGP(0:NEXP-1,N+2,0:NUMDEG),& YPLEGM(0:NEXP-1,N+2,0:NUMDEG) ) !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN OR LOBATTO QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) ELSE IF (DIMENSIONS_RESET .OR. QUAD_RESET) THEN !print*,'at 4' DEALLOCATE ( MU,W,YPLEGP,YPLEGM ) ALLOCATE ( MU(N),W(N),YPLEGP(0:NEXP-1,N+2,0:NUMDEG),& YPLEGM(0:NEXP-1,N+2,0:NUMDEG) ) !OBTAIN MUs AND WEIGHTS FOR GAUSSIAN OR LOBATTO QUADRATURE CALL GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) END IF IF(SURF_IO == 1) THEN PRINT* PRINT*,'THE MUs ARE:' WRITE(*,*) (MU(I),I=1,N) 20 FORMAT(8(1X,F14.8)) PRINT* PRINT*,'THE WEIGHTS ARE:' WRITE(*,*) (W(I),I=1,N) END IF !CREATE LEGENDRE POLYNOMIALS FOR ALL DEGREES TO BE USED !IN CURRENT RUN DO DEG=0,NUMDEG !CREATE LEGENDRE POLYNOMIALS FOR +MUs FOR DEGREE "DEG" CALL PLEG2(MU0,MU,DEG,NEXP,N,YPLEGP(0,1,DEG)) IF (SURF_IO == 1) THEN PRINT* PRINT*,'THE YPLEGP(L,MU,DEG) ARE FOR DEGREE = :', DEG WRITE(*,30) ((YPLEGP(L,J,DEG),J=1,N+2),L=0,NEXP-1) 30 FORMAT(6(1X,F17.13)) END IF !CREATE LEGENDRE POLYNOMIALS FOR -MUs FOR DEGREE "DEG" !(TAKING ADVANTAGE OF A LEGENDRE POLYNOMIAL RELATIONSHIP) DO J=1,N+1 DO L=0,NEXP-1 YPLEGM(L,J,DEG) = ((-1.0D0)**(L+DEG))*YPLEGP(L,J,DEG) END DO END DO !END OF LEGENDRE POLYNOMIAL LOOP END DO !PRINT* !PRINT*,'INSIDE SURF_REF2_2' !PRINT*,'CHI_SURF = ',CHI_SURF !END OF SETUP BLOCK ENDIF ENDIF !PREPARE TO COMPUTE RG AMD GAMMA IF (DEGREE == 0) THEN DELTA_0_M = 1.0D0 ELSE DELTA_0_M = 0.0D0 END IF !CONSTRUCT SURFACE REFLECTION MATRIX (RG) RG = 0.0D0 DO I=1,N DO J=1,N DO L=0,NEXP-1 LEG_PROD = YPLEGM(L,I,DEGREE)*YPLEGP(L,J,DEGREE) RG(I,J) = RG(I,J) + CHI_SURF(L)*LEG_PROD IF (SURF_IO == 1 .AND. I==1 .AND. J==1) THEN PRINT*,'L = ', L PRINT*,'CHI_SURF(L) = ', CHI_SURF(L) PRINT*,'YPLEGM(L,I,DEGREE) = ', YPLEGM(L,I,DEGREE) PRINT*,'YPLEGP(L,J,DEGREE) = ', YPLEGP(L,J,DEGREE) PRINT*,'RG(I,J) = ',RG(I,J) END IF END DO RG(I,J) = W(J)*MU(J)*(2.0D0 - DELTA_0_M)*RG(I,J) END DO END DO !DISPLAY RG IF (SURF_IO == 1) THEN PRINT* PRINT*,'THE RG MATRIX LOOKS LIKE:' WRITE(*,31) ((RG(I,J),J=1,N),I=1,N) 31 FORMAT(9(1X,F9.5)) END IF !CONSTRUCT SURFACE REFLECTION VECTOR FOR DIRECT BEAM (GAMMA) GAMMA = 0.0D0 DO I=1,N DO L=0,NEXP-1 LEG_PROD = YPLEGP(L,I,DEGREE)*YPLEGM(L,N+1,DEGREE) GAMMA(I) = GAMMA(I) + CHI_SURF(L)*LEG_PROD IF (SURF_IO == 1 .AND. I==1) THEN PRINT*,'L = ', L PRINT*,'CHI_SURF(L) = ', CHI_SURF(L) PRINT*,'YPLEGP(L,I,DEGREE) = ', YPLEGP(L,I,DEGREE) PRINT*,'YPLEGM(L,N+1,DEGREE) = ',YPLEGM(L,N+1,DEGREE) PRINT*,'GAMMA(I) = ',GAMMA(I) END IF END DO GAMMA(I) = (2.0D0 - DELTA_0_M)*GAMMA(I) END DO !DISPLAY GAMMA IF (SURF_IO == 1) THEN PRINT* PRINT*,'THE GAMMA VECTOR LOOKS LIKE:' WRITE(*,32) (GAMMA(I),I=1,N) 32 FORMAT(1X,F9.5) END IF !GAUSSIAN WEIGHT TEST AND COMPUTATION OF SI SUM = 0.0D0 SI = 0.0D0 DO J=1,N SUM = SUM + W(J) SI = SI + W(J)*MU(J) END DO IF (SURF_IO == 1) THEN PRINT* WRITE(*,35) SUM 35 FORMAT(1X,'THE SUM OF THE WEIGHTS = ',F14.12) PRINT* WRITE(*,40) SI 40 FORMAT(1X,'SI = ',F14.12) END IF !DIFFUSE REFLECTIVITY NORM TEST FOR M=0 (SHOULD BE EQUAL TO THE ALBEDO !FOR EACH QUADRATURE ANGLE) IF((DEGREE == 0).AND.(SURF_IO == 1)) THEN PRINT* DO I=1,N SUM = 0.0D0 DO J=1,N SUM = SUM + 2.0D0*RG(I,J) END DO WRITE(*,45) I,SUM 45 FORMAT(1X,'THE SUM OF DIFFUSE NORM TEST ',I3,' IS: ',F18.12) END DO END IF !DIRECT BEAM REFLECTIVITY NORM TEST FOR M=0 (SHOULD BE EQUAL TO THE !ALBEDO) IF((DEGREE == 0).AND.(SURF_IO == 1)) THEN PRINT* SUM = 0.0D0 DO J=1,N SUM = SUM + 2.0D0*W(J)*MU(J)*GAMMA(J) END DO WRITE(*,50) SUM 50 FORMAT(1X,'THE SUM OF DIRECT BEAM NORM TEST IS: ',F18.12) END IF !END PROGRAM IF (SURF_IO == 1) THEN PRINT* PRINT*,'LEAVING SURF_REF2_2' END IF END SUBROUTINE SURF_REF2_2 !******************************************************************************* !******************************************************************************* END MODULE RADIANT2S