Context Navigation

callcorrk.F90 @ 537

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

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

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

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

Download in other formats:

Original Format