Context Navigation

optcv.F90 @ 3728

Last change on this file since 3728 was 3696, checked in by emillour, 3 months ago
Generic PCM: Bug fix in the computation of rho() in calc_rayleigh Made calc_rayleigh.F90 a module and commented out some prints which swamped the standard output. GM+EM
Property svn:executable set to ``*
File size: 20.4 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	TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
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_CO2, igas_O2
15	use comcstfi_mod, only: g, r, mugaz
16	use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, &
17	generic_continuum_database,rayleigh
18	use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
19	use tpindex_mod, only: tpindex
20	use interpolate_continuum_mod, only: interpolate_continuum
21	use calc_rayleigh_mod, only: calc_rayleigh
22
23	implicit none
24
25	!==================================================================
26	!
27	! Purpose
28	! -------
29	! Calculates shortwave optical constants at each level.
30	!
31	! Authors
32	! -------
33	! Adapted from the NASA Ames code by R. Wordsworth (2009)
34	!
35	!==================================================================
36	!
37	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
38	! IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
39	! LAYER: WBAR, DTAU, COSBAR
40	! LEVEL: TAU
41	!
42	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
43	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
44	!
45	! TLEV(L) - Temperature at the layer boundary
46	! PLEV(L) - Pressure at the layer boundary (i.e. level)
47	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
48	!
49	!-------------------------------------------------------------------
50
51
52	real*8,intent(out) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
53	real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
54	real*8,intent(out) :: TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
55	real*8,intent(out) :: TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
56	real*8,intent(in) :: PLEV(L_LEVELS)
57	real*8,intent(in) :: TMID(L_LEVELS), PMID(L_LEVELS)
58	real*8,intent(out) :: COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
59	real*8,intent(out) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
60
61	! for aerosols
62	real*8,intent(in) :: QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
63	real*8,intent(in) :: QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
64	real*8,intent(in) :: GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
65	real*8,intent(in) :: TAUAERO(L_LEVELS,NAERKIND)
66
67	! local arrays (saved for convenience as need be allocated)
68	real*8,save,allocatable :: TAUAEROLK(:,:,:)
69	real*8,save,allocatable :: TAEROS(:,:,:)
70	!$OMP THREADPRIVATE(TAUAEROLK,TAEROS)
71
72	integer L, NW, NG, K, LK, IAER
73	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
74	real*8 ANS, TAUGAS
75	real*8 TAURAY(L_LEVELS,L_NSPECTV)
76	real*8 TRAY(L_LEVELS,L_NSPECTV)
77	real*8 DPR(L_LEVELS), U(L_LEVELS)
78	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
79
80	real*8,intent(out) :: taugsurf(L_NSPECTV,L_NGAUSS-1)
81	real*8 DCONT,DAERO
82	real*8 DRAYAER
83	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
84	double precision p_cross
85
86	! variable species mixing ratio variables
87	real*8,intent(in) :: QVAR(L_LEVELS)
88	real*8,intent(in) :: MUVAR(L_LEVELS)
89	real*8,intent(in) :: FRACVAR(ngasmx,L_LEVELS)
90	real*8 :: WRATIO(L_LEVELS)
91	real*8 KCOEF(4)
92	integer NVAR(L_LEVELS)
93
94	! temporary variables to reduce memory access time to gasv
95	real*8 tmpk(2,2)
96	real*8 tmpkvar(2,2,2)
97
98	! temporary variables for multiple aerosol calculation
99	real*8 atemp(L_NLAYRAD,L_NSPECTV)
100	real*8 btemp(L_NLAYRAD,L_NSPECTV)
101	real*8 ctemp(L_NLAYRAD,L_NSPECTV)
102
103	! variables for k in units m^-1
104	real*8 dz(L_LEVELS)
105
106
107	integer igas, jgas
108
109	integer interm
110
111	logical :: firstcall=.true.
112	!$OMP THREADPRIVATE(firstcall)
113
114	if (firstcall) then
115	! allocate local arrays of size "naerkind" (which are also
116	! "saved" so that this is done only once in for all even if
117	! we don't need to store the value from a time step to the next)
118	allocate(TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND))
119	allocate(TAEROS(L_LEVELS,L_NSPECTV,NAERKIND))
120	firstcall=.false.
121	endif ! of if (firstcall)
122
123	!! AS: to save time in computing continuum (see bilinearbig)
124	IF (.not.ALLOCATED(indv)) THEN
125	ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
126	indv = -9999 ! this initial value means "to be calculated"
127	ENDIF
128
129	!=======================================================================
130	! Determine the total gas opacity throughout the column, for each
131	! spectral interval, NW, and each Gauss point, NG.
132	! Calculate the continuum opacities, i.e., those that do not depend on
133	! NG, the Gauss index.
134
135	taugsurf(:,:) = 0.0
136	dpr(:) = 0.0
137	lkcoef(:,:) = 0.0
138
139	do K=2,L_LEVELS
140	DPR(k) = PLEV(K)-PLEV(K-1)
141
142	! if we have continuum opacities, we need dz
143	if(kastprof)then
144	dz(k) = dpr(k)(1000.0d08.314463d0/muvar(k))TMID(K)/(gPMID(K))
145	U(k) = CmkDPR(k)mugaz/muvar(k)
146	else
147	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
148	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
149	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
150	endif
151
152	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
153	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
154
155	do LK=1,4
156	LKCOEF(K,LK) = LCOEF(LK)
157	end do
158	end do ! levels
159
160	! Spectral dependance of aerosol absorption
161	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
162	! but visible does not handle very well diffusion in first layer.
163	! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
164	! in the 4 first semilayers in optcv, but not optci.
165	! This solves random variations of the sw heating at the model top.
166	do iaer=1,naerkind
167	do NW=1,L_NSPECTV
168	TAEROS(1:4,NW,IAER)=0.d0
169	do K=5,L_LEVELS
170	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
171	end do ! levels
172	end do
173	end do
174
175	!=======================================================================
176	! Set up the wavelength independent part of the Rayleigh scattering.
177	! WAVEV is in microns. There is no Rayleigh scattering in the IR.
178
179	if(rayleigh) then
180	call calc_rayleigh(QVAR,MUVAR,PMID,TMID,TAURAY)
181	else
182	print*,'setspv: No Rayleigh scattering, check for NaN in output!'
183	do NW=1,L_NSPECTV
184	TAURAY(:,NW) = 1E-16
185	end do
186	endif
187
188	! Computation of pressure dependant part of Rayleigh scattering
189	do NW=1,L_NSPECTV
190	TRAY(1:4,NW) = 1d-30
191	do K=5,L_LEVELS
192	TRAY(K,NW) = TAURAY(K,NW) * DPR(K)
193	end do ! levels
194	end do
195
196	! we ignore K=1...
197	do K=2,L_LEVELS
198
199	do NW=1,L_NSPECTV
200
201	DRAYAER = TRAY(K,NW)
202	! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
203	do iaer=1,naerkind
204	DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
205	end do
206
207	DCONT = 0.0 ! continuum absorption
208
209	if(continuum.and.(.not.graybody).and.callgasvis)then
210	! include continua if necessary
211
212	if(generic_continuum_database)then
213	T_cont = dble(TMID(k))
214	do igas=1,ngasmx
215
216	if(gfrac(igas).eq.-1)then ! variable
217	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
218	elseif(varspec) then
219	p_cont = dble(PMID(k)scalepFRACVAR(igas,k)*(1.-QVAR(k)))
220	else
221	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
222	endif
223
224	dtemp=0.0
225
226	if (igas .eq. igas_N2) then
227	call interpolate_continuum('',igas_N2,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
228	do jgas=1,ngasmx
229	if (jgas .eq. igas_H2) then
230	call interpolate_continuum('',igas_N2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
231	elseif (jgas .eq. igas_O2) then
232	call interpolate_continuum('',igas_N2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
233	elseif (jgas .eq. igas_CH4) then
234	call interpolate_continuum('',igas_N2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
235	endif
236	enddo
237	elseif (igas .eq. igas_O2) then
238	call interpolate_continuum('',igas_O2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
239	do jgas=1,ngasmx
240	if (jgas .eq. igas_CO2) then
241	call interpolate_continuum('',igas_CO2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
242	endif
243	enddo
244	elseif (igas .eq. igas_H2) then
245	call interpolate_continuum('',igas_H2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
246	do jgas=1,ngasmx
247	if (jgas .eq. igas_CH4) then
248	call interpolate_continuum('',igas_H2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
249	elseif (jgas .eq. igas_He) then
250	call interpolate_continuum('',igas_H2,igas_He,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
251	endif
252	enddo
253	elseif (igas .eq. igas_CH4) then
254	call interpolate_continuum('',igas_CH4,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
255	elseif (igas .eq. igas_H2O) then
256	call interpolate_continuum('',igas_H2O,igas_H2O,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
257	do jgas=1,ngasmx
258	if (jgas .eq. igas_N2) then
259	call interpolate_continuum('',igas_H2O,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
260	elseif (jgas .eq. igas_O2) then
261	call interpolate_continuum('',igas_H2O,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
262	elseif (jgas .eq. igas_CO2) then
263	call interpolate_continuum('',igas_H2O,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
264	endif
265	enddo
266	elseif (igas .eq. igas_CO2) then
267	call interpolate_continuum('',igas_CO2,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
268	do jgas=1,ngasmx
269	if (jgas .eq. igas_H2) then
270	call interpolate_continuum('',igas_CO2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
271	elseif (jgas .eq. igas_CH4) then
272	call interpolate_continuum('',igas_CO2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
273	endif
274	enddo
275	endif
276
277	DCONT = DCONT + dtemp
278
279	enddo ! igas=1,ngasmx
280
281	else ! generic_continuum_database
282
283
284	wn_cont = dble(wnov(nw))
285	T_cont = dble(TMID(k))
286	do igas=1,ngasmx
287
288	if(gfrac(igas).eq.-1)then ! variable
289	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
290	elseif(varspec) then
291	p_cont = dble(PMID(k)scalepFRACVAR(igas,k)*(1.-QVAR(k)))
292	else
293	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
294	endif
295
296	dtemp=0.0
297	if(igas.eq.igas_N2)then
298
299	interm = indv(nw,igas,igas)
300	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
301	indv(nw,igas,igas) = interm
302	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
303
304	elseif(igas.eq.igas_H2)then
305
306	! first do self-induced absorption
307	interm = indv(nw,igas,igas)
308	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
309	indv(nw,igas,igas) = interm
310
311	! then cross-interactions with other gases
312	do jgas=1,ngasmx
313	if(varspec) then
314	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
315	else
316	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
317	endif
318	dtempc = 0.0
319	if(jgas.eq.igas_N2)then
320	interm = indv(nw,igas,jgas)
321	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
322	indv(nw,igas,jgas) = interm
323	! should be irrelevant in the visible
324	elseif(jgas.eq.igas_CO2)then
325	interm = indv(nw,igas,jgas)
326	call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
327	indv(nw,igas,jgas) = interm
328	! might not be relevant in the visible
329	elseif(jgas.eq.igas_He)then
330	interm = indv(nw,igas,jgas)
331	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
332	indv(nw,igas,jgas) = interm
333	endif
334	dtemp = dtemp + dtempc
335	enddo
336
337	elseif(igas.eq.igas_CH4)then
338
339	! first do self-induced absorption
340	interm = indv(nw,igas,igas)
341	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
342	indv(nw,igas,igas) = interm
343
344	! then cross-interactions with other gases
345	do jgas=1,ngasmx
346	if(varspec) then
347	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
348	else
349	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
350	endif
351	dtempc = 0.0d0
352	if(jgas.eq.igas_H2)then
353	interm = indv(nw,igas,jgas)
354	call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
355	indv(nw,igas,jgas) = interm
356	elseif(jgas.eq.igas_CO2)then
357	interm = indv(nw,igas,jgas)
358	call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
359	indv(nw,igas,jgas) = interm
360	! might not be relevant in the visible
361	elseif(jgas.eq.igas_He)then
362	interm = indv(nw,igas,jgas)
363	call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
364	indv(nw,igas,jgas) = interm
365	endif
366	dtemp = dtemp + dtempc
367	enddo
368
369	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
370	! Compute self and foreign (with air) continuum of H2O
371	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
372	interm = indv(nw,igas,igas)
373	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
374	indv(nw,igas,igas) = interm
375
376	endif
377
378	DCONT = DCONT + dtemp
379
380	enddo
381
382	DCONT = DCONT*dz(k)
383
384	endif ! generic_continuum_database
385	endif ! continuum
386
387	do ng=1,L_NGAUSS-1
388
389	! Now compute TAUGAS
390
391	! Interpolate between water mixing ratios
392	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
393	! the water data range
394
395	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
396
397	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
398	! the execution time of optci/v -> ~ factor 2 on the whole radiative
399	! transfer on the tested simulations !
400
401	IF (corrk_recombin) THEN ! Added by JVO
402	tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
403	ELSE
404	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
405	ENDIF
406
407	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
408	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
409	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
410	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
411
412	else
413
414	IF (corrk_recombin) THEN
415	tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
416	ELSE
417	tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
418	ENDIF
419
420	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
421	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
422
423	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
424	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
425
426	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
427	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
428
429	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
430	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
431
432
433	endif
434
435	! Interpolate the gaseous k-coefficients to the requested T,P values
436
437	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
438	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
439
440	TAUGAS = U(k)*ANS
441
442	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
443	DTAUKV(K,nw,ng) = TAUGAS &
444	+ DRAYAER & ! DRAYAER includes all scattering contributions
445	+ DCONT ! For parameterized continuum aborption
446
447	end do
448
449	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
450	! which holds continuum opacity only
451
452	NG = L_NGAUSS
453	DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
454
455	end do
456	end do
457
458
459	!=======================================================================
460	! Now the full treatment for the layers, where besides the opacity
461	! we need to calculate the scattering albedo and asymmetry factors
462
463	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
464	! but not in the visible
465	! The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
466	! This solves random variations of the sw heating at the model top.
467	do iaer=1,naerkind
468	DO NW=1,L_NSPECTV
469	TAUAEROLK(1:4,NW,IAER)=0.d0
470	DO K=5,L_LEVELS
471	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
472	ENDDO
473	ENDDO
474	end do
475
476	DO NW=1,L_NSPECTV
477	DO L=1,L_NLAYRAD-1
478	K = 2*L+1
479	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))
480	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
481	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
482	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
483	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
484	END DO ! L vertical loop
485
486	! Last level
487	L = L_NLAYRAD
488	K = 2*L+1
489	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
490	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
491	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
492	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
493	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
494
495
496	END DO ! NW spectral loop
497
498	DO NG=1,L_NGAUSS
499	DO NW=1,L_NSPECTV
500	DO L=1,L_NLAYRAD-1
501
502	K = 2*L+1
503	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
504	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
505
506	END DO ! L vertical loop
507
508	! Last level
509
510	L = L_NLAYRAD
511	K = 2*L+1
512	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
513
514	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
515
516	END DO ! NW spectral loop
517	END DO ! NG Gauss loop
518
519	! Total extinction optical depths
520
521	DO NG=1,L_NGAUSS ! full gauss loop
522	DO NW=1,L_NSPECTV
523	TAUCUMV(1,NW,NG)=0.0D0
524	DO K=2,L_LEVELS
525	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
526	END DO
527
528	DO L=1,L_NLAYRAD
529	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
530	END DO
531	TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
532	END DO
533	END DO ! end full gauss loop
534
535
536
537
538	end subroutine optcv
539
540	END MODULE optcv_mod

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

Download in other formats:

Original Format