Context Navigation

callcorrk.F90 @ 747

Last change on this file since 747 was 728, checked in by jleconte, 12 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

Property svn:executable set to *

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

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

Download in other formats:

Original Format