Context Navigation

callcorrk.F90 @ 1660

Last change on this file since 1660 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

Property svn:executable set to *

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

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

Download in other formats:

Original Format