Context Navigation

callcorrk.F90 @ 3572

Last change on this file since 3572 was 3572, checked in by debatzbr, 3 weeks ago
Remove generic_aerosols and generic_condensation, along with their related variables (useless). RENAME THE VARIABLE AEROHAZE TO OPTICHAZE.
Property svn:executable set to ``*
File size: 48.0 KB

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

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

Download in other formats:

Original Format