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

Last change on this file since 3709 was 1461, checked in by emillour, 10 years ago

Titan GCM:
Turned the common block "tgmdat.F" into a module "tgmdat_mod.F90".
This fixes issues in "debug" mode with common variables which seemed to not be correctly shared between routines.
EM
File size: 9.9 KB
Line 
1      SUBROUTINE COOLING(NG,NL,PRESS,TEMP,Z,Q0,zlwup,zlwdn,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      zlwup(nl)         up fluxes,   (+) upward
39c      zlwdn(nl)         down fluxes, (+) downward
40c      pfluxi          IR descendant a la surface (+ vers le bas)
41c
42c   Commons:
43c   --------
44c
45c     COMMON/IRTAUS/dtaui(nlayer,nspeci)
46c     infrared opacities of the differents layers for differents
47c     spectral ranges. This common is initialized by radtitan.
48c   
49c     COMMON /PLANT/ CSUBP,F0PI
50c     This common is initialized by tgmdat.
51c
52c=======================================================================
53c-----------------------------------------------------------------------
54c   Declarations:
55c   ------------
56
57      use dimphy
58      use tgmdat_mod, only: CSUBP,F0PI
59      IMPLICIT NONE
60#include "dimensions.h"
61#include "YOMCST.h"
62#include "clesphys.h"
63      INTEGER NLAYER,NSPECI,NSPC1I
64      PARAMETER(NLAYER=llm)
65      PARAMETER (NSPECI=46,NSPC1I=47)
66
67c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
68      INTEGER   ngrid
69      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
70c
71c   Arguments:
72c   ----------
73
74      INTEGER NG,NL,icld
75      REAL PRESS(NG,NL),TEMP(NG,NL)
76      REAL Z(NG,NL),Q0(NG,NL-1)
77      REAL zlwup(NG,NL),zlwdn(NG,NL),UBARI2
78      real pfluxi(NG)
79
80
81c   Common:
82c   -------
83
84C DTAU IS PASSED EN-MASS, SO ITS DEMENSIONS ARE CRITICAL
85      REAL dtaui(ngrid,NLAYER,NSPECI)
86      REAL dtauip(ngrid,NLAYER,NSPECI)
87      COMMON /IRTAUS/ dtaui,dtauip
88
89c   Local:
90c   ------
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
99c     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
106c   Data:
107c   -----
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      REAL effg    ! effg est une fonction(z en m)
139
140c-----------------------------------------------------------------------
141
142c   Initialisations:
143c   ----------------
144
145      UBARI2=1./1.66
146      UBARI2=UBARI
147
148C ZERO THE FLUXES
149         Q0    = 0.0
150         zlwup = 0.0
151         zlwdn = 0.0
152
153c-----------------------------------------------------------------------
154C WE NOW ENTER A MAJOR LOOP OVER SPECRAL INTERVALS IN THE INFRARED
155C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL
156c-----------------------------------------------------------------------
157
158       DO 2000 K=1,NSPECI    ! *** START OF SPECTRAL LOOP
159
160c-----------------------------------------------------------------------
161C 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
169          DO J=1,NL-1
170             DO ig=1,NG
171                TJ=TEMP(ig,J)
172
173   
174C  Modif: in-lining de la fonction planck pour vectorisation
175C               B0(ig,J)=PLNCK(WAVEN,TJ,DW)
176C     FUNCTION PLNCK(WAV,T,DW)
177C* PLNCK FUNCTION RETURNS B IN CGS UNITS, ERGS CM-2 WAVENUMBER-1
178C* WAVNUM IS WAVENUMBER IN CM-1
179C* 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=WAVNUM*WAVNUM*WAVNUM
185                   PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2)
186                ENDDO
187                B0(ig,J)=.2*PLNCK
188             ENDDO
189
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
201          ENDDO
202
203c-----------------------------------------------------------------------
204C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL
205C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY
206
207           FDI  =0.
208           FDIS =0.
209           FUPI =0.
210           FUPIS=0.
211
212        DO 2220 J=1,NL-1
213           DO 2230 ig=1,NG
214              FDI(ig,J+1) = FDI(ig,J)*EM(ig,J) + 2.*RPI*UBARI*
215     &        B0(ig,J)*(1.-EM(ig,J))
216              FDIS(ig,J+1,K) = FDIS(ig,J,K)*EM(ig,J) + 2.*RPI*UBARI*
217     &        B0(ig,J)*(1.-EM(ig,J))
2182230       CONTINUE
2192220    CONTINUE
220c     write(*,*)
221c     write(*,*) 'cooling : EM  =' ,
222c    & ((EM(i,l),l=1,nl),i=1,ngrid)
223c     write(*,*)
224c     write(*,*) 'cooling : B0  =' ,
225c    & ((B0(i,l),l=1,nl),i=1,ngrid)
226c     write(*,*)
227c     write(*,*) 'cooling : FDI =' ,
228c    & ((FDI(i,l),l=1,nl),i=1,ngrid)
229
230c-----------------------------------------------------------------------
231C UPWARD FLUXES: SURFACE EMISSIONS
232
233        DO 2310 ig=1,NG
234          PLNCK=0.
235          DO I=-2,2,1
236             WAVNUM=WAVEN + I*zz1
237             zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NL))
238             zz3=WAVNUM*WAVNUM*WAVNUM
239             PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2)
240          ENDDO
241c          BSURF=PLNCK( WAVEN, TEMP(ig,NL), DW)
242           BSURF=.2*PLNCK*emis
243        FUPI(ig,NL)   =BSURF*2.*RPI*UBARI+(1-emis)*FDI(ig,NL)
244        FUPIS(ig,NL,K)=BSURF*2.*RPI*UBARI+(1-emis)*FDIS(ig,NL,K)
2452310    CONTINUE
246c     write(*,*)
247c     write(*,*) 'cooling : FUPI/NL =' ,
248c    & ((FUPI(i,l),l=nl,nl),i=1,NG)
249c     write(*,*)
250c     write(*,*) 'cooling : FDI/NL =' ,
251c    & ((FDI(i,l),l=nl,nl),i=1,NG)
252
253        DO 2320 J=NL-1,1,-1
254           DO 2330 ig=1,NG
255              FUPI(ig,J) = FUPI(ig,J+1)*EM(ig,J) + 2.*RPI*UBARI*
256     &        B0(ig,J)*(1.-EM(ig,J))
257              FUPIS(ig,J,K) = FUPIS(ig,J+1,K)*EM(ig,J)+2.*RPI*UBARI*
258     &        B0(ig,J)*(1.-EM(ig,J))
2592330       CONTINUE
2602320    CONTINUE
261c     write(*,*)
262c     write(*,*) 'cooling : EM  =' ,
263c    & ((EM(i,l),l=1,nl),i=1,ngrid)
264c     write(*,*)
265c     write(*,*) 'cooling : B0  =' ,
266c    & ((B0(i,l),l=1,nl),i=1,ngrid)
267c     write(*,*)
268c     write(*,*) 'cooling : FUPI =' ,
269c    & ((FUPI(i,l),l=1,nl),i=1,ngrid)
270
271c compute the downward IR flux at the surface:
272c
273          DO 3520 ig=1,NG
274             pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NL)
2753520      CONTINUE
276
277c compute the up (+ upward) and down (+ downward) IR fluxes:
278c
279          DO J=1,NL
280          DO ig=1,NG
281             zlwup(ig,J)= zlwup(ig,J)+ DWNI(K)*FUPI(ig,J)
282             zlwdn(ig,J)= zlwdn(ig,J)+ DWNI(K)*FDI(ig,J)
283          ENDDO
284          ENDDO
285         
286          DO 3210 J=1,NL-1
287             DO 3220 ig=1,NG
288                QOUT=FUPI(ig,J) + FDI(ig,J+1)   ! OUT OF LAYER
289                QIN =FDI(ig,J)  + FUPI(ig,J+1)  ! INTO LAYER
290                Q0(ig,J)=Q0(ig,J)+(QOUT-QIN)*DWNI(K)
2913220         CONTINUE
2923210      CONTINUE
293
294c     write(*,*)
295c     write(*,*) 'cooling/loop : FUPI =' ,
296c    & ((FUPI(i,l),l=1,nl),i=1,ngrid)
297c     write(*,*)
298c     write(*,*) 'cooling : FDI  =' ,
299c    & ((FDI(i,l),l=1,nl),i=1,ngrid)
300c     write(*,*)
301c     write(*,*) 'cooling : Q0 =' ,
302c    & ((Q0(i,l),l=1,nl-1),i=1,ngrid)
303
304
305c-----------------------------------------------------------------------
306
3072000  CONTINUE ! *** END SPECTRAL INTERVAL COMPUTATIONS
308
309
310c-----------------------------------------------------------------------
311
312c   convertion erg/cm2 -> J/m2
313      DO 3550 ig=1,NG
314         pfluxi(ig)  = 1.e-3*pfluxi(ig)
315         zlwup(ig,:) = 1.e-3*zlwup(ig,:)
316         zlwdn(ig,:) = 1.e-3*zlwdn(ig,:)
3173550  CONTINUE
318
319c     PRINT*,'flux IR'
320c     WRITE(*,'(8e10.2)') pfluxi
321
322C COMPUTE THE BASELINE COOLING RATE
323
324       DO 3000 J=1,NL-1
325C TURN THE Q'S INTO TIMESCALES.....
326          DO 3300 ig=1,NG
327          COLDEN = RHOP*(PRESS(ig,J+1)-PRESS(ig,J))/effg(Z(ig,J))
328c            Q0(J) = (COLDEN * CSUBP )/Q0(J)
329             Q0(ig,J) = Q0(ig,J) / (COLDEN*CSUBP)
3303300      CONTINUE
3313000  CONTINUE
332
333c     write(*,*)
334c     write(*,*) 'cooling/end : Q0 ='
335c     write(*,*) ((Q0(k,l)*1e7,l=1,nl-1),k=1,ngrid)
336c-----------------------------------------------------------------------
337
338      RETURN
339      END
Note: See TracBrowser for help on using the repository browser.