| [3] | 1 | SUBROUTINE GFLUXV(NTODO,WAVEN,DTDEL,TDEL,WDEL,CDEL |
|---|
| 2 | & , F0PI,RSF,BTOP,BSURF,FP,FM,FMIDP,FMIDM,IPRINT) |
|---|
| 3 | c PARAMETER (NL=101) |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | INTEGER IPRINT,I,J,K,IDELTA,NAYER,NL,NTODO |
|---|
| 6 | PARAMETER (NL=401) |
|---|
| 7 | C ??FLAG? THIS VALUE (101) MUST BE .GE. NTODO |
|---|
| 8 | C* THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS |
|---|
| 9 | C FOR THE VISIBLE FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT |
|---|
| 10 | C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS |
|---|
| 11 | C MEASURED FROM THE TOP OF EACH LAEYER. (DTAU) TOP OF EACH LAYER HAS |
|---|
| 12 | C OPTICAL DEPTH TAU(N).IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N |
|---|
| 13 | C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES |
|---|
| 14 | C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES. |
|---|
| 15 | C THIS SUBROUTINE DIFFERS FROM ITS IR CONTERPART IN THAT HERE WE SOLVE FOR |
|---|
| 16 | C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR |
|---|
| 17 | C J.A.S., 37, 630-642, 1980. |
|---|
| 18 | REAL B81,B82,R81 |
|---|
| 19 | REAL AP,AM,F0PI,DENOM,UBARI,UBARV,UBAR0,CMMID,TAUMID,RSF,BTOP |
|---|
| 20 | REAL EM,EP,PI,WAVEN,CPMID,BSURF,G4 |
|---|
| 21 | COMMON /UBARED/ UBARI,UBARV,UBAR0 |
|---|
| 22 | C THIS NEXT ROW OF VARIABLES ARE THOSE ACTUALLY USED IN THE |
|---|
| 23 | C ROUTINE |
|---|
| 24 | REAL W0(NL), COSBAR(NL),DTAU(NL),TAU(NL) |
|---|
| 25 | REAL wdel(ntodo-1),cdel(ntodo-1),dtdel(ntodo-1) |
|---|
| 26 | REAL tdel(ntodo),fm(ntodo),fp(ntodo) |
|---|
| 27 | REAL fmidm(ntodo),fmidp(ntodo) |
|---|
| 28 | REAL ALPHA(NL),LAMDA(NL),XK1(NL),XK2(NL) |
|---|
| 29 | REAL G1(NL),G2(NL),G3(NL) |
|---|
| 30 | REAL GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL) |
|---|
| 31 | &,E1(NL),E2(NL),E3(NL),E4(NL),EXPTRM(NL) |
|---|
| 32 | DATA PI/3.14159265358979323846/ |
|---|
| 33 | DATA IDELTA/1/ |
|---|
| 34 | NAYER=NTODO-1 |
|---|
| 35 | 505 CONTINUE |
|---|
| 36 | C TURN ON THE DELTA-FUNCTION IF REQUIRED HERE |
|---|
| 37 | c PRINT*,' UBAR0 ',UBARI,UBARV,UBAR0 |
|---|
| 38 | |
|---|
| 39 | IF (IDELTA .EQ. 0) THEN |
|---|
| 40 | DO 101 J=1,NAYER |
|---|
| 41 | W0(J)=WDEL(J) |
|---|
| 42 | COSBAR(J)=CDEL(J) |
|---|
| 43 | DTAU(J)=DTDEL(J) |
|---|
| 44 | TAU(J)=TDEL(J) |
|---|
| 45 | 101 CONTINUE |
|---|
| 46 | TAU(NTODO)=TDEL(NTODO) |
|---|
| 47 | ELSE |
|---|
| 48 | C FOR THE DELTA FUNCTION HERE... |
|---|
| 49 | TAU(1)=TDEL(1)*(1.-WDEL(1)*CDEL(1)**2) |
|---|
| 50 | DO 102 J=1,NAYER |
|---|
| 51 | W0(J)=WDEL(J)*(1.-CDEL(J)**2)/(1.-WDEL(J)*CDEL(J)**2) |
|---|
| 52 | COSBAR(J)=CDEL(J)/(1.+CDEL(J)) |
|---|
| 53 | DTAU(J)=DTDEL(J)*(1.-WDEL(J)*CDEL(J)**2) |
|---|
| 54 | TAU(J+1)=TAU(J)+DTAU(J) |
|---|
| 55 | c print*,'W0=',W0(J),' COSBAR=',COSBAR(J) |
|---|
| 56 | c print*,'DTAU=',DTAU(J),' TAU=',TAU(J) |
|---|
| 57 | c print*,'WDEL=',WDEL(J),' CDEL=',CDEL(J),' DTDEL=',DTDEL(J) |
|---|
| 58 | 102 CONTINUE |
|---|
| 59 | ENDIF |
|---|
| 60 | C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH |
|---|
| 61 | C AS DEFINED BY M&W - THIS IS THE WAY THE IR IS DONE |
|---|
| 62 | DO 20 J=1,NAYER |
|---|
| 63 | ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) ) |
|---|
| 64 | C THIS SET OF G'S IS FOR THE TWO STREAM HEMI-CONSTANT |
|---|
| 65 | C G1(J)=(1.-W0(J)*0.5*(1.+COSBAR(J)))/UBARV |
|---|
| 66 | C G2(J)=W0(J)*0.5*(1.-COSBAR(J))/UBARV |
|---|
| 67 | C G3(J)=0.5*(1.-COSBAR(J)) |
|---|
| 68 | C SET OF CONSTANTS DETERMINED BY DOM |
|---|
| 69 | G1(J)= (SQRT(3.)*0.5)*(2. - W0(J)*(1.+COSBAR(J))) |
|---|
| 70 | G2(J)= (SQRT(3.)*W0(J)*0.5)*(1.-COSBAR(J)) |
|---|
| 71 | G3(J)=0.5*(1.-SQRT(3.)*COSBAR(J)*UBAR0) |
|---|
| 72 | LAMDA(J)=SQRT(G1(J)**2 - G2(J)**2) |
|---|
| 73 | GAMA(J)=(G1(J)-LAMDA(J))/G2(J) |
|---|
| 74 | |
|---|
| 75 | 20 CONTINUE |
|---|
| 76 | DO 7 J=1,NAYER |
|---|
| 77 | G4=1.-G3(J) |
|---|
| 78 | DENOM=LAMDA(J)**2 - 1./UBAR0**2 |
|---|
| 79 | c print*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 |
|---|
| 80 | c print*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 |
|---|
| 81 | c & - G2(J)**2) |
|---|
| 82 | c print*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 |
|---|
| 83 | C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0 |
|---|
| 84 | C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN |
|---|
| 85 | C THE SCATTERING GOES TO ZERO |
|---|
| 86 | C PREVENT THIS WITH A IF STATEMENT |
|---|
| 87 | IF ( DENOM .EQ. 0.) THEN |
|---|
| 88 | PRINT*,UBAR0 |
|---|
| 89 | DENOM=1.E-4 |
|---|
| 90 | PRINT*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 |
|---|
| 91 | PRINT*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 |
|---|
| 92 | & - G2(J)**2) |
|---|
| 93 | PRINT*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 |
|---|
| 94 | WRITE (6,99) |
|---|
| 95 | 99 FORMAT (' DENOM ZERO; RESET # 1IN GFLUXV, W0=0?') |
|---|
| 96 | END IF |
|---|
| 97 | AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM |
|---|
| 98 | AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM |
|---|
| 99 | C* CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED |
|---|
| 100 | C AT THE TOP OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J) |
|---|
| 101 | CPM1(J)=AP*EXP(-TAU(J)/UBAR0) |
|---|
| 102 | CMM1(J)=AM*EXP(-TAU(J)/UBAR0) |
|---|
| 103 | C* CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
|---|
| 104 | C BOTTOM OF THE LAYER. THAT IS AT HIGHER OPTICAL DEPTH TAU(J+1) |
|---|
| 105 | CP(J)=AP*EXP(-TAU(J+1)/UBAR0) |
|---|
| 106 | CM(J)=AM*EXP(-TAU(J+1)/UBAR0) |
|---|
| 107 | c print*,'AM=',AM,' AP=',AP |
|---|
| 108 | c print*,'CPM1(J)=',CPM1(J),' CMM1(J)=',CMM1(J) |
|---|
| 109 | c print*,'CP(J)=',CP(J),' CM(J)=',CM(J) |
|---|
| 110 | 7 CONTINUE |
|---|
| 111 | C* |
|---|
| 112 | C* NOW CALCULATE THE EXPONENTIAL TERMS NEEDED |
|---|
| 113 | C* FOR THE TRIDIAGONAL ROTATED LAYERED METHOD |
|---|
| 114 | C* WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 |
|---|
| 115 | C* WE CLIPP IT TO AVOID OVERFLOW. |
|---|
| 116 | C* EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS |
|---|
| 117 | C* CORRECTED IN THE DSOLVER ROUTINE. ??FLAG? |
|---|
| 118 | DO 103 J=1,NAYER |
|---|
| 119 | EXPTRM(J)=35. |
|---|
| 120 | IF ( LAMDA(J)*DTAU(J) .LT. 35. ) EXPTRM(J)=LAMDA(J)*DTAU(J) |
|---|
| 121 | 103 CONTINUE |
|---|
| 122 | C* NEED TO CLIPP THE EXPONENTIAL HERE. |
|---|
| 123 | DO 8 J=1,NAYER |
|---|
| 124 | EP=EXP(EXPTRM(J)) |
|---|
| 125 | EM=1.0/EP |
|---|
| 126 | E1(J)=EP+GAMA(J)*EM |
|---|
| 127 | E2(J)=EP-GAMA(J)*EM |
|---|
| 128 | E3(J)=GAMA(J)*EP+EM |
|---|
| 129 | E4(J)=GAMA(J)*EP-EM |
|---|
| 130 | 8 CONTINUE |
|---|
| 131 | C |
|---|
| 132 | B81=BTOP |
|---|
| 133 | B82=BSURF |
|---|
| 134 | R81=RSF |
|---|
| 135 | CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1 |
|---|
| 136 | & ,E1,E2,E3,E4,B81,B82,R81,XK1,XK2) |
|---|
| 137 | C |
|---|
| 138 | C EVALUATE THE NTODO FLUXES THROUGH THE NAYER LAYERS |
|---|
| 139 | C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM |
|---|
| 140 | C AT THE THE TOP OF EACH LAYER,J = LEVEL J |
|---|
| 141 | DO 46 J=1,NAYER |
|---|
| 142 | FP(J)= XK1(J) + GAMA(J)*XK2(J) + CPM1(J) |
|---|
| 143 | FM(J)=GAMA(J)*XK1(J) + XK2(J) + CMM1(J) |
|---|
| 144 | 46 CONTINUE |
|---|
| 145 | C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NTODO |
|---|
| 146 | J=NAYER |
|---|
| 147 | EP=EXP(EXPTRM(J)) |
|---|
| 148 | EM=1.0/EP |
|---|
| 149 | FP(J+1)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CP(J) |
|---|
| 150 | FM(J+1)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CM(J) |
|---|
| 151 | C NOTE THAT WE HAVE SOLVED FOR THE FLUXES DIRECTLY AND NO |
|---|
| 152 | C FURTHER INTEGRATION IS NEEDED, THE UBARV TERM IS ABSORBED |
|---|
| 153 | C INTO THE DEFINITION OF G1 THU G4 |
|---|
| 154 | C UBARV = .5 IS HEMISPHERIC CONSTANT |
|---|
| 155 | C UBARV = SQRT(3) IS GAUSS QUADRATURE |
|---|
| 156 | C AND OTHER CASES AS PER MEADOR AND WEAVOR, JAS, 37, 630-643,1980. |
|---|
| 157 | C |
|---|
| 158 | C ADD THE DIRECT FLUX TERM TO THE DOWNWELLING RADIATION, LIOU 182 |
|---|
| 159 | DO 80 J=1,NTODO |
|---|
| 160 | FM(J)=FM(J)+UBAR0*F0PI*EXP(-TAU(J)/UBAR0) |
|---|
| 161 | 80 CONTINUE |
|---|
| 162 | C |
|---|
| 163 | C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS. |
|---|
| 164 | C EXACLTY ANALOGOUS TO THE ABOVE COMPUTATION |
|---|
| 165 | C |
|---|
| 166 | DO 1982 J=1,NAYER |
|---|
| 167 | EP=EXP(0.5*EXPTRM(J)) |
|---|
| 168 | EM=1.0/EP |
|---|
| 169 | G4=1.-G3(J) |
|---|
| 170 | DENOM=LAMDA(J)**2 - 1./UBAR0**2 |
|---|
| 171 | C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0 |
|---|
| 172 | C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN |
|---|
| 173 | C THE SCATTERING GOES TO ZERO |
|---|
| 174 | C PREVENT THIS WITH A IF STATEMENT |
|---|
| 175 | IF ( DENOM .EQ. 0.) THEN |
|---|
| 176 | DENOM=1.E-4 |
|---|
| 177 | PRINT*,'G1(J)**2- G2(J)**2',G1(J)**2 - G2(J)**2 |
|---|
| 178 | PRINT*,'SQRT(G1(J)**2 - G2(J)**2)',SQRT(G1(J)**2 |
|---|
| 179 | & - G2(J)**2) |
|---|
| 180 | PRINT*,'DENOM=',LAMDA(J)**2,1./UBAR0**2 |
|---|
| 181 | WRITE (6,78) |
|---|
| 182 | 78 FORMAT (' DENOM ZERO; RESET # 2 IN GFLUXV, W0=0?') |
|---|
| 183 | END IF |
|---|
| 184 | AM=F0PI*W0(J)*(G4 *(G1(J)+1./UBAR0) +G2(J)*G3(J) )/DENOM |
|---|
| 185 | AP=F0PI*W0(J)*(G3(J)*(G1(J)-1./UBAR0) +G2(J)*G4 )/DENOM |
|---|
| 186 | C* CPMID AND CMMID ARE THE CPLUS AND CMINUS TERMS EVALUATED |
|---|
| 187 | C AT THE MIDDLE OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)+ |
|---|
| 188 | TAUMID= (TAU(J)+0.5*DTAU(J) ) |
|---|
| 189 | CPMID=AP*EXP(-TAUMID/UBAR0) |
|---|
| 190 | CMMID=AM*EXP(-TAUMID/UBAR0) |
|---|
| 191 | FMIDP(J)=XK1(J)*EP + GAMA(J)*XK2(J)*EM + CPMID |
|---|
| 192 | FMIDM(J)=XK1(J)*EP*GAMA(J) + XK2(J)*EM + CMMID |
|---|
| 193 | C ADD THE DIRECT FLUX TO THE DOWNWELLING TERM |
|---|
| 194 | FMIDM(J)= FMIDM(J) +UBAR0*F0PI*EXP(-TAUMID/UBAR0) |
|---|
| 195 | 1982 CONTINUE |
|---|
| 196 | C ** NOW PRINTOUT IF NECESSARY |
|---|
| 197 | IF (IPRINT .LT. 9) RETURN |
|---|
| 198 | WRITE(6,601) F0PI,RSF,BTOP,BSURF |
|---|
| 199 | DO 120 J=1,NAYER |
|---|
| 200 | WRITE (6,301) TAU(J),FP(J),FM(J),DTAU(J),W0(J),COSBAR(J) |
|---|
| 201 | & ,ALPHA(J), LAMDA(J),G1(J),G2(J),G3(J) |
|---|
| 202 | 120 CONTINUE |
|---|
| 203 | J=NTODO |
|---|
| 204 | WRITE (6,301) TAU(J),FP(J),FM(J) |
|---|
| 205 | WRITE (6,602) |
|---|
| 206 | DO 130 J=1,NAYER |
|---|
| 207 | WRITE (6,301) CP(J),CM(J),E1(J),E2(J),E3(J),E4(J),CPM1(J) |
|---|
| 208 | & ,CMM1(J),GAMA(J) |
|---|
| 209 | 130 CONTINUE |
|---|
| 210 | 301 FORMAT(1X,1P13E10.3) |
|---|
| 211 | 602 FORMAT(' CP(J),CM(J),E1(J),E2(J),EM1(J),EM2(J),CPM1(J),CMM1(J)' |
|---|
| 212 | & ,',GAMA(J)') |
|---|
| 213 | 601 FORMAT(1X,'F0PI,RSF,BTOP,BSURF= ',1P4E10.3,/ |
|---|
| 214 | & ' TAU,FUP,FDOWN,DTAU(J),W0(J),COSBAR(J)' |
|---|
| 215 | & ,',ALPHA(J),LAMDA(J),G1,G2,G3') |
|---|
| 216 | C *******************************************************************012 |
|---|
| 217 | RETURN |
|---|
| 218 | END |
|---|
