source: trunk/LMDZ.TITAN/libf/phytitan/optci.F @ 201

Last change on this file since 201 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: 10.1 KB
RevLine 
[104]1      SUBROUTINE OPTCI(ykim,qaer,nmicro,IPRINT)
[102]2      use dimphy
[119]3      use infotrac
[3]4#include "dimensions.h"
5#include "microtab.h"
6#include "numchimrad.h"
7#include "clesphys.h"
8
9c   Arguments:
10c   ---------
[97]11      REAL    ykim(klon,klev,nqtot)
[119]12      real    qaer(klon,klev,nqtot)
[3]13      integer nmicro
14c   ---------
15
[104]16c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
17      INTEGER   ngrid
18      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
19c
[3]20      PARAMETER(NLAYER=llm,NLEVEL=NLAYER+1)
21      PARAMETER (NSPECI=46,NSPC1I=47,NSPECV=24,NSPC1V=25)
22
23      COMMON /ATM/ Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
24
25      COMMON /GASS/ CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
26     & ,XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
27
28      COMMON /STRATO/ C2H2(NLAYER),C2H6(NLAYER)
29      COMMON /STRAT2/ HCN(NLAYER)
30
31      COMMON /AERSOL/ RADIUS(NLAYER), XNUMB(NLAYER)
32     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
33
[175]34      COMMON /CLOUD/
35     &               RCLDI(NSPECI), XICLDI(NSPECI)
36     &             , RCLDV(NSPECV), XICLDV(NSPECV)
37     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
38     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
[3]39
[104]40      COMMON /TAUS/   TAUHI(ngrid,NSPECI),TAUCI(ngrid,NSPECI),
41     &                TAUGI(ngrid,NSPECI),TAURV(ngrid,NSPECV),
42     &                TAUHV(ngrid,NSPECV),TAUCV(ngrid,NSPECV),
43     &                TAUGV(ngrid,NSPECV)
[3]44
[104]45      COMMON /TAUD/   TAUHID(ngrid,NLAYER,NSPECI)
[175]46     &               ,TAUCID(ngrid,NLAYER,NSPECI)
[104]47     &               ,TAUGID(ngrid,NLAYER,NSPECI)
48     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
[175]49     &               ,TAUCVD(ngrid,NLAYER,NSPECV)
[104]50     &               ,TAUGVD(ngrid,NLAYER,NSPECV)
[3]51
52
[104]53      COMMON /OPTICI/ DTAUI(ngrid,NLAYER,NSPECI)
54     &               ,TAUI (ngrid,NLEVEL,NSPECI)
55     &               ,WBARI(ngrid,NLAYER,NSPECI)
56     &               ,COSBI(ngrid,NLAYER,NSPECI)
[175]57     &               ,DTAUIP(ngrid,NLAYER,NSPECI)
58     &               ,TAUIP(ngrid,NLEVEL,NSPECI)
59     &               ,WBARIP(ngrid,NLAYER,NSPECI)
60     &               ,COSBIP(ngrid,NLAYER,NSPECI)
[3]61
62      COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
63     &                DWNI(NSPECI), WLNI(NSPECI)
64
[175]65      REAL DTAUP(ngrid,NLAYER,NSPECI),DTAUPP(ngrid,NLAYER,NSPECI)
66      COMMON /IRTAUS/ DTAUP,DTAUPP
67
[3]68      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
69      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
70      COMMON /CONST/RGAS,RHOP,PI,SIGMA
71      COMMON /part/v,rayon,vrat,dr,dv
72
[175]73c-----Rayons nuages et "composition" de la goutte
74c     sur la grille ...
75      integer ncount(ngrid,NLAYER)
76      real    rmcbar(ngrid,NLAYER)
77      real    xfbar(ngrid,NLAYER,4)
78      COMMON/rnuabar/ncount,rmcbar,xfbar
79
[3]80      DIMENSION PROD(NLEVEL)
81* nrad dans microtab.h
82      real v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
83      real xv1(klev,nspeci),xv2(klev,nspeci)
84      real xv3(klev,nspeci)
85      REAL QF1(nrad,NSPECI),QF2(nrad,NSPECI)
86      REAL QF3(nrad,NSPECI),QF4(nrad,NSPECI)
87      REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
88      REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
89      real emu
90      REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
91     
92      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
93
94      integer iopti,iwarning     ! iopti: premier appel, une seule boucle sur les l.d'o.
[130]95      integer ig,seulmtunpt
[3]96      save iopti,iwarning,seulmtunpt
97      data iopti,iwarning,seulmtunpt/0,0,0/
98
[175]99      real   zqaer_1pt(NLAYER,2*nrad)
100#include "optci_1pt.h"
101
[3]102      character*100 dummy
103      real   dummy2,dummy3
104
105C THE PRESSURE INDUCED TRANSITIONS ARE FROM REGIS
106C THE LAST SEVENTEEN INTERVALS ARE THE BANDS FROM GNF.
107C
108C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE INFRARED
109C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE IR
110C LAYER: WBAR, DTAU, COSBAR
111C LEVEL: TAU
112C
113       print*,'START OPTCI'
114
115c Diagnostic eventuellement:
116c      if (nmicro.gt.0) then
117c      sum=0.
118c      do nng=2,klon
119c        do i=1,klev
120c         do j=1,nmicro
121c          print*,'j,rj',j,rayon(j)
122c          print*,'paer',qaer(nng,i,j)
123c           sum=sum+qaer(nng,i,j)*rayon(j)**3.*1.3333*3.1415*1000.
124c         enddo
125c        enddo
126c        enddo
127c      print*,sum/(klon-1),'SOMME COLONNE/OPTCI'
128c      endif
129
130
131c      do inq=1,nrad
132c          print*,inq,rayon(inq),vrat,qaer(12,25,inq)
133c      enddo
134             
135C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
136c INITIALISATIONS UNE SEULE FOIS
137C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
138
139      if (iopti.eq.0) then
140
[97]141c verif pour taille zqaer_1pt, sachant que si microfi=0 et nqtot=1,
[3]142c il faut quand meme qu on lise la look-up table de dim nrad=10
143c et si microfi=1, on doit avoir nmicro=nrad (dans microtab.h)
[175]144c
145c Nouvelle verif pour nuages :
146c La condition ci-dessus n'est plus realisable !
147c nmicro comprend maintenant aussi des glaces
148c Donc on teste juste que nmicro soit > 2*nrad (ou nrad si on ne fait pas de nuages)
149       if (microfi.ge.1) then
150         if ((clouds.eq.1).and.(nmicro.lt.2*nrad)) then
151           print*,"OPTCI :"
152           print*,"clouds = 1 MAIS nmicro < 2*nrad"
153           print*,"Probleme pour zqaer_1pt dans optci."
154           stop
155         endif
156         if ((clouds.eq.0).and.(nmicro.lt.nrad)) then
157           print*,"OPTCI :"
158           print*,"nmicro < nrad"
159           print*,"Probleme pour zqaer_1pt dans optci."
160           stop
161         endif
[3]162       endif
163
164      DO 420 K=1,NSPECI
165C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE.
166          CALL THOLIN(WLNI(K),TNR,TNI)
167          REALI(K)=TNR
168          XIMGI(K)=TNI*FHIR
169C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
[175]170          CALL LIQCH4(WLNI(K),TNR,TNI)
171          RCLDI(K)=TNR
172          XICLDI(K)=TNI
173          CALL LIQC2H6(WLNI(K),TNR,TNI)
174          RCLDI2(K)=TNR
175          XICLDI2(K)=TNI
[3]176 420  CONTINUE
177
178C
179C ZERO ALL OPTICAL DEPTHS.
180C ??FLAG? FOR SOME APPLCIATIONS THE TOP OPACITY MAY NOT VANISH
181
182c      open  (unit=1,file='xsetupi')
183c      do i=1,klev
184c       read(1,*) a
185c        do j=1,nspeci
186c            read(1,*) xv1(i,j),xv2(i,j),xv3(i,j)
187c        enddo
188c       enddo
189c       close(1)
190
191      endif    ! fin initialisations premier appel
192
193c************************************************************************
194c************************************************************************
195      DO 79 ig=1,klon      ! BOUCLE SUR GRILLE HORIZONTALE   
196c      print*,'ig NEW optci',ig
197c************************************************************************
198c************************************************************************
199
200        if (.not.ylellouch) then
201       
202            XN2(1) = ykim(ig,1,iradn2)
203            CH4(1) = ykim(ig,1,iradch4)
204             H2(1) = ykim(ig,1,iradh2)
205            do j=2,nlayer
206               XN2(j) = (ykim(ig,j,iradn2)+ykim(ig,j-1,iradn2))/2.
207               CH4(j) = (ykim(ig,j,iradch4)+ykim(ig,j-1,iradch4))/2.
208                H2(j) = (ykim(ig,j,iradh2)+ykim(ig,j-1,iradh2))/2.
209            enddo
210            XN2(nlevel) = ykim(ig,nlayer,iradn2)
211            CH4(nlevel) = ykim(ig,nlayer,iradch4)
212             H2(nlevel) = ykim(ig,nlayer,iradh2)     
213
214            do j=1,nlayer
215               emu = ( xmu(j) + xmu(j+1) )/2.
216               C2H2(j) = ykim(ig,j,iradc2h2) * 26./emu
217               C2H6(j) = ykim(ig,j,iradc2h6) * 30./emu
218                HCN(j) = ykim(ig,j,iradhcn ) * 27./emu
219            enddo
220                 
221        endif
222
223c     if ((.not.ylellouch).and.(ig.eq.klon/2)) then
224c        print*,' LAYER      C2H2         C2H6       HCN masmix ratios'
225c        do j=1,nlayer
226c            print*,j,C2H2(j),C2H6(j),HCN(j)
227c        enddo
228c     endif   
229
[175]230        if (microfi.ge.1) then
231          do iq=1,2*nrad
232c           si on ne fait pas de nuages on ne copie que les nrad premieres valeurs.
233            if (clouds.eq.0.and.iq.gt.nrad) then
234              zqaer_1pt(:,iq)=0.
235            else
236              do j=1,NLAYER
237                zqaer_1pt(j,iq)=qaer(ig,j,iq)
238              enddo
239            endif
240          enddo
[3]241        else
242         if (ig.eq.1)  then
[175]243           zqaer_1pt = 0.
[3]244c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
245c boucle sur nrad=10 (dans microtab.h)
246           open(10,file="qaer_eq_1d.dat")
247           do iq=1,15
248             read(10,'(A100)') dummy
249           enddo
250           do j=NLAYER,1,-1
251             read(10,*) dummy2,dummy3,(zqaer_1pt(j,iq),iq=1,nrad)
252           enddo
253           close(10)
254c ici, les tableaux definissant la structure des aerosols sont
255c remplis: rf,df(nq),rayon(nq,)v(nq)......
[175]256           call rdf()
[3]257         endif
258        endif
259
260c        if ((ig.eq.klon/2).or.(microfi.eq.0))  then
261c       print*,"Q01=",zqaer_1pt(:,1)
262c       print*,"Q05=",zqaer_1pt(:,5)
263c       print*,"Q10=",zqaer_1pt(:,10)
264c       stop
265c        endif
266       
267        if (seulmtunpt.eq.0) then
[175]268          call optci_1pt2(zqaer_1pt,rmcbar(ig,:),xfbar(ig,:,:),
269     &                   iopti,IPRINT)
[3]270           iopti = 1
271        endif
272
273c Pas de microphysique, ni de composition variable: un seul passage
274c dans optci_1pt.
275        if ((microfi.eq.0).and.(ylellouch)) then
276           seulmtunpt = 1
277        endif
278       
[175]279        COSBI(ig,:,:)  = COSBI_1pt(:,:)
280        WBARI(ig,:,:)  = WBARI_1pt(:,:)
281        DTAUI(ig,:,:)  = DTAUI_1pt(:,:)
282        TAUI(ig,:,:)   = TAUI_1pt(:,:)
[3]283
[175]284        COSBIP(ig,:,:)  = COSBIP_1pt(:,:)
285        WBARIP(ig,:,:)  = WBARIP_1pt(:,:)
286        DTAUIP(ig,:,:)  = DTAUIP_1pt(:,:)
287        TAUIP(ig,:,:)   = TAUIP_1pt(:,:)
288
289        TAUHI(ig,:)    = TAUHI_1pt(:)
290        TAUCI(ig,:)    = TAUCI_1pt(:)
291        TAUGI(ig,:)    = TAUGI_1pt(:)
292
293        TAUHID(ig,:,:) = TAUHID_1pt(:,:)
294        TAUCID(ig,:,:) = TAUCID_1pt(:,:)
295        TAUGID(ig,:,:) = TAUGID_1pt(:,:)
296
[3]297c************************************************************************
298c************************************************************************
299  79  CONTINUE   ! FIN BOUCLE GRILLE HORIZONTALE
300c************************************************************************
301c************************************************************************
[175]302C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
303        DO 225 IG=1,klon
304         DO 220 J=1,NLAYER
305          DO 230 K=1,NSPECI
306              DTAUP(IG,J,K)=DTAUI(IG,J,K)
307              DTAUPP(IG,J,K)=DTAUIP(IG,J,K)
308230       CONTINUE
309220      CONTINUE
310225     CONTINUE
[3]311
312      print*, 'FIN OPTCI'
313
314      RETURN
315      END
[119]316
Note: See TracBrowser for help on using the repository browser.