Context Navigation

callcorrk.F90 @ 867

Last change on this file since 867 was 861, checked in by aslmd, 12 years ago

LMDZ.GENERIC

Fixed an allocating bug which arises from previous modifications owing to the double use of callcorrk with CLF_varying
Fixed a small bug in a diagnostic in the end of calc_rayleigh. Some picky compilers complain.
Fixed a small bug with the array noms which is not allocated when tracer is false. But still need in physdem1.

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 sugas_corrk ! set up gaseous absorption properties
226	call setspi ! basic infrared properties
227	call setspv ! basic visible 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