Context Navigation

surfprop.F90 @ 3957

Last change on this file since 3957 was 3698, checked in by emillour, 9 months ago

Pluto PCM:
Some fixes to enable runnnig with dynamico:

add "strictboundcorrk" flag in callcorrk_pluto to enable running even if outside of kmatrix temperatures (when strictboundcorrk=.true.)
add premature exiting of writediagsoil if not with lon-lat grid
while at it, turned surfprop.F90 into a module
in physiq, enforce the possibility to output subsurface-related field in most cases, not only when "fast=.true."
adapt reference xml files: subsurface quantities need to be defined on a dedicated grid, otherwise XIOS will generate misleading garbage values. Updated files are put in "deftank/dynamico" for now.

File size: 21.2 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 DEFAULT
346	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
347	stop
348	END SELECT
349
350	!---------------------------------
351	! 1.D. THOLINS
352	!---------------------------------
353
354	else
355
356	SELECT CASE (mode_tholins)
357
358	CASE (0) ! Default, 2 albedos, one for the tropics, one for the poles
359
360	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
361	albedo(ig,:)=alb_tho_eq
362	emis(ig)=emis_tho_eq
363	else
364	albedo(ig,:)=alb_tho ! default = 0.1
365	emis(ig)=emis_tho
366	endif
367
368	CASE (1) ! Special mode one region with a different albedo
369
370	if (latitude(ig)180./pi.le.tholateq.and.latitude(ig)180./pi.gt.-tholateq) then
371	albedo(ig,:)=alb_tho_eq
372	emis(ig)=emis_tho_eq
373	else
374	albedo(ig,:)=alb_tho ! default = 0.1
375	emis(ig)=emis_tho
376	endif
377
378	if (latitude(ig)180./pi.le.tholatn.and.latitude(ig)180./pi.ge.tholats &
379	.and.longitude(ig)180./pi.ge.tholone.and.longitude(ig)180./pi.le.tholonw) then
380	albedo(ig,:)=alb_tho_spe
381	emis(ig)=emis_tho_spe
382	endif
383
384	CASE (2) ! Depends on the albedo map
385
386	albedo(ig,:)=albedodat(ig)
387	if (albedo(ig,1).gt.alb_ch4) then
388	emis(ig)=emis_ch4
389	else
390	emis(ig)=emis_tho
391	endif
392
393	CASE DEFAULT
394	write(,) 'STOP in surfprop.F90:mod_ch4 not found'
395	stop
396	END SELECT
397
398	!---------------------------------
399	! 1.E. Tholins read from file
400	!---------------------------------
401
402	if (specalb) then
403	albedo(ig,:)=albedodat(ig) ! Specific ground properties
404	!if (albedodat(ig).lt.0.25) then
405	! albedo(ig,:)=alb_tho_eq
406	!else if (albedodat(ig).lt.0.40) then
407	! albedo(ig,:)=0.35
408	!else if (albedodat(ig).lt.0.70) then
409	! albedo(ig,:)=0.51
410	!endif
411	endif
412
413	endif ! ice_case
414
415	ENDDO ! ig ngrid
416
417	!-----------------------------------------------------------------------
418	! 2) Changes of Thermal inertia with time
419	!-----------------------------------------------------------------------
420
421	IF (changeti) then ! change of seasonal TI
422	facls=8.
423	DO ig=1,ngrid
424
425	! get depth of the different layers
426	! limit diurnal / seasonal
427	if (changetid.and.methane) then
428	if (qsurf(ig,igcm_n2)>1.e-3) then
429	emin=faclsITN2d/volcapa(daysec/pi)**0.5
430	tid=ITN2d
431	else if (qsurf(ig,igcm_ch4_ice)>1.e-3) then
432	emin=faclsITCH4d/volcapa(daysec/pi)**0.5
433	tid=ITCH4d
434	else
435	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
436	tid=ITH2Od
437	endif
438	else if (changetid) then
439	if (qsurf(ig,igcm_n2)>1.e-3) then
440	emin=faclsITN2d/volcapa(daysec/pi)**0.5
441	tid=ITN2d
442	else
443	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
444	tid=ITH2Od
445	endif
446	else if (ngrid.ne.1) then
447	emin=faclsinertiedat(ig,1)/volcapa(daysec/pi)**0.5
448	tid=inertiedat(ig,1)
449	else
450	emin=faclsITH2Od/volcapa(daysec/pi)**0.5
451	tid=ITH2Od
452	endif
453
454	! limit for N2
455	emax=emin+qsurf(ig,igcm_n2)/1000.
456
457	! limit for CH4
458	if (methane) then
459	emax2=emax+qsurf(ig,igcm_ch4_ice)/1000.
460	else
461	emax2=0.
462	endif
463
464	do isoil=0,nsoilmx-1
465	if (mlayer(isoil).le.emin) then ! diurnal TI
466	therm_inertia(ig,isoil+1)=tid
467	else if (isoil.gt.0.and.(mlayer(isoil).gt.emin).and.(mlayer(isoil-1).lt.emin)) then ! still diurnal TI
468	therm_inertia(ig,isoil+1)=tid
469	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax)) then ! TI N2
470	therm_inertia(ig,isoil+1)=ITN2
471	else if ((mlayer(isoil).gt.emin).and.(mlayer(isoil).le.emax2)) then
472	therm_inertia(ig,isoil+1)=ITCH4
473	else
474	therm_inertia(ig,isoil+1)=ITH2O
475	endif
476
477	enddo
478	ENDDO
479
480	ELSE
481
482	DO ig=1,ngrid
483	therm_inertia(ig,:)=inertiedat(ig,:)
484	ENDDO
485
486	ENDIF
487
488
489	!-----------------------------------------------------------------------
490	! 3) Tests changements emis
491	!-----------------------------------------------------------------------
492	!n2cover=0.
493	!n2coversun=0.
494	!DO ig=1,ngrid
495	! if (qsurf(ig,igcm_n2).ge.0.001) then
496	! n2cover=n2cover+area(ig)
497	! if (pfra(ig).gt.0.) then
498	! n2coversun=n2coversun+area(ig)
499	! endif
500	! endif
501	!enddo
502	!gamm=n2cover/n2coversun
503	!gamm=1.
504	!tsb=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2bstephan))*(1/4.)
505	!tsa=(1/gammFat1AU/dist2(1.-alb_n2b)/(emis_n2astephan))*(1/4.)
506	!frac=max((35.6-tsb)/abs(tsa-tsb),0.)
507	!write(,) 'gamm=',gamm,n2cover,n2coversun
508	!write(,) 'tsb=',tsb
509	!write(,) 'tsa=',tsa
510	!write(55,*) 'frac=',frac,tsb,tsa
511
512
513	end subroutine surfprop
514
515	end module surfprop_mod

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

Download in other formats:

Original Format