Context Navigation

surfprop.F90 @ 3421

Last change on this file since 3421 was 3377, checked in by afalco, 19 months ago
Pluto PCM: initialize some variables to 0. 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	emis(:) = 1 !default to 1
137
138	! Loop on all points
139	DO ig=1,ngrid
140
141	! Looking for dominant ice:
142	ice_case=0 ! Tholins
143	if (methane) then
144	if (qsurf(ig,igcm_ch4_ice).ge.thres_ch4ice) then
145	ice_case=1
146	endif
147	endif
148	if (carbox) then
149	if (qsurf(ig,igcm_co_ice).ge.thres_coice) then
150	ice_case=2
151	endif
152	endif
153	if (qsurf(ig,igcm_n2).ge.thres_n2ice) then
154	ice_case=3
155	endif
156
157	!---------------------------------
158	! 1.A. N2
159	!---------------------------------
160
161	if (ice_case.eq.3) then
162
163	! EMISSIVITY N2
164	! change emis if phase alpha different de phase beta
165	frac=0.5(1.+tanh(6.(35.6-tsurface(ig))/0.5))
166	if (frac.gt.1.) frac=1.
167	emis(ig)=fracemis_n2a+(1.-frac)emis_n2b
168
169	! ALBEDO N2 and emissivity changes
170	SELECT CASE (mode_n2)
171
172	CASE (0) ! fixed albedo
173	albedo(ig,:)=min(max(alb_n2b+adjust,0.),1.)
174
175	CASE (1) ! Albedo decreases with thickness
176	albedo(ig,:)=(alb_n2b-deltab)/(1.-10000.)*qsurf(ig,igcm_n2)+alb_n2b
177	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
178	CASE (2) ! Special Sputnik differences of albedo
179	if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)180./pi.le.-161.or.longitude(ig)180./pi.ge.146)) then
180	if ( (latitude(ig)*180./pi.ge.25.).or. &
181	((longitude(ig)*180./pi.gt.140.).and. &
182	(longitude(ig)*180./pi.lt.165.)) ) then
183	albedo(ig,:)=alb_n2b-deltab
184	else
185	albedo(ig,:)=alb_n2b+deltab
186	endif
187	else
188	albedo(ig,:)=alb_n2b
189	endif
190
191	CASE (3) ! Albedo increases (delta neg) or decreases (delta pos) with sublimation rates
192	albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
193	albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
194	!! Triton:
195	!albedo(ig,:)=deltab/1.e-4*fluold(ig,igcm_n2)+albedo(ig,:)
196	!albedo(ig,:)=min(max(alb_n2b-deltab,albedo(ig,:)),alb_n2b+deltab)
197
198	CASE (4) ! Albedo Difference in N/S
199	if (latitude(ig)*180./pi.ge.0.) then
200	albedo(ig,:)=min(max(alb_n2b-deltab+adjust,0.),1.)
201	else
202	albedo(ig,:)=min(max(alb_n2b+deltab+adjust,0.),1.)
203	endif
204
205	CASE (5) ! Special Sputnik differences of albedo in patches
206	!! Patches Nogcm
207	!if (qsurf(ig,igcm_n2).ge.1.e4.and.(longitude(ig)180./pi.le.180.).and.(longitude(ig)180./pi.ge.174.) ) then
208	! if (((latitude(ig)180./pi.le.46.).and.(latitude(ig)180./pi.ge.42.)) &
209	! .or. ((latitude(ig)180./pi.le.36.).and.(latitude(ig)180./pi.ge.32.)) &
210	! .or. ((latitude(ig)180./pi.le.26.).and.(latitude(ig)180./pi.ge.22.)) &
211	! .or. ((latitude(ig)180./pi.le.16.).and.(latitude(ig)180./pi.ge.12.)) &
212	! .or. ((latitude(ig)180./pi.le.6.).and.(latitude(ig)180./pi.ge.2.)) &
213	! ) then
214	! albedo(ig,:)=alb_n2b-deltab
215
216	!! Patches GCM
217	if (qsurf(ig,igcm_n2).ge.1.e4) then
218	if (((latitude(ig)180./pi.lt.33.).and.(latitude(ig)180./pi.gt.32.).and. &
219	(longitude(ig)180./pi.lt.165.5).and.(longitude(ig)180./pi.gt.164.5)) &
220	.or. ((latitude(ig)180./pi.lt.30.5).and.(latitude(ig)180./pi.gt.29.5).and. &
221	(longitude(ig)180./pi.lt.169.).and.(longitude(ig)180./pi.gt.168.))) then
222	albedo(ig,:)=alb_n2b-deltab
223	else if (((latitude(ig)180./pi.lt.30.5).and.(latitude(ig)180./pi.gt.29.).and. &
224	(longitude(ig)180./pi.lt.165.5).and.(longitude(ig)180./pi.gt.164.5)) &
225	.or. ((latitude(ig)180./pi.lt.33.).and.(latitude(ig)180./pi.gt.32.).and. &
226	(longitude(ig)180./pi.lt.169.).and.(longitude(ig)180./pi.gt.168.))) then
227	albedo(ig,:)=alb_n2b+deltab
228	else
229	albedo(ig,:)=alb_n2b
230	endif
231	else
232	albedo(ig,:)=alb_n2b
233	endif
234
235	CASE (7) ! Albedo from albedodat and adjusted emissivity
236	albedo(ig,:)=albedodat(ig)
237	! adjust emissivity if convergeps = true
238	emis(ig)=min(max(emis(ig)+adjust,0.),1.)
239
240	CASE DEFAULT
241	write(,) 'STOP in surfprop.F90: mod_n2 not found'
242	stop
243	END SELECT
244
245
246	!---------------------------------
247	! 1.B. CO
248	!---------------------------------
249
250	else if (ice_case.eq.2) then
251	albedo(ig,:)=alb_co
252	emis(ig)=emis_co
253
254	!---------------------------------
255	! 1.C. CH4
256	!---------------------------------
257
258	else if (ice_case.eq.1) then
259
260	SELECT CASE (mode_ch4)
261
262	CASE (0) ! 2 albedos, one for the tropics, one for the poles
263	emis(ig)=emis_ch4
264	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
265	albedo(ig,:)=alb_ch4_eq
266	else
267	albedo(ig,:)=alb_ch4
268	endif
269
270	CASE (1) ! 3 albedos, one for the tropics, 2 for the poles (north and south)
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 if (latitude(ig)*180./pi.le.-metlateq) then
275	albedo(ig,:)=alb_ch4_s
276	else
277	albedo(ig,:)=alb_ch4
278	endif
279
280	CASE (2) ! 2 albedos + albedo feedback
281
282	emis(ig)=emis_ch4
283	albedo(ig,:)=alb_ch4
284
285	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
286	albedo(ig,:)=alb_ch4_eq
287	endif
288
289	!! Albedo feedback
290	if (latitude(ig)180./pi.le.fdch4_latn.and.latitude(ig)180./pi.gt.fdch4_lats) then
291	if (longitude(ig)180./pi.le.fdch4_lone.and.longitude(ig)180./pi.gt.fdch4_lonw) then
292	if (qsurf(ig,igcm_ch4_ice).lt.fdch4_maxice) then ! fd would not apply on BTD
293	SELECT CASE (feedback_met)
294	CASE (0) ! Default (linear from alb_ch4_eq)
295	albedo(ig,:)=(1.-alb_ch4_eq)/0.002*qsurf(ig,igcm_ch4_ice)+alb_ch4_eq
296	!emis(ig)=(1.-emis_ch4)/0.002*qsurf(ig,igcm_ch4_ice)+emis_ch4
297	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
298	!if (emis(ig).gt.1.) emis(ig)=1.
299	CASE (1) ! Hyperbolic tangent old
300	albedo(ig,:)=0.60.5(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.3 !
301	!emis(ig)=0.20.5(1.+tanh(6.*(0.001+qsurf(ig,igcm_ch4_ice))/0.5))+0.8 !
302	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
303	!if (emis(ig).gt.1.) emis(ig)=1.
304	CASE (2) ! hyperbolic tangent old
305	albedo(ig,:)=0.50.6(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.3 !
306	!emis(ig)=0.20.5(1.+tanh(1000.*(qsurf(ig,igcm_ch4_ice)-0.002)))+0.8 !
307	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
308	!if (emis(ig).gt.1.) emis(ig)=1.
309	CASE (3) ! hyperbolic tangent equation with parameters
310	albedo(ig,:)=aa(-1+tanh(fdch4_amplqsurf(ig,igcm_ch4_ice)))+bb
311	if (albedo(ig,1).gt.fdch4_maxalb) albedo(ig,:)=fdch4_maxalb
312	CASE DEFAULT
313	write(,) 'STOP surfprop.F90: fd wrong'
314	stop
315	END SELECT
316	endif
317	endif
318	endif
319
320	CASE (3) ! Eq, poles N, pole S + depending on Ls
321	emis(ig)=emis_ch4
322	if (latitude(ig)180./pi.le.metlateq.and.latitude(ig)180./pi.gt.-metlateq) then
323	albedo(ig,:)=alb_ch4_eq ! Pure methane ice
324	else if (latitude(ig)*180./pi.le.-metlateq) then
325	albedo(ig,:)=alb_ch4_s ! Pure methane ice
326	if (pls.le.metls2.and.pls.gt.metls1) then
327	albedo(ig,:)=alb_ch4_s+deltab ! Also used for N2, careful
328	endif
329	else
330	albedo(ig,:)=alb_ch4 ! Pure methane ice
331	endif
332
333	CASE (4) ! Albedo from albedodat
334	emis(ig)=emis_ch4
335	albedo(ig,:)=albedodat(ig)
336
337	CASE DEFAULT
338	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
339	stop
340	END SELECT
341
342	!---------------------------------
343	! 1.D. THOLINS
344	!---------------------------------
345
346	else
347
348	SELECT CASE (mode_tholins)
349
350	CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
351
352	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
353	albedo(ig,:)=alb_tho_eq
354	emis(ig)=emis_tho_eq
355	else
356	albedo(ig,:)=alb_tho ! default = 0.1
357	emis(ig)=emis_tho
358	endif
359
360	CASE (1) ! Special mode one region with a different albedo
361
362	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
363	albedo(ig,:)=alb_tho_eq
364	emis(ig)=emis_tho_eq
365	else
366	albedo(ig,:)=alb_tho ! default = 0.1
367	emis(ig)=emis_tho
368	endif
369
370	if (latitude(ig)180./pi.le.tholatn.and.latitude(ig)180./pi.ge.tholats &
371	.and.longitude(ig)180./pi.ge.tholone.and.longitude(ig)180./pi.le.tholonw) then
372	albedo(ig,:)=alb_tho_spe
373	emis(ig)=emis_tho_spe
374	endif
375
376	CASE (2) ! Depends on the albedo map
377
378	albedo(ig,:)=albedodat(ig)
379	if (albedo(ig,1).gt.alb_ch4) then
380	emis(ig)=emis_ch4
381	else
382	emis(ig)=emis_tho
383	endif
384
385	CASE DEFAULT
386	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
387	stop
388	END SELECT
389
390	!---------------------------------
391	! 1.E. Tholins read from file
392	!---------------------------------
393
394	if (specalb) then
395	albedo(ig,:)=albedodat(ig) ! Specific ground properties
396	!if (albedodat(ig).lt.0.25) then
397	! albedo(ig,:)=alb_tho_eq
398	!else if (albedodat(ig).lt.0.40) then
399	! albedo(ig,:)=0.35
400	!else if (albedodat(ig).lt.0.70) then
401	! albedo(ig,:)=0.51
402	!endif
403	endif
404
405	endif ! ice_case
406
407	ENDDO ! ig ngrid
408
409	!-----------------------------------------------------------------------
410	! 2) Changes of Thermal inertia with time
411	!-----------------------------------------------------------------------
412
413	IF (changeti) then ! change of seasonal TI
414	facls=8.
415	DO ig=1,ngrid
416
417	! get depth of the different layers
418	! limit diurnal / seasonal
419	if (changetid) then
420	if (qsurf(ig,igcm_n2)>1.e-3) then
421	emin=faclsITN2d/volcapa(daysec/pi)**0.5
422	tid=ITN2d
423	else if (qsurf(ig,igcm_ch4_ice)>1.e-3) then
424	emin=faclsITCH4d/volcapa(daysec/pi)**0.5
425	tid=ITCH4d
426	else
427	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
428	tid=ITH2Od
429	endif
430	else if (ngrid.ne.1) then
431	emin=faclsinertiedat(ig,1)/volcapa(daysec/pi)**0.5
432	tid=inertiedat(ig,1)
433	else
434	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
435	tid=ITH2Od
436	endif
437
438	! limit for N2
439	emax=emin+qsurf(ig,igcm_n2)/1000.
440
441	! limit for CH4
442	if (methane) then
443	emax2=emax+qsurf(ig,igcm_ch4_ice)/1000.
444	else
445	emax2=0.
446	endif
447
448	do isoil=0,nsoilmx-1
449	if (mlayer(isoil).le.emin) then ! diurnal TI
450	tidat(ig,isoil+1)=tid
451	else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
452	tidat(ig,isoil+1)=tid
453	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
454	tidat(ig,isoil+1)=ITN2
455	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
456	tidat(ig,isoil+1)=ITCH4
457	else
458	tidat(ig,isoil+1)=ITH2O
459	endif
460
461	enddo
462	ENDDO
463
464	ELSE
465
466	DO ig=1,ngrid
467	tidat(ig,:)=inertiedat(ig,:)
468	ENDDO
469
470	ENDIF
471
472
473	!-----------------------------------------------------------------------
474	! 3) Tests changements emis
475	!-----------------------------------------------------------------------
476	!n2cover=0.
477	!n2coversun=0.
478	!DO ig=1,ngrid
479	! if (qsurf(ig,igcm_n2).ge.0.001) then
480	! n2cover=n2cover+area(ig)
481	! if (pfra(ig).gt.0.) then
482	! n2coversun=n2coversun+area(ig)
483	! endif
484	! endif
485	!enddo
486	!gamm=n2cover/n2coversun
487	!gamm=1.
488	!tsb=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2bstephan))*(1/4.)
489	!tsa=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2astephan))*(1/4.)
490	!frac=max((35.6-tsb)/abs(tsa-tsb),0.)
491	!write(,) 'gamm=',gamm,n2cover,n2coversun
492	!write(,) 'tsb=',tsb
493	!write(,) 'tsa=',tsa
494	!write(55,*) 'frac=',frac,tsb,tsa
495
496
497	end subroutine surfprop
498

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

Download in other formats:

Original Format