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

Last change on this file since 3708 was 3693, checked in by emillour, 3 months ago

Generic PCM:
Minor corrections:

  • about OpenMP in rad_common_h.F90 : (unclosed bracket in the THREADPRIVATE statement)
  • in interpolate_continuum.F90 : "filename" character string size must be specified as '*' (i.e. variable) and not a fixed number

Took the opportunity to also clean-up interpolate_continuum.F90 by
making it a module, adding some intent() statements, specifying
why saved variables are not threadprivate, get rid of useless "return"
statement at the end of routines, etc.
EM
File size: 21.4 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               dtemp=0.0
206               
207               if (igas .eq. igas_N2) then
208                call interpolate_continuum('',igas_N2,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
209                do jgas=1,ngasmx
210                 if (jgas .eq. igas_H2) then
211                  call interpolate_continuum('',igas_N2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
212                 elseif (jgas .eq. igas_O2) then
213                  call interpolate_continuum('',igas_N2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
214                 elseif (jgas .eq. igas_CH4) then
215                  call interpolate_continuum('',igas_N2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
216                 endif
217                enddo
218               elseif (igas .eq. igas_O2) then
219                call interpolate_continuum('',igas_O2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
220                do jgas=1,ngasmx
221                 if (jgas .eq. igas_CO2) then
222                  call interpolate_continuum('',igas_CO2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
223                 endif
224                enddo
225               elseif (igas .eq. igas_H2) then
226                call interpolate_continuum('',igas_H2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
227                do jgas=1,ngasmx
228                 if (jgas .eq. igas_CH4) then
229                  call interpolate_continuum('',igas_H2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
230                 elseif (jgas .eq. igas_He) then
231                  call interpolate_continuum('',igas_H2,igas_He,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
232                 endif
233                enddo   
234               elseif (igas .eq. igas_CH4) then
235                call interpolate_continuum('',igas_CH4,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
236               elseif (igas .eq. igas_H2O) then
237                call interpolate_continuum('',igas_H2O,igas_H2O,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
238                do jgas=1,ngasmx
239                 if (jgas .eq. igas_N2) then
240                  call interpolate_continuum('',igas_H2O,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
241                 elseif (jgas .eq. igas_O2) then
242                  call interpolate_continuum('',igas_H2O,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
243                 elseif (jgas .eq. igas_CO2) then
244                  call interpolate_continuum('',igas_H2O,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
245                 endif
246                enddo   
247               elseif (igas .eq. igas_CO2) then
248                call interpolate_continuum('',igas_CO2,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
249                do jgas=1,ngasmx
250                 if (jgas .eq. igas_H2) then
251                  call interpolate_continuum('',igas_CO2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
252                 elseif (jgas .eq. igas_CH4) then
253                  call interpolate_continuum('',igas_CO2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
254                 endif
255                enddo
256               endif
257               
258               DCONT = DCONT + dtemp
259               
260             enddo ! igas=1,ngasmx
261           else ! generic_continuum_database
262            wn_cont = dble(wnoi(nw))
263            T_cont  = dble(TMID(k))
264            do igas=1,ngasmx
265
266              if(gfrac(igas).eq.-1)then ! variable
267                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
268              elseif(varspec) then
269                 p_cont  = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k)))
270              else
271                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
272              endif
273
274              dtemp=0.0d0
275              if(igas.eq.igas_N2)then
276
277                 interm = indi(nw,igas,igas)
278                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
279                 indi(nw,igas,igas) = interm
280
281              elseif(igas.eq.igas_H2)then
282
283                 ! first do self-induced absorption
284                 interm = indi(nw,igas,igas)
285                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
286                 indi(nw,igas,igas) = interm
287
288                 ! then cross-interactions with other gases
289                 do jgas=1,ngasmx
290                    if(varspec) then
291                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
292                    else
293                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
294                    endif
295                    dtempc  = 0.0d0
296                    if(jgas.eq.igas_N2)then
297                       interm = indi(nw,igas,jgas)
298                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
299                       indi(nw,igas,jgas) = interm
300                    elseif(jgas.eq.igas_CO2)then
301                       interm = indi(nw,igas,jgas)
302                       call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
303                       indi(nw,igas,jgas) = interm
304                    elseif(jgas.eq.igas_He)then
305                       interm = indi(nw,igas,jgas)
306                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
307                       indi(nw,igas,jgas) = interm
308                    endif
309                    dtemp = dtemp + dtempc
310                 enddo
311
312              elseif(igas.eq.igas_CH4)then
313
314                 ! first do self-induced absorption
315                 interm = indi(nw,igas,igas)
316                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
317                 indi(nw,igas,igas) = interm
318
319                 ! then cross-interactions with other gases
320                 do jgas=1,ngasmx
321                    if(varspec) then
322                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
323                    else
324                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
325                    endif
326                    dtempc  = 0.0d0
327                    if(jgas.eq.igas_H2)then
328                       interm = indi(nw,igas,jgas)
329                       call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
330                       indi(nw,igas,jgas) = interm
331                    elseif(jgas.eq.igas_CO2)then
332                       interm = indi(nw,igas,jgas)
333                       call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
334                       indi(nw,igas,jgas) = interm
335                    elseif(jgas.eq.igas_He)then
336                       interm = indi(nw,igas,jgas)
337                       call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
338                       indi(nw,igas,jgas) = interm
339                    endif
340                    dtemp = dtemp + dtempc
341                 enddo
342
343              elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
344                 ! Compute self and foreign (with air) continuum of H2O
345                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
346                 interm = indi(nw,igas,igas)
347                 call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
348                 indi(nw,igas,igas) = interm
349
350              endif
351
352              DCONT = DCONT + dtemp
353
354           enddo
355
356           ! Oobleck test
357           !rho = PMID(k)*scalep / (TMID(k)*286.99)
358           !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
359           !   DCONT = rho * 0.125 * 4.6e-4
360           !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
361           !   DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g
362           !   DCONT = rho * 1.0 * 4.6e-4
363           !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
364           !   DCONT = rho * 0.125 * 4.6e-4
365           !endif
366
367           DCONT = DCONT*dz(k)
368
369        endif ! generic_continuum_database
370       
371        endif ! continuum
372
373        do ng=1,L_NGAUSS-1
374
375           ! Now compute TAUGAS
376
377           ! Interpolate between water mixing ratios
378           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
379           ! the water data range
380
381           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
382           
383              ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
384              ! the execution time of optci/v -> ~ factor 2 on the whole radiative
385              ! transfer on the tested simulations !
386
387              IF (corrk_recombin) THEN ! added by JVO
388                tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
389              ELSE
390                tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
391              ENDIF
392
393              KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
394              KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
395              KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
396              KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
397
398           else
399
400              IF (corrk_recombin) THEN ! added by JVO
401                tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
402              ELSE
403                tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
404              ENDIF
405
406              KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) *  &
407                        ( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
408
409              KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) *  &
410                        ( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
411
412              KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) *  &
413                        ( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
414             
415              KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) *  &
416                        ( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
417
418           endif
419
420           ! Interpolate the gaseous k-coefficients to the requested T,P values
421
422           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
423                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
424
425           TAUGAS  = U(k)*ANS
426
427           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
428           DTAUKI(K,nw,ng) = TAUGAS    &
429                             + DCONT   & ! For parameterized continuum absorption
430                             + DAERO     ! For aerosol absorption
431
432        end do
433
434        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
435        ! which holds continuum opacity only
436
437        NG              = L_NGAUSS
438        DTAUKI(K,nw,ng) = 0.d0      &
439                          + DCONT   & ! For parameterized continuum absorption
440                          + DAERO     ! For aerosol absorption
441
442     end do
443  end do
444
445  !=======================================================================
446  !     Now the full treatment for the layers, where besides the opacity
447  !     we need to calculate the scattering albedo and asymmetry factors
448
449  do iaer=1,naerkind
450    DO NW=1,L_NSPECTI
451     DO K=2,L_LEVELS
452           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
453     ENDDO
454    ENDDO
455  end do
456 
457  DO NW=1,L_NSPECTI
458     DO L=1,L_NLAYRAD-1
459        K              = 2*L+1
460        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
461     END DO ! L vertical loop
462     
463     ! Last level
464     L           = L_NLAYRAD
465     K           = 2*L+1   
466     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
467     
468  END DO                    ! NW spectral loop
469 
470
471  DO NW=1,L_NSPECTI
472     NG = L_NGAUSS
473     DO L=1,L_NLAYRAD-1
474
475        K              = 2*L+1
476        DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
477
478        atemp = 0.
479        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
480           do iaer=1,naerkind
481              atemp = atemp +                                     &
482                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
483                   GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
484           end do
485           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
486        else
487           WBARI(L,nw,ng) = 0.0D0
488           DTAUI(L,NW,NG) = 1.0D-9
489        endif
490
491        if(btemp(L,nw) .GT. 0.0d0) then
492           cosbi(L,NW,NG) = atemp/btemp(L,nw)
493        else
494           cosbi(L,NW,NG) = 0.0D0
495        end if
496
497     END DO ! L vertical loop
498     
499     ! Last level
500     
501     L              = L_NLAYRAD
502     K              = 2*L+1
503     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
504
505     atemp = 0.
506     if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
507        do iaer=1,naerkind
508           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
509        end do
510        WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
511     else
512        WBARI(L,nw,ng) = 0.0D0
513        DTAUI(L,NW,NG) = 1.0D-9
514     endif
515
516     if(btemp(L,nw) .GT. 0.0d0) then
517        cosbi(L,NW,NG) = atemp/btemp(L,nw)
518     else
519        cosbi(L,NW,NG) = 0.0D0
520     end if
521     
522
523     ! Now the other Gauss points, if needed.
524
525     DO NG=1,L_NGAUSS-1
526        IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
527
528           DO L=1,L_NLAYRAD-1
529              K              = 2*L+1
530              DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
531
532              if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
533
534                 WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
535
536              else
537                 WBARI(L,nw,ng) = 0.0D0
538                 DTAUI(L,NW,NG) = 1.0D-9
539              endif
540
541              cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
542           END DO ! L vertical loop
543           
544           ! Last level
545           L              = L_NLAYRAD
546           K              = 2*L+1
547           DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
548
549           if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
550
551              WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
552
553           else
554              WBARI(L,nw,ng) = 0.0D0
555              DTAUI(L,NW,NG) = 1.0D-9
556           endif
557
558           cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
559           
560        END IF
561
562     END DO                 ! NG Gauss loop
563  END DO                    ! NW spectral loop
564
565  ! Total extinction optical depths
566
567  DO NG=1,L_NGAUSS       ! full gauss loop
568     DO NW=1,L_NSPECTI       
569        TAUCUMI(1,NW,NG)=0.0D0
570        DO K=2,L_LEVELS
571           TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
572        END DO
573     END DO                 ! end full gauss loop
574  END DO
575
576  ! be aware when comparing with textbook results
577  ! (e.g. Pierrehumbert p. 218) that
578  ! taucumi does not take the <cos theta>=0.5 factor into
579  ! account. It is the optical depth for a vertically
580  ! ascending ray with angle theta = 0.
581
582  !open(127,file='taucum.out')
583  !do nw=1,L_NSPECTI
584  !   write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
585  !enddo
586  !close(127)
587 
588!  print*,'WBARI'
589!  print*,WBARI
590!  print*,'DTAUI'
591!  print*,DTAUI
592!  call abort
593
594end subroutine optci
595
596END MODULE optci_mod
597
Note: See TracBrowser for help on using the repository browser.