Context Navigation

callcorrk.F90 @ 3380

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

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

Download in other formats:

Original Format