[135] | 1 | SUBROUTINE GFLUXI(NLL,TLEV,NW,DW,DTAU,TAUCUM,W0,COSBAR,UBARI, |
---|
| 2 | * RSF,BTOP,BSURF,FTOPUP,FMIDP,FMIDM) |
---|
| 3 | |
---|
| 4 | C THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS |
---|
| 5 | C FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT |
---|
| 6 | C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS |
---|
| 7 | C MEASURED FROM THE TOP OF EACH LAYER. THE TOP OF EACH LAYER HAS |
---|
| 8 | C OPTICAL DEPTH ZERO. IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N |
---|
| 9 | C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES |
---|
| 10 | C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES. |
---|
| 11 | C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY |
---|
| 12 | C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER |
---|
| 13 | C |
---|
| 14 | C NLL = NUMBER OF LEVELS (NLAYERS + 1) MUST BE LESS THAT NL (101) |
---|
| 15 | C TLEV(L_LEVELS) = ARRAY OF TEMPERATURES AT GCM LEVELS |
---|
| 16 | C WAVEN = WAVELENGTH FOR THE COMPUTATION |
---|
| 17 | C DW = WAVENUMBER INTERVAL |
---|
| 18 | C DTAU(NLAYER) = ARRAY OPTICAL DEPTH OF THE LAYERS |
---|
| 19 | C W0(NLEVEL) = SINGLE SCATTERING ALBEDO |
---|
| 20 | C COSBAR(NLEVEL) = ASYMMETRY FACTORS, 0=ISOTROPIC |
---|
| 21 | C UBARI = AVERAGE ANGLE, MUST BE EQUAL TO 0.5 IN IR |
---|
| 22 | C RSF = SURFACE REFLECTANCE |
---|
| 23 | C BTOP = UPPER BOUNDARY CONDITION ON IR INTENSITY (NOT FLUX) |
---|
| 24 | C BSURF = SURFACE EMISSION = (1-RSFI)*PLANCK, INTENSITY (NOT FLUX) |
---|
| 25 | C FP(NLEVEL) = UPWARD FLUX AT LEVELS |
---|
| 26 | C FM(NLEVEL) = DOWNWARD FLUX AT LEVELS |
---|
| 27 | C FMIDP(NLAYER) = UPWARD FLUX AT LAYER MIDPOINTS |
---|
| 28 | C FMIDM(NLAYER) = DOWNWARD FLUX AT LAYER MIDPOINTS |
---|
| 29 | C |
---|
| 30 | C----------------------------------------------------------------------C |
---|
| 31 | |
---|
| 32 | use radinc_h |
---|
| 33 | use radcommon_h, only: planckir |
---|
| 34 | |
---|
| 35 | implicit none |
---|
| 36 | |
---|
| 37 | #include "comcstfi.h" |
---|
| 38 | |
---|
| 39 | INTEGER NLP |
---|
| 40 | PARAMETER (NLP=101) ! MUST BE LARGER THAN NLEVEL |
---|
| 41 | |
---|
| 42 | INTEGER NLL, NLAYER, L, NW, NT, NT2 |
---|
| 43 | REAL*8 TERM, CPMID, CMMID |
---|
| 44 | REAL*8 PLANCK |
---|
| 45 | REAL*8 EM,EP |
---|
| 46 | REAL*8 COSBAR(L_NLAYRAD), W0(L_NLAYRAD), DTAU(L_NLAYRAD) |
---|
| 47 | REAL*8 TAUCUM(L_LEVELS), DTAUK |
---|
| 48 | REAL*8 TLEV(L_LEVELS) |
---|
| 49 | REAL*8 WAVEN, DW, UBARI, RSF |
---|
| 50 | REAL*8 BTOP, BSURF, FMIDP(L_NLAYRAD), FMIDM(L_NLAYRAD) |
---|
| 51 | REAL*8 B0(NLP),B1(NLP),ALPHA(NLP),LAMDA(NLP),XK1(NLP),XK2(NLP) |
---|
| 52 | REAL*8 GAMA(NLP),CP(NLP),CM(NLP),CPM1(NLP),CMM1(NLP),E1(NLP) |
---|
| 53 | REAL*8 E2(NLP),E3(NLP),E4(NLP) |
---|
| 54 | |
---|
| 55 | REAL*8 FTOPUP, FLUXUP, FLUXDN |
---|
| 56 | |
---|
| 57 | real*8 :: TAUMAX = L_TAUMAX |
---|
| 58 | |
---|
| 59 | C======================================================================C |
---|
| 60 | |
---|
| 61 | C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH IN THE INFRARED |
---|
| 62 | |
---|
| 63 | |
---|
[543] | 64 | IF (NLL .GT. NLP) STOP 'PARAMETER NL TOO SMALL IN GFLUXI' |
---|
[135] | 65 | |
---|
| 66 | NLAYER = L_NLAYRAD |
---|
| 67 | |
---|
| 68 | DO L=1,L_NLAYRAD-1 |
---|
[804] | 69 | |
---|
| 70 | !---------------------------------------------------- |
---|
| 71 | ! There is a problem when W0 = 1 |
---|
| 72 | ! open(888,file='W0') |
---|
| 73 | ! if ((W0(L).eq.0.).or.(W0(L).eq.1.)) then |
---|
| 74 | ! write(888,*) W0(L), L, 'gfluxi' |
---|
| 75 | ! endif |
---|
| 76 | ! Prevent this with an if statement: |
---|
[959] | 77 | if (W0(L).eq.1.D0) then |
---|
| 78 | W0(L) = 0.99999D0 |
---|
[804] | 79 | endif |
---|
| 80 | !---------------------------------------------------- |
---|
| 81 | |
---|
[959] | 82 | ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) ) |
---|
| 83 | LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI |
---|
[135] | 84 | |
---|
[543] | 85 | NT = int(TLEV(2*L)*NTfac) - NTstar+1 |
---|
| 86 | NT2 = int(TLEV(2*L+2)*NTfac) - NTstar+1 |
---|
[135] | 87 | |
---|
| 88 | B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L) |
---|
| 89 | B0(L) = PLANCKIR(NW,NT) |
---|
| 90 | END DO |
---|
| 91 | |
---|
| 92 | C Take care of special lower layer |
---|
| 93 | |
---|
| 94 | L = L_NLAYRAD |
---|
[804] | 95 | |
---|
| 96 | if (W0(L).eq.1.) then |
---|
[959] | 97 | W0(L) = 0.99999D0 |
---|
[804] | 98 | end if |
---|
| 99 | |
---|
[959] | 100 | ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) ) |
---|
| 101 | LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI |
---|
[135] | 102 | |
---|
[995] | 103 | ! Tsurf is used for 1st layer source function |
---|
| 104 | ! -- same results for most thin atmospheres |
---|
| 105 | ! -- and stabilizes integrations |
---|
[543] | 106 | NT = int(TLEV(2*L+1)*NTfac) - NTstar+1 |
---|
[995] | 107 | |
---|
| 108 | !! For deep, opaque, thick first layers (e.g. Saturn) |
---|
| 109 | !! what is below works much better, not unstable, ... |
---|
| 110 | !! ... and actually fully accurate because 1st layer temp (JL) |
---|
| 111 | !NT = int(TLEV(2*L)*NTfac) - NTstar+1 |
---|
| 112 | !! (or this one yields same results |
---|
| 113 | !NT = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstar+1 |
---|
| 114 | |
---|
[543] | 115 | NT2 = int(TLEV(2*L)*NTfac) - NTstar+1 |
---|
[135] | 116 | B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L) |
---|
| 117 | B0(L) = PLANCKIR(NW,NT2) |
---|
| 118 | |
---|
| 119 | DO L=1,L_NLAYRAD |
---|
[959] | 120 | GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L)) |
---|
| 121 | TERM = UBARI/(1.0D0-W0(L)*COSBAR(L)) |
---|
[135] | 122 | |
---|
| 123 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
---|
| 124 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
---|
| 125 | |
---|
| 126 | CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM |
---|
| 127 | CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM |
---|
| 128 | |
---|
| 129 | C CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED |
---|
| 130 | C AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH |
---|
| 131 | |
---|
| 132 | CPM1(L) = B0(L)+B1(L)*TERM |
---|
| 133 | CMM1(L) = B0(L)-B1(L)*TERM |
---|
| 134 | END DO |
---|
| 135 | |
---|
| 136 | C NOW CALCULATE THE EXPONENTIAL TERMS NEEDED |
---|
| 137 | C FOR THE TRIDIAGONAL ROTATED LAYERED METHOD |
---|
| 138 | C WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX) |
---|
| 139 | C WE CLIP IT TO AVOID OVERFLOW. |
---|
| 140 | |
---|
| 141 | DO L=1,L_NLAYRAD |
---|
| 142 | |
---|
| 143 | C CLIP THE EXPONENTIAL HERE. |
---|
| 144 | |
---|
| 145 | EP = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX)) |
---|
[959] | 146 | EM = 1.0D0/EP |
---|
[135] | 147 | E1(L) = EP+GAMA(L)*EM |
---|
| 148 | E2(L) = EP-GAMA(L)*EM |
---|
| 149 | E3(L) = GAMA(L)*EP+EM |
---|
| 150 | E4(L) = GAMA(L)*EP-EM |
---|
| 151 | END DO |
---|
| 152 | |
---|
| 153 | c B81=BTOP ! RENAME BEFORE CALLING DSOLVER - used to be to set |
---|
| 154 | c B82=BSURF ! them to real*8 - but now everything is real*8 |
---|
| 155 | c R81=RSF ! so this may not be necessary |
---|
| 156 | |
---|
| 157 | C Double precision tridiagonal solver |
---|
| 158 | |
---|
| 159 | CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP, |
---|
| 160 | * BSURF,RSF,XK1,XK2) |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS. |
---|
| 164 | |
---|
| 165 | DO L=1,L_NLAYRAD-1 |
---|
| 166 | DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L) |
---|
| 167 | EP = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL |
---|
[959] | 168 | EM = 1.0D0/EP |
---|
[961] | 169 | TERM = UBARI/(1.D0-W0(L)*COSBAR(L)) |
---|
[135] | 170 | |
---|
| 171 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
---|
| 172 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
---|
| 173 | |
---|
| 174 | CPMID = B0(L)+B1(L)*DTAUK +B1(L)*TERM |
---|
| 175 | CMMID = B0(L)+B1(L)*DTAUK -B1(L)*TERM |
---|
| 176 | FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID |
---|
| 177 | FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID |
---|
| 178 | |
---|
| 179 | C FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT |
---|
| 180 | |
---|
| 181 | FMIDP(L) = FMIDP(L)*PI |
---|
| 182 | FMIDM(L) = FMIDM(L)*PI |
---|
| 183 | END DO |
---|
| 184 | |
---|
| 185 | C And now, for the special bottom layer |
---|
| 186 | |
---|
| 187 | L = L_NLAYRAD |
---|
| 188 | |
---|
| 189 | EP = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL |
---|
[959] | 190 | EM = 1.0D0/EP |
---|
| 191 | TERM = UBARI/(1.D0-W0(L)*COSBAR(L)) |
---|
[135] | 192 | |
---|
| 193 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
---|
| 194 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
---|
| 195 | |
---|
| 196 | CPMID = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM |
---|
| 197 | CMMID = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM |
---|
| 198 | FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID |
---|
| 199 | FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID |
---|
| 200 | |
---|
| 201 | C FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT |
---|
| 202 | |
---|
| 203 | FMIDP(L) = FMIDP(L)*PI |
---|
| 204 | FMIDM(L) = FMIDM(L)*PI |
---|
| 205 | |
---|
| 206 | C FLUX AT THE PTOP LEVEL |
---|
| 207 | |
---|
[959] | 208 | EP = 1.0D0 |
---|
| 209 | EM = 1.0D0 |
---|
| 210 | TERM = UBARI/(1.0D0-W0(1)*COSBAR(1)) |
---|
[135] | 211 | |
---|
| 212 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
---|
| 213 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
---|
| 214 | |
---|
| 215 | CPMID = B0(1)+B1(1)*TERM |
---|
| 216 | CMMID = B0(1)-B1(1)*TERM |
---|
| 217 | |
---|
| 218 | FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID |
---|
| 219 | FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID |
---|
| 220 | |
---|
| 221 | C FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT |
---|
| 222 | |
---|
| 223 | FTOPUP = (FLUXUP-FLUXDN)*PI |
---|
| 224 | |
---|
| 225 | |
---|
| 226 | RETURN |
---|
| 227 | END |
---|