SUBROUTINE COOLING(NG,NL,PRESS,TEMP,Z,Q0,lwnet,pfluxi,icld) c======================================================================= c c Author : C. P. Mc Kay 01/02/91 c ------ c c Object : c -------- c C THIS SUBROUTINE RETURNS THE COOLING RATE IN TITAN'S ATMOSPHERE C INPUTS ARE PRESS(BARS), TEMP(K), Z(KM) C OUTPUT IS: Q(K/SEC)C C C COOLING RATE COMPUTED NEGLECTING SCATTERING. C THE TRICK OF THIS ROUTINE IS THAT IT READS IN THE OPACITIES C FOR EACH LAYER AT EACH WAVENUMBER IN THE SPECTRAL DOMAIN C THESE OPACITIES ARE HELD CONSTANT WITH TEMPERATURE AND TIME. c c Interface: c ---------- c c Arguments: c ---------- c c input: c ------ c c nl number of levels c press(nl) pressure levels (layers) c temp(nl) temperature (layers) c z(nl) altitude (m, levels) c c output: c ------- c c q0(nl-1) radiative cooling in K/sec c lwnet(nl) net fluxes, (+) upward c pfluxi IR descendant a la surface (+ vers le bas) c c Commons: c -------- c c COMMON/IRTAUS/dtaui(nlayer,nspeci) c infrared opacities of the differents layers for differents c spectral ranges. This common is initialized by radtitan. c c COMMON /PLANT/ CSUBP,F0PI c This common is initialized by tgmdat. c c======================================================================= c----------------------------------------------------------------------- c Declarations: c ------------ use dimphy IMPLICIT NONE #include "dimensions.h" #include "YOMCST.h" #include "clesphys.h" INTEGER NLAYER,NSPECI,NSPC1I PARAMETER(NLAYER=llm) PARAMETER (NSPECI=46,NSPC1I=47) c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX INTEGER ngrid PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon c c Arguments: c ---------- INTEGER NG,NL,icld REAL PRESS(NG,NL),TEMP(NG,NL) REAL Z(NG,NL),Q0(NG,NL-1) REAL lwnet(NG,NL),UBARI2 real pfluxi(NG) c Common: c ------- C DTAU IS PASSED EN-MASS, SO ITS DEMENSIONS ARE CRITICAL REAL dtaui(ngrid,NLAYER,NSPECI) REAL dtauip(ngrid,NLAYER,NSPECI) COMMON /IRTAUS/ dtaui,dtauip COMMON /PLANT/ CSUBP,F0PI REAL CSUBP,F0PI c Local: c ------ REAL WNOI(NSPECI),DWNI(NSPECI) ! SPECTAL INTERVALS REAL B0(ngrid,llm+1) REAL EM(ngrid,llm+1) REAL DW,WAVEN,TJ,BSURF,QOUT,QIN,eff_g,COLDEN INTEGER ig,K,J,I,L c EXTERNAL PLNCK REAL PLNCK,zz1,zz2,zz3,zz4,WAVNUM,Xtest REAL FNETIS(ngrid,llm+1),FNETI(ngrid,llm+1) REAL FDIS(ngrid,llm+1,nspeci),FUPIS(ngrid,llm+1,nspeci) REAL FDI(ngrid,llm+1), FUPI(ngrid,llm+1) c Data: c ----- REAL RHOP,UBARI DATA RHOP/1.E4/ ! CONVERSION FROM PRESSURE TO MASS DATA UBARI/0.5/ ! MEAN COSINE FOR 2-STREAM DATA WNOI/ & 11.500, 20.000, 31.250, 50.000, 75.000, & 100.000, 125.000, 150.000, 175.000, 200.000, & 225.000, 250.000, 275.000, 300.000, 325.000, & 350.000, 375.000, 400.000, 425.000, 450.000, & 475.000, 500.000, 525.000, 550.000, 575.000, & 600.000, 628.750, 662.838, 681.757, 683.919, & 686.541, 689.623, 692.704, 695.786, 715.141, & 733.836, 735.597, 737.358, 739.119, 742.720, & 748.160, 753.600, 834.560, 917.333, 926.400, & 935.466/ DATA DWNI/ & 7.000, 10.000, 12.500, 25.000, 25.000, & 25.000, 25.000, 25.000, 25.000, 25.000, & 25.000, 25.000, 25.000, 25.000, 25.000, & 25.000, 25.000, 25.000, 25.000, 25.000, & 25.000, 25.000, 25.000, 25.000, 25.000, & 25.000, 32.500, 35.676, 2.162, 2.162, & 3.082, 3.082, 3.082, 3.082, 35.629, & 1.761, 1.761, 1.761, 1.761, 5.440, & 5.440, 5.440, 156.480, 9.067, 9.067, & 9.067/ save RHOP,UBARI,WNOI,DWNI REAL effg ! effg est une fonction(z en m) c----------------------------------------------------------------------- c Initialisations: c ---------------- UBARI2=1./1.66 UBARI2=UBARI C ZERO THE NET FLUXES Q0 = 0.0 lwnet = 0.0 c----------------------------------------------------------------------- C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE INFRARED C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL c----------------------------------------------------------------------- DO 2000 K=1,NSPECI ! *** START OF SPECTRAL LOOP c----------------------------------------------------------------------- C SET UP ALTITIDUE PARAMETERS WAVEN=WNOI(K) DW=DWNI(K) zz1=DW/(2.*2) EM = 0. B0 = 0. DO J=1,NL-1 DO ig=1,NG TJ=TEMP(ig,J) C Modif: in-lining de la fonction planck pour vectorisation C B0(ig,J)=PLNCK(WAVEN,TJ,DW) C FUNCTION PLNCK(WAV,T,DW) C* PLNCK FUNCTION RETURNS B IN CGS UNITS, ERGS CM-2 WAVENUMBER-1 C* WAVNUM IS WAVENUMBER IN CM-1 C* T IS IN KELVIN PLNCK=0. DO I=-2,2,1 WAVNUM=WAVEN + I*zz1 zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,J)) zz3=WAVNUM*WAVNUM*WAVNUM PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2) ENDDO B0(ig,J)=.2*PLNCK ENDDO IF (ICLD.EQ.1) THEN DO ig=1,NG zz4=EXP(-DTAUI(ig,J,K)/UBARI2) EM(ig,J)=zz4 ENDDO ELSE DO ig=1,NG zz4=EXP(-DTAUIP(ig,J,K)/UBARI2) EM(ig,J)=zz4 ENDDO ENDIF ENDDO c----------------------------------------------------------------------- C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY FDI =0. FDIS =0. FUPI =0. FUPIS=0. DO 2220 J=1,NL-1 DO 2230 ig=1,NG FDI(ig,J+1) = FDI(ig,J)*EM(ig,J) + 2.*RPI*UBARI* & B0(ig,J)*(1.-EM(ig,J)) FDIS(ig,J+1,K) = FDIS(ig,J,K)*EM(ig,J) + 2.*RPI*UBARI* & B0(ig,J)*(1.-EM(ig,J)) 2230 CONTINUE 2220 CONTINUE c write(*,*) c write(*,*) 'cooling : EM =' , c & ((EM(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : B0 =' , c & ((B0(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : FDI =' , c & ((FDI(i,l),l=1,nl),i=1,ngrid) c----------------------------------------------------------------------- C UPWARD FLUXES: SURFACE EMISSIONS DO 2310 ig=1,NG PLNCK=0. DO I=-2,2,1 WAVNUM=WAVEN + I*zz1 zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NL)) zz3=WAVNUM*WAVNUM*WAVNUM PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2) ENDDO c BSURF=PLNCK( WAVEN, TEMP(ig,NL), DW) BSURF=.2*PLNCK*emis FUPI(ig,NL) =BSURF*2.*RPI*UBARI+(1-emis)*FDI(ig,NL) FUPIS(ig,NL,K)=BSURF*2.*RPI*UBARI+(1-emis)*FDIS(ig,NL,K) 2310 CONTINUE c write(*,*) c write(*,*) 'cooling : FUPI/NL =' , c & ((FUPI(i,l),l=nl,nl),i=1,NG) c write(*,*) c write(*,*) 'cooling : FDI/NL =' , c & ((FDI(i,l),l=nl,nl),i=1,NG) DO 2320 J=NL-1,1,-1 DO 2330 ig=1,NG FUPI(ig,J) = FUPI(ig,J+1)*EM(ig,J) + 2.*RPI*UBARI* & B0(ig,J)*(1.-EM(ig,J)) FUPIS(ig,J,K) = FUPIS(ig,J+1,K)*EM(ig,J)+2.*RPI*UBARI* & B0(ig,J)*(1.-EM(ig,J)) 2330 CONTINUE 2320 CONTINUE c write(*,*) c write(*,*) 'cooling : EM =' , c & ((EM(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : B0 =' , c & ((B0(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : FUPI =' , c & ((FUPI(i,l),l=1,nl),i=1,ngrid) c compute the downward IR flux at the surface: c DO 3520 ig=1,NG pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NL) 3520 CONTINUE c compute the net IR flux, (+) upward: c DO J=1,NL DO ig=1,NG lwnet(ig,J)= lwnet(ig,J)+ DWNI(K)*(FUPI(ig,J)-FDI(ig,J)) ENDDO ENDDO DO 3210 J=1,NL-1 DO 3220 ig=1,NG QOUT=FUPI(ig,J) + FDI(ig,J+1) ! OUT OF LAYER QIN =FDI(ig,J) + FUPI(ig,J+1) ! INTO LAYER Q0(ig,J)=Q0(ig,J)+(QOUT-QIN)*DWNI(K) 3220 CONTINUE 3210 CONTINUE c write(*,*) c write(*,*) 'cooling/loop : FUPI =' , c & ((FUPI(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : FDI =' , c & ((FDI(i,l),l=1,nl),i=1,ngrid) c write(*,*) c write(*,*) 'cooling : Q0 =' , c & ((Q0(i,l),l=1,nl-1),i=1,ngrid) c----------------------------------------------------------------------- 2000 CONTINUE ! *** END SPECTRAL INTERVAL COMPUTATIONS c----------------------------------------------------------------------- c convertion erg/cm2 -> J/m2 DO 3550 ig=1,NG pfluxi(ig) = 1.e-3*pfluxi(ig) lwnet(ig,:) = 1.e-3*lwnet(ig,:) 3550 CONTINUE c PRINT*,'flux IR' c WRITE(*,'(8e10.2)') pfluxi C COMPUTE THE BASELINE COOLING RATE DO 3000 J=1,NL-1 C TURN THE Q'S INTO TIMESCALES..... DO 3300 ig=1,NG COLDEN = RHOP*(PRESS(ig,J+1)-PRESS(ig,J))/effg(Z(ig,J)) c Q0(J) = (COLDEN * CSUBP )/Q0(J) Q0(ig,J) = Q0(ig,J) / (COLDEN*CSUBP) 3300 CONTINUE 3000 CONTINUE c write(*,*) c write(*,*) 'cooling/end : Q0 =' c write(*,*) ((Q0(k,l)*1e7,l=1,nl-1),k=1,ngrid) c----------------------------------------------------------------------- RETURN END