Context Navigation

callcorrk.F90 @ 2937

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

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

Download in other formats:

Original Format