Context Navigation

surfprop.F90 @ 3986

Last change on this file since 3986 was 3986, checked in by vbelissa, 2 weeks ago
PLUTO PCM : Adding option for CH4 albedo VB
File size: 21.9 KB

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

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

Download in other formats:

Original Format