source: trunk/LMDZ.GENERIC/libf/phystd/optci.F90 @ 3809

Last change on this file since 3809 was 3797, checked in by mturbet, 4 weeks ago

additional corrrection for the new continuum routine
File size: 20.3 KB
Line 
1MODULE optci_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
8     QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
9     TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
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_CO2, igas_O2
16  use comcstfi_mod, only: g, r, mugaz
17  use callkeys_mod, only: kastprof,continuum,graybody,varspec, &
18                          generic_continuum_database
19  use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
20  use tpindex_mod, only: tpindex
21  use interpolate_continuum_mod, only: interpolate_continuum
22
23  implicit none
24
25  !==================================================================
26  !     
27  !     Purpose
28  !     -------
29  !     Calculates longwave optical constants at each level. For each
30  !     layer and spectral interval in the IR it calculates WBAR, DTAU
31  !     and COSBAR. For each level it calculates TAU.
32  !     
33  !     TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
34  !     at the *bottom* of layer L), LW is the spectral wavelength interval.
35  !     
36  !     TLEV(L) - Temperature at the layer boundary (i.e., level)
37  !     PLEV(L) - Pressure at the layer boundary (i.e., level)
38  !
39  !     Authors
40  !     -------
41  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
42  !     
43  !==================================================================
44
45
46  real*8,intent(out) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
47  real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
48  real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
49  real*8,intent(out) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
50  real*8,intent(in) :: PLEV(L_LEVELS)
51  real*8,intent(in) :: TLEV(L_LEVELS) ! not used
52  real*8,intent(in) :: TMID(L_LEVELS)
53  real*8,intent(in) :: PMID(L_LEVELS)
54  real*8,intent(out) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
55  real*8,intent(out) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
56
57  ! for aerosols
58  real*8,intent(in) ::  QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
59  real*8,intent(in) ::  QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
60  real*8,intent(in) ::  GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
61  real*8,intent(in) ::  TAUAERO(L_LEVELS,NAERKIND)
62
63  ! local variables (saved for convenience as need be allocated)
64  real*8,save,allocatable :: TAUAEROLK(:,:,:)
65  real*8,save,allocatable :: TAEROS(:,:,:)
66!$OMP THREADPRIVATE(TAUAEROLK,TAEROS)
67
68  integer L, NW, NG, K, LK, IAER
69  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
70  real*8  ANS, TAUGAS
71  real*8  DPR(L_LEVELS), U(L_LEVELS)
72  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
73
74  real*8,intent(out) :: taugsurf(L_NSPECTI,L_NGAUSS-1)
75  real*8 DCONT,DAERO
76  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
77  double precision p_cross
78
79  ! variable species mixing ratio variables
80  real*8,intent(in) :: QVAR(L_LEVELS)
81  real*8,intent(in) :: MUVAR(L_LEVELS)
82  real*8,intent(in) ::  FRACVAR(ngasmx,L_LEVELS)
83  real*8  WRATIO(L_LEVELS)
84  real*8  KCOEF(4)
85  integer NVAR(L_LEVELS)
86 
87  ! temporary variables to reduce memory access time to gasi
88  real*8 tmpk(2,2)
89  real*8 tmpkvar(2,2,2)
90
91  ! temporary variables for multiple aerosol calculation
92  real*8 atemp
93  real*8 btemp(L_NLAYRAD,L_NSPECTI)
94
95  ! variables for k in units m^-1
96  real*8 dz(L_LEVELS)
97  !real*8 rho !! see test below
98
99  integer igas, jgas
100
101  integer interm
102 
103  logical :: firstcall=.true.
104!$OMP THREADPRIVATE(firstcall)
105
106  !--- Kasting's CIA ----------------------------------------
107  !real*8, parameter :: Ci(L_NSPECTI)=[                         &
108  !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
109  !     5.4E-7, 1.6E-6, 0.0,                               &
110  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            &
111  !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
112  !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
113  !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
114  !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
115  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
116  !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
117  !----------------------------------------------------------
118
119  if (firstcall) then
120    ! allocate local arrays of size "naerkind" (which are also
121    ! "saved" so that this is done only once in for all even if
122    ! we don't need to store the value from a time step to the next)
123    allocate(TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND))
124    allocate(TAEROS(L_LEVELS,L_NSPECTI,NAERKIND))
125    firstcall=.false.
126  endif ! of if (firstcall)
127
128  !! AS: to save time in computing continuum (see bilinearbig)
129  IF (.not.ALLOCATED(indi)) THEN
130      ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
131      indi = -9999  ! this initial value means "to be calculated"
132  ENDIF
133
134  !=======================================================================
135  !     Determine the total gas opacity throughout the column, for each
136  !     spectral interval, NW, and each Gauss point, NG.
137
138  taugsurf(:,:) = 0.0
139  dpr(:)        = 0.0
140  lkcoef(:,:)   = 0.0
141
142  do K=2,L_LEVELS
143     DPR(k) = PLEV(K)-PLEV(K-1)
144
145     !--- Kasting's CIA ----------------------------------------
146     !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
147     ! this is CO2 path length (in cm) as written by Francois
148     ! delta_z = delta_p * R_specific * T / (g * P)
149     ! But Kasting states that W is in units of _atmosphere_ cm
150     ! So we do
151     !dz(k)=dz(k)*(PMID(K)/1013.25)
152     !dz(k)=dz(k)/100.0 ! in m for SI calc
153     !----------------------------------------------------------
154
155     ! if we have continuum opacities, we need dz
156     if(kastprof)then
157        dz(k) = dpr(k)*(1000.0d0*8.314463d0/muvar(k))*TMID(K)/(g*PMID(K))
158        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)
159     else
160        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
161        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
162            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
163     endif
164
165     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
166          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
167
168     do LK=1,4
169        LKCOEF(K,LK) = LCOEF(LK)
170     end do
171  end do                    ! levels
172
173  ! Spectral dependance of aerosol absorption
174  do iaer=1,naerkind
175     DO NW=1,L_NSPECTI
176        do K=2,L_LEVELS
177           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
178        end do                    ! levels
179     END DO
180  end do
181
182  do NW=1,L_NSPECTI
183
184     do K=2,L_LEVELS
185     
186        DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
187
188        DCONT = 0.0d0 ! continuum absorption
189
190        if(continuum.and.(.not.graybody))then
191           ! include continua if necessary
192           
193           if(generic_continuum_database)then
194             T_cont  = dble(TMID(k))
195             do igas=1,ngasmx
196             
197               if(gfrac(igas).eq.-1)then ! variable
198                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
199               elseif(varspec) then
200                 p_cont  = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k)))
201               else
202                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
203               endif
204             
205               do jgas=1,ngasmx
206                 if(gfrac(jgas).eq.-1)then ! variable
207                   p_cross  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
208                 elseif(varspec) then
209                   p_cross  = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
210                 else
211                   p_cross  = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
212                 endif
213               
214                 dtemp=0.0
215
216                 if ( ((igas .eq. igas_N2) .and. (jgas .eq. igas_N2)) .or.   &
217                     ((igas .eq. igas_N2) .and. (jgas .eq. igas_H2)) .or.    &
218                     ((igas .eq. igas_N2) .and. (jgas .eq. igas_O2)) .or.    &
219                     ((igas .eq. igas_N2) .and. (jgas .eq. igas_CH4)) .or.   &
220                     ((igas .eq. igas_O2) .and. (jgas .eq. igas_O2)) .or.    &
221                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_O2)) .or.   &
222                     ((igas .eq. igas_H2) .and. (jgas .eq. igas_H2)) .or.    &
223                     ((igas .eq. igas_H2) .and. (jgas .eq. igas_CH4)) .or.   &
224                     ((igas .eq. igas_H2) .and. (jgas .eq. igas_He)) .or.    &
225                     ((igas .eq. igas_CH4) .and. (jgas .eq. igas_CH4)) .or.  &
226                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_H2O)) .or.  &
227                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_N2)) .or.   &
228                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_O2)) .or.   &
229                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_CO2)) .or.  &
230                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_CO2)) .or.  &
231                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_H2)) .or.   &
232                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_CH4))  ) then
233
234                   call interpolate_continuum('',igas,jgas,'IR',nw,T_cont,p_cont,p_cross,dtemp,.false.)
235
236                 endif
237               
238                 DCONT = DCONT + dtemp
239                 
240               enddo ! jgas=1,ngasmx
241               
242             enddo ! igas=1,ngasmx
243             
244           else ! generic_continuum_database
245           
246            wn_cont = dble(wnoi(nw))
247            T_cont  = dble(TMID(k))
248            do igas=1,ngasmx
249
250              if(gfrac(igas).eq.-1)then ! variable
251                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
252              elseif(varspec) then
253                 p_cont  = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k)))
254              else
255                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
256              endif
257
258              dtemp=0.0d0
259              if(igas.eq.igas_N2)then
260
261                 interm = indi(nw,igas,igas)
262                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
263                 indi(nw,igas,igas) = interm
264
265              elseif(igas.eq.igas_H2)then
266
267                 ! first do self-induced absorption
268                 interm = indi(nw,igas,igas)
269                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
270                 indi(nw,igas,igas) = interm
271
272                 ! then cross-interactions with other gases
273                 do jgas=1,ngasmx
274                    if(varspec) then
275                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
276                    else
277                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
278                    endif
279                    dtempc  = 0.0d0
280                    if(jgas.eq.igas_N2)then
281                       interm = indi(nw,igas,jgas)
282                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
283                       indi(nw,igas,jgas) = interm
284                    elseif(jgas.eq.igas_CO2)then
285                       interm = indi(nw,igas,jgas)
286                       call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
287                       indi(nw,igas,jgas) = interm
288                    elseif(jgas.eq.igas_He)then
289                       interm = indi(nw,igas,jgas)
290                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
291                       indi(nw,igas,jgas) = interm
292                    endif
293                    dtemp = dtemp + dtempc
294                 enddo
295
296              elseif(igas.eq.igas_CH4)then
297
298                 ! first do self-induced absorption
299                 interm = indi(nw,igas,igas)
300                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
301                 indi(nw,igas,igas) = interm
302
303                 ! then cross-interactions with other gases
304                 do jgas=1,ngasmx
305                    if(varspec) then
306                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
307                    else
308                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
309                    endif
310                    dtempc  = 0.0d0
311                    if(jgas.eq.igas_H2)then
312                       interm = indi(nw,igas,jgas)
313                       call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
314                       indi(nw,igas,jgas) = interm
315                    elseif(jgas.eq.igas_CO2)then
316                       interm = indi(nw,igas,jgas)
317                       call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
318                       indi(nw,igas,jgas) = interm
319                    elseif(jgas.eq.igas_He)then
320                       interm = indi(nw,igas,jgas)
321                       call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
322                       indi(nw,igas,jgas) = interm
323                    endif
324                    dtemp = dtemp + dtempc
325                 enddo
326
327              elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
328                 ! Compute self and foreign (with air) continuum of H2O
329                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
330                 interm = indi(nw,igas,igas)
331                 call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
332                 indi(nw,igas,igas) = interm
333
334              endif
335
336              DCONT = DCONT + dtemp
337
338           enddo
339
340           ! Oobleck test
341           !rho = PMID(k)*scalep / (TMID(k)*286.99)
342           !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
343           !   DCONT = rho * 0.125 * 4.6e-4
344           !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
345           !   DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g
346           !   DCONT = rho * 1.0 * 4.6e-4
347           !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
348           !   DCONT = rho * 0.125 * 4.6e-4
349           !endif
350
351           endif ! generic_continuum_database
352           
353           DCONT = DCONT*dz(k)
354       
355        endif ! continuum
356
357        do ng=1,L_NGAUSS-1
358
359           ! Now compute TAUGAS
360
361           ! Interpolate between water mixing ratios
362           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
363           ! the water data range
364
365           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
366           
367              ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
368              ! the execution time of optci/v -> ~ factor 2 on the whole radiative
369              ! transfer on the tested simulations !
370
371              IF (corrk_recombin) THEN ! added by JVO
372                tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
373              ELSE
374                tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
375              ENDIF
376
377              KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
378              KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
379              KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
380              KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
381
382           else
383
384              IF (corrk_recombin) THEN ! added by JVO
385                tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
386              ELSE
387                tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
388              ENDIF
389
390              KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) *  &
391                        ( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
392
393              KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) *  &
394                        ( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
395
396              KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) *  &
397                        ( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
398             
399              KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) *  &
400                        ( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
401
402           endif
403
404           ! Interpolate the gaseous k-coefficients to the requested T,P values
405
406           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
407                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
408
409           TAUGAS  = U(k)*ANS
410
411           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
412           DTAUKI(K,nw,ng) = TAUGAS    &
413                             + DCONT   & ! For parameterized continuum absorption
414                             + DAERO     ! For aerosol absorption
415
416        end do
417
418        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
419        ! which holds continuum opacity only
420
421        NG              = L_NGAUSS
422        DTAUKI(K,nw,ng) = 0.d0      &
423                          + DCONT   & ! For parameterized continuum absorption
424                          + DAERO     ! For aerosol absorption
425
426     end do
427  end do
428
429  !=======================================================================
430  !     Now the full treatment for the layers, where besides the opacity
431  !     we need to calculate the scattering albedo and asymmetry factors
432
433  do iaer=1,naerkind
434    DO NW=1,L_NSPECTI
435     DO K=2,L_LEVELS
436           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
437     ENDDO
438    ENDDO
439  end do
440 
441  DO NW=1,L_NSPECTI
442     DO L=1,L_NLAYRAD-1
443        K              = 2*L+1
444        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
445     END DO ! L vertical loop
446     
447     ! Last level
448     L           = L_NLAYRAD
449     K           = 2*L+1   
450     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
451     
452  END DO                    ! NW spectral loop
453 
454
455  DO NW=1,L_NSPECTI
456     NG = L_NGAUSS
457     DO L=1,L_NLAYRAD-1
458
459        K              = 2*L+1
460        DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
461
462        atemp = 0.
463        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
464           do iaer=1,naerkind
465              atemp = atemp +                                     &
466                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
467                   GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
468           end do
469           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
470        else
471           WBARI(L,nw,ng) = 0.0D0
472           DTAUI(L,NW,NG) = 1.0D-9
473        endif
474
475        if(btemp(L,nw) .GT. 0.0d0) then
476           cosbi(L,NW,NG) = atemp/btemp(L,nw)
477        else
478           cosbi(L,NW,NG) = 0.0D0
479        end if
480
481     END DO ! L vertical loop
482     
483     ! Last level
484     
485     L              = L_NLAYRAD
486     K              = 2*L+1
487     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
488
489     atemp = 0.
490     if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
491        do iaer=1,naerkind
492           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
493        end do
494        WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
495     else
496        WBARI(L,nw,ng) = 0.0D0
497        DTAUI(L,NW,NG) = 1.0D-9
498     endif
499
500     if(btemp(L,nw) .GT. 0.0d0) then
501        cosbi(L,NW,NG) = atemp/btemp(L,nw)
502     else
503        cosbi(L,NW,NG) = 0.0D0
504     end if
505     
506
507     ! Now the other Gauss points, if needed.
508
509     DO NG=1,L_NGAUSS-1
510        IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
511
512           DO L=1,L_NLAYRAD-1
513              K              = 2*L+1
514              DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
515
516              if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
517
518                 WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
519
520              else
521                 WBARI(L,nw,ng) = 0.0D0
522                 DTAUI(L,NW,NG) = 1.0D-9
523              endif
524
525              cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
526           END DO ! L vertical loop
527           
528           ! Last level
529           L              = L_NLAYRAD
530           K              = 2*L+1
531           DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
532
533           if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
534
535              WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
536
537           else
538              WBARI(L,nw,ng) = 0.0D0
539              DTAUI(L,NW,NG) = 1.0D-9
540           endif
541
542           cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
543           
544        END IF
545
546     END DO                 ! NG Gauss loop
547  END DO                    ! NW spectral loop
548
549  ! Total extinction optical depths
550
551  DO NG=1,L_NGAUSS       ! full gauss loop
552     DO NW=1,L_NSPECTI       
553        TAUCUMI(1,NW,NG)=0.0D0
554        DO K=2,L_LEVELS
555           TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
556        END DO
557     END DO                 ! end full gauss loop
558  END DO
559
560  ! be aware when comparing with textbook results
561  ! (e.g. Pierrehumbert p. 218) that
562  ! taucumi does not take the <cos theta>=0.5 factor into
563  ! account. It is the optical depth for a vertically
564  ! ascending ray with angle theta = 0.
565
566  !open(127,file='taucum.out')
567  !do nw=1,L_NSPECTI
568  !   write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
569  !enddo
570  !close(127)
571 
572!  print*,'WBARI'
573!  print*,WBARI
574!  print*,'DTAUI'
575!  print*,DTAUI
576!  call abort
577
578end subroutine optci
579
580END MODULE optci_mod
581
Note: See TracBrowser for help on using the repository browser.