| [3] | 1 | SUBROUTINE COOLING(ngrid,NL,PRESS,TEMP,Z,Q0,lwnet,pfluxi) |
|---|
| 2 | IMPLICIT NONE |
|---|
| 3 | |
|---|
| 4 | c======================================================================= |
|---|
| 5 | c |
|---|
| 6 | c Author : C. P. Mc Kay 01/02/91 |
|---|
| 7 | c ------ |
|---|
| 8 | c |
|---|
| 9 | c Object : |
|---|
| 10 | c -------- |
|---|
| 11 | c |
|---|
| 12 | C THIS SUBROUTINE RETURNS THE COOLING RATE IN TITAN'S ATMOSPHERE |
|---|
| 13 | C INPUTS ARE PRESS(BARS), TEMP(K), Z(KM) |
|---|
| 14 | C OUTPUT IS: Q(K/SEC)C |
|---|
| 15 | C |
|---|
| 16 | C COOLING RATE COMPUTED NEGLECTING SCATTERING. |
|---|
| 17 | C THE TRICK OF THIS ROUTINE IS THAT IT READS IN THE OPACITIES |
|---|
| 18 | C FOR EACH LAYER AT EACH WAVENUMBER IN THE SPECTRAL DOMAIN |
|---|
| 19 | C THESE OPACITIES ARE HELD CONSTANT WITH TEMPERATURE AND TIME. |
|---|
| 20 | c |
|---|
| 21 | c Interface: |
|---|
| 22 | c ---------- |
|---|
| 23 | c |
|---|
| 24 | c Arguments: |
|---|
| 25 | c ---------- |
|---|
| 26 | c |
|---|
| 27 | c input: |
|---|
| 28 | c ------ |
|---|
| 29 | c |
|---|
| 30 | c nl number of levels |
|---|
| 31 | c press(nl) pressure levels (layers) |
|---|
| 32 | c temp(nl) temperature (layers) |
|---|
| 33 | c z(nl) altitude (m, levels) |
|---|
| 34 | c |
|---|
| 35 | c output: |
|---|
| 36 | c ------- |
|---|
| 37 | c |
|---|
| 38 | c q0(nl-1) radiative cooling in K/sec |
|---|
| 39 | c lwnet(nl) net fluxes, (+) upward |
|---|
| 40 | c pfluxi IR descendant a la surface (+ vers le bas) |
|---|
| 41 | c |
|---|
| 42 | c Commons: |
|---|
| 43 | c -------- |
|---|
| 44 | c |
|---|
| 45 | c COMMON/IRTAUS/dtaui(nlayer,nspeci) |
|---|
| 46 | c infrared opacities of the differents layers for differents |
|---|
| 47 | c spectral ranges. This common is initialized by radtitan. |
|---|
| 48 | c |
|---|
| 49 | c======================================================================= |
|---|
| 50 | c----------------------------------------------------------------------- |
|---|
| 51 | c Declarations: |
|---|
| 52 | c ------------ |
|---|
| 53 | |
|---|
| 54 | #include "dimensions.h" |
|---|
| 55 | #include "dimphy.h" |
|---|
| 56 | #include "YOMCST.h" |
|---|
| 57 | INTEGER NLAYER,NSPECI,NSPC1I |
|---|
| 58 | PARAMETER(NLAYER=llm) |
|---|
| 59 | PARAMETER (NSPECI=46,NSPC1I=47) |
|---|
| 60 | |
|---|
| 61 | c Arguments: |
|---|
| 62 | c ---------- |
|---|
| 63 | |
|---|
| 64 | INTEGER NL,ngrid |
|---|
| 65 | REAL PRESS(ngrid,nl),TEMP(ngrid,nl) |
|---|
| 66 | REAL Z(ngrid,nl),Q0(ngrid,nl-1) |
|---|
| 67 | REAL lwnet(ngrid,nl),UBARI2 |
|---|
| 68 | real pfluxi(ngrid) |
|---|
| 69 | |
|---|
| 70 | |
|---|
| 71 | c Common: |
|---|
| 72 | c ------- |
|---|
| 73 | |
|---|
| 74 | C DTAU IS PASSED EN-MASS, SO ITS DEMENSIONS ARE CRITICAL |
|---|
| 75 | COMMON /IRTAUS/ DTAUI(klon,NLAYER,NSPECI) |
|---|
| 76 | REAL dtaui |
|---|
| 77 | |
|---|
| 78 | c Local: |
|---|
| 79 | c ------ |
|---|
| 80 | |
|---|
| 81 | REAL WNOI(NSPECI),DWNI(NSPECI) ! SPECTAL INTERVALS |
|---|
| 82 | REAL B0(klon,llm+1) |
|---|
| 83 | REAL EM(klon,llm+1) |
|---|
| 84 | REAL DW,WAVEN,TJ,BSURF,QOUT,QIN,eff_g,COLDEN |
|---|
| 85 | |
|---|
| 86 | INTEGER ig,K,J,NLEVEL,I,L |
|---|
| 87 | |
|---|
| 88 | c EXTERNAL PLNCK |
|---|
| 89 | REAL PLNCK,zz1,zz2,zz3,zz4,WAVNUM,Xtest |
|---|
| 90 | |
|---|
| 91 | REAL FNETIS(ngrid,llm+1),FNETI(ngrid,llm+1) |
|---|
| 92 | REAL FDIS(ngrid,llm+1,nspeci),FUPIS(ngrid,llm+1,nspeci) |
|---|
| 93 | REAL FDI(ngrid,llm+1), FUPI(ngrid,llm+1) |
|---|
| 94 | |
|---|
| 95 | c Data: |
|---|
| 96 | c ----- |
|---|
| 97 | |
|---|
| 98 | REAL RHOP,CSUBP,UBARI,RSFI |
|---|
| 99 | DATA RHOP/1.E4/ ! CONVERSION FROM PRESSURE TO MASS |
|---|
| 100 | DATA UBARI/0.5/ ! MEAN COSINE FOR 2-STREAM |
|---|
| 101 | DATA RSFI/0.0/ ! SURFACE ALBEDO |
|---|
| 102 | DATA WNOI/ |
|---|
| 103 | & 11.500, 20.000, 31.250, 50.000, 75.000, |
|---|
| 104 | & 100.000, 125.000, 150.000, 175.000, 200.000, |
|---|
| 105 | & 225.000, 250.000, 275.000, 300.000, 325.000, |
|---|
| 106 | & 350.000, 375.000, 400.000, 425.000, 450.000, |
|---|
| 107 | & 475.000, 500.000, 525.000, 550.000, 575.000, |
|---|
| 108 | & 600.000, 628.750, 662.838, 681.757, 683.919, |
|---|
| 109 | & 686.541, 689.623, 692.704, 695.786, 715.141, |
|---|
| 110 | & 733.836, 735.597, 737.358, 739.119, 742.720, |
|---|
| 111 | & 748.160, 753.600, 834.560, 917.333, 926.400, |
|---|
| 112 | & 935.466/ |
|---|
| 113 | DATA DWNI/ |
|---|
| 114 | & 7.000, 10.000, 12.500, 25.000, 25.000, |
|---|
| 115 | & 25.000, 25.000, 25.000, 25.000, 25.000, |
|---|
| 116 | & 25.000, 25.000, 25.000, 25.000, 25.000, |
|---|
| 117 | & 25.000, 25.000, 25.000, 25.000, 25.000, |
|---|
| 118 | & 25.000, 25.000, 25.000, 25.000, 25.000, |
|---|
| 119 | & 25.000, 32.500, 35.676, 2.162, 2.162, |
|---|
| 120 | & 3.082, 3.082, 3.082, 3.082, 35.629, |
|---|
| 121 | & 1.761, 1.761, 1.761, 1.761, 5.440, |
|---|
| 122 | & 5.440, 5.440, 156.480, 9.067, 9.067, |
|---|
| 123 | & 9.067/ |
|---|
| 124 | |
|---|
| 125 | |
|---|
| 126 | save RHOP,UBARI,RSFI,WNOI,DWNI |
|---|
| 127 | |
|---|
| 128 | c----------------------------------------------------------------------- |
|---|
| 129 | |
|---|
| 130 | c Initialisations: |
|---|
| 131 | c ---------------- |
|---|
| 132 | |
|---|
| 133 | CSUBP = RCPD*1e4 ! HEAT CAPACITY OF TITANS ATMOSPHERE in CGS |
|---|
| 134 | |
|---|
| 135 | UBARI2=1./1.66 |
|---|
| 136 | UBARI2=UBARI |
|---|
| 137 | NLEVEL=NL |
|---|
| 138 | |
|---|
| 139 | C ZERO THE NET FLUXES |
|---|
| 140 | Q0 = 0.0 |
|---|
| 141 | lwnet = 0.0 |
|---|
| 142 | |
|---|
| 143 | c----------------------------------------------------------------------- |
|---|
| 144 | C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE INFRARED |
|---|
| 145 | C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL |
|---|
| 146 | c----------------------------------------------------------------------- |
|---|
| 147 | |
|---|
| 148 | DO 2000 K=1,NSPECI ! *** START OF SPECTRAL LOOP |
|---|
| 149 | |
|---|
| 150 | c----------------------------------------------------------------------- |
|---|
| 151 | C SET UP ALTITIDUE PARAMETERS |
|---|
| 152 | |
|---|
| 153 | WAVEN=WNOI(K) |
|---|
| 154 | DW=DWNI(K) |
|---|
| 155 | zz1=DW/(2.*2) |
|---|
| 156 | EM = 0. |
|---|
| 157 | B0 = 0. |
|---|
| 158 | |
|---|
| 159 | DO J=1,NLEVEL-1 |
|---|
| 160 | DO ig=1,ngrid |
|---|
| 161 | TJ=TEMP(ig,J) |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | C Modif: in-lining de la fonction planck pour vectorisation |
|---|
| 165 | C B0(ig,J)=PLNCK(WAVEN,TJ,DW) |
|---|
| 166 | C FUNCTION PLNCK(WAV,T,DW) |
|---|
| 167 | C* PLNCK FUNCTION RETURNS B IN CGS UNITS, ERGS CM-2 WAVENUMBER-1 |
|---|
| 168 | C* WAVNUM IS WAVENUMBER IN CM-1 |
|---|
| 169 | C* T IS IN KELVIN |
|---|
| 170 | PLNCK=0. |
|---|
| 171 | DO I=-2,2,1 |
|---|
| 172 | WAVNUM=WAVEN + I*zz1 |
|---|
| 173 | zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,J)) |
|---|
| 174 | zz3=WAVNUM*WAVNUM*WAVNUM |
|---|
| 175 | PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2) |
|---|
| 176 | ENDDO |
|---|
| 177 | B0(ig,J)=.2*PLNCK |
|---|
| 178 | ENDDO |
|---|
| 179 | |
|---|
| 180 | DO ig=1,ngrid |
|---|
| 181 | zz4=EXP(-DTAUI(ig,J,K)/UBARI2) |
|---|
| 182 | EM(ig,J)=zz4 |
|---|
| 183 | ENDDO |
|---|
| 184 | ENDDO |
|---|
| 185 | |
|---|
| 186 | c----------------------------------------------------------------------- |
|---|
| 187 | C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL |
|---|
| 188 | C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY |
|---|
| 189 | |
|---|
| 190 | FDI =0. |
|---|
| 191 | FDIS =0. |
|---|
| 192 | FUPI =0. |
|---|
| 193 | FUPIS=0. |
|---|
| 194 | |
|---|
| 195 | DO 2220 J=1,NLEVEL-1 |
|---|
| 196 | DO 2230 ig=1,ngrid |
|---|
| 197 | FDI(ig,J+1) = FDI(ig,J)*EM(ig,J) + 2.*RPI*UBARI* |
|---|
| 198 | & B0(ig,J)*(1.-EM(ig,J)) |
|---|
| 199 | FDIS(ig,J+1,K) = FDIS(ig,J,K)*EM(ig,J) + 2.*RPI*UBARI* |
|---|
| 200 | & B0(ig,J)*(1.-EM(ig,J)) |
|---|
| 201 | 2230 CONTINUE |
|---|
| 202 | 2220 CONTINUE |
|---|
| 203 | c write(*,*) |
|---|
| 204 | c write(*,*) 'cooling : EM =' , |
|---|
| 205 | c & ((EM(i,l),l=1,nl),i=1,ngrid) |
|---|
| 206 | c write(*,*) |
|---|
| 207 | c write(*,*) 'cooling : B0 =' , |
|---|
| 208 | c & ((B0(i,l),l=1,nl),i=1,ngrid) |
|---|
| 209 | c write(*,*) |
|---|
| 210 | c write(*,*) 'cooling : FDI =' , |
|---|
| 211 | c & ((FDI(i,l),l=1,nl),i=1,ngrid) |
|---|
| 212 | |
|---|
| 213 | c----------------------------------------------------------------------- |
|---|
| 214 | C UPWARD FLUXES: SURFACE EMISSIONS |
|---|
| 215 | |
|---|
| 216 | DO 2310 ig=1,ngrid |
|---|
| 217 | PLNCK=0. |
|---|
| 218 | DO I=-2,2,1 |
|---|
| 219 | WAVNUM=WAVEN + I*zz1 |
|---|
| 220 | zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NLEVEL)) |
|---|
| 221 | zz3=WAVNUM*WAVNUM*WAVNUM |
|---|
| 222 | PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2) |
|---|
| 223 | ENDDO |
|---|
| 224 | c BSURF=PLNCK( WAVEN, TEMP(ig,NLEVEL), DW) |
|---|
| 225 | BSURF=.2*PLNCK |
|---|
| 226 | FUPI(ig,NLEVEL)=BSURF * 2.*RPI*UBARI + RSFI*FDI(ig,NLEVEL) |
|---|
| 227 | FUPIS(ig,NLEVEL,K)=BSURF*2.*RPI*UBARI+RSFI*FDIS(ig,NLEVEL,K) |
|---|
| 228 | 2310 CONTINUE |
|---|
| 229 | c write(*,*) |
|---|
| 230 | c write(*,*) 'cooling : FUPI/NLEVEL =' , |
|---|
| 231 | c & ((FUPI(i,l),l=nl,nl),i=1,ngrid) |
|---|
| 232 | c write(*,*) |
|---|
| 233 | c write(*,*) 'cooling : FDI/NLEVEL =' , |
|---|
| 234 | c & ((FDI(i,l),l=nl,nl),i=1,ngrid) |
|---|
| 235 | |
|---|
| 236 | DO 2320 J=NLEVEL-1,1,-1 |
|---|
| 237 | DO 2330 ig=1,ngrid |
|---|
| 238 | FUPI(ig,J) = FUPI(ig,J+1)*EM(ig,J) + 2.*RPI*UBARI* |
|---|
| 239 | & B0(ig,J)*(1.-EM(ig,J)) |
|---|
| 240 | FUPIS(ig,J,K) = FUPIS(ig,J+1,K)*EM(ig,J)+2.*RPI*UBARI* |
|---|
| 241 | & B0(ig,J)*(1.-EM(ig,J)) |
|---|
| 242 | 2330 CONTINUE |
|---|
| 243 | 2320 CONTINUE |
|---|
| 244 | c write(*,*) |
|---|
| 245 | c write(*,*) 'cooling : EM =' , |
|---|
| 246 | c & ((EM(i,l),l=1,nl),i=1,ngrid) |
|---|
| 247 | c write(*,*) |
|---|
| 248 | c write(*,*) 'cooling : B0 =' , |
|---|
| 249 | c & ((B0(i,l),l=1,nl),i=1,ngrid) |
|---|
| 250 | c write(*,*) |
|---|
| 251 | c write(*,*) 'cooling : FUPI =' , |
|---|
| 252 | c & ((FUPI(i,l),l=1,nl),i=1,ngrid) |
|---|
| 253 | |
|---|
| 254 | c compute the downward IR flux at the surface: |
|---|
| 255 | c |
|---|
| 256 | DO 3520 ig=1,ngrid |
|---|
| 257 | pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NLEVEL) |
|---|
| 258 | 3520 CONTINUE |
|---|
| 259 | |
|---|
| 260 | c compute the net IR flux, (+) upward: |
|---|
| 261 | c |
|---|
| 262 | DO J=1,NLEVEL |
|---|
| 263 | DO ig=1,ngrid |
|---|
| 264 | lwnet(ig,J)= lwnet(ig,J)+ DWNI(K)*(FUPI(ig,J)-FDI(ig,J)) |
|---|
| 265 | ENDDO |
|---|
| 266 | ENDDO |
|---|
| 267 | |
|---|
| 268 | DO 3210 J=1,NLEVEL-1 |
|---|
| 269 | DO 3220 ig=1,ngrid |
|---|
| 270 | QOUT=FUPI(ig,J) + FDI(ig,J+1) ! OUT OF LAYER |
|---|
| 271 | QIN =FDI(ig,J) + FUPI(ig,J+1) ! INTO LAYER |
|---|
| 272 | Q0(ig,J)=Q0(ig,J)+(QOUT-QIN)*DWNI(K) |
|---|
| 273 | 3220 CONTINUE |
|---|
| 274 | 3210 CONTINUE |
|---|
| 275 | |
|---|
| 276 | c write(*,*) |
|---|
| 277 | c write(*,*) 'cooling/loop : FUPI =' , |
|---|
| 278 | c & ((FUPI(i,l),l=1,nl),i=1,ngrid) |
|---|
| 279 | c write(*,*) |
|---|
| 280 | c write(*,*) 'cooling : FDI =' , |
|---|
| 281 | c & ((FDI(i,l),l=1,nl),i=1,ngrid) |
|---|
| 282 | c write(*,*) |
|---|
| 283 | c write(*,*) 'cooling : Q0 =' , |
|---|
| 284 | c & ((Q0(i,l),l=1,nl-1),i=1,ngrid) |
|---|
| 285 | |
|---|
| 286 | |
|---|
| 287 | c----------------------------------------------------------------------- |
|---|
| 288 | |
|---|
| 289 | 2000 CONTINUE ! *** END SPECTRAL INTERVAL COMPUTATIONS |
|---|
| 290 | |
|---|
| 291 | |
|---|
| 292 | c----------------------------------------------------------------------- |
|---|
| 293 | |
|---|
| 294 | c convertion erg/cm2 -> J/m2 |
|---|
| 295 | DO 3550 ig=1,ngrid |
|---|
| 296 | pfluxi(ig) = 1.e-3*pfluxi(ig) |
|---|
| 297 | lwnet(ig,:) = 1.e-3*lwnet(ig,:) |
|---|
| 298 | 3550 CONTINUE |
|---|
| 299 | |
|---|
| 300 | c PRINT*,'flux IR' |
|---|
| 301 | c WRITE(*,'(8e10.2)') pfluxi |
|---|
| 302 | |
|---|
| 303 | C COMPUTE THE BASELINE COOLING RATE |
|---|
| 304 | |
|---|
| 305 | DO 3000 J=1,NLEVEL-1 |
|---|
| 306 | C TURN THE Q'S INTO TIMESCALES..... |
|---|
| 307 | DO 3300 ig=1,ngrid |
|---|
| 308 | eff_g = RG*(RA/(RA+Z(ig,J)))**2 ! 10% DIFF AT 1 MBAR |
|---|
| 309 | COLDEN = RHOP*(PRESS(ig,J+1)-PRESS(ig,J))/eff_g |
|---|
| 310 | c Q0(J) = (COLDEN * CSUBP )/Q0(J) |
|---|
| 311 | Q0(ig,J) = Q0(ig,J) / (COLDEN*CSUBP) |
|---|
| 312 | 3300 CONTINUE |
|---|
| 313 | 3000 CONTINUE |
|---|
| 314 | |
|---|
| 315 | c write(*,*) |
|---|
| 316 | c write(*,*) 'cooling/end : Q0 =' |
|---|
| 317 | c write(*,*) ((Q0(k,l)*1e7,l=1,nl-1),k=1,ngrid) |
|---|
| 318 | c----------------------------------------------------------------------- |
|---|
| 319 | |
|---|
| 320 | RETURN |
|---|
| 321 | END |
|---|
