Context Navigation

callcorrk.F90 @ 919

Last change on this file since 919 was 918, checked in by jleconte, 12 years ago

28/03/2013 == JL

optimization of optci and optcv routines. 15to 25% gain on these routines.

around 10% on the whole code with 1 scatterer.

No changes on output (byte to byte)

corrected bug in gray case in callcorrk.
added profiling option in makegcm_ifort. See the file for details

Property svn:executable set to *

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

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

Download in other formats:

Original Format