source: trunk/LMDZ.MARS/libf/phymars/geticecover.F90 @ 4016

Last change on this file since 4016 was 3902, checked in by emillour, 9 months ago

Mars PCM:
More code tidying: turn geticecover, interp_line, orbite, simpleclouds, tabfi
and tcondco2 into modules.
EM
File size: 14.2 KB
Line 
1MODULE geticecover_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7SUBROUTINE geticecover( ngrid, ls, lontab, lattab, icecover )
8  !!**********************************************!!
9  ! A. Spiga (aymeric.spiga@upmc.fr) October 2011  !
10  !!**********************************************!!
11  IMPLICIT NONE
12  !! INPUTS
13  REAL, DIMENSION(ngrid), INTENT(IN)  :: lontab,lattab
14  REAL,                   INTENT(IN)  :: ls
15  INTEGER,                INTENT(IN)  :: ngrid
16  !! OUTPUTS
17  REAL, DIMENSION(ngrid), INTENT(OUT) :: icecover
18  !! LOCAL
19  INTEGER :: ig
20  icecover(:) = 0.
21  DO ig=1,ngrid
22    icecover(ig) = isitice( ls, lontab(ig), lattab(ig) )
23  ENDDO
24END SUBROUTINE geticecover
25
26REAL FUNCTION isitice( ls, lon, lat )
27  !!**********************************************!!
28  ! A. Spiga (aymeric.spiga@upmc.fr) October 2011  !
29  !!**********************************************!!
30  IMPLICIT NONE
31  REAL, INTENT(IN) :: ls,lon,lat
32  isitice = 0.
33  !! for speedup purposes:
34  !! useless to call functions for lats between -40. and 40.
35  IF (lat .le. -40.) THEN
36    IF ( lat .le. splatcrocus(ls,lon) ) isitice = 1.
37  ENDIF
38  IF (lat .ge. 40.) THEN
39    IF ( lat .ge. nplatcrocus(ls,lon) ) isitice = 1.
40  ENDIF
41END FUNCTION isitice
42
43REAL FUNCTION nplatcrocus( ls, lon )
44  !;Purpose: To return the areocentric latitude of the IR cap edge for TES year "1".
45  !;
46  !;Inputs:
47  !;  ls: Season in degrees of L_s.
48  !;  lo: East Longitude of interest.
49  !;
50  !;Dates of TES data used:  1999 FEB 28 to  2001 JAN 15
51  !;Dates of TES data used:  2001 JAN 15 to  2002 DEC 2
52  !;Dates of TES data used:  2002 DEC 2 to  1980 JAN 1
53  !;
54  !;Written by T.N. Titus Thu Dec 15 11:54:26 2005
55  !;U.S. Geological Survey Astrogeology Team.
56  !;2255 North Gemini Drive
57  !;Flagstaff, AZ 86001 USA
58  !;http://www.mars-ice.org
59  !;
60  !;If using this function for research or publication, please cite:
61  !; T. N. Titus (2005), Mars Polar Cap Edges Tracked over 3 Full Mars Years,
62  !;   36th Annual Lunar and Planetary Science Conference, March 14-18, 2005,
63  !;   in League City, Texas, abstract no.1993
64  !!**********************************************************!!
65  ! ADAPTATION: A. Spiga (aymeric.spiga@upmc.fr) October 2011  !
66  !!**********************************************************!!
67  IMPLICIT NONE
68  !! INPUT
69  REAL, INTENT(IN) :: ls,lon
70  !! LOCAL
71  REAL, DIMENSION(40) :: gen1, gen2, quan, c, s
72  REAL, DIMENSION(40,2,3) :: tes
73  INTEGER :: year
74  REAL, DIMENSION(3) :: const, line
75  !!*****************************************************************************************
76  !! CONSTANT
77  const(:) = (/ 71.4648,70.4607,71.5339 /)
78  !! GENERIC TABS
79  gen1(:) = (/ &
80  & 0.00000,1.00000,2.00000,3.00000,4.00000, &
81  & 0.00000,1.00000,2.00000,3.00000,4.00000, &
82  & 0.00000,1.00000,2.00000,3.00000,4.00000, &
83  & 0.00000,1.00000,2.00000,3.00000,4.00000, &
84  & 1.00000,2.00000,3.00000,4.00000, &
85  & 1.00000,2.00000,3.00000,4.00000, &
86  & 1.00000,2.00000,3.00000,4.00000, &
87  & 1.00000,2.00000,3.00000,4.00000, &
88  & 1.00000,2.00000,3.00000,4.00000 &
89  /)
90  gen2(:) = (/ &
91  & 1.00000,1.00000,1.00000,1.00000,1.00000, &
92  & 2.00000,2.00000,2.00000,2.00000,2.00000, &
93  & 3.00000,3.00000,3.00000,3.00000,3.00000, &
94  & 4.00000,4.00000,4.00000,4.00000,4.00000, &
95  &  0.00000, 0.00000, 0.00000, 0.00000, &
96  & -1.00000,-1.00000,-1.00000,-1.00000, &
97  & -2.00000,-2.00000,-2.00000,-2.00000, &
98  & -3.00000,-3.00000,-3.00000,-3.00000, &
99  & -4.00000,-4.00000,-4.00000,-4.00000 &
100  /)
101  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NORTHERN CAP
102  !!! FIRST YEAR TES
103  tes(:,1,1) = (/ &
104  & -1.04928,  -0.217329,  0.101174, -0.0374879, 0.00576928, &
105  & -0.467708,  0.213426, -0.322449,  0.248625, -0.0879115, &
106  & -0.0852889,-0.0280091,-0.185435,  0.0160633, 0.0485029, &
107  & -0.0402631, 0.163700, -0.115792, -0.0468377,-0.00756858, &
108  & -11.1382,   0.918513,  0.754372, -0.288945,  0.0435517, &
109  & 0.143201,  -0.195657,  0.0378341,-0.640119,  0.429567, &
110  & -0.114784, -0.0829540,-0.179708,  0.0328001, 0.0342086, &
111  & -0.0193232, 0.0279848, 0.180155, -0.154807,  0.0971645 &
112  /)
113  tes(:,2,1) = (/ &
114  &  0.673412,   0.194463,  0.0213440, -0.219141,   0.124593, &
115  &  0.0526366, -0.604908,  0.113683,   0.0372814, -0.0714123, &
116  & -0.00710610,-0.0933308, 0.0241633,  0.0559299, -0.0924497, &
117  &  0.0747056,  0.0400294,-0.127268,  -0.0222301,  0.0272872, &
118  & -19.2296,    4.86860,  -0.454119,  -0.375147,  -0.119044, &
119  &  0.209839,  -0.271209,  0.367698,   0.171835,   0.350654, &
120  & -0.288149,  -0.0152939, 0.0195934,  0.184573,  -0.0272326, &
121  & -0.0466995,  0.299228,- 0.0593913, -0.0660705,  0.0282595 &
122  /)
123  !!! SECOND YEAR TES
124  tes(:,1,2) = (/ &
125  & -0.728762,-0.671618,0.154561,0.112470,-0.152911,-0.633290, &
126  &  0.372959,-0.241842,0.100173,0.0558815,-0.0258014,0.0138332, &
127  & -0.113160,0.0965581,-0.0478080,0.0199076,0.0360519,-0.119565, &
128  & -0.0564988,0.00166589,-10.8664,1.34367,0.334286,-0.0889369, &
129  & 0.170191,-0.241918,0.0406942,-0.0312268,-0.542561,0.440734,&
130  & -0.181151,-0.178009,-0.216867,-0.00154585,0.0381691,0.0333462,&
131  & 0.0558955,0.0512524,-0.0389391,0.0705681 &
132  /)
133  tes(:,2,2) = (/ &
134  & 0.133511,0.301999,0.354979,-0.471594,0.151863,0.0575986, &
135  & -0.478402,-0.0168485,0.148232,-0.128787,-0.0914331,-0.0603194, &
136  & 0.0611823,0.0668322,-0.0833629,0.0176089,0.0570124,-0.0641502, &
137  & 0.0196993,-0.0272954,-18.8258,4.08857,-0.252963,-0.0642321, &
138  & -0.463212,0.380846,-0.00712758,0.0374067,0.230256,0.154810, &
139  & -0.235321,0.0838983,0.0737168,0.115381,-0.00787243,-0.0289239, &
140  & 0.169424,-0.00376776,-0.0134735,0.00663328 &
141  /)
142  !!! THIRD YEAR TES
143  tes(:,1,3) = (/ &
144  & -0.170075,-0.706807,-0.0725778,0.266425,-0.166670,-0.408371, &
145  & 0.106833,-0.318539,0.319422,-0.0361432,-0.0366784,-0.0231430, &
146  & -0.129665,0.0478013,0.0283050,0.399130,-0.177498,-0.113099, &
147  & 0.0885522,-0.0503254,-12.2470,1.40947,0.401761,-0.149156, &
148  & -0.256222,-0.0307999,0.173418,-0.112992,-0.628345,0.286791, &
149  & -0.0395388,-0.0634619,-0.263644,0.0727097,0.0679133,0.0350731, &
150  & -0.140744,0.0673651,0.0457041,-0.0454068 &
151  /)
152  tes(:,2,3) = (/ &
153  & 0.447200,-0.251112,0.597279,-0.374734,-0.0290520,-0.0208475, &
154  & -0.682811,0.146231,0.0122089,-0.0917365,0.00708168,-0.143992, &
155  & 0.0724996,0.102962,-0.0703368,-0.0429708,-0.173888,0.189929, &
156  & -0.107449,-0.0505178,-19.3010,4.83420,-1.09834,0.430944, &
157  & -0.643534,0.710448,-0.299251,0.135935,0.134539,0.320505, &
158  & -0.200302,-0.0106718,0.0166824,0.260126,-0.0992902,-0.0182585, &
159  & -0.0627401,0.232410,-0.126134,-0.0100335 &
160  /)
161  !!! MAIN CALCULATIONS
162  quan = (acos(-1.)/180.)*ls*gen1 + (acos(-1.)/180.)*lon*gen2
163  c = cos(quan)
164  s = sin(quan)
165  line(:) = 0.
166  DO year=1,3
167    line(year) = DOT_PRODUCT(c,tes(:,1,year)) - DOT_PRODUCT(s,tes(:,2,year)) + const(year)
168    !!! this is taken into account in isitice
169    !if (line(year) < -90.) line(year) = -90.
170    !if (line(year) >  90.) line(year) =  90.
171  ENDDO
172  !!! THIS IS TEMPORARY, BUT WHY NOT FOR MODELING
173  nplatcrocus = (line(1)+line(2)+line(3))/3.
174END FUNCTION nplatcrocus
175
176REAL FUNCTION splatcrocus( ls, lon )
177  !;Purpose: To return the areocentric latitude of the IR cap edge for TES year "1".
178  !;
179  !;Inputs:
180  !;  ls: Season in degrees of L_s.
181  !;  lo: East Longitude of interest.
182  !;
183  !;Dates of TES data used:  1999 FEB 28 to  2001 JAN 15
184  !;Dates of TES data used:  2001 JAN 15 to  2002 DEC 2
185  !;Dates of TES data used:  2002 DEC 2 to  1980 JAN 1
186  !;
187  !;Written by T.N. Titus Thu Dec 15 11:54:26 2005
188  !;U.S. Geological Survey Astrogeology Team.
189  !;2255 North Gemini Drive
190  !;Flagstaff, AZ 86001 USA
191  !;http://www.mars-ice.org
192  !;
193  !;If using this function for research or publication, please cite:
194  !; T. N. Titus (2005), Mars Polar Cap Edges Tracked over 3 Full Mars Years,
195  !;   36th Annual Lunar and Planetary Science Conference, March 14-18, 2005,
196  !;   in League City, Texas, abstract no.1993
197  !!**********************************************************!!
198  ! ADAPTATION: A. Spiga (aymeric.spiga@upmc.fr) October 2011  !
199  !!**********************************************************!!
200  IMPLICIT NONE
201  !! INPUT
202  REAL, INTENT(IN) :: ls,lon
203  !! LOCAL
204  REAL, DIMENSION(60) :: gen1, gen2, quan, c, s
205  REAL, DIMENSION(60,2,2) :: tes
206  INTEGER :: year
207  REAL, DIMENSION(2) :: const, line
208  !!*****************************************************************************************
209  !! CONSTANT
210  const(:) = (/ -69.6039, -68.8420 /)
211  !! GENERIC TABS
212  gen1(:) = (/ &
213  & 0.00000,1.00000,2.00000,3.00000,4.00000, &
214  & 5.00000,0.00000,1.00000,2.00000,3.00000, &
215  & 4.00000,5.00000,0.00000,1.00000,2.00000, &
216  & 3.00000,4.00000,5.00000,0.00000,1.00000, &
217  & 2.00000,3.00000,4.00000,5.00000,0.00000, &
218  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
219  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
220  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
221  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
222  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
223  & 1.00000,2.00000,3.00000,4.00000,5.00000, &
224  & 1.00000,2.00000,3.00000,4.00000,5.00000 &
225  /)
226  gen2(:) = (/ &
227  &  1.00000, 1.00000, 1.00000, 1.00000, 1.00000, &
228  &  1.00000, 2.00000, 2.00000, 2.00000, 2.00000, &
229  &  2.00000, 2.00000, 3.00000, 3.00000, 3.00000, &
230  &  3.00000, 3.00000, 3.00000, 4.00000, 4.00000, &
231  &  4.00000, 4.00000, 4.00000, 4.00000, 5.00000, &
232  &  5.00000, 5.00000, 5.00000, 5.00000, 5.00000, &
233  &  0.00000, 0.00000, 0.00000, 0.00000, 0.00000, &
234  & -1.00000,-1.00000,-1.00000,-1.00000,-1.00000, &
235  & -2.00000,-2.00000,-2.00000,-2.00000,-2.00000, &
236  & -3.00000,-3.00000,-3.00000,-3.00000,-3.00000, &
237  & -4.00000,-4.00000,-4.00000,-4.00000,-4.00000, &
238  & -5.00000,-5.00000,-5.00000,-5.00000,-5.00000 &
239  /)
240  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOUTHERN CAP
241  !!! FIRST YEAR TES
242  tes(:,1,1) = (/ &
243  &  3.53762,  -1.57584,  -2.97332,  -1.07670,    0.272623,  &
244  &  0.311819,  0.842387,  1.11375,   0.388617,  -0.575892,  &
245  & -0.502105, -0.0357803,-0.689386,  0.743998,   0.796186,  &
246  &  0.304092, -0.172464, -0.290923, -0.542742,  -0.577187,  &
247  & -0.439117,  0.266552,  0.393525,  0.0654230,  0.296123,  &
248  &  0.0949434,-0.162267, -0.175361,  0.000162221,0.105111,  &
249  & -12.4148,   0.177898,  1.45822,   0.309320,  -0.148393, 3.51232, &
250  &  0.355009, -0.782760, -0.427043, -0.178686,  -0.429962,  &
251  & -0.649747, -0.487551,  0.259459,  0.139186,  -0.964653,  &
252  & -0.409536,  0.137911,  0.293971,  0.0248885,  0.238855,  &
253  &  0.612519,  0.172469, -0.0462983,-0.0129877,  0.00532444,&
254  & -0.189464, -0.139035, -0.0231797, 0.0903127 &
255  /)
256  tes(:,2,1) = (/ &
257  &  3.21178,   3.81330,    0.357977,  -1.85466,   -1.01477,  &
258  & -0.209586, -0.797844,   0.334659,   0.828534,   0.519493, &
259  & -0.0512511,-0.204263,  -0.924137,  -0.694012, -0.0666508, &
260  &  0.632104,  0.397041,  -0.00521332, 0.933057,  0.0136105, &
261  & -0.646719, -0.480330,  -0.0122939,  0.214389,  0.0880269, &
262  &  0.0896852, 0.158164,  -0.00717643,-0.130884, -0.0425004, &
263  & -21.5000,  -5.88472,   -0.767056,   0.261383,  0.232822,  &
264  &  1.18087,   2.13123,    1.00744,   -0.0527458,-0.0417051, &
265  &  1.05940,   0.00149598,-0.501423,  -0.290405, -0.0562713, &
266  &  0.335254, -0.480024,  -0.604500,   0.0411040, 0.0403603, &
267  & -0.489831,  0.109007,   0.217500,   0.109769, -0.000845153, &
268  &  0.152139,  0.206418,  -0.0980599, -0.180992, -0.00685824 &
269  /)
270  !!! SECOND YEAR TES
271  tes(:,1,2) = (/ &
272  &  3.10733,  -0.0926926, -1.38769,   -0.638410,  0.0598460, &
273  &  0.0212056,-0.574512,  -0.687073,  -0.193427,  0.161192,  &
274  &  0.242137,  0.0178444, -0.427617,   0.892971,  0.650940,  &
275  & -0.0372970,-0.311147,  -0.123042,  -0.221688, -0.467191,  &
276  & -0.659134,  0.0246518,  0.257168,   0.124850,  0.174360,  &
277  &  0.320583,  0.136864,  -0.0606024, -0.0941828, 0.00372431,&
278  & -11.1678,  -0.376356,   0.316746,  -0.100600,  -0.119510, &
279  &  1.35807,  -1.53000,   -0.867922,   0.301309,   0.258290, &
280  &  0.0813258, 0.601762,   0.268272,   0.0989900, -0.0737015,&
281  & -0.968944, -0.599203,   0.217550,   0.466813,   0.0202996,&
282  &  0.693900,  0.653567,  -0.0888332, -0.219429,  -0.0298150,&
283  & -0.186738, -0.329258,  -0.155475,   0.0241376,  0.0613455 &
284  /)
285  tes(:,2,2) = (/ &
286  &  0.745574,  2.34362,    0.617967,  -0.847462,  -0.440483, &
287  & -0.135479,  0.546528,  -0.0659467, -0.627697,  -0.358895, &
288  &  0.0292012, 0.0797757, -0.953241,  -0.386098,   0.325419, &
289  &  0.543846,  0.108148,  -0.0596165,  1.15557,    0.436356, &
290  & -0.392191, -0.500867,  -0.205739,   0.194082,  -0.239364, &
291  & -0.136630,  0.135228,   0.213587,   0.0299697, -0.0101247,&
292  & -19.9097,  -5.34634,    0.103437,  -0.343399,  -0.348643, &
293  &  2.20957,   1.02330,   -0.722428,  -0.574345,   0.149871, &
294  & -0.635328, -0.317596,   0.154507,   0.110691,   0.0265673,&
295  &  0.593630, -0.516898  ,-0.689954,   0.00465230, 0.120493, &
296  & -0.323229,  0.491035,   0.508041,   0.0373988, -0.106504, &
297  &  0.241758,  0.0470212, -0.236231,  -0.186565,   0.0105776 &
298  /)
299  !!! MAIN CALCULATIONS
300  quan = (acos(-1.)/180.)*ls*gen1 + (acos(-1.)/180.)*lon*gen2
301  c = cos(quan)
302  s = sin(quan)
303  line(:) = 0.
304  DO year=1,2
305    line(year) = DOT_PRODUCT(c,tes(:,1,year)) - DOT_PRODUCT(s,tes(:,2,year)) + const(year)
306    !!! this is taken into account in isitice
307    !if (line(year) < -90.) line(year) = -90.
308    !if (line(year) >  90.) line(year) =  90.
309  ENDDO
310  !!! THIS IS TEMPORARY, BUT WHY NOT FOR MODELING
311  splatcrocus = (line(1)+line(2))/2.
312END FUNCTION splatcrocus
313
314END MODULE geticecover_mod
315
316!PROGRAM main
317!
318!implicit none
319!INTEGER, PARAMETER :: ngrid = 12
320!REAL,DIMENSION(ngrid) :: lontab,lattab,icecover
321!REAL :: ls,lo,outputs,isitice
322!REAL :: nplatcrocus,splatcrocus
323!INTEGER :: i,j,zels
324!
325!DO i=0,360,30
326!print *, 'ls,np,sp ', float(i), nplatcrocus(float(i),0.), splatcrocus(float(i),0.)
327!ENDDO
328!
329!print *, 'isitice(90.,0.,50.)',isitice(90.,0.,50.)
330!print *, 'isitice(90.,0.,87.)',isitice(90.,0.,87.)
331!print *, 'isitice(90.,0.,-50.)',isitice(90.,0.,-50.)
332!print *, 'isitice(90.,0.,-87.)',isitice(90.,0.,-87.)
333!print *, 'isitice(290.,0.,50.)',isitice(290.,0.,50.)
334!print *, 'isitice(290.,0.,87.)',isitice(290.,0.,87.)
335!print *, 'isitice(290.,0.,-50.)',isitice(290.,0.,-50.)
336!print *, 'isitice(290.,0.,-87.)',isitice(290.,0.,-87.)
337!print *, 'isitice(120.,0.,90.)',isitice(120.,0.,90.)
338!
339!lontab(:) = 0.
340!lattab = (/ -90., -75., -60., -45., -30., -15.,  15.,  30.,  45.,  60.,  75.,  90. /)
341!print *,'ls. lat:', lattab
342!DO zels=0,360,30
343!  CALL geticecover( ngrid, float(zels), lontab, lattab, icecover )
344!  print *, zels, icecover
345!ENDDO
346!
347!END PROGRAM main
Note: See TracBrowser for help on using the repository browser.