MODULE RADCALV! Module wrapper for RadCal subroutineUSE PRECISION_PARAMETERSUSE GLOBAL_CONSTANTS, ONLY: PI,SQRTPIIMPLICIT NONEPRIVATECHARACTER(255), PARAMETER :: iradid='$Id: irad.f90 710 2007-09-28 19:52:12Z drjfloyd $'CHARACTER(255), PARAMETER :: iradrev='$Revision: 710 $'CHARACTER(255), PARAMETER :: iraddate='$Date: 2007-09-28 15:52:12 -0400 (Fri, 28 Sep 2007) $'PUBLIC OMMAX,OMMIN,DD,SPECIE,SVF,PLANCK,P,RCT,RCALLOC,INIT_RADCAL,RADCAL,RCDEALLOC,GET_REV_iradREAL(EB), ALLOCATABLE, DIMENSION(:,:) :: P,X,UUU,GC,GAMMA,SD15, SD,SD7,SD3REAL(EB), ALLOCATABLE, DIMENSION(:) :: SVF,XPART,RCT,SPECIE,TAU,TAUS,AC,AD,QW,TTAU,XTOT,XT,XSTAR,AB,PKPA,AMBDA,ATOT,BCNT,DDREAL(EB) OMMIN,OMMAX,TWALLINTEGER NPT,NPRINT,NOMCONTAINSSUBROUTINE GET_REV_irad(MODULE_REV,MODULE_DATE)INTEGER,INTENT(INOUT) :: MODULE_REVCHARACTER(255),INTENT(INOUT) :: MODULE_DATEWRITE(MODULE_DATE,'(A)') iradrev(INDEX(iradrev,':')+1:LEN_TRIM(iradrev)-2)READ (MODULE_DATE,'(I5)') MODULE_REVWRITE(MODULE_DATE,'(A)') iraddateEND SUBROUTINE GET_REV_iradSUBROUTINE INIT_RADCALNPT=1TWALL=0._EB!      TWALL=293.15IF(OMMAX<1100._EB) THEN    NOM=INT((OMMAX-OMMIN)/5._EB)ELSEIF(OMMIN>5000._EB) THEN   NOM=INT((OMMAX-OMMIN)/50._EB)ELSEIF(OMMIN<1100._EB.AND.OMMAX>5000._EB) THEN   NOM=INT((1100._EB-OMMIN)/5._EB)+INT((5000._EB-1100._EB)/25._EB) +INT((OMMAX-5000._EB)/50._EB)ELSEIF(OMMIN<1100._EB) THEN   NOM=INT((1100._EB-OMMIN)/5._EB)+INT((OMMAX-1100._EB)/25._EB)ELSEIF(OMMAX>5000._EB) THEN   NOM=INT((5000._EB-OMMIN)/25._EB)+INT((OMMAX-5000._EB)/50._EB)ELSE   NOM=INT((OMMAX-OMMIN)/25._EB)ENDIFNPRINT=1END SUBROUTINE INIT_RADCAL!**************************************************************************SUBROUTINE RADCAL(AMEAN,AP0)!INTEGER I,II,J,K,L,KK,NM,N,MM,KMAX,LMAX,KMINREAL (EB) DOM,ABGAS,PTOT,TEMP,UK,GKD,GKDD,XD,YD,XX,ENN,ARG,ARGNEW, &           RSL,RSS,ABLONG,ABSHRT,ABIL,ABIS,OMEGA,WL,DAMBDA,AP0, &           SDWEAK,GDINV,GDDINV,YC,Y,AIWALL,AMEAN,XC,AOM,Q, &           LTERM!  [NOTE: THE TOTAL INTENSITY CALCULATED IS THAT WHICH LEAVES INTERVAL J=1.!  P(I,J) IS PARTIAL PRESSURE, ATM, OF SPECIES I IN  INTERVAL J.!  I=1,2,3,4,5, OR 6 IMPLIES SPECIES IS CO2, H2O, CH4, CO, O2, OR N2, RESP.]DOM=5.0_EBOMEGA=OMMIN-DOMNM=NOM-1!     LOOP 1000 COMPUTES EACH SPECTRAL CONTRIBUTION!     *********************************************L1000: DO KK=1,NOM   OMEGA=OMEGA+DOM   IF(OMEGA>1100._EB) OMEGA=OMEGA+20._EB   IF(OMEGA>5000._EB) OMEGA=OMEGA+25._EB   AMBDA(KK)=10000._EB/OMEGA   ABGAS=0._EB!!          LOOP 200 COMPUTES THE CONTRIBUTION OF EACH SPECIES TO TAU!	   *********************************************************   L200: DO I=1,4!     IF SPECIE(I) IS SET TO 0., THAT PARTICULAR RADIATING SPECIES IS!     NOT PRESENT.  THE SPECIES CONSIDERED ARE!          I   SPECIES!          1     CO2!          2     H2O!          3     CH4!          4     CO!          5     PARTICULATES      IF(SPECIE(I)==0._EB) CYCLE L200!!               LOOP 100 IS FOR EACH ELEMENT ALONG PATH!               ***************************************      L100: DO J=1,NPT!     (CALCULATION PROCEEDS IN ACCORDANCE WITH THE SLG MODEL, TABLE 5-18!     IN NASA SP-3080._EB)         IF(KK<=1) THEN            UUU(I,J)=273./RCT(J)*P(I,J)*100.*DD(J)            GC(I,J)=0._EB            PTOT=0._EB            DO II=1,6               PTOT=P(II,J)+PTOT               GC(I,J)=GC(I,J)+GAMMA(I,II)*P(II,J)*(273./RCT(J))**.5            ENDDO            GC(I,J)=GC(I,J)+GAMMA(I,7)*P(I,J)*273./RCT(J)         ENDIF         IF(P(I,J)==0._EB) THEN            IF(J>1) THEN               XSTAR(J)=XSTAR(J-1)               AC(J)=AC(J-1)               AD(J)=AD(J-1)            ELSE               XSTAR(1)=1.E-34_EB               AC(1)=1._EB               AD(1)=1._EB            ENDIF            X(I,J)=XSTAR(J)         ELSE            TEMP=RCT(J)            SELECT CASE (I)            CASE (1)            CALL CO2(OMEGA,TEMP,GC(1,J),SDWEAK,GDINV,GDDINV)            CASE (2)            CALL H2O(OMEGA,TEMP,GC(2,J),SDWEAK,GDINV,GDDINV)            CASE (3)            CALL FUEL(OMEGA,TEMP,P(3,J),PTOT,GC(3,J),SDWEAK,GDINV,GDDINV)            CASE (4)            CALL CO(OMEGA,TEMP,GC(4,J),SDWEAK,GDINV,GDDINV)               END SELECT            UK=SDWEAK*UUU(I,J)            IF(J==1) THEN               XSTAR(1)=UK+1.E-34_EB               ABGAS=UK/DD(1)+ABGAS               AD(1)=GDDINV               AC(1)=GDINV            ELSE               GKD=UK*GDINV               GKDD=UK*GDDINV               XSTAR(J)=XSTAR(J-1)+UK               AD(J)=(XSTAR(J-1)*AD(J-1)+GKDD)/XSTAR(J)               AC(J)=(XSTAR(J-1)*AC(J-1)+GKD)/XSTAR(J)            ENDIF            IF(XSTAR(J)>=1.E-6_EB) THEN               XD=1.7_EB*AD(J)*(DLOG(1._EB+(XSTAR(J)/1.7_EB/AD(J))**2))**.5_EB               YD=1._EB-(XD/XSTAR(J))**2               XC=XSTAR(J)/(1._EB+XSTAR(J)/4._EB/AC(J))**.5_EB!!   THE FOLLOWING LOOP COMPUTES THE OPTICAL THICKNESS, XC, FOR METHANE USING!   THE GODSON EQUATION AND AN APPROXIMATION TO THE LADENBERG-REICHE!   FUNCTION AS RECOMMENDED BY BROSMER AND TIEN (JQSRT 33,P 521).  THE!   ERROR FUNCTION IS FOUND FROM ITS SERIES EXPANSION.!               IF((I==3._EB).AND.(XC<=10)) THEN                  AOM=XC                  XX=.5_EB*SQRTPI*XC                  IF(XX<=3._EB) THEN                     ENN=1._EB                     DO N=1,30                        ENN=ENN*REAL(N,EB)                        MM=2*N+1                        ARG=1.128379_EB*(-1._EB)**N*((.88622693_EB*XC)**MM)/(REAL(MM,EB)*ENN)                        ARGNEW=ARG+AOM!     IF(ABS(ARG/ARGNEW)<.000001)N=30                        AOM=ARGNEW                     ENDDO                  ELSE                     AOM=1._EB-EXP(-XX**2)/(SQRTPI*XX)                  ENDIF                  IF(AOM>=1._EB)AOM=.9999999_EB                  XC=-DLOG(1._EB-AOM)               ENDIF                    YC=1._EB-(XC/XSTAR(J))**2               Y=MAX(1._EB/YC**2+1._EB/YD**2-1._EB,1._EB)               X(I,J)=XSTAR(J)*((1._EB-(Y**(-.5_EB)))**.5_EB)            ELSE               X(I,J)=XSTAR(J)            ENDIF         ENDIF      ENDDO L100   ENDDO L200!  DETERMINE OPTICAL DEPTH OF SOOT   IF(SPECIE(5)==0._EB) THEN      DO J=1,NPT         XPART(J)=0._EB      ENDDO   ELSE      CALL POD(OMEGA)   ENDIFAB(KK)=ABGAS+XPART(1)/DD(1)     !     EVALUATE THE COMBINED SPECTRAL TRANSMITTANCE AND RADIANCE!     *********************************************************   L500: DO J=1,NPT      XTOT(J)=0._EB      DO I=1,4         IF(SPECIE(I)==0._EB) X(I,J)=0._EB         XTOT(J)=X(I,J)+XTOT(J)      ENDDO      XTOT(J)=XTOT(J)+XPART(J)      IF(XTOT(J)>=99._EB) THEN         TAU(J)=0._EB      ELSE         TAU(J)=EXP(-XTOT(J))      ENDIF      IF(J==1) THEN         QW(KK)=-(TAU(1)-1._EB)*PLANCK(RCT(1),AMBDA(KK))      ELSE         QW(KK)=QW(KK)-(TAU(J)-TAU(J-1))*PLANCK(RCT(J),AMBDA(KK))      ENDIF   ENDDO L500   XT(KK)=XTOT(NPT)   TTAU(KK)=TAU(NPT)   QW(KK)=QW(KK)+TTAU(KK)*PLANCK(TWALL,AMBDA(KK))ENDDO L1000     !     INTEGRATE THE RADIANCE OVER THE SPECTRUM     Q=QW(1)*(AMBDA(1)-AMBDA(2))DO KK=2,NM   Q=Q+QW(KK)*(AMBDA(KK-1)-AMBDA(KK+1))/2._EBENDDOQ=Q+QW(NOM)*(AMBDA(NOM-1)-AMBDA(NOM))     !     DETERMINE SOOT RADIANCE FOR SHORT AND LONG WAVELENGTHS.     RSL=0._EBRSS=0._EBABLONG=0._EBABSHRT=0._EBABIL=0._EBABIS=0._EB   IF(.NOT.(SPECIE(5)==0._EB .AND. TWALL==0._EB)) THEN   KMAX=OMMIN/5*5   DO KK=5,KMAX,5      OMEGA=FLOAT(KK)      WL=10000._EB/OMEGA      DAMBDA=10000._EB/(OMEGA-2.5_EB)-10000._EB/(OMEGA+2.5_EB)      CALL POD(OMEGA)         DO J=1,NPT         IF(XPART(J)>=33._EB) THEN            TAUS(J)=0._EB         ELSE            TAUS(J)=EXP(-XPART(J))         ENDIF         IF(J==1) THEN            RSL=RSL-(TAUS(1)-1._EB)*PLANCK(RCT(1),WL)*DAMBDA            ABLONG=ABLONG+XPART(1)/DD(1)*PLANCK(RCT(1),WL)*DAMBDA*5.5411E7_EB/(RCT(1))**4            ABIL=ABIL+XPART(1)/DD(1)*PLANCK(TWALL,WL)*DAMBDA*5.5411E7_EB /(TWALL+.000001_EB)**4         ELSE            RSL=RSL-(TAUS(J)-TAUS(J-1))*PLANCK(RCT(J),WL)*DAMBDA         ENDIF      ENDDO      RSL=RSL+TAUS(NPT)*PLANCK(TWALL,WL)*DAMBDA   ENDDO   KMIN=OMMAX/100*100   DO KK=KMIN,25000,100      OMEGA=FLOAT(KK)      WL=10000._EB/OMEGA      DAMBDA=10000._EB/(OMEGA-50._EB)-10000._EB/(OMEGA+50._EB)      CALL POD(OMEGA)      DO J=1,NPT         IF(XPART(J)>=33._EB) THEN            TAUS(J)=0._EB         ELSE            TAUS(J)=EXP(-XPART(J))         ENDIF         IF(J/=1) THEN            RSS=RSS-(TAUS(J)-TAUS(J-1))*PLANCK(RCT(J),WL)*DAMBDA         ELSE            RSS=RSS-(TAUS(1)-1._EB)*PLANCK(RCT(1),WL)*DAMBDA            ABSHRT=ABSHRT+XPART(1)/DD(1)*PLANCK(RCT(1),WL)*DAMBDA*5.5411E7_EB/(RCT(1))**4            ABIS=ABIS+XPART(1)/DD(1)*PLANCK(TWALL,WL)*DAMBDA*5.5411E7_EB/(TWALL+.000001_EB)**4         ENDIF      ENDDO      RSS=RSS+TAUS(NPT)*PLANCK(TWALL,WL)*DAMBDA   ENDDOENDIFQ=Q+RSS+RSLIF (NPRINT/=0) THEN   IF(NPRINT==1) THEN      NPRINT=2      DO J=1,NPT         DO I=1,6            PKPA(I)=P(I,J)*101._EB         ENDDO      ENDDO