Context Navigation

callcorrk.F90 @ 3586

Last change on this file since 3586 was 3585, checked in by debatzbr, 12 days ago
Connecting microphysics to radiative transfer + miscellaneous cleans
Property svn:executable set to ``*
File size: 49.0 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	zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar, &
10	dtlw,dtsw,fluxsurf_lw, &
11	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
12	fluxabs_sw,fluxtop_dn, &
13	OLR_nu,OSR_nu,GSR_nu, &
14	int_dtaui,int_dtauv, &
15	tau_col,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 datafile_mod, only: datadir
25	use ioipsl_getin_p_mod, only: getin_p
26	use gases_h, only: ngasmx
27	use radii_mod, only : su_aer_radii, haze_reffrad_fix
28	use aerosol_mod, only : iaero_haze
29	use aeropacity_mod, only: aeropacity
30	use aeroptproperties_mod, only: aeroptproperties
31	use tracer_h, only: igcm_ch4_gas,igcm_n2,mmol
32	use comcstfi_mod, only: pi, mugaz, cpp, r, g
33	use callkeys_mod, only: diurnal,tracer,varfixed,satval, &
34	diagdtau,kastprof,strictboundcorrk,specOLR, &
35	tplanckmin,tplanckmax,global1d, &
36	optichaze,haze_radproffix,&
37	methane,carbox,cooling,nlte,strobel,&
38	ch4fix,vmrch4_proffix,vmrch4fix,&
39	callmufi
40	use optcv_mod, only: optcv
41	use optci_mod, only: optci
42	use sfluxi_mod, only: sfluxi
43	use sfluxv_mod, only: sfluxv
44	use recombin_corrk_mod, only: corrk_recombin, call_recombin
45	use planetwide_mod, only: planetwide_maxval, planetwide_minval
46	use radcommon_h, only: wavev,wavei
47	use mp2m_diagnostics
48	implicit none
49
50	!==================================================================
51	!
52	! Purpose
53	! -------
54	! Solve the radiative transfer using the correlated-k method for
55	! the gaseous absorption and the Toon et al. (1989) method for
56	! scatttering due to aerosols.
57	!
58	! Authors
59	! -------
60	! Emmanuel 01/2001, Forget 09/2001
61	! Robin Wordsworth (2009)
62	!
63	!==================================================================
64
65	!-----------------------------------------------------------------------
66	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
67	! Layer #1 is the layer near the ground.
68	! Layer #nlayer is the layer at the top.
69	!-----------------------------------------------------------------------
70
71
72	! INPUT
73	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
74	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
75	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
76	INTEGER,INTENT(IN) :: nq ! Number of tracers.
77	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
78	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
79	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
80	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
81	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
82	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
83	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
84	REAL,INTENT(IN) :: zzlay(ngrid,nlayer) ! Mid-layer altitude
85	REAL,INTENT(IN) :: zzlev(ngrid,nlayer) ! Altitude at the layer boundaries.
86	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
87	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
88	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
89	REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
90	logical,intent(in) :: firstcall ! Signals first call to physics.
91	logical,intent(in) :: lastcall ! Signals last call to physics.
92
93	! OUTPUT
94	REAL,INTENT(OUT) :: dtau_aer(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
95	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
96	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
97	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
98	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
99	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
100	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
101	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
102	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
103	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
104	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
105	REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV) ! Surface SW radiation in each band (Normalized to the band width (W/m2/cm-1).
106	REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity.
107	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
108	REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
109	REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
110
111
112
113
114
115	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
116	! made "save" variables so they are allocated once in for all, not because
117	! the values need be saved from a time step to the next
118	REAL,SAVE,ALLOCATABLE :: QVISsQREF3d(:,:,:,:)
119	REAL,SAVE,ALLOCATABLE :: omegaVIS3d(:,:,:,:)
120	REAL,SAVE,ALLOCATABLE :: gVIS3d(:,:,:,:)
121	!$OMP THREADPRIVATE(QVISsQREF3d,omegaVIS3d,gVIS3d)
122	REAL,SAVE,ALLOCATABLE :: QIRsQREF3d(:,:,:,:)
123	REAL,SAVE,ALLOCATABLE :: omegaIR3d(:,:,:,:)
124	REAL,SAVE,ALLOCATABLE :: gIR3d(:,:,:,:)
125	!$OMP THREADPRIVATE(QIRsQREF3d,omegaIR3d,gIR3d)
126
127	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
128	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
129
130	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
131	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
132	!$OMP THREADPRIVATE(reffrad,nueffrad)
133
134	!-----------------------------------------------------------------------
135	! Declaration of the variables required by correlated-k subroutines
136	! Numbered from top to bottom (unlike in the GCM)
137	!-----------------------------------------------------------------------
138
139	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
140	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
141
142	! Optical values for the optci/cv subroutines
143	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
144	! NB: Arrays below are "save" to avoid reallocating them at every call
145	! not because their content needs be reused from call to the next
146	REAL*8,allocatable,save :: dtaui(:,:,:)
147	REAL*8,allocatable,save :: dtauv(:,:,:)
148	REAL*8,allocatable,save :: cosbv(:,:,:)
149	REAL*8,allocatable,save :: cosbi(:,:,:)
150	REAL*8,allocatable,save :: wbari(:,:,:)
151	REAL*8,allocatable,save :: wbarv(:,:,:)
152	!$OMP THREADPRIVATE(dtaui,dtauv,cosbv,cosbi,wbari,wbarv)
153	REAL*8,allocatable,save :: tauv(:,:,:)
154	REAL*8,allocatable,save :: taucumv(:,:,:)
155	REAL*8,allocatable,save :: taucumi(:,:,:)
156	!$OMP THREADPRIVATE(tauv,taucumv,taucumi)
157	REAL*8,allocatable,save :: tauaero(:,:)
158	!$OMP THREADPRIVATE(tauaero)
159	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
160	real*8 nfluxtopv_nu(L_NSPECTV)
161	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
162	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
163	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
164	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
165	real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
166	real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
167	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
168	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
169	REAL*8 albi,acosz
170	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
171
172	INTEGER ig,l,k,nw,iaer,iq
173
174	real*8,allocatable,save :: taugsurf(:,:)
175	real*8,allocatable,save :: taugsurfi(:,:)
176	!$OMP THREADPRIVATE(taugsurf,taugsurfi)
177	real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol). index 1 is the top of the atmosphere, index L_LEVELS is the bottom
178
179	! Local aerosol optical properties for each column on RADIATIVE grid.
180	real8,save,allocatable :: QXVAER(:,:,:) ! Extinction coeff (QVISsQREFQREFvis)
181	real*8,save,allocatable :: QSVAER(:,:,:)
182	real*8,save,allocatable :: GVAER(:,:,:)
183	real8,save,allocatable :: QXIAER(:,:,:) ! Extinction coeff (QIRsQREFQREFir)
184	real*8,save,allocatable :: QSIAER(:,:,:)
185	real*8,save,allocatable :: GIAER(:,:,:)
186	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER)
187	real, dimension(:,:,:), save, allocatable :: QREFvis3d
188	real, dimension(:,:,:), save, allocatable :: QREFir3d
189	!$OMP THREADPRIVATE(QREFvis3d,QREFir3d)
190
191	! Miscellaneous :
192	real*8 temp,temp1,temp2,pweight,sig
193	character(len=10) :: tmp1
194	character(len=10) :: tmp2
195	character(len=100) :: message
196	character(len=10),parameter :: subname="callcorrk"
197
198	logical OLRz
199	real*8 NFLUXGNDV_nu(L_NSPECTV)
200
201	! Included by RW for runaway greenhouse 1D study.
202	real vtmp(nlayer)
203	REAL*8 muvarrad(L_LEVELS)
204
205	! Included by MT for albedo calculations.
206	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
207	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
208
209	! NLTE factor for CH4
210	real eps_nlte_sw23(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
211	real eps_nlte_sw33(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
212	real eps_nlte_lw(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw
213	integer Nfine,ifine
214	parameter(Nfine=701)
215	real,save :: levdat(Nfine),vmrdat(Nfine)
216	REAL dtlw_hcn_c2h2(ngrid, nlayer) ! cooling rate (K/s) due to C2H2/HCN (diagnostic)
217	real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix
218
219	REAL dtlw_nu(nlayer,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands
220	REAL dtsw_nu(nlayer,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands
221
222	! local variable
223	REAL dpp ! intermediate
224
225	integer ok ! status (returned by NetCDF functions)
226
227	REAL :: maxvalue,minvalue
228
229	!===============================================================
230	! I.a Initialization on first call
231	!===============================================================
232
233
234	if(firstcall) then
235
236	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
237	if(.not.allocated(QVISsQREF3d)) then
238	allocate(QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind))
239	endif
240	if(.not.allocated(omegaVIS3d)) then
241	allocate(omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind))
242	endif
243	if(.not.allocated(gVIS3d)) then
244	allocate(gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind))
245	endif
246	if (.not.allocated(QIRsQREF3d)) then
247	allocate(QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind))
248	endif
249	if (.not.allocated(omegaIR3d)) then
250	allocate(omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind))
251	endif
252	if (.not.allocated(gIR3d)) then
253	allocate(gIR3d(ngrid,nlayer,L_NSPECTI,naerkind))
254	endif
255	if (.not.allocated(tauaero)) then
256	allocate(tauaero(L_LEVELS,naerkind))
257	endif
258
259	if(.not.allocated(QXVAER)) then
260	allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
261	if (ok /= 0) then
262	write(,) "memory allocation failed for QXVAER!"
263	call abort_physic(subname,'allocation failure for QXVAER',1)
264	endif
265	endif
266	if(.not.allocated(QSVAER)) then
267	allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
268	if (ok /= 0) then
269	write(,) "memory allocation failed for QSVAER!"
270	call abort_physic(subname,'allocation failure for QSVAER',1)
271	endif
272	endif
273	if(.not.allocated(GVAER)) then
274	allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
275	if (ok /= 0) then
276	write(,) "memory allocation failed for GVAER!"
277	call abort_physic(subname,'allocation failure for GVAER',1)
278	endif
279	endif
280	if(.not.allocated(QXIAER)) then
281	allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
282	if (ok /= 0) then
283	write(,) "memory allocation failed for QXIAER!"
284	call abort_physic(subname,'allocation failure for QXIAER',1)
285	endif
286	endif
287	if(.not.allocated(QSIAER)) then
288	allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
289	if (ok /= 0) then
290	write(,) "memory allocation failed for QSIAER!"
291	call abort_physic(subname,'allocation failure for QSIAER',1)
292	endif
293	endif
294	if(.not.allocated(GIAER)) then
295	allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
296	if (ok /= 0) then
297	write(,) "memory allocation failed for GIAER!"
298	call abort_physic(subname,'allocation failure for GIAER',1)
299	endif
300	endif
301
302	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
303	IF(.not.ALLOCATED(QREFvis3d))THEN
304	ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind), stat=ok)
305	IF (ok/=0) THEN
306	write(,) "memory allocation failed for QREFvis3d!"
307	call abort_physic(subname,'allocation failure for QREFvis3d',1)
308	ENDIF
309	ENDIF
310	IF(.not.ALLOCATED(QREFir3d)) THEN
311	ALLOCATE(QREFir3d(ngrid,nlayer,naerkind), stat=ok)
312	IF (ok/=0) THEN
313	write(,) "memory allocation failed for QREFir3d!"
314	call abort_physic(subname,'allocation failure for QREFir3d',1)
315	ENDIF
316	ENDIF
317	! Effective radius and variance of the aerosols
318	IF(.not.ALLOCATED(reffrad)) THEN
319	allocate(reffrad(ngrid,nlayer,naerkind), stat=ok)
320	IF (ok/=0) THEN
321	write(,) "memory allocation failed for reffrad!"
322	call abort_physic(subname,'allocation failure for reffrad',1)
323	ENDIF
324	ENDIF
325	IF(.not.ALLOCATED(nueffrad)) THEN
326	allocate(nueffrad(ngrid,nlayer,naerkind), stat=ok)
327	IF (ok/=0) THEN
328	write(,) "memory allocation failed for nueffrad!"
329	call abort_physic(subname,'allocation failure for nueffrad',1)
330	ENDIF
331	ENDIF
332
333	if (is_master) call system('rm -f surf_vals_long.out')
334
335	!--------------------------------------------------
336	! Set up correlated k
337	!--------------------------------------------------
338
339	!this block is now done at firstcall of physiq_mod
340	! print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
341	! call getin_p("corrkdir",corrkdir)
342	! print*, "corrkdir = ",corrkdir
343	! write( tmp1, '(i3)' ) L_NSPECTI
344	! write( tmp2, '(i3)' ) L_NSPECTV
345	! banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
346	! banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
347
348	! call setspi ! Basic infrared properties.
349	! call setspv ! Basic visible properties.
350	! call sugas_corrk ! Set up gaseous absorption properties.
351	! call suaer_corrk ! Set up aerosol optical properties.
352
353
354	! now that L_NGAUSS has been initialized (by sugas_corrk)
355	! allocate related arrays
356	if(.not.allocated(dtaui)) then
357	ALLOCATE(dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok)
358	if (ok/=0) then
359	write(,) "memory allocation failed for dtaui!"
360	call abort_physic(subname,'allocation failure for dtaui',1)
361	endif
362	endif
363	if(.not.allocated(dtauv)) then
364	ALLOCATE(dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok)
365	if (ok/=0) then
366	write(,) "memory allocation failed for dtauv!"
367	call abort_physic(subname,'allocation failure for dtauv',1)
368	endif
369	endif
370	if(.not.allocated(cosbv)) then
371	ALLOCATE(cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok)
372	if (ok/=0) then
373	write(,) "memory allocation failed for cosbv!"
374	call abort_physic(subname,'allocation failure for cobsv',1)
375	endif
376	endif
377	if(.not.allocated(cosbi)) then
378	ALLOCATE(cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok)
379	if (ok/=0) then
380	write(,) "memory allocation failed for cosbi!"
381	call abort_physic(subname,'allocation failure for cobsi',1)
382	endif
383	endif
384	if(.not.allocated(wbari)) then
385	ALLOCATE(wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok)
386	if (ok/=0) then
387	write(,) "memory allocation failed for wbari!"
388	call abort_physic(subname,'allocation failure for wbari',1)
389	endif
390	endif
391	if(.not.allocated(wbarv)) then
392	ALLOCATE(wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok)
393	if (ok/=0) then
394	write(,) "memory allocation failed for wbarv!"
395	call abort_physic(subname,'allocation failure for wbarv',1)
396	endif
397	endif
398	if(.not.allocated(tauv)) then
399	ALLOCATE(tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS), stat=ok)
400	if (ok/=0) then
401	write(,) "memory allocation failed for tauv!"
402	call abort_physic(subname,'allocation failure for tauv',1)
403	endif
404	endif
405	if(.not.allocated(taucumv)) then
406	ALLOCATE(taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS), stat=ok)
407	if (ok/=0) then
408	write(,) "memory allocation failed for taucumv!"
409	call abort_physic(subname,'allocation failure for taucumv',1)
410	endif
411	endif
412	if(.not.allocated(taucumi)) then
413	ALLOCATE(taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS), stat=ok)
414	if (ok/=0) then
415	write(,) "memory allocation failed for taucumi!"
416	call abort_physic(subname,'allocation failure for taucumi',1)
417	endif
418	endif
419	if(.not.allocated(taugsurf)) then
420	ALLOCATE(taugsurf(L_NSPECTV,L_NGAUSS-1), stat=ok)
421	if (ok/=0) then
422	write(,) "memory allocation failed for taugsurf!"
423	call abort_physic(subname,'allocation failure for taugsurf',1)
424	endif
425	endif
426	if(.not.allocated(taugsurfi)) then
427	ALLOCATE(taugsurfi(L_NSPECTI,L_NGAUSS-1), stat=ok)
428	if (ok/=0) then
429	write(,) "memory allocation failed for taugsurfi!"
430	call abort_physic(subname,'allocation failure for taugsurfi',1)
431	endif
432	endif
433
434	end if ! of if (firstcall)
435
436	!=======================================================================
437	! I.b Initialization on every call
438	!=======================================================================
439
440	qxvaer(:,:,:)=0.0
441	qsvaer(:,:,:)=0.0
442	gvaer(:,:,:) =0.0
443
444	qxiaer(:,:,:)=0.0
445	qsiaer(:,:,:)=0.0
446	giaer(:,:,:) =0.0
447
448	OLR_nu(:,:) = 0.
449	OSR_nu(:,:) = 0.
450	GSR_nu(:,:) = 0.
451
452	!--------------------------------------------------
453	! Effective radius and variance of the aerosols
454	!--------------------------------------------------
455	! Radiative Hazes
456	if (optichaze) then
457	if (callmufi) then
458	! Spherical aerosols
459	sig = 0.2
460	where (mp2m_rc_sph(:,:) > 1e-10)
461	reffrad(:,:,1) = mp2m_rc_sph(:,:) * exp(5.sig*2 / 2.)
462	elsewhere
463	reffrad(:,:,1) = 0d0
464	endwhere
465	if (exp(sig**2) - 1 > 0.1) then
466	nueffrad(:,:,1) = exp(sig**2) - 1
467	else
468	nueffrad(:,:,1) = 0.1
469	endif
470	! Fractal aerosols
471	sig = 0.35
472	where (mp2m_rc_fra(:,:) > 1e-10)
473	reffrad(:,:,2) = mp2m_rc_fra(:,:) * exp(5.sig*2 / 2.)
474	elsewhere
475	reffrad(:,:,2) = 0d0
476	endwhere
477	if (exp(sig**2) - 1 > 0.1) then
478	nueffrad(:,:,2) = exp(sig**2) - 1
479	else
480	nueffrad(:,:,2) = 0.1
481	endif
482
483	else
484	do iaer=1,naerkind
485	if ((iaer.eq.iaero_haze)) then
486	call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer),nueffrad(1,1,iaer))
487	endif
488	end do ! iaer = 1, naerkind.
489	if (haze_radproffix) then
490	call haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad)
491	if (is_master) print*, 'haze_radproffix=T : fixed profile for haze rad'
492	else
493	if (is_master) print*,'reffrad haze:',reffrad(1,1,iaero_haze)
494	if (is_master) print*,'nueff haze',nueffrad(1,1,iaero_haze)
495	endif ! end haze_radproffix
496	endif ! end callmufi
497	endif ! end radiative haze
498
499	! How much light do we get ?
500	do nw=1,L_NSPECTV
501	stel(nw)=stellarf(nw)/(dist_star**2)
502	end do
503
504	if (optichaze) then
505
506	! Get 3D aerosol optical properties.
507	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
508	QVISsQREF3d,omegaVIS3d,gVIS3d, &
509	QIRsQREF3d,omegaIR3d,gIR3d, &
510	QREFvis3d,QREFir3d)
511
512	! Get aerosol optical depths.
513	call aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq,dtau_aer, &
514	reffrad,nueffrad,QREFvis3d,QREFir3d, &
515	tau_col)
516	endif
517
518
519	!-----------------------------------------------------------------------
520	! Prepare CH4 mixing ratio for radiative transfer
521	IF (methane) then
522	vmrch4(:,:)=0.
523
524	if (ch4fix) then
525	if (vmrch4_proffix) then
526	!! Interpolate on the model vertical grid
527	do ig=1,ngrid
528	CALL interp_line(levdat,vmrdat,Nfine, &
529	zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
530	enddo
531	else
532	vmrch4(:,:)=vmrch4fix
533	endif
534	else
535	vmrch4(:,:)=pq(:,:,igcm_ch4_gas)100. &
536	mmol(igcm_n2)/mmol(igcm_ch4_gas)
537	endif
538	ENDIF
539
540	! Prepare NON LTE correction in Pluto atmosphere
541	IF (nlte) then
542	CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,&
543	eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
544	ENDIF
545	! Net atmospheric radiative cooling rate from C2H2 (K.s-1):
546	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547	! dtlw_hcn_c2h2=0.
548	if (cooling) then
549	CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,&
550	pt,dtlw_hcn_c2h2)
551	endif
552
553
554	!-----------------------------------------------------------------------
555	do ig=1,ngrid ! Starting Big Loop over every GCM column
556	!-----------------------------------------------------------------------
557
558
559	!=======================================================================
560	! II. Transformation of the GCM variables
561	!=======================================================================
562
563
564	!-----------------------------------------------------------------------
565	! Aerosol optical properties Qext, Qscat and g.
566	! The transformation in the vertical is the same as for temperature.
567	!-----------------------------------------------------------------------
568
569
570	! AF24: for now only consider one aerosol (=haze)
571	if (optichaze) then
572	do iaer=1,naerkind
573	! Shortwave.
574	do nw=1,L_NSPECTV
575
576	do l=1,nlayer
577
578	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
579	*QREFvis3d(ig,nlayer+1-l,iaer)
580
581	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
582	*QREFvis3d(ig,max(nlayer-l,1),iaer)
583
584	qxvaer(2*l,nw,iaer) = temp1
585	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
586
587	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
588	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
589
590	qsvaer(2*l,nw,iaer) = temp1
591	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
592
593	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
594	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
595
596	gvaer(2*l,nw,iaer) = temp1
597	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
598
599	end do ! nlayer
600
601	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
602	qxvaer(2*nlayer+1,nw,iaer)=0.
603
604	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
605	qsvaer(2*nlayer+1,nw,iaer)=0.
606
607	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
608	gvaer(2*nlayer+1,nw,iaer)=0.
609
610	end do ! L_NSPECTV
611
612	do nw=1,L_NSPECTI
613	! Longwave
614	do l=1,nlayer
615
616	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
617	*QREFir3d(ig,nlayer+1-l,iaer)
618
619	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
620	*QREFir3d(ig,max(nlayer-l,1),iaer)
621
622	qxiaer(2*l,nw,iaer) = temp1
623	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
624
625	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
626	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
627
628	qsiaer(2*l,nw,iaer) = temp1
629	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
630
631	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
632	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
633
634	giaer(2*l,nw,iaer) = temp1
635	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
636
637	end do ! nlayer
638
639	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
640	qxiaer(2*nlayer+1,nw,iaer)=0.
641
642	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
643	qsiaer(2*nlayer+1,nw,iaer)=0.
644
645	giaer(1,nw,iaer)=giaer(2,nw,iaer)
646	giaer(2*nlayer+1,nw,iaer)=0.
647
648	end do ! L_NSPECTI
649
650	end do ! naerkind
651
652	! Test / Correct for freaky s. s. albedo values.
653	do iaer=1,naerkind
654	do k=1,L_LEVELS
655
656	do nw=1,L_NSPECTV
657	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
658	message='Serious problems with qsvaer values'
659	call abort_physic(subname,message,1)
660	endif
661	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
662	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
663	endif
664	end do
665
666	do nw=1,L_NSPECTI
667	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
668	message='Serious problems with qsvaer values'
669	call abort_physic(subname,message,1)
670	endif
671	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
672	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
673	endif
674	end do
675
676	end do ! L_LEVELS
677	end do ! naerkind
678	end if ! optichaze
679
680	!-----------------------------------------------------------------------
681	! Aerosol optical depths
682	!-----------------------------------------------------------------------
683	if (optichaze) then
684	do iaer=1,naerkind ! a bug was here
685	do k=0,nlayer-1
686
687	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
688	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
689	! As 'dtau_aer' is at reference (visible) wavelenght we scale it as
690	! it will be multplied by qxi/v in optci/v
691	temp=dtau_aer(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
692	tauaero(2k+2,iaer)=max(temppweight,0.d0)
693	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
694
695	end do
696	! boundary conditions
697	tauaero(1,iaer) = tauaero(2,iaer)
698	!tauaero(1,iaer) = 0.
699	!JL18 at time of testing, the two above conditions gave the same results bit for bit.
700
701	end do ! naerkind
702	else
703	tauaero(:,:)=0
704	end if ! optichaze
705
706	! Albedo and Emissivity.
707	albi=1-emis(ig) ! Long Wave.
708	DO nw=1,L_NSPECTV ! Short Wave loop.
709	albv(nw)=albedo(ig,nw)
710	ENDDO
711
712	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
713
714	!-----------------------------------------------------------------------
715	! Methane Vapor
716	!-----------------------------------------------------------------------
717	if (methane) then
718	do l=1,nlayer
719	qvar(2l) = vmrch4(ig,nlayer+1-l)/100. &
720	mmol(igcm_ch4_gas)/mmol(igcm_n2)
721	qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig, &
722	max(nlayer-l,1)))/2.)/100.* &
723	mmol(igcm_ch4_gas)/mmol(igcm_n2)
724	end do
725	qvar(1)=qvar(2)
726
727	else
728	do k=1,L_LEVELS
729	qvar(k) = 1.0D-7
730	end do
731	end if ! end methane
732
733	!-----------------------------------------------------------------------
734	! kcm mode only !
735	!-----------------------------------------------------------------------
736
737	DO l=1,nlayer
738	muvarrad(2*l) = muvar(ig,nlayer+2-l)
739	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
740	END DO
741
742	muvarrad(1) = muvarrad(2)
743	muvarrad(2*nlayer+1)=muvar(ig,1)
744
745	! Keep values inside limits for which we have radiative transfer coefficients !!!
746	if(L_REFVAR.gt.1)then ! (there was a bug here)
747	do k=1,L_LEVELS
748	if(qvar(k).lt.wrefvar(1))then
749	qvar(k)=wrefvar(1)+1.0e-8
750	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
751	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
752	endif
753	end do
754	endif
755
756	!-----------------------------------------------------------------------
757	! Pressure and temperature
758	!-----------------------------------------------------------------------
759
760	DO l=1,nlayer
761	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
762	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
763	tlevrad(2*l) = pt(ig,nlayer+1-l)
764	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
765	END DO
766
767	plevrad(1) = 0.
768	!!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.
769
770	tlevrad(1) = tlevrad(2)
771	tlevrad(2*nlayer+1)=tsurf(ig)
772
773	pmid(1) = pplay(ig,nlayer)/scalep
774	pmid(2) = pmid(1)
775
776	tmid(1) = tlevrad(2)
777	tmid(2) = tmid(1)
778
779	DO l=1,L_NLAYRAD-1
780	tmid(2l+1) = tlevrad(2l+1)
781	tmid(2l+2) = tlevrad(2l+1)
782	pmid(2l+1) = plevrad(2l+1)
783	pmid(2l+2) = plevrad(2l+1)
784	END DO
785	pmid(L_LEVELS) = plevrad(L_LEVELS)
786	tmid(L_LEVELS) = tlevrad(L_LEVELS)
787
788	!!Alternative interpolation:
789	! pmid(3) = pmid(1)
790	! pmid(4) = pmid(1)
791	! tmid(3) = tmid(1)
792	! tmid(4) = tmid(1)
793	! DO l=2,L_NLAYRAD-1
794	! tmid(2l+1) = tlevrad(2l)
795	! tmid(2l+2) = tlevrad(2l)
796	! pmid(2l+1) = plevrad(2l)
797	! pmid(2l+2) = plevrad(2l)
798	! END DO
799	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
800	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
801
802	! Test for out-of-bounds pressure.
803	if(plevrad(3).lt.pgasmin)then
804	print*,'Warning: minimum pressure is outside the radiative'
805	print*,'transfer kmatrix bounds, exiting.'
806	print*,'Pressure:', plevrad(3), 'Pa'
807	message="Minimum pressure outside of kmatrix bounds"
808	!call abort_physic(subname,message,1)
809	elseif(plevrad(L_LEVELS).gt.pgasmax)then
810	print*,'Maximum pressure is outside the radiative'
811	print*,'transfer kmatrix bounds, exiting.'
812	message="Minimum pressure outside of kmatrix bounds"
813	call abort_physic(subname,message,1)
814	endif
815
816	! Test for out-of-bounds temperature.
817	! -- JVO 20 : Also add a sanity test checking that tlevrad is
818	! within Planck function temperature boundaries,
819	! which would cause gfluxi/sfluxi to crash.
820	do k=1,L_LEVELS
821
822	if(tlevrad(k).lt.tgasmin)then
823	print*,'Minimum temperature is outside the radiative'
824	print*,'transfer kmatrix bounds'
825	print*,"k=",k," tlevrad(k)=",tlevrad(k)
826	print*,"tgasmin=",tgasmin
827	if (strictboundcorrk) then
828	message="Minimum temperature outside of kmatrix bounds"
829	call abort_physic(subname,message,1)
830	else
831	print,'**********************************************'
832	print*,'we allow model to continue with tlevrad<tgasmin'
833	print*,' ... we assume we know what you are doing ... '
834	print*,' ... but do not let this happen too often ... '
835	print,'**********************************************'
836	!tlevrad(k)=tgasmin ! Used in the source function !
837	endif
838	elseif(tlevrad(k).gt.tgasmax)then
839	print*,'Maximum temperature is outside the radiative'
840	print*,'transfer kmatrix bounds, exiting.'
841	print*,"k=",k," tlevrad(k)=",tlevrad(k)
842	print*,"tgasmax=",tgasmax
843	if (strictboundcorrk) then
844	message="Maximum temperature outside of kmatrix bounds"
845	call abort_physic(subname,message,1)
846	else
847	print,'**********************************************'
848	print*,'we allow model to continue with tlevrad>tgasmax'
849	print*,' ... we assume we know what you are doing ... '
850	print*,' ... but do not let this happen too often ... '
851	print,'**********************************************'
852	!tlevrad(k)=tgasmax ! Used in the source function !
853	endif
854	endif
855
856	if (tlevrad(k).lt.tplanckmin) then
857	print*,'Minimum temperature is outside the boundaries for'
858	print*,'Planck function integration set in callphys.def, aborting.'
859	print*,"k=",k," tlevrad(k)=",tlevrad(k)
860	print*,"tplanckmin=",tplanckmin
861	message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
862	call abort_physic(subname,message,1)
863	else if (tlevrad(k).gt.tplanckmax) then
864	print*,'Maximum temperature is outside the boundaries for'
865	print*,'Planck function integration set in callphys.def, aborting.'
866	print*,"k=",k," tlevrad(k)=",tlevrad(k)
867	print*,"tplanckmax=",tplanckmax
868	message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
869	call abort_physic(subname,message,1)
870	endif
871
872	enddo
873
874	do k=1,L_NLAYRAD+1
875	if(tmid(k).lt.tgasmin)then
876	print*,'Minimum temperature is outside the radiative'
877	print*,'transfer kmatrix bounds, exiting.'
878	print*,"k=",k," tmid(k)=",tmid(k)
879	print*,"tgasmin=",tgasmin
880	if (strictboundcorrk) then
881	message="Minimum temperature outside of kmatrix bounds"
882	call abort_physic(subname,message,1)
883	else
884	print,'**********************************************'
885	print*,'we allow model to continue but with tmid=tgasmin'
886	print*,' ... we assume we know what you are doing ... '
887	print*,' ... but do not let this happen too often ... '
888	print,'**********************************************'
889	tmid(k)=tgasmin
890	endif
891	elseif(tmid(k).gt.tgasmax)then
892	print*,'Maximum temperature is outside the radiative'
893	print*,'transfer kmatrix bounds, exiting.'
894	print*,"k=",k," tmid(k)=",tmid(k)
895	print*,"tgasmax=",tgasmax
896	if (strictboundcorrk) then
897	message="Maximum temperature outside of kmatrix bounds"
898	call abort_physic(subname,message,1)
899	else
900	print,'**********************************************'
901	print*,'we allow model to continue but with tmid=tgasmax'
902	print*,' ... we assume we know what you are doing ... '
903	print*,' ... but do not let this happen too often ... '
904	print,'**********************************************'
905	tmid(k)=tgasmax
906	endif
907	endif
908	enddo
909
910	!=======================================================================
911	! III. Calling the main radiative transfer subroutines
912	!=======================================================================
913
914	! ----------------------------------------------------------------
915	! Recombine reference corrk tables if needed - Added by JVO, 2020.
916	if (corrk_recombin) then
917	call call_recombin(ig,nlayer,pq(ig,:,:),pplay(ig,:),pt(ig,:),qvar(:),tmid(:),pmid(:))
918	endif
919	! ----------------------------------------------------------------
920
921	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
922	glat_ig=glat(ig)
923
924	!-----------------------------------------------------------------------
925	! Short Wave Part
926	!-----------------------------------------------------------------------
927
928	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
929	if((ngrid.eq.1).and.(global1d))then
930	do nw=1,L_NSPECTV
931	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
932	end do
933	else
934	do nw=1,L_NSPECTV
935	stel_fract(nw)= stel(nw) * fract(ig)
936	end do
937	endif
938
939	call optcv(dtauv,tauv,taucumv,plevrad, &
940	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
941	tmid,pmid,taugsurf,qvar,muvarrad)
942
943	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
944	acosz,stel_fract,nfluxtopv,fluxtopvdn,nfluxoutv_nu,&
945	nfluxgndv_nu,nfluxtopv_nu, &
946	fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf)
947
948	else ! During the night, fluxes = 0.
949	nfluxtopv = 0.0d0
950	fluxtopvdn = 0.0d0
951	nfluxoutv_nu(:) = 0.0d0
952	nfluxgndv_nu(:) = 0.0d0
953	do l=1,L_NLAYRAD
954	fmnetv(l)=0.0d0
955	fmnetv_nu(l,:)=0.0d0
956	fluxupv(l)=0.0d0
957	fluxdnv(l)=0.0d0
958	end do
959	end if
960
961
962	! Equivalent Albedo Calculation (for OUTPUT). MT2015
963	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
964	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
965	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
966	DO nw=1,L_NSPECTV
967	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
968	ENDDO
969	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
970	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
971	else
972	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
973	endif
974	else
975	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
976	endif
977
978
979	!-----------------------------------------------------------------------
980	! Long Wave Part
981	!-----------------------------------------------------------------------
982
983	call optci(plevrad,tlevrad,dtaui,taucumi, &
984	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
985	taugsurfi,qvar,muvarrad)
986
987	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
988	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
989	fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
990
991	!-----------------------------------------------------------------------
992	! Transformation of the correlated-k code outputs
993	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
994
995	! Flux incident at the top of the atmosphere
996	fluxtop_dn(ig)=fluxtopvdn
997
998	fluxtop_lw(ig) = real(nfluxtopi)
999	fluxabs_sw(ig) = real(-nfluxtopv)
1000	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
1001	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
1002
1003	! Flux absorbed by the surface. By MT2015.
1004	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
1005
1006	if(fluxtop_dn(ig).lt.0.0)then
1007	print*,'Achtung! fluxtop_dn has lost the plot!'
1008	print*,'fluxtop_dn=',fluxtop_dn(ig)
1009	print*,'acosz=',acosz
1010	print*,'dtau_aer=',dtau_aer(ig,:,:)
1011	print*,'temp= ',pt(ig,:)
1012	print*,'pplay= ',pplay(ig,:)
1013	message="Achtung! fluxtop_dn has lost the plot!"
1014	call abort_physic(subname,message,1)
1015	endif
1016
1017	! Spectral output, for exoplanet observational comparison
1018	if(specOLR)then
1019	do nw=1,L_NSPECTI
1020	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
1021	end do
1022	do nw=1,L_NSPECTV
1023	GSR_nu(ig,nw)=nfluxgndv_nu(nw)/DWNV(nw)
1024	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
1025	end do
1026	endif
1027
1028	! Finally, the heating rates
1029	DO l=2,L_NLAYRAD
1030	! dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
1031	! glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
1032	dpp = glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2l-1)))
1033	do nw=1,L_NSPECTV
1034	dtsw_nu(L_NLAYRAD+1-l,nw)= &
1035	(fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
1036	end do
1037
1038	! dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
1039	! glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
1040	do nw=1,L_NSPECTI
1041	dtlw_nu(L_NLAYRAD+1-l,nw)= &
1042	(fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
1043	end do
1044	END DO
1045
1046	! These are values at top of atmosphere
1047	! dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
1048	! glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
1049	! dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
1050	! glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
1051	dpp = g/(cppscalep(plevrad(3)-plevrad(1)))
1052	do nw=1,L_NSPECTV
1053	dtsw_nu(L_NLAYRAD,nw)= &
1054	(fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp
1055	end do
1056	do nw=1,L_NSPECTI
1057	dtlw_nu(L_NLAYRAD,nw)= &
1058	(fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp
1059	end do
1060
1061	! Optical thickness diagnostics (added by JVO)
1062	if (diagdtau) then
1063	do l=1,L_NLAYRAD
1064	do nw=1,L_NSPECTV
1065	int_dtauv(ig,l,nw) = 0.0d0
1066	DO k=1,L_NGAUSS
1067	! Output exp(-tau) because gweight ponderates exp and not tau itself
1068	int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
1069	ENDDO
1070	enddo
1071	do nw=1,L_NSPECTI
1072	int_dtaui(ig,l,nw) = 0.0d0
1073	DO k=1,L_NGAUSS
1074	! Output exp(-tau) because gweight ponderates exp and not tau itself
1075	int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
1076	ENDDO
1077	enddo
1078	enddo
1079	endif
1080
1081	! **********************************************************
1082	! NON NLTE correction in Pluto atmosphere
1083	! And conversion of LW spectral heating rates to total rates
1084	! **********************************************************
1085
1086	if (.not.nlte) then
1087	eps_nlte_sw23(ig,:) =1. ! IF no NLTE
1088	eps_nlte_sw33(ig,:) =1. ! IF no NLTE
1089	eps_nlte_lw(ig,:) =1. ! IF no NLTE
1090	endif
1091
1092	do l=1,nlayer
1093
1094	!LW
1095	dtlw(ig,l) =0.
1096	! dtlw_co(ig,l) =0. ! only for diagnostic
1097	do nw=1,L_NSPECTI
1098	! wewei : wavelength in micrometer
1099	if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
1100	dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
1101	else
1102	!dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
1103	dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
1104	! dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
1105	end if
1106	dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
1107	end do
1108	! adding c2h2 if cooling active
1109	! dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l)
1110
1111	!SW
1112	dtsw(ig,l) =0.
1113
1114	if (strobel) then
1115
1116	do nw=1,L_NSPECTV
1117	if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
1118	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
1119	elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
1120	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
1121	else
1122	dtsw_nu(l,nw)=dtsw_nu(l,nw)
1123	end if
1124	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
1125	end do
1126
1127	else ! total heating rates multiplied by nlte coef
1128
1129	do nw=1,L_NSPECTV
1130	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
1131	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
1132	enddo
1133
1134	endif
1135
1136
1137	end do
1138	! **********************************************************
1139
1140
1141	!-----------------------------------------------------------------------
1142	end do ! End of big loop over every GCM column.
1143	!-----------------------------------------------------------------------
1144
1145
1146	!-----------------------------------------------------------------------
1147	! Additional diagnostics
1148	!-----------------------------------------------------------------------
1149
1150	! IR spectral output, for exoplanet observational comparison
1151	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
1152
1153	print*,'Saving scalar quantities in surf_vals.out...'
1154	print*,'psurf = ', pplev(1,1),' Pa'
1155	open(116,file='surf_vals.out')
1156	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
1157	real(-nfluxtopv),real(nfluxtopi)
1158	close(116)
1159
1160
1161	! USEFUL COMMENT - Do Not Remove.
1162	!
1163	! if(specOLR)then
1164	! open(117,file='OLRnu.out')
1165	! do nw=1,L_NSPECTI
1166	! write(117,*) OLR_nu(1,nw)
1167	! enddo
1168	! close(117)
1169	!
1170	! open(127,file='OSRnu.out')
1171	! do nw=1,L_NSPECTV
1172	! write(127,*) OSR_nu(1,nw)
1173	! enddo
1174	! close(127)
1175	! endif
1176
1177	! OLR vs altitude: do it as a .txt file.
1178	OLRz=.false.
1179	if(OLRz)then
1180	print*,'saving IR vertical flux for OLRz...'
1181	open(118,file='OLRz_plevs.out')
1182	open(119,file='OLRz.out')
1183	do l=1,L_NLAYRAD
1184	write(118,) plevrad(2l)
1185	do nw=1,L_NSPECTI
1186	write(119,*) fluxupi_nu(l,nw)
1187	enddo
1188	enddo
1189	close(118)
1190	close(119)
1191	endif
1192
1193	endif
1194
1195	! See physiq.F for explanations about CLFvarying. This is temporary.
1196	if (lastcall) then
1197	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
1198	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
1199	!$OMP BARRIER
1200	!$OMP MASTER
1201	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
1202	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
1203	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
1204	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
1205	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
1206	!$OMP END MASTER
1207	!$OMP BARRIER
1208	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
1209	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
1210	endif
1211
1212
1213	end subroutine callcorrk
1214
1215	END MODULE callcorrk_mod

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

Download in other formats:

Original Format