source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/optcv.f90 @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 16.6 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!-----------------------------------------------------------------------
38! INCLUDE "comcstfi.h"
39
40      common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
41      common/comcstfi/avocado!,molrad,visc
42     
43      real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
44      real avocado!,molrad,visc
45
46!
47! For Fortran 77/Fortran 90 compliance always use line continuation
48! symbols '&' in columns 73 and 6
49!
50! Group commons according to their type for minimal performance impact
51
52      COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
53     &   , co2cond,callsoil                                             &
54     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
55     &   , callstats,calleofdump                                        &
56     &   , enertest                                                     &
57     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
58     &   , radfixed                                                     &
59     &   , meanOLR, specOLR                                             &
60     &   , kastprof                                                     &
61     &   , nosurf, oblate                                               &     
62     &   , newtonian, testradtimes                                      &
63     &   , check_cpp_match, force_cpp                                   &
64     &   , rayleigh                                                     &
65     &   , stelbbody                                                    &
66     &   , nearco2cond                                                  &
67     &   , tracer, mass_redistrib, varactive, varfixed                  &
68     &   , sedimentation,water,watercond,waterrain                      &
69     &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                       &
70     &   , aerofixco2,aerofixh2o                                        &
71     &   , hydrology, sourceevol                                        &
72     &   , CLFvarying                                                   &
73     &   , strictboundcorrk                                             &                                       
74     &   , ok_slab_ocean                                                &
75     &   , ok_slab_sic                                                  &
76     &   , ok_slab_heat_transp                                         
77
78
79      COMMON/callkeys_i/iaervar,iddist,iradia,startype
80     
81      COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb,     &
82     &                  Tstrat,tplanet,obs_tau_col_tropo,               &
83     &                  obs_tau_col_strato,pres_bottom_tropo,           &
84     &                  pres_top_tropo,pres_bottom_strato,              &
85     &                  pres_top_strato,size_tropo,size_strato,satval,  &
86     &                  CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
87     &                  maxicethick,Tsaldiff,tau_relax,cloudlvl,        &
88     &                  icetstep,intheat,flatten,Rmean,J2,MassPlanet
89     
90      logical callrad,corrk,calldifv,UseTurbDiff                        &
91     &   , calladj,co2cond,callsoil                                     &
92     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
93     &   , callstats,calleofdump                                        &
94     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
95     &   , strictboundcorrk                                             
96
97      logical enertest
98      logical nonideal
99      logical meanOLR
100      logical specOLR
101      logical kastprof
102      logical newtonian
103      logical check_cpp_match
104      logical force_cpp
105      logical testradtimes
106      logical rayleigh
107      logical stelbbody
108      logical ozone
109      logical nearco2cond
110      logical tracer
111      logical mass_redistrib
112      logical varactive
113      logical varfixed
114      logical radfixed
115      logical sedimentation
116      logical water,watercond,waterrain
117      logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
118      logical aerofixco2,aerofixh2o
119      logical hydrology
120      logical sourceevol
121      logical CLFvarying
122      logical nosurf
123      logical oblate
124      logical ok_slab_ocean
125      logical ok_slab_sic
126      logical ok_slab_heat_transp
127
128      integer iddist
129      integer iaervar
130      integer iradia
131      integer startype
132
133      real topdustref
134      real Nmix_co2
135      real dusttau
136      real Fat1AU
137      real stelTbb
138      real Tstrat
139      real tplanet
140      real obs_tau_col_tropo
141      real obs_tau_col_strato
142      real pres_bottom_tropo
143      real pres_top_tropo
144      real pres_bottom_strato
145      real pres_top_strato
146      real size_tropo
147      real size_strato
148      real satval
149      real CLFfixval
150      real n2mixratio
151      real co2supsat
152      real pceil
153      real albedosnow
154      real maxicethick
155      real Tsaldiff
156      real tau_relax
157      real cloudlvl
158      real icetstep
159      real intheat
160      real flatten
161      real Rmean
162      real J2
163      real MassPlanet
164
165
166  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
167  real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
168  real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
169  real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
170  real*8 PLEV(L_LEVELS)
171  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
172  real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
173  real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
174
175  ! for aerosols
176  real*8  QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
177  real*8  QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
178  real*8  GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
179  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
180  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
181  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
182
183  integer L, NW, NG, K, LK, IAER
184  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
185  real*8  ANS, TAUGAS
186  real*8  TAURAY(L_NSPECTV)
187  real*8  TRAY(L_LEVELS,L_NSPECTV)
188  real*8  TRAYAER
189  real*8  DPR(L_LEVELS), U(L_LEVELS)
190  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
191
192  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
193  real*8 DCONT,DAERO
194  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
195  double precision p_cross
196
197  ! variable species mixing ratio variables
198  real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
199  real*8  KCOEF(4)
200  integer NVAR(L_LEVELS)
201
202  ! temporary variables for multiple aerosol calculation
203  real*8 atemp(L_NLAYRAD,L_NSPECTV)
204  real*8 btemp(L_NLAYRAD,L_NSPECTV)
205  real*8 ctemp(L_NLAYRAD,L_NSPECTV)
206
207  ! variables for k in units m^-1
208  real*8 dz(L_LEVELS)
209
210  integer igas, jgas
211
212  integer interm
213
214  !! AS: to save time in computing continuum (see bilinearbig)
215  IF (.not.ALLOCATED(indv)) THEN
216      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
217      indv = -9999 ! this initial value means "to be calculated"
218  ENDIF
219
220  !=======================================================================
221  !     Determine the total gas opacity throughout the column, for each
222  !     spectral interval, NW, and each Gauss point, NG.
223  !     Calculate the continuum opacities, i.e., those that do not depend on
224  !     NG, the Gauss index.
225
226  taugsurf(:,:) = 0.0
227  dpr(:)        = 0.0
228  lkcoef(:,:)   = 0.0
229
230  do K=2,L_LEVELS
231     DPR(k) = PLEV(K)-PLEV(K-1)
232
233     ! if we have continuum opacities, we need dz
234     if(kastprof)then
235        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
236        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
237     else
238        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
239        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
240            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
241     endif
242
243     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
244          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
245
246     do LK=1,4
247        LKCOEF(K,LK) = LCOEF(LK)
248     end do
249  end do                    ! levels
250
251
252  do iaer=1,naerkind
253     do NW=1,L_NSPECTV
254        do K=2,L_LEVELS
255           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
256        end do                    ! levels
257     end do
258  end do
259  do NW=1,L_NSPECTV
260     do K=2,L_LEVELS
261        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
262     end do                    ! levels
263  end do
264 
265  !     we ignore K=1...
266  do K=2,L_LEVELS
267
268     do NW=1,L_NSPECTV
269
270        TRAYAER = TRAY(K,NW)
271        !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
272        do iaer=1,naerkind
273           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
274        end do
275
276        DCONT = 0.0 ! continuum absorption
277
278        if(continuum.and.(.not.graybody).and.callgasvis)then
279           ! include continua if necessary
280           wn_cont = dble(wnov(nw))
281           T_cont  = dble(TMID(k))
282           do igas=1,ngasmx
283
284              if(gfrac(igas).eq.-1)then ! variable
285                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
286              else
287                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
288              endif
289
290              dtemp=0.0
291              if(igas.eq.igas_N2)then
292
293                 interm = indv(nw,igas,igas)
294!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
295                 indv(nw,igas,igas) = interm
296                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
297
298              elseif(igas.eq.igas_H2)then
299
300                 ! first do self-induced absorption
301                 interm = indv(nw,igas,igas)
302                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
303                 indv(nw,igas,igas) = interm
304
305                 ! then cross-interactions with other gases
306                 do jgas=1,ngasmx
307                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
308                    dtempc  = 0.0
309                    if(jgas.eq.igas_N2)then
310                       interm = indv(nw,igas,jgas)
311                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
312                       indv(nw,igas,jgas) = interm
313                       ! should be irrelevant in the visible
314                    elseif(jgas.eq.igas_He)then
315                       interm = indv(nw,igas,jgas)
316                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
317                       indv(nw,igas,jgas) = interm
318                    endif
319                    dtemp = dtemp + dtempc
320                 enddo
321
322              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
323
324                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
325                 if(H2Ocont_simple)then
326                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
327                 else
328                    interm = indv(nw,igas,igas)
329                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
330                    indv(nw,igas,igas) = interm
331                 endif
332
333              endif
334
335              DCONT = DCONT + dtemp
336
337           enddo
338
339           DCONT = DCONT*dz(k)
340
341        endif
342
343        do ng=1,L_NGAUSS-1
344
345           ! Now compute TAUGAS
346
347           ! Interpolate between water mixing ratios
348           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
349           ! the water data range
350
351           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
352              KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
353              KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
354              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
355              KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
356           else
357
358              KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
359                   (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
360                   GASV(MT(K),MP(K),NVAR(K),NW,NG))
361
362              KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
363                   (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
364                   GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
365
366              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
367                   (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
368                   GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
369
370              KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
371                   (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
372                   GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
373
374           endif
375
376           ! Interpolate the gaseous k-coefficients to the requested T,P values
377
378           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
379                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
380
381           TAUGAS  = U(k)*ANS
382
383           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
384           DTAUKV(K,nw,ng) = TAUGAS & 
385                             + TRAYAER & ! TRAYAER includes all scattering contributions
386                             + DCONT ! For parameterized continuum aborption
387
388        end do
389
390        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
391        ! which holds continuum opacity only
392
393        NG              = L_NGAUSS
394        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
395
396        do iaer=1,naerkind
397           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) +  TAEROS(K,NW,IAER)
398        end do ! a bug was here!
399
400     end do
401  end do
402
403
404  !=======================================================================
405  !     Now the full treatment for the layers, where besides the opacity
406  !     we need to calculate the scattering albedo and asymmetry factors
407
408  do iaer=1,naerkind
409    DO NW=1,L_NSPECTV
410      DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
411           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
412      ENDDO
413    ENDDO
414  end do
415
416  DO NW=1,L_NSPECTV
417     DO L=1,L_NLAYRAD-1
418        K              = 2*L+1
419        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))
420        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
421        ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
422        btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
423        COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
424     END DO ! L vertical loop
425     
426     !last level
427     L              = L_NLAYRAD
428     K              = 2*L+1
429     atemp(L,NW)    = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
430     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
431     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
432     btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
433     COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
434     
435     
436  END DO                    ! NW spectral loop
437
438  DO NG=1,L_NGAUSS
439    DO NW=1,L_NSPECTV
440     DO L=1,L_NLAYRAD-1
441
442        K              = 2*L+1
443        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
444        WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
445
446      END DO ! L vertical loop
447
448        !     No vertical averaging on bottom layer
449
450        L              = L_NLAYRAD
451        K              = 2*L+1
452        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
453
454        WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
455     END DO                 ! NW spectral loop
456  END DO                    ! NG Gauss loop
457
458  ! Total extinction optical depths
459
460  DO NG=1,L_NGAUSS       ! full gauss loop
461     DO NW=1,L_NSPECTV       
462        TAUV(1,NW,NG)=0.0D0
463        DO L=1,L_NLAYRAD
464           TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
465        END DO
466
467        TAUCUMV(1,NW,NG)=0.0D0
468        DO K=2,L_LEVELS
469           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
470        END DO
471     END DO           
472  END DO                 ! end full gauss loop
473
474
475  return
476
477
478end subroutine optcv
Note: See TracBrowser for help on using the repository browser.