Context Navigation

callcorrk.F90 @ 1157

Last change on this file since 1157 was 1147, checked in by aslmd, 11 years ago
LMDZ.GENERIC. a suggestion by JL. do not modify tlevrad (if outside tgasmin or tgasmax) because used in source function.
Property svn:executable set to ``*
File size: 33.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	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,back2lay_reffrad
17	use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay
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,fluxtopvdn
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 szangle
128	logical global1d
129	save szangle,global1d
130	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
131	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
132
133	real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol)
134
135	! Local aerosol optical properties for each column on RADIATIVE grid
136	real*8,save :: QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
137	real*8,save :: QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
138	real*8,save :: GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
139	real*8,save :: QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
140	real*8,save :: QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
141	real*8,save :: GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
142
143	!REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
144	!REAL :: QREFir3d(ngrid,nlayermx,naerkind)
145	!save QREFvis3d, QREFir3d
146	real, dimension(:,:,:), save, allocatable :: QREFvis3d
147	real, dimension(:,:,:), save, allocatable :: QREFir3d
148
149
150	! Misc.
151	logical nantest
152	real*8 tempv(L_NSPECTV)
153	real*8 tempi(L_NSPECTI)
154	real*8 temp,temp1,temp2,pweight
155	character(len=10) :: tmp1
156	character(len=10) :: tmp2
157
158	! for fixed water vapour profiles
159	integer i_var
160	real RH
161	real*8 pq_temp(nlayer)
162	real ptemp, Ttemp, qsat
163
164	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
165
166	!real ptime, pday
167	logical OLRz
168	real*8 NFLUXGNDV_nu(L_NSPECTV)
169
170
171	! for weird cloud test
172	real pqtest(ngrid,nlayer,nq)
173
174	real maxrad, minrad
175
176	real,external :: CBRT
177
178	! included by RW for runaway greenhouse 1D study
179	real vtmp(nlayermx)
180	REAL*8 muvarrad(L_LEVELS)
181
182
183	!===============================================================
184	! Initialization on first call
185
186	qxvaer(:,:,:)=0.0
187	qsvaer(:,:,:)=0.0
188	gvaer(:,:,:) =0.0
189
190	qxiaer(:,:,:)=0.0
191	qsiaer(:,:,:)=0.0
192	giaer(:,:,:) =0.0
193
194	if(firstcall) then
195
196	!!! ALLOCATED instances are necessary because of CLFvarying
197	!!! strategy to call callcorrk twice in physiq...
198	IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayermx,naerkind))
199	IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayermx,naerkind))
200	! Effective radius and variance of the aerosols
201	IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
202	IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
203
204	call system('rm -f surf_vals_long.out')
205
206	if(naerkind.gt.4)then
207	print*,'Code not general enough to deal with naerkind > 4 yet.'
208	call abort
209	endif
210	call su_aer_radii(ngrid,reffrad,nueffrad)
211
212
213
214	!--------------------------------------------------
215	! set up correlated k
216	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
217	call getin("corrkdir",corrkdir)
218	print*, "corrkdir = ",corrkdir
219	write( tmp1, '(i3)' ) L_NSPECTI
220	write( tmp2, '(i3)' ) L_NSPECTV
221	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
222	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
223
224	call setspi ! basic infrared properties
225	call setspv ! basic visible properties
226	call sugas_corrk ! set up gaseous absorption properties
227	call suaer_corrk ! set up aerosol optical properties
228
229	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
230
231	if((igcm_h2o_vap.eq.0) .and. varactive)then
232	print*,'varactive in callcorrk but no h2o_vap tracer.'
233	stop
234	endif
235
236	OLR_nu(:,:) = 0.
237	OSR_nu(:,:) = 0.
238
239	if (ngrid.eq.1) then
240	PRINT*, 'Simulate global averaged conditions ?'
241	global1d = .false. ! default value
242	call getin("global1d",global1d)
243	write(,) "global1d = ",global1d
244	! Test of incompatibility:
245	! if global1d is true, there should not be any diurnal cycle
246	if (global1d.and.diurnal) then
247	print*,'if global1d is true, diurnal must be set to false'
248	stop
249	endif
250
251	if (global1d) then
252	PRINT *,'Solar Zenith angle (deg.) ?'
253	PRINT *,'(assumed for averaged solar flux S/4)'
254	szangle=60.0 ! default value
255	call getin("szangle",szangle)
256	write(,) "szangle = ",szangle
257	endif
258	endif
259
260	end if ! of if (firstcall)
261
262	!=======================================================================
263	! Initialization on every call
264
265	!--------------------------------------------------
266	! Effective radius and variance of the aerosols
267	do iaer=1,naerkind
268
269	if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
270	call co2_reffrad(ngrid,nq,pq,reffrad(1,1,iaero_co2))
271	print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
272	print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
273	end if
274	if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
275	call h2o_reffrad(ngrid,pq(1,1,igcm_h2o_ice),pt, &
276	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
277	print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
278	print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
279	endif
280	if(iaer.eq.iaero_dust)then
281	call dust_reffrad(ngrid,reffrad(1,1,iaero_dust))
282	print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
283	endif
284	if(iaer.eq.iaero_h2so4)then
285	call h2so4_reffrad(ngrid,reffrad(1,1,iaero_h2so4))
286	print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
287	endif
288	if(iaer.eq.iaero_back2lay)then
289	call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
290	endif
291	end do !iaer=1,naerkind
292
293
294	! how much light we get
295	do nw=1,L_NSPECTV
296	stel(nw)=stellarf(nw)/(dist_star**2)
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(nosurf.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	!!! JL13: in the following, I changed some indices in the interpolations so that the model rsults are less dependent on the number of layers
466	!!! the older verions are commented with the commetn !JL13index
467
468
469	!-----------------------------------------------------------------------
470	! Water vapour (to be generalised for other gases eventually)
471
472	if(varactive)then
473
474	i_var=igcm_h2o_vap
475	do l=1,nlayer
476	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
477	qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)
478	!JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
479	!JL13index ! Average approximation as for temperature...
480	end do
481	qvar(1)=qvar(2)
482
483	elseif(varfixed)then
484
485	do l=1,nlayermx ! here we will assign fixed water vapour profiles globally
486	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
487	if(RH.lt.0.0) RH=0.0
488
489	ptemp=pplay(ig,l)
490	Ttemp=pt(ig,l)
491	call watersat(Ttemp,ptemp,qsat)
492
493	!pq_temp(l) = qsat ! fully saturated everywhere
494	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
495	end do
496
497	do l=1,nlayer
498	qvar(2*l) = pq_temp(nlayer+1-l)
499	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
500	end do
501	qvar(1)=qvar(2)
502
503	! Lowest layer of atmosphere
504	RH = satval * (1 - 0.02) / 0.98
505	if(RH.lt.0.0) RH=0.0
506
507	! ptemp = pplev(ig,1)
508	! Ttemp = tsurf(ig)
509	! call watersat(Ttemp,ptemp,qsat)
510
511	qvar(2nlayermx+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
512
513	else
514	do k=1,L_LEVELS
515	qvar(k) = 1.0D-7
516	end do
517	end if
518
519	if(.not.kastprof)then
520	! IMPORTANT: Now convert from kg/kg to mol/mol
521	do k=1,L_LEVELS
522	qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
523	end do
524	end if
525
526	!-----------------------------------------------------------------------
527	! kcm mode only
528	if(kastprof)then
529
530	! initial values equivalent to mugaz
531	DO l=1,nlayer
532	muvarrad(2*l) = mugaz
533	muvarrad(2*l+1) = mugaz
534	END DO
535
536	if(ngasmx.gt.1)then
537
538	DO l=1,nlayer
539	muvarrad(2*l) = muvar(ig,nlayer+2-l)
540	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
541	muvar(ig,max(nlayer+1-l,1)))/2
542	END DO
543
544	muvarrad(1) = muvarrad(2)
545	muvarrad(2*nlayermx+1)=muvar(ig,1)
546
547	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
548	print*,'Assumes that the variable gas is H2O!!!'
549	print*,'Assumes that there is only one tracer'
550	!i_var=igcm_h2o_vap
551	i_var=1
552	if(nq.gt.1)then
553	print*,'Need 1 tracer only to run kcm1d.e'
554	stop
555	endif
556	do l=1,nlayer
557	vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi))
558	!vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
559	end do
560
561	do l=1,nlayer
562	qvar(2*l) = vtmp(nlayer+1-l)
563	qvar(2*l+1) = vtmp(nlayer+1-l)
564	! qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
565	end do
566	qvar(1)=qvar(2)
567
568	print*,'Warning: reducing qvar in callcorrk.'
569	print*,'Temperature profile no longer consistent ', &
570	'with saturated H2O. qsat=',satval
571	do k=1,L_LEVELS
572	qvar(k) = qvar(k)*satval
573	end do
574
575	endif
576	else ! if kastprof
577	DO l=1,nlayer
578	muvarrad(2*l) = muvar(ig,nlayer+2-l)
579	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
580	END DO
581
582	muvarrad(1) = muvarrad(2)
583	muvarrad(2*nlayermx+1)=muvar(ig,1)
584	endif
585
586	! Keep values inside limits for which we have radiative transfer coefficients
587	if(L_REFVAR.gt.1)then ! there was a bug here!
588	do k=1,L_LEVELS
589	if(qvar(k).lt.wrefvar(1))then
590	qvar(k)=wrefvar(1)+1.0e-8
591	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
592	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
593	endif
594	end do
595	endif
596
597	!-----------------------------------------------------------------------
598	! Pressure and temperature
599
600	DO l=1,nlayer
601	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
602	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
603	tlevrad(2*l) = pt(ig,nlayer+1-l)
604	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
605	END DO
606
607	plevrad(1) = 0.
608	plevrad(2) = max(pgasmin,0.0001*plevrad(3))
609
610	tlevrad(1) = tlevrad(2)
611	tlevrad(2*nlayermx+1)=tsurf(ig)
612
613	tmid(1) = tlevrad(2)
614	tmid(2) = tlevrad(2)
615	pmid(1) = plevrad(2)
616	pmid(2) = plevrad(2)
617
618	DO l=1,L_NLAYRAD-1
619	tmid(2l+1) = tlevrad(2l)
620	tmid(2l+2) = tlevrad(2l)
621	pmid(2l+1) = plevrad(2l)
622	pmid(2l+2) = plevrad(2l)
623	!JL13index tmid(2l+1) = tlevrad(2l+1)
624	!JL13index tmid(2l+2) = tlevrad(2l+1)
625	!JL13index pmid(2l+1) = plevrad(2l+1)
626	!JL13index pmid(2l+2) = plevrad(2l+1)
627	END DO
628	pmid(L_LEVELS) = plevrad(L_LEVELS-1)
629	tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
630	!JL13index pmid(L_LEVELS) = plevrad(L_LEVELS)
631	!JL13index tmid(L_LEVELS) = tlevrad(L_LEVELS)
632
633	! test for out-of-bounds pressure
634	if(plevrad(3).lt.pgasmin)then
635	print*,'Minimum pressure is outside the radiative'
636	print*,'transfer kmatrix bounds, exiting.'
637	call abort
638	elseif(plevrad(L_LEVELS).gt.pgasmax)then
639	print*,'Maximum pressure is outside the radiative'
640	print*,'transfer kmatrix bounds, exiting.'
641	call abort
642	endif
643
644	! test for out-of-bounds temperature
645	do k=1,L_LEVELS
646	if(tlevrad(k).lt.tgasmin)then
647	print*,'Minimum temperature is outside the radiative'
648	print*,'transfer kmatrix bounds'
649	print*,"k=",k," tlevrad(k)=",tlevrad(k)
650	print*,"tgasmin=",tgasmin
651	if (strictboundcorrk) then
652	call abort
653	else
654	print,'**********************************************'
655	print*,'we allow model to continue with tlevrad=tgasmin'
656	print*,' ... we assume we know what you are doing ... '
657	print*,' ... but do not let this happen too often ... '
658	print,'**********************************************'
659	!tlevrad(k)=tgasmin
660	endif
661	elseif(tlevrad(k).gt.tgasmax)then
662	print*,'Maximum temperature is outside the radiative'
663	print*,'transfer kmatrix bounds, exiting.'
664	print*,"k=",k," tlevrad(k)=",tlevrad(k)
665	print*,"tgasmax=",tgasmax
666	if (strictboundcorrk) then
667	call abort
668	else
669	print,'**********************************************'
670	print*,'we allow model to continue with tlevrad=tgasmax'
671	print*,' ... we assume we know what you are doing ... '
672	print*,' ... but do not let this happen too often ... '
673	print,'**********************************************'
674	!tlevrad(k)=tgasmax
675	endif
676	endif
677	enddo
678	do k=1,L_NLAYRAD+1
679	if(tmid(k).lt.tgasmin)then
680	print*,'Minimum temperature is outside the radiative'
681	print*,'transfer kmatrix bounds, exiting.'
682	print*,"k=",k," tmid(k)=",tmid(k)
683	print*,"tgasmin=",tgasmin
684	if (strictboundcorrk) then
685	call abort
686	else
687	print,'**********************************************'
688	print*,'we allow model to continue with tmid=tgasmin'
689	print*,' ... we assume we know what you are doing ... '
690	print*,' ... but do not let this happen too often ... '
691	print,'**********************************************'
692	tmid(k)=tgasmin
693	endif
694	elseif(tmid(k).gt.tgasmax)then
695	print*,'Maximum temperature is outside the radiative'
696	print*,'transfer kmatrix bounds, exiting.'
697	print*,"k=",k," tmid(k)=",tmid(k)
698	print*,"tgasmax=",tgasmax
699	if (strictboundcorrk) then
700	call abort
701	else
702	print,'**********************************************'
703	print*,'we allow model to continue with tmid=tgasmin'
704	print*,' ... we assume we know what you are doing ... '
705	print*,' ... but do not let this happen too often ... '
706	print,'**********************************************'
707	tmid(k)=tgasmax
708	endif
709	endif
710	enddo
711
712	!=======================================================================
713	! Calling the main radiative transfer subroutines
714
715
716	!-----------------------------------------------------------------------
717	! Shortwave
718
719	if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
720
721	if((ngrid.eq.1).and.(global1d))then
722	do nw=1,L_NSPECTV
723	stel_fract(nw)= stel(nw) * 0.25 / acosz
724	! globally averaged = divide by 4
725	! but we correct for solar zenith angle
726	end do
727	else
728	do nw=1,L_NSPECTV
729	stel_fract(nw)= stel(nw) * fract(ig)*(1.-eclipse(ig))
730	end do
731	endif
732
733	call optcv(dtauv,tauv,taucumv,plevrad, &
734	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
735	tmid,pmid,taugsurf,qvar,muvarrad)
736
737	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
738	acosz,stel_fract,gweight, &
739	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
740	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
741
742	else ! during the night, fluxes = 0
743	nfluxtopv = 0.0d0
744	fluxtopvdn = 0.0d0
745	nfluxoutv_nu(:) = 0.0d0
746	nfluxgndv_nu(:) = 0.0d0
747	do l=1,L_NLAYRAD
748	fmnetv(l)=0.0d0
749	fluxupv(l)=0.0d0
750	fluxdnv(l)=0.0d0
751	end do
752	end if
753
754	!-----------------------------------------------------------------------
755	! Longwave
756
757	call optci(plevrad,tlevrad,dtaui,taucumi, &
758	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
759	taugsurfi,qvar,muvarrad)
760
761	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
762	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
763	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
764
765	!-----------------------------------------------------------------------
766	! Transformation of the correlated-k code outputs
767	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
768
769	! Flux incident at the top of the atmosphere
770	fluxtop_dn(ig)=fluxtopvdn
771
772	fluxtop_lw(ig) = real(nfluxtopi)
773	fluxabs_sw(ig) = real(-nfluxtopv)
774	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
775	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
776
777	if(fluxtop_dn(ig).lt.0.0)then
778	print*,'Achtung! fluxtop_dn has lost the plot!'
779	print*,'fluxtop_dn=',fluxtop_dn(ig)
780	print*,'acosz=',acosz
781	print*,'aerosol=',aerosol(ig,:,:)
782	print*,'temp= ',pt(ig,:)
783	print*,'pplay= ',pplay(ig,:)
784	call abort
785	endif
786
787	! Spectral output, for exoplanet observational comparison
788	if(specOLR)then
789	do nw=1,L_NSPECTI
790	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
791	end do
792	do nw=1,L_NSPECTV
793	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
794	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
795	end do
796	endif
797
798	! Finally, the heating rates
799
800	DO l=2,L_NLAYRAD
801	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
802	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
803	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
804	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
805	END DO
806
807	! These are values at top of atmosphere
808	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
809	g/(cppscalep*(plevrad(3)-plevrad(1)))
810	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
811	g/(cppscalep*(plevrad(3)-plevrad(1)))
812
813	! ---------------------------------------------------------------
814	end do ! end of big loop over every GCM column (ig = 1:ngrid)
815
816
817	!-----------------------------------------------------------------------
818	! Additional diagnostics
819
820	! IR spectral output, for exoplanet observational comparison
821
822
823	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
824
825	print*,'Saving scalar quantities in surf_vals.out...'
826	print*,'psurf = ', pplev(1,1),' Pa'
827	open(116,file='surf_vals.out')
828	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
829	real(-nfluxtopv),real(nfluxtopi)
830	close(116)
831
832	! I am useful, please don`t remove me!
833	! if(specOLR)then
834	! open(117,file='OLRnu.out')
835	! do nw=1,L_NSPECTI
836	! write(117,*) OLR_nu(1,nw)
837	! enddo
838	! close(117)
839	!
840	! open(127,file='OSRnu.out')
841	! do nw=1,L_NSPECTV
842	! write(127,*) OSR_nu(1,nw)
843	! enddo
844	! close(127)
845	! endif
846
847	! OLR vs altitude: do it as a .txt file
848	OLRz=.false.
849	if(OLRz)then
850	print*,'saving IR vertical flux for OLRz...'
851	open(118,file='OLRz_plevs.out')
852	open(119,file='OLRz.out')
853	do l=1,L_NLAYRAD
854	write(118,) plevrad(2l)
855	do nw=1,L_NSPECTI
856	write(119,*) fluxupi_nu(l,nw)
857	enddo
858	enddo
859	close(118)
860	close(119)
861	endif
862
863	endif
864
865	! see physiq.F for explanations about CLFvarying. This is temporary.
866	if (lastcall .and. .not.CLFvarying) then
867	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
868	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
869	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
870	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
871	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
872	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
873	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
874	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
875	endif
876
877
878	end subroutine callcorrk

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

Download in other formats:

Original Format