Context Navigation

surfprop.F90 @ 3607

Last change on this file since 3607 was 3539, checked in by tbertrand, 2 months ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
File size: 21.2 KB

Rev	Line
[3483]	1	subroutine surfprop(ngrid,nq,pfra,qsurf,tsurface, &
	2	capcal,adjust,dist,fluold,ptimestep,zls, &
	3	albedo,emis,therm_inertia)
[3275]	4
	5	! use comgeomfi_h, only:
[3329]	6	use radinc_h, only: L_NSPECTV, L_NSPECTI
[3275]	7	use comcstfi_mod, only: pi
	8	use comsoil_h, only: nsoilmx,mlayer,volcapa,inertiedat
	9	use geometry_mod, only: latitude, longitude
	10	! use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
	11	use callkeys_mod, only: methane,carbox,mode_n2,mode_ch4,mode_tholins, &
	12	feedback_met,fdch4_depalb,fdch4_finalb, &
	13	thres_ch4ice,thres_coice,thres_n2ice, &
	14	changeti,changetid,deltab,specalb, &
	15	fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_latn, &
	16	fdch4_maxalb,fdch4_maxice,fdch4_ampl, &
	17	tholateq,tholatn,tholats,tholone,tholonw, &
	18	metlateq,metls1,metls2
	19	use surfdat_h, only: alb_n2a,alb_n2b,alb_ch4,alb_ch4_eq,alb_ch4_s, &
	20	alb_co,alb_tho,alb_tho_eq,alb_tho_spe,albedodat, &
	21	emis_n2b,emis_n2a,emis_ch4,emis_co, &
	22	emis_tho,emis_tho_eq,emis_tho_spe, &
	23	itch4,itch4d,ith2o,ith2od,itn2,itn2d
	24	USE tracer_h, only: igcm_ch4_ice,igcm_co_ice,igcm_n2
	25	use time_phylmdz_mod, only: daysec
	26
	27	implicit none
	28
	29	!==================================================================
	30	! Purpose
	31	! -------
	32	! set up the properties of pluto's surface
	33	!
	34	! Inputs
	35	! ------
	36	! ngrid Number of vertical columns
	37	! nlayer Number of layers
	38	! qsurf(ngrid,iq) Amount of ice on surface kg/m2
	39	! tsurface(ngrid) surface temperature
	40	!
	41	! Outputs
	42	! -------
	43	! albedo(ngrid)
	44	! emis(ngrid)
[3483]	45	! therm_inertia(ngrid,nsoilmx)
[3275]	46	!
	47	! Both
	48	! ----
	49	!
	50	! Authors
	51	! -------
	52	! Tanguy Bertrand
	53	!
	54	!==================================================================
	55
	56	!-----------------------------------------------------------------------
	57	! Arguments
	58
[3483]	59
[3275]	60	INTEGER ngrid, nq
[3483]	61	REAL,INTENT(IN) :: pfra(ngrid)
[3275]	62	REAL,INTENT(IN) :: qsurf(ngrid,nq)
[3483]	63	REAL,INTENT(IN) :: tsurface(ngrid)
[3275]	64	REAL,INTENT(IN) :: fluold(ngrid,nq)
	65	REAL,INTENT(IN) :: ptimestep
	66	REAL,INTENT(IN) :: zls
	67	REAL,INTENT(IN) :: capcal(ngrid)
	68	REAL,INTENT(IN) :: adjust
	69	REAL,INTENT(IN) :: dist
[3329]	70	REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
[3275]	71	REAL,INTENT(OUT) :: emis(ngrid)
[3539]	72	REAL,INTENT(OUT) :: therm_inertia(ngrid,nsoilmx)
[3275]	73	!-----------------------------------------------------------------------
	74	! Local variables
	75	REAL stephan
	76	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	77	SAVE stephan
	78	REAL tsa,tsb,coef,aa,bb
	79	REAL emin,emax,emax2,tid
	80	REAL n2cover,n2coversun,gamm,frac,pls,facls
	81	INTEGER ig,isoil,ice_case
	82	LOGICAL firstcall
	83	SAVE firstcall
	84	DATA firstcall/.true./
	85	!-----------------------------------------------------------------------
	86	! Special aging for methane (not implemented)
	87	! REAL albstep ! Time constant for albedo change (age)
	88	! DATA albstep/1.e7/
	89	! REAL albmet(ngrid)
	90	! SAVE albmet
	91	! REAL albmetmin ! Min alb for methane ice
	92	! DATA albmetmin/0.5/
	93	! REAL albmetmax ! Max alb for methane ice
	94	! DATA albmetmax/0.8/
	95
	96	!-----------------------------------------------------------------------
	97	! 1) ALBEDOS and EMISSIVITY
	98	! A. N2
	99	! CASE (0) ! fixed albedo
	100	! CASE (1) ! Albedo decreases with thickness
	101	! CASE (2) ! Special Sputnik differences of albedo
	102	! CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimationi rates
	103	! CASE (4) ! Albedo Difference in N/S (e.g. used for Triton)
	104	! CASE (5) ! Special Sputnik differences of albedo in small (1 pixel) patches (e.g. simulating dark patches / plumes)
	105	! --> EMISSIVITY N2: based on the alpha/beta transition
	106	! B. CO
	107	! C. CH4
	108	! CASE (0) ! 2 albedos, one for the tropics, one for the poles
	109	! CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south)
	110	! CASE (2) ! 2 albedos + albedo feedback
	111	! SELECT CASE (feedback_met)
	112	! CASE (0) ! Default (linear from alb_ch4_eq)
	113	! CASE (1) ! Hyperbolic tangent old
	114	! CASE (2) ! hyperbolic tangent old
	115	! CASE (3) ! hyperbolic tangent equation with parameters
	116	! CASE (3) ! Eq, poles N, pole S + depending on Ls
	117	! D. Tholins
	118	! CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
	119	! CASE (1) ! Special mode one region with a different albedo
	120	! E. Tholins read from file
	121	! specalb
	122
	123	! 2) Changes of Thermal inertia with time
	124	! if (changeti) : change of seasonal TI
	125	! if (changetid) : change of diurnal TI
	126
	127	! 3) OTHER TESTS
	128	!-----------------------------------------------------------------------
	129	! 1) ALBEDOS and EMISSIVITY
	130	!-----------------------------------------------------------------------
	131	! Some input parameters:
	132	pls=zls*180./pi
	133	! for equation of feedbackalbedo
	134	if (methane.and.mode_ch4.eq.2.and.feedback_met.eq.3) then
	135	aa=fdch4_finalb-fdch4_depalb
	136	bb=fdch4_finalb
	137	endif
[3377]	138	emis(:) = 1 !default to 1
[3275]	139
	140	! Loop on all points
	141	DO ig=1,ngrid
	142
	143	! Looking for dominant ice:
	144	ice_case=0 ! Tholins
	145	if (methane) then
	146	if (qsurf(ig,igcm_ch4_ice).ge.thres_ch4ice) then
	147	ice_case=1
	148	endif
	149	endif
	150	if (carbox) then
	151	if (qsurf(ig,igcm_co_ice).ge.thres_coice) then
	152	ice_case=2
	153	endif
	154	endif
	155	if (qsurf(ig,igcm_n2).ge.thres_n2ice) then
	156	ice_case=3
	157	endif
	158
	159	!---------------------------------
	160	! 1.A. N2
	161	!---------------------------------
	162
	163	if (ice_case.eq.3) then
	164
	165	! EMISSIVITY N2
	166	! change emis if phase alpha different de phase beta
	167	frac=0.5(1.+tanh(6.(35.6-tsurface(ig))/0.5))
	168	if (frac.gt.1.) frac=1.
	169	emis(ig)=fracemis_n2a+(1.-frac)emis_n2b
	170
	171	! ALBEDO N2 and emissivity changes
	172	SELECT CASE (mode_n2)
	173
	174	CASE (0) ! fixed albedo
[3329]	175	albedo(ig,:)=min(max(alb_n2b+adjust,0.),1.)
[3275]	176
	177	CASE (1) ! Albedo decreases with thickness
[3329]	178	albedo(ig,:)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b
	179	albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b) ! should not be higher than alb_n2b and not lower than alb_n2b-deltab
[3275]	180	CASE (2) ! Special Sputnik differences of albedo
	181	if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)180./pi.le.-161.or.longitude(ig)180./pi.ge.146)) then
	182	if ( (latitude(ig)*180./pi.ge.25.).or. &
	183	((longitude(ig)*180./pi.gt.140.).and. &
	184	(longitude(ig)*180./pi.lt.165.)) ) then
[3329]	185	albedo(ig,:)=alb_n2b-deltab
[3275]	186	else
[3329]	187	albedo(ig,:)=alb_n2b+deltab
[3275]	188	endif
	189	else
[3329]	190	albedo(ig,:)=alb_n2b
[3275]	191	endif
	192
	193	CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates
[3329]	194	albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
	195	albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
[3275]	196	!! Triton:
[3329]	197	!albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
	198	!albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
[3275]	199
	200	CASE (4) ! Albedo Difference in N/S
	201	if (latitude(ig)*180./pi.ge.0.) then
[3329]	202	albedo(ig,:)=min(max(alb_n2b-deltab+adjust,0.),1.)
[3275]	203	else
[3329]	204	albedo(ig,:)=min(max(alb_n2b+deltab+adjust,0.),1.)
[3275]	205	endif
	206
	207	CASE (5) ! Special Sputnik differences of albedo in patches
	208	!! Patches Nogcm
	209	!if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)180./pi.le.180.).and.(longitude(ig)180./pi.ge.174.) ) then
	210	! if (((latitude(ig)180./pi.le.46.).and.(latitude(ig)180./pi.ge.42.)) &
	211	! .or. ((latitude(ig)180./pi.le.36.).and.(latitude(ig)180./pi.ge.32.)) &
	212	! .or. ((latitude(ig)180./pi.le.26.).and.(latitude(ig)180./pi.ge.22.)) &
	213	! .or. ((latitude(ig)180./pi.le.16.).and.(latitude(ig)180./pi.ge.12.)) &
	214	! .or. ((latitude(ig)180./pi.le.6.).and.(latitude(ig)180./pi.ge.2.)) &
	215	! ) then
[3329]	216	! albedo(ig,:)=alb_n2b-deltab
[3275]	217
	218	!! Patches GCM
	219	if (qsurf(ig,igcm_n2).ge.1.e4) then
	220	if (((latitude(ig)180./pi.lt.33.).and.(latitude(ig)180./pi.gt.32.).and. &
	221	(longitude(ig)180./pi.lt.165.5).and.(longitude(ig)180./pi.gt.164.5)) &
	222	.or. ((latitude(ig)180./pi.lt.30.5).and.(latitude(ig)180./pi.gt.29.5).and. &
	223	(longitude(ig)180./pi.lt.169.).and.(longitude(ig)180./pi.gt.168.))) then
[3329]	224	albedo(ig,:)=alb_n2b-deltab
[3275]	225	else if (((latitude(ig)180./pi.lt.30.5).and.(latitude(ig)180./pi.gt.29.).and. &
	226	(longitude(ig)180./pi.lt.165.5).and.(longitude(ig)180./pi.gt.164.5)) &
	227	.or. ((latitude(ig)180./pi.lt.33.).and.(latitude(ig)180./pi.gt.32.).and. &
	228	(longitude(ig)180./pi.lt.169.).and.(longitude(ig)180./pi.gt.168.))) then
[3329]	229	albedo(ig,:)=alb_n2b+deltab
[3275]	230	else
[3329]	231	albedo(ig,:)=alb_n2b
[3275]	232	endif
	233	else
[3329]	234	albedo(ig,:)=alb_n2b
[3275]	235	endif
	236
	237	CASE (7) ! Albedo from albedodat and adjusted emissivity
[3329]	238	albedo(ig,:)=albedodat(ig)
[3275]	239	! adjust emissivity if convergeps = true
	240	emis(ig)=min(max(emis(ig)+adjust,0.),1.)
	241
	242	CASE DEFAULT
	243	write(,) 'STOP in surfprop.F90: mod_n2 not found'
	244	stop
	245	END SELECT
	246
	247
	248	!---------------------------------
	249	! 1.B. CO
	250	!---------------------------------
	251
	252	else if (ice_case.eq.2) then
[3329]	253	albedo(ig,:)=alb_co
[3275]	254	emis(ig)=emis_co
	255
	256	!---------------------------------
	257	! 1.C. CH4
	258	!---------------------------------
	259
	260	else if (ice_case.eq.1) then
	261
	262	SELECT CASE (mode_ch4)
	263
	264	CASE (0) ! 2 albedos, one for the tropics, one for the poles
	265	emis(ig)=emis_ch4
	266	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
[3329]	267	albedo(ig,:)=alb_ch4_eq
[3275]	268	else
[3329]	269	albedo(ig,:)=alb_ch4
[3275]	270	endif
	271
	272	CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south)
	273	emis(ig)=emis_ch4
	274	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
[3329]	275	albedo(ig,:)=alb_ch4_eq
[3275]	276	else if (latitude(ig)*180./pi.le.-metlateq) then
[3329]	277	albedo(ig,:)=alb_ch4_s
[3275]	278	else
[3329]	279	albedo(ig,:)=alb_ch4
[3275]	280	endif
	281
	282	CASE (2) ! 2 albedos + albedo feedback
	283
	284	emis(ig)=emis_ch4
[3329]	285	albedo(ig,:)=alb_ch4
[3275]	286
	287	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
[3329]	288	albedo(ig,:)=alb_ch4_eq
[3275]	289	endif
	290
	291	!! Albedo feedback
	292	if (latitude(ig)180./pi.le.fdch4_latn.and.latitude(ig)180./pi.gt.fdch4_lats) then
	293	if (longitude(ig)180./pi.le.fdch4_lone.and.longitude(ig)180./pi.gt.fdch4_lonw) then
	294	if (qsurf(ig,igcm_ch4_ice).lt.fdch4_maxice) then ! fd would not apply on BTD
	295	SELECT CASE (feedback_met)
	296	CASE (0) ! Default (linear from alb_ch4_eq)
[3329]	297	albedo(ig,:)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq
[3275]	298	!emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4
[3329]	299	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
[3275]	300	!if (emis(ig).gt.1.) emis(ig)=1.
	301	CASE (1) ! Hyperbolic tangent old
[3329]	302	albedo(ig,:)=0.60.5(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !
[3275]	303	!emis(ig)=0.20.5(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 !
[3329]	304	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
[3275]	305	!if (emis(ig).gt.1.) emis(ig)=1.
	306	CASE (2) ! hyperbolic tangent old
[3329]	307	albedo(ig,:)=0.50.6(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !
[3275]	308	!emis(ig)=0.20.5(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 !
[3329]	309	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
[3275]	310	!if (emis(ig).gt.1.) emis(ig)=1.
	311	CASE (3) ! hyperbolic tangent equation with parameters
[3329]	312	albedo(ig,:)=aa(-1+tanh(fdch4_amplqsurf(ig,igcm_ch4_ice)))+bb
	313	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
[3275]	314	CASE DEFAULT
	315	write(,) 'STOP surfprop.F90: fd wrong'
	316	stop
	317	END SELECT
	318	endif
	319	endif
	320	endif
	321
	322	CASE (3) ! Eq, poles N, pole S + depending on Ls
	323	emis(ig)=emis_ch4
	324	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
[3329]	325	albedo(ig,:)=alb_ch4_eq ! Pure methane ice
[3275]	326	else if (latitude(ig)*180./pi.le.-metlateq) then
[3329]	327	albedo(ig,:)=alb_ch4_s ! Pure methane ice
[3275]	328	if (pls.le.metls2.and.pls.gt.metls1) then
[3329]	329	albedo(ig,:)=alb_ch4_s+deltab ! Also used for N2, careful
[3275]	330	endif
	331	else
[3329]	332	albedo(ig,:)=alb_ch4 ! Pure methane ice
[3275]	333	endif
	334
	335	CASE (4) ! Albedo from albedodat
	336	emis(ig)=emis_ch4
[3329]	337	albedo(ig,:)=albedodat(ig)
[3275]	338
	339	CASE DEFAULT
	340	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
	341	stop
	342	END SELECT
	343
	344	!---------------------------------
	345	! 1.D. THOLINS
	346	!---------------------------------
	347
	348	else
	349
	350	SELECT CASE (mode_tholins)
	351
	352	CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
	353
	354	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
[3329]	355	albedo(ig,:)=alb_tho_eq
[3275]	356	emis(ig)=emis_tho_eq
	357	else
[3329]	358	albedo(ig,:)=alb_tho ! default = 0.1
[3275]	359	emis(ig)=emis_tho
	360	endif
	361
	362	CASE (1) ! Special mode one region with a different albedo
	363
	364	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
[3329]	365	albedo(ig,:)=alb_tho_eq
[3275]	366	emis(ig)=emis_tho_eq
	367	else
[3329]	368	albedo(ig,:)=alb_tho ! default = 0.1
[3275]	369	emis(ig)=emis_tho
	370	endif
	371
	372	if (latitude(ig)180./pi.le.tholatn.and.latitude(ig)180./pi.ge.tholats &
	373	.and.longitude(ig)180./pi.ge.tholone.and.longitude(ig)180./pi.le.tholonw) then
[3329]	374	albedo(ig,:)=alb_tho_spe
[3275]	375	emis(ig)=emis_tho_spe
	376	endif
	377
	378	CASE (2) ! Depends on the albedo map
	379
[3329]	380	albedo(ig,:)=albedodat(ig)
	381	if (albedo(ig,1).gt.alb_ch4) then
[3275]	382	emis(ig)=emis_ch4
	383	else
	384	emis(ig)=emis_tho
	385	endif
	386
	387	CASE DEFAULT
	388	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
	389	stop
	390	END SELECT
	391
	392	!---------------------------------
	393	! 1.E. Tholins read from file
	394	!---------------------------------
	395
	396	if (specalb) then
[3329]	397	albedo(ig,:)=albedodat(ig) ! Specific ground properties
[3275]	398	!if (albedodat(ig).lt.0.25) then
[3329]	399	! albedo(ig,:)=alb_tho_eq
[3275]	400	!else if (albedodat(ig).lt.0.40) then
[3329]	401	! albedo(ig,:)=0.35
[3275]	402	!else if (albedodat(ig).lt.0.70) then
[3329]	403	! albedo(ig,:)=0.51
[3275]	404	!endif
	405	endif
	406
	407	endif ! ice_case
	408
	409	ENDDO ! ig ngrid
	410
	411	!-----------------------------------------------------------------------
	412	! 2) Changes of Thermal inertia with time
	413	!-----------------------------------------------------------------------
	414
	415	IF (changeti) then ! change of seasonal TI
	416	facls=8.
	417	DO ig=1,ngrid
	418
	419	! get depth of the different layers
	420	! limit diurnal / seasonal
[3539]	421	if (changetid.and.methane) then
[3275]	422	if (qsurf(ig,igcm_n2)>1.e-3) then
	423	emin=faclsITN2d/volcapa(daysec/pi)**0.5
	424	tid=ITN2d
	425	else if (qsurf(ig,igcm_ch4_ice)>1.e-3) then
	426	emin=faclsITCH4d/volcapa(daysec/pi)**0.5
	427	tid=ITCH4d
	428	else
	429	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
	430	tid=ITH2Od
	431	endif
[3539]	432	else if (changetid) then
	433	if (qsurf(ig,igcm_n2)>1.e-3) then
	434	emin=faclsITN2d/volcapa(daysec/pi)**0.5
	435	tid=ITN2d
	436	else
	437	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
	438	tid=ITH2Od
	439	endif
[3275]	440	else if (ngrid.ne.1) then
	441	emin=faclsinertiedat(ig,1)/volcapa(daysec/pi)**0.5
	442	tid=inertiedat(ig,1)
	443	else
	444	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
	445	tid=ITH2Od
	446	endif
	447
	448	! limit for N2
	449	emax=emin+qsurf(ig,igcm_n2)/1000.
	450
	451	! limit for CH4
	452	if (methane) then
	453	emax2=emax+qsurf(ig,igcm_ch4_ice)/1000.
	454	else
	455	emax2=0.
	456	endif
	457
	458	do isoil=0,nsoilmx-1
	459	if (mlayer(isoil).le.emin) then ! diurnal TI
[3483]	460	therm_inertia(ig,isoil+1)=tid
[3275]	461	else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
[3483]	462	therm_inertia(ig,isoil+1)=tid
[3275]	463	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
[3483]	464	therm_inertia(ig,isoil+1)=ITN2
[3275]	465	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
[3483]	466	therm_inertia(ig,isoil+1)=ITCH4
[3275]	467	else
[3483]	468	therm_inertia(ig,isoil+1)=ITH2O
[3275]	469	endif
	470
	471	enddo
	472	ENDDO
	473
	474	ELSE
	475
	476	DO ig=1,ngrid
[3483]	477	therm_inertia(ig,:)=inertiedat(ig,:)
[3275]	478	ENDDO
	479
	480	ENDIF
	481
	482
	483	!-----------------------------------------------------------------------
	484	! 3) Tests changements emis
	485	!-----------------------------------------------------------------------
	486	!n2cover=0.
	487	!n2coversun=0.
	488	!DO ig=1,ngrid
	489	! if (qsurf(ig,igcm_n2).ge.0.001) then
	490	! n2cover=n2cover+area(ig)
	491	! if (pfra(ig).gt.0.) then
	492	! n2coversun=n2coversun+area(ig)
	493	! endif
	494	! endif
	495	!enddo
	496	!gamm=n2cover/n2coversun
	497	!gamm=1.
	498	!tsb=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2bstephan))*(1/4.)
	499	!tsa=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2astephan))*(1/4.)
	500	!frac=max((35.6-tsb)/abs(tsa-tsb),0.)
	501	!write(,) 'gamm=',gamm,n2cover,n2coversun
	502	!write(,) 'tsb=',tsb
	503	!write(,) 'tsa=',tsa
	504	!write(55,*) 'frac=',frac,tsb,tsa
	505
	506
	507	end subroutine surfprop
	508

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

Download in other formats:

Original Format