Context Navigation

callcorrk.F90 @ 1902

Last change on this file since 1902 was 1897, checked in by jvatant, 8 years ago
Making Titan's hazy again - part II + Major updates of J.Burgalat YAMMS library and optical coupling, including : ++ Added the routines for haze optics inside YAMMS ++ Calling rad. transf. with interactive haze is plugged in but should stay unactive as long as the microphysics is in test phase : cf "uncoupl_optic_haze" flag : true for now ! ++ Also some sanity checks for negative tendencies and some others upkeep of YAMMS model + Also added a temporary CPP key USE_QTEST in physiq_mod that enables to have microphysical tendencies separated from dynamics for debugging and test phases -- JVO and JB
Property svn:executable set to ``*
File size: 19.2 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
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	lastcall)
9
10	use mod_phys_lmdz_para, only : is_master
11	use radinc_h
12	use radcommon_h
13	use gases_h
14	USE tracer_h
15	use callkeys_mod, only: global1d, szangle
16	use comcstfi_mod, only: pi, mugaz, cpp
17	use callkeys_mod, only: diurnal,tracer, &
18	strictboundcorrk,specOLR
19
20	implicit none
21
22	!==================================================================
23	!
24	! Purpose
25	! -------
26	! Solve the radiative transfer using the correlated-k method for
27	! the gaseous absorption and the Toon et al. (1989) method for
28	! scatttering due to aerosols.
29	!
30	! Authors
31	! -------
32	! Emmanuel 01/2001, Forget 09/2001
33	! Robin Wordsworth (2009)
34	!
35	!==================================================================
36
37	!-----------------------------------------------------------------------
38	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
39	! Layer #1 is the layer near the ground.
40	! Layer #nlayer is the layer at the top.
41	!-----------------------------------------------------------------------
42
43
44	! INPUT
45	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
46	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
47	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/m2).
48	INTEGER,INTENT(IN) :: nq ! Number of tracers.
49	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
50	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
51	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
52	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
53	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
54	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
55	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
56	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
57	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
58	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
59	logical,intent(in) :: lastcall ! Signals last call to physics.
60
61	! OUTPUT
62	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
63	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
64	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
65	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
66	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
67	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
68	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
69	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
70	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
71	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
72	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
73
74
75	!-----------------------------------------------------------------------
76	! Declaration of the variables required by correlated-k subroutines
77	! Numbered from top to bottom (unlike in the GCM)
78	!-----------------------------------------------------------------------
79
80	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
81	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
82
83	! Optical values for the optci/cv subroutines
84	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
85	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
86	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
87	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
88	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
89	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
90	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
91	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
92	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
93	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
94
95	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
96	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
97	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
98	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
99	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
100	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
101	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
102	REAL*8 albi,acosz
103	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
104
105	INTEGER ig,l,k,nw
106
107	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
108	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
109
110	logical OLRz
111	real*8 NFLUXGNDV_nu(L_NSPECTV)
112
113	! Included by MT for albedo calculations.
114	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
115	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
116
117
118	!=======================================================================
119	! I. Initialization on every call
120	!=======================================================================
121
122
123	! How much light do we get ?
124	do nw=1,L_NSPECTV
125	stel(nw)=stellarf(nw)/(dist_star**2)
126	end do
127
128
129	!-----------------------------------------------------------------------
130	do ig=1,ngrid ! Starting Big Loop over every GCM column
131	!-----------------------------------------------------------------------
132
133
134	!=======================================================================
135	! II. Transformation of the GCM variables
136	!=======================================================================
137
138
139	! Albedo and Emissivity.
140	albi=1-emis(ig) ! Long Wave.
141	DO nw=1,L_NSPECTV ! Short Wave loop.
142	albv(nw)=albedo(ig,nw)
143	ENDDO
144
145	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
146	acosz = cos(pi*szangle/180.0)
147	print*,'acosz=',acosz,', szangle=',szangle
148	else
149	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
150	endif
151
152	!-----------------------------------------------------------------------
153	! Pressure and temperature
154	!-----------------------------------------------------------------------
155
156	DO l=1,nlayer
157	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
158	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
159	tlevrad(2*l) = pt(ig,nlayer+1-l)
160	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
161	END DO
162
163	plevrad(1) = 0.
164	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.
165
166	tlevrad(1) = tlevrad(2)
167	tlevrad(2*nlayer+1)=tsurf(ig)
168
169	pmid(1) = max(pgasmin,0.0001*plevrad(3))
170	pmid(2) = pmid(1)
171
172	tmid(1) = tlevrad(2)
173	tmid(2) = tmid(1)
174
175	DO l=1,L_NLAYRAD-1
176	tmid(2l+1) = tlevrad(2l+1)
177	tmid(2l+2) = tlevrad(2l+1)
178	pmid(2l+1) = plevrad(2l+1)
179	pmid(2l+2) = plevrad(2l+1)
180	END DO
181	pmid(L_LEVELS) = plevrad(L_LEVELS)
182	tmid(L_LEVELS) = tlevrad(L_LEVELS)
183
184	!!Alternative interpolation:
185	! pmid(3) = pmid(1)
186	! pmid(4) = pmid(1)
187	! tmid(3) = tmid(1)
188	! tmid(4) = tmid(1)
189	! DO l=2,L_NLAYRAD-1
190	! tmid(2l+1) = tlevrad(2l)
191	! tmid(2l+2) = tlevrad(2l)
192	! pmid(2l+1) = plevrad(2l)
193	! pmid(2l+2) = plevrad(2l)
194	! END DO
195	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
196	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
197
198	! Test for out-of-bounds pressure.
199	if(plevrad(3).lt.pgasmin)then
200	print*,'Minimum pressure is outside the radiative'
201	print*,'transfer kmatrix bounds, exiting.'
202	call abort
203	elseif(plevrad(L_LEVELS).gt.pgasmax)then
204	print*,'Maximum pressure is outside the radiative'
205	print*,'transfer kmatrix bounds, exiting.'
206	call abort
207	endif
208
209	! Test for out-of-bounds temperature.
210	do k=1,L_LEVELS
211	if(tlevrad(k).lt.tgasmin)then
212	print*,'Minimum temperature is outside the radiative'
213	print*,'transfer kmatrix bounds'
214	print*,"k=",k," tlevrad(k)=",tlevrad(k)
215	print*,"tgasmin=",tgasmin
216	if (strictboundcorrk) then
217	call abort
218	else
219	print,'**********************************************'
220	print*,'we allow model to continue with tlevrad=tgasmin'
221	print*,' ... we assume we know what you are doing ... '
222	print*,' ... but do not let this happen too often ... '
223	print,'**********************************************'
224	!tlevrad(k)=tgasmin
225	endif
226	elseif(tlevrad(k).gt.tgasmax)then
227	! print*,'Maximum temperature is outside the radiative'
228	! print*,'transfer kmatrix bounds, exiting.'
229	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
230	! print*,"tgasmax=",tgasmax
231	if (strictboundcorrk) then
232	call abort
233	else
234	! print,'**********************************************'
235	! print*,'we allow model to continue with tlevrad=tgasmax'
236	! print*,' ... we assume we know what you are doing ... '
237	! print*,' ... but do not let this happen too often ... '
238	! print,'**********************************************'
239	!tlevrad(k)=tgasmax
240	endif
241	endif
242	enddo
243	do k=1,L_NLAYRAD+1
244	if(tmid(k).lt.tgasmin)then
245	print*,'Minimum temperature is outside the radiative'
246	print*,'transfer kmatrix bounds, exiting.'
247	print*,"k=",k," tmid(k)=",tmid(k)
248	print*,"tgasmin=",tgasmin
249	if (strictboundcorrk) then
250	call abort
251	else
252	print,'**********************************************'
253	print*,'we allow model to continue with tmid=tgasmin'
254	print*,' ... we assume we know what you are doing ... '
255	print*,' ... but do not let this happen too often ... '
256	print,'**********************************************'
257	tmid(k)=tgasmin
258	endif
259	elseif(tmid(k).gt.tgasmax)then
260	! print*,'Maximum temperature is outside the radiative'
261	! print*,'transfer kmatrix bounds, exiting.'
262	! print*,"k=",k," tmid(k)=",tmid(k)
263	! print*,"tgasmax=",tgasmax
264	if (strictboundcorrk) then
265	call abort
266	else
267	! print,'**********************************************'
268	! print*,'we allow model to continue with tmid=tgasmin'
269	! print*,' ... we assume we know what you are doing ... '
270	! print*,' ... but do not let this happen too often ... '
271	! print,'**********************************************'
272	tmid(k)=tgasmax
273	endif
274	endif
275	enddo
276
277	!=======================================================================
278	! III. Calling the main radiative transfer subroutines
279	!=======================================================================
280
281
282	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
283	glat_ig=glat(ig)
284
285	!-----------------------------------------------------------------------
286	! Short Wave Part
287	!-----------------------------------------------------------------------
288
289	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
290	if((ngrid.eq.1).and.(global1d))then
291	do nw=1,L_NSPECTV
292	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
293	end do
294	else
295	do nw=1,L_NSPECTV
296	stel_fract(nw)= stel(nw) * fract(ig)
297	end do
298	endif
299
300	call optcv(pq(ig,:,1:nmicro),nlayer,plevrad,tmid,pmid, &
301	dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf)
302
303	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
304	acosz,stel_fract, &
305	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
306	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
307
308	else ! During the night, fluxes = 0.
309	nfluxtopv = 0.0d0
310	fluxtopvdn = 0.0d0
311	nfluxoutv_nu(:) = 0.0d0
312	nfluxgndv_nu(:) = 0.0d0
313	do l=1,L_NLAYRAD
314	fmnetv(l)=0.0d0
315	fluxupv(l)=0.0d0
316	fluxdnv(l)=0.0d0
317	end do
318	end if
319
320
321	! Equivalent Albedo Calculation (for OUTPUT). MT2015
322	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
323	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
324	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
325	DO nw=1,L_NSPECTV
326	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
327	ENDDO
328	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
329	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
330	else
331	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
332	endif
333	else
334	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
335	endif
336
337
338	!-----------------------------------------------------------------------
339	! Long Wave Part
340	!-----------------------------------------------------------------------
341
342	call optci(pq(ig,:,1:nmicro),nlayer,plevrad,tlevrad,tmid,pmid, &
343	dtaui,taucumi,cosbi,wbari,taugsurfi)
344
345	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
346	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
347	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
348
349	!-----------------------------------------------------------------------
350	! Transformation of the correlated-k code outputs
351	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
352
353	! Flux incident at the top of the atmosphere
354	fluxtop_dn(ig)=fluxtopvdn
355
356	fluxtop_lw(ig) = real(nfluxtopi)
357	fluxabs_sw(ig) = real(-nfluxtopv)
358	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
359	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
360
361	! Flux absorbed by the surface. By MT2015.
362	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
363
364	if(fluxtop_dn(ig).lt.0.0)then
365	print*,'Achtung! fluxtop_dn has lost the plot!'
366	print*,'fluxtop_dn=',fluxtop_dn(ig)
367	print*,'acosz=',acosz
368	print*,'temp= ',pt(ig,:)
369	print*,'pplay= ',pplay(ig,:)
370	call abort
371	endif
372
373	! Spectral output, for exoplanet observational comparison
374	if(specOLR)then
375	do nw=1,L_NSPECTI
376	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
377	end do
378	do nw=1,L_NSPECTV
379	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
380	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
381	end do
382	endif
383
384	! Finally, the heating rates
385
386	DO l=2,L_NLAYRAD
387	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
388	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
389	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
390	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
391	END DO
392
393	! These are values at top of atmosphere
394	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
395	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
396	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
397	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
398
399
400	!-----------------------------------------------------------------------
401	end do ! End of big loop over every GCM column.
402	!-----------------------------------------------------------------------
403
404
405
406	!-----------------------------------------------------------------------
407	! Additional diagnostics
408	!-----------------------------------------------------------------------
409
410	! IR spectral output, for exoplanet observational comparison
411	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
412
413	print*,'Saving scalar quantities in surf_vals.out...'
414	print*,'psurf = ', pplev(1,1),' Pa'
415	open(116,file='surf_vals.out')
416	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
417	real(-nfluxtopv),real(nfluxtopi)
418	close(116)
419
420
421	! USEFUL COMMENT - Do Not Remove.
422	!
423	! if(specOLR)then
424	! open(117,file='OLRnu.out')
425	! do nw=1,L_NSPECTI
426	! write(117,*) OLR_nu(1,nw)
427	! enddo
428	! close(117)
429	!
430	! open(127,file='OSRnu.out')
431	! do nw=1,L_NSPECTV
432	! write(127,*) OSR_nu(1,nw)
433	! enddo
434	! close(127)
435	! endif
436
437	! OLR vs altitude: do it as a .txt file.
438	OLRz=.false.
439	if(OLRz)then
440	print*,'saving IR vertical flux for OLRz...'
441	open(118,file='OLRz_plevs.out')
442	open(119,file='OLRz.out')
443	do l=1,L_NLAYRAD
444	write(118,) plevrad(2l)
445	do nw=1,L_NSPECTI
446	write(119,*) fluxupi_nu(l,nw)
447	enddo
448	enddo
449	close(118)
450	close(119)
451	endif
452
453	endif
454
455	! See physiq.F for explanations about CLFvarying. This is temporary.
456	if (lastcall) then
457	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
458	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
459	!$OMP BARRIER
460	!$OMP MASTER
461	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
462	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
463	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
464	!$OMP END MASTER
465	!$OMP BARRIER
466	endif
467
468
469	end subroutine callcorrk

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

Download in other formats:

Original Format