source: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90 @ 2655

Last change on this file since 2655 was 2655, checked in by gmilcareck, 3 years ago

Major changes to CIA interpolation:
1) Add contribution from CH4 (H2-CH4,He-CH4,CH4-CH4) ;
2) Add some tests before interpolation for H2, He and CH4 ;
3) Add the possibility to choose between a normal or equilibrium ortho:para
fraction for CIA H2;
4) Change "strictboundH2H2cia" to the generic "strictboundcia" for H2,He,CH4.
It can be added for others CIA (N2,H,CO2...) if you want.
  • Property svn:executable set to *
File size: 14.4 KB
Line 
1MODULE optcv_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
8     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
9     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
10
11  use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
12  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
13  use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4
14  use comcstfi_mod, only: g, r, mugaz
15  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
16  use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
17  use tpindex_mod, only: tpindex
18
19  implicit none
20
21  !==================================================================
22  !     
23  !     Purpose
24  !     -------
25  !     Calculates shortwave optical constants at each level.
26  !     
27  !     Authors
28  !     -------
29  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
30  !     
31  !==================================================================
32  !     
33  !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL 
34  !     IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
35  !     LAYER: WBAR, DTAU, COSBAR
36  !     LEVEL: TAU
37  !     
38  !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
39  !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
40  !     
41  !     TLEV(L) - Temperature at the layer boundary
42  !     PLEV(L) - Pressure at the layer boundary (i.e. level)
43  !     GASV(NT,NPS,NW,NG) - Visible k-coefficients
44  !     
45  !-------------------------------------------------------------------
46
47
48  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
49  real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
50  real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
51  real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
52  real*8 PLEV(L_LEVELS)
53  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
54  real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
55  real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
56
57  ! for aerosols
58  real*8  QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
59  real*8  QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
60  real*8  GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
61  real*8  TAUAERO(L_LEVELS,NAERKIND)
62  real*8  TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND)
63  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
64
65  integer L, NW, NG, K, LK, IAER
66  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
67  real*8  ANS, TAUGAS
68  real*8  TAURAY(L_NSPECTV)
69  real*8  TRAY(L_LEVELS,L_NSPECTV)
70  real*8  DPR(L_LEVELS), U(L_LEVELS)
71  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
72
73  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
74  real*8 DCONT,DAERO
75  real*8 DRAYAER
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  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
81  real*8  KCOEF(4)
82  integer NVAR(L_LEVELS)
83 
84  ! temporary variables to reduce memory access time to gasv
85  real*8 tmpk(2,2)
86  real*8 tmpkvar(2,2,2)
87
88  ! temporary variables for multiple aerosol calculation
89  real*8 atemp(L_NLAYRAD,L_NSPECTV)
90  real*8 btemp(L_NLAYRAD,L_NSPECTV)
91  real*8 ctemp(L_NLAYRAD,L_NSPECTV)
92
93  ! variables for k in units m^-1
94  real*8 dz(L_LEVELS)
95
96
97  integer igas, jgas
98
99  integer interm
100
101  !! AS: to save time in computing continuum (see bilinearbig)
102  IF (.not.ALLOCATED(indv)) THEN
103      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
104      indv = -9999 ! this initial value means "to be calculated"
105  ENDIF
106
107  !=======================================================================
108  !     Determine the total gas opacity throughout the column, for each
109  !     spectral interval, NW, and each Gauss point, NG.
110  !     Calculate the continuum opacities, i.e., those that do not depend on
111  !     NG, the Gauss index.
112
113  taugsurf(:,:) = 0.0
114  dpr(:)        = 0.0
115  lkcoef(:,:)   = 0.0
116
117  do K=2,L_LEVELS
118     DPR(k) = PLEV(K)-PLEV(K-1)
119
120     ! if we have continuum opacities, we need dz
121     if(kastprof)then
122        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
123        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)
124     else
125        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
126        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
127            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
128     endif
129
130     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
131          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
132
133     do LK=1,4
134        LKCOEF(K,LK) = LCOEF(LK)
135     end do
136  end do                    ! levels
137
138  ! Spectral dependance of aerosol absorption
139            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
140            !   but visible does not handle very well diffusion in first layer.
141            !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
142            !   in the 4 first semilayers in optcv, but not optci.
143            !   This solves random variations of the sw heating at the model top.
144  do iaer=1,naerkind
145     do NW=1,L_NSPECTV
146        TAEROS(1:4,NW,IAER)=0.d0
147        do K=5,L_LEVELS
148           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
149        end do                    ! levels
150     end do
151  end do
152 
153  ! Rayleigh scattering
154  do NW=1,L_NSPECTV
155     TRAY(1:4,NW)   = 1d-30
156     do K=5,L_LEVELS
157        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
158     end do                    ! levels
159  end do
160 
161  !     we ignore K=1...
162  do K=2,L_LEVELS
163
164     do NW=1,L_NSPECTV
165
166        DRAYAER = TRAY(K,NW)
167        !     DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
168        do iaer=1,naerkind
169           DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
170        end do
171
172        DCONT = 0.0 ! continuum absorption
173
174        if(continuum.and.(.not.graybody).and.callgasvis)then
175           ! include continua if necessary
176           wn_cont = dble(wnov(nw))
177           T_cont  = dble(TMID(k))
178           do igas=1,ngasmx
179
180              if(gfrac(igas).eq.-1)then ! variable
181                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
182              else
183                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
184              endif
185
186              dtemp=0.0
187              if(igas.eq.igas_N2)then
188
189                 interm = indv(nw,igas,igas)
190!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
191                 indv(nw,igas,igas) = interm
192                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
193
194              elseif(igas.eq.igas_H2)then
195
196                 ! first do self-induced absorption
197                 interm = indv(nw,igas,igas)
198                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
199                 indv(nw,igas,igas) = interm
200
201                 ! then cross-interactions with other gases
202                 do jgas=1,ngasmx
203                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
204                    dtempc  = 0.0
205                    if(jgas.eq.igas_N2)then
206                       interm = indv(nw,igas,jgas)
207                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
208                       indv(nw,igas,jgas) = interm
209                       ! should be irrelevant in the visible
210                    elseif(jgas.eq.igas_He)then
211                       interm = indv(nw,igas,jgas)
212                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
213                       indv(nw,igas,jgas) = interm
214                    endif
215                    dtemp = dtemp + dtempc
216                 enddo
217                 
218              elseif(igas.eq.igas_CH4)then
219
220                 ! first do self-induced absorption
221                 interm = indv(nw,igas,igas)
222                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
223                 indv(nw,igas,igas) = interm
224
225                 ! then cross-interactions with other gases
226                 do jgas=1,ngasmx
227                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
228                    dtempc  = 0.0d0
229                    if(jgas.eq.igas_H2)then
230                       interm = indv(nw,igas,jgas)
231                       call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
232                       indv(nw,igas,jgas) = interm
233                    elseif(jgas.eq.igas_He)then
234                       interm = indv(nw,igas,jgas)
235                       call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
236                       indv(nw,igas,jgas) = interm
237                    endif
238                    dtemp = dtemp + dtempc
239                 enddo
240
241              elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
242                 ! Compute self and foreign (with air) continuum of H2O
243                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
244                 interm = indv(nw,igas,igas)
245                 call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
246                 indv(nw,igas,igas) = interm
247
248              endif
249
250              DCONT = DCONT + dtemp
251
252           enddo
253
254           DCONT = DCONT*dz(k)
255
256        endif
257
258        do ng=1,L_NGAUSS-1
259
260           ! Now compute TAUGAS
261
262           ! Interpolate between water mixing ratios
263           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
264           ! the water data range
265
266           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
267           
268              ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
269              ! the execution time of optci/v -> ~ factor 2 on the whole radiative
270              ! transfer on the tested simulations !
271
272              IF (corrk_recombin) THEN ! Added by JVO
273                tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
274              ELSE
275                tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
276              ENDIF
277             
278              KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
279              KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
280              KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
281              KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
282
283           else
284
285              IF (corrk_recombin) THEN
286                tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
287              ELSE
288                tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
289              ENDIF
290
291              KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) *  &
292                        ( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
293
294              KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) *  &
295                        ( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
296
297              KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) *  &
298                        ( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
299             
300              KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) *  &
301                        ( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
302
303
304           endif
305
306           ! Interpolate the gaseous k-coefficients to the requested T,P values
307
308           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
309                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
310
311           TAUGAS  = U(k)*ANS
312
313           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
314           DTAUKV(K,nw,ng) = TAUGAS &
315                             + DRAYAER & ! DRAYAER includes all scattering contributions
316                             + DCONT ! For parameterized continuum aborption
317
318        end do
319
320        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
321        ! which holds continuum opacity only
322
323        NG              = L_NGAUSS
324        DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
325
326     end do
327  end do
328
329
330  !=======================================================================
331  !     Now the full treatment for the layers, where besides the opacity
332  !     we need to calculate the scattering albedo and asymmetry factors
333
334            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
335            !   but not in the visible
336            !   The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
337            !   This solves random variations of the sw heating at the model top.
338  do iaer=1,naerkind
339    DO NW=1,L_NSPECTV
340      TAUAEROLK(1:4,NW,IAER)=0.d0
341      DO K=5,L_LEVELS
342           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
343      ENDDO
344    ENDDO
345  end do
346
347  DO NW=1,L_NSPECTV
348     DO L=1,L_NLAYRAD-1
349        K              = 2*L+1
350        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))
351        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
352        ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))  ! JVO 2017 : does this 0.999 is really meaningful ?
353        btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
354        COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
355     END DO ! L vertical loop
356     
357     ! Last level
358     L           = L_NLAYRAD
359     K           = 2*L+1
360     atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
361     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
362     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
363     btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
364     COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
365     
366     
367  END DO                    ! NW spectral loop
368
369  DO NG=1,L_NGAUSS
370    DO NW=1,L_NSPECTV
371     DO L=1,L_NLAYRAD-1
372
373        K              = 2*L+1
374        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
375        WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
376
377      END DO ! L vertical loop
378
379        ! Last level
380
381        L              = L_NLAYRAD
382        K              = 2*L+1
383        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
384
385        WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
386
387     END DO                 ! NW spectral loop
388  END DO                    ! NG Gauss loop
389
390  ! Total extinction optical depths
391
392  DO NG=1,L_NGAUSS       ! full gauss loop
393     DO NW=1,L_NSPECTV       
394        TAUCUMV(1,NW,NG)=0.0D0
395        DO K=2,L_LEVELS
396           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
397        END DO
398
399        DO L=1,L_NLAYRAD
400           TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
401        END DO
402        TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
403     END DO           
404  END DO                 ! end full gauss loop
405
406
407
408
409end subroutine optcv
410
411END MODULE optcv_mod
Note: See TracBrowser for help on using the repository browser.