Context Navigation

cooling.F @ 1243

Last change on this file since 1243 was 888, checked in by slebonnois, 12 years ago
SL: small modifications to the tools, to Venus default .def files and to outputs (including forgotten modifications linked to the 1D); + bug corrections in phytitan
File size: 9.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format