Context Navigation

callcorrk.F90 @ 2987

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

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

Download in other formats:

Original Format