Context Navigation

callcorrk.F90 @ 613

Last change on this file since 613 was 600, checked in by jleconte, 13 years ago

Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

OLR now equal to LW net heating/cooling at equilibrium!
As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

Property svn:executable set to *

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

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

Download in other formats:

Original Format