Context Navigation

callcorrk.F90 @ 3614

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

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

Download in other formats:

Original Format