source: trunk/libf/phytitan/cooling.F @ 119

Last change on this file since 119 was 104, checked in by slebonnois, 14 years ago

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

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