Context Navigation

callcorrk.F90 @ 1351

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

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

Download in other formats:

Original Format