source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/optcv.F90 @ 285

Last change on this file since 285 was 285, checked in by aslmd, 10 years ago

modification proposed by Sandrine Guerlet to overcome the problem in the penultimate layer when calculating visible heating rates
  • Property svn:executable set to *
File size: 11.7 KB
Line 
1SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
2     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
3     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
4
5  use radinc_h
6  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
7  use gases_h
8
9  implicit none
10
11  !==================================================================
12  !     
13  !     Purpose
14  !     -------
15  !     Calculates shortwave optical constants at each level.
16  !     
17  !     Authors
18  !     -------
19  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
20  !     
21  !==================================================================
22  !     
23  !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL 
24  !     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
25  !     LAYER: WBAR, DTAU, COSBAR
26  !     LEVEL: TAU
27  !     
28  !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
29  !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
30  !     
31  !     TLEV(L) - Temperature at the layer boundary
32  !     PLEV(L) - Pressure at the layer boundary (i.e. level)
33  !     GASV(NT,NPS,NW,NG) - Visible k-coefficients
34  !     
35  !-------------------------------------------------------------------
36
37#include "comcstfi.h"
38#include "callkeys.h"
39
40
41  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
42  real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
43  real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
44  real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
45  real*8 PLEV(L_LEVELS)
46  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
47  real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
48  real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
49
50  ! for aerosols
51  real*8  QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
52  real*8  QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
53  real*8  GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
54  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
55  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
56  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
57
58  integer L, NW, NG, K, LK, IAER
59  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
60  real*8  ANS, TAUGAS
61  real*8  TAURAY(L_NSPECTV)
62  real*8  TRAY(L_LEVELS,L_NSPECTV)
63  real*8  TRAYAER
64  real*8  DPR(L_LEVELS), U(L_LEVELS)
65  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
66
67  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
68  real*8 DCONT,DAERO
69  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
70  double precision p_cross
71
72  ! variable species mixing ratio variables
73  real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
74  real*8  KCOEF(4)
75  integer NVAR(L_LEVELS)
76
77  ! temporary variables for multiple aerosol calculation
78  real*8 atemp(L_NLAYRAD,L_NSPECTV)
79  real*8 btemp(L_NLAYRAD,L_NSPECTV)
80  real*8 ctemp(L_NLAYRAD,L_NSPECTV)
81
82  ! variables for k in units m^-1
83  real*8 dz(L_LEVELS)
84
85  integer igas, jgas
86
87  integer interm
88
89  !! AS: to save time in computing continuum (see bilinearbig)
90  IF (.not.ALLOCATED(indv)) THEN
91      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
92      indv = -9999 ! this initial value means "to be calculated"
93  ENDIF
94
95  !=======================================================================
96  !     Determine the total gas opacity throughout the column, for each
97  !     spectral interval, NW, and each Gauss point, NG.
98  !     Calculate the continuum opacities, i.e., those that do not depend on
99  !     NG, the Gauss index.
100
101  taugsurf(:,:) = 0.0
102  dpr(:)        = 0.0
103  lkcoef(:,:)   = 0.0
104
105  do K=2,L_LEVELS
106     DPR(k) = PLEV(K)-PLEV(K-1)
107
108     ! if we have continuum opacities, we need dz
109     if(kastprof)then
110        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
111        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
112     else
113        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
114        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
115            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
116     endif
117
118     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
119          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
120
121     do LK=1,4
122        LKCOEF(K,LK) = LCOEF(LK)
123     end do
124  end do                    ! levels
125
126
127  do iaer=1,naerkind
128     do NW=1,L_NSPECTV
129        do K=2,L_LEVELS
130           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
131        end do                    ! levels
132     end do
133  end do
134  do NW=1,L_NSPECTV
135     do K=2,L_LEVELS
136        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
137     end do                    ! levels
138  end do
139 
140  !     we ignore K=1...
141  do K=2,L_LEVELS
142
143     do NW=1,L_NSPECTV
144
145        TRAYAER = TRAY(K,NW)
146        !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
147        do iaer=1,naerkind
148           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
149        end do
150
151        DCONT = 0.0 ! continuum absorption
152
153        if(continuum.and.(.not.graybody).and.callgasvis)then
154           ! include continua if necessary
155           wn_cont = dble(wnov(nw))
156           T_cont  = dble(TMID(k))
157           do igas=1,ngasmx
158
159              if(gfrac(igas).eq.-1)then ! variable
160                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
161              else
162                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
163              endif
164
165              dtemp=0.0
166              if(igas.eq.igas_N2)then
167
168                 interm = indv(nw,igas,igas)
169!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
170                 indv(nw,igas,igas) = interm
171                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
172
173              elseif(igas.eq.igas_H2)then
174
175                 ! first do self-induced absorption
176                 interm = indv(nw,igas,igas)
177                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
178                 indv(nw,igas,igas) = interm
179
180                 ! then cross-interactions with other gases
181                 do jgas=1,ngasmx
182                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
183                    dtempc  = 0.0
184                    if(jgas.eq.igas_N2)then
185                       interm = indv(nw,igas,jgas)
186                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
187                       indv(nw,igas,jgas) = interm
188                       ! should be irrelevant in the visible
189                    elseif(jgas.eq.igas_He)then
190                       interm = indv(nw,igas,jgas)
191                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
192                       indv(nw,igas,jgas) = interm
193                    endif
194                    dtemp = dtemp + dtempc
195                 enddo
196
197              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
198
199                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
200                 if(H2Ocont_simple)then
201                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
202                 else
203                    interm = indv(nw,igas,igas)
204                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
205                    indv(nw,igas,igas) = interm
206                 endif
207
208              endif
209
210              DCONT = DCONT + dtemp
211
212           enddo
213
214           DCONT = DCONT*dz(k)
215
216        endif
217
218        do ng=1,L_NGAUSS-1
219
220           ! Now compute TAUGAS
221
222           ! Interpolate between water mixing ratios
223           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
224           ! the water data range
225
226           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
227              KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
228              KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
229              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
230              KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
231           else
232
233              KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
234                   (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
235                   GASV(MT(K),MP(K),NVAR(K),NW,NG))
236
237              KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
238                   (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
239                   GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
240
241              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
242                   (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
243                   GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
244
245              KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
246                   (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
247                   GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
248
249           endif
250
251           ! Interpolate the gaseous k-coefficients to the requested T,P values
252
253           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
254                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
255
256           TAUGAS  = U(k)*ANS
257
258           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
259           DTAUKV(K,nw,ng) = TAUGAS & 
260                             + TRAYAER & ! TRAYAER includes all scattering contributions
261                             + DCONT ! For parameterized continuum aborption
262
263        end do
264
265        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
266        ! which holds continuum opacity only
267
268        NG              = L_NGAUSS
269        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
270
271        do iaer=1,naerkind
272           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) +  TAEROS(K,NW,IAER)
273        end do ! a bug was here!
274
275     end do
276  end do
277
278
279  !=======================================================================
280  !     Now the full treatment for the layers, where besides the opacity
281  !     we need to calculate the scattering albedo and asymmetry factors
282
283  do iaer=1,naerkind
284    DO NW=1,L_NSPECTV
285      DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
286           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
287      ENDDO
288    ENDDO
289  end do
290
291  DO NW=1,L_NSPECTV
292     DO L=1,L_NLAYRAD-1
293        K              = 2*L+1
294        atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
295        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
296        ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
297        btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
298        COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
299     END DO ! L vertical loop
300     
301     !last level
302     L              = L_NLAYRAD
303     K              = 2*L+1
304     atemp(L,NW)    = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
305     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
306     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
307     btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
308     COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
309     
310     
311  END DO                    ! NW spectral loop
312
313  DO NG=1,L_NGAUSS
314    DO NW=1,L_NSPECTV
315     DO L=1,L_NLAYRAD-1
316
317        K              = 2*L+1
318        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
319        WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
320
321      END DO ! L vertical loop
322
323        !     No vertical averaging on bottom layer
324
325        L              = L_NLAYRAD
326        K              = 2*L+1
327        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
328
329        WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
330     END DO                 ! NW spectral loop
331  END DO                    ! NG Gauss loop
332
333  ! Total extinction optical depths
334
335  DO NG=1,L_NGAUSS       ! full gauss loop
336     DO NW=1,L_NSPECTV       
337        TAUV(1,NW,NG)=DTAUKV(2,NW,NG)
338        DO L=1,L_NLAYRAD
339           TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
340        END DO
341
342        TAUCUMV(1,NW,NG)=0.0D0
343        DO K=2,L_LEVELS
344           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
345        END DO
346     END DO           
347  END DO                 ! end full gauss loop
348
349
350  return
351
352
353end subroutine optcv
Note: See TracBrowser for help on using the repository browser.