Context Navigation

callcorrk.F90 @ 726

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

17/07/2012 == JL for LK

Generalization of aerosol scheme:
- any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
- addition of a module with the id numbers for aerosols (aerosol_mod.F90).
- initialization of aerosols id numbers in iniaerosol.F90
- compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

If starting from an old start file, recreate start file with the q=0 option in newstart.e.
update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

To add other aerosols, see Laura Kerber.

Property svn:executable set to *

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

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

Download in other formats:

Original Format