Context Navigation

callcorrk.F90 @ 3469

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

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

Download in other formats:

Original Format