Context Navigation

callcorrk.F90 @ 2245

Last change on this file since 2245 was 2138, checked in by jvatant, 7 years ago
Correct outptuts of dtaui/v with correct ponderations of weights, but now it is attenuation exp(-tau) --JVO
Property svn:executable set to ``*
File size: 26.5 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
3	tsurf,fract,dist_star, &
4	dtlw,dtsw,fluxsurf_lw, &
5	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
6	fluxabs_sw,fluxtop_dn, &
7	OLR_nu,OSR_nu, &
8	int_dtaui,int_dtauv, &
9	lastcall)
10
11	use mod_phys_lmdz_para, only : is_master
12	use radinc_h
13	use radcommon_h
14	use gases_h
15	USE tracer_h
16	use callkeys_mod, only: global1d, szangle
17	use comcstfi_mod, only: pi, mugaz, cpp
18	use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin, &
19	strictboundcorrk,specOLR,diagdtau
20	use geometry_mod, only: latitude
21
22	implicit none
23
24	!==================================================================
25	!
26	! Purpose
27	! -------
28	! Solve the radiative transfer using the correlated-k method for
29	! the gaseous absorption and the Toon et al. (1989) method for
30	! scatttering due to aerosols.
31	!
32	! Authors
33	! -------
34	! Emmanuel 01/2001, Forget 09/2001
35	! Robin Wordsworth (2009)
36	! Jan Vatant d'Ollone (2018) -> corrk recombining case
37	!
38	!==================================================================
39
40	!-----------------------------------------------------------------------
41	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
42	! Layer #1 is the layer near the ground.
43	! Layer #nlayer is the layer at the top.
44	!-----------------------------------------------------------------------
45
46
47	! INPUT
48	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
49	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
50	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/kg).
51	INTEGER,INTENT(IN) :: nq ! Number of tracers.
52	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
53	REAL,INTENT(IN) :: zday ! Time elapsed since Ls=0 (sols).
54	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
55	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
56	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
57	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
58	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
59	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
60	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
61	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
62	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
63	logical,intent(in) :: lastcall ! Signals last call to physics.
64
65	! OUTPUT
66	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
67	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
68	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
69	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
70	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
71	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
72	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
73	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
74	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
75	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
76	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
77	REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
78	REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
79
80
81	!-----------------------------------------------------------------------
82	! Declaration of the variables required by correlated-k subroutines
83	! Numbered from top to bottom (unlike in the GCM)
84	!-----------------------------------------------------------------------
85
86	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
87	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
88
89	! Optical values for the optci/cv subroutines
90	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
91	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
92	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
93	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
94	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
95	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
96	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
97	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
98	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
99	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
100
101	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
102	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
103	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
104	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
105	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
106	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
107	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
108	REAL*8 albi,acosz
109	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
110
111	INTEGER ig,l,k,nw,iq,ip,ilay,it
112
113	LOGICAL found
114
115	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
116	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
117
118	logical OLRz
119	real*8 NFLUXGNDV_nu(L_NSPECTV)
120
121	! Included by MT for albedo calculations.
122	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
123	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
124
125	! For variable haze
126	REAL*8 seashazefact(L_LEVELS)
127
128	! For muphys optics
129	REAL*8 pqmo(ngrid,nlayer,nmicro) ! Tracers for microphysics optics (X/m2).
130	REAL*8 i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics )
131
132	! For corr-k recombining
133	REAL*8 pqr(ngrid,L_PINT,L_REFVAR) ! Tracers for corr-k recombining (mol/mol).
134	REAL*8 fact, tmin, tmax
135
136	LOGICAL usept(L_PINT,L_NTREF) ! mask if pfref grid point will be used
137	INTEGER inflay(L_PINT) ! nearest inferior GCM layer for pfgasref grid points
138
139	!=======================================================================
140	! I. Initialization on every call
141	!=======================================================================
142
143
144	! How much light do we get ?
145	do nw=1,L_NSPECTV
146	stel(nw)=stellarf(nw)/(dist_star**2)
147	end do
148
149	! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2
150	! NOTE: it should be moved somewhere else: calmufi performs the same kind of
151	! computations... waste of time...
152	i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
153	pqmo(:,:,:) = 0.0
154	DO iq=1,nmicro
155	pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:)
156	ENDDO
157
158	! Default value for fixed species for whom vmr has been taken
159	! into account while computing high-resolution spectra
160	if (corrk_recombin) pqr(:,:,:) = 1.0
161
162	!-----------------------------------------------------------------------
163	do ig=1,ngrid ! Starting Big Loop over every GCM column
164	!-----------------------------------------------------------------------
165
166	! Recombine reference corr-k if needed
167	if (corrk_recombin) then
168
169	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
170	! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values
171	! who match the atmospheric conditions ) which is then processed as a standard pre-mix in
172	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
173
174	! Extract tracers for variable radiative species
175	! Also find the nearest GCM layer under each ref pressure
176	do ip=1,L_PINT
177
178	ilay=0
179	found = .false.
180	do l=1,nlayer
181	if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar])
182	found=.true.
183	ilay=l
184	endif
185	enddo
186
187	if (.not. found ) then ! set to min
188	do iq=1,L_REFVAR
189	if ( radvar_mask(iq) ) then
190	pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
191	endif
192	enddo
193	else
194	if (ilay==nlayer) then ! set to max
195	do iq=1,L_REFVAR
196	if ( radvar_mask(iq) ) then
197	pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
198	endif
199	enddo
200	else ! standard
201	fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(1,ilay+1) ) / ( pplay(1,ilay) - pplay(1,ilay+1) ) ! pfgasref=log(p[mbar])
202	do iq=1,L_REFVAR
203	if ( radvar_mask(iq) ) then
204	pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))*fact pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact)
205	pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
206	endif
207	enddo
208	endif ! if ilay==nlayer
209	endif ! if not found
210
211	inflay(ip) = ilay
212
213	enddo ! ip=1,L_PINT
214
215	! NB : The following usept is a trick to call recombine only for the reference T-P
216	! grid points that are useful given the temperature range at this altitude
217	! It saves a looot of time - JVO 18
218	usept(:,:) = .true.
219	do ip=1,L_PINT-1
220	if ( inflay(ip+1)==nlayer ) then
221	usept(ip,:) = .false.
222	endif
223	if ( inflay(ip) == 0 ) then
224	usept(ip+1:,:) = .false.
225	endif
226	if ( usept(ip,1) ) then ! if not all false
227	tmin = minval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
228	tmax = maxval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
229	do it=1,L_NTREF-1
230	if ( tgasref(it+1) .lt. tmin ) then
231	usept(ip,it) = .false.
232	endif
233	enddo
234	do it=2,L_NTREF
235	if ( tgasref(it-1) .gt. tmax ) then
236	usept(ip,it) = .false.
237	endif
238	enddo
239	! in case of out-of-bounds
240	if ( tgasref(1) .lt. tmin ) usept(ip,1) = .true.
241	if ( tgasref(L_NTREF) .gt. tmax ) usept(ip,L_NTREF) = .true.
242	endif
243	enddo ! ip=1,L_PINT-1
244	! deal with last bound
245	if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:)
246
247
248	do ip=1,L_PINT
249
250	! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied
251	do it=1,L_NTREF
252
253	if ( usept(ip,it) .eqv. .false. ) cycle
254
255	do l=1,L_REFVAR
256	if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01 & ! +- 1%
257	.or. ( useptold(ip,it) .eqv. .false. ) ) then ! in case T change but not the tracers
258	call recombin_corrk( pqr(ig,ip,:),ip,it )
259	exit ! one is enough to trigger the update
260	endif
261	enddo
262
263	enddo
264
265	enddo ! ip=1,L_PINT
266
267	useptold(:,:)=usept(:,:)
268
269	endif ! if corrk_recombin
270
271	!=======================================================================
272	! II. Transformation of the GCM variables
273	!=======================================================================
274
275
276	! Albedo and Emissivity.
277	albi=1-emis(ig) ! Long Wave.
278	DO nw=1,L_NSPECTV ! Short Wave loop.
279	albv(nw)=albedo(ig,nw)
280	ENDDO
281
282	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
283	acosz = cos(pi*szangle/180.0)
284	print*,'acosz=',acosz,', szangle=',szangle
285	else
286	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
287	endif
288
289	!-----------------------------------------------------------------------
290	! Pressure and temperature
291	!-----------------------------------------------------------------------
292
293	DO l=1,nlayer
294	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
295	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
296	tlevrad(2*l) = pt(ig,nlayer+1-l)
297	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
298	END DO
299
300	plevrad(1) = 0.
301	plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
302
303	tlevrad(1) = tlevrad(2)
304	tlevrad(2*nlayer+1)=tsurf(ig)
305
306	pmid(1) = max(pgasmin,0.0001*plevrad(3))
307	pmid(2) = pmid(1)
308
309	tmid(1) = tlevrad(2)
310	tmid(2) = tmid(1)
311
312	DO l=1,L_NLAYRAD-1
313	tmid(2l+1) = tlevrad(2l+1)
314	tmid(2l+2) = tlevrad(2l+1)
315	pmid(2l+1) = plevrad(2l+1)
316	pmid(2l+2) = plevrad(2l+1)
317	END DO
318	pmid(L_LEVELS) = plevrad(L_LEVELS)
319	tmid(L_LEVELS) = tlevrad(L_LEVELS)
320
321	!!Alternative interpolation:
322	! pmid(3) = pmid(1)
323	! pmid(4) = pmid(1)
324	! tmid(3) = tmid(1)
325	! tmid(4) = tmid(1)
326	! DO l=2,L_NLAYRAD-1
327	! tmid(2l+1) = tlevrad(2l)
328	! tmid(2l+2) = tlevrad(2l)
329	! pmid(2l+1) = plevrad(2l)
330	! pmid(2l+2) = plevrad(2l)
331	! END DO
332	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
333	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
334
335	! Test for out-of-bounds pressure.
336	if(plevrad(3).lt.pgasmin)then
337	print*,'Minimum pressure is outside the radiative'
338	print*,'transfer kmatrix bounds, exiting.'
339	call abort
340	elseif(plevrad(L_LEVELS).gt.pgasmax)then
341	print*,'Maximum pressure is outside the radiative'
342	print*,'transfer kmatrix bounds, exiting.'
343	call abort
344	endif
345
346	! Test for out-of-bounds temperature.
347	do k=1,L_LEVELS
348	if(tlevrad(k).lt.tgasmin)then
349	print*,'Minimum temperature is outside the radiative'
350	print*,'transfer kmatrix bounds'
351	print*,"k=",k," tlevrad(k)=",tlevrad(k)
352	print*,"tgasmin=",tgasmin
353	if (strictboundcorrk) then
354	call abort
355	else
356	print,'**********************************************'
357	print*,'we allow model to continue with tlevrad=tgasmin'
358	print*,' ... we assume we know what you are doing ... '
359	print*,' ... but do not let this happen too often ... '
360	print,'**********************************************'
361	!tlevrad(k)=tgasmin
362	endif
363	elseif(tlevrad(k).gt.tgasmax)then
364	! print*,'Maximum temperature is outside the radiative'
365	! print*,'transfer kmatrix bounds, exiting.'
366	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
367	! print*,"tgasmax=",tgasmax
368	if (strictboundcorrk) then
369	call abort
370	else
371	! print,'**********************************************'
372	! print*,'we allow model to continue with tlevrad=tgasmax'
373	! print*,' ... we assume we know what you are doing ... '
374	! print*,' ... but do not let this happen too often ... '
375	! print,'**********************************************'
376	!tlevrad(k)=tgasmax
377	endif
378	endif
379	enddo
380	do k=1,L_NLAYRAD+1
381	if(tmid(k).lt.tgasmin)then
382	print*,'Minimum temperature is outside the radiative'
383	print*,'transfer kmatrix bounds, exiting.'
384	print*,"k=",k," tmid(k)=",tmid(k)
385	print*,"tgasmin=",tgasmin
386	if (strictboundcorrk) then
387	call abort
388	else
389	print,'**********************************************'
390	print*,'we allow model to continue with tmid=tgasmin'
391	print*,' ... we assume we know what you are doing ... '
392	print*,' ... but do not let this happen too often ... '
393	print,'**********************************************'
394	tmid(k)=tgasmin
395	endif
396	elseif(tmid(k).gt.tgasmax)then
397	! print*,'Maximum temperature is outside the radiative'
398	! print*,'transfer kmatrix bounds, exiting.'
399	! print*,"k=",k," tmid(k)=",tmid(k)
400	! print*,"tgasmax=",tgasmax
401	if (strictboundcorrk) then
402	call abort
403	else
404	! print,'**********************************************'
405	! print*,'we allow model to continue with tmid=tgasmin'
406	! print*,' ... we assume we know what you are doing ... '
407	! print*,' ... but do not let this happen too often ... '
408	! print,'**********************************************'
409	tmid(k)=tgasmax
410	endif
411	endif
412	enddo
413
414	!=======================================================================
415	! III. Calling the main radiative transfer subroutines
416	!=======================================================================
417
418	Cmk(:) = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
419	gzlat_ig(:) = gzlat(ig,:)
420
421	! Compute the haze seasonal modulation factor
422	if (seashaze) call season_haze(zday,latitude(ig),plevrad,seashazefact)
423
424	!-----------------------------------------------------------------------
425	! Short Wave Part
426	!-----------------------------------------------------------------------
427
428	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
429	if((ngrid.eq.1).and.(global1d))then
430	do nw=1,L_NSPECTV
431	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
432	end do
433	else
434	do nw=1,L_NSPECTV
435	stel_fract(nw)= stel(nw) * fract(ig)
436	end do
437	endif
438
439	call optcv(pqmo(ig,:,:),nlayer,plevrad,tmid,pmid, &
440	dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf,seashazefact)
441
442	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
443	acosz,stel_fract, &
444	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
445	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
446
447	else ! During the night, fluxes = 0.
448	nfluxtopv = 0.0d0
449	fluxtopvdn = 0.0d0
450	nfluxoutv_nu(:) = 0.0d0
451	nfluxgndv_nu(:) = 0.0d0
452	do l=1,L_NLAYRAD
453	fmnetv(l)=0.0d0
454	fluxupv(l)=0.0d0
455	fluxdnv(l)=0.0d0
456	end do
457	end if
458
459
460	! Equivalent Albedo Calculation (for OUTPUT). MT2015
461	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
462	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
463	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
464	DO nw=1,L_NSPECTV
465	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
466	ENDDO
467	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
468	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
469	else
470	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
471	endif
472	else
473	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
474	endif
475
476
477	!-----------------------------------------------------------------------
478	! Long Wave Part
479	!-----------------------------------------------------------------------
480
481	call optci(pqmo(ig,:,:),nlayer,plevrad,tlevrad,tmid,pmid, &
482	dtaui,taucumi,cosbi,wbari,taugsurfi,seashazefact)
483
484	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
485	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
486	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
487
488	!-----------------------------------------------------------------------
489	! Transformation of the correlated-k code outputs
490	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
491
492	! Flux incident at the top of the atmosphere
493	fluxtop_dn(ig)=fluxtopvdn
494
495	fluxtop_lw(ig) = real(nfluxtopi)
496	fluxabs_sw(ig) = real(-nfluxtopv)
497	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
498	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
499
500	! Flux absorbed by the surface. By MT2015.
501	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
502
503	if(fluxtop_dn(ig).lt.0.0)then
504	print*,'Achtung! fluxtop_dn has lost the plot!'
505	print*,'fluxtop_dn=',fluxtop_dn(ig)
506	print*,'acosz=',acosz
507	print*,'temp= ',pt(ig,:)
508	print*,'pplay= ',pplay(ig,:)
509	call abort
510	endif
511
512	! Spectral output, for exoplanet observational comparison
513	if(specOLR)then
514	do nw=1,L_NSPECTI
515	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
516	end do
517	do nw=1,L_NSPECTV
518	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
519	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
520	end do
521	endif
522
523	! Finally, the heating rates
524
525	DO l=2,L_NLAYRAD
526	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
527	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
528	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
529	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
530	END DO
531
532	! These are values at top of atmosphere
533	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
534	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
535	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
536	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
537
538
539	! Optical thickness diagnostics (added by JVO)
540	if (diagdtau) then
541	do l=1,L_NLAYRAD
542	do nw=1,L_NSPECTV
543	int_dtauv(ig,l,nw) = 0.0d0
544	DO k=1,L_NGAUSS
545	! Output exp(-tau) because gweight ponderates exp and not tau itself
546	int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
547	ENDDO
548	enddo
549	do nw=1,L_NSPECTI
550	int_dtaui(ig,l,nw) = 0.0d0
551	DO k=1,L_NGAUSS
552	! Output exp(-tau) because gweight ponderates exp and not tau itself
553	int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
554	ENDDO
555	enddo
556	enddo
557	endif
558
559
560	!-----------------------------------------------------------------------
561	end do ! End of big loop over every GCM column.
562	!-----------------------------------------------------------------------
563
564
565	!-----------------------------------------------------------------------
566	! Additional diagnostics
567	!-----------------------------------------------------------------------
568
569	! IR spectral output, for exoplanet observational comparison
570	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
571
572	print*,'Saving scalar quantities in surf_vals.out...'
573	print*,'psurf = ', pplev(1,1),' Pa'
574	open(116,file='surf_vals.out')
575	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
576	real(-nfluxtopv),real(nfluxtopi)
577	close(116)
578
579
580	! USEFUL COMMENT - Do Not Remove.
581	!
582	! if(specOLR)then
583	! open(117,file='OLRnu.out')
584	! do nw=1,L_NSPECTI
585	! write(117,*) OLR_nu(1,nw)
586	! enddo
587	! close(117)
588	!
589	! open(127,file='OSRnu.out')
590	! do nw=1,L_NSPECTV
591	! write(127,*) OSR_nu(1,nw)
592	! enddo
593	! close(127)
594	! endif
595
596	! OLR vs altitude: do it as a .txt file.
597	OLRz=.false.
598	if(OLRz)then
599	print*,'saving IR vertical flux for OLRz...'
600	open(118,file='OLRz_plevs.out')
601	open(119,file='OLRz.out')
602	do l=1,L_NLAYRAD
603	write(118,) plevrad(2l)
604	do nw=1,L_NSPECTI
605	write(119,*) fluxupi_nu(l,nw)
606	enddo
607	enddo
608	close(118)
609	close(119)
610	endif
611
612	endif
613
614	if (lastcall) then
615	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
616	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
617	IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb )
618	IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb )
619	IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold )
620	IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold )
621	!$OMP BARRIER
622	!$OMP MASTER
623	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
624	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
625	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
626	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
627	!$OMP END MASTER
628	!$OMP BARRIER
629	endif
630
631
632	end subroutine callcorrk

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

Download in other formats:

Original Format