Context Navigation

surfprop.F90 @ 3542

Last change on this file since 3542 was 3539, checked in by tbertrand, 13 months ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
File size: 21.2 KB

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

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

Download in other formats:

Original Format