Context Navigation

callcorrk.F90 @ 1990

Last change on this file since 1990 was 1988, checked in by jleconte, 7 years ago

28/08/2018 == JL

We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv.
This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
Additional speedup corrections have been made in gfluxi that change nothing to the result.
Finally, if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv and changes in last commit.)
This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.

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

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

Download in other formats:

Original Format