source: trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90 @ 3595

Last change on this file since 3595 was 3184, checked in by afalco, 13 months ago

Pluto PCM:
Added LMDZ.PLUTO, a copy of the generic model,
cleaned from some unnecessary modules (water, ...)
AF
  • Property svn:executable set to *
File size: 13.2 KB
RevLine 
[3184]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                     igas_CH4, igas_N2
15  use comcstfi_mod, only: g, r, mugaz
16  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
17  use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
18  use tpindex_mod, only: tpindex
19
20  implicit none
21
22  !==================================================================
23  !     
24  !     Purpose
25  !     -------
26  !     Calculates shortwave optical constants at each level.
27  !     
28  !     Authors
29  !     -------
30  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
31  !     
32  !==================================================================
33  !     
34  !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL 
35  !     IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
36  !     LAYER: WBAR, DTAU, COSBAR
37  !     LEVEL: TAU
38  !     
39  !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
40  !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
41  !     
42  !     TLEV(L) - Temperature at the layer boundary
43  !     PLEV(L) - Pressure at the layer boundary (i.e. level)
44  !     GASV(NT,NPS,NW,NG) - Visible k-coefficients
45  !     
46  !-------------------------------------------------------------------
47
48
49  real*8,intent(out) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
50  real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
51  real*8,intent(out) :: TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
52  real*8,intent(out) :: TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
53  real*8,intent(in) :: PLEV(L_LEVELS)
54  real*8,intent(in) :: TMID(L_LEVELS), PMID(L_LEVELS)
55  real*8,intent(out) :: COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
56  real*8,intent(out) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
57
58  ! for aerosols
59  real*8,intent(in) :: QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
60  real*8,intent(in) :: QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
61  real*8,intent(in) :: GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
62  real*8,intent(in) :: TAUAERO(L_LEVELS,NAERKIND)
63 
64  ! local arrays (saved for convenience as need be allocated)
65  real*8,save,allocatable :: TAUAEROLK(:,:,:)
66  real*8,save,allocatable :: TAEROS(:,:,:)
67!$OMP THREADPRIVATE(TAUAEROLK,TAEROS)
68
69  integer L, NW, NG, K, LK, IAER
70  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
71  real*8  ANS, TAUGAS
72  real*8,intent(in) :: TAURAY(L_NSPECTV)
73  real*8  TRAY(L_LEVELS,L_NSPECTV)
74  real*8  DPR(L_LEVELS), U(L_LEVELS)
75  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
76
77  real*8,intent(out) :: taugsurf(L_NSPECTV,L_NGAUSS-1)
78  real*8 DCONT,DAERO
79  real*8 DRAYAER
80  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
81  double precision p_cross
82
83  ! variable species mixing ratio variables
84  real*8,intent(in) :: QVAR(L_LEVELS)
85  real*8,intent(in) :: MUVAR(L_LEVELS)
86  real*8 :: WRATIO(L_LEVELS)
87  real*8  KCOEF(4)
88  integer NVAR(L_LEVELS)
89 
90  ! temporary variables to reduce memory access time to gasv
91  real*8 tmpk(2,2)
92  real*8 tmpkvar(2,2,2)
93
94  ! temporary variables for multiple aerosol calculation
95  real*8 atemp(L_NLAYRAD,L_NSPECTV)
96  real*8 btemp(L_NLAYRAD,L_NSPECTV)
97  real*8 ctemp(L_NLAYRAD,L_NSPECTV)
98
99  ! variables for k in units m^-1
100  real*8 dz(L_LEVELS)
101
102
103  integer igas, jgas
104
105  integer interm
106
107  logical :: firstcall=.true.
108!$OMP THREADPRIVATE(firstcall)
109
110  if (firstcall) then
111    ! allocate local arrays of size "naerkind" (which are also
112    ! "saved" so that this is done only once in for all even if
113    ! we don't need to store the value from a time step to the next)
114    allocate(TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND))
115    allocate(TAEROS(L_LEVELS,L_NSPECTV,NAERKIND))
116    firstcall=.false.
117  endif ! of if (firstcall)
118
119  !! AS: to save time in computing continuum (see bilinearbig)
120  IF (.not.ALLOCATED(indv)) THEN
121      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
122      indv = -9999 ! this initial value means "to be calculated"
123  ENDIF
124
125  !=======================================================================
126  !     Determine the total gas opacity throughout the column, for each
127  !     spectral interval, NW, and each Gauss point, NG.
128  !     Calculate the continuum opacities, i.e., those that do not depend on
129  !     NG, the Gauss index.
130
131  taugsurf(:,:) = 0.0
132  dpr(:)        = 0.0
133  lkcoef(:,:)   = 0.0
134
135  do K=2,L_LEVELS
136     DPR(k) = PLEV(K)-PLEV(K-1)
137
138     ! if we have continuum opacities, we need dz
139     if(kastprof)then
140        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
141        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)
142     else
143        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
144        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
145            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
146     endif
147
148     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
149          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
150
151     do LK=1,4
152        LKCOEF(K,LK) = LCOEF(LK)
153     end do
154  end do                    ! levels
155
156  ! Spectral dependance of aerosol absorption
157            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
158            !   but visible does not handle very well diffusion in first layer.
159            !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
160            !   in the 4 first semilayers in optcv, but not optci.
161            !   This solves random variations of the sw heating at the model top.
162  do iaer=1,naerkind
163     do NW=1,L_NSPECTV
164        TAEROS(1:4,NW,IAER)=0.d0
165        do K=5,L_LEVELS
166           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
167        end do                    ! levels
168     end do
169  end do
170 
171  ! Rayleigh scattering
172  do NW=1,L_NSPECTV
173     TRAY(1:4,NW)   = 1d-30
174     do K=5,L_LEVELS
175        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
176     end do                    ! levels
177  end do
178 
179  !     we ignore K=1...
180  do K=2,L_LEVELS
181
182     do NW=1,L_NSPECTV
183
184        DRAYAER = TRAY(K,NW)
185        !     DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
186        do iaer=1,naerkind
187           DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
188        end do
189
190        DCONT = 0.0 ! continuum absorption
191
192        if(continuum.and.(.not.graybody).and.callgasvis)then
193           ! include continua if necessary
194           wn_cont = dble(wnov(nw))
195           T_cont  = dble(TMID(k))
196           do igas=1,ngasmx
197
198              if(gfrac(igas).eq.-1)then ! variable
199                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
200              else
201                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
202              endif
203
204              dtemp=0.0
205              if(igas.eq.igas_N2)then
206
207                 interm = indv(nw,igas,igas)
208!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
209                 indv(nw,igas,igas) = interm
210                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
211
212               ! elseif(igas.eq.igas_H2)then !AF24: removed
213
214               elseif(igas.eq.igas_CH4)then
215
216                 ! first do self-induced absorption
217                 interm = indv(nw,igas,igas)
218                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
219                 indv(nw,igas,igas) = interm
220
221                 ! then cross-interactions with other gases  !AF24: removed
222                 
223            !   elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then  !AF24: removed
224                 
225              endif
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.