Context Navigation

callcorrk.F90 @ 2719

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

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

Download in other formats:

Original Format