source: trunk/LMDZ.TITAN/libf/phytitan/radtitan.F @ 778

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

S.LEBONNOIS:

  • Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
  • Quelques ajouts dans la doc.
File size: 8.7 KB
RevLine 
[104]1       SUBROUTINE RADTITAN(p,nq,nmicro,ycomp,qaer)
[3]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
[102]33      USE infotrac
34      use dimphy
35      USE comgeomphy
36      IMPLICIT NONE
[3]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
[104]51
[3]52c
[104]53c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
54      INTEGER   ngrid
55      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
56c
[3]57
58c   Arguments:
59c   ----------
60
61      INTEGER nq,nmicro 
62
63      REAL p(klon,nlevel)
64      REAL ycomp(klon,nlayer,nq)
[104]65      REAL qaer(klon,klev,nq)
[3]66
67c   Local:
68c   ------
69
70      INTEGER I,J,IG,K,IPRINT
71      INTEGER IPREM
72      LOGICAL notfirstcall
73      SAVE IPREM,notfirstcall
74      data notfirstcall/.false./
75
[104]76      REAL emu,somcoslat,coslat(ngrid)
[3]77 
78      REAL PCH4, effg,FH2L,RHCH4L,SSUM    ! effg est une fonction(z)
79
80c   COMMONS for interface with local subroutines:
81c   ---------------------------------------------
82
83      REAL DZED(NLAYER)
84      REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
85      REAL  CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
86      REAL  XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
87      REAL  C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER)
88      REAL  RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
[175]89      REAL RGAS,RHOP,PI,SIGMA
[3]90
91      COMMON /VERTICAL/ DZED
92
93      COMMON /ATM/ Z,PRESS
94     &            ,DEN,TEMP
95
96
97      COMMON /GASS/ CH4,XN2
98     &              ,H2,AR
99     &              ,XMU,GAS1
100     &              ,COLDEN
101
102      COMMON /STRATO/ C2H2,C2H6
103      COMMON /STRAT2/ HCN
104
105      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
[175]106      COMMON /CONST/RGAS,RHOP,PI,SIGMA
[3]107
108c-----------------------------------------------------------------------
109c   1. Initialisations:
110c   -------------------
111
112       REAL xpoub,kkk,xvis,xir
113
114
115C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
116C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
117      IPRINT=1
118
119C MODIFY ADJUSTABLE NUMBERS HERE -- NOT IN COMMON
120C&&
121      FHAZE=0.3
122C&&
123       if(iprem.eq.0) then
124         TAUFAC=0
125         FHVIS=2.0
126         FHIR=.2
127       print*,'ouverture du fichier initpar'
128       open (unit=1,file='initpar')
129       read(1,*) xpoub,kkk,xvis,xir
130       close(1)
131         FHVIS= xvis
132         FHIR = xir
133       print*,'ouverture du fichier initpar ok'
134       print*,'DANS RADTITAN'
135       print*,'-------------'
136       print*,'FHVIS  = ',FHVIS
137       print*,'FHIR   = ',FHIR
138       iprem=1
139       endif
140
141c-----------------------------------------------------------------------
142c   2. Calcul of the atmospheric profile:
143c   -------------------------------------
144
145       print*,'dans radtitan ',klon
146       print*,notfirstcall
147       IF(notfirstcall) GOTO 300           !F au premier appel!
148       print*,notfirstcall
149
150      DO 210 J=1,NLEVEL
151         PRESS(J)=SSUM(klon,p(1,j),1)/FLOAT(klon)
152210   CONTINUE
153
[106]154c  a cause du tableau predefini dans lell.F (et lell_light.F)
155c     IF(press(nlevel-1).GE.1.44) then
156      IF(press(nlevel-1).GE.1.48) then
[3]157           STOP'pression au sol trop grande'
158          PRINT*,'pression au sol trop grande'
159      endif
160
[106]161c      PRESS(nlevel)=1.48
162c      XCORR=1.48/PRESS(nlevel)
[3]163c     DO 211 J=1,NLEVEL
164c        PRESS(J)=XCORR*PRESS(J)
165c11   CONTINUE
166
167c *********************************************************
168c + 20/1/00: S.Lebonnois: model with chemistry
169c   ++ 22/07/02: ajout HCN ++
170c *********************************************************
171      if (ylellouch) then
172c------------------------------------------------------
173c  initialisation de l'atmosphere et de la composition
174c------------------------------------------------------
175          CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
176     &              CH4,H2,XN2,AR,IPRINT)
177 
178          print*,'LELLOUCH'
179          do i=1,55
180             print*,z(i),PRESS(i)
181          enddo
182C
183C
184C    NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
185          CALL GASSES(IPRINT)
186         
187      else
188c------------------------------------------------------
189c  initialisation seulement de l'atmosphere
190c------------------------------------------------------
191          CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
192     &              CH4,H2,XN2,AR,IPRINT)
193 
194          print*,'LELLOUCH LIGHT'
195          do i=1,55
196             print*,z(i),PRESS(i)
197          enddo
198
199c ++ remplace gasses.F ++
200
201          do i=1,nq
[102]202             if (tname(i).eq."CH4") then
[3]203                iradch4=i
[102]204             elseif (tname(i).eq."C2H2") then
[3]205                iradc2h2=i
[102]206             elseif (tname(i).eq."C2H6") then
[3]207                iradc2h6=i
[102]208             elseif (tname(i).eq."HCN") then
[3]209                iradhcn=i
[102]210             elseif (tname(i).eq."N2") then
[3]211                iradn2=i
[102]212             elseif (tname(i).eq."H2") then
[3]213                iradh2=i
214             endif
215          enddo
216         
217c          print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
218         
219          print*,' ALT   CH4 mass mixing ratio '
220         
221          somcoslat=0.
222          do j=1,klon
223            coslat(j) = cos(rlatd(j)*RPI/180.)
224            somcoslat=somcoslat+coslat(j)
225          enddo
226          do i=1,nlayer
227             colden(i)=rhop*(press(i+1)-press(i))/effg(z(i))
228             gas1(i)=0.
229             emu=(xmu(i+1)+xmu(i))/2.
230             do j=1,klon
231               gas1(i) = gas1(i) +
232     $            coslat(j)/somcoslat*ycomp(j,i,iradch4)*(16./emu)
233             enddo
234             print*,z(i),gas1(i)
235          enddo
236         
237          RHCH4=0.
238          do j=1,klon
239             RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
240          enddo
241          RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
242          print*,'RHCH4 = ',RHCH4
243
244      endif
245
246c *********************************************************
247
248C
249C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
250      CALL SETSPI(IPRINT)
251      CALL SETSPV(IPRINT)
252C SET UP PIA COEFFICIENTS
253      CALL SETPIA(IPRINT,1)
254
255      IF (TAUFAC .GT. 0.)  CALL  CLD(IPRINT)
256
257C
258C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
259C  INFRARED. AND THEN IN THE VISIBLE.
260
261C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
262C  AND AT EACH CALL OF THE PHYSICS
263
264      print*,'aerosol/gas/cloud properties'
265
[104]266      CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
[3]267      print*,'On sort de optci'
268
269C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
270C  INFRARED. AND THEN IN THE VISIBLE.
271
[104]272        CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
[3]273
274        do j=1,NLAYER
275        DZED(j)=Z(J)-Z(J+1)
276        enddo
277       
278c      print*,wlnv
279c      print*,""
280c      print*,wlni
281c      stop
282
283300   CONTINUE  ! fin notfirstcall
284       
285       
286c -----------------------------
287c on ne recalcule pas optci si microfi=0 et compo lellouch
288c -----------------------------
[175]289      IF ((MICROFI.ge.1).or.(.not.ylellouch)) THEN 
[3]290      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
291       print*,'aerosol/gas/cloud properties'
[104]292       CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
[3]293      ENDIF
294      ENDIF
295     
296c ni optcv si microfi=0
297
[175]298      IF (MICROFI.ge.1) THEN 
[3]299      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
300       print*,'aerosol/gas/cloud properties'
[104]301       CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
[3]302      ENDIF
303      ENDIF
304     
305c -----------------------------
306         if (klon.eq.1) then
307           ig=1
308         else
309           ig=klon/2
310         endif
311c       print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
312c       print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
313c       print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
314c       print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
315c       print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
316c       print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
317c       stop
318
319      notfirstcall=.true.
320
321      RETURN
322 191  FORMAT(F8.2,1P10E10.2)
323 192  FORMAT(a8,1P10E10.2)
324      END
Note: See TracBrowser for help on using the repository browser.