Context Navigation

callcorrk.F90 @ 3210

Last change on this file since 3210 was 3197, checked in by afalco, 22 months ago
Pluto PCM: Fixed some flags for aerosols AF
Property svn:executable set to ``*
File size: 47.4 KB

Line
1	MODULE callcorrk_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
8	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
9	tsurf,fract,dist_star,aerosol,muvar, &
10	dtlw,dtsw,fluxsurf_lw, &
11	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
12	fluxabs_sw,fluxtop_dn, &
13	OLR_nu,OSR_nu,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
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
33	use comcstfi_mod, only: pi, mugaz, cpp
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	use optcv_mod, only: optcv
39	use optci_mod, only: optci
40	use sfluxi_mod, only: sfluxi
41	use sfluxv_mod, only: sfluxv
42	use recombin_corrk_mod, only: corrk_recombin, call_recombin
43	use generic_cloud_common_h, only: Psat_generic, epsi_generic
44	use generic_tracer_index_mod, only: generic_tracer_index
45	use planetwide_mod, only: planetwide_maxval, planetwide_minval
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) :: tsurf(ngrid) ! Surface temperature (K).
83	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
84	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
85	REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
86	REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer) ! Fraction of clouds (%).
87	logical,intent(in) :: clearsky
88	logical,intent(in) :: firstcall ! Signals first call to physics.
89	logical,intent(in) :: lastcall ! Signals last call to physics.
90
91	! OUTPUT
92	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
93	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
94	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
95	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
96	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
97	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
98	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
99	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
100	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
101	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
102	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
103	REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV) ! Surface SW radiation in each band (Normalized to the band width (W/m2/cm-1).
104	REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity.
105	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
106	REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! Column Fraction of clouds (%).
107	REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
108	REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
109
110
111
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 :: cosbv(:,:,:)
148	REAL*8,allocatable,save :: cosbi(:,:,:)
149	REAL*8,allocatable,save :: wbari(:,:,:)
150	REAL*8,allocatable,save :: wbarv(:,:,:)
151	!$OMP THREADPRIVATE(dtaui,dtauv,cosbv,cosbi,wbari,wbarv)
152	REAL*8,allocatable,save :: tauv(:,:,:)
153	REAL*8,allocatable,save :: taucumv(:,:,:)
154	REAL*8,allocatable,save :: taucumi(:,:,:)
155	!$OMP THREADPRIVATE(tauv,taucumv,taucumi)
156	REAL*8,allocatable,save :: tauaero(:,:)
157	!$OMP THREADPRIVATE(tauaero)
158	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
159	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
160	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
161	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
162	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
163	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
164	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
165	REAL*8 albi,acosz
166	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
167
168	INTEGER ig,l,k,nw,iaer,iq
169
170	real*8,allocatable,save :: taugsurf(:,:)
171	real*8,allocatable,save :: taugsurfi(:,:)
172	!$OMP THREADPRIVATE(taugsurf,taugsurfi)
173	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
174
175	! Local aerosol optical properties for each column on RADIATIVE grid.
176	real8,save,allocatable :: QXVAER(:,:,:) ! Extinction coeff (QVISsQREFQREFvis)
177	real*8,save,allocatable :: QSVAER(:,:,:)
178	real*8,save,allocatable :: GVAER(:,:,:)
179	real8,save,allocatable :: QXIAER(:,:,:) ! Extinction coeff (QIRsQREFQREFir)
180	real*8,save,allocatable :: QSIAER(:,:,:)
181	real*8,save,allocatable :: GIAER(:,:,:)
182	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER)
183	real, dimension(:,:,:), save, allocatable :: QREFvis3d
184	real, dimension(:,:,:), save, allocatable :: QREFir3d
185	!$OMP THREADPRIVATE(QREFvis3d,QREFir3d)
186
187
188	! Miscellaneous :
189	real*8 temp,temp1,temp2,pweight
190	character(len=10) :: tmp1
191	character(len=10) :: tmp2
192	character(len=100) :: message
193	character(len=10),parameter :: subname="callcorrk"
194
195	! For fixed water vapour profiles.
196	integer i_var
197	real RH
198	real*8 pq_temp(nlayer)
199	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
200	real psat,qsat
201
202	logical OLRz
203	real*8 NFLUXGNDV_nu(L_NSPECTV)
204
205	! Included by RW for runaway greenhouse 1D study.
206	real vtmp(nlayer)
207	REAL*8 muvarrad(L_LEVELS)
208
209	! Included by MT for albedo calculations.
210	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
211	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
212
213	! local variable
214	integer ok ! status (returned by NetCDF functions)
215
216	integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
217	logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
218	real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
219	!$OMP THREADPRIVATE(metallicity)
220	REAL, SAVE :: qvap_deep ! deep mixing ratio of water vapor when simulating bottom less planets
221	!$OMP THREADPRIVATE(qvap_deep)
222
223	REAL :: maxvalue,minvalue
224
225	!===============================================================
226	! I.a Initialization on first call
227	!===============================================================
228
229
230	if(firstcall) then
231
232	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
233	if(.not.allocated(QVISsQREF3d)) then
234	allocate(QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind))
235	endif
236	if(.not.allocated(omegaVIS3d)) then
237	allocate(omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind))
238	endif
239	if(.not.allocated(gVIS3d)) then
240	allocate(gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind))
241	endif
242	if (.not.allocated(QIRsQREF3d)) then
243	allocate(QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind))
244	endif
245	if (.not.allocated(omegaIR3d)) then
246	allocate(omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind))
247	endif
248	if (.not.allocated(gIR3d)) then
249	allocate(gIR3d(ngrid,nlayer,L_NSPECTI,naerkind))
250	endif
251	if (.not.allocated(tauaero)) then
252	allocate(tauaero(L_LEVELS,naerkind))
253	endif
254
255	if(.not.allocated(QXVAER)) then
256	allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
257	if (ok /= 0) then
258	write(,) "memory allocation failed for QXVAER!"
259	call abort_physic(subname,'allocation failure for QXVAER',1)
260	endif
261	endif
262	if(.not.allocated(QSVAER)) then
263	allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
264	if (ok /= 0) then
265	write(,) "memory allocation failed for QSVAER!"
266	call abort_physic(subname,'allocation failure for QSVAER',1)
267	endif
268	endif
269	if(.not.allocated(GVAER)) then
270	allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
271	if (ok /= 0) then
272	write(,) "memory allocation failed for GVAER!"
273	call abort_physic(subname,'allocation failure for GVAER',1)
274	endif
275	endif
276	if(.not.allocated(QXIAER)) then
277	allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
278	if (ok /= 0) then
279	write(,) "memory allocation failed for QXIAER!"
280	call abort_physic(subname,'allocation failure for QXIAER',1)
281	endif
282	endif
283	if(.not.allocated(QSIAER)) then
284	allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
285	if (ok /= 0) then
286	write(,) "memory allocation failed for QSIAER!"
287	call abort_physic(subname,'allocation failure for QSIAER',1)
288	endif
289	endif
290	if(.not.allocated(GIAER)) then
291	allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok)
292	if (ok /= 0) then
293	write(,) "memory allocation failed for GIAER!"
294	call abort_physic(subname,'allocation failure for GIAER',1)
295	endif
296	endif
297
298	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
299	IF(.not.ALLOCATED(QREFvis3d))THEN
300	ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind), stat=ok)
301	IF (ok/=0) THEN
302	write(,) "memory allocation failed for QREFvis3d!"
303	call abort_physic(subname,'allocation failure for QREFvis3d',1)
304	ENDIF
305	ENDIF
306	IF(.not.ALLOCATED(QREFir3d)) THEN
307	ALLOCATE(QREFir3d(ngrid,nlayer,naerkind), stat=ok)
308	IF (ok/=0) THEN
309	write(,) "memory allocation failed for QREFir3d!"
310	call abort_physic(subname,'allocation failure for QREFir3d',1)
311	ENDIF
312	ENDIF
313	! Effective radius and variance of the aerosols
314	IF(.not.ALLOCATED(reffrad)) THEN
315	allocate(reffrad(ngrid,nlayer,naerkind), stat=ok)
316	IF (ok/=0) THEN
317	write(,) "memory allocation failed for reffrad!"
318	call abort_physic(subname,'allocation failure for reffrad',1)
319	ENDIF
320	ENDIF
321	IF(.not.ALLOCATED(nueffrad)) THEN
322	allocate(nueffrad(ngrid,nlayer,naerkind), stat=ok)
323	IF (ok/=0) THEN
324	write(,) "memory allocation failed for nueffrad!"
325	call abort_physic(subname,'allocation failure for nueffrad',1)
326	ENDIF
327	ENDIF
328
329	#ifndef MESOSCALE
330	if (is_master) call system('rm -f surf_vals_long.out')
331	#endif
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
436	if(varfixed .and. generic_condensation)then
437	write(,) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
438	qvap_deep=-1. ! default value
439	call getin_p("qvap_deep",qvap_deep)
440	write(,) " qvap_deep = ",qvap_deep
441
442	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
443	call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
444	endif
445
446	end if ! of if (firstcall)
447
448	!=======================================================================
449	! I.b Initialization on every call
450	!=======================================================================
451
452	qxvaer(:,:,:)=0.0
453	qsvaer(:,:,:)=0.0
454	gvaer(:,:,:) =0.0
455
456	qxiaer(:,:,:)=0.0
457	qsiaer(:,:,:)=0.0
458	giaer(:,:,:) =0.0
459
460	OLR_nu(:,:) = 0.
461	OSR_nu(:,:) = 0.
462	GSR_nu(:,:) = 0.
463
464	!--------------------------------------------------
465	! Effective radius and variance of the aerosols
466	!--------------------------------------------------
467	if (aerohaze) then
468	do iaer=1,naerkind
469	if ((iaer.eq.iaero_haze)) then
470	call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), &
471	nueffrad(1,1,iaer))
472	endif
473	end do !iaer=1,naerkind.
474	if (haze_radproffix) then
475	print*, 'haze_radproffix=T : fixed profile for haze rad'
476	else
477	print*,'reffrad haze:',reffrad(1,1,iaero_haze)
478	print*,'nueff haze',nueffrad(1,1,iaero_haze)
479	endif
480	endif
481
482
483	! How much light do we get ?
484	do nw=1,L_NSPECTV
485	stel(nw)=stellarf(nw)/(dist_star**2)
486	end do
487
488	if (aerohaze) then
489	! if (haze_radproffix) then
490	! call haze_reffrad_fix(ngrid,nlayer,zzlay,&
491	! reffrad,nueffrad)
492	! endif
493
494	! Get 3D aerosol optical properties.
495	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
496	QVISsQREF3d,omegaVIS3d,gVIS3d, &
497	QIRsQREF3d,omegaIR3d,gIR3d, &
498	QREFvis3d,QREFir3d)
499
500	! Get aerosol optical depths.
501	call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol, &
502	reffrad,nueffrad,QREFvis3d,QREFir3d, &
503	tau_col,cloudfrac,totcloudfrac,clearsky)
504	endif
505
506	!-----------------------------------------------------------------------
507	do ig=1,ngrid ! Starting Big Loop over every GCM column
508	!-----------------------------------------------------------------------
509
510
511	!=======================================================================
512	! II. Transformation of the GCM variables
513	!=======================================================================
514
515
516	!-----------------------------------------------------------------------
517	! Aerosol optical properties Qext, Qscat and g.
518	! The transformation in the vertical is the same as for temperature.
519	!-----------------------------------------------------------------------
520
521
522	if (aerohaze) then
523	do iaer=1,naerkind
524	! Shortwave.
525	do nw=1,L_NSPECTV
526
527	do l=1,nlayer
528
529	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
530	*QREFvis3d(ig,nlayer+1-l,iaer)
531
532	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
533	*QREFvis3d(ig,max(nlayer-l,1),iaer)
534
535	qxvaer(2*l,nw,iaer) = temp1
536	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
537
538	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
539	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
540
541	qsvaer(2*l,nw,iaer) = temp1
542	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
543
544	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
545	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
546
547	gvaer(2*l,nw,iaer) = temp1
548	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
549
550	end do ! nlayer
551
552	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
553	qxvaer(2*nlayer+1,nw,iaer)=0.
554
555	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
556	qsvaer(2*nlayer+1,nw,iaer)=0.
557
558	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
559	gvaer(2*nlayer+1,nw,iaer)=0.
560
561	end do ! L_NSPECTV
562
563	do nw=1,L_NSPECTI
564	! Longwave
565	do l=1,nlayer
566
567	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
568	*QREFir3d(ig,nlayer+1-l,iaer)
569
570	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
571	*QREFir3d(ig,max(nlayer-l,1),iaer)
572
573	qxiaer(2*l,nw,iaer) = temp1
574	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
575
576	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
577	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
578
579	qsiaer(2*l,nw,iaer) = temp1
580	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
581
582	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
583	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
584
585	giaer(2*l,nw,iaer) = temp1
586	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
587
588	end do ! nlayer
589
590	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
591	qxiaer(2*nlayer+1,nw,iaer)=0.
592
593	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
594	qsiaer(2*nlayer+1,nw,iaer)=0.
595
596	giaer(1,nw,iaer)=giaer(2,nw,iaer)
597	giaer(2*nlayer+1,nw,iaer)=0.
598
599	end do ! L_NSPECTI
600
601	end do ! naerkind
602
603	! Test / Correct for freaky s. s. albedo values.
604	do iaer=1,naerkind
605	do k=1,L_LEVELS
606
607	do nw=1,L_NSPECTV
608	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
609	message='Serious problems with qsvaer values'
610	call abort_physic(subname,message,1)
611	endif
612	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
613	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
614	endif
615	end do
616
617	do nw=1,L_NSPECTI
618	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
619	message='Serious problems with qsvaer values'
620	call abort_physic(subname,message,1)
621	endif
622	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
623	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
624	endif
625	end do
626
627	end do ! L_LEVELS
628	end do ! naerkind
629
630	endif ! aerohaze
631	!-----------------------------------------------------------------------
632	! Aerosol optical depths
633	!-----------------------------------------------------------------------
634	IF (aerohaze) THEN
635	do iaer=1,naerkind ! a bug was here
636	do k=0,nlayer-1
637
638	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
639	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
640	! As 'aerosol' is at reference (visible) wavelenght we scale it as
641	! it will be multplied by qxi/v in optci/v
642	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
643	tauaero(2k+2,iaer)=max(temppweight,0.d0)
644	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
645
646	end do
647	! boundary conditions
648	tauaero(1,iaer) = tauaero(2,iaer)
649	!tauaero(1,iaer) = 0.
650	!JL18 at time of testing, the two above conditions gave the same results bit for bit.
651
652	end do ! naerkind
653	ELSE
654	tauaero(:,:)=0
655	ENDIF ! aerohaze
656
657	! Albedo and Emissivity.
658	albi=1-emis(ig) ! Long Wave.
659	DO nw=1,L_NSPECTV ! Short Wave loop.
660	albv(nw)=albedo(ig,nw)
661	ENDDO
662
663	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
664
665
666	!-----------------------------------------------------------------------
667	! GCS (Generic Condensable Specie) Vapor
668	! If you have GCS tracers and they are : variable & radiatively active
669	!
670	! NC22
671	!-----------------------------------------------------------------------
672
673	if (generic_condensation) then
674
675	! For now, only one GCS tracer can be both variable and radiatively active
676	! If you set two GCS tracers, that are variable and radiatively active,
677	! the last one in tracer.def will be chosen as the one that will be vadiatively active
678
679	do iq=1,nq
680
681	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
682
683	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
684
685	if(varactive)then
686
687	i_var=igcm_generic_vap
688	do l=1,nlayer
689	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
690	qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)
691	!JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
692	!JL13index ! Average approximation as for temperature...
693	end do
694	qvar(1)=qvar(2)
695
696	elseif(varfixed .and. (qvap_deep .ge. 0))then
697
698	do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
699
700	call Psat_generic(pt(ig,l),pplay(ig,l),metallicity,psat,qsat)
701
702	if (qsat .lt. qvap_deep) then
703	pq_temp(l) = qsat ! fully saturated everywhere
704	else
705	pq_temp(l) = qvap_deep
706	end if
707
708	end do
709
710	do l=1,nlayer
711	qvar(2*l) = pq_temp(nlayer+1-l)
712	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
713	end do
714
715	qvar(1)=qvar(2)
716
717	else
718	do k=1,L_LEVELS
719	qvar(k) = 1.0D-7
720	end do
721	end if ! varactive/varfixed
722
723	endif
724
725	end do ! do iq=1,nq loop on tracers
726
727	end if ! if (generic_condensation)
728
729	!-----------------------------------------------------------------------
730	! No Water vapor and No GCS (Generic Condensable Specie) vapor
731	!-----------------------------------------------------------------------
732
733	if (.not. generic_condensation) then
734	do k=1,L_LEVELS
735	qvar(k) = 1.0D-7
736	end do
737	end if ! if (.not. generic_condensation)
738
739
740	if(.not.kastprof)then
741	! IMPORTANT: Now convert from kg/kg to mol/mol.
742	do k=1,L_LEVELS
743	if (generic_condensation) then
744	do iq=1,nq
745	call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
746	if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
747
748	epsi_generic=constants_epsi_generic(iq)
749
750	qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic))
751
752	endif
753	end do ! do iq=1,nq loop on tracers
754	endif
755	end do
756	end if
757
758	!-----------------------------------------------------------------------
759	! kcm mode only !
760	!-----------------------------------------------------------------------
761
762	DO l=1,nlayer
763	muvarrad(2*l) = muvar(ig,nlayer+2-l)
764	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
765	END DO
766
767	muvarrad(1) = muvarrad(2)
768	muvarrad(2*nlayer+1)=muvar(ig,1)
769
770	! Keep values inside limits for which we have radiative transfer coefficients !!!
771	if(L_REFVAR.gt.1)then ! (there was a bug here)
772	do k=1,L_LEVELS
773	if(qvar(k).lt.wrefvar(1))then
774	qvar(k)=wrefvar(1)+1.0e-8
775	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
776	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
777	endif
778	end do
779	endif
780
781	!-----------------------------------------------------------------------
782	! Pressure and temperature
783	!-----------------------------------------------------------------------
784
785	DO l=1,nlayer
786	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
787	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
788	tlevrad(2*l) = pt(ig,nlayer+1-l)
789	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
790	END DO
791
792	plevrad(1) = 0.
793	! 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.
794
795	tlevrad(1) = tlevrad(2)
796	tlevrad(2*nlayer+1)=tsurf(ig)
797
798	pmid(1) = pplay(ig,nlayer)/scalep
799	pmid(2) = pmid(1)
800
801	tmid(1) = tlevrad(2)
802	tmid(2) = tmid(1)
803
804	DO l=1,L_NLAYRAD-1
805	tmid(2l+1) = tlevrad(2l+1)
806	tmid(2l+2) = tlevrad(2l+1)
807	pmid(2l+1) = plevrad(2l+1)
808	pmid(2l+2) = plevrad(2l+1)
809	END DO
810	pmid(L_LEVELS) = plevrad(L_LEVELS)
811	tmid(L_LEVELS) = tlevrad(L_LEVELS)
812
813	!!Alternative interpolation:
814	! pmid(3) = pmid(1)
815	! pmid(4) = pmid(1)
816	! tmid(3) = tmid(1)
817	! tmid(4) = tmid(1)
818	! DO l=2,L_NLAYRAD-1
819	! tmid(2l+1) = tlevrad(2l)
820	! tmid(2l+2) = tlevrad(2l)
821	! pmid(2l+1) = plevrad(2l)
822	! pmid(2l+2) = plevrad(2l)
823	! END DO
824	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
825	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
826
827	! Test for out-of-bounds pressure.
828	if(plevrad(3).lt.pgasmin)then
829	print*,'Minimum pressure is outside the radiative'
830	print*,'transfer kmatrix bounds, exiting.'
831	message="Minimum pressure outside of kmatrix bounds"
832	call abort_physic(subname,message,1)
833	elseif(plevrad(L_LEVELS).gt.pgasmax)then
834	print*,'Maximum pressure is outside the radiative'
835	print*,'transfer kmatrix bounds, exiting.'
836	message="Minimum pressure outside of kmatrix bounds"
837	call abort_physic(subname,message,1)
838	endif
839
840	! Test for out-of-bounds temperature.
841	! -- JVO 20 : Also add a sanity test checking that tlevrad is
842	! within Planck function temperature boundaries,
843	! which would cause gfluxi/sfluxi to crash.
844	do k=1,L_LEVELS
845
846	if(tlevrad(k).lt.tgasmin)then
847	print*,'Minimum temperature is outside the radiative'
848	print*,'transfer kmatrix bounds'
849	print*,"k=",k," tlevrad(k)=",tlevrad(k)
850	print*,"tgasmin=",tgasmin
851	if (strictboundcorrk) then
852	message="Minimum temperature outside of kmatrix bounds"
853	call abort_physic(subname,message,1)
854	else
855	print,'**********************************************'
856	print*,'we allow model to continue with tlevrad<tgasmin'
857	print*,' ... we assume we know what you are doing ... '
858	print*,' ... but do not let this happen too often ... '
859	print,'**********************************************'
860	!tlevrad(k)=tgasmin ! Used in the source function !
861	endif
862	elseif(tlevrad(k).gt.tgasmax)then
863	print*,'Maximum temperature is outside the radiative'
864	print*,'transfer kmatrix bounds, exiting.'
865	print*,"k=",k," tlevrad(k)=",tlevrad(k)
866	print*,"tgasmax=",tgasmax
867	if (strictboundcorrk) then
868	message="Maximum temperature outside of kmatrix bounds"
869	call abort_physic(subname,message,1)
870	else
871	print,'**********************************************'
872	print*,'we allow model to continue with tlevrad>tgasmax'
873	print*,' ... we assume we know what you are doing ... '
874	print*,' ... but do not let this happen too often ... '
875	print,'**********************************************'
876	!tlevrad(k)=tgasmax ! Used in the source function !
877	endif
878	endif
879
880	if (tlevrad(k).lt.tplanckmin) then
881	print*,'Minimum temperature is outside the boundaries for'
882	print*,'Planck function integration set in callphys.def, aborting.'
883	print*,"k=",k," tlevrad(k)=",tlevrad(k)
884	print*,"tplanckmin=",tplanckmin
885	message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
886	call abort_physic(subname,message,1)
887	else if (tlevrad(k).gt.tplanckmax) then
888	print*,'Maximum temperature is outside the boundaries for'
889	print*,'Planck function integration set in callphys.def, aborting.'
890	print*,"k=",k," tlevrad(k)=",tlevrad(k)
891	print*,"tplanckmax=",tplanckmax
892	message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
893	call abort_physic(subname,message,1)
894	endif
895
896	enddo
897
898	do k=1,L_NLAYRAD+1
899	if(tmid(k).lt.tgasmin)then
900	print*,'Minimum temperature is outside the radiative'
901	print*,'transfer kmatrix bounds, exiting.'
902	print*,"k=",k," tmid(k)=",tmid(k)
903	print*,"tgasmin=",tgasmin
904	if (strictboundcorrk) then
905	message="Minimum temperature outside of kmatrix bounds"
906	call abort_physic(subname,message,1)
907	else
908	print,'**********************************************'
909	print*,'we allow model to continue but with tmid=tgasmin'
910	print*,' ... we assume we know what you are doing ... '
911	print*,' ... but do not let this happen too often ... '
912	print,'**********************************************'
913	tmid(k)=tgasmin
914	endif
915	elseif(tmid(k).gt.tgasmax)then
916	print*,'Maximum temperature is outside the radiative'
917	print*,'transfer kmatrix bounds, exiting.'
918	print*,"k=",k," tmid(k)=",tmid(k)
919	print*,"tgasmax=",tgasmax
920	if (strictboundcorrk) then
921	message="Maximum temperature outside of kmatrix bounds"
922	call abort_physic(subname,message,1)
923	else
924	print,'**********************************************'
925	print*,'we allow model to continue but with tmid=tgasmax'
926	print*,' ... we assume we know what you are doing ... '
927	print*,' ... but do not let this happen too often ... '
928	print,'**********************************************'
929	tmid(k)=tgasmax
930	endif
931	endif
932	enddo
933
934	!=======================================================================
935	! III. Calling the main radiative transfer subroutines
936	!=======================================================================
937
938	! ----------------------------------------------------------------
939	! Recombine reference corrk tables if needed - Added by JVO, 2020.
940	if (corrk_recombin) then
941	call call_recombin(ig,nlayer,pq(ig,:,:),pplay(ig,:),pt(ig,:),qvar(:),tmid(:),pmid(:))
942	endif
943	! ----------------------------------------------------------------
944
945	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
946	glat_ig=glat(ig)
947
948	!-----------------------------------------------------------------------
949	! Short Wave Part
950	!-----------------------------------------------------------------------
951
952	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
953	if((ngrid.eq.1).and.(global1d))then
954	do nw=1,L_NSPECTV
955	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
956	end do
957	else
958	do nw=1,L_NSPECTV
959	stel_fract(nw)= stel(nw) * fract(ig)
960	end do
961	endif
962
963	call optcv(dtauv,tauv,taucumv,plevrad, &
964	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
965	tmid,pmid,taugsurf,qvar,muvarrad)
966
967	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
968	acosz,stel_fract, &
969	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
970	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
971
972	else ! During the night, fluxes = 0.
973	nfluxtopv = 0.0d0
974	fluxtopvdn = 0.0d0
975	nfluxoutv_nu(:) = 0.0d0
976	nfluxgndv_nu(:) = 0.0d0
977	do l=1,L_NLAYRAD
978	fmnetv(l)=0.0d0
979	fluxupv(l)=0.0d0
980	fluxdnv(l)=0.0d0
981	end do
982	end if
983
984
985	! Equivalent Albedo Calculation (for OUTPUT). MT2015
986	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
987	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
988	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
989	DO nw=1,L_NSPECTV
990	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
991	ENDDO
992	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
993	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
994	else
995	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
996	endif
997	else
998	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
999	endif
1000
1001
1002	!-----------------------------------------------------------------------
1003	! Long Wave Part
1004	!-----------------------------------------------------------------------
1005
1006	call optci(plevrad,tlevrad,dtaui,taucumi, &
1007	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
1008	taugsurfi,qvar,muvarrad)
1009
1010	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
1011	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
1012	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
1013
1014	!-----------------------------------------------------------------------
1015	! Transformation of the correlated-k code outputs
1016	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
1017
1018	! Flux incident at the top of the atmosphere
1019	fluxtop_dn(ig)=fluxtopvdn
1020
1021	fluxtop_lw(ig) = real(nfluxtopi)
1022	fluxabs_sw(ig) = real(-nfluxtopv)
1023	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
1024	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
1025
1026	! Flux absorbed by the surface. By MT2015.
1027	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
1028
1029	if(fluxtop_dn(ig).lt.0.0)then
1030	print*,'Achtung! fluxtop_dn has lost the plot!'
1031	print*,'fluxtop_dn=',fluxtop_dn(ig)
1032	print*,'acosz=',acosz
1033	print*,'aerosol=',aerosol(ig,:,:)
1034	print*,'temp= ',pt(ig,:)
1035	print*,'pplay= ',pplay(ig,:)
1036	message="Achtung! fluxtop_dn has lost the plot!"
1037	call abort_physic(subname,message,1)
1038	endif
1039
1040	! Spectral output, for exoplanet observational comparison
1041	if(specOLR)then
1042	do nw=1,L_NSPECTI
1043	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
1044	end do
1045	do nw=1,L_NSPECTV
1046	GSR_nu(ig,nw)=nfluxgndv_nu(nw)/DWNV(nw)
1047	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
1048	end do
1049	endif
1050
1051	! Finally, the heating rates
1052
1053	DO l=2,L_NLAYRAD
1054	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
1055	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
1056	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
1057	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
1058	END DO
1059
1060	! These are values at top of atmosphere
1061	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
1062	glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
1063	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
1064	glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
1065
1066	! Optical thickness diagnostics (added by JVO)
1067	if (diagdtau) then
1068	do l=1,L_NLAYRAD
1069	do nw=1,L_NSPECTV
1070	int_dtauv(ig,l,nw) = 0.0d0
1071	DO k=1,L_NGAUSS
1072	! Output exp(-tau) because gweight ponderates exp and not tau itself
1073	int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
1074	ENDDO
1075	enddo
1076	do nw=1,L_NSPECTI
1077	int_dtaui(ig,l,nw) = 0.0d0
1078	DO k=1,L_NGAUSS
1079	! Output exp(-tau) because gweight ponderates exp and not tau itself
1080	int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
1081	ENDDO
1082	enddo
1083	enddo
1084	endif
1085
1086
1087	!-----------------------------------------------------------------------
1088	end do ! End of big loop over every GCM column.
1089	!-----------------------------------------------------------------------
1090
1091
1092
1093	!-----------------------------------------------------------------------
1094	! Additional diagnostics
1095	!-----------------------------------------------------------------------
1096
1097	! IR spectral output, for exoplanet observational comparison
1098	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
1099
1100	print*,'Saving scalar quantities in surf_vals.out...'
1101	print*,'psurf = ', pplev(1,1),' Pa'
1102	open(116,file='surf_vals.out')
1103	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
1104	real(-nfluxtopv),real(nfluxtopi)
1105	close(116)
1106
1107
1108	! USEFUL COMMENT - Do Not Remove.
1109	!
1110	! if(specOLR)then
1111	! open(117,file='OLRnu.out')
1112	! do nw=1,L_NSPECTI
1113	! write(117,*) OLR_nu(1,nw)
1114	! enddo
1115	! close(117)
1116	!
1117	! open(127,file='OSRnu.out')
1118	! do nw=1,L_NSPECTV
1119	! write(127,*) OSR_nu(1,nw)
1120	! enddo
1121	! close(127)
1122	! endif
1123
1124	! OLR vs altitude: do it as a .txt file.
1125	OLRz=.false.
1126	if(OLRz)then
1127	print*,'saving IR vertical flux for OLRz...'
1128	open(118,file='OLRz_plevs.out')
1129	open(119,file='OLRz.out')
1130	do l=1,L_NLAYRAD
1131	write(118,) plevrad(2l)
1132	do nw=1,L_NSPECTI
1133	write(119,*) fluxupi_nu(l,nw)
1134	enddo
1135	enddo
1136	close(118)
1137	close(119)
1138	endif
1139
1140	endif
1141
1142	! See physiq.F for explanations about CLFvarying. This is temporary.
1143	if (lastcall) then
1144	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
1145	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
1146	!$OMP BARRIER
1147	!$OMP MASTER
1148	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
1149	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
1150	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
1151	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
1152	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
1153	!$OMP END MASTER
1154	!$OMP BARRIER
1155	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
1156	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
1157	endif
1158
1159
1160	end subroutine callcorrk
1161
1162	END MODULE callcorrk_mod

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

Download in other formats:

Original Format