Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
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

Line
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