Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

surfprop.F90 @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 19.5 KB

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

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

Download in other formats:

Original Format