Context Navigation

callcorrk.F90 @ 2340

Last change on this file since 2340 was 2297, checked in by jvatant, 6 years ago

Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.

It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)

*aeronlay_tauref (Optical depth of aerosol layer at ref wavelenght)
*aeronlay_lamref (Ref wavelenght (m))
*aeronlay_choice (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
*aeronlay_pbot (Bottom pressure (Pa))
*aeronlay_ptop (Top pressure (Pa) - useful only if choice=1)
*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height - useful only if choice=2 )
*aeronlay_size (Particle size (m))
*optprop_aeronlay_vis File for VIS opt properties.
*optprop_aeronlay_ir File for IR opt properties.

+Extra info :

+ In addition of solving the bug from 2-layer it enables different optical properties.
+ The former scheme are left for retrocompatibility (for now) but you should use the new one.
+ See aeropacity.F90 for the calculations

+ Each layer can have different optical properties, size of particle ...
+ You have different choices for vertical profile of the aerosol layers :

aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).

In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.

+ The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+ Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..

Property svn:executable set to *

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

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

Download in other formats:

Original Format