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

Last change on this file since 1525 was 1461, checked in by emillour, 10 years ago

Titan GCM:
Turned the common block "tgmdat.F" into a module "tgmdat_mod.F90".
This fixes issues in "debug" mode with common variables which seemed to not be correctly shared between routines.
EM
File size: 8.6 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      USE optcld, only : iniqcld
37      use moyzon_mod, only:plevmoy
38      USE TGMDAT_MOD, ONLY: RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,
39     &                      RCLOUD,FARGON
40      USE TGMDAT_MOD, ONLY: RHOP
41      IMPLICIT NONE
42#include "dimensions.h"
43#include "clesphys.h"
44#include "microtab.h"
45#include "numchimrad.h"
46#include "YOMCST.h"
47
48c Pour le CRAY, les block data doivent etre declares external
49c pour etre pris en compte
50      EXTERNAL TGMDAT
51
52      INTEGER NLEVEL,NLAYER,NSPECI,NSPC1I,NSPECV,NSPC1V,NSPV
53      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
54      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
55      PARAMETER (NSPV=21)  ! LDO POUR CALCUL ALBEDO
56
57c
58c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
59      INTEGER   ngrid
60      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
61c
62
63c   Arguments:
64c   ----------
65
66      INTEGER nq,nmicro 
67
68      REAL p(klon,nlevel)
69      REAL ycomp(klon,nlayer,nq)
70      REAL qaer(klon,klev,nq)
71
72c   Local:
73c   ------
74
75      INTEGER I,J,IG,K,IPRINT
76      INTEGER IPREM
77      LOGICAL notfirstcall
78      SAVE IPREM,notfirstcall
79      data notfirstcall/.false./
80
81      REAL emu,somcoslat,coslat(ngrid)
82 
83      REAL PCH4, effg,FH2L,RHCH4L,SSUM    ! effg est une fonction(z)
84
85c   COMMONS for interface with local subroutines:
86c   ---------------------------------------------
87
88      REAL DZED(NLAYER)
89      REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
90      REAL  CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
91      REAL  XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
92      REAL  C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER)
93
94      COMMON /VERTICAL/ DZED
95
96      COMMON /ATM/ Z,PRESS
97     &            ,DEN,TEMP
98
99
100      COMMON /GASS/ CH4,XN2
101     &              ,H2,AR
102     &              ,XMU,GAS1
103     &              ,COLDEN
104
105      COMMON /STRATO/ C2H2,C2H6
106      COMMON /STRAT2/ HCN
107
108c-----------------------------------------------------------------------
109c   1. Initialisations:
110c   -------------------
111
112
113C IPRINT CONTOLS OUTPUT AMOUNT:0=IRREDUCIBLE OUTPUT,LESS THAN 1 PAGE
114C PER RUN, 0=MINIMAL OUTPUT, 1=BACKGROUND ATM AND SPEC; 10=FULL DEBUG
115      IPRINT=1
116
117C&&
118      FHAZE=0.3
119C&&
120       if(iprem.eq.0) then
121         TAUFAC=0
122c xvis et xir lus dans physiq.def  (ancien fichier initpar)
123       FHVIS= xvis
124       FHIR = xir
125c      on initialise le paquet optcld
126       if (clouds.eq.1) call iniqcld()
127       iprem=1
128       endif
129
130c-----------------------------------------------------------------------
131c   2. Calcul of the atmospheric profile:
132c   -------------------------------------
133
134       print*,'dans radtitan ',klon
135       print*,notfirstcall
136       IF(notfirstcall) GOTO 300           !F au premier appel!
137       print*,notfirstcall
138
139c   pression moyenne globale
140c   passage au pressions en bar avec indice 1 au sommet.
141c   (similaire zp dans radlwsw)
142      DO 210 J=2,NLEVEL
143         PRESS(J)=plevmoy(NLEVEL+1-j)*1.e-5
144210   CONTINUE
145      PRESS(1) = PRESS(2)*0.001
146
147c  a cause du tableau predefini dans lell.F (et lell_light.F)
148c     IF(press(nlevel-1).GE.1.44) then
149      IF(press(nlevel-1).GE.1.48) then
150           STOP'pression au sol trop grande'
151          PRINT*,'pression au sol trop grande'
152      endif
153
154c      PRESS(nlevel)=1.48
155c      XCORR=1.48/PRESS(nlevel)
156c     DO 211 J=1,NLEVEL
157c        PRESS(J)=XCORR*PRESS(J)
158c11   CONTINUE
159
160c *********************************************************
161c + 20/1/00: S.Lebonnois: model with chemistry
162c   ++ 22/07/02: ajout HCN ++
163c *********************************************************
164      if (ylellouch) then
165c------------------------------------------------------
166c  initialisation de l'atmosphere et de la composition
167c------------------------------------------------------
168          CALL LELL(NLEVEL,Z,RHCH4L,FH2L,FARGON,TEMP,PRESS,DEN,XMU,
169     &              CH4,H2,XN2,AR,IPRINT)
170 
171          print*,'LELLOUCH'
172          do i=1,55
173             print*,z(i),PRESS(i)
174          enddo
175C
176C
177C    NOW CALCULATE THE LAYER AVERAGE GAS MIXING RATIOS.
178          CALL GASSES(IPRINT)
179
180      else
181c------------------------------------------------------
182c  initialisation seulement de l'atmosphere
183c------------------------------------------------------
184          CALL LELL_LIGHT(NLEVEL,Z,FARGON,TEMP,PRESS,DEN,XMU,
185     &              CH4,H2,XN2,AR,IPRINT)
186 
187          print*,'LELLOUCH LIGHT'
188          do i=1,55
189             print*,z(i),PRESS(i)
190          enddo
191
192c ++ remplace gasses.F ++
193
194          do i=1,nq
195             if (tname(i).eq."CH4") then
196                iradch4=i
197             elseif (tname(i).eq."C2H2") then
198                iradc2h2=i
199             elseif (tname(i).eq."C2H6") then
200                iradc2h6=i
201             elseif (tname(i).eq."HCN") then
202                iradhcn=i
203             elseif (tname(i).eq."N2") then
204                iradn2=i
205             elseif (tname(i).eq."H2") then
206                iradh2=i
207             endif
208          enddo
209         
210c          print*,iradch4,iradc2h2,iradc2h6,iradhcn,iradn2,iradh2
211         
212          print*,' ALT   CH4 mass mixing ratio '
213         
214          somcoslat=0.
215          do j=1,klon
216            coslat(j) = cos(rlatd(j)*RPI/180.)
217            somcoslat=somcoslat+coslat(j)
218          enddo
219          do i=1,nlayer
220c attention ici, Z en km doit etre passe en m
221             colden(i)=rhop*(press(i+1)-press(i))/effg(z(i)*1000.)
222             gas1(i)=0.
223             emu=(xmu(i+1)+xmu(i))/2.
224             do j=1,klon
225               gas1(i) = gas1(i) +
226     $            coslat(j)/somcoslat*ycomp(j,i,iradch4)*(16./emu)
227             enddo
228             print*,z(i),gas1(i)
229          enddo
230         
231          RHCH4=0.
232          do j=1,klon
233             RHCH4 = RHCH4 + coslat(j)/somcoslat*ycomp(j,nlayer,iradch4)
234          enddo
235          RHCH4 = RHCH4*press(nlevel)/PCH4(temp(nlevel))
236          print*,'RHCH4 = ',RHCH4
237
238      endif
239
240c *********************************************************
241
242C
243C CALL A ROUTINE THAT SETS UP THE IR SPECTRAL INTERVALS
244      CALL SETSPI(IPRINT)
245      CALL SETSPV(IPRINT)
246C SET UP PIA COEFFICIENTS
247      CALL SETPIA(IPRINT,1)
248
249      IF (TAUFAC .GT. 0.)  CALL  CLD(IPRINT)
250
251C
252C CALL A SUBROUTINE THAT SETS UP THE OPTICAL PROPERTIES IN THE
253C  INFRARED. AND THEN IN THE VISIBLE.
254
255C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
256C  AND AT EACH CALL OF THE PHYSICS
257
258      print*,'aerosol/gas/cloud properties'
259
260      CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
261      print*,'On sort de optci'
262
263C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
264C  INFRARED. AND THEN IN THE VISIBLE.
265
266        CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
267
268        do j=1,NLAYER
269        DZED(j)=Z(J)-Z(J+1)
270        enddo
271       
272c      print*,wlnv
273c      print*,""
274c      print*,wlni
275c      stop
276
277300   CONTINUE  ! fin notfirstcall
278       
279       
280c -----------------------------
281c on ne recalcule pas optci si microfi=0 et compo lellouch
282c -----------------------------
283      IF ((MICROFI.ge.1).or.(.not.ylellouch)) THEN 
284      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
285       print*,'aerosol/gas/cloud properties'
286       CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
287      ENDIF
288      ENDIF
289     
290c ni optcv si microfi=0
291
292      IF (MICROFI.ge.1) THEN 
293      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
294       print*,'aerosol/gas/cloud properties'
295       CALL OPTCV(qaer,nmicro,IPRINT)        ! #2
296      ENDIF
297      ENDIF
298     
299c -----------------------------
300         if (klon.eq.1) then
301           ig=1
302         else
303           ig=klon/2
304         endif
305c       print*,"DTAUI(equateur,:,1)=",DTAUI(ig,:,1)
306c       print*,"DTAUI(equateur,:,10)=",DTAUI(ig,:,10)
307c       print*,"DTAUI(equateur,:,NSPECI)=",DTAUI(ig,:,NSPECI)
308c       print*,"DTAUV(equateur,:,1,2)=",DTAUV(ig,:,1,2)
309c       print*,"DTAUV(equateur,:,10,2)=",DTAUV(ig,:,10,2)
310c       print*,"DTAUV(equateur,:,NSPECV,2)=",DTAUV(ig,:,NSPECV,2)
311c       stop
312
313      notfirstcall=.true.
314
315      RETURN
316 191  FORMAT(F8.2,1P10E10.2)
317 192  FORMAT(a8,1P10E10.2)
318      END
Note: See TracBrowser for help on using the repository browser.