Context Navigation

callcorrk.F90 @ 2099

Last change on this file since 2099 was 2051, checked in by jleconte, 7 years ago

12/12/2018 == JL

Correct a bug from commit 2032 in callcorrk: automatic arrays must be allocated

only if they haven't been before (callcorrk is called twice in physiq when cloud cover is not uniform).

Property svn:executable set to *

File size: 40.4 KB

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

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

Download in other formats:

Original Format