Context Navigation

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

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

Small fix for openMP

YM

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

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

Download in other formats: