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