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

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

SL : corrections et modifications dans phytitan correspondant a celles
faites apres compilation Venus. Titan pas encore compile.

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