Context Navigation

callcorrk.F90 @ 2276

Last change on this file since 2276 was 2269, checked in by emillour, 5 years ago

Generic GCM:

Cleanup OpenMP statements in callcorrk.F90 and also use "call abort_physic" instead of "stop" or "call abort"

Property svn:executable set to *

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

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

Download in other formats:

Original Format