Context Navigation

callcorrk.F90 @ 2870

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

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

Download in other formats:

Original Format