Context Navigation

callcorrk.F90 @ 437

Last change on this file since 437 was 374, checked in by emillour, 14 years ago
Generic GCM: Upgrade: The location of the 'datagcm' directory can now be given in the callphys.def file ( datadir = /absolute/path/to/datagcm ). Changed "datafile.h" into a F90 module "datafile_mod.F90" and spread this change to all routines that used to use "datafile.h". EM
Property svn:executable set to ``*
File size: 29.7 KB

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

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

Download in other formats:

Original Format