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

Last change on this file since 2520 was 2520, checked in by Guillaume Chaverot, 4 years ago

update of the H2O continuum

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