Context Navigation

callcorrk.F90 @ 724

Last change on this file since 724 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

Property svn:executable set to *

File size: 32.0 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
2	albedo,emis,mu0,pplev,pplay,pt, &
3	tsurf,fract,dist_star,aerosol,muvar, &
4	dtlw,dtsw,fluxsurf_lw, &
5	fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn, &
6	OLR_nu,OSR_nu, &
7	reffrad,tau_col,cloudfrac,totcloudfrac, &
8	clearsky,firstcall,lastcall)
9
10	use radinc_h
11	use radcommon_h
12	use watercommon_h
13	use datafile_mod, only: datadir
14	use ioipsl_getincom
15	use gases_h
16
17	implicit none
18
19	!==================================================================
20	!
21	! Purpose
22	! -------
23	! Solve the radiative transfer using the correlated-k method for
24	! the gaseous absorption and the Toon et al. (1989) method for
25	! scatttering due to aerosols.
26	!
27	! Authors
28	! -------
29	! Emmanuel 01/2001, Forget 09/2001
30	! Robin Wordsworth (2009)
31	!
32	!==================================================================
33
34	#include "dimphys.h"
35	#include "comcstfi.h"
36	#include "callkeys.h"
37	#include "tracer.h"
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 #nlayermx is the layer at the top.
43
44	! INPUT
45	INTEGER icount
46	INTEGER ngrid,nlayer
47	REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
48	REAL albedo(ngrid) ! SW albedo
49	REAL emis(ngrid) ! LW emissivity
50	REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer
51	REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries
52
53	REAL pt(ngrid,nlayermx) ! air temperature (K)
54	REAL tsurf(ngrid) ! surface temperature (K)
55	REAL dist_star,mu0(ngrid) ! distance star-planet (AU)
56	REAL fract(ngrid) ! fraction of day
57
58	! Globally varying aerosol optical properties on GCM grid
59	! Not needed everywhere so not in radcommon_h
60	REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
61	REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
62	REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
63
64	REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
65	REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
66	REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
67
68	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
69	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
70
71	! REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
72	! REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
73
74	REAL reffrad(ngrid,nlayer,naerkind)
75	REAL nueffrad(ngrid,nlayer,naerkind)
76
77	! OUTPUT
78	REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
79	REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
80	REAL fluxsurf_lw(ngridmx) ! incident LW flux to surf (W/m2)
81	REAL fluxtop_lw(ngridmx) ! outgoing LW flux to space (W/m2)
82	REAL fluxsurf_sw(ngridmx) ! incident SW flux to surf (W/m2)
83	REAL fluxabs_sw(ngridmx) ! SW flux absorbed by planet (W/m2)
84	REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2)
85	REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
86	REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
87	!-----------------------------------------------------------------------
88	! Declaration of the variables required by correlated-k subroutines
89	! Numbered from top to bottom unlike in the GCM!
90
91	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
92	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
93
94	! Optical values for the optci/cv subroutines
95	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
96	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
97	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
98	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
99	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
100	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
101	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
102	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
103	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
104	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
105
106	REAL*8 tauaero(L_LEVELS+1,naerkind)
107	REAL*8 nfluxtopv,nfluxtopi,nfluxtop
108	real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
109	real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
110	real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
111	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
112	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
113	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
114	REAL*8 albi,albv,acosz
115
116	INTEGER ig,l,k,nw,iaer,irad
117
118	real fluxtoplanet
119	real szangle
120	logical global1d
121	save szangle,global1d
122	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
123	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
124
125	real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol)
126	REAL pq(ngridmx,nlayer,nq)
127	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
128	integer nq
129
130	! Local aerosol optical properties for each column on RADIATIVE grid
131	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
132	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
133	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
134	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
135	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
136	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
137
138	save qxvaer, qsvaer, gvaer
139	save qxiaer, qsiaer, giaer
140	save QREFvis3d, QREFir3d
141
142	REAL tau_col(ngrid) ! diagnostic from aeropacity
143
144	! Misc.
145	logical firstcall, lastcall, nantest
146	real*8 tempv(L_NSPECTV)
147	real*8 tempi(L_NSPECTI)
148	real*8 temp,temp1,temp2,pweight
149	character(len=10) :: tmp1
150	character(len=10) :: tmp2
151
152	! for fixed water vapour profiles
153	integer i_var
154	real RH
155	real*8 pq_temp(nlayer)
156	real ptemp, Ttemp, qsat
157
158	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
159
160	!real ptime, pday
161	logical OLRz
162	real*8 NFLUXGNDV_nu(L_NSPECTV)
163
164	! for H2O cloud fraction in aeropacity
165	real cloudfrac(ngridmx,nlayermx)
166	real totcloudfrac(ngridmx)
167	logical clearsky
168
169	! for weird cloud test
170	real pqtest(ngridmx,nlayer,nq)
171
172	! are particle radii fixed?
173	logical, save :: radfixed
174	real maxrad, minrad
175
176	! Local water cloud optical properties
177	real, parameter :: rad_froid=35.0e-6
178	real, parameter :: rad_chaud=10.0e-6
179	real, parameter :: coef_chaud=0.13
180	real, parameter :: coef_froid=0.09
181	real zfice
182
183	real CBRT
184	external CBRT
185
186	! included by RW for runaway greenhouse 1D study
187	real muvar(ngridmx,nlayermx+1)
188	real vtmp(nlayermx)
189	REAL*8 muvarrad(L_LEVELS)
190
191	!=======================================================================
192	! Initialization on first call
193
194	qxvaer(:,:,:)=0.0
195	qsvaer(:,:,:)=0.0
196	gvaer(:,:,:) =0.0
197
198	qxiaer(:,:,:)=0.0
199	qsiaer(:,:,:)=0.0
200	giaer(:,:,:) =0.0
201	radfixed=.false.
202
203	if(firstcall) then
204	if(kastprof)then
205	radfixed=.true.
206	endif
207
208	call system('rm -f surf_vals_long.out')
209
210	!--------------------------------------------------
211	! Effective radius and variance of the aerosols
212
213	do iaer=1,naerkind
214	! these values will change once the microphysics gets to work
215	! UNLESS tracer=.false., in which case we should be working with
216	! a fixed aerosol layer, and be able to define reffrad in a
217	! .def file. To be improved!
218
219	if(iaer.eq.1)then ! CO2 ice
220	do l=1,nlayer
221	do ig=1,ngrid
222	reffrad(ig,l,iaer) = 1.e-4
223	nueffrad(ig,l,iaer) = 0.1
224	enddo
225	enddo
226	endif
227
228	if(iaer.eq.2)then ! H2O ice
229	do l=1,nlayer
230	do ig=1,ngrid
231	reffrad(ig,l,iaer) = 1.e-5
232	nueffrad(ig,l,iaer) = 0.1
233	enddo
234	enddo
235	endif
236
237	if(iaer.eq.3)then ! dust
238	do l=1,nlayer
239	do ig=1,ngrid
240	reffrad(ig,l,iaer) = 1.e-5
241	nueffrad(ig,l,iaer) = 0.1
242	enddo
243	enddo
244	endif
245
246	if(iaer.gt.3)then
247	print*,'Error in callcorrk, naerkind is too high.'
248	print*,'The code still needs generalisation to arbitrary'
249	print*,'aerosol kinds and number.'
250	call abort
251	endif
252
253	enddo
254
255	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
256	call getin("corrkdir",corrkdir)
257	print*, "corrkdir = ",corrkdir
258	write( tmp1, '(i3)' ) L_NSPECTI
259	write( tmp2, '(i3)' ) L_NSPECTV
260	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
261	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
262
263	call sugas_corrk ! set up gaseous absorption properties
264	call setspi ! basic infrared properties
265	call setspv ! basic visible properties
266	call suaer_corrk ! set up aerosol optical properties
267
268	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
269
270	if((igcm_h2o_vap.eq.0) .and. varactive)then
271	print*,'varactive in callcorrk but no h2o_vap tracer.'
272	stop
273	endif
274
275	OLR_nu(:,:) = 0.
276	OSR_nu(:,:) = 0.
277
278	if (ngridmx.eq.1) then
279	PRINT*, 'Simulate global averaged conditions ?'
280	global1d = .false. ! default value
281	call getin("global1d",global1d)
282	write(,) "global1d = ",global1d
283	! Test of incompatibility:
284	! if global1d is true, there should not be any diurnal cycle
285	if (global1d.and.diurnal) then
286	print*,'if global1d is true, diurnal must be set to false'
287	stop
288	endif
289
290	if (global1d) then
291	PRINT *,'Solar Zenith angle (deg.) ?'
292	PRINT *,'(assumed for averaged solar flux S/4)'
293	szangle=60.0 ! default value
294	call getin("szangle",szangle)
295	write(,) "szangle = ",szangle
296	endif
297	endif
298
299	firstcall=.false.
300
301	end if
302
303	!=======================================================================
304	! Initialization on every call
305
306	do l=1,nlayer
307	do ig=1,ngrid
308	do iaer=1,naerkind
309	nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
310	enddo
311	enddo
312	enddo
313
314
315
316	if(radfixed)then
317	do l=1,nlayer
318	do ig=1,ngrid
319	reffrad(ig,l,1) = 5.e-5 ! CO2 ice
320	enddo
321	enddo
322	print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
323	if(naerkind.ge.2)then
324	do l=1,nlayer
325	do ig=1,ngrid
326	!reffrad(ig,l,2) = 2.e-5 ! H2O ice
327	reffrad(ig,l,2) = rad_chaud ! H2O ice
328
329	zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
330	zfice = MIN(MAX(zfice,0.0),1.0)
331	reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
332	nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
333	enddo
334	enddo
335	print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
336	endif
337	if(naerkind.eq.3)then
338	do l=1,nlayer
339	do ig=1,ngrid
340	reffrad(ig,l,naerkind) = 2.e-6 ! dust
341	enddo
342	enddo
343	print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
344	endif
345	if(naerkind.gt.3)then
346	print*,'Code not general enough to deal with naerkind > 3 yet.'
347	call abort
348	endif
349
350	else
351
352	maxrad=0.0
353	!averad=0.0
354	minrad=1.0
355	do l=1,nlayer
356
357	!masse = (pplev(ig,l) - pplev(ig,l+1))/g
358
359	do ig=1,ngrid
360	!if(tracer)then
361	if(tracer.and.igcm_co2_ice.gt.0)then
362
363	if(igcm_co2_ice.lt.1)then
364	print*,'Tracers but no CO2 ice still seems to be a problem...'
365	print*,'Aborting in callcorrk.'
366	stop
367	endif
368
369	reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
370	(4Nmix_co2pi*rho_co2) )
371	endif
372	reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
373	reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
374
375	!averad = averad + reffrad(ig,l,1)*area(ig)
376	maxrad = max(reffrad(ig,l,1),maxrad)
377	minrad = min(reffrad(ig,l,1),minrad)
378	enddo
379	enddo
380	if(igcm_co2_ice.gt.0)then
381	print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
382	print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
383	endif
384
385	if((naerkind.ge.2).and.water)then
386	maxrad=0.0
387	minrad=1.0
388	do l=1,nlayer
389	do ig=1,ngrid
390	reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
391	(4Nmix_h2opi*rho_ice) )
392	reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
393	reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)
394
395	maxrad = max(reffrad(ig,l,2),maxrad)
396	minrad = min(reffrad(ig,l,2),minrad)
397	enddo
398	enddo
399	print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
400	print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
401	endif
402
403	if(naerkind.eq.3)then
404	do l=1,nlayer
405	do ig=1,ngrid
406	reffrad(ig,l,naerkind) = 2.e-6 ! dust
407	enddo
408	enddo
409	endif
410
411	endif
412
413
414	! how much light we get
415	fluxtoplanet=0
416	do nw=1,L_NSPECTV
417	stel(nw)=stellarf(nw)/(dist_star**2)
418	fluxtoplanet=fluxtoplanet + stel(nw)
419	end do
420
421	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
422	QVISsQREF3d,omegaVIS3d,gVIS3d, &
423	QIRsQREF3d,omegaIR3d,gIR3d, &
424	QREFvis3d,QREFir3d) ! get 3D aerosol optical properties
425
426	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
427	reffrad,QREFvis3d,QREFir3d, &
428	tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths
429
430	!-----------------------------------------------------------------------
431	! Starting Big Loop over every GCM column
432	do ig=1,ngridmx
433
434	!=======================================================================
435	! Transformation of the GCM variables
436
437	!-----------------------------------------------------------------------
438	! Aerosol optical properties Qext, Qscat and g
439	! The transformation in the vertical is the same as for temperature
440
441	! shortwave
442	do iaer=1,naerkind
443	DO nw=1,L_NSPECTV
444	do l=1,nlayermx
445
446	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) &
447	*QREFvis3d(ig,nlayermx+1-l,iaer)
448
449	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
450	*QREFvis3d(ig,max(nlayermx-l,1),iaer)
451
452	qxvaer(2*l,nw,iaer) = temp1
453	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
454
455	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
456	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
457
458	qsvaer(2*l,nw,iaer) = temp1
459	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
460
461	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
462	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
463
464	gvaer(2*l,nw,iaer) = temp1
465	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
466
467	end do
468
469	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
470	qxvaer(2*nlayermx+1,nw,iaer)=0.
471
472	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
473	qsvaer(2*nlayermx+1,nw,iaer)=0.
474
475	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
476	gvaer(2*nlayermx+1,nw,iaer)=0.
477
478	end do
479
480	! longwave
481	DO nw=1,L_NSPECTI
482	do l=1,nlayermx
483
484	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) &
485	*QREFir3d(ig,nlayermx+1-l,iaer)
486
487	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
488	*QREFir3d(ig,max(nlayermx-l,1),iaer)
489
490	qxiaer(2*l,nw,iaer) = temp1
491	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
492
493	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
494	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
495
496	qsiaer(2*l,nw,iaer) = temp1
497	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
498
499	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
500	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
501
502	giaer(2*l,nw,iaer) = temp1
503	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
504
505	end do
506
507	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
508	qxiaer(2*nlayermx+1,nw,iaer)=0.
509
510	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
511	qsiaer(2*nlayermx+1,nw,iaer)=0.
512
513	giaer(1,nw,iaer)=giaer(2,nw,iaer)
514	giaer(2*nlayermx+1,nw,iaer)=0.
515
516	end do
517	end do
518
519	! test / correct for freaky s. s. albedo values
520	do iaer=1,naerkind
521	do k=1,L_LEVELS+1
522
523	do nw=1,L_NSPECTV
524	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
525	print*,'Serious problems with qsvaer values in callcorrk'
526	call abort
527	endif
528	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
529	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
530	endif
531	end do
532
533	do nw=1,L_NSPECTI
534	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
535	print*,'Serious problems with qsiaer values in callcorrk'
536	call abort
537	endif
538	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
539	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
540	endif
541	end do
542
543	end do
544	end do
545
546	!-----------------------------------------------------------------------
547	! Aerosol optical depths
548
549	do iaer=1,naerkind ! a bug was here
550	do k=0,nlayer-1
551
552	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
553	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
554
555	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
556
557	tauaero(2k+2,iaer)=max(temppweight,0.d0)
558	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
559	!
560	end do
561	! boundary conditions
562	tauaero(1,iaer) = tauaero(2,iaer)
563	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
564	!tauaero(1,iaer) = 0.
565	!tauaero(L_LEVELS+1,iaer) = 0.
566	end do
567
568	! Albedo and emissivity
569	albi=1-emis(ig) ! longwave
570	albv=albedo(ig) ! shortwave
571
572	if(noradsurf.and.(albv.gt.0.0))then
573	print*,'For open lower boundary in callcorrk must'
574	print*,'have surface albedo set to zero!'
575	call abort
576	endif
577
578	if ((ngridmx.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
579	acosz = cos(pi*szangle/180.0)
580	print*,'acosz=',acosz,', szangle=',szangle
581	else
582	acosz=mu0(ig) ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude
583	endif
584
585	!-----------------------------------------------------------------------
586	! Water vapour (to be generalised for other gases eventually)
587
588	if(varactive)then
589
590	i_var=igcm_h2o_vap
591	do l=1,nlayer
592	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
593	qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
594	! Average approximation as for temperature...
595	end do
596	qvar(1)=qvar(2)
597
598	elseif(varfixed)then
599
600	do l=1,nlayermx ! here we will assign fixed water vapour profiles globally
601	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
602	if(RH.lt.0.0) RH=0.0
603
604	ptemp=pplay(ig,l)
605	Ttemp=pt(ig,l)
606	call watersat(Ttemp,ptemp,qsat)
607
608	!pq_temp(l) = qsat ! fully saturated everywhere
609	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
610	end do
611
612	do l=1,nlayer
613	qvar(2*l) = pq_temp(nlayer+1-l)
614	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
615	end do
616	qvar(1)=qvar(2)
617
618	! Lowest layer of atmosphere
619	RH = satval * (1 - 0.02) / 0.98
620	if(RH.lt.0.0) RH=0.0
621
622	ptemp = pplev(ig,1)
623	Ttemp = tsurf(ig)
624	call watersat(Ttemp,ptemp,qsat)
625
626	!qvar(2*nlayermx+1)=qsat ! fully saturated everywhere
627	qvar(2nlayermx+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
628	!qvar=0.005 ! completely constant profile (JL)
629
630	else
631	do k=1,L_LEVELS
632	qvar(k) = 1.0D-7
633	end do
634	end if
635
636	if(.not.kastprof)then
637	! IMPORTANT: Now convert from kg/kg to mol/mol
638	do k=1,L_LEVELS
639	qvar(k) = qvar(k)/epsi
640	end do
641	end if
642
643	!-----------------------------------------------------------------------
644	! kcm mode only
645	if(kastprof)then
646
647	! initial values equivalent to mugaz
648	DO l=1,nlayer
649	muvarrad(2*l) = mugaz
650	muvarrad(2*l+1) = mugaz
651	END DO
652
653	!do k=1,L_LEVELS
654	! qvar(k) = 0.0
655	!end do
656	!print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
657	endif
658
659
660	if(kastprof.and.(ngasmx.gt.1))then
661
662	DO l=1,nlayer
663	muvarrad(2*l) = muvar(ig,nlayer+2-l)
664	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
665	muvar(ig,max(nlayer+1-l,1)))/2
666	END DO
667
668	muvarrad(1) = muvarrad(2)
669	muvarrad(2*nlayermx+1)=muvar(ig,1)
670
671	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
672	print*,'Assumes that the variable gas is H2O!!!'
673	print*,'Assumes that there is only one tracer'
674	!i_var=igcm_h2o_vap
675	i_var=1
676	if(nqmx.gt.1)then
677	print*,'Need 1 tracer only to run kcm1d.e'
678	stop
679	endif
680	do l=1,nlayer
681	vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
682	end do
683
684	do l=1,nlayer
685	qvar(2*l) = vtmp(nlayer+1-l)
686	qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
687	end do
688	qvar(1)=qvar(2)
689
690	print*,'Warning: reducing qvar in callcorrk.'
691	print*,'Temperature profile no longer consistent ', &
692	'with saturated H2O.'
693	do k=1,L_LEVELS
694	qvar(k) = qvar(k)*satval
695	end do
696
697	endif
698
699	! Keep values inside limits for which we have radiative transfer coefficients
700	if(L_REFVAR.gt.1)then ! there was a bug here!
701	do k=1,L_LEVELS
702	if(qvar(k).lt.wrefvar(1))then
703	qvar(k)=wrefvar(1)+1.0e-8
704	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
705	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
706	endif
707	end do
708	endif
709
710	!-----------------------------------------------------------------------
711	! Pressure and temperature
712
713	DO l=1,nlayer
714	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
715	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
716	tlevrad(2*l) = pt(ig,nlayer+1-l)
717	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
718	END DO
719
720	plevrad(1) = 0.
721	plevrad(2) = max(pgasmin,0.0001*plevrad(3))
722
723	tlevrad(1) = tlevrad(2)
724	tlevrad(2*nlayermx+1)=tsurf(ig)
725
726	tmid(1) = tlevrad(2)
727	tmid(2) = tlevrad(2)
728	pmid(1) = plevrad(2)
729	pmid(2) = plevrad(2)
730
731	DO l=1,L_NLAYRAD-1
732	tmid(2l+1) = tlevrad(2l+1)
733	tmid(2l+2) = tlevrad(2l+1)
734	pmid(2l+1) = plevrad(2l+1)
735	pmid(2l+2) = plevrad(2l+1)
736	END DO
737	pmid(L_LEVELS) = plevrad(L_LEVELS)
738	tmid(L_LEVELS) = tlevrad(L_LEVELS)
739
740	! test for out-of-bounds pressure
741	if(plevrad(3).lt.pgasmin)then
742	print*,'Minimum pressure is outside the radiative'
743	print*,'transfer kmatrix bounds, exiting.'
744	call abort
745	elseif(plevrad(L_LEVELS).gt.pgasmax)then
746	print*,'Maximum pressure is outside the radiative'
747	print*,'transfer kmatrix bounds, exiting.'
748	call abort
749	endif
750
751	! test for out-of-bounds temperature
752	do k=1,L_LEVELS
753	if(tlevrad(k).lt.tgasmin)then
754	print*,'Minimum temperature is outside the radiative'
755	print*,'transfer kmatrix bounds, exiting.'
756	!print*,'WARNING, OVERRIDING FOR TEST'
757	call abort
758	elseif(tlevrad(k).gt.tgasmax)then
759	print*,'Maximum temperature is outside the radiative'
760	print*,'transfer kmatrix bounds, exiting.'
761	!print*,'WARNING, OVERRIDING FOR TEST'
762	call abort
763	endif
764	enddo
765
766	!=======================================================================
767	! Calling the main radiative transfer subroutines
768
769
770	!-----------------------------------------------------------------------
771	! Shortwave
772
773	if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
774
775	fluxtoplanet=0.
776
777	if((ngridmx.eq.1).and.(global1d))then
778	do nw=1,L_NSPECTV
779	stel_fract(nw)= stel(nw) * 0.25 / acosz
780	fluxtoplanet=fluxtoplanet + stel_fract(nw)
781	! globally averaged = divide by 4
782	! but we correct for solar zenith angle
783	end do
784	else
785	do nw=1,L_NSPECTV
786	stel_fract(nw)= stel(nw) * fract(ig)
787	fluxtoplanet=fluxtoplanet + stel_fract(nw)
788	end do
789	endif
790
791	call optcv(dtauv,tauv,taucumv,plevrad, &
792	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
793	tmid,pmid,taugsurf,qvar,muvarrad)
794
795	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
796	acosz,stel_fract,gweight, &
797	nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, &
798	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
799
800	else ! during the night, fluxes = 0
801	nfluxtopv = 0.0
802	nfluxoutv_nu(:) = 0.0
803	nfluxgndv_nu(:) = 0.0
804	do l=1,L_NLAYRAD
805	fmnetv(l)=0.0
806	fluxupv(l)=0.0
807	fluxdnv(l)=0.0
808	end do
809	end if
810
811	!-----------------------------------------------------------------------
812	! Longwave
813
814	call optci(plevrad,tlevrad,dtaui,taucumi, &
815	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
816	taugsurfi,qvar,muvarrad)
817
818	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
819	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
820	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
821
822	!-----------------------------------------------------------------------
823	! Transformation of the correlated-k code outputs
824	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
825
826	! Flux incident at the top of the atmosphere
827	fluxtop_dn(ig)=fluxdnv(1)
828
829	fluxtop_lw(ig) = real(nfluxtopi)
830	fluxabs_sw(ig) = real(-nfluxtopv)
831	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
832	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
833
834	if(fluxtop_dn(ig).lt.0.0)then
835	print*,'Achtung! fluxtop_dn has lost the plot!'
836	print*,'fluxtop_dn=',fluxtop_dn(ig)
837	print*,'acosz=',acosz
838	print*,'aerosol=',aerosol(ig,:,:)
839	print*,'temp= ',pt(ig,:)
840	print*,'pplay= ',pplay(ig,:)
841	call abort
842	endif
843
844	! Spectral output, for exoplanet observational comparison
845	if(specOLR)then
846	do nw=1,L_NSPECTI
847	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
848	end do
849	do nw=1,L_NSPECTV
850	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
851	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
852	end do
853	endif
854
855	! Finally, the heating rates
856
857	DO l=2,L_NLAYRAD
858	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
859	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
860	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
861	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
862	END DO
863
864	! These are values at top of atmosphere
865	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
866	g/(cppscalep*(plevrad(3)-plevrad(1)))
867	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
868	g/(cppscalep*(plevrad(3)-plevrad(1)))
869
870	! ---------------------------------------------------------------
871	end do ! end of big loop over every GCM column (ig = 1:ngrid)
872
873
874	!-----------------------------------------------------------------------
875	! Additional diagnostics
876
877	! IR spectral output, for exoplanet observational comparison
878
879
880	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
881
882	print*,'Saving scalar quantities in surf_vals.out...'
883	print*,'psurf = ', pplev(1,1),' Pa'
884	open(116,file='surf_vals.out')
885	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
886	real(-nfluxtopv),real(nfluxtopi)
887	close(116)
888
889	! I am useful, please don`t remove me!
890	! if(specOLR)then
891	! open(117,file='OLRnu.out')
892	! do nw=1,L_NSPECTI
893	! write(117,*) OLR_nu(1,nw)
894	! enddo
895	! close(117)
896	!
897	! open(127,file='OSRnu.out')
898	! do nw=1,L_NSPECTV
899	! write(127,*) OSR_nu(1,nw)
900	! enddo
901	! close(127)
902	! endif
903
904	! OLR vs altitude: do it as a .txt file
905	OLRz=.false.
906	if(OLRz)then
907	print*,'saving IR vertical flux for OLRz...'
908	open(118,file='OLRz_plevs.out')
909	open(119,file='OLRz.out')
910	do l=1,L_NLAYRAD
911	write(118,) plevrad(2l)
912	do nw=1,L_NSPECTI
913	write(119,*) fluxupi_nu(l,nw)
914	enddo
915	enddo
916	close(118)
917	close(119)
918	endif
919
920	endif
921
922	! see physiq.F for explanations about CLFvarying. This is temporary.
923	if (lastcall .and. .not.CLFvarying) then
924	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
925	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
926	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
927	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
928	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
929	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
930	endif
931
932
933	end subroutine callcorrk

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

Download in other formats:

Original Format