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

Last change on this file since 2057 was 2032, checked in by emillour, 7 years ago

Generic GCM:

  • correct a bug introduced in revision 2026; now that L_NGAUSS is a parameter read in via sugas_corrk (called at first call by callcorrk), automatic arrays of size L_NGAUSS cannot be declared in callcorrk and must be allocated once the value of L_NGAUSS has been set.
  • turned optci, optcv and callcorrk into modules in the process.

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