Context Navigation

callcorrk.F @ 3558

Last change on this file since 3558 was 3373, checked in by afalco, 12 months ago
Pluto.old PCM: Corrected some runtime issues with gfortran. AF
File size: 31.6 KB

Rev	Line
[3175]	1	subroutine callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
	2	& albedo,emis,mu0,pplev,pplay,pt,
	3	& tsurf,fract,dist_star,igout,aerosol,cpp3D,
	4	& dtlw,dtsw,fluxsurf_lw,
	5	& fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn,
	6	& reffrad,tau_col,ptime,pday,firstcall,lastcall,zzlay)
[3275]	7
[3175]	8	use radinc_h
	9	use radcommon_h
[3275]	10	use ioipsl_getincom
[3175]	11	use radii_mod
	12	use aerosol_mod
	13	use datafile_mod, only: datadir
	14
	15	implicit none
	16
	17	!==================================================================
	18	!
	19	! Purpose
	20	! -------
	21	! Solve the radiative transfer using the correlated-k method for
	22	! the gaseous absorption and the Toon et al. (1989) method for
	23	! scatttering due to aerosols.
	24	!
	25	! Authors
[3275]	26	! -------
[3175]	27	! Emmanuel 01/2001, Forget 09/2001
	28	! Robin Wordsworth (2009)
	29	! Modif Pluton Tanguy Bertrand 2017
	30	!==================================================================
	31
	32	#include "dimphys.h"
	33	#include "comcstfi.h"
	34	#include "callkeys.h"
	35	#include "tracer.h"
	36
	37	!-----------------------------------------------------------------------
	38	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
[3275]	39	! Layer #1 is the layer near the ground.
[3175]	40	! Layer #nlayermx is the layer at the top.
	41
	42	! INPUT
	43	INTEGER icount
	44	INTEGER ngrid,nlayer
	45	INTEGER igout
[3275]	46	REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau
[3175]	47	REAL albedo(ngrid) ! SW albedo
	48	REAL emis(ngrid) ! LW emissivity
	49	REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer
	50	REAL zzlay(ngrid,nlayermx) ! altitude at the middle of the layers
	51	REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries
	52
	53	REAL pt(ngrid,nlayermx) ! air temperature (K)
	54	REAL tsurf(ngrid) ! surface temperature (K)
	55	REAL dist_star,mu0(ngrid) ! distance star-planet (AU)
	56	REAL fract(ngrid) ! fraction of day
	57
	58	! Globally varying aerosol optical properties on GCM grid
	59	! Not needed everywhere so not in radcommon_h
	60	REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	61	REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	62	REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	63
	64	REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	65	REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	66	REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	67
	68	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
	69	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
	70
	71	! REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
	72	! REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
	73
	74	REAL reffrad(ngrid,nlayer,naerkind)
	75	REAL nueffrad(ngrid,nlayer,naerkind)
	76	REAL profhaze(ngrid,nlayer) ! TB17 fixed profile of haze mmr
	77
	78	! OUTPUT
	79	REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
	80	REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
	81	REAL fluxsurf_lw(ngridmx) ! incident LW flux to surf (W/m2)
	82	REAL fluxtop_lw(ngridmx) ! outgoing LW flux to space (W/m2)
	83	REAL fluxsurf_sw(ngridmx) ! incident SW flux to surf (W/m2)
	84	REAL fluxtop_sw(ngridmx) ! outgoing LW flux to space (W/m2)
	85	REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2)
	86
	87	!-----------------------------------------------------------------------
	88	! Declaration of the variables required by correlated-k subroutines
	89	! Numbered from top to bottom unlike in the GCM!
	90
	91	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	92	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	93
	94	! Optical values for the optci/cv subroutines
	95	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
	96	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	97	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	98	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	99	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	100	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	101	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	102	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	103	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	104	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	105
	106	REAL*8 tauaero(L_LEVELS+1,naerkind)
	107	REAL*8 nfluxtopv,nfluxtopi,nfluxtop
	108	real*8 NFLUXTOPV_nu(L_NSPECTV)
	109	real*8 NFLUXTOPI_nu(L_NSPECTI)
	110	real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
	111	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
[3275]	112	real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
	113	real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
[3175]	114	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	115	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
	116	REAL*8 albi,albv,acosz
	117
	118	INTEGER ig,l,k,nw,iaer,irad
	119
	120	real fluxtoplanet
	121	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	122	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	123
	124	real*8 qvar(L_LEVELS) ! mixing ratio of variable component
	125	REAL pq(ngridmx,nlayer,nq)
	126	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	127	integer nq
	128
	129	! Local aerosol optical properties for each column on RADIATIVE grid
	130	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	131	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	132	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	133	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	134	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	135	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	136
	137	save qxvaer, qsvaer, gvaer
	138	save qxiaer, qsiaer, giaer
[3275]	139	save QREFvis3d, QREFir3d
[3175]	140
	141	REAL tau_col(ngrid) ! diagnostic from aeropacity
	142
	143	! Misc.
	144	logical firstcall, lastcall, nantest
	145	real*8 tempv(L_NSPECTV)
	146	real*8 tempi(L_NSPECTI)
	147	real*8 temp,temp1,temp2,pweight
	148	character(len=10) :: tmp1
	149	character(len=10) :: tmp2
	150
	151	! for fixed vapour profiles
	152	real RH
	153	real*8 pq_temp(nlayer)
	154	real ptemp, Ttemp, qsat
	155
	156	! for OLR spec
	157	integer OLRcount
	158	save OLRcount
	159	integer OLRcount2
	160	save OLRcount2
	161	character(len=2) :: tempOLR
	162	character(len=30) :: filenomOLR
	163	real ptime, pday
	164
	165	REAL epsi_ch4
	166	SAVE epsi_ch4
	167
	168	logical diagrad_OLRz,diagrad_OLR,diagrad_surf,diagrad_rates
	169	real OLR_nu(ngrid,L_NSPECTI)
	170
	171	! Allow variations in cp with location
	172	REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
	173
	174	! NLTE factor for CH4
	175	real eps_nlte_sw23(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
	176	real eps_nlte_sw33(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
	177	real eps_nlte_lw(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw
	178	REAL dtlw_nu(nlayermx,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands
	179	REAL dtsw_nu(nlayermx,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands
	180
	181	REAL dpp ! intermediate
	182
	183	REAL dtlw_co(ngridmx, nlayermx) ! cooling rate (K/s) due to CO (diagnostic)
	184	REAL dtlw_hcn_c2h2(ngridmx, nlayermx) ! cooling rate (K/s) due to C2H2/HCN (diagnostic)
	185
	186	!!read altitudes and vmrch4
	187	integer Nfine,ifine
	188	parameter(Nfine=701)
	189	character(len=100) :: file_path
	190	real,save :: levdat(Nfine),vmrdat(Nfine)
	191	real :: vmrch4(ngridmx,nlayermx) ! vmr ch4 from vmrch4_proffix
	192
	193	!=======================================================================
	194	! Initialization on first call
	195
	196	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,qxvaer)
	197	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,qsvaer)
	198	CALL zerophys((L_LEVELS+1)L_NSPECTVnaerkind,gvaer)
	199
	200	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,qxiaer)
	201	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,qsiaer)
	202	CALL zerophys((L_LEVELS+1)L_NSPECTInaerkind,giaer)
	203
	204
	205	if(firstcall) then
	206
	207	print*, "callcorrk: Correlated-k data folder:",trim(datadir)
	208	call getin("corrkdir",corrkdir)
	209	print*, "corrkdir = ",corrkdir
	210	write( tmp1, '(i3)' ) L_NSPECTI
	211	write( tmp2, '(i3)' ) L_NSPECTV
	212	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	213	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	214
	215	print*,'starting sugas'
	216	call sugas_corrk ! set up gaseous absorption properties
	217	print*,'starting setspi'
	218	call setspi ! basic infrared properties
	219	print*,'starting setspv'
	220	call setspv ! basic visible properties
	221
[3275]	222	! Radiative Hazes
[3175]	223	if (aerohaze) then
	224
	225	print*,'aerohaze: starting suaer_corrk'
	226	call suaer_corrk ! set up aerosol optical properties
	227	print*,'ending suaer_corrk'
	228
	229	!--------------------------------------------------
	230	! Effective radius and variance of the aerosols
	231	!--------------------------------------------------
	232	do iaer=1,naerkind
[3275]	233	if ((iaer.eq.iaero_haze)) then
	234	call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer),
[3175]	235	& nueffrad(1,1,iaer))
	236	endif
	237	end do !iaer=1,naerkind.
	238	if (haze_radproffix) then
	239	print*, 'haze_radproffix=T : fixed profile for haze rad'
[3275]	240	else
[3175]	241	print*,'reffrad haze:',reffrad(1,1,iaero_haze)
	242	print*,'nueff haze',nueffrad(1,1,iaero_haze)
	243	endif
	244	endif ! radiative haze
	245
	246	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
	247
	248	epsi_ch4=mmol(igcm_ch4_gas)/mmol(igcm_n2)
	249
	250	! If fixed profile of CH4 gas
	251	IF (vmrch4_proffix) then
	252	file_path=trim(datadir)//'/gas_prop/vmr_ch4.txt'
	253	open(115,file=file_path,form='formatted')
	254	do ifine=1,Nfine
	255	read(115,*) levdat(ifine), vmrdat(ifine)
	256	enddo
	257	close(115)
	258	ENDIF
	259
	260	end if ! firstcall
	261
	262	!=======================================================================
	263
	264	! L_NSPECTV is the number of Visual(or Solar) spectral intervals
	265	! how much light we get
	266	fluxtoplanet=0
	267	DO nw=1,L_NSPECTV
	268	stel(nw)=stellarf(nw)/(dist_star**2) !flux
	269	fluxtoplanet=fluxtoplanet + stel(nw)
	270	END DO
	271
	272	!-----------------------------------------------------------------------
	273	! Get 3D aerosol optical properties.
[3275]	274	! ici on selectionne les proprietes opt correspondant a reffrad
[3175]	275	if (aerohaze) then
	276	!--------------------------------------------------
	277	! Effective radius and variance of the aerosols if profil non
	278	! uniform. Called only if aerohaze=true.
	279	!--------------------------------------------------
	280	if (haze_radproffix) then
	281	call haze_reffrad_fix(ngrid,nlayer,zzlay,
	282	& reffrad,nueffrad)
	283	endif
	284
[3275]	285	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
	286	& QVISsQREF3d,omegaVIS3d,gVIS3d,
	287	& QIRsQREF3d,omegaIR3d,gIR3d,
[3175]	288	& QREFvis3d,QREFir3d)
	289
	290	! Get aerosol optical depths.
[3275]	291	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,
	292	& reffrad,QREFvis3d,QREFir3d,
[3175]	293	& tau_col)
	294	endif
	295
	296	!-----------------------------------------------------------------------
	297	! Prepare CH4 mixing ratio for radiative transfer
	298	IF (methane) then
	299	vmrch4(:,:)=0.
[3275]	300
[3175]	301	if (ch4fix) then
	302	if (vmrch4_proffix) then
	303	!! Interpolate on the model vertical grid
	304	do ig=1,ngridmx
[3275]	305	! CALL interp_line(levdat,vmrdat,Nfine,
	306	! & zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
[3175]	307	enddo
	308	else
	309	vmrch4(:,:)=vmrch4fix
	310	endif
	311	else
	312	vmrch4(:,:)=pq(:,:,igcm_ch4_gas)100.
	313	& mmol(igcm_n2)/mmol(igcm_ch4_gas)
	314	endif
	315	ENDIF
	316
	317	! Prepare NON LTE correction in Pluto atmosphere
	318	IF (nlte) then
[3373]	319	CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
	320	& eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
[3175]	321	ENDIF
	322	c Net atmospheric radiative cooling rate from C2H2 (K.s-1):
	323	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	324	dtlw_hcn_c2h2=0.
	325	if (cooling) then
	326	CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,
	327	& pt,dtlw_hcn_c2h2)
	328	endif
	329
	330	!-----------------------------------------------------------------------
	331	!=======================================================================
	332	!-----------------------------------------------------------------------
	333	! starting big loop over every GCM column
	334	do ig=1,ngridmx
	335
	336	!=======================================================================
	337	! Transformation of the GCM variables
	338
	339	!-----------------------------------------------------------------------
	340	! Aerosol optical properties Qext, Qscat and g on each band
	341	! The transformation in the vertical is the same as for temperature
	342	if (aerohaze) then
	343	! shortwave
	344	do iaer=1,naerkind
[3275]	345	DO nw=1,L_NSPECTV
[3175]	346	do l=1,nlayermx
	347
[3275]	348	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)
[3175]	349	$ *QREFvis3d(ig,nlayermx+1-l,iaer)
	350
[3275]	351	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
[3175]	352	$ *QREFvis3d(ig,max(nlayermx-l,1),iaer)
	353	qxvaer(2*l,nw,iaer) = temp1
	354	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	355
	356	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
	357	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
	358
	359	qsvaer(2*l,nw,iaer) = temp1
	360	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	361
	362	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
	363	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
	364
	365	gvaer(2*l,nw,iaer) = temp1
	366	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	367
	368	end do
	369	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
	370	qxvaer(2*nlayermx+1,nw,iaer)=0.
	371
	372	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
	373	qsvaer(2*nlayermx+1,nw,iaer)=0.
	374
	375	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
	376	gvaer(2*nlayermx+1,nw,iaer)=0.
	377	end do
	378
	379	! longwave
[3275]	380	DO nw=1,L_NSPECTI
[3175]	381	do l=1,nlayermx
	382
[3275]	383	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)
[3175]	384	$ *QREFir3d(ig,nlayermx+1-l,iaer)
	385
[3275]	386	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
[3175]	387	$ *QREFir3d(ig,max(nlayermx-l,1),iaer)
	388
	389	qxiaer(2*l,nw,iaer) = temp1
	390	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	391
	392	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
	393	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
	394
	395	qsiaer(2*l,nw,iaer) = temp1
	396	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	397
	398	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
	399	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
	400
	401	giaer(2*l,nw,iaer) = temp1
	402	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	403
	404	end do
	405
	406	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
	407	qxiaer(2*nlayermx+1,nw,iaer)=0.
	408
	409	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
	410	qsiaer(2*nlayermx+1,nw,iaer)=0.
	411
	412	giaer(1,nw,iaer)=giaer(2,nw,iaer)
	413	giaer(2*nlayermx+1,nw,iaer)=0.
	414	end do
	415	end do ! naerkind
	416
	417	! Test / Correct for freaky s. s. albedo values.
	418	do iaer=1,naerkind
	419	do k=1,L_LEVELS
	420
	421	do nw=1,L_NSPECTV
	422	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
[3275]	423	print*,'Serious problems with qsvaer values'
[3175]	424	print*,'in callcorrk'
	425	call abort
	426	endif
	427	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
	428	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
	429	endif
	430	end do
	431
[3275]	432	do nw=1,L_NSPECTI
[3175]	433	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
	434	print*,'Serious problems with qsiaer values'
	435	print*,'in callcorrk'
	436	call abort
	437	endif
	438	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
	439	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
	440	endif
	441	end do
	442	end do ! levels
	443
	444	end do ! naerkind
	445
	446	endif ! aerohaze
	447
	448	!-----------------------------------------------------------------------
	449	! Aerosol optical depths
[3275]	450	IF (aerohaze) THEN
	451	do iaer=1,naerkind ! heritage generic
[3175]	452	do k=0,nlayer-1
	453	pweight=
	454	$ (pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/
	455	$ (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
[3275]	456	if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then
[3175]	457	temp=aerosol(ig,L_NLAYRAD-k,iaer)/
	458	$ QREFvis3d(ig,L_NLAYRAD-k,iaer)
	459	else
	460	print*, 'stop corrk',k,QREFvis3d(ig,L_NLAYRAD-k,iaer)
	461	stop
	462	end if
	463	tauaero(2k+2,iaer)=max(temppweight,0.d0)
	464	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0) ! tauaero en L_LEVELS soit deux fois plus que nlayer
	465	end do
	466
	467	! generic New boundary conditions
	468	tauaero(1,iaer) = tauaero(2,iaer)
	469	!tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
	470	!tauaero(1,iaer) = 0.
	471	!tauaero(L_LEVELS+1,iaer) = 0.
	472
	473	end do ! naerkind
	474	ELSE
	475	tauaero(:,:)=0
	476	ENDIF
	477	!-----------------------------------------------------------------------
	478
	479	! Albedo and emissivity
	480	albi=1-emis(ig) ! longwave
[3275]	481	albv=albedo(ig) ! shortwave
[3175]	482	acosz=mu0(ig) ! cosine of sun incident angle
	483
	484	!-----------------------------------------------------------------------
[3275]	485	! Methane vapour
[3175]	486
	487	c qvar = mixing ratio
	488	c L_LEVELS (51) different de GCM levels (25) . L_LEVELS = 2(llm-1)+3=2(ngridmx-1)+3
	489	c L_REFVAR The number of different mixing ratio values in
	490	c datagcm/composition.in for the k-coefficients.
	491	qvar(:)=0.
[3275]	492	IF (methane) then
[3175]	493
	494	do l=1,nlayer
	495	qvar(2l) = vmrch4(ig,nlayer+1-l)/100.
	496	& mmol(igcm_ch4_gas)/mmol(igcm_n2)
	497	qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig,
	498	& max(nlayer-l,1)))/2.)/100.*
	499	& mmol(igcm_ch4_gas)/mmol(igcm_n2)
	500	end do
	501	qvar(1)=qvar(2)
	502
	503
	504	! Keep values inside limits for which we have radiative transfer coefficients
	505	if(L_REFVAR.gt.1)then ! there was a bug here!
	506	do k=1,L_LEVELS
	507	if(qvar(k).lt.wrefvar(1))then
	508	qvar(k)=wrefvar(1)+1.0e-8
	509	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
	510	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
	511	endif
	512	end do
	513	endif
	514
	515	! IMPORTANT: Now convert from kg/kg to mol/mol
	516	do k=1,L_LEVELS
	517	qvar(k) = qvar(k)/(epsi_ch4+qvar(k)*(1.-epsi_ch4))
	518	end do
	519	ENDIF ! methane
	520
	521	!-----------------------------------------------------------------------
	522	! Pressure and temperature
	523
	524	! generic updated:
	525	DO l=1,nlayer
	526	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	527	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	528	tlevrad(2*l) = pt(ig,nlayer+1-l)
	529	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	530	END DO
	531
	532	plevrad(1) = 0.
	533	plevrad(2) = 0. !! Trick to have correct calculations of fluxes
	534	! in gflux(i/v).F, but the pmid levels are not impacted by this change.
	535
	536	tlevrad(1) = tlevrad(2)
	537	tlevrad(2*nlayer+1)=tsurf(ig)
	538
	539	pmid(1) = max(pgasmin,0.0001*plevrad(3))
	540	pmid(2) = pmid(1)
	541
	542	tmid(1) = tlevrad(2)
	543	tmid(2) = tmid(1)
	544
	545	! INI
	546	! DO l=1,nlayer
	547	! plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	548	! plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	549	! tlevrad(2*l) = pt(ig,nlayer+1-l)
	550	! tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,
	551	! $ max(nlayer-l,1)))/2
	552	! END DO
	553
	554	!! following lines changed in 03/2015 to solve upper atmosphere bug
	555	! plevrad(1) = 0.
[3275]	556	! plevrad(2) = max(pgasmin,0.0001*plevrad(3))
[3175]	557	!
	558	! tlevrad(1) = tlevrad(2)
	559	! tlevrad(2*nlayermx+1)=tsurf(ig)
	560	!
	561	! tmid(1) = tlevrad(2)
	562	! tmid(2) = tlevrad(2)
	563	!
	564	! pmid(1) = plevrad(2)
[3275]	565	! pmid(2) = plevrad(2)
[3175]	566
	567	DO l=1,L_NLAYRAD-1
	568	tmid(2l+1) = tlevrad(2l+1)
	569	tmid(2l+2) = tlevrad(2l+1)
	570	pmid(2l+1) = plevrad(2l+1)
	571	pmid(2l+2) = plevrad(2l+1)
	572	END DO
	573	! end of changes
	574	pmid(L_LEVELS) = plevrad(L_LEVELS)
	575	tmid(L_LEVELS) = tlevrad(L_LEVELS)
[3275]	576
[3175]	577	!TB
	578	if ((PMID(2).le.1.e-5).and.(ig.eq.1)) then
	579	if ((TMID(2).le.30.).and.(ig.eq.1)) then
	580	write(,) 'Caution! Pres/temp of upper levels lower than'
	581	write(,) 'ref pressure/temp: kcoef fixed for upper levels'
	582	!! cf tpindex.F
	583	endif
	584	endif
	585
	586	! test for out-of-bounds pressure
	587	if(plevrad(3).lt.pgasmin)then
	588	print*,'Minimum pressure is outside the radiative'
	589	print*,'transfer kmatrix bounds, exiting.'
	590	! call abort
	591	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	592	print*,'Maximum pressure is outside the radiative'
	593	print*,'transfer kmatrix bounds, exiting.'
	594	! call abort
	595	endif
	596
	597	! test for out-of-bounds temperature
	598	do k=1,L_LEVELS
	599	if(tlevrad(k).lt.tgasmin)then
	600	print*,'Minimum temperature is outside the radiative'
	601	print*,'transfer kmatrix bounds, exiting.'
	602	print*,tlevrad(k),k,tgasmin
	603	! call abort
	604	elseif(tlevrad(k).gt.tgasmax)then
	605	print*,'Maximum temperature is outside the radiative'
	606	print*,'transfer kmatrix bounds, exiting.'
	607	! call abort
	608	endif
	609	enddo
[3275]	610
[3175]	611	!=======================================================================
	612	! Calling the main radiative transfer subroutines
	613
	614	!-----------------------------------------------------------------------
	615	! Shortwave
[3275]	616
[3175]	617	IF(fract(ig) .GE. 1.0e-4) THEN ! only during daylight IPM?! flux UV...
	618
	619	fluxtoplanet=0.
	620	DO nw=1,L_NSPECTV
	621	stel_fract(nw)= stel(nw) * fract(ig)
	622	fluxtoplanet=fluxtoplanet + stel_fract(nw)
	623	END DO
	624
[3275]	625	!print*, 'starting optcv'
[3175]	626	call optcv(dtauv,tauv,taucumv,plevrad,
	627	$ qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,
	628	$ tmid,pmid,taugsurf,qvar)
	629
	630	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,
	631	$ acosz,stel_fract,gweight,nfluxtopv,nfluxtopv_nu,
	632	$ fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf)
	633
	634	ELSE ! during the night, fluxes = 0
	635	nfluxtopv=0.0
	636	DO l=1,L_NLAYRAD
	637	fmnetv(l)=0.0
	638	fluxupv(l)=0.0
	639	fluxdnv(l)=0.0
	640	END DO
	641	END IF
	642
	643	!-----------------------------------------------------------------------
	644	! Longwave
	645
	646	call optci(plevrad,tlevrad,dtaui,taucumi,
	647	$ qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,
	648	$ taugsurfi,qvar)
	649	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,
	650	$ wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu,
	651	$ fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	652
	653
	654	!-----------------------------------------------------------------------
	655	! Transformation of the correlated-k code outputs
	656	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	657
	658	fluxtop_lw(ig) = fluxupi(1)
	659	fluxsurf_lw(ig) = fluxdni(L_NLAYRAD)
	660	fluxtop_sw(ig) = fluxupv(1)
	661	fluxsurf_sw(ig) = fluxdnv(L_NLAYRAD)
	662
	663	! Flux incident at the top of the atmosphere
	664	fluxtop_dn(ig)=fluxdnv(1)
	665
	666	! IR spectral output from top of the atmosphere
	667	if(specOLR)then
[3275]	668	do nw=1,L_NSPECTI
[3175]	669	OLR_nu(ig,nw)=nfluxtopi_nu(nw)
	670	end do
	671	endif
	672
	673	! **********************************************************
	674	! Finally, the heating rates
[3275]	675	! g/cp*DF/DP
[3175]	676	! **********************************************************
	677
	678	DO l=2,L_NLAYRAD
	679	dpp = g/(cppscalep(plevrad(2l+1)-plevrad(2l-1)))
	680
	681	! DTSW :
	682	!dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))*dpp !averaged dtlw on each wavelength
	683	do nw=1,L_NSPECTV
[3275]	684	dtsw_nu(L_NLAYRAD+1-l,nw)=
[3175]	685	& (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
	686	end do
	687
	688	! DTLW : detailed with wavelength so that we can apply NLTE
	689	!dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))*dpp !averaged dtlw on each wavelength
	690	do nw=1,L_NSPECTI
[3275]	691	dtlw_nu(L_NLAYRAD+1-l,nw)=
	692	& (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
[3175]	693	end do
[3275]	694	END DO
	695
[3175]	696	! values at top of atmosphere
	697	dpp = g/(cppscalep(plevrad(3)-plevrad(1)))
	698
[3275]	699	! SW
[3175]	700	!dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)*dpp
	701	do nw=1,L_NSPECTV
	702	dtsw_nu(L_NLAYRAD,nw)=
	703	& (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp
	704	end do
	705
[3275]	706	! LW
	707	c dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp
[3175]	708	do nw=1,L_NSPECTI
	709	dtlw_nu(L_NLAYRAD,nw)=
	710	& (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp
	711	end do
	712
	713	! **********************************************************
	714	! NON NLTE correction in Pluto atmosphere
	715	! And conversion of LW spectral heating rates to total rates
	716	! **********************************************************
	717
	718	if (.not.nlte) then
[3275]	719	eps_nlte_sw23(ig,:) =1. ! IF no NLTE
	720	eps_nlte_sw33(ig,:) =1. ! IF no NLTE
	721	eps_nlte_lw(ig,:) =1. ! IF no NLTE
[3175]	722	endif
[3275]	723
[3175]	724	do l=1,nlayer
	725
	726	!LW
[3275]	727	dtlw(ig,l) =0.
[3175]	728	! dtlw_co(ig,l) =0. ! only for diagnostic
	729	do nw=1,L_NSPECTI
	730	! wewei : wavelength in micrometer
[3275]	731	if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
[3175]	732	dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
[3275]	733	else
[3175]	734	!dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
	735	dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
	736	! dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
	737	end if
[3275]	738	dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
[3175]	739	end do
	740	! adding c2h2 if cooling active
	741	dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l)
	742
	743	!SW
	744	dtsw(ig,l) =0.
[3275]	745
[3175]	746	if (strobel) then
	747
	748	do nw=1,L_NSPECTV
[3275]	749	if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
[3175]	750	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
[3275]	751	elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
[3175]	752	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
	753	else
	754	dtsw_nu(l,nw)=dtsw_nu(l,nw)
	755	end if
[3275]	756	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
[3175]	757	end do
	758
	759	else ! total heating rates multiplied by nlte coef
	760
	761	do nw=1,L_NSPECTV
	762	dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
	763	dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
	764	enddo
	765
[3275]	766	endif
[3175]	767
	768
	769	end do
	770	! **********************************************************
	771
	772	! Diagnotics for last call for each grid point
[3275]	773	!if (lastcall) then
[3175]	774
	775	!print*,'albedi vis=',albv
	776	!print*,'albedo ir=',albi
	777	!print*,'fluxup ir (:)=',fluxupi
	778	!print*,'flux ir net (:)=',fluxdni-fluxupi
	779	!print*,'cumulative flux net ir (:)=',fmneti
	780	!print*,'cumulative flux net vis (:)=',fmnetv
	781	!print*,'fluxdn vis (:)=',fluxdnv
	782	!print*,'fluxtop vis=',fluxtop_sw
	783	!print*,'fluxsurf vis=',fluxsurf_sw
	784	!print*,'fluxtop ir=',fluxtop_lw
	785	!print*,'fluxsurf ir=',fluxsurf_lw
	786
	787	c write(,) 'pressure, eps_nlte_sw, eps_nlte_lw'
	788	c do l=1,nlayer
	789	c write(,)pplay(1,l),eps_nlte_sw(1,l),eps_nlte_lw(1,l)
	790	c end do
	791
	792	!endif
	793
	794	! ---------------------------------------------------------------
	795	end do ! end of big loop over every GCM column (ig = 1:ngrid)
	796
	797	!-----------------------------------------------------------------------
	798	! Additional diagnostics
	799
	800	! IR spectral output, for exoplanet observational comparison
	801	if(specOLR)then
	802	if(ngrid.ne.1)then
	803	call writediagspec(ngrid,"OLR3D",
	804	& "OLR(lon,lat,band)","W m^-2",3,OLR_nu)
	805	endif
	806	endif
	807
[3275]	808	if(lastcall)then
[3175]	809
	810	! 1D Output
	811	if(ngrid.eq.1)then
	812
	813	! surface diagnotics
	814	diagrad_surf=.true.
	815	if(diagrad_surf)then
	816	open(116,file='surf_vals.out')
	817	write(116,*) tsurf(1),pplev(1,1),
[3275]	818	& fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1)
[3175]	819	do nw=1,L_NSPECTV
	820	write(116,*) wavev(nw),fmnetv_nu(L_NLAYRAD,nw)
	821	enddo
	822	do nw=1,L_NSPECTI
	823	write(116,*) wavei(nw),fmneti_nu(L_NLAYRAD,nw)
	824	enddo
	825	close(116)
	826	endif
	827
	828	! OLR by band
	829	diagrad_OLR=.true.
	830	if(diagrad_OLR)then
	831	open(117,file='OLRnu.out')
[3275]	832	write(117,*) 'IR wavel - band width - OLR'
[3175]	833	do nw=1,L_NSPECTI
	834	write(117,*) wavei(nw),
	835	& abs(1.e4/bwnv(nw)-1.e4/bwnv(nw+1)),OLR_nu(1,nw)
[3275]	836	enddo
[3175]	837	close(117)
	838	endif
	839
	840	! OLR vs altitude: in a .txt file
	841	diagrad_OLRz=.true.
	842	if(diagrad_OLRz)then
	843	open(118,file='OLRz_plevs.out')
	844	open(119,file='OLRz.out')
	845	do l=1,L_NLAYRAD
	846	write(118,) plevrad(2l)
	847	do nw=1,L_NSPECTI
[3275]	848	write(119,*) fluxupi_nu(l,nw)
[3175]	849	enddo
[3275]	850	enddo
[3175]	851	close(118)
	852	close(119)
	853	endif
	854
	855	! Heating rates vs altitude in a .txt file
	856	diagrad_rates=.true.
	857	if(diagrad_rates)then
	858	open(120,file='heating_rates.out')
	859	write(120,*) "Pressure - Alt - HR tot - Rates (wavel SW)"
	860	do l=1,nlayermx
	861	write(120,*) pplay(1,l),zzlay(1,l),dtsw(1,l),dtsw_nu(l,:)
	862	enddo
	863	close(120)
	864
	865	open(121,file='cooling_rates.out')
	866	write(121,*) "Pressure - Alt - CR tot - Rates (wavel LW)"
	867	do l=1,nlayermx
	868	write(121,*) pplay(1,l),zzlay(1,l),dtlw(1,l),dtlw_nu(l,:)
	869	enddo
	870	close(121)
	871
	872	open(122,file='bands.out')
	873	write(122,*) "wavel - bands boundaries (microns)"
	874	do nw=1,L_NSPECTV
	875	write(122,*) wavev(nw),1.e4/bwnv(nw+1),1.e4/bwnv(nw)
	876	enddo
	877	do nw=1,L_NSPECTI
	878	write(122,*) wavei(nw),1.e4/bwni(nw+1),1.e4/bwni(nw)
	879	enddo
	880	close(122)
	881
	882	open(123,file='c2h2_rates.out')
	883	write(123,*) "Pressure - Alt - CR c2h2"
	884	do l=1,nlayermx
	885	write(123,*) pplay(1,l),zzlay(1,l),dtlw_hcn_c2h2(1,l)
	886	enddo
	887	close(123)
	888
	889	endif
	890
	891	endif ! ngrid.eq.1
	892
	893	endif ! lastcall
	894
	895	end

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

Download in other formats:

Original Format