| 1 | !*************************************************************************** |
|---|
| 2 | ! * |
|---|
| 3 | ! RRTM : RAPID RADIATIVE TRANSFER MODEL * |
|---|
| 4 | ! * |
|---|
| 5 | ! ATMOSPHERIC AND ENVIRONMENTAL RESEARCH, INC. * |
|---|
| 6 | ! 840 MEMORIAL DRIVE * |
|---|
| 7 | ! CAMBRIDGE, MA 02139 * |
|---|
| 8 | ! * |
|---|
| 9 | ! ELI J. MLAWER * |
|---|
| 10 | ! STEVEN J. TAUBMAN~ * |
|---|
| 11 | ! SHEPARD A. CLOUGH * |
|---|
| 12 | ! * |
|---|
| 13 | ! ~currently at GFDL * |
|---|
| 14 | ! * |
|---|
| 15 | ! email: mlawer@aer.com * |
|---|
| 16 | ! * |
|---|
| 17 | ! The authors wish to acknowledge the contributions of the * |
|---|
| 18 | ! following people: Patrick D. Brown, Michael J. Iacono, * |
|---|
| 19 | ! Ronald E. Farren, Luke Chen, Robert Bergstrom. * |
|---|
| 20 | ! * |
|---|
| 21 | !*************************************************************************** |
|---|
| 22 | ! Reformatted for F90 by JJMorcrette, ECMWF, 980714 * |
|---|
| 23 | ! * |
|---|
| 24 | !*************************************************************************** |
|---|
| 25 | ! *** mji *** |
|---|
| 26 | ! *** This version of RRTM has been altered to interface with either |
|---|
| 27 | ! the ECMWF numerical weather prediction model or the ECMWF column |
|---|
| 28 | ! radiation model (ECRT) package. |
|---|
| 29 | |
|---|
| 30 | ! Revised, April, 1997; Michael J. Iacono, AER, Inc. |
|---|
| 31 | ! - initial implementation of RRTM in ECRT code |
|---|
| 32 | ! Revised, June, 1999; Michael J. Iacono and Eli J. Mlawer, AER, Inc. |
|---|
| 33 | ! - to implement generalized maximum/random cloud overlap |
|---|
| 34 | |
|---|
| 35 | SUBROUTINE RRTM_RRTM_140GP & |
|---|
| 36 | &( KIDIA , KFDIA , KLON , KLEV & |
|---|
| 37 | &, PAER , PAPH , PAP & |
|---|
| 38 | &, PTS , PTH , PT & |
|---|
| 39 | &, ZEMIS , ZEMIW & |
|---|
| 40 | &, PQ , PCCO2 , POZN & |
|---|
| 41 | &, PCLDF , PTAUCLD & |
|---|
| 42 | &, PEMIT , PFLUX , PFLUC, PTCLEAR & |
|---|
| 43 | &) |
|---|
| 44 | |
|---|
| 45 | ! *** This program is the driver for RRTM, the AER rapid model. |
|---|
| 46 | ! For each atmosphere the user wishes to analyze, this routine |
|---|
| 47 | ! a) calls ECRTATM to read in the atmospheric profile |
|---|
| 48 | ! b) calls SETCOEF to calculate various quantities needed for |
|---|
| 49 | ! the radiative transfer algorithm |
|---|
| 50 | ! c) calls RTRN to do the radiative transfer calculation for |
|---|
| 51 | ! clear or cloudy sky |
|---|
| 52 | ! d) writes out the upward, downward, and net flux for each |
|---|
| 53 | ! level and the heating rate for each layer |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | #include "tsmbkind.h" |
|---|
| 57 | |
|---|
| 58 | USE PARRRTM , ONLY : JPBAND ,JPG ,JPXSEC ,JPGPT ,JPLAY ,& |
|---|
| 59 | &JPINPX |
|---|
| 60 | USE YOERRTWN , ONLY : WAVENUM1 ,WAVENUM2 ,DELWAVE ,NG ,NSPA ,NSPB |
|---|
| 61 | |
|---|
| 62 | !------------------------------Arguments-------------------------------- |
|---|
| 63 | |
|---|
| 64 | ! Input arguments |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | IMPLICIT NONE |
|---|
| 68 | INTEGER_M :: kidia ! First atmosphere index |
|---|
| 69 | INTEGER_M :: kfdia ! Last atmosphere index |
|---|
| 70 | INTEGER_M :: klon ! Number of atmospheres (longitudes) |
|---|
| 71 | INTEGER_M :: klev ! Number of atmospheric layers |
|---|
| 72 | REAL_B :: paer(klon,6,klev) ! Aerosol optical thickness |
|---|
| 73 | REAL_B :: pap(klon,klev) ! Layer pressures (Pa) |
|---|
| 74 | REAL_B :: paph(klon,klev+1) ! Interface pressures (Pa) |
|---|
| 75 | REAL_B :: pts(klon) ! Surface temperature (K) |
|---|
| 76 | REAL_B :: pt(klon,klev) ! Layer temperature (K) |
|---|
| 77 | REAL_B :: zemis(klon) ! Non-window surface emissivity |
|---|
| 78 | REAL_B :: zemiw(klon) ! Window surface emissivity |
|---|
| 79 | REAL_B :: pth(klon,klev+1) ! Interface temperatures (K) |
|---|
| 80 | REAL_B :: pq(klon,klev) ! H2O specific humidity (mmr) |
|---|
| 81 | REAL_B :: pozn(klon,klev) ! O3 mass mixing ratio |
|---|
| 82 | REAL_B :: pcco2 ! CO2 mass mixing ratio |
|---|
| 83 | REAL_B :: rch4 ! CH4 mass mixing ratio |
|---|
| 84 | REAL_B :: rn2o ! N2O mass mixing ratio |
|---|
| 85 | REAL_B :: rcfc11 ! CFC11 mass mixing ratio |
|---|
| 86 | REAL_B :: rcfc12 ! CFC12 mass mixing ratio |
|---|
| 87 | REAL_B :: pcldf(klon,klev) ! Cloud fraction |
|---|
| 88 | REAL_B :: ptaucld(klon,klev,JPBAND)! Cloud optical depth |
|---|
| 89 | REAL_B :: PFLUX(klon,2,klev+1) ! LW total sky flux (1=up, 2=down) |
|---|
| 90 | REAL_B :: PFLUC(klon,2,klev+1) ! LW clear sky flux (1=up, 2=down) |
|---|
| 91 | REAL_B :: PEMIT(klon) ! Surface LW emissivity |
|---|
| 92 | REAL_B :: PTCLEAR(klon) ! clear-sky fraction of column |
|---|
| 93 | |
|---|
| 94 | INTEGER_M :: ICLDLYR(JPLAY) ! Cloud indicator |
|---|
| 95 | REAL_B :: CLDFRAC(JPLAY) ! Cloud fraction |
|---|
| 96 | REAL_B :: TAUCLD(JPLAY,JPBAND) ! Spectral optical thickness |
|---|
| 97 | |
|---|
| 98 | REAL_B :: ABSS1 (JPGPT*JPLAY) |
|---|
| 99 | REAL_B :: ATR1 (JPGPT,JPLAY) |
|---|
| 100 | EQUIVALENCE (ABSS1(1),ATR1(1,1)) |
|---|
| 101 | |
|---|
| 102 | REAL_B :: OD (JPGPT,JPLAY) |
|---|
| 103 | |
|---|
| 104 | REAL_B :: TAUSF1(JPGPT*JPLAY) |
|---|
| 105 | REAL_B :: TF1 (JPGPT,JPLAY) |
|---|
| 106 | EQUIVALENCE (TAUSF1(1),TF1(1,1)) |
|---|
| 107 | |
|---|
| 108 | REAL_B :: COLDRY(JPLAY) |
|---|
| 109 | REAL_B :: WKL(JPINPX,JPLAY) |
|---|
| 110 | |
|---|
| 111 | REAL_B :: WX(JPXSEC,JPLAY) ! Amount of trace gases |
|---|
| 112 | |
|---|
| 113 | REAL_B :: CLFNET (0:JPLAY) |
|---|
| 114 | REAL_B :: CLHTR (0:JPLAY) |
|---|
| 115 | REAL_B :: FNET (0:JPLAY) |
|---|
| 116 | REAL_B :: HTR (0:JPLAY) |
|---|
| 117 | REAL_B :: TOTDFLUC(0:JPLAY) |
|---|
| 118 | REAL_B :: TOTDFLUX(0:JPLAY) |
|---|
| 119 | REAL_B :: TOTUFLUC(0:JPLAY) |
|---|
| 120 | REAL_B :: TOTUFLUX(0:JPLAY) |
|---|
| 121 | |
|---|
| 122 | ! LOCAL INTEGER SCALARS |
|---|
| 123 | INTEGER_M :: i, icld, iplon, K |
|---|
| 124 | INTEGER_M :: ISTART |
|---|
| 125 | INTEGER_M :: IEND |
|---|
| 126 | |
|---|
| 127 | ! LOCAL REAL SCALARS |
|---|
| 128 | REAL_B :: FLUXFAC, HEATFAC, PI, ZEPSEC, ZTCLEAR |
|---|
| 129 | |
|---|
| 130 | !- from AER |
|---|
| 131 | REAL_B :: TAUAERL(JPLAY,JPBAND) |
|---|
| 132 | |
|---|
| 133 | !- from INTFAC |
|---|
| 134 | REAL_B :: FAC00(JPLAY) |
|---|
| 135 | REAL_B :: FAC01(JPLAY) |
|---|
| 136 | REAL_B :: FAC10(JPLAY) |
|---|
| 137 | REAL_B :: FAC11(JPLAY) |
|---|
| 138 | REAL_B :: FORFAC(JPLAY) |
|---|
| 139 | |
|---|
| 140 | !- from INTIND |
|---|
| 141 | INTEGER_M :: JP(JPLAY) |
|---|
| 142 | INTEGER_M :: JT(JPLAY) |
|---|
| 143 | INTEGER_M :: JT1(JPLAY) |
|---|
| 144 | |
|---|
| 145 | !- from PRECISE |
|---|
| 146 | REAL_B :: ONEMINUS |
|---|
| 147 | |
|---|
| 148 | !- from PROFDATA |
|---|
| 149 | REAL_B :: COLH2O(JPLAY) |
|---|
| 150 | REAL_B :: COLCO2(JPLAY) |
|---|
| 151 | REAL_B :: COLO3 (JPLAY) |
|---|
| 152 | REAL_B :: COLN2O(JPLAY) |
|---|
| 153 | REAL_B :: COLCH4(JPLAY) |
|---|
| 154 | REAL_B :: COLO2 (JPLAY) |
|---|
| 155 | REAL_B :: CO2MULT(JPLAY) |
|---|
| 156 | INTEGER_M :: LAYTROP |
|---|
| 157 | INTEGER_M :: LAYSWTCH |
|---|
| 158 | INTEGER_M :: LAYLOW |
|---|
| 159 | |
|---|
| 160 | !- from PROFILE |
|---|
| 161 | REAL_B :: PAVEL(JPLAY) |
|---|
| 162 | REAL_B :: TAVEL(JPLAY) |
|---|
| 163 | REAL_B :: PZ(0:JPLAY) |
|---|
| 164 | REAL_B :: TZ(0:JPLAY) |
|---|
| 165 | REAL_B :: TBOUND |
|---|
| 166 | INTEGER_M :: NLAYERS |
|---|
| 167 | |
|---|
| 168 | !- from SELF |
|---|
| 169 | REAL_B :: SELFFAC(JPLAY) |
|---|
| 170 | REAL_B :: SELFFRAC(JPLAY) |
|---|
| 171 | INTEGER_M :: INDSELF(JPLAY) |
|---|
| 172 | |
|---|
| 173 | !- from SP |
|---|
| 174 | REAL_B :: PFRAC(JPGPT,JPLAY) |
|---|
| 175 | |
|---|
| 176 | !- from SURFACE |
|---|
| 177 | REAL_B :: SEMISS(JPBAND) |
|---|
| 178 | REAL_B :: SEMISLW |
|---|
| 179 | INTEGER_M :: IREFLECT |
|---|
| 180 | |
|---|
| 181 | |
|---|
| 182 | ! HEATFAC is the factor by which one must multiply delta-flux/ |
|---|
| 183 | ! delta-pressure, with flux in w/m-2 and pressure in mbar, to get |
|---|
| 184 | ! the heating rate in units of degrees/day. It is equal to |
|---|
| 185 | ! (g)x(#sec/day)x(1e-5)/(specific heat of air at const. p) |
|---|
| 186 | ! = (9.8066)(86400)(1e-5)/(1.004) |
|---|
| 187 | |
|---|
| 188 | ZEPSEC = 1.E-06_JPRB |
|---|
| 189 | ONEMINUS = _ONE_ - ZEPSEC |
|---|
| 190 | PI = _TWO_*ASIN(_ONE_) |
|---|
| 191 | FLUXFAC = PI * 2.D4 |
|---|
| 192 | HEATFAC = 8.4391_JPRB |
|---|
| 193 | |
|---|
| 194 | ! *** mji *** |
|---|
| 195 | ! For use with ECRT, this loop is over atmospheres (or longitudes) |
|---|
| 196 | DO iplon = kidia,kfdia |
|---|
| 197 | |
|---|
| 198 | ! *** mji *** |
|---|
| 199 | !- Prepare atmospheric profile from ECRT for use in RRTM, and define |
|---|
| 200 | ! other RRTM input parameters. Arrays are passed back through the |
|---|
| 201 | ! existing RRTM commons and arrays. |
|---|
| 202 | ZTCLEAR=_ONE_ |
|---|
| 203 | |
|---|
| 204 | ! print *,'before RRTM_ECRT_140GP' |
|---|
| 205 | |
|---|
| 206 | CALL RRTM_ECRT_140GP & |
|---|
| 207 | &( iplon, klon , klev, icld & |
|---|
| 208 | &, paer , paph , pap & |
|---|
| 209 | &, pts , pth , pt & |
|---|
| 210 | &, zemis, zemiw & |
|---|
| 211 | &, pq , pcco2, pozn, pcldf, ptaucld, ztclear & |
|---|
| 212 | &, CLDFRAC,TAUCLD,COLDRY,WKL,WX & |
|---|
| 213 | &, TAUAERL,PAVEL,TAVEL,PZ,TZ,TBOUND,NLAYERS,SEMISS,IREFLECT) |
|---|
| 214 | |
|---|
| 215 | PTCLEAR(iplon)=ztclear |
|---|
| 216 | |
|---|
| 217 | ISTART = 1 |
|---|
| 218 | IEND = 16 |
|---|
| 219 | |
|---|
| 220 | ! Calculate information needed by the radiative transfer routine |
|---|
| 221 | ! that is specific to this atmosphere, especially some of the |
|---|
| 222 | ! coefficients and indices needed to compute the optical depths |
|---|
| 223 | ! by interpolating data from stored reference atmospheres. |
|---|
| 224 | |
|---|
| 225 | ! print *,'before RRTM_SETCOEF_140GP' |
|---|
| 226 | |
|---|
| 227 | CALL RRTM_SETCOEF_140GP (KLEV,COLDRY,WKL & |
|---|
| 228 | &, FAC00,FAC01,FAC10,FAC11,FORFAC,JP,JT,JT1 & |
|---|
| 229 | &, COLH2O,COLCO2,COLO3,COLN2O,COLCH4,COLO2,CO2MULT & |
|---|
| 230 | &, LAYTROP,LAYSWTCH,LAYLOW,PAVEL,TAVEL,SELFFAC,SELFFRAC,INDSELF) |
|---|
| 231 | |
|---|
| 232 | ! print *,'before RRTM_GASABS1A_140GP' |
|---|
| 233 | |
|---|
| 234 | CALL RRTM_GASABS1A_140GP (KLEV,ATR1,OD,TF1,COLDRY,WX & |
|---|
| 235 | &, TAUAERL,FAC00,FAC01,FAC10,FAC11,FORFAC,JP,JT,JT1,ONEMINUS & |
|---|
| 236 | &, COLH2O,COLCO2,COLO3,COLN2O,COLCH4,COLO2,CO2MULT & |
|---|
| 237 | &, LAYTROP,LAYSWTCH,LAYLOW,SELFFAC,SELFFRAC,INDSELF,PFRAC) |
|---|
| 238 | |
|---|
| 239 | !- Call the radiative transfer routine. |
|---|
| 240 | |
|---|
| 241 | ! *** mji *** |
|---|
| 242 | ! Check for cloud in column. Use ECRT threshold set as flag icld in |
|---|
| 243 | ! routine ECRTATM. If icld=1 then column is cloudy, otherwise it is |
|---|
| 244 | ! clear. Also, set up flag array, icldlyr, for use in radiative |
|---|
| 245 | ! transfer. Set icldlyr to one for each layer with non-zero cloud |
|---|
| 246 | ! fraction. |
|---|
| 247 | |
|---|
| 248 | DO K = 1, KLEV |
|---|
| 249 | IF (ICLD == 1.AND.CLDFRAC(K) > ZEPSEC) THEN |
|---|
| 250 | ICLDLYR(K) = 1 |
|---|
| 251 | ELSE |
|---|
| 252 | ICLDLYR(K) = 0 |
|---|
| 253 | ENDIF |
|---|
| 254 | ENDDO |
|---|
| 255 | |
|---|
| 256 | ! Clear and cloudy parts of column are treated together in RTRN. |
|---|
| 257 | ! Clear radiative transfer is done for clear layers and cloudy radiative |
|---|
| 258 | ! transfer is done for cloudy layers as identified by icldlyr. |
|---|
| 259 | |
|---|
| 260 | ! print *,'before RRTM_RTRN1A_140GP' |
|---|
| 261 | |
|---|
| 262 | CALL RRTM_RTRN1A_140GP (KLEV,ISTART,IEND,ICLDLYR,CLDFRAC,TAUCLD,ABSS1 & |
|---|
| 263 | &, OD,TAUSF1,CLFNET,CLHTR,FNET,HTR,TOTDFLUC,TOTDFLUX,TOTUFLUC,TOTUFLUX & |
|---|
| 264 | &, TAVEL,PZ,TZ,TBOUND,PFRAC,SEMISS,SEMISLW,IREFLECT) |
|---|
| 265 | |
|---|
| 266 | ! *** Pass clear sky and total sky up and down flux profiles to ECRT |
|---|
| 267 | ! output arrays (zflux, zfluc). Array indexing from bottom to top |
|---|
| 268 | ! is preserved for ECRT. |
|---|
| 269 | ! Invert down flux arrays for consistency with ECRT sign conventions. |
|---|
| 270 | |
|---|
| 271 | pemit(iplon) = SEMISLW |
|---|
| 272 | DO i = 0, KLEV |
|---|
| 273 | pfluc(iplon,1,i+1) = TOTUFLUC(i)*FLUXFAC |
|---|
| 274 | pfluc(iplon,2,i+1) = -TOTDFLUC(i)*FLUXFAC |
|---|
| 275 | pflux(iplon,1,i+1) = TOTUFLUX(i)*FLUXFAC |
|---|
| 276 | pflux(iplon,2,i+1) = -TOTDFLUX(i)*FLUXFAC |
|---|
| 277 | ENDDO |
|---|
| 278 | ENDDO |
|---|
| 279 | |
|---|
| 280 | RETURN |
|---|
| 281 | END SUBROUTINE RRTM_RRTM_140GP |
|---|