source: trunk/LMDZ.PLUTO/libf/phypluto/optci_pluto_mod.F90 @ 3585

Last change on this file since 3585 was 3329, checked in by afalco, 9 months ago

Pluto PCM:
Include cooling, hazes in radiative module
AF
File size: 10.1 KB
RevLine 
[3329]1MODULE optci_pluto_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7      subroutine optci_pluto(PLEV,TLEV,DTAUI,TAUCUMI,      &
8           QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
9           TMID,PMID,TAUGSURF,QVAR)
10
11  use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
12                      L_NLEVRAD, L_REFVAR, naerkind
13  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
14  use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, &
15                     igas_CH4, igas_N2
16  use comcstfi_mod, only: g, r, mugaz
17  use callkeys_mod, only: kastprof,continuum,graybody
18  use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
19  use tpindex_mod, only: tpindex
20      implicit none
21
22!==================================================================
23!
24!     Purpose
25!     -------
26!     Calculates longwave optical constants at each level. For each
27!     layer and spectral interval in the IR it calculates WBAR, DTAU
28!     and COSBAR. For each level it calculates TAU.
29!
30!     TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
31!     at the *bottom* of layer L), LW is the spectral wavelength interval.
32!
33!     TLEV(L) - Temperature at the layer boundary (i.e., level)
34!     PLEV(L) - Pressure at the layer boundary (i.e., level)
35!
36!     Authors
37!     -------
38!     Adapted from the NASA Ames code by R. Wordsworth (2009)
39!
40!==================================================================
41
42
43
44      real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
45      real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
46      real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
47      real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
48      real*8 PLEV(L_LEVELS)
49      real*8 TLEV(L_LEVELS)
50      real*8 TMID(L_LEVELS), PMID(L_LEVELS)
51      real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
52      real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
53
54!     For aerosols
55      real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
56      real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
57      real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
58      real*8  TAUAERO(L_LEVELS+1,NAERKIND)
59      real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
60      real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
61
62      integer L, NW, NG, K, LK, IAER
63      integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
64      real*8  ANS, TAUGAS
65      real*8  DPR(L_LEVELS), U(L_LEVELS)
66      real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
67
68      real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
69      real*8 dco2
70
71!     mixing ratio variables
72      real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS)
73      real*8  KCOEF(4)
74      integer NVAR(L_LEVELS)
75
76!     temporary variables for multiple aerosol calculation
77      real*8 atemp, btemp
78
79!     variables for k in units m^-1
80      real*8 rho, dz
81
82!=======================================================================
83!     Determine the total gas opacity throughout the column, for each
84!     spectral interval, NW, and each Gauss point, NG.
85
86!      write(*,*)'L_LEVELS',L_LEVELS
87!      write(*,*)'L_NSPECTI',L_NSPECTI
88      DTAUI(:,:,:)=0.
89      DTAUKI(:,:,:)=0.
90
91      DO NG=1,L_NGAUSS-1
92         do NW=1,L_NSPECTI
93            TAUGSURF(NW,NG) = 0.0D0
94         end do
95      end do
96      do K=2,L_LEVELS
97         DPR(k) = PLEV(K)-PLEV(K-1)
98
99!         rho  = PLEV(K)/(R*TMID(K))
100         rho  = PMID(K)/(R*TMID(K))
101         dz   = -DPR(k)/(g*rho)
102         !print*,'rho=',rho
103         !print*,'dz=',dz
104
105         U(k)   = Cmk*DPR(k)    ! only Cmk line in optci_pluto.F
106         ! soon to be replaced by m^-1 !!!
107
108         call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
109              LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
110
111         do LK=1,4
112            LKCOEF(K,LK) = LCOEF(LK)
113         end do
114
115         DO NW=1,L_NSPECTI
116            do iaer=1,naerkind
117               TAEROS(K,NW,IAER) = TAUAERO(K,IAER)  * QXIAER(K,NW,IAER)
118!               write(22,*) 'TB17 Taero IR:',K,NW,IAER,TAEROS(K,NW,IAER)
119            end do
120         END DO
121      end do                    ! levels
122
123
124      do K=2,L_LEVELS
125         do nw=1,L_NSPECTI
126
127             DCO2   = 0.0 ! continuum absorption (no longer used)
128
129           do ng=1,L_NGAUSS-1
130
131!     Now compute TAUGAS
132!     Interpolate between mixing ratios
133!     WRATIO = 0.0 if the requested amount is equal to, or outside the
134!     the data range
135
136
137            if (L_REFVAR.eq.1)then ! added by RW for special no variable case
138                  KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
139                  KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
140                  KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
141                  KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
142            else
143
144              KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
145                   (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                &
146                   GASI(MT(K),MP(K),NVAR(K),NW,NG))
147
148              KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)+ WRATIO(K)* &
149                   (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -              &
150                   GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
151
152              KCOEF(3)=GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)+WRATIO(K)*  &
153                   (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -            &
154                   GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
155
156              KCOEF(4) =GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
157                   (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -              &
158                   GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
159            endif
160
161!     Interpolate the gaseous k-coefficients to the requested T,P values
162
163               ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +    &
164                    LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
165
166               TAUGAS  = U(k)*ANS
167
168               TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
169
170               DTAUKI(K,nw,ng) = TAUGAS
171               do iaer=1,naerkind
172                 DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
173                       + DCO2 ! For Kasting CIA
174               end do
175
176           end do
177
178!     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
179!     which holds continuum opacity only
180
181           NG              = L_NGAUSS
182           DTAUKI(K,nw,ng) = 0.0
183           do iaer=1,naerkind
184               DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) +  TAEROS(K,NW,IAER) &
185                    + DCO2 ! For parameterized continuum absorption
186           end do ! a bug was found here!!
187
188         end do
189      end do
190
191
192!=======================================================================
193!     Now the full treatment for the layers, where besides the opacity
194!     we need to calculate the scattering albedo and asymmetry factors
195
196      DO NW=1,L_NSPECTI
197         DO K=2,L_LEVELS
198            do iaer=1,naerkind
199               TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
200            end do
201         ENDDO
202      ENDDO
203
204      DO NW=1,L_NSPECTI
205         NG = L_NGAUSS
206         DO L=1,L_NLAYRAD-1
207
208            K              = 2*L+1
209            DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
210
211            atemp = 0.
212            btemp = 0.
213            if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
214               do iaer=1,naerkind
215                  atemp = atemp +                                     &
216                      GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
217                      GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
218                  btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
219!     *                    + 1.e-10
220               end do
221               WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
222            else
223               WBARI(L,nw,ng) = 0.0D0
224               DTAUI(L,NW,NG) = 1.0E-9
225            endif
226
227            if(btemp .GT. 0.0) then
228               cosbi(L,NW,NG) = atemp/btemp
229            else
230               cosbi(L,NW,NG) = 0.0D0
231            end if
232
233         END DO ! L vertical loop
234
235     ! Last level
236
237     L              = L_NLAYRAD
238     K              = 2*L+1
239     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
240
241     atemp = 0.
242     if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
243        do iaer=1,naerkind
244           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
245        end do
246        WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
247     else
248        WBARI(L,nw,ng) = 0.0D0
249        DTAUI(L,NW,NG) = 1.0D-9
250     endif
251
252     if(btemp .GT. 0.0d0) then
253        cosbi(L,NW,NG) = atemp/btemp
254     else
255        cosbi(L,NW,NG) = 0.0D0
256     end if
257
258!     Now the other Gauss points, if needed.
259
260         DO NG=1,L_NGAUSS-1
261            IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
262
263               DO L=1,L_NLAYRAD-1
264                  K              = 2*L+1
265                  DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
266
267                  btemp = 0.
268                  if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
269
270                     do iaer=1,naerkind
271                        btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
272                     end do
273                     WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
274
275                  else
276                     WBARI(L,nw,ng) = 0.0D0
277                     DTAUI(L,NW,NG) = 1.0E-9
278                  endif
279
280                  cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
281               END DO ! L vertical loop
282
283               ! Last level
284               L           = L_NLAYRAD
285               K           = 2*L+1
286               btemp = SUM(TAUAEROLK(K,NW,1:naerkind))
287               DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
288
289               if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
290
291                    WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
292
293               else
294                    WBARI(L,nw,ng) = 0.0D0
295                    DTAUI(L,NW,NG) = 1.0D-9
296               endif
297
298               cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
299
300            END IF
301
302         END DO                 ! NG Gauss loop
303      END DO                    ! NW spectral loop
304
305!     Total extinction optical depths
306
307      DO NW=1,L_NSPECTI
308         DO NG=1,L_NGAUSS       ! full gauss loop
309            TAUI(1,NW,NG)=0.0D0
310            DO L=1,L_NLAYRAD
311               TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
312            END DO
313
314            TAUCUMI(1,NW,NG)=0.0D0
315            DO K=2,L_LEVELS
316               TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
317            END DO
318         END DO                 ! end full gauss loop
319      END DO
320
321
322end subroutine optci_pluto
323
324END MODULE optci_pluto_mod
325
Note: See TracBrowser for help on using the repository browser.