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

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

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

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