| 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=201) ! 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 | |
|---|
| 64 | IF (NLL .GT. NLP) STOP 'PARAMETER NL TOO SMALL IN GFLUXV' |
|---|
| 65 | |
|---|
| 66 | NLAYER = L_NLAYRAD |
|---|
| 67 | |
|---|
| 68 | DO L=1,L_NLAYRAD-1 |
|---|
| 69 | ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) ) |
|---|
| 70 | LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI |
|---|
| 71 | |
|---|
| 72 | NT2 = int(TLEV(2*L+2)*10.0D0)-NTstar +1 |
|---|
| 73 | NT = int(TLEV(2*L)*10.0D0)-NTstar + 1 |
|---|
| 74 | ! if temperatures superior to 50K |
|---|
| 75 | ! NT2 = TLEV(2*L+2)*10.0D0-499 |
|---|
| 76 | ! NT = TLEV(2*L)*10.0D0-499 |
|---|
| 77 | |
|---|
| 78 | B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L) |
|---|
| 79 | B0(L) = PLANCKIR(NW,NT) |
|---|
| 80 | END DO |
|---|
| 81 | |
|---|
| 82 | C Take care of special lower layer |
|---|
| 83 | |
|---|
| 84 | L = L_NLAYRAD |
|---|
| 85 | ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) ) |
|---|
| 86 | LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI |
|---|
| 87 | |
|---|
| 88 | |
|---|
| 89 | NT2 = TLEV(2*L+1)*10.0D0-NTstar +1 |
|---|
| 90 | NT = TLEV(2*L)*10.0D0-NTstar + 1 |
|---|
| 91 | ! if temperatures superior to 50K |
|---|
| 92 | ! NT = TLEV(2*L+1)*10.0D0-499 |
|---|
| 93 | ! NT2 = TLEV(2*L)*10.0D0-499 |
|---|
| 94 | B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L) |
|---|
| 95 | B0(L) = PLANCKIR(NW,NT2) |
|---|
| 96 | |
|---|
| 97 | DO L=1,L_NLAYRAD |
|---|
| 98 | GAMA(L) = (1.0-ALPHA(L))/(1.0+ALPHA(L)) |
|---|
| 99 | TERM = UBARI/(1.0-W0(L)*COSBAR(L)) |
|---|
| 100 | |
|---|
| 101 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
|---|
| 102 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
|---|
| 103 | |
|---|
| 104 | CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM |
|---|
| 105 | CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM |
|---|
| 106 | |
|---|
| 107 | C CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED |
|---|
| 108 | C AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH |
|---|
| 109 | |
|---|
| 110 | CPM1(L) = B0(L)+B1(L)*TERM |
|---|
| 111 | CMM1(L) = B0(L)-B1(L)*TERM |
|---|
| 112 | END DO |
|---|
| 113 | |
|---|
| 114 | C NOW CALCULATE THE EXPONENTIAL TERMS NEEDED |
|---|
| 115 | C FOR THE TRIDIAGONAL ROTATED LAYERED METHOD |
|---|
| 116 | C WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX) |
|---|
| 117 | C WE CLIP IT TO AVOID OVERFLOW. |
|---|
| 118 | |
|---|
| 119 | DO L=1,L_NLAYRAD |
|---|
| 120 | |
|---|
| 121 | C CLIP THE EXPONENTIAL HERE. |
|---|
| 122 | |
|---|
| 123 | EP = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX)) |
|---|
| 124 | EM = 1.0/EP |
|---|
| 125 | E1(L) = EP+GAMA(L)*EM |
|---|
| 126 | E2(L) = EP-GAMA(L)*EM |
|---|
| 127 | E3(L) = GAMA(L)*EP+EM |
|---|
| 128 | E4(L) = GAMA(L)*EP-EM |
|---|
| 129 | END DO |
|---|
| 130 | |
|---|
| 131 | c B81=BTOP ! RENAME BEFORE CALLING DSOLVER - used to be to set |
|---|
| 132 | c B82=BSURF ! them to real*8 - but now everything is real*8 |
|---|
| 133 | c R81=RSF ! so this may not be necessary |
|---|
| 134 | |
|---|
| 135 | C Double precision tridiagonal solver |
|---|
| 136 | |
|---|
| 137 | CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP, |
|---|
| 138 | * BSURF,RSF,XK1,XK2) |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | C NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS. |
|---|
| 142 | |
|---|
| 143 | DO L=1,L_NLAYRAD-1 |
|---|
| 144 | DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L) |
|---|
| 145 | EP = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL |
|---|
| 146 | EM = 1.0/EP |
|---|
| 147 | TERM = UBARI/(1.-W0(L)*COSBAR(L)) |
|---|
| 148 | |
|---|
| 149 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
|---|
| 150 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
|---|
| 151 | |
|---|
| 152 | CPMID = B0(L)+B1(L)*DTAUK +B1(L)*TERM |
|---|
| 153 | CMMID = B0(L)+B1(L)*DTAUK -B1(L)*TERM |
|---|
| 154 | FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID |
|---|
| 155 | FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID |
|---|
| 156 | |
|---|
| 157 | C FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT |
|---|
| 158 | |
|---|
| 159 | FMIDP(L) = FMIDP(L)*PI |
|---|
| 160 | FMIDM(L) = FMIDM(L)*PI |
|---|
| 161 | END DO |
|---|
| 162 | |
|---|
| 163 | C And now, for the special bottom layer |
|---|
| 164 | |
|---|
| 165 | L = L_NLAYRAD |
|---|
| 166 | |
|---|
| 167 | EP = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL |
|---|
| 168 | EM = 1.0/EP |
|---|
| 169 | TERM = UBARI/(1.-W0(L)*COSBAR(L)) |
|---|
| 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)*DTAU(L) +B1(L)*TERM |
|---|
| 175 | CMMID = B0(L)+B1(L)*DTAU(L) -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 | |
|---|
| 184 | C FLUX AT THE PTOP LEVEL |
|---|
| 185 | |
|---|
| 186 | EP = 1.0 |
|---|
| 187 | EM = 1.0 |
|---|
| 188 | TERM = UBARI/(1.0-W0(1)*COSBAR(1)) |
|---|
| 189 | |
|---|
| 190 | C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE |
|---|
| 191 | C BOTTOM OF THE LAYER. THAT IS AT DTAU OPTICAL DEPTH |
|---|
| 192 | |
|---|
| 193 | CPMID = B0(1)+B1(1)*TERM |
|---|
| 194 | CMMID = B0(1)-B1(1)*TERM |
|---|
| 195 | |
|---|
| 196 | FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID |
|---|
| 197 | FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID |
|---|
| 198 | |
|---|
| 199 | C FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT |
|---|
| 200 | |
|---|
| 201 | FTOPUP = (FLUXUP-FLUXDN)*PI |
|---|
| 202 | |
|---|
| 203 | |
|---|
| 204 | RETURN |
|---|
| 205 | END |
|---|
