c ******************************************************************* c Test version of the vector DOM program for c homogeneous layer c c ******************************************************************* Program VDOM_HOM Implicit none INTEGER NLEG, NS, M, I, J, K, L, Jc, Jr, K1, NS_US INTEGER KK, KS, N_SCA, NLEG_MAX c Variables for post=processing INTEGER N_T0, N_MU0 DOUBLE PRECISION TAU_US, XMU, AZIM DOUBLE PRECISION E_t0m0, C_tm0m, C_tnjm, T_t0t_njm, RELAZM, D_M, $ SCAT, SZA, D_COS, D_SIN DOUBLE PRECISION B_L(4,4), S1(4), S2(4), STOC_D(4), QPU(4), $ F_INC(4,1), STOC_U(4) CHARACTER ( LEN = 6 ) KEY_OUT CHARACTER ( LEN = 80 ) OUT_NAME_SS, OUT_NAME_MS DOUBLE PRECISION, ALLOCATABLE :: Greek(:,:), $ DWF_P2U(:,:),DWF_M2U(:,:),ALP_U(:,:), CF(:,:), $ DWF_P2_M(:,:),DWF_M2_M(:,:), $ DWF_P2(:,:), DWF_M2(:,:), T0_LAY(:), SCAT_ANG(:), $ B_LU(:,:,:), B_LD(:,:,:), ALP(:,:), $ B_LUM(:,:,:), B_LDM(:,:,:), $ ALP0(:,:), TTT(:,:), PI_LM(:,:,:), GLM(:,:,:), $ XM(:), AG(:), PIE(:,:), SS_SOL_U(:,:,:,:,:), $ MH(:,:), RHS(:,:), AT_COF(:,:), BT_COF(:,:), $ Part_p(:), Part_m(:), aj_as(:,:), bj_as(:,:), $ QP1_AS(:,:), QM1_AS(:,:), $ FIP(:,:), FIM(:,:), HEL(:,:), H4(:), H3(:), $ PSIP(:,:), PSIM(:,:), Norm(:), XM_MU0(:), $ SUM_E(:,:), SUM_F(:,:), SV(:,:), GLM_1(:,:), GLM_2(:,:), $ RR(:), RI(:), VR(:,:), VI(:,:), WORK(:), $ XM_US(:), TT(:), FIN_UP_PP(:,:,:,:), FIN_DOWN_PP(:,:,:,:), $ COF_11(:) DOUBLE PRECISION, ALLOCATABLE :: NJc(:), NJr(:), IPIV(:) INTEGER FAIL, INFO, SIGN_M DOUBLE PRECISION OMEG, MU0, TAU, TAU_L, RE_EXP, ALBEDO, IM_EXP DOUBLE PRECISION Raj, Rbj, Raj_1, Raj_2, Rbj_1, Rbj_2 DOUBLE PRECISION T1, T2, T3, T4 DOUBLE PRECISION M_FO(4,4), Z_matr(4,4) DOUBLE PRECISION SS_SOL_D(4,1) LOGICAL SS_COR COMPLEX (kind(0.0d0)) :: C_t1, C_t2 F_INC = 0.d0 ; F_INC(1,1) = 1.d0 CALL FILEREAD c CALL PH_MATRIX c stop IF( KEY_OUT .EQ. 'REFLEC' ) TAU_US = 0.0d0 CALL GAULEG ( 0.D0, 1.D0, XM, AG, NS ) c CALL GAULEG ( 0.D0, 1.D0, XM, AG, NS-1 ) c write( *,* ) ' ... Gauss angles ...' c XM(NS) = 1.0 ; AG(NS) = 0.0 c AG(:) = 0.0 c DO I = 1, NS c WRITE( *,'(2x,i3,2x,e12.5,2x,f7.2)') c $ I, XM(I), DACOS(XM(I))/3.14159*180. c ENDDO C Factorial matrix CF: c for m = 0..N and l = m..N: ==> CF(l+1,m+1) = (l-m)! / (l+m)! CF = 0.0d0 ; CF(1:NLEG,1) = 1.d0 DO J = 1, NLEG DO I = J+1, NLEG CF(I,J+1) = CF(I,J) / DBLE( (I+J-1)*(I-J) ) ENDDO ENDDO CF(:,:) = DSQRT( CF(:,:) ) FIN_DOWN_PP = 0.d0 ; FIN_UP_PP = 0.d0 DO M = 0, NLEG-1 WRITE( *,* ) ' ************************************ ' WRITE( *,* ) ' ********* M = ',M,' ******' WRITE( *,* ) ' ************************************ ' SIGN_M = (1 - 2*MODULO(M,2)) IF( M .EQ. 0 ) D_M = 0.5d0 ; IF( M > 0 ) D_M = 1.d0 D_COS = D_M*DCOS( DBLE(M)*RELAZM ) D_SIN = D_M*DSIN( DBLE(M)*RELAZM ) CALL WIGNER_FUNCTION( XM, NS, M, NLEG, DWF_P2, DWF_M2 ) c DO J = M+1, NLEG c ALP0(J,1) = PLGNDR( J-1,M,MU0 )*CF(J,M+1)*SIGN_M c DO I = 1, NS c ALP(J,I) = PLGNDR( J-1,M,XM(I) )*CF(J,M+1)*SIGN_M c ENDDO c ENDDO CALL WIGNER_PLG( XM_MU0, N_MU0, M, NLEG, ALP0 ) CALL WIGNER_PLG( XM, NS, M, NLEG, ALP ) c write(*,*) ' After wigner ' c DO J = M+1, NLEG c WRITE(*,'(2x,I3,3(2xe15.8))') J-1, c $ ALP(J,5), TTT(J,5),ALP(J,5) - TTT(J,5) c ENDDO c write(*,*) ' After do ' c stop CALL WIGNER_FUNCTION( XM_US, NS_US, M, NLEG, DWF_P2U, DWF_M2U ) CALL WIGNER_PLG( XM_US, NS_US, M, NLEG, ALP_U ) c DO J = M+1, NLEG c DO I = 1, NS_US c ALP_U(J,I) = PLGNDR( J-1,M,XM_US(I) )*CF(J,M+1)*SIGN_M c ENDDO c ENDDO CALL MATRIX_QU SV = MATMUL( SUM_F,SUM_E ) CALL DGEEV ('V','V', 4*NS,SV,4*NS, RR,RI, VI,4*NS, VR,4*NS, $ WORK, 20*4*NS, FAIL) IF (FAIL.NE.0) WRITE(*,*) 'EV problem --> ifail = ', fail CALL FI_AND_PSI c Particularly solution: aj and bj coefficients for part. sol DO KS = 1,N_MU0 CALL COEF_FOR_PART_SOL $ ( QM1_AS(:,KS),QP1_AS(:,KS), aj_as(:,KS),bj_as(:,KS) ) ENDDO DO I = 1,NS FIP(1:4*NS, 4*I-3:4*I) = FIP(1:4*NS,4*I-3:4*I)/XM(I)*0.5 FIM(1:4*NS, 4*I-3:4*I) = FIM(1:4*NS,4*I-3:4*I)/XM(I)*0.5 ENDDO c ... Matrix GLM for post-processing GLM = 0.d0 DO L = M+1,NLEG DO J = 1,4*NS DO I = 1,NS H3(4*I-3:4*I) = FIP(J,4*I-3:4*I)*AG(I) H4(4*I-3:4*I) = FIM(J,4*I-3:4*I)*AG(I) ENDDO GLM(L,J,:) = MATMUL( TRANSPOSE( PI_LM(L,:,:)),H3 ) S1(:) = MATMUL( TRANSPOSE( PI_LM(L,:,:)),H4 ) S1(3:4) = -S1(3:4) GLM(L,J,:) = GLM(L,J,:) + (1-2*MODULO(L-1-M,2))*S1(:) ENDDO ENDDO DO KK = 1,N_T0 TAU_L = T0_LAY(KK) IF( KEY_OUT .EQ. 'TRANSM' ) TAU_US = TAU_L DO KS = 1,N_MU0 MU0 = XM_MU0(KS) c Particularly solution: Part_p and Part_m for RHS calculation CALL PART_SOLUTION ( 0.0d0, aj_as(:,KS),bj_as(:,KS) ) c ... Righthand site for coefficients matrix (upper boun. cond.) RHS(1:4*NS,KS) = -Part_p(:) CALL PART_SOLUTION ( TAU_L, aj_as(:,KS),bj_as(:,KS) ) FORALL( K = 1:NS ) H3(K) = Part_p(4*K-3)*AG(K)*XM(K) T1 = 2.0*ALBEDO*( MU0*DEXP(-TAU_L/MU0) + SUM(H3(1:NS)) ) c ... See comment in the PART_SOLUTION !!! c FORALL( I = 1:NS ) Part_m(4*I-1:4*I) = -Part_m(4*I-1:4*I) FORALL( K = 1:NS ) Part_m(4*K-3) = Part_m(4*K-3) - T1 c .... Righthand site for coefficients matrix (bottom bound. cond.) RHS(4*NS+1:8*NS,KS) = -Part_m(:) ENDDO c Matrix for coefficients of the homogeneous solution CALL MAIN_MATRIX c LU factorization of matrix MH CALL DGETRF( 8*NS, 8*NS, MH, 8*NS, IPIV, INFO ) IF (INFO.NE.0) WRITE(*,*) 'DGETRF problem --> INFO = ', INFO c Solution of linear equation for single RHS CALL DGETRS( 'n', 8*NS, N_MU0, MH, 8*NS,IPIV,RHS,8*NS,INFO ) IF (INFO.NE.0) WRITE(*,*) 'DGETRS problem --> INFO = ', INFO c ... A and B coefficients for homogeneous solution I = 0 DO J = 1,8*NS - 1,2 ; I = I + 1 DO KS = 1,N_MU0 AT_COF(I,KS) = RHS(J,KS) ; BT_COF(I,KS) = RHS(J+1,KS) ENDDO ENDDO DO K1 = 1,NS_US DO KS = 1,N_MU0 MU0 = XM_MU0(KS) c .... STOC_U as a output for each harmonic CALL PP_UP (XM_US(K1), TAU_US, K1, AT_COF(:,KS), $ BT_COF(:,KS), aj_as(:,KS), bj_as(:,KS) ) FIN_UP_PP(K1,KS,KK,1:2) = FIN_UP_PP(K1,KS,KK,1:2) $ + D_COS*STOC_U(1:2) FIN_UP_PP(K1,KS,KK,3:4) = FIN_UP_PP(K1,KS,KK,3:4) $ + D_SIN*STOC_U(3:4) c .... STOC_D as a output for each harmonic CALL PP_DOWN(XM_US(K1), TAU_US, K1, AT_COF(:,KS), $ BT_COF(:,KS), aj_as(:,KS), bj_as(:,KS) ) FIN_DOWN_PP(K1,KS,KK,1:2) = FIN_DOWN_PP(K1,KS,KK,1:2) $ + D_COS*STOC_D(1:2) FIN_DOWN_PP(K1,KS,KK,3:4) = FIN_DOWN_PP(K1,KS,KK,3:4) $ + D_SIN*STOC_D(3:4) ENDDO ! End for solar zenith angles loop ENDDO ENDDO ! End for tau_layer loop ENDDO ! End for M-loop XMU = XM_US(1) WRITE( *,* ) ' ........... DOWN ............... ' write( *,* ) ' Stocks down tau = ', tau_l write( *,* ) ' Mu = ', DACOS(XMU)/DACOS(-1.0d0)*180.0d0 write( *,* ) ' Mu0 = ', DACOS(MU0)/DACOS(-1.0d0)*180.0d0 write( *,* ) ' Relazm = ', Relazm/DACOS(-1.0d0)*180.0d0 WRITE( *,* ) $ ' Numeric SS-analytic' c ----- Analytical solution for SS radiation down T1 = OMEG*MU0/4.d0/(XMU - MU0)* $ ( DEXP(-TAU_L/XMU) - DEXP(-TAU_L/MU0) ) CALL Z_MATR_SS ( XMU, MU0, RELAZM ) SS_SOL_D(:,:) = T1*MATMUL( Z_matr,F_INC ) DO I = 1,4 WRITE( *,'(2x,i2,3(2x,ES15.6))') $ I, FIN_DOWN_PP(1,1,1,I), SS_SOL_D(I,1) ENDDO c ..... Analytical solution for SS up-going radiation c for all solar and line-of-sight angles DO KK = 1,N_T0 TAU_L = T0_LAY(KK) DO K1 = 1,NS_US XMU = -XM_US(K1) DO KS = 1,N_MU0 MU0 = XM_MU0(KS) T1 = OMEG*MU0/4.d0/(-XMU + MU0)* $ (1.d0 - DEXP( -TAU_L*(1.d0/MU0 - 1.d0/XMU)) ) CALL Z_MATR_SS ( XMU, MU0, RELAZM ) SS_SOL_U( K1,KS,KK,:,: ) = T1*MATMUL( Z_matr,F_INC ) ENDDO ENDDO ENDDO IF( SS_COR ) THEN FIN_UP_PP(:,:,:,1:4) = FIN_UP_PP(:,:,:,1:4) + $ SS_SOL_U( :,:,:,1:4,1 ) ENDIF WRITE( *,* ) $ ' ********* Output for all user angles ' WRITE( *,* ) $' MU I Q U V' DO KK = 1, N_T0 DO KS = 1,N_MU0 WRITE( *,* ) ' ------ Solar zenith angle : --------' WRITE( *,'(10x,f6.2)') DACOS(XM_MU0(KS))/3.14159*180. DO I = 1,NS_US WRITE( *,'(2x,f6.2,4(2x,ES15.5),2x,es12.4)') $ -XM_US(I),(FIN_UP_PP(I,KS,KK,J), J = 1,4),T0_LAY(KK) ENDDO DO I = NS_US,1,-1 WRITE( *,'(2x,f6.2,4(2x,ES15.5))') $ XM_US(I), (FIN_DOWN_PP(I,KS,KK,J), J = 1,4) ENDDO WRITE( *,* ) ' ' ENDDO ENDDO c ********* Output results ********** IF( KEY_OUT .EQ. 'REFLEC' ) THEN OPEN( 45,FILE = TRIM(OUT_NAME_MS) ) DO KK = 1, N_T0 DO I = 1,NS_US DO KS = 1,N_MU0 T1 = DACOS(XM_MU0(KS))/3.14159*180. S1(:) = FIN_UP_PP(I,KS,KK,:)/XM_MU0(KS) S2(1) = 100.*S1(2)/S1(1) S2(2) =100.* DSQRT( S1(2)**2 + S1(3)**2 )/S1(1) S2(3) = 100.*S1(4)/S1(1) WRITE( 45,'(1x,e12.3,7(1x,ES15.5),3(1x,f6.2))') $ T0_LAY(KK),(S1(J), J = 1,4), (S2(J), J = 1,3), $ -XM_US(I), T1, AZIM ENDDO ENDDO ENDDO close(45) OPEN( 45,FILE = TRIM(OUT_NAME_SS)) DO KK = 1, N_T0 DO I = 1,NS_US DO KS = 1,N_MU0 T1 = DACOS(XM_MU0(KS))/3.14159*180. S1(:) = SS_SOL_U(I,KS,KK,:,1)/XM_MU0(KS) S2(1) = 100.*S1(2)/S1(1) S2(2) =100.* DSQRT( S1(2)**2 + S1(3)**2 )/S1(1) S2(3) = 100.*S1(4)/S1(1) WRITE( 45,'(1x,e12.3,7(1x,ES15.5),3(1x,f6.2))') $ T0_LAY(KK),(S1(J), J = 1,4), (S2(J), J = 1,3), $ -XM_US(I), T1, AZIM ENDDO ENDDO ENDDO close(45) ENDIF IF( KEY_OUT .EQ. 'TRANSM' ) THEN OPEN( 45,FILE='Transm.dat') DO KK = 1, N_T0 DO I = 1,NS_US DO KS = 1,N_MU0 T1 = DACOS(XM_MU0(KS))/3.14159*180. S1(:) = FIN_DOWN_PP(I,KS,KK,:)/XM_MU0(KS) S2(1) = 100.*S1(2)/S1(1) S2(2) =100.* DSQRT( S1(2)**2 + S1(3)**2 )/S1(1) S2(3) = 100.*S1(4)/S1(1) WRITE( 45,'(1x,e12.3,7(1x,ES15.5),3(1x,f6.2))') $ T0_LAY(KK),(S1(J), J = 1,4), (S2(J), J = 1,3), $ XM_US(I), T1, AZIM ENDDO ENDDO ENDDO close(45) ENDIF CONTAINS c ************************************************************************** SUBROUTINE FILEREAD INTEGER NNS, NLEG_REAL, NNS_US, I_DUM INTEGER KEY_SOL CHARACTER (LEN = 80) INP_NAME CHARACTER (LEN = 1) DIM_ANG DOUBLE PRECISION P_I, S_a, S_e, STEP P_I = DACOS(-1.0d0) OPEN( 44,file = 'vector.inp', status = 'old' ) CALL MOVE_POINTER READ( 44,* ) INP_NAME WRITE( *,* ) ' ---->> Filename for Legendre coefficients' WRITE( *,* ) TRIM(INP_NAME) CALL MOVE_POINTER READ(44,*) NS WRITE( *,* ) ' ---->> Numbe of Gauss angles: ' WRITE( *,'(5x,I3)') NS CALL MOVE_POINTER READ(44,*) NLEG_REAL WRITE( *,* ) ' ---->> Numbe of Legendre coefficients: ' WRITE( *,'(5x,I3)') NLEG_REAL CALL MOVE_POINTER READ(44,*) SS_COR CALL MOVE_POINTER READ(44,*) DIM_ANG WRITE( *,* ) ' ---->> Input angles dimension: ' WRITE( *,'(7x,a1)') DIM_ANG READ(44,*) NS_US WRITE( *,* ) ' ---->> Number of user angles' ALLOCATE( XM_US(NS_US) ) READ(44,*) ( XM_US(I), I = 1,NS_US ) IF( DIM_ANG .EQ. 'a' ) THEN WRITE( *,* ) ' **** Angles are recalculated in cosines ' XM_US(:) = DCOS(XM_US(:)/180.0d0*P_I) ENDIF DO I = 1, NS_US WRITE(*,'(5x,I3,2x,e12.4,2x,f7.3)') $ I, XM_US(I), DACOS(XM_US(I))/P_I*180. ENDDO CALL MOVE_POINTER READ(44,*) AZIM WRITE( *,* ) ' ---->> Relative azimuth :' WRITE( *,'(5x,f7.2)' ) AZIM CALL MOVE_POINTER READ(44,*) OMEG WRITE( *,* ) ' ---->> Single-scattering albedo :' WRITE( *,'(5x,E15.7)' ) OMEG CALL MOVE_POINTER READ(44,*) N_T0 ALLOCATE( T0_LAY(N_T0) ) READ(44,*) (T0_LAY(I),I = 1,N_T0) WRITE( *,* ) ' ---->> Optical thickness of the layer :' DO I = 1,N_T0 WRITE( *,'(5x,I3,2x,ES15.7)' ) I,T0_LAY(I) ENDDO CALL MOVE_POINTER READ(44,*) KEY_SOL READ(44,*) DIM_ANG WRITE( *,* ) ' ---->> Input solar angles dimension: ' WRITE( *,'(7x,a1)') DIM_ANG IF( KEY_SOL .EQ. 1 ) THEN WRITE(*,*) ' ---->> Solar angle input scheme 1 ' READ(44,*) N_MU0 ALLOCATE( XM_MU0(N_MU0) ) READ(44,*) (XM_MU0(I), I = 1,N_MU0) ENDIF IF( KEY_SOL .EQ. 2 ) THEN WRITE(*,*) ' ---->> Solar angle input scheme 2 ' READ(44,*) S_a, S_e, STEP N_MU0 = INT( (S_e - S_a)/STEP ) + 1 ALLOCATE( XM_MU0(N_MU0) ) XM_MU0(1) = S_a DO I = 1, N_MU0 - 1 XM_MU0(I+1) = XM_MU0(I) + STEP ENDDO ENDIF IF( DIM_ANG .EQ. 'a' ) THEN WRITE( *,* ) ' **** Angles are recalculated in cosines ' XM_MU0(:) = DCOS(XM_MU0(:)/180.0d0*P_I) ENDIF WRITE( *,* ) ' ---->> Cosines and the solar angles:' DO I = 1, N_MU0 WRITE(*,'(5x,I3,2x,e12.4,2x,f6.3)') $ I, XM_MU0(I), DACOS(XM_MU0(I))/P_I*180. ENDDO CALL MOVE_POINTER READ(44,*) ALBEDO WRITE( *,* ) ' ---->> Albedo :' WRITE( *,'(5x,E15.7)' ) ALBEDO CALL MOVE_POINTER READ(44,*) KEY_OUT WRITE( *,* ) ' ---->> Output mode :' WRITE( *,'(5x,a6)') KEY_OUT CALL MOVE_POINTER READ(44,*) TAU_US WRITE( *,* ) ' ---->> Output optical depth :' IF( KEY_OUT .EQ. 'INSIDE' ) $ WRITE( *,'(5x,E15.7)' ) TAU_US IF( KEY_OUT .EQ. 'TRANSM' ) $ WRITE( *,* ) ' bottom of the layer' IF( KEY_OUT .EQ. 'REFLEC' ) $ WRITE( *,* ) ' top of the layer' CALL MOVE_POINTER READ( 44,* ) OUT_NAME_MS WRITE( *,* ) ' ---->> Output filename for MS results ' WRITE( *,* ) TRIM(OUT_NAME_MS) CALL MOVE_POINTER READ( 44,* ) OUT_NAME_SS WRITE( *,* ) ' ---->> Output filename for SS results ' WRITE( *,* ) TRIM(OUT_NAME_SS) CLOSE( 44 ) IF( KEY_OUT .NE. 'REFLEC' ) SS_COR = .FALSE. WRITE( *,* ) ' ---->> Single-scattering correction ' WRITE( *,* ) SS_COR RELAZM = AZIM/180.0d0 * P_I c ... Number of scattering angles for single-scattering correction N_SCA = N_MU0 + NS_US - 1 ALLOCATE( Greek(5000,6) ) OPEN( 44,file = TRIM(INP_NAME), status = 'old' ) DO I = 1,12 ; READ(44,*) ; ENDDO I = 0 100 I = I + 1 READ( 44,*,END = 200 ) I_DUM, (Greek(I,J), J = 1,6) IF( Greek(I,1) < 1.0d-06 ) GOTO 200 GO TO 100 200 CONTINUE WRITE( *,* ) ' **** Number of LP from file:', I-1 NLEG = I - 1 ; NLEG_MAX = NLEG CLOSE( 44 ) IF( NLEG_REAL < NLEG ) NLEG = NLEG_REAL ALLOCATE( COF_11(1:700), $ DWF_P2U(NLEG,NS_US), DWF_M2U(NLEG,NS_US), $ DWF_P2_M(NLEG_MAX,NS_US), DWF_M2_M(NLEG_MAX,NS_US), $ ALP_U(NLEG,NS_US), CF(NLEG,NLEG), $ DWF_P2(NLEG,NS), DWF_M2(NLEG,NS), $ B_LU(NLEG,4,4), B_LD(NLEG,4,4), ALP(NLEG,NS), $ B_LUM(NLEG_MAX,4,4), B_LDM(NLEG_MAX,4,4), $ ALP0(NLEG,N_MU0), TTT(NLEG,NS), $ PI_LM(NLEG,4*NS,4), GLM(NLEG,4*NS,4), $ XM(NS), AG(NS), PIE(4*NS,4), $ MH(8*NS,8*NS), RHS(8*NS,N_MU0), $ AT_COF(4*NS,N_MU0), BT_COF(4*NS,N_MU0), $ Part_p(4*NS), Part_m(4*NS), $ aj_as(4*NS,N_MU0), bj_as(4*NS,N_MU0), $ H4(4*NS), H3(4*NS), SCAT_ANG(N_SCA), $ QP1_AS(4*NS,N_MU0), QM1_AS(4*NS,N_MU0), $ FIP(4*NS,4*NS), FIM(4*NS,4*NS), HEL(4,4*NS), $ PSIP(4*NS,4*NS), PSIM(4*NS,4*NS), Norm(4*NS), $ SUM_E(4*NS,4*NS), SUM_F(4*NS,4*NS), SV(4*NS,4*NS), $ GLM_1(4*NS,4), GLM_2(4*NS,4), $ NJc(4*NS), NJr(4*NS), RR(4*NS), RI(4*NS), $ VR( 4*NS,4*NS ), VI( 4*NS,4*NS ), WORK(20*4*NS), $ IPIV(8*NS), TT(NS_US), $ SS_SOL_U ( NS_US,N_MU0,N_T0,4,1), $ FIN_UP_PP( NS_US,N_MU0,N_T0,4 ), $ FIN_DOWN_PP( NS_US,N_MU0,N_T0,4 ) ) B_LU = 0.0 ; B_LD = 0.0 DO I = 1, NLEG B_LU(I,1,1) = Greek(I,1) B_LU(I,2,2) = Greek(I,2) B_LD(I,3,3) = Greek(I,3) B_LD(I,4,4) = Greek(I,4) B_LU(I,1,2) = Greek(I,5) B_LU(I,2,1) = Greek(I,5) B_LD(I,3,4) = Greek(I,6) B_LD(I,4,3) = -Greek(I,6) ENDDO c DO K = 1,3 c write( *,* ) ' K = ',K c do i = 1,2 c write( *,'(2x,2(2x,e12.4))') ( B_LU(K,i,j), j=1,2 ) c enddo c write( *,* ) ' ' c do i = 3,4 c write( *,'(2x,2(2x,e12.4))') ( B_LD(K,i,j), j=3,4 ) c enddo c ENDDO B_LUM = 0.0 ; B_LDM = 0.0 DO I = 1, NLEG_MAX B_LUM(I,1,1) = Greek(I,1) B_LUM(I,2,2) = Greek(I,2) B_LDM(I,3,3) = Greek(I,3) B_LDM(I,4,4) = Greek(I,4) B_LUM(I,1,2) = Greek(I,5) B_LUM(I,2,1) = Greek(I,5) B_LDM(I,3,4) = Greek(I,6) B_LDM(I,4,3) = -Greek(I,6) ENDDO DEALLOCATE( Greek ) END SUBROUTINE FILEREAD c ********************************************************************* SUBROUTINE MOVE_POINTER CHARACTER (LEN = 80) TXT 1 read(44,*) TXT IF( TXT(1:1) .EQ. '#' ) GOTO 1 BACKSPACE( 44 ) RETURN END SUBROUTINE MOVE_POINTER c ******************************************************************** SUBROUTINE FORM_4T4 c Matrix M_FO (4x4) is for each angle XM(i), i = 1,2,..,NS, c and each L, L = M,M+1,...,NLEG-1 M_FO(1,1) = ALP(L,I) M_FO(4,4) = M_FO(1,1) M_FO(2,2) = 0.5*SIGN_M*( DWF_P2(L,I) + DWF_M2(L,I) ) M_FO(3,3) = M_FO(2,2) M_FO(2,3) =-0.5*SIGN_M*( DWF_P2(L,I) - DWF_M2(L,I) ) M_FO(3,2) = M_FO(2,3) END SUBROUTINE FORM_4T4 c ******************************************************************* SUBROUTINE MATRIX_QU M_FO = 0.0 ; SUM_E = 0.0 ; SUM_F = 0.0 QM1_AS = 0.0 ; QP1_AS = 0.0 DO L = M+1,NLEG,2 H4 = 0.0 DO I = 1,NS CALL FORM_4T4 PI_LM(L,4*I-3:4*I,:) = M_FO(:,:) PIE (4*I-3:4*I,:) = AG(I) * M_FO H4(4*I-3) = B_LU(L,1,1)*M_FO(1,1) H4(4*I-2) = B_LU(L,1,2)*M_FO(2,2) H4(4*I-1) = B_LU(L,1,2)*M_FO(3,2) ENDDO FORALL( K = 1:N_MU0 ) $ QP1_AS(:,K) = QP1_AS(:,K) + H4(:)*ALP0(L,K) FORALL( I = 1:NS ) H4(4*I-1) = -H4(4*I-1) FORALL( K = 1:N_MU0 ) QM1_AS(:,K) = QM1_AS(:,K) $ + (1-2*MODULO(L-1-M,2))*H4(:)*ALP0(L,K) HEL(:,:) = MATMUL( B_LU(L,:,:),TRANSPOSE(PIE) ) SUM_E = SUM_E + MATMUL( PI_LM(L,:,:),HEL ) HEL(:,:) = MATMUL( B_LD(L,:,:),TRANSPOSE(PIE) ) SUM_F = SUM_F + MATMUL( PI_LM(L,:,:),HEL ) ENDDO c WRITE( *,* ) ' NEXT ' DO L = M+2,NLEG,2 H4 = 0.0 DO I = 1,NS CALL FORM_4T4 PI_LM(L,4*I-3:4*I,:) = M_FO(:,:) PIE (4*I-3:4*I,:) = AG(I) * M_FO H4(4*I-3) = B_LU(L,1,1)*M_FO(1,1) H4(4*I-2) = B_LU(L,1,2)*M_FO(2,2) H4(4*I-1) = B_LU(L,1,2)*M_FO(3,2) ENDDO FORALL( K = 1:N_MU0 ) $ QP1_AS(:,K) = QP1_AS(:,K) + H4(:)*ALP0(L,K) FORALL( I = 1:NS ) H4(4*I-1) = -H4(4*I-1) FORALL( K = 1:N_MU0 ) QM1_AS(:,K) = QM1_AS(:,K) $ + (1-2*MODULO(L-1-M,2))*H4(:)*ALP0(L,K) HEL(:,:) = MATMUL( B_LD(L,:,:),TRANSPOSE(PIE) ) SUM_E = SUM_E + MATMUL( PI_LM(L,:,:),HEL ) HEL(:,:) = MATMUL( B_LU(L,:,:),TRANSPOSE(PIE) ) SUM_F = SUM_F + MATMUL( PI_LM(L,:,:),HEL ) ENDDO QP1_AS = 0.5*OMEG * QP1_AS QM1_AS = 0.5*OMEG * QM1_AS SUM_E = -OMEG*SUM_E SUM_F = -OMEG*SUM_F FORALL( I = 1:4*NS ) SUM_E(I,I) = 1.0 + SUM_E(I,I) SUM_F(I,I) = 1.0 + SUM_F(I,I) END FORALL DO I = 1,NS SUM_E(1:4*NS,4*I-3:4*I) = SUM_E(1:4*NS,4*I-3:4*I)/XM(I) SUM_F(1:4*NS,4*I-3:4*I) = SUM_F(1:4*NS,4*I-3:4*I)/XM(I) ENDDO RETURN END SUBROUTINE MATRIX_QU c ******************************************************************** SUBROUTINE FI_AND_PSI Jc = 0 ; Jr = 0 DO J = 1,4*NS IF( RI(J) > 0.0 ) THEN Jc = Jc + 1 ; NJc(Jc) = J c write(*,'(2(2x,i3),2(2x,e15.7))') Jc, J, RR(J), RI(J) ENDIF IF( RI(J) .EQ. 0.0 ) THEN Jr = Jr + 1 ; NJr(Jr) = J RR(J) = 1./DSQRT(RR(J)) H4 = MATMUL( SUM_E,VR(:,J) )*RR(J) FIP(J,:) = VR(:,J) + H4(:) FIM(J,:) = VR(:,J) - H4(:) H4 = MATMUL( TRANSPOSE(SUM_F),VI(:,J) )*RR(J) PSIP(J,:) = VI(:,J) + H4(:) PSIM(J,:) = -VI(:,J) + H4(:) Norm(J) = DOT_PRODUCT( VI(:,J),VR(:,J) ) ENDIF ENDDO c WRITE( *,* ) '... Number of complex-pairs :', Jc c WRITE( *,* ) '... Number of real-pairs :', Jr c IF( Jc .NE. 0 ) STOP c WRITE( *,* ) '... COMPLEX NORM : ' DO I = 1,Jc J = NJc(I) C_t1 = CMPLX( RR(J),RI(J) ) ; C_t2 = 1.0d0/sqrt(C_t1) RR(J) = DBLE( C_t2 ) ; RR(J+1) = AIMAG( C_t2 ) c T1 = DSQRT( RR(J)**2 + RI(J)**2 ) cc Imaginary and real part of the separation constant for complex cc conjugate eigenvalue J and J+1 are stored in the form: cc Nu = x + iy c RR(J+1) = -DSQRT( (T1 - RR(J))/2.d0 )/T1 c write(*,'(3(2x,e15.7))') RR(J), RI(J), T1 - RR(J) c RR(J) = DSQRT( (T1 + RR(J))/2.d0 )/T1 c write(*,'(2(2x,i3),2(2x,e15.7))') I, J, RR(J), RR(J+1) H3 = MATMUL( SUM_E,VR(:,J) ) H4 = MATMUL( SUM_E,VR(:,J+1) ) FIP(J,:) = VR(:,J) + (H3*RR(J) - H4*RR(J+1)) ! real part FIP(J+1,:) = VR(:,J+1) + (H4*RR(J) + H3*RR(J+1)) ! imaginary FIM(J,:) = VR(:,J) - (H3*RR(J) - H4*RR(J+1)) FIM(J+1,:) = VR(:,J+1) - (H4*RR(J) + H3*RR(J+1)) c PSI - function H3 = MATMUL( TRANSPOSE(SUM_F),VI(:,J) ) H4 = -MATMUL( TRANSPOSE(SUM_F),VI(:,J+1) ) PSIP(J,:) = VI(:,J) + (H3*RR(J) - H4*RR(J+1)) ! real part PSIP(J+1,:) = -VI(:,J+1) + (H4*RR(J) + H3*RR(J+1)) ! imaginary PSIM(J,:) = -VI(:,J) + (H3*RR(J) - H4*RR(J+1)) PSIM(J+1,:) = VI(:,J+1) + (H4*RR(J) + H3*RR(J+1)) Norm(J) = DOT_PRODUCT( VI(:,J), VR(:,J) ) $ + DOT_PRODUCT( VI(:,J+1),VR(:,J+1) ) Norm(J+1) = DOT_PRODUCT( VI(:,J), VR(:,J+1) ) $ - DOT_PRODUCT( VI(:,J+1),VR(:,J) ) ENDDO END SUBROUTINE FI_AND_PSI c ********************************************************************* SUBROUTINE COEF_FOR_PART_SOL( QM1,QP1, aj,bj ) DOUBLE PRECISION QM1(4*NS), QP1(4*NS), aj(4*NS), bj(4*NS) c Each element 4*I (I=1..NS) for QP1 and QM1 equals zero H3 = QM1 FORALL ( I = 1:NS ) H3( 4*I-1 ) = -H3( 4*I-1 ) DO I = 1,Jr ; J = NJr(I) aj(J) = 0.5*( DOT_PRODUCT( PSIP(J,:),QP1 ) + $ DOT_PRODUCT( PSIM(J,:),H3 ) )/Norm(J) bj(J) = 0.5*( DOT_PRODUCT( PSIM(J,:),QP1 ) + $ DOT_PRODUCT( PSIP(J,:),H3 ) )/Norm(J) ENDDO DO I = 1,Jc ; J = NJc(I) T1= DOT_PRODUCT(PSIP(J,:) ,QP1)+DOT_PRODUCT(PSIM(J,:) ,H3) T2= DOT_PRODUCT(PSIP(J+1,:),QP1)+DOT_PRODUCT(PSIM(J+1,:),H3) C_t2 = CMPLX( NORM(J),NORM(J+1) ) ; C_t1 = CMPLX( T1,T2 ) aj(J) = DBLE(C_t1/C_t2)*0.5 ; aj(J+1) = AIMAG(C_t1/C_t2)*0.5 T1= DOT_PRODUCT(PSIM(J,:) ,QP1)+DOT_PRODUCT(PSIP(J,:), H3) T2= DOT_PRODUCT(PSIM(J+1,:),QP1)+DOT_PRODUCT(PSIP(J+1,:),H3) C_t1 = CMPLX( T1,T2 ) bj(J) = DBLE(C_t1/C_t2)*0.5 ; bj(J+1) = AIMAG(C_t1/C_t2)*0.5 ENDDO END SUBROUTINE COEF_FOR_PART_SOL c ************************************************************************ SUBROUTINE PART_SOLUTION ( XT, aj, bj ) Implicit none DOUBLE PRECISION XT, aj(4*NS), bj(4*NS) DOUBLE PRECISION E_t0tmu0, E_tmu0 COMPLEX (kind(0.0d0)) :: C_rr, C_aj, C_bj, C_res c WRITE( *,* ) ' TAU from PART_SOLUTION : ', XT E_tmu0 = DEXP( -XT/MU0 ) E_t0tmu0 = DEXP( -(TAU_L-XT)/MU0 ) Part_p = 0.0 ; Part_m = 0.0 c Real part of the particularly solution DO I = 1,Jr ; J = NJr(I) T3 = aj(J)*( DEXP(-XT/RR(J)) - E_tmu0 )/( RR(J) - MU0 ) T4 = bj(J)*(1.0 - DEXP( -(TAU_L-XT)/RR(J) )*E_t0tmu0 ) $ /(RR(J)+MU0)*E_tmu0 Part_p = Part_p + RR(J)*( T3*FIP(J,:) + T4*FIM(J,:) ) Part_m = Part_m + RR(J)*( T3*FIM(J,:) + T4*FIP(J,:) ) ENDDO c Imaginary part of the particularly solution DO I = 1,Jc ; J = NJc(I) C_aj = CMPLX( aj(J),aj(J+1) ) ; C_rr = CMPLX( RR(J),RR(J+1) ) C_res = C_rr*C_aj*( EXP(-CMPLX(XT)/C_rr) - E_tmu0 )/(C_rr-MU0) Part_p = Part_p + $ 2.0*( DBLE(C_res)*FIP(J,:) - AIMAG(C_res)*FIP(J+1,:) ) Part_m = Part_m + $ 2.0*( DBLE(C_res)*FIM(J,:) - AIMAG(C_res)*FIM(J+1,:) ) C_bj = CMPLX( bj(J),bj(J+1) ) C_res = C_bj*C_rr* $ ( 1.0 - EXP( -CMPLX(TAU_L-XT)/C_rr )*E_t0tmu0 ) $ /(C_rr + MU0)*E_tmu0 Part_p = Part_p + $ 2.0*( DBLE(C_res)*FIM(J,:) - AIMAG(C_res)*FIM(J+1,:) ) Part_m = Part_m + $ 2.0*( DBLE(C_res)*FIP(J,:) - AIMAG(C_res)*FIP(J+1,:) ) ENDDO Part_p(:) = Part_p(:)*MU0 ; Part_m(:) = Part_m(:)*MU0 c ... No change of sing if using for R.h.s. of the main matrix only !! c FORALL( I = 1:NS ) Part_m(4*I-1:4*I) = -Part_m(4*I-1:4*I) END SUBROUTINE PART_SOLUTION c ********************************************************************* SUBROUTINE MAIN_MATRIX DO I = 1,Jr ; J = NJr(I) RE_EXP = DEXP(-TAU_L/RR(J)) FORALL( K = 1:NS ) H3(K) = FIP(J,4*K-3)*AG(K)*XM(K) Raj = 2.0*ALBEDO*SUM( H3(1:NS) ) FORALL( K = 1:NS ) H3(K) = FIM(J,4*K-3)*AG(K)*XM(K) Rbj = 2.0*ALBEDO*SUM( H3(1:NS) ) H3(:) = FIM(J,:) FORALL( K = 1:NS ) H3(4*K-3) = H3(4*K-3) - Raj MH(1 :4*NS,2*J-1) = FIP(J,:) ! *Aj*FIP MH(4*NS+1:8*NS,2*J-1) = H3(:)*RE_EXP ! *Aj*(FIM - Raj) H3(:) = FIP(J,:) FORALL( K = 1:NS ) H3(4*K-3) = H3(4*K-3) - Rbj MH(1 :4*NS,2*J) = FIM(J,:)*RE_EXP ! *Bj*FIM MH(4*NS+1:8*NS,2*J) = H3(:) ! *Bj*(FIP - Rbj) ENDDO DO I = 1,Jc ; J = NJc(I) C_t1 = CMPLX(RR(J),RR(J+1)) ; C_t2 = EXP(-CMPLX(TAU_L)/C_t1) RE_EXP = DBLE( C_t2 ) ; IM_EXP = AIMAG( C_t2 ) FORALL( K = 1:NS ) H3(K) = $ (FIP(J,4*K-3)*RE_EXP - FIP(J+1,4*K-3)*IM_EXP)*AG(K)*XM(K) Raj_1 = 2.0*ALBEDO*SUM( H3(1:NS) ) FORALL( K = 1:NS ) H4(K) = $ (FIP(J,4*K-3)*IM_EXP + FIP(J+1,4*K-3)*RE_EXP)*AG(K)*XM(K) Raj_2 = 2.0*ALBEDO*SUM( H4(1:NS) ) FORALL( K = 1:NS ) $ H3(K) = (FIM(J,4*K-3) - FIM(J+1,4*K-3))*AG(K)*XM(K) Rbj_1 = 2.0*ALBEDO*SUM( H3(1:NS) ) FORALL( K = 1:NS ) H3(K) = $ (FIM(J,4*K-3) + FIM(J+1,4*K-3))*AG(K)*XM(K) Rbj_2 = 2.0*ALBEDO*SUM( H3(1:NS) ) H4(:) = FIM(J,:)*RE_EXP - FIM(J+1,:)*IM_EXP MH(1 :4*NS,2*J ) = H4(:) ! *Bj(1)*FM(1) FORALL( K = 1:NS ) H4(4*K-3) = H4(4*K-3) - Raj_1 MH(4*NS+1:8*NS,2*J-1) = H4(:) ! *Aj(1)*(FM(1)-Raj(1)) H3(:) = FIP(J,:) MH(1 :4*NS,2*J-1) = H3(:) ! *Aj(1)*FP(1) FORALL( K = 1:NS ) H3(4*K-3) = H3(4*K-3) - Rbj_1 MH(4*NS+1:8*NS,2*J ) = H3(:) ! *Bj(1)*(FP(1)-Rbj(1)) H4(:) = FIM(J,:)*IM_EXP + FIM(J+1,:)*RE_EXP MH(1 :4*NS,2*J+2) = H4(:) ! *Bj(2)*FM(2) FORALL( K = 1:NS ) H4(4*K-3) = H4(4*K-3) - Raj_2 MH(4*NS+1:8*NS,2*J+1) = H4(:) ! *Aj(2)*(FM(2)-Raj(2)) H3(:) = FIP(J+1,:) MH(1 :4*NS,2*J+1) = H3(:) ! *Aj(2)*FP(2) FORALL( K = 1:NS ) H3(4*K-3) = H3(4*K-3) - Rbj_2 MH(4*NS+1:8*NS,2*J+2) = H3(:) ! *Bj(2)*(FP(2)-Rbj(2)) ENDDO END SUBROUTINE MAIN_MATRIX c ********************************************************************* SUBROUTINE PP_UP ( XMU, TAU, NU, A_COF, B_COF, aj, bj ) IMPLICIT NONE INTEGER NU DOUBLE PRECISION XMU, TAU DOUBLE PRECISION A_COF(4*NS), B_COF(4*NS), aj(4*NS), bj(4*NS) INTEGER SIGN_ML DOUBLE PRECISION B_1(4), B_2(4), BM_FO(NLEG,4,4), BB_L(NLEG,4,4) DOUBLE PRECISION SH1(4,1), SH2(4,1) DOUBLE PRECISION T_t0_tm0m, C_t0_tnjm, T_t0_tnjm COMPLEX (kind(0.0d0)) :: C_bj, C_t2, C_t3, C_rr, C_aj, $ C_nuC, C_nuT, C_nuW, C_nuZ c Albedo = 0 only !!! E_t0m0 = DEXP(-TAU_L/MU0) IF( XMU .EQ. 0.0d0 ) THEN T1 = 0.0d0 T_t0_tm0m = 1.0d0/MU0*DEXP(-TAU/MU0) ELSE T1 = DEXP( -(TAU_L-TAU)/XMU ) T_t0_tm0m = ( 1.d0 - DEXP(-(TAU_L-TAU)/MU0) $ *DEXP(-(TAU_L-TAU)/XMU) )/( XMU + MU0 )*DEXP(-TAU/MU0) ENDIF B_L = 0.d0 ; M_FO = 0.d0 BM_FO = 0.0d0 ; BB_L = 0.0d0 DO L = M+1,NLEG BM_FO(L,1,1) = ALP_U(L,NU) BM_FO(L,4,4) = ALP_U(L,NU) BM_FO(L,2,2) = 0.5*SIGN_M*( DWF_P2U(L,NU) + DWF_M2U(L,NU) ) BM_FO(L,3,3) = 0.5*SIGN_M*( DWF_P2U(L,NU) + DWF_M2U(L,NU) ) BM_FO(L,2,3) =-0.5*SIGN_M*( DWF_P2U(L,NU) - DWF_M2U(L,NU) ) BM_FO(L,3,2) =-0.5*SIGN_M*( DWF_P2U(L,NU) - DWF_M2U(L,NU) ) BB_L(L,1:2,1:2) = B_LU(L,1:2,1:2) BB_L(L,3:4,3:4) = B_LD(L,3:4,3:4) ENDDO B_1 = 0.d0 DO J = 1,4*NS S1(:) = 0.0d0 ; S2(:) = 0.d0 DO L = M+1,NLEG SIGN_ML = ( 1 - 2*MODULO(L-1-M,2) ) FORALL( I = 1:4 ) $ B_1(I) = DOT_PRODUCT( BB_L(L,i,:),GLM(L,J,:) ) c B_1(:) = MATMUL( BB_L(L,:,:),GLM(L,J,:) ) FORALL( I = 1:4 ) $ B_2(I) = DOT_PRODUCT( BM_FO(L,i,:),B_1(:) ) S1(:) = S1(:) + B_2(:) c S1(:) = S1(:) + MATMUL( BM_FO(L,1:4,1:4),B_1(1:4) ) B_1(3:4) = -B_1(3:4) FORALL( I = 1:4 ) $ B_2(I) = DOT_PRODUCT( BM_FO(L,i,:),B_1(:) ) S2(:) = S2(:) + SIGN_ML*B_2(:) c S2 = S2 + SIGN_ML*MATMUL( BM_FO(L,:,:),B_1 ) ENDDO GLM_1(J,:) = S1(:) ; GLM_2(J,:) = S2(:) ENDDO IF( AG(1) .EQ. 0.0d0 ) GO TO 12 S1 = 0.0 ; S2 = 0.0 DO I = 1,Jr ; J = NJr(I) T2 = DEXP( -(TAU_L-TAU)/RR(J) ) IF( (RR(J) - XMU) .EQ. 0.d0 ) THEN C_t0_tnjm = T2*(TAU_L-TAU)/RR(J)/XMU ELSE C_t0_tnjm = (T2 - T1)/(RR(J) - XMU) ENDIF T_t0_tnjm = (1.0 - T2*T1 )/(RR(J) + XMU)*DEXP(-TAU/RR(J)) S1 = S1 + GLM_1(J,:)*RR(J)*( B_COF(J)*C_t0_tnjm $ + bj(J)*MU0*(MU0*T_t0_tm0m - RR(J)*C_t0_tnjm*E_t0m0) $ /(RR(J) + MU0) ) S2 = S2 + GLM_2(J,:)*RR(J)*( A_COF(J)*T_t0_tnjm $ + aj(J)*MU0*(RR(J)*T_t0_tnjm - MU0*T_t0_tm0m) $ /(RR(J) - MU0) ) ENDDO c ... Complex part of the solution DO I = 1,Jc ; J = NJc(I) C_rr = CMPLX( RR(J),RR(J+1) ) C_t2 = EXP( -CMPLX(TAU_L-TAU)/C_rr ) C_nuC = C_rr*(C_t2 - T1)/(C_rr - XMU) B_1(:) = B_COF(J) *GLM_1(J,:) + B_COF(J+1)*GLM_1(J+1,:) B_2(:) = B_COF(J+1)*GLM_1(J,:) - B_COF(J) *GLM_1(J+1,:) C_bj = CMPLX( bj(J),bj(J+1) ) C_nuW = MU0*C_rr*C_bj* $ (MU0*T_t0_tm0m - C_nuC*E_t0m0)/(C_rr + MU0) S1 = S1 + DBLE(C_nuC)*B_1 + AIMAG(C_nuC)*B_2 + 2.0* $ (DBLE(C_nuW)*GLM_1(J,:) - AIMAG(C_nuW)*GLM_1(J+1,:)) C_t3 = EXP( -CMPLX(TAU)/C_rr ) C_nuT = C_rr*(1.0 - C_t2*T1 )/(C_rr + XMU)*C_t3 B_1(:) = A_COF(J) *GLM_2(J,:) + A_COF(J+1)*GLM_2(J+1,:) B_2(:) = A_COF(J+1)*GLM_2(J,:) - A_COF(J) *GLM_2(J+1,:) C_aj = CMPLX( aj(J),aj(J+1) ) C_nuZ = MU0*C_rr*C_aj*(C_nuT - MU0*T_t0_tm0m) $ /(C_rr - MU0) S2 = S2 + DBLE(C_nuT)*B_1 + AIMAG(C_nuT)*B_2 + 2.0* $ (DBLE(C_nuZ)*GLM_2(J,:) - AIMAG(C_nuZ)*GLM_2(J+1,:)) ENDDO STOC_U = S1 + S2 c ... Test for single scattering 12 continue QPU = 0.d0 IF( SS_COR ) GOTO 13 S1(4) = 0.d0 DO L = M+1,NLEG SIGN_ML = ( 1 - 2*MODULO(L-1-M,2) ) S1(1) = B_LU(L,1,1)*ALP_U(L,NU) S1(2) = B_LU(L,1,2)*0.5*SIGN_M*(DWF_P2U(L,NU)+DWF_M2U(L,NU)) S1(3) = B_LU(L,1,2)*0.5*SIGN_M*(DWF_P2U(L,NU)-DWF_M2U(L,NU)) QPU = QPU + SIGN_ML*S1*ALP0(L,KS) ENDDO 13 CONTINUE S1(:) = STOC_U(:) ; S1(3:4) = -S1(3:4) STOC_U(:) = 0.5*OMEG*( MU0*T_t0_tm0m*QPU(:) + S1(:) ) END SUBROUTINE PP_UP c ********************************************************************* SUBROUTINE PP_DOWN ( XMU, TAU, NU, A_COF, B_COF, aj, bj ) IMPLICIT NONE INTEGER NU DOUBLE PRECISION XMU, TAU DOUBLE PRECISION A_COF(4*NS), B_COF(4*NS), aj(4*NS), bj(4*NS) DOUBLE PRECISION A_1(4), A_2(4) COMPLEX (kind(0.0d0)) :: C_nuY, C_bj, C_t2, C_t3, C_rr, C_aj, $ C_nuC, C_nuX, C_nuT IF( XMU .EQ. 0.0d0 ) THEN T1 = 0.0d0 ELSE T1 = DEXP( -TAU/XMU ) ENDIF E_t0m0 = DEXP(-TAU_L/MU0) IF( (MU0 - XMU) .EQ. 0.d0 ) THEN C_tm0m = DEXP(-TAU/MU0)*TAU/MU0/XMU ELSE C_tm0m = ( DEXP(-TAU/MU0) - T1 )/( MU0-XMU ) ENDIF STOC_D = 0.0 S1 = 0.0 ; S2 = 0.0 IF( AG(1) .EQ. 0.0d0 ) GO TO 12 DO I = 1,Jr ; J = NJr(I) T2 = DEXP( -TAU/RR(J) ) ; T3 = DEXP( -(TAU_L-TAU)/RR(J) ) IF( (RR(J) - XMU) .EQ. 0.d0 ) THEN C_tnjm = T2*TAU/RR(J)/XMU ELSE C_tnjm = (T2 - T1)/(RR(J) - XMU) ENDIF S1 = S1 + GLM_1(J,:)*RR(J)*( A_COF(J)*C_tnjm + $ aj(J)*MU0*(RR(J)*C_tnjm - MU0*C_tm0m)/(RR(J) - MU0) ) T_t0t_njm = (1.0 - T2*T1 )/(RR(J) + XMU)*T3 S2 = S2 + GLM_2(J,:)*RR(J)*( B_COF(J)*T_t0t_njm + $ bj(J)*MU0*(MU0*C_tm0m - RR(J)*E_t0m0*T_t0t_njm) $ /(RR(J) + MU0) ) ENDDO c ... Complex part of the solution DO I = 1,Jc ; J = NJc(I) C_rr = CMPLX( RR(J),RR(J+1) ) ; C_t2 = EXP(-CMPLX(TAU)/C_rr) C_nuC = C_rr*(C_t2 - T1)/(C_rr - XMU) A_1(:) = A_COF(J) *GLM_1(J,:) + A_COF(J+1)*GLM_1(J+1,:) A_2(:) = A_COF(J+1)*GLM_1(J,:) - A_COF(J) *GLM_1(J+1,:) C_aj = CMPLX( aj(J),aj(J+1) ) C_nuX = MU0*C_rr*C_aj*(C_nuC - MU0*C_tm0m)/(C_rr - MU0) S1 = S1 + DBLE(C_nuC)*A_1 + AIMAG(C_nuC)*A_2 + 2.0* $ (DBLE(C_nuX)*GLM_1(J,:) - AIMAG(C_nuX)*GLM_1(J+1,:)) C_t3 = EXP( -CMPLX(TAU_L-TAU)/C_rr ) C_nuT = C_rr*(1.0 - C_t2*T1 )/(C_rr + XMU)*C_t3 A_1(:) = B_COF(J) *GLM_2(J,:) + B_COF(J+1)*GLM_2(J+1,:) A_2(:) = B_COF(J+1)*GLM_2(J,:) - B_COF(J) *GLM_2(J+1,:) C_bj = CMPLX( bj(J),bj(J+1) ) C_nuY = MU0*C_rr*C_bj*(MU0*C_tm0m - E_t0m0*C_nuT) $ /(C_rr + MU0) S2 = S2 + DBLE(C_nuT)*A_1 + AIMAG(C_nuT)*A_2 + 2.0* $ (DBLE(C_nuY)*GLM_2(J,:) - AIMAG(C_nuY)*GLM_2(J+1,:)) ENDDO STOC_D = S1 + S2 12 continue QPU = 0.d0 ; S1(4) = 0.d0 DO L = M+1,NLEG S1(1) = B_LU(L,1,1)*ALP_U(L,NU) S1(2) = B_LU(L,1,2)*0.5*SIGN_M*(DWF_P2U(L,NU)+DWF_M2U(L,NU)) S1(3) =-B_LU(L,1,2)*0.5*SIGN_M*(DWF_P2U(L,NU)-DWF_M2U(L,NU)) QPU = QPU + S1*ALP0(L,KS) ENDDO STOC_D(:) = 0.5*OMEG*STOC_D(:) STOC_D(:) = STOC_D(:) + 0.5*OMEG*MU0*C_tm0m*QPU(:) END SUBROUTINE PP_DOWN c ********************************************************************** SUBROUTINE Z_MATR_SS ( XMU, MU0, RELAZM ) Implicit none DOUBLE PRECISION XMU, MU0, RELAZM DOUBLE PRECISION SCAT, T1, T2, PL_RES(NLEG_MAX,NS_US) DOUBLE PRECISION Rot_r(4,4), Rot_l(4,4), Ph_matr(4,4) LOGICAL No_Rotate SCAT = XMU*MU0 + DSQRT( (1.d0-XMU**2)*(1.d0-MU0**2) )*DCOS(RELAZM) Rot_r = 0.d0 ; Rot_l = 0.d0 No_Rotate = .FALSE. IF( DABS(xmu) .EQ. 1.d0 .OR. mu0 .EQ. 1.d0 .OR. RELAZM .EQ. 0.d0 ) $ No_Rotate = .TRUE. IF( No_Rotate ) THEN c WRITE( *,* ) ' ------- No rotation --------- ' FORALL( I = 1:4 ) Rot_r(I,I) = 1.d0 ; Rot_l(I,I) = 1.d0 ENDFORALL ELSE c ... Rotation matrixes t1 = ( xmu - mu0*scat )/ dsqrt( (1.d0-mu0**2)*(1.d0-scat**2) ) t2 = ( mu0 - xmu*scat )/ dsqrt( (1.d0-xmu**2)*(1.d0-scat**2) ) Rot_r(1,1) = 1.d0 ; Rot_r(4,4) = 1.d0 Rot_r(2,2) = 2.d0*t1**2 - 1.d0 ; Rot_r(3,3) = Rot_r(2,2) Rot_r(3,2) = 2.d0*dsqrt( DABS(1.d0 - t1**2) )*t1 Rot_r(2,3) = -Rot_r(3,2) Rot_l(1,1) = 1.d0 ; Rot_l(4,4) = 1.d0 Rot_l(2,2) = 2.d0*t2**2 - 1.d0 ; Rot_l(3,3) = Rot_l(2,2) Rot_l(3,2) = 2.d0*dsqrt( DABS(1.d0 - t2**2) )*t2 Rot_l(2,3) = -Rot_l(3,2) ENDIF TT(:) = scat CALL WIGNER_PLG( TT, NS_US, 0, NLEG_MAX, PL_RES ) Ph_matr = 0.d0 Ph_matr(1,1) = DOT_PRODUCT( PL_RES(:,1),B_LUM(:,1,1) ) CALL WIGNER_FUNCTION( TT, NS_US, 2, NLEG_MAX, DWF_P2_M, DWF_M2_M ) t1 = 0.d0 ; t2 = 0.d0 DO J = 1,NLEG_MAX t1 = t1 + (B_LUM(J,2,2) + B_LDM(J,3,3))*DWF_P2_M(J,1) t2 = t2 + (B_LUM(J,2,2) - B_LDM(J,3,3))*DWF_M2_M(J,1) ENDDO Ph_matr(2,2) = 0.5*( t1 + t2 ) Ph_matr(3,3) = 0.5*( t1 - t2 ) Ph_matr(4,4) = DOT_PRODUCT( PL_RES(:,1),B_LDM(:,4,4) ) CALL WIGNER_FUNCTION( TT, NS_US, 0, NLEG_MAX, DWF_P2_M, DWF_M2_M ) T1 = 0.d0 DO J = 1,NLEG_MAX T1 = T1 + B_LUM(J,1,2)*DWF_P2_M(J,1) ENDDO Ph_matr(1,2) = T1 Ph_matr(2,1) = T1 Z_matr = MATMUL( Rot_l,MATMUL( Ph_matr,Rot_r ) ) RETURN END SUBROUTINE Z_MATR_SS c **************************************************************** Subroutine GAULEG ( X1, X2, X, W, N ) C Taken from Numerical Recipes chapter 4 C Returns Gauss Legendre abscissa and weights for range [X1,X2] INTEGER N DOUBLE PRECISION X(N), W(N), X1, X2 C Local variables DOUBLE PRECISION YEPS PARAMETER ( YEPS = 3.0D-14 ) INTEGER M, I, J DOUBLE PRECISION YXM, YXL, YZ, YP1, YP2, YP3, YPP, YZ1 C ************************************* M = (N+1)/2 YXM = 0.5D0 * (X2+X1) YXL = 0.5D0 * (X2-X1) DO 10 I = 1, M YZ = DCOS( 3.141592654D0 * (I-0.25D0) / (N+0.5D0) ) 1 CONTINUE YP1 = 1.0D0 YP2 = 0.0D0 DO 20 J = 1, N YP3 = YP2 YP2 = YP1 YP1 = ( (2.0D0*J-1.0D0) * YZ * YP2 - (J-1.0D0) * YP3 ) / J 20 CONTINUE YPP = N * (YZ*YP1-YP2) / (YZ*YZ-1.0D0) YZ1 = YZ YZ = YZ1 - YP1 / YPP IF ( DABS(YZ-YZ1) .GT. YEPS ) GOTO 1 X(I) = YXM - YZ * YXL X(N+1-I) = YXM + YXL * YZ W(I) = 2.0D0 * YXL / ( (1.0D0-YZ*YZ) * YPP * YPP ) W(N+1-I) = W(I) 10 CONTINUE END SUBROUTINE END ! End of program c ******************************************************************** SUBROUTINE WIGNER_FUNCTION( XM, NS, M, NLEG, DWF_P2, DWF_M2 ) Implicit none c Input variables c --------------- INTEGER NS,M,NLEG DOUBLE PRECISION XM(NS) c Output varable c --------------- DOUBLE PRECISION DWF_P2(NLEG,NS),DWF_M2(NLEG,NS) c Local variables c --------------- INTEGER S,IS,S0, I, L DOUBLE PRECISION NFAC DOUBLE PRECISION C1, C2 DOUBLE PRECISION HM(NS), HP(NS) HM(:) = 1.0d0 - XM(:) ; HP(:) = 1.0d0 + XM(:) S0 = MAX(M,2) DWF_P2(1:S0,:) = 0.0d0 ; DWF_M2(1:S0,:) = 0.0d0 IF( M .EQ. 0 ) THEN DWF_P2(S0+1,:) = DSQRT(6.D0)/4.*HM(:)*HP(:) DWF_M2(S0+1,:) = DWF_P2(S0+1,:) ENDIF IF( M .EQ. 1 ) THEN DWF_P2(S0+1,:) = 0.5*DSQRT(HM(:)*HP(:))*HP(:) DWF_M2(S0+1,:) = 0.5*DSQRT(HP(:)*HM(:))*HM(:) ENDIF IF( M .EQ. 2 ) THEN DWF_P2(S0+1,:) = HP(:)*HP(:)*0.25d0 DWF_M2(S0+1,:) = HM(:)*HM(:)*0.25d0 ENDIF IF( M .GT. 2 ) THEN NFAC = 0.5d0*DBLE(M*(M-1)) DO I = 3,M ; NFAC = NFAC*(M+I)/I ; ENDDO C1 = (1-2*MODULO(M,2))*DSQRT(NFAC)/2.0d0**M HM(:) = dsqrt(HM(:)) ; HP(:) = dsqrt(HP(:)) DWF_P2(S0+1,:) = C1 * HM(:)**(M-2)*HP(:)**(M+2) DWF_M2(S0+1,:) = C1 * HP(:)**(M-2)*HM(:)**(M+2) ENDIF DO L = S0+1,NLEG - 1 C1 = 1.0d0/(L-1)/DSQRT( DBLE((L-M)*(L+M)*(L+2)*(L-2)) ) C2 = L*DSQRT( DBLE((L-1-M)*(L-1+M)*(L-3)*(L+1)) ) DWF_P2(L+1,:) = C1* $ ((2*L-1)*(L*(L-1)*XM(:)-2*M)*DWF_P2(L,:)-C2*DWF_P2(L-1,:)) DWF_M2(L+1,:) = -C1* $ ((2*L-1)*(L*(L-1)*XM(:)+2*M)*DWF_M2(L,:)+C2*DWF_M2(L-1,:)) ENDDO DO L = 0,NLEG-1 DWF_M2(L+1,:) = (1 - 2*MODULO(L+M,2)) * DWF_M2(L+1,:) ENDDO RETURN END c **************************************************************** SUBROUTINE WIGNER_PLG( XM, NS, M, NLEG, RES ) Implicit none c Input variables c --------------- INTEGER NS,M,NLEG DOUBLE PRECISION XM(NS) c Output varable c --------------- DOUBLE PRECISION RES(NLEG,NS) c Local variables c --------------- INTEGER I, L DOUBLE PRECISION MFAC DOUBLE PRECISION HM(NS) IF( M .EQ. 0 ) THEN RES(1,:) = 1.0d0 ; RES(2,:) = XM(:) DO L = 2,NLEG - 1 RES(L+1,:) = ((2*L-1)*XM(:)*RES(L,:) - (L-1)*RES(L-1,:))/L ENDDO RETURN ENDIF RES(1:M,:) = 0.0d0 ; MFAC = 1.d0 DO I = 1,M ; MFAC = MFAC*(M+I)/I ; ENDDO HM(:) = DSQRT( (1.0d0 - XM(:))*(1.0d0 + XM(:)) )/2.0d0 c ... Without (-1)**M as has been used by Siewert RES(M+1,:) = DSQRT(MFAC)*HM(:)**M DO L = M+1,NLEG - 1 RES(L+1,:) = ( (2*L-1)*XM(:)*RES(L,:) $ - DSQRT(DBLE((L-1-M)*(L-1+M)))*RES(L-1,:) ) $ /DSQRT(DBLE((L-M)*(L+M))) ENDDO RETURN END c ****************************************************************