Context Navigation

callcorrk.F90 @ 3469

Last change on this file since 3469 was 3259, checked in by gmilcareck, 9 months ago
Fix a bug with global1d. When using global1d, the latitude used when global1d=false was still taken into account at global1d=true in the conditions for the presence of solar flux.
Property svn:executable set to ``*
File size: 55.6 KB

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

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

Download in other formats:

Original Format