Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callcorrk.F90 @ 224

Last change on this file since 224 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

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

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

Download in other formats: