Context Navigation

callcorrk.F90 @ 1680

Last change on this file since 1680 was 1648, checked in by jvatant, 8 years ago

Modifications to custom radiative transfer to Titan
+ Enables an altitude dependant gfrac for CIA computations

-> many radical changes in su_gases and co ..
-> read vertical CH4 profile with call_profilgases
-> Now you need a 'profile.def' that I will add in the deftank

+ Added interpolate CIA routines for CH4
+ Added temporary mean aerosol profile opacity routine (disr_haze)

Property svn:executable set to *

File size: 31.0 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
3	tsurf,fract,dist_star,aerosol, &
4	dtlw,dtsw,fluxsurf_lw, &
5	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
6	fluxabs_sw,fluxtop_dn, &
7	OLR_nu,OSR_nu, &
8	tau_col,firstcall,lastcall)
9
10	use radinc_h
11	use radcommon_h
12	use datafile_mod, only: datadir
13	use ioipsl_getin_p_mod, only: getin_p
14	use gases_h
15	use radii_mod, only : su_aer_radii,back2lay_reffrad
16	use aerosol_mod, only : iaero_back2lay
17	USE tracer_h
18	use comcstfi_mod, only: pi, mugaz, cpp
19	use callkeys_mod, only: diurnal,tracer,nosurf, &
20	strictboundcorrk,specOLR
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	!
37	!==================================================================
38
39	!-----------------------------------------------------------------------
40	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
41	! Layer #1 is the layer near the ground.
42	! Layer #nlayer is the layer at the top.
43	!-----------------------------------------------------------------------
44
45
46	! INPUT
47	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
48	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
49	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
50	INTEGER,INTENT(IN) :: nq ! Number of tracers.
51	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
52	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
53	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
54	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
55	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
56	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
57	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
58	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
59	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
60	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
61	logical,intent(in) :: firstcall ! Signals first call to physics.
62	logical,intent(in) :: lastcall ! Signals last call to physics.
63
64	! OUTPUT
65	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
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) :: tau_col(ngrid) ! Diagnostic from aeropacity.
77	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
78
79
80	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
81	REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
82	REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
83	REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
84	REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
85	REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
86	REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
87
88	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
89	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
90
91	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
92	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
93	!$OMP THREADPRIVATE(reffrad,nueffrad)
94
95	!-----------------------------------------------------------------------
96	! Declaration of the variables required by correlated-k subroutines
97	! Numbered from top to bottom (unlike in the GCM)
98	!-----------------------------------------------------------------------
99
100	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
101	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
102
103	! Optical values for the optci/cv subroutines
104	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
105	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
106	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
107	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
108	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
109	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
110	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
111	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
112	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
113	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
114
115	REAL*8 tauaero(L_LEVELS+1,naerkind)
116	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
117	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
118	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
119	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
120	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
121	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
122	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
123	REAL*8 albi,acosz
124	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
125
126	INTEGER ig,l,k,nw,iaer
127
128	real szangle
129	logical global1d
130	save szangle,global1d
131	!$OMP THREADPRIVATE(szangle,global1d)
132	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
133	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
134
135	! Local aerosol optical properties for each column on RADIATIVE grid.
136	real*8,save,allocatable :: QXVAER(:,:,:)
137	real*8,save,allocatable :: QSVAER(:,:,:)
138	real*8,save,allocatable :: GVAER(:,:,:)
139	real*8,save,allocatable :: QXIAER(:,:,:)
140	real*8,save,allocatable :: QSIAER(:,:,:)
141	real*8,save,allocatable :: GIAER(:,:,:)
142
143	real, dimension(:,:,:), save, allocatable :: QREFvis3d
144	real, dimension(:,:,:), save, allocatable :: QREFir3d
145	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
146
147
148	! Miscellaneous :
149	real*8 temp,temp1,temp2,pweight
150	character(len=10) :: tmp1
151	character(len=10) :: tmp2
152
153	logical OLRz
154	real*8 NFLUXGNDV_nu(L_NSPECTV)
155
156	! Included by MT for albedo calculations.
157	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
158	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
159
160
161	!===============================================================
162	! I.a Initialization on first call
163	!===============================================================
164
165
166	if(firstcall) then
167
168	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
169	if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
170	if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
171	if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
172	if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
173	if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
174	if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
175
176	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
177	IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
178	IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
179	! Effective radius and variance of the aerosols
180	IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
181	IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
182
183	call system('rm -f surf_vals_long.out')
184
185	if(naerkind.gt.4)then
186	print*,'Code not general enough to deal with naerkind > 4 yet.'
187	call abort
188	endif
189	call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
190
191
192	!--------------------------------------------------
193	! Set up correlated k
194	!--------------------------------------------------
195
196
197	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
198	call getin_p("corrkdir",corrkdir)
199	print*, "corrkdir = ",corrkdir
200	write( tmp1, '(i3)' ) L_NSPECTI
201	write( tmp2, '(i3)' ) L_NSPECTV
202	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
203	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
204
205	call setspi ! Basic infrared properties.
206	call setspv ! Basic visible properties.
207	call sugas_corrk ! Set up gaseous absorption properties.
208	call suaer_corrk ! Set up aerosol optical properties.
209
210
211	OLR_nu(:,:) = 0.
212	OSR_nu(:,:) = 0.
213
214	if (ngrid.eq.1) then
215	PRINT*, 'Simulate global averaged conditions ?'
216	global1d = .false. ! default value
217	call getin_p("global1d",global1d)
218	write(,) "global1d = ",global1d
219
220	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
221	if (global1d.and.diurnal) then
222	print*,'if global1d is true, diurnal must be set to false'
223	stop
224	endif
225
226	if (global1d) then
227	PRINT *,'Solar Zenith angle (deg.) ?'
228	PRINT *,'(assumed for averaged solar flux S/4)'
229	szangle=60.0 ! default value
230	call getin_p("szangle",szangle)
231	write(,) "szangle = ",szangle
232	endif
233	endif
234
235	end if ! of if (firstcall)
236
237	!=======================================================================
238	! I.b Initialization on every call
239	!=======================================================================
240
241	qxvaer(:,:,:)=0.0
242	qsvaer(:,:,:)=0.0
243	gvaer(:,:,:) =0.0
244
245	qxiaer(:,:,:)=0.0
246	qsiaer(:,:,:)=0.0
247	giaer(:,:,:) =0.0
248
249	!--------------------------------------------------
250	! Effective radius and variance of the aerosols
251	!--------------------------------------------------
252
253	do iaer=1,naerkind
254
255	if(iaer.eq.iaero_back2lay)then
256	call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
257	endif
258
259	end do !iaer=1,naerkind.
260
261
262	! How much light do we get ?
263	do nw=1,L_NSPECTV
264	stel(nw)=stellarf(nw)/(dist_star**2)
265	end do
266
267	! Get 3D aerosol optical properties.
268	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
269	QVISsQREF3d,omegaVIS3d,gVIS3d, &
270	QIRsQREF3d,omegaIR3d,gIR3d, &
271	QREFvis3d,QREFir3d)
272
273	! Get aerosol optical depths.
274	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
275	reffrad,QREFvis3d,QREFir3d, &
276	tau_col)
277
278
279
280	!-----------------------------------------------------------------------
281	do ig=1,ngrid ! Starting Big Loop over every GCM column
282	!-----------------------------------------------------------------------
283
284
285	!=======================================================================
286	! II. Transformation of the GCM variables
287	!=======================================================================
288
289
290	!-----------------------------------------------------------------------
291	! Aerosol optical properties Qext, Qscat and g.
292	! The transformation in the vertical is the same as for temperature.
293	!-----------------------------------------------------------------------
294
295
296	do iaer=1,naerkind
297	! Shortwave.
298	do nw=1,L_NSPECTV
299
300	do l=1,nlayer
301
302	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
303	*QREFvis3d(ig,nlayer+1-l,iaer)
304
305	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
306	*QREFvis3d(ig,max(nlayer-l,1),iaer)
307
308	qxvaer(2*l,nw,iaer) = temp1
309	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
310
311	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
312	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
313
314	qsvaer(2*l,nw,iaer) = temp1
315	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
316
317	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
318	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
319
320	gvaer(2*l,nw,iaer) = temp1
321	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
322
323	end do ! nlayer
324
325	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
326	qxvaer(2*nlayer+1,nw,iaer)=0.
327
328	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
329	qsvaer(2*nlayer+1,nw,iaer)=0.
330
331	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
332	gvaer(2*nlayer+1,nw,iaer)=0.
333
334	end do ! L_NSPECTV
335
336	do nw=1,L_NSPECTI
337	! Longwave
338	do l=1,nlayer
339
340	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
341	*QREFir3d(ig,nlayer+1-l,iaer)
342
343	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
344	*QREFir3d(ig,max(nlayer-l,1),iaer)
345
346	qxiaer(2*l,nw,iaer) = temp1
347	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
348
349	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
350	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
351
352	qsiaer(2*l,nw,iaer) = temp1
353	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
354
355	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
356	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
357
358	giaer(2*l,nw,iaer) = temp1
359	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
360
361	end do ! nlayer
362
363	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
364	qxiaer(2*nlayer+1,nw,iaer)=0.
365
366	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
367	qsiaer(2*nlayer+1,nw,iaer)=0.
368
369	giaer(1,nw,iaer)=giaer(2,nw,iaer)
370	giaer(2*nlayer+1,nw,iaer)=0.
371
372	end do ! L_NSPECTI
373
374	end do ! naerkind
375
376	! Test / Correct for freaky s. s. albedo values.
377	do iaer=1,naerkind
378	do k=1,L_LEVELS+1
379
380	do nw=1,L_NSPECTV
381	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
382	print*,'Serious problems with qsvaer values'
383	print*,'in callcorrk'
384	call abort
385	endif
386	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
387	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
388	endif
389	end do
390
391	do nw=1,L_NSPECTI
392	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
393	print*,'Serious problems with qsiaer values'
394	print*,'in callcorrk'
395	call abort
396	endif
397	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
398	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
399	endif
400	end do
401
402	end do ! L_LEVELS
403	end do ! naerkind
404
405	!-----------------------------------------------------------------------
406	! Aerosol optical depths
407	!-----------------------------------------------------------------------
408
409	do iaer=1,naerkind ! a bug was here
410	do k=0,nlayer-1
411
412	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
413	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
414	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
415	tauaero(2k+2,iaer)=max(temppweight,0.d0)
416	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
417
418	end do
419	! boundary conditions
420	tauaero(1,iaer) = tauaero(2,iaer)
421	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
422	!tauaero(1,iaer) = 0.
423	!tauaero(L_LEVELS+1,iaer) = 0.
424
425	end do ! naerkind
426
427	! Albedo and Emissivity.
428	albi=1-emis(ig) ! Long Wave.
429	DO nw=1,L_NSPECTV ! Short Wave loop.
430	albv(nw)=albedo(ig,nw)
431	ENDDO
432
433	if (nosurf) then ! Case with no surface.
434	DO nw=1,L_NSPECTV
435	if(albv(nw).gt.0.0) then
436	print*,'For open lower boundary in callcorrk must'
437	print*,'have spectral surface band albedos all set to zero!'
438	call abort
439	endif
440	ENDDO
441	endif
442
443	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
444	acosz = cos(pi*szangle/180.0)
445	print*,'acosz=',acosz,', szangle=',szangle
446	else
447	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
448	endif
449
450	!-----------------------------------------------------------------------
451	! Pressure and temperature
452	!-----------------------------------------------------------------------
453
454	DO l=1,nlayer
455	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
456	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
457	tlevrad(2*l) = pt(ig,nlayer+1-l)
458	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
459	END DO
460
461	plevrad(1) = 0.
462	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.
463
464	tlevrad(1) = tlevrad(2)
465	tlevrad(2*nlayer+1)=tsurf(ig)
466
467	pmid(1) = max(pgasmin,0.0001*plevrad(3))
468	pmid(2) = pmid(1)
469
470	tmid(1) = tlevrad(2)
471	tmid(2) = tmid(1)
472
473	DO l=1,L_NLAYRAD-1
474	tmid(2l+1) = tlevrad(2l+1)
475	tmid(2l+2) = tlevrad(2l+1)
476	pmid(2l+1) = plevrad(2l+1)
477	pmid(2l+2) = plevrad(2l+1)
478	END DO
479	pmid(L_LEVELS) = plevrad(L_LEVELS)
480	tmid(L_LEVELS) = tlevrad(L_LEVELS)
481
482	!!Alternative interpolation:
483	! pmid(3) = pmid(1)
484	! pmid(4) = pmid(1)
485	! tmid(3) = tmid(1)
486	! tmid(4) = tmid(1)
487	! DO l=2,L_NLAYRAD-1
488	! tmid(2l+1) = tlevrad(2l)
489	! tmid(2l+2) = tlevrad(2l)
490	! pmid(2l+1) = plevrad(2l)
491	! pmid(2l+2) = plevrad(2l)
492	! END DO
493	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
494	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
495
496	! Test for out-of-bounds pressure.
497	if(plevrad(3).lt.pgasmin)then
498	print*,'Minimum pressure is outside the radiative'
499	print*,'transfer kmatrix bounds, exiting.'
500	call abort
501	elseif(plevrad(L_LEVELS).gt.pgasmax)then
502	print*,'Maximum pressure is outside the radiative'
503	print*,'transfer kmatrix bounds, exiting.'
504	call abort
505	endif
506
507	! Test for out-of-bounds temperature.
508	do k=1,L_LEVELS
509	if(tlevrad(k).lt.tgasmin)then
510	print*,'Minimum temperature is outside the radiative'
511	print*,'transfer kmatrix bounds'
512	print*,"k=",k," tlevrad(k)=",tlevrad(k)
513	print*,"tgasmin=",tgasmin
514	if (strictboundcorrk) then
515	call abort
516	else
517	print,'**********************************************'
518	print*,'we allow model to continue with tlevrad=tgasmin'
519	print*,' ... we assume we know what you are doing ... '
520	print*,' ... but do not let this happen too often ... '
521	print,'**********************************************'
522	!tlevrad(k)=tgasmin
523	endif
524	elseif(tlevrad(k).gt.tgasmax)then
525	! print*,'Maximum temperature is outside the radiative'
526	! print*,'transfer kmatrix bounds, exiting.'
527	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
528	! print*,"tgasmax=",tgasmax
529	if (strictboundcorrk) then
530	call abort
531	else
532	! print,'**********************************************'
533	! print*,'we allow model to continue with tlevrad=tgasmax'
534	! print*,' ... we assume we know what you are doing ... '
535	! print*,' ... but do not let this happen too often ... '
536	! print,'**********************************************'
537	!tlevrad(k)=tgasmax
538	endif
539	endif
540	enddo
541	do k=1,L_NLAYRAD+1
542	if(tmid(k).lt.tgasmin)then
543	print*,'Minimum temperature is outside the radiative'
544	print*,'transfer kmatrix bounds, exiting.'
545	print*,"k=",k," tmid(k)=",tmid(k)
546	print*,"tgasmin=",tgasmin
547	if (strictboundcorrk) then
548	call abort
549	else
550	print,'**********************************************'
551	print*,'we allow model to continue with tmid=tgasmin'
552	print*,' ... we assume we know what you are doing ... '
553	print*,' ... but do not let this happen too often ... '
554	print,'**********************************************'
555	tmid(k)=tgasmin
556	endif
557	elseif(tmid(k).gt.tgasmax)then
558	! print*,'Maximum temperature is outside the radiative'
559	! print*,'transfer kmatrix bounds, exiting.'
560	! print*,"k=",k," tmid(k)=",tmid(k)
561	! print*,"tgasmax=",tgasmax
562	if (strictboundcorrk) then
563	call abort
564	else
565	! print,'**********************************************'
566	! print*,'we allow model to continue with tmid=tgasmin'
567	! print*,' ... we assume we know what you are doing ... '
568	! print*,' ... but do not let this happen too often ... '
569	! print,'**********************************************'
570	tmid(k)=tgasmax
571	endif
572	endif
573	enddo
574
575	!=======================================================================
576	! III. Calling the main radiative transfer subroutines
577	!=======================================================================
578
579
580	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
581	glat_ig=glat(ig)
582
583	!-----------------------------------------------------------------------
584	! Short Wave Part
585	!-----------------------------------------------------------------------
586
587	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
588	if((ngrid.eq.1).and.(global1d))then
589	do nw=1,L_NSPECTV
590	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
591	end do
592	else
593	do nw=1,L_NSPECTV
594	stel_fract(nw)= stel(nw) * fract(ig)
595	end do
596	endif
597
598	call optcv(dtauv,tauv,taucumv,plevrad, &
599	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
600	tmid,pmid,taugsurf,gweight)
601
602	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
603	acosz,stel_fract,gweight, &
604	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
605	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
606
607	else ! During the night, fluxes = 0.
608	nfluxtopv = 0.0d0
609	fluxtopvdn = 0.0d0
610	nfluxoutv_nu(:) = 0.0d0
611	nfluxgndv_nu(:) = 0.0d0
612	do l=1,L_NLAYRAD
613	fmnetv(l)=0.0d0
614	fluxupv(l)=0.0d0
615	fluxdnv(l)=0.0d0
616	end do
617	end if
618
619
620	! Equivalent Albedo Calculation (for OUTPUT). MT2015
621	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
622	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
623	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
624	DO nw=1,L_NSPECTV
625	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
626	ENDDO
627	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
628	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
629	else
630	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
631	endif
632	else
633	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
634	endif
635
636
637	!-----------------------------------------------------------------------
638	! Long Wave Part
639	!-----------------------------------------------------------------------
640
641	call optci(plevrad,tlevrad,dtaui,taucumi, &
642	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
643	taugsurfi,gweight)
644
645	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
646	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
647	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
648
649	!-----------------------------------------------------------------------
650	! Transformation of the correlated-k code outputs
651	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
652
653	! Flux incident at the top of the atmosphere
654	fluxtop_dn(ig)=fluxtopvdn
655
656	fluxtop_lw(ig) = real(nfluxtopi)
657	fluxabs_sw(ig) = real(-nfluxtopv)
658	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
659	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
660
661	! Flux absorbed by the surface. By MT2015.
662	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
663
664	if(fluxtop_dn(ig).lt.0.0)then
665	print*,'Achtung! fluxtop_dn has lost the plot!'
666	print*,'fluxtop_dn=',fluxtop_dn(ig)
667	print*,'acosz=',acosz
668	print*,'aerosol=',aerosol(ig,:,:)
669	print*,'temp= ',pt(ig,:)
670	print*,'pplay= ',pplay(ig,:)
671	call abort
672	endif
673
674	! Spectral output, for exoplanet observational comparison
675	if(specOLR)then
676	do nw=1,L_NSPECTI
677	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
678	end do
679	do nw=1,L_NSPECTV
680	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
681	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
682	end do
683	endif
684
685	! Finally, the heating rates
686
687	DO l=2,L_NLAYRAD
688	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
689	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
690	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
691	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
692	END DO
693
694	! These are values at top of atmosphere
695	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
696	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
697	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
698	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
699
700
701	!-----------------------------------------------------------------------
702	end do ! End of big loop over every GCM column.
703	!-----------------------------------------------------------------------
704
705
706
707	!-----------------------------------------------------------------------
708	! Additional diagnostics
709	!-----------------------------------------------------------------------
710
711	! IR spectral output, for exoplanet observational comparison
712	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
713
714	print*,'Saving scalar quantities in surf_vals.out...'
715	print*,'psurf = ', pplev(1,1),' Pa'
716	open(116,file='surf_vals.out')
717	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
718	real(-nfluxtopv),real(nfluxtopi)
719	close(116)
720
721
722	! USEFUL COMMENT - Do Not Remove.
723	!
724	! if(specOLR)then
725	! open(117,file='OLRnu.out')
726	! do nw=1,L_NSPECTI
727	! write(117,*) OLR_nu(1,nw)
728	! enddo
729	! close(117)
730	!
731	! open(127,file='OSRnu.out')
732	! do nw=1,L_NSPECTV
733	! write(127,*) OSR_nu(1,nw)
734	! enddo
735	! close(127)
736	! endif
737
738	! OLR vs altitude: do it as a .txt file.
739	OLRz=.false.
740	if(OLRz)then
741	print*,'saving IR vertical flux for OLRz...'
742	open(118,file='OLRz_plevs.out')
743	open(119,file='OLRz.out')
744	do l=1,L_NLAYRAD
745	write(118,) plevrad(2l)
746	do nw=1,L_NSPECTI
747	write(119,*) fluxupi_nu(l,nw)
748	enddo
749	enddo
750	close(118)
751	close(119)
752	endif
753
754	endif
755
756	! See physiq.F for explanations about CLFvarying. This is temporary.
757	if (lastcall) then
758	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
759	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
760	!$OMP BARRIER
761	!$OMP MASTER
762	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
763	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
764	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
765	!$OMP END MASTER
766	!$OMP BARRIER
767	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
768	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
769	endif
770
771
772	end subroutine callcorrk

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

Download in other formats:

Original Format