Context Navigation

callcorrk.F90 @ 3661

Last change on this file since 3661 was 3659, checked in by emoisan, 9 months ago
Titan CRM: Deal with cloud fraction Fcloudy=1. EMo
Property svn:executable set to ``*
File size: 40.9 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, &
3	pt,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,popthi,popthv,poptti,popttv, &
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	tplanckmin,tplanckmax,callclouds,Fcloudy
21	use geometry_mod, only: latitude
22
23	implicit none
24
25	!==================================================================
26	!
27	! Purpose
28	! -------
29	! Solve the radiative transfer using the correlated-k method for
30	! the gaseous absorption and the Toon et al. (1989) method for
31	! scatttering due to aerosols.
32	!
33	! Authors
34	! -------
35	! Emmanuel 01/2001, Forget 09/2001
36	! Robin Wordsworth (2009)
37	!
38	! Modified
39	! --------
40	! Jan Vatant d'Ollone (2018)
41	! --> corrk recombining case
42	! B. de Batz de Trenquelléon (2023)
43	! --> Titan's haze and coulds optics
44	! --> clear/dark columns method
45	! --> optical diagnostics
46	!
47	!==================================================================
48
49	!-----------------------------------------------------------------------
50	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
51	! Layer #1 is the layer near the ground.
52	! Layer #nlayer is the layer at the top.
53	!-----------------------------------------------------------------------
54
55
56	! INPUT
57	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
58	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
59	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/kg).
60	INTEGER,INTENT(IN) :: nq ! Number of tracers.
61	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
62	REAL,INTENT(IN) :: zday ! Time elapsed since Ls=0 (sols).
63	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
64	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
65	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
66	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
67	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
68	REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf).
69	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
70	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
71	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
72	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
73	logical,intent(in) :: lastcall ! Signals last call to physics.
74
75	! OUTPUT
76	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
77	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
78	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
79	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
80	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
81	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
82	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
83	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
84	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
85	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
86	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
87	REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags ().
88	REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags ().
89	! Diagnostics
90	REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
91	REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
92	REAL,INTENT(OUT) :: poptti(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
93	REAL,INTENT(OUT) :: popttv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
94
95
96
97	!-----------------------------------------------------------------------
98	! Declaration of the variables required by correlated-k subroutines
99	! Numbered from top to bottom (unlike in the GCM)
100	!-----------------------------------------------------------------------
101
102	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
103	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
104
105	! Optical values for the optci/cv subroutines
106	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
107	! IR
108	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
109	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
110	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
111	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
112	! VI
113	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
114	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
115	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
116	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
117	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
118
119	! Temporary optical values for the optci/cv subroutines (clear column)
120	REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
121	REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
122	REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
123	REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
124	REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
125	REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
126	REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
127	REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
128	REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
129	! Temporary optical values for the optci/cv subroutines (dark column)
130	REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
131	REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
132	REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
133	REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
134	REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
135	REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
136	REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
137	REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
138	REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
139
140	! Optical diagnostics
141	! Haze
142	REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6)
143	REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6)
144	! Clouds
145	REAL*8 diag_optti(L_LEVELS,L_NSPECTI,6)
146	REAL*8 diag_opttv(L_LEVELS,L_NSPECTV,6)
147	! Temporary optical diagnostics (clear column)
148	REAL*8 diag_optti_cc(L_LEVELS,L_NSPECTI,6)
149	REAL*8 diag_opttv_cc(L_LEVELS,L_NSPECTV,6)
150	! Temporary optical diagnostics (dark column)
151	REAL*8 diag_optti_dc(L_LEVELS,L_NSPECTI,6)
152	REAL*8 diag_opttv_dc(L_LEVELS,L_NSPECTV,6)
153
154
155	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
156	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
157	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
158	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
159	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
160	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
161	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
162	REAL*8 albi,acosz
163	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
164
165	INTEGER ig,l,k,nw,iq,ip,ilay,it,lev2lay,cdcolumn,ng
166
167	LOGICAL found
168
169	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
170	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
171
172	! Miscellaneous :
173	character(len=100) :: message
174	character(len=10),parameter :: subname="callcorrk"
175
176	logical OLRz
177	real*8 NFLUXGNDV_nu(L_NSPECTV)
178
179	! Included by MT for albedo calculations.
180	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
181	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
182
183	! For variable haze
184	REAL*8 seashazefact(L_LEVELS)
185
186	! For muphys optics
187	REAL*8 pqmo(ngrid,nlayer,nmicro) ! Tracers for microphysics optics (X/m2).
188	REAL*8 i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics )
189
190	! For corr-k recombining
191	REAL*8 pqr(ngrid,L_PINT,L_REFVAR) ! Tracers for corr-k recombining (mol/mol).
192	REAL*8 fact, tmin, tmax
193
194	LOGICAL usept(L_PINT,L_NTREF) ! mask if pfref grid point will be used
195	INTEGER inflay(L_PINT) ! nearest inferior GCM layer for pfgasref grid points
196
197	!=======================================================================
198	! I. Initialization on every call
199	!=======================================================================
200
201
202	! How much light do we get ?
203	do nw=1,L_NSPECTV
204	stel(nw)=stellarf(nw)/(dist_star**2)
205	end do
206
207	! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2
208	! NOTE: it should be moved somewhere else: calmufi performs the same kind of
209	! computations... waste of time...
210	i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
211	pqmo(:,:,:) = 0.0
212	DO iq=1,nmicro
213	pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:)
214	ENDDO
215
216	! Default value for fixed species for whom vmr has been taken
217	! into account while computing high-resolution spectra
218	if (corrk_recombin) pqr(:,:,:) = 1.0
219
220	!-----------------------------------------------------------------------
221	do ig=1,ngrid ! Starting Big Loop over every GCM column
222	!-----------------------------------------------------------------------
223
224	! Recombine reference corr-k if needed
225	if (corrk_recombin) then
226
227	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
228	! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values
229	! who match the atmospheric conditions ) which is then processed as a standard pre-mix in
230	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
231
232	! Extract tracers for variable radiative species
233	! Also find the nearest GCM layer under each ref pressure
234	do ip=1,L_PINT
235
236	ilay=0
237	found = .false.
238	do l=1,nlayer
239	if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar])
240	found=.true.
241	ilay=l
242	endif
243	enddo
244
245	if (.not. found ) then ! set to min
246	do iq=1,L_REFVAR
247	if ( radvar_mask(iq) ) then
248	pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
249	endif
250	enddo
251	else
252	if (ilay==nlayer) then ! set to max
253	do iq=1,L_REFVAR
254	if ( radvar_mask(iq) ) then
255	pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
256	endif
257	enddo
258	else ! standard
259	fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(ig,ilay+1) ) / ( pplay(ig,ilay) - pplay(ig,ilay+1) ) ! pfgasref=log(p[mbar])
260	do iq=1,L_REFVAR
261	if ( radvar_mask(iq) ) then
262	pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))*fact pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact)
263	pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
264	endif
265	enddo
266	endif ! if ilay==nlayer
267	endif ! if not found
268
269	inflay(ip) = ilay
270
271	enddo ! ip=1,L_PINT
272
273	! NB : The following usept is a trick to call recombine only for the reference T-P
274	! grid points that are useful given the temperature range at this altitude
275	! It saves a looot of time - JVO 18
276	usept(:,:) = .true.
277	do ip=1,L_PINT-1
278	if ( inflay(ip+1)==nlayer ) then
279	usept(ip,:) = .false.
280	endif
281	if ( inflay(ip) == 0 ) then
282	usept(ip+1:,:) = .false.
283	endif
284	if ( usept(ip,1) ) then ! if not all false
285	tmin = minval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer)))
286	tmax = maxval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer)))
287	do it=1,L_NTREF-1
288	if ( tgasref(it+1) .lt. tmin ) then
289	usept(ip,it) = .false.
290	endif
291	enddo
292	do it=2,L_NTREF
293	if ( tgasref(it-1) .gt. tmax ) then
294	usept(ip,it) = .false.
295	endif
296	enddo
297	! in case of out-of-bounds
298	if ( tgasref(1) .lt. tmin ) usept(ip,1) = .true.
299	if ( tgasref(L_NTREF) .gt. tmax ) usept(ip,L_NTREF) = .true.
300	endif
301	enddo ! ip=1,L_PINT-1
302	! deal with last bound
303	if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:)
304
305
306	do ip=1,L_PINT
307
308	! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied
309	do it=1,L_NTREF
310
311	if ( usept(ip,it) .eqv. .false. ) cycle
312
313	do l=1,L_REFVAR
314	if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01 & ! +- 1%
315	.or. ( useptold(ip,it) .eqv. .false. ) ) then ! in case T change but not the tracers
316	call recombin_corrk( pqr(ig,ip,:),ip,it )
317	exit ! one is enough to trigger the update
318	endif
319	enddo
320
321	enddo
322
323	enddo ! ip=1,L_PINT
324
325	useptold(:,:)=usept(:,:)
326
327	endif ! if corrk_recombin
328
329	!=======================================================================
330	! II. Transformation of the GCM variables
331	!=======================================================================
332
333
334	! Albedo and Emissivity.
335	albi=1-emis(ig) ! Long Wave.
336	DO nw=1,L_NSPECTV ! Short Wave loop.
337	albv(nw)=albedo(ig,nw)
338	ENDDO
339
340	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
341	acosz = cos(pi*szangle/180.0)
342	print*,'acosz=',acosz,', szangle=',szangle
343	else
344	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
345	endif
346
347	!-----------------------------------------------------------------------
348	! Pressure and temperature
349	!-----------------------------------------------------------------------
350
351	DO l=1,nlayer
352	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
353	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
354	tlevrad(2*l) = pt(ig,nlayer+1-l)
355	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
356	END DO
357
358	plevrad(1) = 0.
359	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.
360
361	tlevrad(1) = tlevrad(2)
362	tlevrad(2*nlayer+1)=tsurf(ig)
363
364	pmid(1) = max(pgasmin,0.0001*plevrad(3))
365	pmid(2) = pmid(1)
366
367	tmid(1) = tlevrad(2)
368	tmid(2) = tmid(1)
369
370	DO l=1,L_NLAYRAD-1
371	tmid(2l+1) = tlevrad(2l+1)
372	tmid(2l+2) = tlevrad(2l+1)
373	pmid(2l+1) = plevrad(2l+1)
374	pmid(2l+2) = plevrad(2l+1)
375	END DO
376	pmid(L_LEVELS) = plevrad(L_LEVELS)
377	tmid(L_LEVELS) = tlevrad(L_LEVELS)
378
379	!!Alternative interpolation:
380	! pmid(3) = pmid(1)
381	! pmid(4) = pmid(1)
382	! tmid(3) = tmid(1)
383	! tmid(4) = tmid(1)
384	! DO l=2,L_NLAYRAD-1
385	! tmid(2l+1) = tlevrad(2l)
386	! tmid(2l+2) = tlevrad(2l)
387	! pmid(2l+1) = plevrad(2l)
388	! pmid(2l+2) = plevrad(2l)
389	! END DO
390	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
391	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
392
393	! Test for out-of-bounds pressure.
394	if(plevrad(3).lt.pgasmin)then
395	print*,'Minimum pressure is outside the radiative'
396	print*,'transfer kmatrix bounds, exiting.'
397	call abort
398	elseif(plevrad(L_LEVELS).gt.pgasmax)then
399	print*,'Maximum pressure is outside the radiative'
400	print*,'transfer kmatrix bounds, exiting.'
401	call abort
402	endif
403
404	! Test for out-of-bounds temperature.
405	do k=1,L_LEVELS
406	if(tlevrad(k).lt.tgasmin)then
407	print*,'Minimum temperature is outside the radiative'
408	print*,'transfer kmatrix bounds'
409	print*,"k=",k," tlevrad(k)=",tlevrad(k)
410	print*,"tgasmin=",tgasmin
411	if (strictboundcorrk) then
412	call abort
413	else
414	print,'**********************************************'
415	print*,'we allow model to continue with tlevrad=tgasmin'
416	print*,' ... we assume we know what you are doing ... '
417	print*,' ... but do not let this happen too often ... '
418	print,'**********************************************'
419	!tlevrad(k)=tgasmin
420	endif
421	elseif(tlevrad(k).gt.tgasmax)then
422	! print*,'Maximum temperature is outside the radiative'
423	! print*,'transfer kmatrix bounds, exiting.'
424	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
425	! print*,"tgasmax=",tgasmax
426	if (strictboundcorrk) then
427	call abort
428	else
429	! print,'**********************************************'
430	! print*,'we allow model to continue with tlevrad=tgasmax'
431	! print*,' ... we assume we know what you are doing ... '
432	! print*,' ... but do not let this happen too often ... '
433	! print,'**********************************************'
434	!tlevrad(k)=tgasmax
435	endif
436	endif
437
438	if (tlevrad(k).lt.tplanckmin) then
439	print*,'Minimum temperature is outside the boundaries for'
440	print*,'Planck function integration set in callphys.def,'
441	print*,' aborting.'
442	print*,"k=",k," tlevrad(k)=",tlevrad(k)
443	print*,"tplanckmin=",tplanckmin
444	message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
445	call abort_physic(subname,message,1)
446	else if (tlevrad(k).gt.tplanckmax) then
447	print*,'Maximum temperature is outside the boundaries for'
448	print*,'Planck function integration set in callphys.def, aborting.'
449	print*,"k=",k," tlevrad(k)=",tlevrad(k)
450	print*,"tplanckmax=",tplanckmax
451	message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
452	call abort_physic(subname,message,1)
453	endif
454
455	enddo
456
457	do k=1,L_NLAYRAD+1
458	if(tmid(k).lt.tgasmin)then
459	print*,'Minimum temperature is outside the radiative'
460	print*,'transfer kmatrix bounds, exiting.'
461	print*,"k=",k," tmid(k)=",tmid(k)
462	print*,"tgasmin=",tgasmin
463	if (strictboundcorrk) then
464	call abort
465	else
466	print,'**********************************************'
467	print*,'we allow model to continue with tmid=tgasmin'
468	print*,' ... we assume we know what you are doing ... '
469	print*,' ... but do not let this happen too often ... '
470	print,'**********************************************'
471	tmid(k)=tgasmin
472	endif
473	elseif(tmid(k).gt.tgasmax)then
474	! print*,'Maximum temperature is outside the radiative'
475	! print*,'transfer kmatrix bounds, exiting.'
476	! print*,"k=",k," tmid(k)=",tmid(k)
477	! print*,"tgasmax=",tgasmax
478	if (strictboundcorrk) then
479	call abort
480	else
481	! print,'**********************************************'
482	! print*,'we allow model to continue with tmid=tgasmin'
483	! print*,' ... we assume we know what you are doing ... '
484	! print*,' ... but do not let this happen too often ... '
485	! print,'**********************************************'
486	tmid(k)=tgasmax
487	endif
488	endif
489	enddo
490
491	!=======================================================================
492	! III. Calling the main radiative transfer subroutines
493	!=======================================================================
494
495	Cmk(:) = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
496	gzlat_ig(:) = gzlat(ig,:)
497
498	! Compute the haze seasonal modulation factor
499	if (seashaze) call season_haze(zday,latitude(ig),plevrad,seashazefact)
500
501	!-----------------------------------------------------------------------
502	! Short Wave Part
503	!-----------------------------------------------------------------------
504
505	! Clear column :
506	cdcolumn = 0
507	call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, &
508	dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
509	diag_opthv,diag_opttv_cc,cdcolumn)
510	! Dark column :
511	cdcolumn = 1
512	call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, &
513	dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
514	diag_opthv,diag_opttv_dc,cdcolumn)
515
516	! Mean opacity, ssa and asf :
517	where ((exp(-dtauv_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtauv_dc(:,:,:)).ge.1.d-40))
518	dtauv(:,:,:) = -log((1-Fcloudy)exp(-dtauv_cc(:,:,:)) + Fcloudyexp(-dtauv_dc(:,:,:)))
519	elsewhere
520	dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
521	endwhere
522	do ng = 1, L_NGAUSS
523	do nw = 1, L_NSPECTV
524	taucumv(1,nw,ng) = 0.0d0
525	do k = 2, L_LEVELS
526	if ((exp(-taucumv_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumv_dc(k,nw,ng)).ge.1.d-40)) then
527	taucumv(k,nw,ng) = taucumv(k-1,nw,ng) - log((1-Fcloudy)exp(-taucumv_cc(k,nw,ng)) + Fcloudyexp(-taucumv_dc(k,nw,ng)))
528	else
529	taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
530	end if
531	end do
532	do l = 1, L_NLAYRAD
533	tauv(l,nw,ng) = taucumv(2*l,nw,ng)
534	end do
535	tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng)
536	end do
537	end do
538	wbarv = (1-Fcloudy) * wbarv_cc + Fcloudy * wbarv_dc
539	cosbv = (1-Fcloudy) * cosbv_cc + Fcloudy * cosbv_dc
540
541	! Diagnostics for clouds :
542	if (callclouds) then
543	where (diag_opttv_cc(:,:,1) .lt. 1.d-30)
544	diag_opttv_cc(:,:,1) = 1.d-30
545	endwhere
546	where (diag_opttv_dc(:,:,1) .lt. 1.d-30)
547	diag_opttv_dc(:,:,1) = 1.d-30
548	endwhere
549	if (Fcloudy .eq. 1) then !avoid log(0) errors
550	diag_opttv(:,:,1) = diag_opttv_dc(:,:,1)
551	diag_opttv(:,:,2) = diag_opttv_dc(:,:,2)
552	diag_opttv(:,:,3) = diag_opttv_dc(:,:,3)
553	diag_opttv(:,:,4) = diag_opttv_dc(:,:,4)
554	diag_opttv(:,:,5) = diag_opttv_dc(:,:,5)
555	diag_opttv(:,:,6) = diag_opttv_dc(:,:,6)
556	else !correct for cloud fraction
557	diag_opttv(:,:,1) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,1)) + Fcloudyexp(-diag_opttv_dc(:,:,1)) )
558	diag_opttv(:,:,2) = (1-Fcloudy) * diag_opttv_cc(:,:,2) + Fcloudy * diag_opttv_dc(:,:,2)
559	diag_opttv(:,:,3) = (1-Fcloudy) * diag_opttv_cc(:,:,3) + Fcloudy * diag_opttv_dc(:,:,3)
560	diag_opttv(:,:,4) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,4)) + Fcloudyexp(-diag_opttv_dc(:,:,4)) )
561	diag_opttv(:,:,5) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,5)) + Fcloudyexp(-diag_opttv_dc(:,:,5)) )
562	diag_opttv(:,:,6) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,6)) + Fcloudyexp(-diag_opttv_dc(:,:,6)) )
563	endif
564	endif
565
566	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
567	if((ngrid.eq.1).and.(global1d))then
568	do nw=1,L_NSPECTV
569	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
570	end do
571	else
572	do nw=1,L_NSPECTV
573	stel_fract(nw)= stel(nw) * fract(ig)
574	end do
575	endif
576	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
577	acosz,stel_fract, &
578	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
579	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
580
581	else ! During the night, fluxes = 0.
582	nfluxtopv = 0.0d0
583	fluxtopvdn = 0.0d0
584	nfluxoutv_nu(:) = 0.0d0
585	nfluxgndv_nu(:) = 0.0d0
586	do l=1,L_NLAYRAD
587	fmnetv(l)=0.0d0
588	fluxupv(l)=0.0d0
589	fluxdnv(l)=0.0d0
590	end do
591	end if
592
593
594	! Equivalent Albedo Calculation (for OUTPUT). MT2015
595	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
596	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
597	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
598	DO nw=1,L_NSPECTV
599	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
600	ENDDO
601	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
602	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
603	else
604	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
605	endif
606	else
607	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
608	endif
609
610
611	!-----------------------------------------------------------------------
612	! Long Wave Part
613	!-----------------------------------------------------------------------
614
615
616	! Clear column :
617	cdcolumn = 0
618	call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
619	dtaui_cc,taucumi_cc,cosbi_cc,wbari_cc,taugsurfi,seashazefact, &
620	diag_opthi,diag_optti_cc,cdcolumn)
621	! Dark column :
622	cdcolumn = 1
623	call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
624	dtaui_dc,taucumi_dc,cosbi_dc,wbari_dc,taugsurfi,seashazefact, &
625	diag_opthi,diag_optti_dc,cdcolumn)
626
627	! Mean opacity, ssa and asf :
628	where ((exp(-dtaui_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtaui_dc(:,:,:)).ge.1.d-40))
629	dtaui(:,:,:) = -log((1-Fcloudy)exp(-dtaui_cc(:,:,:)) + Fcloudyexp(-dtaui_dc(:,:,:)))
630	elsewhere
631	dtaui(:,:,:) = dtaui_dc(:,:,:) ! No need to average...
632	endwhere
633	do ng = 1, L_NGAUSS
634	do nw = 1, L_NSPECTI
635	taucumi(1,nw,ng) = 0.0d0
636	do k = 2, L_LEVELS
637	if ((exp(-taucumi_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumi_dc(k,nw,ng)).ge.1.d-40)) then
638	taucumi(k,nw,ng) = taucumi(k-1,nw,ng) - log((1-Fcloudy)exp(-taucumi_cc(k,nw,ng)) + Fcloudyexp(-taucumi_dc(k,nw,ng)))
639	else
640	taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + taucumi_dc(k,nw,ng) ! No need to average...
641	end if
642	end do
643	end do
644	end do
645	wbari = (1-Fcloudy) * wbari_cc + Fcloudy * wbari_dc
646	cosbi = (1-Fcloudy) * cosbi_cc + Fcloudy * cosbi_dc
647
648	! Diagnostics for clouds :
649	if (callclouds) then
650	where (diag_optti_cc(:,:,1) .lt. 1.d-30)
651	diag_optti_cc(:,:,1) = 1.d-30
652	endwhere
653	where (diag_optti_dc(:,:,1) .lt. 1.d-30)
654	diag_optti_dc(:,:,1) = 1.d-30
655	endwhere
656	if (Fcloudy .eq. 1) then !avoid log(0) errors
657	diag_optti(:,:,1) = diag_optti_dc(:,:,1)
658	diag_optti(:,:,2) = diag_optti_dc(:,:,2)
659	diag_optti(:,:,3) = diag_optti_dc(:,:,3)
660	diag_optti(:,:,4) = diag_optti_dc(:,:,4)
661	diag_optti(:,:,5) = diag_optti_dc(:,:,5)
662	diag_optti(:,:,6) = diag_optti_dc(:,:,6)
663	else !correct for cloud fraction
664	diag_optti(:,:,1) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,1)) + Fcloudyexp(-diag_optti_dc(:,:,1)) )
665	diag_optti(:,:,2) = (1-Fcloudy) * diag_optti_cc(:,:,2) + Fcloudy * diag_optti_dc(:,:,2)
666	diag_optti(:,:,3) = (1-Fcloudy) * diag_optti_cc(:,:,3) + Fcloudy * diag_optti_dc(:,:,3)
667	diag_optti(:,:,4) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,4)) + Fcloudyexp(-diag_optti_dc(:,:,4)) )
668	diag_optti(:,:,5) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,5)) + Fcloudyexp(-diag_optti_dc(:,:,5)) )
669	diag_optti(:,:,6) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,6)) + Fcloudyexp(-diag_optti_dc(:,:,6)) )
670	endif
671	endif
672
673	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
674	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
675	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
676
677	!-----------------------------------------------------------------------
678	! Transformation of the correlated-k code outputs
679	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
680
681	! Flux incident at the top of the atmosphere
682	fluxtop_dn(ig)=fluxtopvdn
683
684	fluxtop_lw(ig) = real(nfluxtopi)
685	fluxabs_sw(ig) = real(-nfluxtopv)
686	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
687	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
688
689	! Flux absorbed by the surface. By MT2015.
690	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
691
692	if(fluxtop_dn(ig).lt.0.0)then
693	print*,'Achtung! fluxtop_dn has lost the plot!'
694	print*,'fluxtop_dn=',fluxtop_dn(ig)
695	print*,'acosz=',acosz
696	print*,'temp= ',pt(ig,:)
697	print*,'pplay= ',pplay(ig,:)
698	call abort
699	endif
700
701	! Spectral output, for exoplanet observational comparison
702	if(specOLR)then
703	do nw=1,L_NSPECTI
704	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
705	end do
706	do nw=1,L_NSPECTV
707	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
708	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
709	end do
710	endif
711
712	! Finally, the heating rates
713
714	DO l=2,L_NLAYRAD
715	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
716	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
717	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
718	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
719	END DO
720
721	! These are values at top of atmosphere
722	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
723	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
724	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
725	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
726
727
728	! Optical thickness diagnostics (added by JVO)
729	if (diagdtau) then
730	do l=1,L_NLAYRAD
731	do nw=1,L_NSPECTV
732	int_dtauv(ig,l,nw) = 0.0d0
733	DO k=1,L_NGAUSS
734	! Output exp(-tau) because gweight ponderates exp and not tau itself
735	int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
736	ENDDO
737	enddo
738	do nw=1,L_NSPECTI
739	int_dtaui(ig,l,nw) = 0.0d0
740	DO k=1,L_NGAUSS
741	! Output exp(-tau) because gweight ponderates exp and not tau itself
742	int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
743	ENDDO
744	enddo
745	enddo
746	endif
747
748
749	! Diagnostics : optical properties of haze and clouds (dtau, tau, k, w, g) :
750	do l = 2, L_LEVELS, 2
751	lev2lay = L_NLAYRAD+1 - l/2
752	! Visible :
753	do nw = 1, L_NSPECTV
754	popthv(ig,lev2lay,nw,:) = 0.0d0
755	popttv(ig,lev2lay,nw,:) = 0.0d0
756	! Optical thickness (dtau) :
757	popthv(ig,lev2lay,nw,1) = (diag_opthv(l,nw,1) + diag_opthv(l+1,nw,1)) / 2.0
758	if (callclouds) then
759	popttv(ig,lev2lay,nw,1) = (diag_opttv(l,nw,1) + diag_opttv(l+1,nw,1)) / 2.0
760	endif
761	! Opacity (tau) :
762	do k = L_NLAYRAD, lev2lay, -1
763	popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
764	if (callclouds) then
765	popttv(ig,lev2lay,nw,2) = popttv(ig,lev2lay,nw,2) + popttv(ig,k,nw,1)
766	endif
767	enddo
768	! Extinction (k) :
769	popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
770	if (callclouds) then
771	popttv(ig,lev2lay,nw,3) = popttv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
772	endif
773	! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
774	popthv(ig,lev2lay,nw,4) = (diag_opthv(l,nw,2) + diag_opthv(l+1,nw,2)) / 2.0
775	popthv(ig,lev2lay,nw,5) = (diag_opthv(l,nw,3) + diag_opthv(l+1,nw,3)) / 2.0
776	if (callclouds) then
777	popttv(ig,lev2lay,nw,4) = (diag_opttv(l,nw,2) + diag_opttv(l+1,nw,2)) / 2.0
778	popttv(ig,lev2lay,nw,5) = (diag_opttv(l,nw,3) + diag_opttv(l+1,nw,3)) / 2.0
779	endif
780	! DRAYAER, TAUGAS, DCONT :
781	popthv(ig,lev2lay,nw,6) = (diag_opthv(l,nw,4) + diag_opthv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
782	popthv(ig,lev2lay,nw,7) = (diag_opthv(l,nw,5) + diag_opthv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
783	popthv(ig,lev2lay,nw,8) = (diag_opthv(l,nw,6) + diag_opthv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
784	if (callclouds) then
785	popttv(ig,lev2lay,nw,6) = (diag_opttv(l,nw,4) + diag_opttv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
786	popttv(ig,lev2lay,nw,7) = (diag_opttv(l,nw,5) + diag_opttv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
787	popttv(ig,lev2lay,nw,8) = (diag_opttv(l,nw,6) + diag_opttv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
788	endif
789	enddo
790	! Infra-Red
791	do nw=1,L_NSPECTI
792	popthi(ig,lev2lay,nw,:) = 0.0d0
793	poptti(ig,lev2lay,nw,:) = 0.0d0
794	! Optical thickness (dtau) :
795	popthi(ig,lev2lay,nw,1) = (diag_opthi(l,nw,1) + diag_opthi(l+1,nw,1)) / 2.0
796	if (callclouds) then
797	poptti(ig,lev2lay,nw,1) = (diag_optti(l,nw,1) + diag_optti(l+1,nw,1)) / 2.0
798	endif
799	! Opacity (tau) :
800	do k = L_NLAYRAD, lev2lay, -1
801	popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
802	if (callclouds) then
803	poptti(ig,lev2lay,nw,2) = poptti(ig,lev2lay,nw,2) + poptti(ig,k,nw,1)
804	endif
805	enddo
806	! Extinction (k) :
807	popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
808	if (callclouds) then
809	poptti(ig,lev2lay,nw,3) = poptti(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
810	endif
811	! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
812	popthi(ig,lev2lay,nw,4) = (diag_opthi(l,nw,2) + diag_opthi(l+1,nw,2)) / 2.0
813	popthi(ig,lev2lay,nw,5) = (diag_opthi(l,nw,3) + diag_opthi(l+1,nw,3)) / 2.0
814	if (callclouds) then
815	poptti(ig,lev2lay,nw,4) = (diag_optti(l,nw,2) + diag_optti(l+1,nw,2)) / 2.0
816	poptti(ig,lev2lay,nw,5) = (diag_optti(l,nw,3) + diag_optti(l+1,nw,3)) / 2.0
817	endif
818	! DRAYAER, TAUGAS, DCONT :
819	popthi(ig,lev2lay,nw,6) = (diag_opthi(l,nw,4) + diag_opthi(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
820	popthi(ig,lev2lay,nw,7) = (diag_opthi(l,nw,5) + diag_opthi(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
821	popthi(ig,lev2lay,nw,8) = (diag_opthi(l,nw,6) + diag_opthi(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
822	if (callclouds) then
823	poptti(ig,lev2lay,nw,6) = (diag_optti(l,nw,4) + diag_optti(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
824	poptti(ig,lev2lay,nw,7) = (diag_optti(l,nw,5) + diag_optti(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
825	poptti(ig,lev2lay,nw,8) = (diag_optti(l,nw,6) + diag_optti(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
826	endif
827	enddo
828	enddo
829
830
831	!-----------------------------------------------------------------------
832	end do ! End of big loop over every GCM column.
833	!-----------------------------------------------------------------------
834
835
836	!-----------------------------------------------------------------------
837	! Additional diagnostics
838	!-----------------------------------------------------------------------
839
840	! IR spectral output, for exoplanet observational comparison
841	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
842
843	print*,'Saving scalar quantities in surf_vals.out...'
844	print*,'psurf = ', pplev(1,1),' Pa'
845	open(116,file='surf_vals.out')
846	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
847	real(-nfluxtopv),real(nfluxtopi)
848	close(116)
849
850	! USEFUL COMMENT - Do Not Remove.
851	!
852	! if(specOLR)then
853	! open(117,file='OLRnu.out')
854	! do nw=1,L_NSPECTI
855	! write(117,*) OLR_nu(1,nw)
856	! enddo
857	! close(117)
858	!
859	! open(127,file='OSRnu.out')
860	! do nw=1,L_NSPECTV
861	! write(127,*) OSR_nu(1,nw)
862	! enddo
863	! close(127)
864	! endif
865
866	! OLR vs altitude: do it as a .txt file.
867	OLRz=.false.
868	if(OLRz)then
869	print*,'saving IR vertical flux for OLRz...'
870	open(118,file='OLRz_plevs.out')
871	open(119,file='OLRz.out')
872	do l=1,L_NLAYRAD
873	write(118,) plevrad(2l)
874	do nw=1,L_NSPECTI
875	write(119,*) fluxupi_nu(l,nw)
876	enddo
877	enddo
878	close(118)
879	close(119)
880	endif
881
882	endif
883
884	if (lastcall) then
885	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
886	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
887	IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb )
888	IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb )
889	IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold )
890	IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold )
891	!$OMP BARRIER
892	!$OMP MASTER
893	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
894	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
895	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
896	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
897	!$OMP END MASTER
898	!$OMP BARRIER
899	endif
900
901
902	end subroutine callcorrk

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

Download in other formats:

Original Format