C Last change: GL 5 Oct 2004 4:54 pm PROGRAM DOUBLE PARAMETER(NZM=32,NA=4*NZM,NZ=NZM,KOM=140,NZ2=NZ+2) REAL W(NZ),ANGLE(NZ),OM(KOM),EWORK(NZ),WKAREA(NA),P(NZ,NZ), 1 R(NZ,NZ),M(NZ,NZ),B(NZ,NZ),FGLOBE(2),MIE_COEF(129), 2 T(NZ,NZ),C(NZ,NZ),VP(NZ),DOWNIN(NZ),UPINT(NZ),VLN(NZ,2), 3 VLP(NZ,2),FU(2),FD(2),VN(NZ),PO(NZ),BO(NZ),VPOLD(NZ), 5 RSURF(NZ,NZ),PL(KOM,NZ2,NA),VLAZ(19,32) INTEGER PHI c c DATA ANGLE/0.0483077,0.144472,0.2392874,0.3318686,0.4213513, c A 0.5068999,0.5877158,0.6630443,0.7321821,0.7944838, cc B 0.8493676,0.8963212,0.9349061,0.9647623,0.9856115,0.9972639/ c DATA W/0.0965401,0.0956387,0.0938444,0.0911739,0.0876521, c A 0.0833119,0.0781939,0.0723458,0.0658222,0.0586841,0.0509981, c B 0.0428359,0.0342739,0.0253921,0.0162744,0.0070186/ c DATA ANGLE/0.0950125098376,0.281603550779,0.45801677766, c 1 0.617876244403,0.75540440836,0.8656312023878, c 2 0.94457502307323,0.989400934992/ c DATA W/0.1894506104551,0.1826034150449,0.16915651940, c 1 0.1495959888166,0.1246289712555,0.0951585116825, c 1 0.06225352393865,0.0271524594118/ C DATA W/0.2955242247,0.2692667193,0.2190863625,0.1494513492, C 1 0.0666713443/ C DATA ANGLE/0.148874339,0.4333953941,0.6794095683,0.8650633667, C 1 0.9739065285/ DATA ALPHA/1.00/ DATA FO/1./ * C PRINT *, 'ANGLEO = ,G= ,SURF ALB=' C READ *, ANGLEO,G,SRALB ANGLEO=1.0 G=0.85 SRALB=0.0 NAA=1 NZM1=NZM+2 NN=KOM PI=ACOS(-1.) G1=0.9 G2=-0.7 G=0.85 ! OPEN(2,file='coef_1_10.dat') ! OPEN(2,file='coef_1_40.dat') ! do k=1,129 ! READ(2,*) MIE_COEF(k) ! end do ! CLOSE(2) DO 10 I=1,KOM c OM(I) = 0.9*G1**(I-1)+0.1*G2**(I-1) OM(I)=G**(I-1) c OM(I)=0. 10 CONTINUE ! f_scale= OM(NA+1) ! WRITE(*,*) f_scale ! do I=2,NA ! OM(I)=(OM(I)-f_scale)/(1.-f_scale) ! end do OM(1)=1 c OM(3)=0.5 CALL QUAD(NZ,ANGLE,W) WRITE(*,*) ANGLE WRITE(*,*) W CALL PLEG(ANGLE,ANGLEO,PL) * * This is a loop over azimuthal expansion: NAA specified * C DOMO=0.04 C DCT=0.15 IAA=1 C CT=-1 c CT=0.7 C NTAU=20 C NOMO=25 C DO 60 ITAU=1,NTAU C DT=10**CT C CT=CT+DCT C COMO=0. c write(*,*) DT,CT C DO 60 IOMO=1,NOMO C OMO=1.-COMO**2 C COMO=COMO+DOMO DT=10. OMO=1. ! DT=(1.-OMO*f_scale)*DT ! OMO=OMO*(1.-f_scale)/(1.-OMO*f_scale) c write(*,*) DT,OMO c if(ITAU.eq.18) write(*,'(''calling vscat'')') CALL VSCAT(ITHICK,PI,D,DT,IAA,W,OM,P,B,PO,BO,PL) c if(ITAU.eq.18) write(*,'(''calling isi'')') CALL ISI(OMO,PI,FO,D,ANGLEO,ANGLE,W, 1 PO,BO,P,B,R,T,VN,VP) c if(ITAU.eq.18) write(*,'(''calling dsi'')') CALL DSI(ITHICK,R,T,VN,VP,D,ANGLEO,DT,SRALB) c if(ITAU.eq.18) write(*,'(''calling asi'')') CALL ASI(R,T,M,EWORK,VN,VP, 2 SRALB,RSURF,ANGLE,W,IAA,VLN,VLP,ANGLEO,FO,DT) * * COMPUTE UPWARD AND DOWNWARD FLUXES FOR EACH LAYER * DO 21 L=1,2 TR = L-1 IF (ANGLEO .EQ. 0) THEN TR = 0 GO TO 11 END IF IF (TR*DT/ANGLEO .LT. 100.) THEN TR = ANGLEO*FO*DEXP(-DBLE(TR*DT/ANGLEO)) ELSE TR = 0 END IF 11 FU(L) = 0. FD(L) = 0. temp=0 DO 20,I=1,NZ temp=temp+2.*PI*W(I)*ANGLE(I)*TR*SRALB/PI FD(L) = FD(L)+2.*PI*W(I)*ANGLE(I)*VLP(I,L) FU(L) = FU(L)+2.*PI*W(I)*ANGLE(I)*VLN(I,L) 20 CONTINUE FGLOBE(L)=FD(L)+TR 21 CONTINUE WRITE(*,*) (PI*VLN(I,1)/FGLOBE(1),I=1,NZ) write(*,*) DT,OMO,FU(1),FU(2),FD(1),FGLOBE(2),TR 60 CONTINUE END SUBROUTINE VSCAT(ITHICK,PI,D,DT,IAA,W,OM,P,B,PO,BO,PL) PARAMETER(NZ=32,KOM=140,NA=4*NZ,NZ2=NZ+2,NA2=NZ*4) REAL P(NZ,NZ),B(NZ,NZ),W(NZ),OM(KOM),PO(NZ),BO(NZ), + PL(KOM,NZ2,NA) c WRITE(*,*) OM N1 = NZ+1 NP2 = NZ+2 PI2=PI*2. IA = IAA-1 DO 10, I=1,NZ DO 10,J=1,NZ SUM = 0. SUM1 = 0. U = 1 DO 9,K=IAA,NA2 KK = K-1 FK = 2*(K-1)+1 CONS = OM(K)*PL(K,I,IA+1)*PL(K,J,IA+1)*FK SUM = SUM+CONS SUM1 = SUM1+CONS*U U = -U 9 CONTINUE P(I,J) = SUM B(I,J) = SUM1 10 CONTINUE DO 12,I=1,NZ SUM = 0. SUM1 = 0. IF (IAA .EQ. 1) FAC = 1. U = 1 DO 11,K=IAA,NA2 FK = 2*(K-1)+1 KK = K-1 IF (IAA .GT. 1) FAC = 2 CONS = OM(K)*PL(K,I,IAA)*PL(K,N1,IAA)*FK*FAC SUM = SUM+CONS SUM1 = SUM1+CONS*U U = -U 11 CONTINUE PO(I) = SUM BO(I) = SUM1 12 CONTINUE ITHICK=ALOG(DT/0.0001)/ALOG(2.) + 1. IF(ITHICK.LE.0) ITHICK=1 D=DT/(2.**ITHICK) WRITE(*,*) D, ITHICK ******************************************************************** * THIS IS A DIAGNOSTIC TEST OF PHASE FUNCTION NORMALIZATION: * THE INTEGRAL OF THE PHASE FUNCTION (SUM BELOW) SHOULD EQUAL UNITY OPEN(2,'PB_matrices.txt') ! do 20 I=1,NZ SUM = 0. DO 15,J=1,NZ SUM = SUM+0.5*W(J)*(P(I,J)+B(I,J)) 15 CONTINUE * WRITE(*,101) SUM 100 FORMAT (1X,8F8.4) 101 FORMAT (1X,F8.6) 20 CONTINUE WRITE (2,100) ((P(I,J),J=1,NZ),I=1,NZ) WRITE (2,100) ((B(I,J),J=1,NZ),I=1,NZ) WRITE (2,100) (PO(J),J=1,NZ) WRITE (2,100) (BO(J),J=1,NZ) close(2) ********************************************************************** RETURN END SUBROUTINE PLEG(ANGLE,ANGLEO,PL) PARAMETER(NZM=32,NA=4*NZM,NZ=NZM,KOM=140,NZ2=NZ+2) ******************************************************************** * This subroutine calculatres the associated Legendre polys * used in VSCAT. Output=array PL ******************************************************************** REAL CM(NA),DM1(NA),ANGLE(NZ),LM1,LM2,MU,PL(KOM,NZ2,NA) CM(1) = 1. DM1(1) = 1. NZ1=NZ+1 NAM2=2*NA-2 DO 20,I=1,NZ1 IF (I .LE. NZ) MU = ANGLE(I) IF (I .EQ. NZ1) MU = ANGLEO DO 20,IM=1,NA M = IM-1 rm=M rm=0.5*rm * * COMPUTE INITIAL VALUES * IF (IM .EQ. 1) THEN PL(1,I,1) = 1. PL(2,I,1) = MU ELSE CM(IM) = -SQRT((2.*M-1.)/(2.*M))*CM(IM-1) PL(IM,I,IM) = CM(IM)*(1-MU**2)**rm DM1(IM) = -SQRT((2.*M+1.)/(2.*M))*DM1(IM-1) PL(IM+1,I,IM) = DM1(IM)*MU*(1-MU**2)**rm END IF DO 10,IL=1,NAM2 L = IL-1 IF (IL .LT. IM) THEN PL(IL,I,IM) = 0. GO TO 10 ELSE * * USE VARYING-DEGREE RECURRENCE FORMULA TO COMPUTE SUCCESSIVE * VALUES FOR L >= M. * LM2 = (L-M+2)*(L+M+2) LM1 = (L+M+1)*(L-M+1) PL(IL+2,I,IM) = ((2*L+3)*MU*PL(IL+1,I,IM)- $ SQRT(LM1)*PL(IL,I,IM))/SQRT(LM2) END IF 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE ISI(OMO,PI,FO,D,ANGLEO,ANGLE,W, 1 PO,BO,P,B,R,T,VN,VP) ****************************************************************** * This subroutine initalizes the integration procedure * and establishes the reflection(R), transmission (T)and * source functions (VN-up and VP-down) ********************************************************************* PARAMETER(NZ=32) REAL P(NZ,NZ),B(NZ,NZ),R(NZ,NZ),T(NZ,NZ),M(NZ,NZ),C(NZ,NZ), 1 PO(NZ),BO(NZ),VN(NZ),VP(NZ),ANGLE(NZ),W(NZ) IF (ANGLEO .EQ. 0) THEN TR = 0 GO TO 1 END IF IF (D/ANGLEO .LT. 100.) THEN TR = DEXP(-DBLE(D/ANGLEO)) ELSE TR = 0 END IF 1 TR=ANGLEO*(1.-TR)/D DO 6 I=1,NZ DO 5 J=1,NZ M(I,J)=0 5 C(I,J)=0 M(I,I)=1./ANGLE(I) 6 C(I,I)=W(I) CALL GMPRD(P,C,T,NZ,NZ,NZ) CALL GMPRD(B,C,R,NZ,NZ,NZ) CALL GMPRD(M,T,P,NZ,NZ,NZ) CALL GMPRD(M,R,B,NZ,NZ,NZ) DO 9 I=1,NZ DO 9 J=1,NZ P(I,J)=P(I,J)*OMO*0.5 B(I,J)=B(I,J)*OMO*0.5 DEL=0 IF(I.EQ.J) DEL=1 P(I,J)=DEL/ANGLE(I)-P(I,J) 9 CONTINUE c WRITE(*,*) ((P(I,J),J=1,NZ),I=1,NZ) c WRITE(*,*) ((B(I,J),J=1,NZ),I=1,NZ) DO 8 I=1,NZ DO 7 J=1,NZ R(I,J)=B(I,J)*D 7 T(I,J)=-P(I,J)*D T(I,I)=1.+T(I,I) VN(I)=.25*FO*D*OMO*BO(I)*M(I,I)*TR/PI 8 VP(I)=.25*FO*D*OMO*PO(I)*M(I,I)*TR/PI c WRITE(*,*) (VN(I),I=1,NZ) c WRITE(*,*) (VP(I),I=1,NZ) RETURN END SUBROUTINE DSI(ITHICK,R,T,VN,VP,D,ANGLEO,DT,SRALB) PARAMETER(NZ=32) ****************************************************************** * This is the doubling algorithm ******************************************************************* REAL R(NZ,NZ),T(NZ,NZ),RWORK(NZ,NZ),TWORK(NZ,NZ),TNEW(NZ,NZ), 1 EWORK(NZ),VN(NZ),VP(NZ),VPNEW(NZ),WKAREA(NZ) OPEN(3,'RT_matrices.txt') SPALB=0 IF (ANGLEO .EQ. 0) THEN TT = 0 GO TO 1 END IF IF (D/ANGLEO .LT. 100.) THEN TT = DEXP(-DBLE(D/ANGLEO)) CC=SPALB*DEXP(DBLE(D-DT)/ANGLEO) ELSE TT = 0 CC=0. END IF 1 CONTINUE DO 20 ILOOP=1,ITHICK CALL GMPRD(R,R,RWORK,NZ,NZ,NZ) DO 10 I=1,NZ DO 5 J=1,NZ 5 RWORK(I,J)=-RWORK(I,J) 10 RWORK(I,I)=1.+RWORK(I,I) CALL INVERT(NZ,RWORK) CALL GMPRD(T,RWORK,TWORK,NZ,NZ,NZ) CALL GMPRD(TWORK,R,RWORK,NZ,NZ,NZ) CALL GMPRD(RWORK,T,TNEW,NZ,NZ,NZ) CALL GMPRD(R,VN,EWORK,NZ,NZ,1) DO 11 I=1,NZ 11 EWORK(I)=EWORK(I)*TT+VP(I) CALL GMPRD(TWORK,EWORK,VPNEW,NZ,NZ,1) CALL GMPRD(R,VP,EWORK,NZ,NZ,1) DO 12 I=1,NZ EWORK(I)=EWORK(I)+VN(I)*TT 12 VPNEW(I)=VPNEW(I)+VP(I)*TT CALL GMPRD(TWORK,EWORK,VP,NZ,NZ,1) CALL GMADD(VP,VN,VN,NZ,1) CALL GMADD(R,TNEW,R,NZ,NZ) CALL GMPRD(TWORK,T,TNEW,NZ,NZ,NZ) DO 15 I=1,NZ VP(I)=VPNEW(I) DO 15 J=1,NZ T(I,J)=TNEW(I,J) 15 CONTINUE TT=TT*TT c WRITE(6,92) ILOOP 92 FORMAT(1H ,' ILOOP = ',I10) c WRITE(6,91) (VP(I),I=1,NZ) c WRITE(6,91) (VN(I),I=1,NZ) 20 CONTINUE IF(TT.eq.0.) then CC=0. else CC=CC/TT endif 91 FORMAT(1H ,4G15.5) write(3,99) ((R(I,J),J=1,NZ),I=1,NZ) WRITE(3,99) ((T(I,J),J=1,NZ),I=1,NZ) DO 13 I=1,NZ VP(I)=VPNEW(I)+VN(I)*CC VN(I)=VN(I)+VPNEW(I)*CC 13 CONTINUE 99 format(8f8.5) close(3) RETURN END SUBROUTINE ASI(R,T,RWORK,EWORK,VN,VP,SRALB, 2 RSURF,ANGLE,W,IAA,VLN,VLP,ANGLEO,FO,DT) PARAMETER(NZ=32) ******************************************************************* * This routine adds a lower reflecting layer to scattering layer * The reflection of this layer is specified by albedo SRALB ********************************************************************* REAL R(NZ,NZ),T(NZ,NZ),RWORK(NZ,NZ),VN(NZ),VP(NZ),EWORK(NZ), + UPINT(NZ),RSURF(NZ,NZ),ANGLE(NZ),W(NZ), + DOWNIN(NZ),VLN(NZ,2),VLP(NZ,2),V0(NZ),V1(NZ),V2(NZ) PI=ACOS(-1.) DO 2 I=1,NZ 2 DOWNIN(I)=0. DO 50 I=1,NZ V0(I)=0 DUMMY=2.*SRALB*ANGLE(I)*W(I) IF(IAA.EQ.1) V0(I)=SRALB*FO*ANGLEO*EXP(-DT/ANGLEO)/PI DO 50 J=1,NZ IF(IAA.GT.1) THEN RSURF(J,I)=0. GO TO 50 ENDIF RSURF(J,I)=DUMMY 50 CONTINUE CALL GMPRD(R,RSURF,RWORK,NZ,NZ,NZ) DO 52 I=1,NZ DO 51 J=1,NZ 51 RWORK(I,J)=-RWORK(I,J) 52 RWORK(I,I)=1.+RWORK(I,I) CALL INVERT(NZ,RWORK) CALL GMPRD(RWORK,VP,EWORK,NZ,NZ,1) CALL GMPRD(R,V0,V2,NZ,NZ,1) CALL GMPRD(RWORK,V2,V1,NZ,NZ,1) CALL GMADD(EWORK,V1,EWORK,NZ,1) CALL GMPRD(RSURF,EWORK,UPINT,NZ,NZ,1) CALL GMADD(UPINT,V0,UPINT,NZ,1) DO 53 I=1,NZ VLN(I,2)=UPINT(I) VLP(I,1)=DOWNIN(I) 53 DOWNIN(I)=EWORK(I) CALL GMPRD(T,UPINT,VP,NZ,NZ,1) CALL GMADD(VP,VN,UPINT,NZ,1) c CALL GMPRD(R,UPINT,DOWNIN,NZ,NZ,1) c CALL GMADD(DOWNIN,VP,DOWNIN,NZ,1) DO 65 I=1,NZ VLP(I,2)=DOWNIN(I) VLN(I,1)=UPINT(I) 65 CONTINUE RETURN END SUBROUTINE GMADD(A,B,R,N,M) DIMENSION A(1),B(1),R(1) NM=N*M DO 10 I=1,NM 10 R(I)=A(I)+B(I) RETURN END SUBROUTINE GMSUB(A,B,R,N,M) DIMENSION A(1),B(1),R(1) NM=N*M DO 10 I=1,NM 10 R(I)=A(I)-B(I) RETURN END SUBROUTINE GMPRD(A,B,R,N,M,L) DIMENSION A(1),B(1),R(1) IR=0 IK=-M DO 10 K=1,L IK=IK+M DO 10 J=1,N IR=IR+1 JI=J-N IB=IK R(IR)=0. DO 10 I=1,M JI=JI+N IB=IB+1 10 R(IR)=R(IR)+A(JI)*B(IB) RETURN END SUBROUTINE INVERT(N,A) * * MATRIX INVERSION USING GAUSS-JORDAN REDUCTION WITH * PARTIAL PIVOTING. * DIMENSION A(N,N),INTER(16,2) * * SEARCH FOR LARGEST PIVOT ELEMENT * DO 14,K=1,N JJ = K IF (K .EQ. N) GO TO 6 KP1 = K+1 BIG = ABS(A(K,K)) DO 5,I=KP1,N AB = ABS(A(I,K)) IF (BIG .GE. AB) GO TO 5 BIG = AB JJ = I 5 CONTINUE * * STORE NUMBERS OF ROWS INTERCHANGED. IF JJ=K, THERE IS * NO INTERCHANGE. * 6 INTER(K,1) = K INTER(K,2) = JJ IF (JJ .EQ. K) GO TO 9 * * ROW INTERCHANGE * DO 8,J=1,N STORK = A(K,J) STORJJ = A(JJ,J) A(K,J) = STORJJ A(JJ,J) = STORK 8 CONTINUE 9 CONTINUE * * CALCULATE NEW ELEMENTS OF PIVOT ROW EXCEPT FOR PIVOT * ELEMENT. * DO 10,J=1,N IF (J .EQ. K) GO TO 10 A(K,J) = A(K,J)/A(K,K) 10 CONTINUE * * CALCULATE NEW ELEMENT REPLACING PIVOT TERM. * A(K,K) = 1./A(K,K) * * CALCULATE NEW ELEMENTS NOT IN PIVOT ROW OR COLUMN. * DO 12,I=1,N IF (I .EQ. K) GO TO 12 DO 11,J=1,N IF (J .EQ. K) GO TO 11 A(I,J) = A(I,J)-A(K,J)*A(I,K) 11 CONTINUE 12 CONTINUE * * CALCULATE NEW ELEMENTS FOR PIVOT COLUMN, EXCEPT FOR * PIVOT ELEMENT. * DO 13,I=1,N IF (I .EQ. K) GO TO 13 A(I,K) = -A(I,K)*A(K,K) 13 CONTINUE 14 CONTINUE * * REARRANGE COLUMNS OF FINAL MATRIX. * DO 16,L=1,N K = N-L+1 KROW = INTER(K,1) IROW = INTER(K,2) IF (KROW .EQ. IROW) GO TO 16 DO 15,I=1,N STORI = A(I,IROW) STORK = A(I,KROW) A(I,IROW) = STORK A(I,KROW) = STORI 15 CONTINUE 16 CONTINUE RETURN END SUBROUTINE QUAD (NZ, ANGLE, W) C Generates the abscissas and weights for an even 2*NZ point C Gauss-Legendre quadrature. Only the NZ positive points are C returned. INTEGER NZ REAL ANGLE(*), W(*) INTEGER N, I, J, K, NUM PARAMETER (NUM=50) REAL*8 ABSCISSAS(NUM), WEIGHTS(NUM) REAL*8 Z, Z1, P1, P2, P3, PP, EPS PARAMETER (EPS=3.0D-14) N=2*NZ DO 12 I=1,NZ Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0)) K=0 1 CONTINUE P1=1.D0 P2=0.D0 DO 11 J=1,N P3=P2 P2=P1 P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J 11 CONTINUE PP=N*(Z*P1-P2)/(Z*Z-1.D0) Z1=Z Z=Z1-P1/PP K=K+1 IF (ABS(Z-Z1).GT.EPS .AND. K.LT.10) GO TO 1 ABSCISSAS(NZ+1-I)=Z WEIGHTS(NZ+1-I)=2.D0/((1.D0-Z*Z)*PP*PP) 12 CONTINUE DO 15 I=1,NZ ANGLE(I)=SNGL(ABSCISSAS(I)) W(I)=SNGL(WEIGHTS(I)) 15 CONTINUE RETURN END