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

Last change on this file since 3663 was 3663, checked in by gmilcareck, 4 months ago

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