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