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

Last change on this file since 2546 was 2543, checked in by aslmd, 4 years ago

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ

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