source: trunk/LMDZ.MARS/libf/phymars/albedocaps.F90 @ 3026

Last change on this file since 3026 was 2999, checked in by llange, 18 months ago

Mars PCM
Include perenial_co2ice (equivalent of watercap) to distinguich between CO2 frost and perenial CO2 ice for paleoclimate studies.
When no frost is present and we dig into perenial ice, the surface albedo is changed. The albedo for seasonal ice is set to 0.65, and the perenial ice albedo can be fixed in the callphys.def. I recommand values between 0.8 and 0.9.
To use this, paleoclimate must be set to True and TESalbedo to false in the callphys.def. Else, it runs as usual with TES albedo
LL
File size: 17.7 KB
Line 
1subroutine albedocaps(zls,ngrid,piceco2,piceco2_peren,psolaralb,emisref)
2
3! routine which changes the albedo (and emissivity) of the surface
4! depending on the presence of CO2 ice on the surface
5
6use ioipsl_getin_p_mod, only: getin_p
7use geometry_mod, only: latitude ! grid point latitudes (rad)
8use surfdat_h, only: TESicealbedo, TESice_Ncoef, TESice_Scoef, &
9                     emisice, albedice, watercaptag, albedo_h2o_cap, &
10                     emissiv, albedodat
11USE mod_phys_lmdz_transfert_para, ONLY: bcast
12USE mod_phys_lmdz_para, ONLY: is_master
13USE paleoclimate_mod, ONLY: paleoclimate,albedo_perenialco2
14
15implicit none
16
17include"callkeys.h"
18
19! arguments:
20real,intent(in) :: zls ! solar longitude (rad)
21integer,intent(in) :: ngrid
22real,intent(in) :: piceco2(ngrid) ! amount of CO2 ice on the surface (kg/m2)
23real,intent(inout) :: piceco2_peren(ngrid) ! amount of perenial co2 ice (kg/m^2)
24real,intent(out) :: psolaralb(ngrid,2) ! albedo of the surface
25real,intent(out) :: emisref(ngrid) ! emissivity of the surface
26
27
28! local variables:
29logical,save :: firstcall=.true.
30!$OMP THREADPRIVATE(firstcall)
31integer :: ig,icap
32
33!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34
35! local variables:
36real,save :: zls_old ! value of zls from a previous call
37real,save :: pi,radeg ! to convert radians to degrees
38
39!$OMP THREADPRIVATE(zls_old,pi,radeg)
40
41! TES datasets: (hard coded fixed length/sizes; for now)
42integer,parameter :: TESlonsize=72
43! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
44real,save :: TESlon(TESlonsize)
45integer,parameter :: TESlatsize=30
46! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
47! to TESlatn(30)=89 ; TESlatn(8)=45
48real,save :: TESlatn(TESlatsize)
49! latitudes (south hemisphere file), in degrees, from TESlats(1)=-89,
50! to TESlats(30)=-31 ; TESlats(23)=-45
51real,save :: TESlats(TESlatsize)
52integer,parameter :: TESlssize=72
53! Solar longitude in TES files, TESls(1)=2.5 to TESls(72)=357.5
54real,save :: TESls(TESlssize)
55! TES North albedo (=-1 for missing values)
56real,save :: TESalbn(TESlonsize,TESlatsize,TESlssize)
57! TES South albedo (=-1 for missing values)
58real,save :: TESalbs(TESlonsize,TESlatsize,TESlssize)
59
60!$OMP THREADPRIVATE(TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
61
62
63!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
64
65! 1. Initializations
66! AS: OK firstcall absolute
67if (firstcall) then
68  ! find out if user wants to use TES cap albedoes or not
69  TESicealbedo=.false. ! default value
70  write(*,*)" albedocaps: Use TES Cap albedoes ?"
71  call getin_p("TESicealbedo",TESicealbedo)
72  write(*,*)" albedocaps: TESicealbedo = ",TESicealbedo
73 
74  ! if using TES albedoes, load coeffcients
75  if (TESicealbedo) then
76    write(*,*)" albedocaps: Coefficient for Northern Cap ?"
77    TESice_Ncoef=1.0 ! default value
78    call getin_p("TESice_Ncoef",TESice_Ncoef)
79    write(*,*)" albedocaps: TESice_Ncoef = ",TESice_Ncoef
80   
81    write(*,*)" albedocaps: Coefficient for Southern Cap ?"
82    TESice_Scoef=1.0 ! default value
83    call getin_p("TESice_Scoef",TESice_Scoef)
84    write(*,*)" albedocaps: TESice_Scoef = ",TESice_Scoef
85  endif
86 
87  call read_TES_icecap_albedo(zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
88 
89  firstcall=.false.
90endif ! of if (firstcall)
91   
92do ig=1,ngrid
93  if (latitude(ig).lt.0.) then
94    icap=2 ! Southern hemisphere
95  else
96    icap=1 ! Northern hemisphere
97  endif
98 
99  if (piceco2(ig).gt.0) then
100    ! set emissivity of surface to be the ice emissivity
101    emisref(ig)=emisice(icap)
102    ! set the surface albedo to be the ice albedo
103    if (TESicealbedo) then
104      call TES_icecap_albedo(zls,ig,psolaralb(ig,1),icap,zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
105      psolaralb(ig,2)=psolaralb(ig,1)
106    else
107      psolaralb(ig,1)=albedice(icap)
108      psolaralb(ig,2)=albedice(icap)
109      if(paleoclimate) then
110         if((piceco2_peren(ig).gt.0.).and.(piceco2(ig).lt.piceco2_peren(ig))) then
111             psolaralb(ig,1) = albedo_perenialco2
112             psolaralb(ig,2) = albedo_perenialco2
113             piceco2_peren(ig) = piceco2(ig)
114         endif
115      endif
116    endif
117  else if (watercaptag(ig) .and. water) then
118    ! there is a water ice cap: set the surface albedo to the water ice one
119    ! to do : emissivity
120    emisref(ig) = 1
121    psolaralb(ig,1)=albedo_h2o_cap
122    psolaralb(ig,2)=albedo_h2o_cap
123  else
124    ! set emissivity of surface to be bare ground emissivity
125    emisref(ig)=emissiv
126    ! set the surface albedo to bare ground albedo
127    psolaralb(ig,1)=albedodat(ig)
128    psolaralb(ig,2)=albedodat(ig)
129  endif ! of if (piceco2(ig).gt.0)
130enddo ! of ig=1,ngrid
131
132end subroutine albedocaps
133
134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135
136subroutine read_TES_icecap_albedo(zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
137
138use datafile_mod, only: datadir
139use netcdf, only: nf90_open, NF90_NOWRITE, NF90_NOERR, &
140                  nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
141USE mod_phys_lmdz_para, ONLY: is_master
142USE mod_phys_lmdz_transfert_para, ONLY: bcast
143                 
144implicit none
145
146! arguments:
147
148! local variables:
149real:: zls_old ! value of zls from a previous call
150real:: pi,radeg ! to convert radians to degrees
151character(len=20),parameter :: modname="TES_icecap_albedo"
152
153! TES datasets: (hard coded fixed length/sizes; for now)
154integer,parameter:: TESlonsize=72
155! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
156real:: TESlon(TESlonsize)
157integer,parameter:: TESlatsize=30
158! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
159! to TESlatn(30)=89 ; TESlatn(8)=45
160real:: TESlatn(TESlatsize)
161! latitudes (south hemisphere file), in degrees, from TESlats(1)=-89,
162! to TESlats(30)=-31 ; TESlats(23)=-45
163real:: TESlats(TESlatsize)
164integer,parameter:: TESlssize=72
165! Solar longitude in TES files, TESls(1)=2.5 to TESls(72)=357.5
166real:: TESls(TESlssize)
167! TES North albedo (=-1 for missing values)
168real:: TESalbn(TESlonsize,TESlatsize,TESlssize)
169! TES South albedo (=-1 for missing values)
170real:: TESalbs(TESlonsize,TESlatsize,TESlssize)
171
172!NetCDF variables:
173integer :: ierr ! NetCDF status
174integer :: nid ! NetCDF file ID
175integer :: nvarid ! NetCDF variable ID
176
177! 0. Preliminary stuff
178
179if(is_master) then
180
181! Load TES albedoes for Northern Hemisphere
182  ierr=nf90_open(trim(datadir)//"/npsc_albedo.nc",NF90_NOWRITE,nid)
183  IF (ierr.NE.NF90_NOERR) THEN
184    write(*,*)'Problem opening npsc_albedo.nc (phymars/albedocaps.F90)'
185    write(*,*)'It should be in :',trim(datadir),'/'
186    write(*,*)'1) You can change this directory address in callfis.def with'
187    write(*,*)'   datadir=/path/to/datafiles'
188    write(*,*)'2) If necessary, npsc_albedo.nc (and other datafiles)'
189    write(*,*)'   can be obtained online on:'
190    write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
191    CALL abort_physic(modname,"missing input file",1)
192  ELSE
193    write(*,*) "albedocaps: using file ",trim(datadir)//"/npsc_albedo.nc"
194  ENDIF
195 
196  ierr=nf90_inq_varid(nid,"longitude",nvarid)
197  if (ierr.ne.NF90_NOERR) then
198    write(*,*) "Failed to find longitude in file!"
199    write(*,*)trim(nf90_strerror(ierr))
200    call abort_physic(modname,"failed finding longitude",1)
201  else
202    ierr=nf90_get_var(nid,nvarid,TESlon)
203    if (ierr.ne.NF90_NOERR) then
204      write(*,*) "Failed loading longitude data from file!"
205      write(*,*)trim(nf90_strerror(ierr))
206      call abort_physic(modname,"failed loading longitude",1)
207    endif
208  endif
209 
210  ierr=nf90_inq_varid(nid,"latitude",nvarid)
211  if (ierr.ne.NF90_NOERR) then
212    write(*,*) "Failed to find latitude in file!"
213    write(*,*)trim(nf90_strerror(ierr))
214    call abort_physic(modname,"failed finding latitude",1)
215  else
216    ierr=nf90_get_var(nid,nvarid,TESlatn)
217    if (ierr.ne.NF90_NOERR) then
218      write(*,*) "Failed loading latitude data from file!"
219      write(*,*)trim(nf90_strerror(ierr))
220      call abort_physic(modname,"failed loading latitude",1)
221    endif
222  endif
223 
224  ierr=nf90_inq_varid(nid,"time",nvarid)
225  if (ierr.ne.NF90_NOERR) then
226    write(*,*) "Failed to find time in file!"
227    write(*,*)trim(nf90_strerror(ierr))
228    call abort_physic(modname,"failed finding time",1)
229  else
230    ierr=nf90_get_var(nid,nvarid,TESls)
231    if (ierr.ne.NF90_NOERR) then
232      write(*,*) "Failed loading time data from file!"
233      write(*,*)trim(nf90_strerror(ierr))
234      call abort_physic(modname,"failed loading time",1)
235    endif
236  endif
237 
238  ierr=nf90_inq_varid(nid,"albedo",nvarid)
239  if (ierr.ne.NF90_NOERR) then
240    write(*,*) "Failed to find albedo in file!"
241    write(*,*)trim(nf90_strerror(ierr))
242    call abort_physic(modname,"failed finding albedo",1)
243  else
244    ierr=nf90_get_var(nid,nvarid,TESalbn)
245    if (ierr.ne.NF90_NOERR) then
246      write(*,*) "Failed loading albedo data from file!"
247      write(*,*)trim(nf90_strerror(ierr))
248      call abort_physic(modname,"failed loading albedo",1)
249    endif
250  endif
251
252  ierr=nf90_close(nid)
253
254! Load albedoes for Southern Hemisphere
255  ierr=nf90_open(trim(datadir)//"/spsc_albedo.nc",NF90_NOWRITE,nid)
256  IF (ierr.NE.NF90_NOERR) THEN
257    write(*,*)'Problem opening spsc_albedo.nc (phymars/albedocaps.F90)'
258    write(*,*)'It should be in :',trim(datadir),'/'
259    write(*,*)'1) You can change this directory address in callfis.def with'
260    write(*,*)'   datadir=/path/to/datafiles'
261    write(*,*)'2) If necessary, spsc_albedo.nc (and other datafiles)'
262    write(*,*)'   can be obtained online on:'
263    write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
264    CALL abort_physic(modname,"missing input file",1)
265  ELSE
266    write(*,*) "albedocaps: using file ",trim(datadir)//"/spsc_albedo.nc"
267  ENDIF
268
269  ierr=nf90_inq_varid(nid,"latitude",nvarid)
270  if (ierr.ne.NF90_NOERR) then
271    write(*,*) "Failed to find latitude in file!"
272    write(*,*)trim(nf90_strerror(ierr))
273    call abort_physic(modname,"failed finding latitude",1)
274  else
275    ierr=nf90_get_var(nid,nvarid,TESlats)
276    if (ierr.ne.NF90_NOERR) then
277      write(*,*) "Failed loading latitude data from file!"
278      write(*,*)trim(nf90_strerror(ierr))
279      call abort_physic(modname,"failed loading latitude",1)
280    endif
281  endif
282
283  ierr=nf90_inq_varid(nid,"albedo",nvarid)
284  if (ierr.ne.NF90_NOERR) then
285    write(*,*) "Failed to find albedo in file!"
286    write(*,*)trim(nf90_strerror(ierr))
287    call abort_physic(modname,"failed finding albedo",1)
288  else
289    ierr=nf90_get_var(nid,nvarid,TESalbs)
290    if (ierr.ne.NF90_NOERR) then
291      write(*,*) "Failed loading albedo data from file!"
292      write(*,*)trim(nf90_strerror(ierr))
293      call abort_physic(modname,"failed loading albedo",1)
294    endif
295  endif
296
297  ierr=nf90_close(nid)
298
299! constants:
300  pi=acos(-1.)
301  radeg=180/pi
302
303  zls_old=-999 ! dummy initialization
304   
305endif !is_master
306   
307    call bcast(TESlon)
308    call bcast(TESlatn)
309    call bcast(TESlats)
310    call bcast(TESls)
311    call bcast(TESalbn)
312    call bcast(TESalbs)
313    call bcast(pi)
314    call bcast(zls_old)
315    call bcast(radeg)
316 
317end subroutine read_TES_icecap_albedo
318
319subroutine TES_icecap_albedo(zls,ig,alb,icap,zls_old,pi,radeg,TESlon,TESlatn,TESlats,TESls,TESalbn,TESalbs)
320
321use geometry_mod, only: latitude, longitude ! in radians
322use surfdat_h, only: albedice, TESice_Ncoef, TESice_Scoef
323use datafile_mod, only: datadir
324use netcdf, only: nf90_open, NF90_NOWRITE, NF90_NOERR, &
325                  nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
326USE mod_phys_lmdz_transfert_para, ONLY: bcast
327                 
328implicit none
329
330! arguments:
331real,intent(in) :: zls ! solar longitude (rad)
332integer,intent(in) :: ig ! grid point index
333real,intent(out) :: alb ! (interpolated) TES ice albedo at that grid point
334integer :: icap ! =1: Northern hemisphere =2: Southern hemisphere
335
336! local variables:
337integer,save :: tinf,tsup ! encompassing time indexes of TES data
338real,save :: reltime ! relative position in-between time indexes (in [0;1])
339integer :: latinf,latsup ! encompassing latitude indexes of TES data
340real :: rellat ! relative position in-between latitude indexes (in[0;1])
341integer :: loninf,lonsup ! encompassing longitude indexes of TES data
342real :: rellon !relative position in-between longitude indexes (in[0;1])
343real :: zlsd ! solar longitude, in degrees
344real :: latd ! latitude, in degrees
345real :: lond ! longitude, in degrees
346integer :: i
347
348!$OMP THREADPRIVATE(tinf,tsup,reltime)
349
350! TES datasets: (hard coded fixed length/sizes; for now)
351real,parameter :: TESdeltalon=5.0 ! step in longitude in TES files
352real,parameter :: TESdeltalat=2.0 ! step in latitude in TES files
353real,parameter :: TESlatnmin=45. ! minimum TES latitude (North hemisphere)
354real,parameter :: TESlatsmax=-45. ! maximum TES latitude (South hemisphere)
355! to TESlats(30)=-31 ; TESlats(23)=-45
356real,parameter :: TESdeltals=5.0 ! step in solar longitude in TES files
357
358real:: zls_old ! value of zls from a previous call
359real:: pi,radeg ! to convert radians to degrees
360
361! TES datasets: (hard coded fixed length/sizes; for now)
362integer,parameter:: TESlonsize=72
363! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
364real:: TESlon(TESlonsize)
365integer,parameter:: TESlatsize=30
366! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
367! to TESlatn(30)=89 ; TESlatn(8)=45
368real:: TESlatn(TESlatsize)
369! latitudes (south hemisphere file), in degrees, from TESlats(1)=-89,
370! to TESlats(30)=-31 ; TESlats(23)=-45
371real:: TESlats(TESlatsize)
372integer,parameter:: TESlssize=72
373! Solar longitude in TES files, TESls(1)=2.5 to TESls(72)=357.5
374real:: TESls(TESlssize)
375! TES North albedo (=-1 for missing values)
376real:: TESalbn(TESlonsize,TESlatsize,TESlssize)
377! TES South albedo (=-1 for missing values)
378real:: TESalbs(TESlonsize,TESlatsize,TESlssize)
379! encompassing nodes arranged as follow : 4 3
380real :: val(4)                          ! 1 2
381
382! 1. Identify encompassing latitudes
383
384! Check that latitude is such that there is TES data to use
385! (ie: latitude 45 deg and poleward) otherwise use 'default' albedoes
386
387latd=latitude(ig)*radeg ! latitude, in degrees
388if (icap.eq.1) then
389 ! North hemisphere
390 if (latd.lt.TESlatnmin) then
391   alb=albedice(1)
392   ! the job is done; quit this routine
393   return
394 else
395   ! find encompassing latitudes
396   if (latd.ge.TESlatn(TESlatsize)) then
397     latinf=TESlatsize
398     latsup=TESlatsize
399     rellat=0.
400   else
401     do i=1,TESlatsize-1
402       if ((latd.ge.TESlatn(i)).and.(latd.lt.TESlatn(i+1))) then
403         latinf=i
404         latsup=i+1
405         rellat=(latd-TESlatn(i))/TESdeltalat
406         exit ! found encompassing indexes; quit loop
407       endif
408     enddo
409   endif
410 endif ! of if (latd.lt.TESlatnmin)
411else ! icap=2
412 ! South hemisphere
413 if (latd.gt.TESlatsmax) then
414   alb=albedice(2)
415   ! the job is done; quit this routine
416   return
417 else
418   ! find encompassing latitudes
419   if (latd.lt.TESlats(1)) then
420     latinf=1
421     latsup=1
422     rellat=0.
423   else
424     do i=1,TESlatsize-1
425       if ((latd.ge.TESlats(i)).and.(latd.lt.TESlats(i+1))) then
426         latinf=i
427         latsup=i+1
428         rellat=(latd-TESlats(i))/TESdeltalat
429         exit ! found encompassing indexes; quit loop
430       endif
431     enddo
432   endif
433 endif ! of if (latd.gt.-45.)
434endif ! of if (icap.eq.1)
435
436! 2. Identify encompassing time indexes
437if (zls.ne.zls_old) then
438  zlsd=zls*radeg ! solar longitude, in degrees
439 
440  if (zlsd.lt.TESls(1)) then
441    tinf=TESlssize
442    tsup=1
443    reltime=0.5+zlsd/TESdeltals
444  else
445    if (zlsd.ge.TESls(TESlssize)) then
446      tinf=TESlssize
447      tsup=1
448      reltime=(360.-zlsd)/TESdeltals
449    else
450      ! look for encompassing indexes
451      do i=1,TESlssize-1
452        if ((zlsd.ge.TESls(i)).and.(zlsd.lt.TESls(i+1))) then
453          tinf=i
454          tsup=i+1
455          reltime=(zlsd-TESls(i))/TESdeltals
456          exit ! quit loop, we found the indexes
457        endif
458      enddo
459    endif
460  endif ! of if (zlsd.lt.TESls(1))
461 
462  zls_old=zls ! store current zls
463endif ! of if (zls.ne.zls_old)
464
465! 3. Identify encompassing longitudes
466lond=longitude(ig)*radeg ! east longitude, in degrees
467if (lond.lt.TESlon(1)) then
468  loninf=TESlonsize
469  lonsup=1
470  rellon=0.5+(180.+lond)/TESdeltalon
471else
472  if (lond.ge.TESlon(TESlonsize)) then
473    loninf=TESlonsize
474    lonsup=1
475    rellon=(180-lond)/TESdeltalon
476  else
477    do i=1,TESlonsize-1
478      if ((lond.ge.TESlon(i)).and.(lond.lt.TESlon(i+1))) then
479        loninf=i
480        lonsup=i+1
481        rellon=(lond-TESlon(i))/TESdeltalon
482        exit ! quit loop, we found the indexes
483      endif
484    enddo
485  endif ! of if (lond.ge.TESlon(TESlonsize))
486endif ! of if (lond.lt.TESlon(1))
487
488! 4. Use linear interpolation in time to build encompassing nodal values
489!    encompassing nodes are arranged as follow : 4 3
490!                                                1 2
491if (icap.eq.1) then
492  ! Northern hemisphere
493  val(1)=(1.-reltime)*TESalbn(loninf,latinf,tinf) &
494         +reltime*TESalbn(loninf,latinf,tsup)
495  val(2)=(1.-reltime)*TESalbn(lonsup,latinf,tinf) &
496         +reltime*TESalbn(lonsup,latinf,tsup)
497  val(3)=(1.-reltime)*TESalbn(lonsup,latsup,tinf) &
498         +reltime*TESalbn(lonsup,latsup,tsup)
499  val(4)=(1.-reltime)*TESalbn(loninf,latsup,tinf) &
500         +reltime*TESalbn(loninf,latsup,tsup)
501else
502  ! Southern hemisphere
503  val(1)=(1.-reltime)*TESalbs(loninf,latinf,tinf) &
504         +reltime*TESalbs(loninf,latinf,tsup)
505  val(2)=(1.-reltime)*TESalbs(lonsup,latinf,tinf) &
506         +reltime*TESalbs(lonsup,latinf,tsup)
507  val(3)=(1.-reltime)*TESalbs(lonsup,latsup,tinf) &
508         +reltime*TESalbs(lonsup,latsup,tsup)
509  val(4)=(1.-reltime)*TESalbs(loninf,latsup,tinf) &
510         +reltime*TESalbs(loninf,latsup,tsup)
511endif ! of if (icap.eq.1)
512
513! 5. Use bilinear interpolation to compute albedo
514alb=(1.-rellon)*(1.-rellat)*val(1) &
515    +rellon*(1.-rellat)*val(2) &
516    +rellon*rellat*val(3) &
517    +(1.-rellon)*rellat*val(4)
518
519! 6. Apply coefficient to interpolated TES albedo
520if (icap.eq.1) then
521  alb=alb*TESice_Ncoef
522else
523  alb=alb*TESice_Scoef
524endif ! of if (icap.eq.1)
525
526! Make sure that returned albedo is never greater than 0.90
527if (alb.gt.0.90) alb=0.90
528
529end subroutine TES_icecap_albedo
530
Note: See TracBrowser for help on using the repository browser.