Context Navigation

callcorrk.F90 @ 3187

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

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

Download in other formats:

Original Format