Context Navigation

callcorrk.F90 @ 3105

Last change on this file since 3105 was 3044, checked in by jleconte, 21 months ago
perfect_vap_profile option
Property svn:executable set to ``*
File size: 54.9 KB

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

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

Download in other formats:

Original Format