Context Navigation

surfini.F @ 288

Last change on this file since 288 was 286, checked in by emillour, 13 years ago
Mars GCM: Fixed problem about undefined tracer names in 'surfini.F' by calling 'initracer' before 'surfini' in physiq.F EM
File size: 9.6 KB

Line
1	SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2	IMPLICIT NONE
3	c=======================================================================
4	c
5	c creation des calottes pour l'etat initial
6	c
7	c=======================================================================
8	c-----------------------------------------------------------------------
9	c Declarations:
10	c -------------
11	#include "dimensions.h"
12	#include "dimphys.h"
13	#include "surfdat.h"
14	#include "callkeys.h"
15	#include "tracer.h"
16	#include "comgeomfi.h"
17	#include "comcstfi.h"
18	c
19	INTEGER ngrid,ig,icap,iq,alternate
20	REAL piceco2(ngrid),psolaralb(ngrid,2)
21	REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
22
23	EXTERNAL ISMIN,ISMAX
24	INTEGER ISMIN,ISMAX
25	c
26	c=======================================================================
27
28	c water ice outliers
29	c ------------------------------------------
30
31	IF ((water) .and. (caps)) THEN
32
33	c Perennial H20 north cap defined by watercaptag=true (allows surface to be
34	c hollowed by sublimation in vdifc).
35
36	c We might not want albedodat to be modified because it is used to write
37	c restart files. Instead, albedo is directly modified when needed (i.e.
38	c if we have watercaptag and no co2 ice), below and in albedocaps.F90
39
40	c "Dryness coefficient" controlling the evaporation and
41	c sublimation from the ground water ice (close to 1)
42	c HERE, the goal is to correct for the fact
43	c that the simulated permanent water ice polar caps
44	c is larger than the actual cap and the atmospheric
45	c opacity not always realistic.
46
47	alternate = 0
48
49	do ig=1,ngridmx
50
51	#ifdef MESOSCALE
52	do iq=1,nqmx
53	write(,) "all qsurf to zero. dirty."
54	qsurf(ig,iq)=0. !! on jette les inputs GCM
55	!! on regle juste watercaptag
56	!! il faudrait garder les inputs GCM
57	!! si elles sont consequentes
58	enddo
59	if ( ( lati(ig)*180./pi .gt. 70. ) .and.
60	. ( albedodat(ig) .ge. 0.26 ) ) then
61	write(,)"outlier ",ig
62	watercaptag(ig) = .true.
63	dryness(ig) = 1.
64	albedodat(ig) = albedo_h2o_ice !! pour output
65	else
66	watercaptag(ig) = .false.
67	dryness(ig) = 1.
68	endif
69	#else
70
71	dryness (ig) = 1
72
73	!!c Towards olympia planitia water caps ('relatively' low latitude ones)
74	!!c---------------- proposition par AS --------------------
75	!!c--------------------------------------------------------
76	!!c if ( ( lati(ig)*180./pi .ge. 75 ) .and.
77	!!c . ( lati(ig)*180./pi .le. 77 ) .and.
78	!!c . ( ( ( long(ig)*180./pi .ge. 000. ) .and.
79	!!c . ( long(ig)*180./pi .le. 120. ) )
80	!!c . .or.
81	!!c . ( ( long(ig)*180./pi .ge. -130. ) .and.
82	!!c . ( long(ig)*180./pi .le. -115. ) ) ) ) then
83	!!c---------------- proposition par TN --------------------
84	!!c---------------- HIGHLY EXPERIMENTAL -------------------
85	!!c--------------------------------------------------------
86	!! if ( ( ( lati(ig)*180./pi .ge. 73. ) .and. ! cap #1
87	!! . ( lati(ig)*180./pi .le. 75.1 ) .and.
88	!! . ( long(ig)*180./pi .ge. 95. ) .and.
89	!! . ( long(ig)*180./pi .le. 110. ) )
90	!! . .or.
91	!! . ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2
92	!! . ( lati(ig)*180./pi .le. 80. ) .and.
93	!! . ( long(ig)*180./pi .ge. 110. ) .and.
94	!! . ( long(ig)*180./pi .le. 140. ) )
95	!! . .or.
96	!! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #3
97	!! . ( lati(ig)*180./pi .le. 78. ) .and.
98	!! . ( long(ig)*180./pi .ge. 155. ) .and.
99	!! . ( long(ig)*180./pi .le. 180. ) )
100	!!c . .or.
101	!!c . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater)
102	!!c . ( lati(ig)*180./pi .le. 72. ) .and.
103	!!c . ( long(ig)*180./pi .ge. 163. ) .and.
104	!!c . ( long(ig)*180./pi .le. 168. ) )
105	!! . .or.
106	!! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
107	!! . ( lati(ig)*180./pi .le. 78. ) .and.
108	!! . ( long(ig)*180./pi .ge. -160.) .and.
109	!! . ( long(ig)*180./pi .le. -120.) ) )
110	!! . then
111	!!
112	!! if (temptag) then
113	!!
114	!! if ((alternate .eq. 0)) then !!! 1/2 en 64x48 sinon trop large en lat
115	!! if (ngridmx.ne.1) watercaptag(ig)=.true.
116	!! write(,) "outliers ", lati(ig)*180./pi,
117	!! . long(ig)*180./pi
118	!! !dryness(ig) = 1.
119	!! alternate = 1
120	!! else
121	!! alternate = 0
122	!! endif !end if alternate = 0
123	!!
124	!! else
125	!!
126	!! if (ngridmx.ne.1) watercaptag(ig)=.true.
127	!! write(,) "outliers ", lati(ig)*180./pi,
128	!! . long(ig)*180./pi
129	!!
130	!! endif ! end if temptag
131	!!
132	!! endif
133	!!
134	!!
135	!!c Opposite olympia planitia water cap
136	!!c---------------- proposition par AS --------------------
137	!!c--------------------------------------------------------
138	!!c if ( ( lati(ig)*180./pi .ge. 82 ) .and.
139	!!c . ( lati(ig)*180./pi .le. 84 ) .and.
140	!!c . ( ( long(ig)*180./pi .gt. -030. ) .and.
141	!!c . ( long(ig)*180./pi .lt. 090. ) ) ) then
142	!!c---------------- proposition par TN --------------------
143	!!c--------------------------------------------------------
144	!! if ( ( lati(ig)*180./pi .ge. 80 ) .and.
145	!! . ( lati(ig)*180./pi .le. 84 ) .and.
146	!! . ( ( long(ig)*180./pi .gt. -030. ) .and.
147	!! . ( long(ig)*180./pi .lt. 090. ) ) ) then
148	!! if (ngridmx.ne.1) then
149	!! watercaptag(ig)=.true.
150	!! write(,) "central cap add ", lati(ig)*180./pi,
151	!! . long(ig)*180./pi
152	!! endif
153	!! !dryness(ig) = 1.
154	!! endif
155
156	c Central cap
157	c---------------- anciens reglages FF -------------------
158	c--------------------------------------------------------
159
160	if (lati(ig)*180./pi.gt.84) then
161	PRINT,'central cap', lati(ig)180./pi,
162	. long(ig)*180./pi
163	if (ngridmx.ne.1) watercaptag(ig)=.true.
164	!dryness(ig) = 1.
165	c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
166	c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
167	c if (ngridmx.ne.1) watercaptag(ig)=.true.
168	c dryness(ig) = 1.
169	c endif
170	c if (lati(ig)*180./pi.ge.85) then
171	c if (ngridmx.ne.1) watercaptag(ig)=.true.
172	c dryness(ig) = 1.
173	c endif
174	endif ! (lati>80 deg)
175	#endif
176	end do ! (ngridmx)
177	ENDIF ! (caps & water)
178
179	c ===============================================================
180	c INITIAL ALBEDO
181	c ===============================================================
182
183	write(,)"surfini: water frost thickness",
184	s frost_albedo_threshold
185	write(,)"surfini: water ice albedo:", albedo_h2o_ice
186	write(,)"surfini: water ice TI:", inert_h2o_ice
187
188	c To start with : Initial albedo = observed dataset
189	c -------------------------------------------------
190	DO ig=1,ngrid
191	psolaralb(ig,1)=albedodat(ig)
192	psolaralb(ig,2)=albedodat(ig)
193	END DO
194	PRINT*,'minimum albedo sans water caps',
195	s albedodat(ISMIN(ngrid,albedodat,1))
196	PRINT*,'maximum albedo sans water caps',
197	s albedodat(ISMAX(ngrid,albedodat,1))
198
199	c initial albedo if permanent H2O ice is present
200	c ------------------------------------------------
201	IF ((water) .and. (caps)) THEN
202	DO ig=1,ngrid
203	IF (watercaptag(ig)) THEN
204	psolaralb(ig,1) = albedo_h2o_ice
205	psolaralb(ig,2) = albedo_h2o_ice
206	ENDIF
207	END DO
208	PRINT*,'minimum albedo avec water caps',
209	s psolaralb(ISMIN(ngrid,psolaralb,1),1)
210	PRINT*,'maximum albedo avec water caps',
211	s psolaralb(ISMAX(ngrid,psolaralb,1),1)
212	ENDIF
213
214	c changing initial albedo if CO2 ice is present
215	c -------------------------------------------
216
217	DO ig=1,ngrid
218	IF (piceco2(ig) .GT. 0.) THEN
219	IF(ig.GT.ngrid/2+1) THEN
220	icap=2
221	ELSE
222	icap=1
223	ENDIF
224	psolaralb(ig,1) = albedice(icap)
225	psolaralb(ig,2) = albedice(icap)
226	END IF
227	END DO
228
229	c changing initial albedo if water ice frost is present
230	c -------------------------------------------
231	IF (water) THEN
232	do iq=1,nqmx
233	c if there is frost and surface albedo is set to albedo_h2o_ice
234	if(noms(iq).eq."h2o_ice") then
235	do ig=1,ngrid
236	if ((piceco2(ig) .eq. 0.).and.
237	& (qsurf(ig,iq).gt.frost_albedo_threshold)) then
238	psolaralb(ig,1) = albedo_h2o_ice
239	psolaralb(ig,2) = albedo_h2o_ice
240	c PRINT*,'surfini.F frost',
241	c & lati(ig)180./pi, long(ig)180./pi
242	endif
243	enddo
244	endif
245	end do
246	PRINT*,'minimum albedo avec givre et co2',
247	s psolaralb(ISMIN(ngrid,psolaralb,1),1)
248	PRINT*,'maximum albedo avec givre et co2',
249	s psolaralb(ISMAX(ngrid,psolaralb,1),1)
250	END IF
251
252
253	RETURN
254	END

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

Download in other formats:

Original Format