Context Navigation

callcorrk.F90 @ 3982

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

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

Download in other formats:

Original Format