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

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

SL : mise a jour de phytitan pour etre conforme aux sources actuelles
utilisees sur gnome.

File size: 11.9 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
255c  a cause du tableau predefini dans lell.F (et lell_light.F)
256c     IF(press(nlevel-1).GE.1.44) then
257      IF(press(nlevel-1).GE.1.48) then
258           STOP'pression au sol trop grande'
259          PRINT*,'pression au sol trop grande'
260      endif
261
262c      PRESS(nlevel)=1.48
263c      XCORR=1.48/PRESS(nlevel)
264c     DO 211 J=1,NLEVEL
265c        PRESS(J)=XCORR*PRESS(J)
266c11   CONTINUE
267
268c *********************************************************
269c + 20/1/00: S.Lebonnois: model with chemistry
270c   ++ 22/07/02: ajout HCN ++
271c *********************************************************
272      if (ylellouch) then
273c------------------------------------------------------
274c  initialisation de l'atmosphere et de la composition
275c------------------------------------------------------
276          CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
277     &              CH4,H2,XN2,AR,IPRINT)
278 
279          print*,'LELLOUCH'
280          do i=1,55
281             print*,z(i),PRESS(i)
282          enddo
283C
284C
285C    NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
286          CALL GASSES(IPRINT)
287         
288      else
289c------------------------------------------------------
290c  initialisation seulement de l'atmosphere
291c------------------------------------------------------
292          CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
293     &              CH4,H2,XN2,AR,IPRINT)
294 
295          print*,'LELLOUCH LIGHT'
296          do i=1,55
297             print*,z(i),PRESS(i)
298          enddo
299
300c ++ remplace gasses.F ++
301
302          do i=1,nq
303             if (tname(i).eq."CH4") then
304                iradch4=i
305             elseif (tname(i).eq."C2H2") then
306                iradc2h2=i
307             elseif (tname(i).eq."C2H6") then
308                iradc2h6=i
309             elseif (tname(i).eq."HCN") then
310                iradhcn=i
311             elseif (tname(i).eq."N2") then
312                iradn2=i
313             elseif (tname(i).eq."H2") then
314                iradh2=i
315             endif
316          enddo
317         
318c          print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
319         
320          print*,' ALT   CH4 mass mixing ratio '
321         
322          somcoslat=0.
323          do j=1,klon
324            coslat(j) = cos(rlatd(j)*RPI/180.)
325            somcoslat=somcoslat+coslat(j)
326          enddo
327          do i=1,nlayer
328             colden(i)=rhop*(press(i+1)-press(i))/effg(z(i))
329             gas1(i)=0.
330             emu=(xmu(i+1)+xmu(i))/2.
331             do j=1,klon
332               gas1(i) = gas1(i) +
333     $            coslat(j)/somcoslat*ycomp(j,i,iradch4)*(16./emu)
334             enddo
335             print*,z(i),gas1(i)
336          enddo
337         
338          RHCH4=0.
339          do j=1,klon
340             RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
341          enddo
342          RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
343          print*,'RHCH4 = ',RHCH4
344
345      endif
346
347c *********************************************************
348
349C
350C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
351      CALL SETSPI(IPRINT)
352      CALL SETSPV(IPRINT)
353C SET UP PIA COEFFICIENTS
354      CALL SETPIA(IPRINT,1)
355
356      IF (TAUFAC .GT. 0.)  CALL  CLD(IPRINT)
357
358C
359C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
360C  INFRARED. AND THEN IN THE VISIBLE.
361
362C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
363C  AND AT EACH CALL OF THE PHYSICS
364
365      print*,'aerosol/gas/cloud properties'
366
367      CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
368      print*,'On sort de optci'
369
370C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
371        DO 225 IG=1,klon
372         DO 220 J=1,NLAYER
373          DO 230 K=1,NSPECI
374              DTAUP(IG,J,K)=DTAUI(IG,J,K)
375230       CONTINUE
376220      CONTINUE
377225     CONTINUE
378
379C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
380C  INFRARED. AND THEN IN THE VISIBLE.
381
382        CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
383
384        do j=1,NLAYER
385        DZED(j)=Z(J)-Z(J+1)
386        enddo
387       
388c      print*,wlnv
389c      print*,""
390c      print*,wlni
391c      stop
392
393300   CONTINUE  ! fin notfirstcall
394       
395       
396c -----------------------------
397c on ne recalcule pas optci si microfi=0 et compo lellouch
398c -----------------------------
399      IF ((MICROFI.eq.1).or.(.not.ylellouch)) THEN 
400      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
401       print*,'aerosol/gas/cloud properties'
402       CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
403        DO  IG=1,klon
404         DO  J=1,NLAYER
405          DO  K=1,NSPECI
406              DTAUP(IG,J,K)=DTAUI(IG,J,K)
407          ENDDO
408         ENDDO
409        ENDDO
410      ENDIF
411      ENDIF
412     
413c ni optcv si microfi=0
414
415      IF (MICROFI.eq.1) THEN 
416      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
417       print*,'aerosol/gas/cloud properties'
418       CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
419      ENDIF
420      ENDIF
421     
422c -----------------------------
423         if (klon.eq.1) then
424           ig=1
425         else
426           ig=klon/2
427         endif
428c       print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
429c       print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
430c       print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
431c       print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
432c       print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
433c       print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
434c       stop
435
436      notfirstcall=.true.
437
438      RETURN
439 191  FORMAT(F8.2,1P10E10.2)
440 192  FORMAT(a8,1P10E10.2)
441      END
Note: See TracBrowser for help on using the repository browser.