Context Navigation

callcorrk.F90 @ 1704

Last change on this file since 1704 was 1699, checked in by jleconte, 8 years ago

21/04/2017 == JL

Add some arch files for a cluster in Bordeaux
added some is_master before printouts in callcorrk and physiq_mod
corrected a bug in bilinear big

Property svn:executable set to *

File size: 38.8 KB

Line
1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
3	tsurf,fract,dist_star,aerosol,muvar, &
4	dtlw,dtsw,fluxsurf_lw, &
5	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
6	fluxabs_sw,fluxtop_dn, &
7	OLR_nu,OSR_nu, &
8	tau_col,cloudfrac,totcloudfrac, &
9	clearsky,firstcall,lastcall)
10
11	use mod_phys_lmdz_para, only : is_master
12	use radinc_h
13	use radcommon_h
14	use watercommon_h
15	use datafile_mod, only: datadir
16	use ioipsl_getin_p_mod, only: getin_p
17	use gases_h
18	use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
19	use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_nh3, iaero_aurora
20	USE tracer_h
21	use comcstfi_mod, only: pi, mugaz, cpp
22	use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval, &
23	kastprof,strictboundcorrk,specOLR,CLFvarying
24
25	implicit none
26
27	!==================================================================
28	!
29	! Purpose
30	! -------
31	! Solve the radiative transfer using the correlated-k method for
32	! the gaseous absorption and the Toon et al. (1989) method for
33	! scatttering due to aerosols.
34	!
35	! Authors
36	! -------
37	! Emmanuel 01/2001, Forget 09/2001
38	! Robin Wordsworth (2009)
39	!
40	!==================================================================
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
48
49	! INPUT
50	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
51	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
52	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
53	INTEGER,INTENT(IN) :: nq ! Number of tracers.
54	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
55	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
56	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
57	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
58	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
59	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
60	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
61	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
62	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
63	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
64	REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
65	REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer) ! Fraction of clouds (%).
66	logical,intent(in) :: clearsky
67	logical,intent(in) :: firstcall ! Signals first call to physics.
68	logical,intent(in) :: lastcall ! Signals last call to physics.
69
70	! OUTPUT
71	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
72	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
73	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
74	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
75	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
76	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
77	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
78	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
79	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
80	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
81	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
82	REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity.
83	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
84	REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! Column Fraction of clouds (%).
85
86
87
88
89
90	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
91	REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
92	REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
93	REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
94	REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
95	REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
96	REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
97
98	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
99	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
100
101	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
102	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
103	!$OMP THREADPRIVATE(reffrad,nueffrad)
104
105	!-----------------------------------------------------------------------
106	! Declaration of the variables required by correlated-k subroutines
107	! Numbered from top to bottom (unlike in the GCM)
108	!-----------------------------------------------------------------------
109
110	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
111	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
112
113	! Optical values for the optci/cv subroutines
114	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
115	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
116	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
117	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
118	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
119	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
120	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
121	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
122	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
123	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
124
125	REAL*8 tauaero(L_LEVELS+1,naerkind)
126	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
127	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
128	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
129	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
130	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
131	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
132	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
133	REAL*8 albi,acosz
134	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
135
136	INTEGER ig,l,k,nw,iaer
137
138	real szangle
139	logical global1d
140	save szangle,global1d
141	!$OMP THREADPRIVATE(szangle,global1d)
142	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
143	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
144	real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol).
145
146	! Local aerosol optical properties for each column on RADIATIVE grid.
147	real*8,save,allocatable :: QXVAER(:,:,:)
148	real*8,save,allocatable :: QSVAER(:,:,:)
149	real*8,save,allocatable :: GVAER(:,:,:)
150	real*8,save,allocatable :: QXIAER(:,:,:)
151	real*8,save,allocatable :: QSIAER(:,:,:)
152	real*8,save,allocatable :: GIAER(:,:,:)
153
154	real, dimension(:,:,:), save, allocatable :: QREFvis3d
155	real, dimension(:,:,:), save, allocatable :: QREFir3d
156	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
157
158
159	! Miscellaneous :
160	real*8 temp,temp1,temp2,pweight
161	character(len=10) :: tmp1
162	character(len=10) :: tmp2
163
164	! For fixed water vapour profiles.
165	integer i_var
166	real RH
167	real*8 pq_temp(nlayer)
168	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
169	real ptemp, Ttemp, qsat
170
171	logical OLRz
172	real*8 NFLUXGNDV_nu(L_NSPECTV)
173
174	! Included by RW for runaway greenhouse 1D study.
175	real vtmp(nlayer)
176	REAL*8 muvarrad(L_LEVELS)
177
178	! Included by MT for albedo calculations.
179	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
180	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
181
182
183	!===============================================================
184	! I.a Initialization on first call
185	!===============================================================
186
187
188	if(firstcall) then
189
190	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
191	if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
192	if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
193	if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
194	if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
195	if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
196	if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
197
198	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
199	IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
200	IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,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,nlayer,reffrad,nueffrad)
212
213
214	!--------------------------------------------------
215	! Set up correlated k
216	!--------------------------------------------------
217
218
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
247	! Test of incompatibility : 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	! I.b Initialization on every call
266	!=======================================================================
267
268	qxvaer(:,:,:)=0.0
269	qsvaer(:,:,:)=0.0
270	gvaer(:,:,:) =0.0
271
272	qxiaer(:,:,:)=0.0
273	qsiaer(:,:,:)=0.0
274	giaer(:,:,:) =0.0
275
276	!--------------------------------------------------
277	! Effective radius and variance of the aerosols
278	!--------------------------------------------------
279
280	do iaer=1,naerkind
281
282	if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
283	call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
284	if (is_master) then
285	print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
286	print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
287	end if
288	end if
289
290	if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
291	call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
292	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
293	if (is_master) then
294	print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
295	print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
296	end if
297	endif
298
299	if(iaer.eq.iaero_dust)then
300	call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
301	if (is_master) then
302	print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
303	end if
304	endif
305
306	if(iaer.eq.iaero_h2so4)then
307	call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
308	if (is_master) then
309	print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
310	end if
311	endif
312
313	if(iaer.eq.iaero_back2lay)then
314	call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
315	endif
316	! if(iaer.eq.iaero_nh3)then
317	! call nh3_reffrad(ngrid,nlayer,reffrad(1,1,iaero_nh3))
318	! endif
319	! if(iaer.eq.iaero_aurora)then
320	! call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora))
321	! endif
322
323	end do !iaer=1,naerkind.
324
325
326	! How much light do we get ?
327	do nw=1,L_NSPECTV
328	stel(nw)=stellarf(nw)/(dist_star**2)
329	end do
330
331	! Get 3D aerosol optical properties.
332	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
333	QVISsQREF3d,omegaVIS3d,gVIS3d, &
334	QIRsQREF3d,omegaIR3d,gIR3d, &
335	QREFvis3d,QREFir3d)
336
337	! Get aerosol optical depths.
338	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
339	reffrad,QREFvis3d,QREFir3d, &
340	tau_col,cloudfrac,totcloudfrac,clearsky)
341
342
343
344	!-----------------------------------------------------------------------
345	do ig=1,ngrid ! Starting Big Loop over every GCM column
346	!-----------------------------------------------------------------------
347
348
349	!=======================================================================
350	! II. Transformation of the GCM variables
351	!=======================================================================
352
353
354	!-----------------------------------------------------------------------
355	! Aerosol optical properties Qext, Qscat and g.
356	! The transformation in the vertical is the same as for temperature.
357	!-----------------------------------------------------------------------
358
359
360	do iaer=1,naerkind
361	! Shortwave.
362	do nw=1,L_NSPECTV
363
364	do l=1,nlayer
365
366	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
367	*QREFvis3d(ig,nlayer+1-l,iaer)
368
369	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
370	*QREFvis3d(ig,max(nlayer-l,1),iaer)
371
372	qxvaer(2*l,nw,iaer) = temp1
373	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
374
375	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
376	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
377
378	qsvaer(2*l,nw,iaer) = temp1
379	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
380
381	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
382	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
383
384	gvaer(2*l,nw,iaer) = temp1
385	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
386
387	end do ! nlayer
388
389	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
390	qxvaer(2*nlayer+1,nw,iaer)=0.
391
392	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
393	qsvaer(2*nlayer+1,nw,iaer)=0.
394
395	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
396	gvaer(2*nlayer+1,nw,iaer)=0.
397
398	end do ! L_NSPECTV
399
400	do nw=1,L_NSPECTI
401	! Longwave
402	do l=1,nlayer
403
404	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
405	*QREFir3d(ig,nlayer+1-l,iaer)
406
407	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
408	*QREFir3d(ig,max(nlayer-l,1),iaer)
409
410	qxiaer(2*l,nw,iaer) = temp1
411	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
412
413	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
414	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
415
416	qsiaer(2*l,nw,iaer) = temp1
417	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
418
419	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
420	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
421
422	giaer(2*l,nw,iaer) = temp1
423	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
424
425	end do ! nlayer
426
427	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
428	qxiaer(2*nlayer+1,nw,iaer)=0.
429
430	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
431	qsiaer(2*nlayer+1,nw,iaer)=0.
432
433	giaer(1,nw,iaer)=giaer(2,nw,iaer)
434	giaer(2*nlayer+1,nw,iaer)=0.
435
436	end do ! L_NSPECTI
437
438	end do ! naerkind
439
440	! Test / Correct for freaky s. s. albedo values.
441	do iaer=1,naerkind
442	do k=1,L_LEVELS+1
443
444	do nw=1,L_NSPECTV
445	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
446	print*,'Serious problems with qsvaer values'
447	print*,'in callcorrk'
448	call abort
449	endif
450	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
451	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
452	endif
453	end do
454
455	do nw=1,L_NSPECTI
456	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
457	print*,'Serious problems with qsiaer values'
458	print*,'in callcorrk'
459	call abort
460	endif
461	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
462	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
463	endif
464	end do
465
466	end do ! L_LEVELS
467	end do ! naerkind
468
469	!-----------------------------------------------------------------------
470	! Aerosol optical depths
471	!-----------------------------------------------------------------------
472
473	do iaer=1,naerkind ! a bug was here
474	do k=0,nlayer-1
475
476	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
477	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
478	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
479	tauaero(2k+2,iaer)=max(temppweight,0.d0)
480	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
481
482	end do
483	! boundary conditions
484	tauaero(1,iaer) = tauaero(2,iaer)
485	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
486	!tauaero(1,iaer) = 0.
487	!tauaero(L_LEVELS+1,iaer) = 0.
488
489	end do ! naerkind
490
491	! Albedo and Emissivity.
492	albi=1-emis(ig) ! Long Wave.
493	DO nw=1,L_NSPECTV ! Short Wave loop.
494	albv(nw)=albedo(ig,nw)
495	ENDDO
496
497	if (nosurf) then ! Case with no surface.
498	DO nw=1,L_NSPECTV
499	if(albv(nw).gt.0.0) then
500	print*,'For open lower boundary in callcorrk must'
501	print*,'have spectral surface band albedos all set to zero!'
502	call abort
503	endif
504	ENDDO
505	endif
506
507	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
508	acosz = cos(pi*szangle/180.0)
509	print*,'acosz=',acosz,', szangle=',szangle
510	else
511	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
512	endif
513
514	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
515	!!! Note by JL13 : In the following, some indices were changed in the interpolations,
516	!!! so that the model results are less dependent on the number of layers !
517	!!!
518	!!! --- The older versions are commented with the comment !JL13index ---
519	!!!
520	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
521
522
523	!-----------------------------------------------------------------------
524	! Water vapour (to be generalised for other gases eventually ...)
525	!-----------------------------------------------------------------------
526
527	if(varactive)then
528
529	i_var=igcm_h2o_vap
530	do l=1,nlayer
531	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
532	qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)
533	!JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
534	!JL13index ! Average approximation as for temperature...
535	end do
536	qvar(1)=qvar(2)
537
538	elseif(varfixed)then
539
540	do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
541	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
542	if(RH.lt.0.0) RH=0.0
543
544	ptemp=pplay(ig,l)
545	Ttemp=pt(ig,l)
546	call watersat(Ttemp,ptemp,qsat)
547
548	!pq_temp(l) = qsat ! fully saturated everywhere
549	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
550	end do
551
552	do l=1,nlayer
553	qvar(2*l) = pq_temp(nlayer+1-l)
554	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
555	end do
556
557	qvar(1)=qvar(2)
558
559	! Lowest layer of atmosphere
560	RH = satval * (1 - 0.02) / 0.98
561	if(RH.lt.0.0) RH=0.0
562
563	! ptemp = pplev(ig,1)
564	! Ttemp = tsurf(ig)
565	! call watersat(Ttemp,ptemp,qsat)
566
567	qvar(2nlayer+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
568
569	else
570	do k=1,L_LEVELS
571	qvar(k) = 1.0D-7
572	end do
573	end if ! varactive/varfixed
574
575	if(.not.kastprof)then
576	! IMPORTANT: Now convert from kg/kg to mol/mol.
577	do k=1,L_LEVELS
578	qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
579	end do
580	end if
581
582	!-----------------------------------------------------------------------
583	! kcm mode only !
584	!-----------------------------------------------------------------------
585
586	if(kastprof)then
587
588	! Initial values equivalent to mugaz.
589	DO l=1,nlayer
590	muvarrad(2*l) = mugaz
591	muvarrad(2*l+1) = mugaz
592	END DO
593
594	if(ngasmx.gt.1)then
595
596	DO l=1,nlayer
597	muvarrad(2*l) = muvar(ig,nlayer+2-l)
598	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
599	muvar(ig,max(nlayer+1-l,1)))/2
600	END DO
601
602	muvarrad(1) = muvarrad(2)
603	muvarrad(2*nlayer+1) = muvar(ig,1)
604
605	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
606	print*,'Assumes that the variable gas is H2O!!!'
607	print*,'Assumes that there is only one tracer'
608
609	!i_var=igcm_h2o_vap
610	i_var=1
611
612	if(nq.gt.1)then
613	print*,'Need 1 tracer only to run kcm1d.e'
614	stop
615	endif
616
617	do l=1,nlayer
618	vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi))
619	!vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
620	end do
621
622	do l=1,nlayer
623	qvar(2*l) = vtmp(nlayer+1-l)
624	qvar(2*l+1) = vtmp(nlayer+1-l)
625	! qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
626	end do
627	qvar(1)=qvar(2)
628
629	print*,'Warning: reducing qvar in callcorrk.'
630	print*,'Temperature profile no longer consistent ', &
631	'with saturated H2O. qsat=',satval
632
633	do k=1,L_LEVELS
634	qvar(k) = qvar(k)*satval
635	end do
636
637	endif
638	else ! if kastprof
639	DO l=1,nlayer
640	muvarrad(2*l) = muvar(ig,nlayer+2-l)
641	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
642	END DO
643
644	muvarrad(1) = muvarrad(2)
645	muvarrad(2*nlayer+1)=muvar(ig,1)
646	endif ! if kastprof
647
648	! Keep values inside limits for which we have radiative transfer coefficients !!!
649	if(L_REFVAR.gt.1)then ! (there was a bug here)
650	do k=1,L_LEVELS
651	if(qvar(k).lt.wrefvar(1))then
652	qvar(k)=wrefvar(1)+1.0e-8
653	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
654	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
655	endif
656	end do
657	endif
658
659	!-----------------------------------------------------------------------
660	! Pressure and temperature
661	!-----------------------------------------------------------------------
662
663	DO l=1,nlayer
664	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
665	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
666	tlevrad(2*l) = pt(ig,nlayer+1-l)
667	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
668	END DO
669
670	plevrad(1) = 0.
671	plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
672
673	tlevrad(1) = tlevrad(2)
674	tlevrad(2*nlayer+1)=tsurf(ig)
675
676	pmid(1) = max(pgasmin,0.0001*plevrad(3))
677	pmid(2) = pmid(1)
678
679	tmid(1) = tlevrad(2)
680	tmid(2) = tmid(1)
681
682	DO l=1,L_NLAYRAD-1
683	tmid(2l+1) = tlevrad(2l+1)
684	tmid(2l+2) = tlevrad(2l+1)
685	pmid(2l+1) = plevrad(2l+1)
686	pmid(2l+2) = plevrad(2l+1)
687	END DO
688	pmid(L_LEVELS) = plevrad(L_LEVELS)
689	tmid(L_LEVELS) = tlevrad(L_LEVELS)
690
691	!!Alternative interpolation:
692	! pmid(3) = pmid(1)
693	! pmid(4) = pmid(1)
694	! tmid(3) = tmid(1)
695	! tmid(4) = tmid(1)
696	! DO l=2,L_NLAYRAD-1
697	! tmid(2l+1) = tlevrad(2l)
698	! tmid(2l+2) = tlevrad(2l)
699	! pmid(2l+1) = plevrad(2l)
700	! pmid(2l+2) = plevrad(2l)
701	! END DO
702	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
703	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
704
705	! Test for out-of-bounds pressure.
706	if(plevrad(3).lt.pgasmin)then
707	print*,'Minimum pressure is outside the radiative'
708	print*,'transfer kmatrix bounds, exiting.'
709	call abort
710	elseif(plevrad(L_LEVELS).gt.pgasmax)then
711	print*,'Maximum pressure is outside the radiative'
712	print*,'transfer kmatrix bounds, exiting.'
713	call abort
714	endif
715
716	! Test for out-of-bounds temperature.
717	do k=1,L_LEVELS
718	if(tlevrad(k).lt.tgasmin)then
719	print*,'Minimum temperature is outside the radiative'
720	print*,'transfer kmatrix bounds'
721	print*,"k=",k," tlevrad(k)=",tlevrad(k)
722	print*,"tgasmin=",tgasmin
723	if (strictboundcorrk) then
724	call abort
725	else
726	print,'**********************************************'
727	print*,'we allow model to continue with tlevrad=tgasmin'
728	print*,' ... we assume we know what you are doing ... '
729	print*,' ... but do not let this happen too often ... '
730	print,'**********************************************'
731	!tlevrad(k)=tgasmin
732	endif
733	elseif(tlevrad(k).gt.tgasmax)then
734	print*,'Maximum temperature is outside the radiative'
735	print*,'transfer kmatrix bounds, exiting.'
736	print*,"k=",k," tlevrad(k)=",tlevrad(k)
737	print*,"tgasmax=",tgasmax
738	if (strictboundcorrk) then
739	call abort
740	else
741	print,'**********************************************'
742	print*,'we allow model to continue with tlevrad=tgasmax'
743	print*,' ... we assume we know what you are doing ... '
744	print*,' ... but do not let this happen too often ... '
745	print,'**********************************************'
746	!tlevrad(k)=tgasmax
747	endif
748	endif
749	enddo
750	do k=1,L_NLAYRAD+1
751	if(tmid(k).lt.tgasmin)then
752	print*,'Minimum temperature is outside the radiative'
753	print*,'transfer kmatrix bounds, exiting.'
754	print*,"k=",k," tmid(k)=",tmid(k)
755	print*,"tgasmin=",tgasmin
756	if (strictboundcorrk) then
757	call abort
758	else
759	print,'**********************************************'
760	print*,'we allow model to continue with tmid=tgasmin'
761	print*,' ... we assume we know what you are doing ... '
762	print*,' ... but do not let this happen too often ... '
763	print,'**********************************************'
764	tmid(k)=tgasmin
765	endif
766	elseif(tmid(k).gt.tgasmax)then
767	print*,'Maximum temperature is outside the radiative'
768	print*,'transfer kmatrix bounds, exiting.'
769	print*,"k=",k," tmid(k)=",tmid(k)
770	print*,"tgasmax=",tgasmax
771	if (strictboundcorrk) then
772	call abort
773	else
774	print,'**********************************************'
775	print*,'we allow model to continue with tmid=tgasmin'
776	print*,' ... we assume we know what you are doing ... '
777	print*,' ... but do not let this happen too often ... '
778	print,'**********************************************'
779	tmid(k)=tgasmax
780	endif
781	endif
782	enddo
783
784	!=======================================================================
785	! III. Calling the main radiative transfer subroutines
786	!=======================================================================
787
788
789	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
790	glat_ig=glat(ig)
791
792	!-----------------------------------------------------------------------
793	! Short Wave Part
794	!-----------------------------------------------------------------------
795
796	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
797	if((ngrid.eq.1).and.(global1d))then
798	do nw=1,L_NSPECTV
799	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
800	end do
801	else
802	do nw=1,L_NSPECTV
803	stel_fract(nw)= stel(nw) * fract(ig)
804	end do
805	endif
806
807	call optcv(dtauv,tauv,taucumv,plevrad, &
808	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
809	tmid,pmid,taugsurf,qvar,muvarrad)
810
811	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
812	acosz,stel_fract,gweight, &
813	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
814	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
815
816	else ! During the night, fluxes = 0.
817	nfluxtopv = 0.0d0
818	fluxtopvdn = 0.0d0
819	nfluxoutv_nu(:) = 0.0d0
820	nfluxgndv_nu(:) = 0.0d0
821	do l=1,L_NLAYRAD
822	fmnetv(l)=0.0d0
823	fluxupv(l)=0.0d0
824	fluxdnv(l)=0.0d0
825	end do
826	end if
827
828
829	! Equivalent Albedo Calculation (for OUTPUT). MT2015
830	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
831	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
832	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
833	DO nw=1,L_NSPECTV
834	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
835	ENDDO
836	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
837	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
838	else
839	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
840	endif
841	else
842	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
843	endif
844
845
846	!-----------------------------------------------------------------------
847	! Long Wave Part
848	!-----------------------------------------------------------------------
849
850	call optci(plevrad,tlevrad,dtaui,taucumi, &
851	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
852	taugsurfi,qvar,muvarrad)
853
854	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
855	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
856	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
857
858	!-----------------------------------------------------------------------
859	! Transformation of the correlated-k code outputs
860	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
861
862	! Flux incident at the top of the atmosphere
863	fluxtop_dn(ig)=fluxtopvdn
864
865	fluxtop_lw(ig) = real(nfluxtopi)
866	fluxabs_sw(ig) = real(-nfluxtopv)
867	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
868	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
869
870	! Flux absorbed by the surface. By MT2015.
871	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
872
873	if(fluxtop_dn(ig).lt.0.0)then
874	print*,'Achtung! fluxtop_dn has lost the plot!'
875	print*,'fluxtop_dn=',fluxtop_dn(ig)
876	print*,'acosz=',acosz
877	print*,'aerosol=',aerosol(ig,:,:)
878	print*,'temp= ',pt(ig,:)
879	print*,'pplay= ',pplay(ig,:)
880	call abort
881	endif
882
883	! Spectral output, for exoplanet observational comparison
884	if(specOLR)then
885	do nw=1,L_NSPECTI
886	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
887	end do
888	do nw=1,L_NSPECTV
889	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
890	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
891	end do
892	endif
893
894	! Finally, the heating rates
895
896	DO l=2,L_NLAYRAD
897	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
898	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
899	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
900	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
901	END DO
902
903	! These are values at top of atmosphere
904	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
905	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
906	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
907	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
908
909
910	!-----------------------------------------------------------------------
911	end do ! End of big loop over every GCM column.
912	!-----------------------------------------------------------------------
913
914
915
916	!-----------------------------------------------------------------------
917	! Additional diagnostics
918	!-----------------------------------------------------------------------
919
920	! IR spectral output, for exoplanet observational comparison
921	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
922
923	print*,'Saving scalar quantities in surf_vals.out...'
924	print*,'psurf = ', pplev(1,1),' Pa'
925	open(116,file='surf_vals.out')
926	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
927	real(-nfluxtopv),real(nfluxtopi)
928	close(116)
929
930
931	! USEFUL COMMENT - Do Not Remove.
932	!
933	! if(specOLR)then
934	! open(117,file='OLRnu.out')
935	! do nw=1,L_NSPECTI
936	! write(117,*) OLR_nu(1,nw)
937	! enddo
938	! close(117)
939	!
940	! open(127,file='OSRnu.out')
941	! do nw=1,L_NSPECTV
942	! write(127,*) OSR_nu(1,nw)
943	! enddo
944	! close(127)
945	! endif
946
947	! OLR vs altitude: do it as a .txt file.
948	OLRz=.false.
949	if(OLRz)then
950	print*,'saving IR vertical flux for OLRz...'
951	open(118,file='OLRz_plevs.out')
952	open(119,file='OLRz.out')
953	do l=1,L_NLAYRAD
954	write(118,) plevrad(2l)
955	do nw=1,L_NSPECTI
956	write(119,*) fluxupi_nu(l,nw)
957	enddo
958	enddo
959	close(118)
960	close(119)
961	endif
962
963	endif
964
965	! See physiq.F for explanations about CLFvarying. This is temporary.
966	if (lastcall .and. .not.CLFvarying) then
967	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
968	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
969	!$OMP BARRIER
970	!$OMP MASTER
971	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
972	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
973	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
974	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
975	!$OMP END MASTER
976	!$OMP BARRIER
977	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
978	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
979	endif
980
981
982	end subroutine callcorrk

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

Download in other formats:

Original Format