source: trunk/LMDZ.TITAN/libf/phytitan/cooling.F @ 828

Last change on this file since 828 was 808, checked in by slebonnois, 12 years ago

SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).
File size: 9.7 KB
Line 
1      SUBROUTINE COOLING(NG,NL,PRESS,TEMP,Z,Q0,lwnet,pfluxi,icld)
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     COMMON /PLANT/ CSUBP,F0PI
49c     This common is initialized by tgmdat.
50c
51c=======================================================================
52c-----------------------------------------------------------------------
53c   Declarations:
54c   ------------
55
56      use dimphy
57      IMPLICIT NONE
58#include "dimensions.h"
59#include "YOMCST.h"
60#include "clesphys.h"
61      INTEGER NLAYER,NSPECI,NSPC1I
62      PARAMETER(NLAYER=llm)
63      PARAMETER (NSPECI=46,NSPC1I=47)
64
65c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
66      INTEGER   ngrid
67      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
68c
69c   Arguments:
70c   ----------
71
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)
77
78
79c   Common:
80c   -------
81
82C DTAU IS PASSED EN-MASS, SO ITS DEMENSIONS ARE CRITICAL
83      REAL dtaui(ngrid,NLAYER,NSPECI)
84      REAL dtauip(ngrid,NLAYER,NSPECI)
85      COMMON /IRTAUS/ dtaui,dtauip
86
87      COMMON /PLANT/ CSUBP,F0PI
88      REAL CSUBP,F0PI
89
90c   Local:
91c   ------
92
93      REAL WNOI(NSPECI),DWNI(NSPECI)   ! SPECTAL INTERVALS
94      REAL B0(ngrid,llm+1)
95      REAL EM(ngrid,llm+1)
96      REAL DW,WAVEN,TJ,BSURF,QOUT,QIN,eff_g,COLDEN
97
98      INTEGER ig,K,J,I,L
99
100c     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
107c   Data:
108c   -----
109
110      REAL RHOP,UBARI
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
137      save RHOP,UBARI,WNOI,DWNI
138
139c-----------------------------------------------------------------------
140
141c   Initialisations:
142c   ----------------
143
144      UBARI2=1./1.66
145      UBARI2=UBARI
146
147C ZERO THE NET FLUXES
148         Q0    = 0.0
149         lwnet = 0.0
150
151c-----------------------------------------------------------------------
152C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE INFRARED
153C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
154c-----------------------------------------------------------------------
155
156       DO 2000 K=1,NSPECI    ! *** START OF SPECTRAL LOOP
157
158c-----------------------------------------------------------------------
159C SET UP ALTITIDUE PARAMETERS
160
161          WAVEN=WNOI(K)
162          DW=DWNI(K)
163          zz1=DW/(2.*2)
164          EM = 0.
165          B0 = 0.
166
167          DO J=1,NL-1
168             DO ig=1,NG
169                TJ=TEMP(ig,J)
170
171   
172C  Modif: in-lining de la fonction planck pour vectorisation
173C               B0(ig,J)=PLNCK(WAVEN,TJ,DW)
174C     FUNCTION PLNCK(WAV,T,DW)
175C* PLNCK FUNCTION RETURNS B IN CGS UNITS, ERGS CM-2 WAVENUMBER-1
176C* WAVNUM IS WAVENUMBER IN CM-1
177C* T IS IN KELVIN
178                PLNCK=0.
179                DO I=-2,2,1
180                   WAVNUM=WAVEN + I*zz1
181                   zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,J))
182                   zz3=WAVNUM*WAVNUM*WAVNUM
183                   PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2)
184                ENDDO
185                B0(ig,J)=.2*PLNCK
186             ENDDO
187
188             IF (ICLD.EQ.1) THEN
189               DO ig=1,NG
190                 zz4=EXP(-DTAUI(ig,J,K)/UBARI2)
191                 EM(ig,J)=zz4
192               ENDDO
193             ELSE
194               DO ig=1,NG
195                 zz4=EXP(-DTAUIP(ig,J,K)/UBARI2)
196                 EM(ig,J)=zz4
197               ENDDO
198             ENDIF
199          ENDDO
200
201c-----------------------------------------------------------------------
202C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
203C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
204
205           FDI  =0.
206           FDIS =0.
207           FUPI =0.
208           FUPIS=0.
209
210        DO 2220 J=1,NL-1
211           DO 2230 ig=1,NG
212              FDI(ig,J+1) = FDI(ig,J)*EM(ig,J) + 2.*RPI*UBARI*
213     &        B0(ig,J)*(1.-EM(ig,J))
214              FDIS(ig,J+1,K) = FDIS(ig,J,K)*EM(ig,J) + 2.*RPI*UBARI*
215     &        B0(ig,J)*(1.-EM(ig,J))
2162230       CONTINUE
2172220    CONTINUE
218c     write(*,*)
219c     write(*,*) 'cooling : EM  =' ,
220c    & ((EM(i,l),l=1,nl),i=1,ngrid)
221c     write(*,*)
222c     write(*,*) 'cooling : B0  =' ,
223c    & ((B0(i,l),l=1,nl),i=1,ngrid)
224c     write(*,*)
225c     write(*,*) 'cooling : FDI =' ,
226c    & ((FDI(i,l),l=1,nl),i=1,ngrid)
227
228c-----------------------------------------------------------------------
229C UPWARD FLUXES: SURFACE EMISSIONS
230
231        DO 2310 ig=1,NG
232          PLNCK=0.
233          DO I=-2,2,1
234             WAVNUM=WAVEN + I*zz1
235             zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NL))
236             zz3=WAVNUM*WAVNUM*WAVNUM
237             PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2)
238          ENDDO
239c          BSURF=PLNCK( WAVEN, TEMP(ig,NL), DW)
240           BSURF=.2*PLNCK*emis
241        FUPI(ig,NL)   =BSURF*2.*RPI*UBARI+(1-emis)*FDI(ig,NL)
242        FUPIS(ig,NL,K)=BSURF*2.*RPI*UBARI+(1-emis)*FDIS(ig,NL,K)
2432310    CONTINUE
244c     write(*,*)
245c     write(*,*) 'cooling : FUPI/NL =' ,
246c    & ((FUPI(i,l),l=nl,nl),i=1,NG)
247c     write(*,*)
248c     write(*,*) 'cooling : FDI/NL =' ,
249c    & ((FDI(i,l),l=nl,nl),i=1,NG)
250
251        DO 2320 J=NL-1,1,-1
252           DO 2330 ig=1,NG
253              FUPI(ig,J) = FUPI(ig,J+1)*EM(ig,J) + 2.*RPI*UBARI*
254     &        B0(ig,J)*(1.-EM(ig,J))
255              FUPIS(ig,J,K) = FUPIS(ig,J+1,K)*EM(ig,J)+2.*RPI*UBARI*
256     &        B0(ig,J)*(1.-EM(ig,J))
2572330       CONTINUE
2582320    CONTINUE
259c     write(*,*)
260c     write(*,*) 'cooling : EM  =' ,
261c    & ((EM(i,l),l=1,nl),i=1,ngrid)
262c     write(*,*)
263c     write(*,*) 'cooling : B0  =' ,
264c    & ((B0(i,l),l=1,nl),i=1,ngrid)
265c     write(*,*)
266c     write(*,*) 'cooling : FUPI =' ,
267c    & ((FUPI(i,l),l=1,nl),i=1,ngrid)
268
269c compute the downward IR flux at the surface:
270c
271          DO 3520 ig=1,NG
272             pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NL)
2733520      CONTINUE
274
275c compute the net IR flux, (+) upward:
276c
277          DO J=1,NL
278          DO ig=1,NG
279             lwnet(ig,J)= lwnet(ig,J)+ DWNI(K)*(FUPI(ig,J)-FDI(ig,J))
280          ENDDO
281          ENDDO
282         
283          DO 3210 J=1,NL-1
284             DO 3220 ig=1,NG
285                QOUT=FUPI(ig,J) + FDI(ig,J+1)   ! OUT OF LAYER
286                QIN =FDI(ig,J)  + FUPI(ig,J+1)  ! INTO LAYER
287                Q0(ig,J)=Q0(ig,J)+(QOUT-QIN)*DWNI(K)
2883220         CONTINUE
2893210      CONTINUE
290
291c     write(*,*)
292c     write(*,*) 'cooling/loop : FUPI =' ,
293c    & ((FUPI(i,l),l=1,nl),i=1,ngrid)
294c     write(*,*)
295c     write(*,*) 'cooling : FDI  =' ,
296c    & ((FDI(i,l),l=1,nl),i=1,ngrid)
297c     write(*,*)
298c     write(*,*) 'cooling : Q0 =' ,
299c    & ((Q0(i,l),l=1,nl-1),i=1,ngrid)
300
301
302c-----------------------------------------------------------------------
303
3042000  CONTINUE ! *** END SPECTRAL INTERVAL COMPUTATIONS
305
306
307c-----------------------------------------------------------------------
308
309c   convertion erg/cm2 -> J/m2
310      DO 3550 ig=1,NG
311         pfluxi(ig)  = 1.e-3*pfluxi(ig)
312         lwnet(ig,:) = 1.e-3*lwnet(ig,:)
3133550  CONTINUE
314
315c     PRINT*,'flux IR'
316c     WRITE(*,'(8e10.2)') pfluxi
317
318C COMPUTE THE BASELINE COOLING RATE
319
320       DO 3000 J=1,NL-1
321C TURN THE Q'S INTO TIMESCALES.....
322          DO 3300 ig=1,NG
323             eff_g = RG*(RA/(RA+Z(ig,J)))**2 ! 10% DIFF AT 1 MBAR
324             COLDEN = RHOP*(PRESS(ig,J+1)-PRESS(ig,J))/eff_g
325c            Q0(J) = (COLDEN * CSUBP )/Q0(J)
326             Q0(ig,J) = Q0(ig,J) / (COLDEN*CSUBP)
3273300      CONTINUE
3283000  CONTINUE
329
330c     write(*,*)
331c     write(*,*) 'cooling/end : Q0 ='
332c     write(*,*) ((Q0(k,l)*1e7,l=1,nl-1),k=1,ngrid)
333c-----------------------------------------------------------------------
334
335      RETURN
336      END
Note: See TracBrowser for help on using the repository browser.