Context Navigation

surfprop.F90 @ 3334

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

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

Download in other formats:

Original Format