Context Navigation

callcorrk.F90 @ 3444

Last change on this file since 3444 was 3444, checked in by afalco, 2 months ago
Pluto PCM: reffrad not written if not initialized. AF
Property svn:executable set to ``*
File size: 52.8 KB

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

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

Download in other formats:

Original Format