Context Navigation

optcv.F90

Last change on this file was 3184, checked in by afalco, 10 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

Line
1	MODULE optcv_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE 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.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
141	U(k) = CmkDPR(k)mugaz/muvar(k)
142	else
143	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
144	U(k) = CmkDPR(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)scalepQVAR(k)) ! qvar = mol/mol
200	else
201	p_cont = dble(PMID(k)scalepgfrac(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
385	end subroutine optcv
386
387	END MODULE optcv_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format