source: trunk/libf/phytitan/radtitan.F @ 6

Last change on this file since 6 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 11.7 KB
Line 
1       SUBROUTINE RADTITAN(p,nq,nmicro,ycomp)
2
3      IMPLICIT NONE
4
5
6c=======================================================================
7c
8c   Authors: C.P. Mc Kay  01/02/91
9c   -------
10c
11c   Object:  Computation of the solar and infra-red
12c   -------  Opacities    (dans des common...)
13c
14c ON TITAN       ADAPTED FROM BEST.FOR FEB 91
15c                           C.P. McKAY
16c
17c   Arguments:
18c   ----------
19c
20c      Input:
21c      ------
22c
23c        p(klon,nl)    pressure (level)
24c        nq       nombre de traceurs
25c        nmicro   nombre de traceurs microphysiques
26c        ycomp(klon,nlayer,nq)
27c
28c      Output:
29c      -------
30c
31c=======================================================================
32c-----------------------------------------------------------------------
33c   Declarations:
34c   -------------
35
36#include "dimensions.h"
37#include "dimphy.h"
38#include "clesphys.h"
39#include "microtab.h"
40#include "numchimrad.h"
41#include "comgeomphy.h"
42#include "YOMCST.h"
43#include "advtrac.h"     !! pour noms des traceurs
44
45c Pour le CRAY, les block data doivent etre declares external
46c pour etre pris en compte
47      EXTERNAL TGMDAT
48
49      INTEGER NLEVEL,NLAYER,NSPECI,NSPC1I,NSPECV,NSPC1V,NSPV
50      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
51      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
52      PARAMETER (NSPV=21)  ! LDO POUR CALCUL ALBEDO
53c
54
55c   Arguments:
56c   ----------
57
58      INTEGER nq,nmicro 
59
60      REAL p(klon,nlevel)
61      REAL ycomp(klon,nlayer,nq)
62
63c   Local:
64c   ------
65
66      INTEGER I,J,IG,K,IPRINT
67      INTEGER IPREM
68      LOGICAL notfirstcall
69      SAVE IPREM,notfirstcall
70      data notfirstcall/.false./
71
72      REAL emu,somcoslat,coslat(klon)
73 
74      REAL PCH4, effg,FH2L,RHCH4L,SSUM    ! effg est une fonction(z)
75
76c   COMMONS for interface with local subroutines:
77c   ---------------------------------------------
78
79      REAL DTAUP(klon,NLAYER,NSPECI)
80      REAL UBARI,UBARV,UBAR0
81      REAL DZED(NLAYER)
82      REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
83      REAL  DTDP(NLAYER),CONVEQ
84      REAL  CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
85      REAL  XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
86      REAL  C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER)
87      REAL  SOLARF(NSPECV),PEXPON(NSPECV),ATERM(4,NSPECV)
88      INTEGER NTERM(NSPECV)
89      REAL  BTERM(4,NSPECV)
90      REAL  RADIUS(NLAYER), XNUMB(NLAYER),REALI(NSPECI), XIMGI(NSPECI)
91      REAL  REALV(NSPECV), XIMGV(NSPECV)
92      REAL  RADCLD(NLAYER), XNCLD(NLAYER),RCLDI(NSPECI),  XICLDI(NSPECI)
93      REAL  RCLDV(NSPECV),  XICLDV(NSPECV)
94      REAL  TAUHI(klon,NSPECI),TAUCI(klon,NSPECI)
95      REAL  TAUGI(klon,NSPECI), TAUGV(klon,NSPECV)
96      REAL  TAURV(klon,NSPECV),TAUHV(klon,NSPECV)
97      REAL  TAUCV(klon,NSPECV)
98c
99      REAL  DTAUI(klon,NLAYER,NSPECI)
100      REAL  TAUI(klon,NLEVEL,NSPECI)
101      REAL  WBARI(klon,NLAYER,NSPECI)
102      REAL  COSBI(klon,NLAYER,NSPECI)
103      REAL  BWNI(NSPC1I),WNOI(NSPECI)
104      REAL  WLNI(NSPECI),DWNI(NSPECI)
105c
106      REAL DTAUV(klon,NLAYER,NSPECV,4)
107      REAL TAUV(klon,NLEVEL,NSPECV,4)
108      REAL WBARV(klon,NLAYER,NSPECV,4)
109      REAL COSBV(klon,NLAYER,NSPECV,4)
110      REAL BWNV(NSPC1V), WNOV(NSPECV),DWNV(NSPECV), WLNV(NSPECV)
111      REAL FNETV(klon,NLEVEL),FUPV(klon,NLEVEL,NSPECV) 
112      REAL FDV(klon,NLEVEL,NSPECV),FMNETV(klon,NLEVEL)
113      REAL FMUPV(NLEVEL),FMDV(NLEVEL)
114      REAL FNET(klon,NLEVEL),FMNET(klon,NLEVEL)
115      REAL THEAT(klon,NLAYER)
116      REAL CSUBP,RSFI,RSFV,F0PI
117      REAL  RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
118      REAL RGAS,RHOP,PII,SIGMA
119      REAL TIDAL
120
121      COMMON /IRTAUS/ DTAUP
122
123      COMMON /VERTICAL/ DZED
124
125      COMMON /ATM/ Z,PRESS
126     &            ,DEN,TEMP
127
128      COMMON /LAPSE/ DTDP,CONVEQ
129      COMMON /UBARED/ UBARI,UBARV,UBAR0
130
131
132
133      COMMON /GASS/ CH4,XN2
134     &              ,H2,AR
135     &              ,XMU,GAS1
136     &              ,COLDEN
137
138      COMMON /STRATO/ C2H2,C2H6
139      COMMON /STRAT2/ HCN
140
141      COMMON /VISGAS/ SOLARF,NTERM
142     &               ,PEXPON,ATERM
143     &               ,BTERM
144
145      COMMON /AERSOL/ RADIUS, XNUMB
146     &               ,REALI, XIMGI
147     &               ,REALV, XIMGV
148
149      COMMON /CLOUD/ RADCLD, XNCLD
150     &             , RCLDI,  XICLDI
151     &             , RCLDV,  XICLDV
152
153      COMMON /TAUS/   TAUHI,TAUCI,
154     &                TAUGI, TAUGV,
155     &                TAURV,TAUHV,
156     &                TAUCV
157
158*     INFRARED  CHARACTERISTICS
159*------------------------------
160
161
162      COMMON /SPECTI/ BWNI,WNOI
163     &               ,DWNI,WLNI
164
165      COMMON /OPTICI/ DTAUI,
166     &                TAUI,
167     &                WBARI,
168     &                COSBI
169
170
171
172*     VISIBLE  CHARACTERISTICS
173*------------------------------
174 
175
176
177      COMMON /OPTICV/ DTAUV
178     &               ,TAUV
179     &               ,WBARV
180     &               ,COSBV
181
182      COMMON /SPECTV/ BWNV, WNOV
183     &               ,DWNV, WLNV
184
185      COMMON /FLUXvV/ FNETV,     
186     &               FUPV,
187     &               FDV,
188     &               FMNETV
189
190      COMMON /FLUX/  FNET,      FMNET
191     &              ,THEAT
192
193
194      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
195      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
196      COMMON /CONST/RGAS,RHOP,PII,SIGMA
197      COMMON /IO/ TIDAL
198
199* common relatifs aux aerosols 
200* nrad dans microtab.h
201      REAL qaer(klon,nlayer,nqmx),volume(nrad),rayon(nrad),vrat,
202     &      drayon(nrad),dvolume(nrad)
203
204      common/traceurs/qaer
205      common/part/volume,rayon,vrat,
206     &      drayon,dvolume
207c-----------------------------------------------------------------------
208c   1. Initialisations:
209c   -------------------
210
211       REAL xpoub,kkk,xvis,xir
212
213
214C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
215C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
216      IPRINT=1
217
218C MODIFY ADJUSTABLE NUMBERS HERE -- NOT IN COMMON
219C&&
220      FHAZE=0.3
221C&&
222       if(iprem.eq.0) then
223         TAUFAC=0
224         FHVIS=2.0
225         FHIR=.2
226       print*,'ouverture du fichier initpar'
227       open (unit=1,file='initpar')
228       read(1,*) xpoub,kkk,xvis,xir
229       close(1)
230         FHVIS= xvis
231         FHIR = xir
232       print*,'ouverture du fichier initpar ok'
233       print*,'DANS RADTITAN'
234       print*,'-------------'
235       print*,'FHVIS  = ',FHVIS
236       print*,'FHIR   = ',FHIR
237       iprem=1
238       endif
239
240c-----------------------------------------------------------------------
241c   2. Calcul of the atmospheric profile:
242c   -------------------------------------
243
244       print*,'dans radtitan ',klon
245       print*,notfirstcall
246       IF(notfirstcall) GOTO 300           !F au premier appel!
247       print*,notfirstcall
248
249      DO 210 J=1,NLEVEL
250         PRESS(J)=SSUM(klon,p(1,j),1)/FLOAT(klon)
251210   CONTINUE
252
253      IF(press(nlevel-1).GE.1.44) then
254           STOP'pression au sol trop grande'
255          PRINT*,'pression au sol trop grande'
256      endif
257
258c      PRESS(nlevel)=1.44
259c      XCORR=1.44/PRESS(nlevel)
260c     DO 211 J=1,NLEVEL
261c        PRESS(J)=XCORR*PRESS(J)
262c11   CONTINUE
263
264c *********************************************************
265c + 20/1/00: S.Lebonnois: model with chemistry
266c   ++ 22/07/02: ajout HCN ++
267c *********************************************************
268      if (ylellouch) then
269c------------------------------------------------------
270c  initialisation de l'atmosphere et de la composition
271c------------------------------------------------------
272          CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
273     &              CH4,H2,XN2,AR,IPRINT)
274 
275          print*,'LELLOUCH'
276          do i=1,55
277             print*,z(i),PRESS(i)
278          enddo
279C
280C
281C    NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
282          CALL GASSES(IPRINT)
283         
284      else
285c------------------------------------------------------
286c  initialisation seulement de l'atmosphere
287c------------------------------------------------------
288          CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
289     &              CH4,H2,XN2,AR,IPRINT)
290 
291          print*,'LELLOUCH LIGHT'
292          do i=1,55
293             print*,z(i),PRESS(i)
294          enddo
295
296c ++ remplace gasses.F ++
297
298          do i=1,nq
299             if (tnom(i).eq."CH4") then
300                iradch4=i
301             elseif (tnom(i).eq."C2H2") then
302                iradc2h2=i
303             elseif (tnom(i).eq."C2H6") then
304                iradc2h6=i
305             elseif (tnom(i).eq."HCN") then
306                iradhcn=i
307             elseif (tnom(i).eq."N2") then
308                iradn2=i
309             elseif (tnom(i).eq."H2") then
310                iradh2=i
311             endif
312          enddo
313         
314c          print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
315         
316          print*,' ALT   CH4 mass mixing ratio '
317         
318          somcoslat=0.
319          do j=1,klon
320            coslat(j) = cos(rlatd(j)*RPI/180.)
321            somcoslat=somcoslat+coslat(j)
322          enddo
323          do i=1,nlayer
324             colden(i)=rhop*(press(i+1)-press(i))/effg(z(i))
325             gas1(i)=0.
326             emu=(xmu(i+1)+xmu(i))/2.
327             do j=1,klon
328               gas1(i) = gas1(i) +
329     $            coslat(j)/somcoslat*ycomp(j,i,iradch4)*(16./emu)
330             enddo
331             print*,z(i),gas1(i)
332          enddo
333         
334          RHCH4=0.
335          do j=1,klon
336             RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
337          enddo
338          RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
339          print*,'RHCH4 = ',RHCH4
340
341      endif
342
343c *********************************************************
344
345C
346C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
347      CALL SETSPI(IPRINT)
348      CALL SETSPV(IPRINT)
349C SET UP PIA COEFFICIENTS
350      CALL SETPIA(IPRINT,1)
351
352      IF (TAUFAC .GT. 0.)  CALL  CLD(IPRINT)
353
354C
355C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
356C  INFRARED. AND THEN IN THE VISIBLE.
357
358C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
359C  AND AT EACH CALL OF THE PHYSICS
360
361      print*,'aerosol/gas/cloud properties'
362
363      CALL OPTCI(ycomp,nmicro,IPRINT)        ! #1
364      print*,'On sort de optci'
365
366C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
367        DO 225 IG=1,klon
368         DO 220 J=1,NLAYER
369          DO 230 K=1,NSPECI
370              DTAUP(IG,J,K)=DTAUI(IG,J,K)
371230       CONTINUE
372220      CONTINUE
373225     CONTINUE
374
375C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
376C  INFRARED. AND THEN IN THE VISIBLE.
377
378        CALL OPTCV(nmicro,IPRINT)        ! #2
379
380        do j=1,NLAYER
381        DZED(j)=Z(J)-Z(J+1)
382        enddo
383       
384c      print*,wlnv
385c      print*,""
386c      print*,wlni
387c      stop
388
389300   CONTINUE  ! fin notfirstcall
390       
391       
392c -----------------------------
393c on ne recalcule pas optci si microfi=0 et compo lellouch
394c -----------------------------
395      IF ((MICROFI.eq.1).or.(.not.ylellouch)) THEN 
396      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
397       print*,'aerosol/gas/cloud properties'
398       CALL OPTCI(ycomp,nmicro,IPRINT)        ! #1
399        DO  IG=1,klon
400         DO  J=1,NLAYER
401          DO  K=1,NSPECI
402              DTAUP(IG,J,K)=DTAUI(IG,J,K)
403          ENDDO
404         ENDDO
405        ENDDO
406      ENDIF
407      ENDIF
408     
409c ni optcv si microfi=0
410
411      IF (MICROFI.eq.1) THEN 
412      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
413       print*,'aerosol/gas/cloud properties'
414       CALL OPTCV(nmicro,IPRINT)        ! #2
415      ENDIF
416      ENDIF
417     
418c -----------------------------
419         if (klon.eq.1) then
420           ig=1
421         else
422           ig=klon/2
423         endif
424c       print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
425c       print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
426c       print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
427c       print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
428c       print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
429c       print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
430c       stop
431
432      notfirstcall=.true.
433
434      RETURN
435 191  FORMAT(F8.2,1P10E10.2)
436 192  FORMAT(a8,1P10E10.2)
437      END
Note: See TracBrowser for help on using the repository browser.