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

Last change on this file since 1268 was 1268, checked in by aslmd, 11 years ago

LMDZ.MARS. related to r1266. forgot to remove a few now-obsolete dimensions.h includes in Mars physics.
File size: 13.4 KB
Line 
1subroutine albedocaps(zls,ngrid,piceco2,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
6! to use the 'getin' routine
7use ioipsl_getincom, only: getin
8use comgeomfi_h, only: lati ! grid point latitudes (rad)
9use surfdat_h, only: TESicealbedo, TESice_Ncoef, TESice_Scoef, &
10                     emisice, albedice, watercaptag, albedo_h2o_ice, &
11                     emissiv, albedodat
12implicit none
13
14#include"callkeys.h"
15
16! arguments:
17real,intent(in) :: zls ! solar longitude (rad)
18integer,intent(in) :: ngrid
19real,intent(in) :: piceco2(ngrid) ! amount of CO2 ice on the surface (kg/m2)
20real,intent(out) :: psolaralb(ngrid,2) ! albedo of the surface
21real,intent(out) :: emisref(ngrid) ! emissivity of the surface
22
23
24! local variables:
25logical,save :: firstcall=.true.
26integer :: ig,icap
27
28! 1. Initializations
29if (firstcall) then
30  ! find out if user wants to use TES cap albedoes or not
31  TESicealbedo=.false. ! default value
32  write(*,*)" albedocaps: Use TES Cap albedoes ?"
33  call getin("TESicealbedo",TESicealbedo)
34  write(*,*)" albedocaps: TESicealbedo = ",TESicealbedo
35 
36  ! if using TES albedoes, load coeffcients
37  if (TESicealbedo) then
38    write(*,*)" albedocaps: Coefficient for Northern Cap ?"
39    TESice_Ncoef=1.0 ! default value
40    call getin("TESice_Ncoef",TESice_Ncoef)
41    write(*,*)" albedocaps: TESice_Ncoef = ",TESice_Ncoef
42   
43    write(*,*)" albedocaps: Coefficient for Southern Cap ?"
44    TESice_Scoef=1.0 ! default value
45    call getin("TESice_Scoef",TESice_Scoef)
46    write(*,*)" albedocaps: TESice_Scoef = ",TESice_Scoef
47  endif
48 
49  firstcall=.false.
50endif ! of if (firstcall)
51
52do ig=1,ngrid
53  if (lati(ig).lt.0.) then
54    icap=2 ! Southern hemisphere
55  else
56    icap=1 ! Northern hemisphere
57  endif
58
59  if (piceco2(ig).gt.0) then
60    ! set emissivity of surface to be the ice emissivity
61    emisref(ig)=emisice(icap)
62    ! set the surface albedo to be the ice albedo
63    if (TESicealbedo) then
64      call TES_icecap_albedo(zls,ig,psolaralb(ig,1),icap)
65      psolaralb(ig,2)=psolaralb(ig,1)
66    else
67      psolaralb(ig,1)=albedice(icap)
68      psolaralb(ig,2)=albedice(icap)
69    endif
70  else if (watercaptag(ig) .and. water) then
71    ! there is a water ice cap: set the surface albedo to the water ice one
72    ! to do : emissivity
73    emisref(ig) = 1
74    psolaralb(ig,1)=albedo_h2o_ice
75    psolaralb(ig,2)=albedo_h2o_ice
76  else
77    ! set emissivity of surface to be bare ground emissivity
78    emisref(ig)=emissiv
79    ! set the surface albedo to bare ground albedo
80    psolaralb(ig,1)=albedodat(ig)
81    psolaralb(ig,2)=albedodat(ig)
82  endif ! of if (piceco2(ig).gt.0)
83enddo ! of ig=1,ngrid
84end subroutine albedocaps
85
86!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87subroutine TES_icecap_albedo(zls,ig,alb,icap)
88
89use comgeomfi_h, only: lati, long
90use surfdat_h, only: albedice, TESice_Ncoef, TESice_Scoef
91use netcdf, only: nf90_open, NF90_NOWRITE, NF90_NOERR, &
92                  nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
93                 
94implicit none
95#include"dimensions.h"
96#include"datafile.h"
97
98! arguments:
99real,intent(in) :: zls ! solar longitude (rad)
100integer,intent(in) :: ig ! grid point index
101real,intent(out) :: alb ! (interpolated) TES ice albedo at that grid point
102integer :: icap ! =1: Northern hemisphere =2: Southern hemisphere
103
104! local variables:
105logical,save :: firstcall=.true.
106real,save :: zls_old ! value of zls from a previous call
107integer,save :: tinf,tsup ! encompassing time indexes of TES data
108real,save :: reltime ! relative position in-between time indexes (in [0;1])
109integer :: latinf,latsup ! encompassing latitude indexes of TES data
110real :: rellat ! relative position in-between latitude indexes (in[0;1])
111integer :: loninf,lonsup ! encompassing longitude indexes of TES data
112real :: rellon !relative position in-between longitude indexes (in[0;1])
113real,save :: pi,radeg ! to convert radians to degrees
114real :: zlsd ! solar longitude, in degrees
115real :: latd ! latitude, in degrees
116real :: lond ! longitude, in degrees
117integer :: i
118
119! TES datasets: (hard coded fixed length/sizes; for now)
120integer,parameter :: TESlonsize=72
121real,parameter :: TESdeltalon=5.0 ! step in longitude in TES files
122! longitudes, in TES files, in degrees, from TESlon(1)=-177.5 to TESlon(72)=177.5
123real,save :: TESlon(TESlonsize)
124integer,parameter :: TESlatsize=30
125real,parameter :: TESdeltalat=2.0 ! step in latitude in TES files
126! latitudes (north hemisphere file), in degrees, from TESlatn(1)=31,
127! to TESlatn(30)=89 ; TESlatn(8)=45
128real,parameter :: TESlatnmin=45. ! minimum TES latitude (North hemisphere)
129real,parameter :: TESlatsmax=-45. ! maximum TES latitude (South hemisphere)
130real,save :: TESlatn(TESlatsize)
131! latitudes (south hemisphere file), in degrees, from TESlats(1)=-89,
132! to TESlats(30)=-31 ; TESlats(23)=-45
133real,save :: TESlats(TESlatsize)
134integer,parameter :: TESlssize=72
135real,parameter :: TESdeltals=5.0 ! step in solar longitude in TES files
136! Solar longitude in TES files, TESls(1)=2.5 to TESls(72)=357.5
137real,save :: TESls(TESlssize)
138! TES North albedo (=-1 for missing values)
139real,save :: TESalbn(TESlonsize,TESlatsize,TESlssize)
140! TES South albedo (=-1 for missing values)
141real,save :: TESalbs(TESlonsize,TESlatsize,TESlssize)
142! encompassing nodes arranged as follow : 4 3
143real :: val(4)                          ! 1 2
144
145!NetCDF variables:
146integer :: ierr ! NetCDF status
147integer :: nid ! NetCDF file ID
148integer :: nvarid ! NetCDF variable ID
149
150! 0. Preliminary stuff
151if (firstcall) then
152! Load TES albedoes for Northern Hemisphere
153  ! Note: datafile() is defined in "datafile.h"
154  ierr=nf90_open(trim(datafile)//"/npsc_albedo.nc",NF90_NOWRITE,nid)
155  IF (ierr.NE.NF90_NOERR) THEN
156    write(*,*)'Problem opening npsc_albedo.nc (phymars/albedocaps.F90)'
157    write(*,*)'It should be in :',trim(datafile),'/'
158    write(*,*)'1) You can change this directory address in callfis.def with'
159    write(*,*)'   datadir=/path/to/datafiles'
160    write(*,*)'2) If necessary, npsc_albedo.nc (and other datafiles)'
161    write(*,*)'   can be obtained online on:'
162    write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
163    CALL ABORT
164  ELSE
165    write(*,*) "albedocaps: using file ",trim(datafile)//"/npsc_albedo.nc"
166  ENDIF
167 
168  ierr=nf90_inq_varid(nid,"longitude",nvarid)
169  if (ierr.ne.NF90_NOERR) then
170    write(*,*) "Failed to find longitude in file!"
171    write(*,*)trim(nf90_strerror(ierr))
172    stop
173  else
174    ierr=nf90_get_var(nid,nvarid,TESlon)
175    if (ierr.ne.NF90_NOERR) then
176      write(*,*) "Failed loading longitude data from file!"
177      write(*,*)trim(nf90_strerror(ierr))
178      stop
179    endif
180  endif
181 
182  ierr=nf90_inq_varid(nid,"latitude",nvarid)
183  if (ierr.ne.NF90_NOERR) then
184    write(*,*) "Failed to find latitude in file!"
185    write(*,*)trim(nf90_strerror(ierr))
186    stop
187  else
188    ierr=nf90_get_var(nid,nvarid,TESlatn)
189    if (ierr.ne.NF90_NOERR) then
190      write(*,*) "Failed loading latitude data from file!"
191      write(*,*)trim(nf90_strerror(ierr))
192      stop
193    endif
194  endif
195 
196  ierr=nf90_inq_varid(nid,"time",nvarid)
197  if (ierr.ne.NF90_NOERR) then
198    write(*,*) "Failed to find time in file!"
199    write(*,*)trim(nf90_strerror(ierr))
200    stop
201  else
202    ierr=nf90_get_var(nid,nvarid,TESls)
203    if (ierr.ne.NF90_NOERR) then
204      write(*,*) "Failed loading time data from file!"
205      write(*,*)trim(nf90_strerror(ierr))
206      stop
207    endif
208  endif
209 
210  ierr=nf90_inq_varid(nid,"albedo",nvarid)
211  if (ierr.ne.NF90_NOERR) then
212    write(*,*) "Failed to find albedo in file!"
213    write(*,*)trim(nf90_strerror(ierr))
214    stop
215  else
216    ierr=nf90_get_var(nid,nvarid,TESalbn)
217    if (ierr.ne.NF90_NOERR) then
218      write(*,*) "Failed loading albedo data from file!"
219      write(*,*)trim(nf90_strerror(ierr))
220      stop
221    endif
222  endif
223
224  ierr=nf90_close(nid)
225
226! Load albedoes for Southern Hemisphere
227  ierr=nf90_open(trim(datafile)//"/spsc_albedo.nc",NF90_NOWRITE,nid)
228  IF (ierr.NE.NF90_NOERR) THEN
229    write(*,*)'Problem opening spsc_albedo.nc (phymars/albedocaps.F90)'
230    write(*,*)'It should be in :',trim(datafile),'/'
231    write(*,*)'1) You can change this directory address in callfis.def with'
232    write(*,*)'   datadir=/path/to/datafiles'
233    write(*,*)'2) If necessary, spsc_albedo.nc (and other datafiles)'
234    write(*,*)'   can be obtained online on:'
235    write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
236    CALL ABORT
237  ELSE
238    write(*,*) "albedocaps: using file ",trim(datafile)//"/spsc_albedo.nc"
239  ENDIF
240
241  ierr=nf90_inq_varid(nid,"latitude",nvarid)
242  if (ierr.ne.NF90_NOERR) then
243    write(*,*) "Failed to find latitude in file!"
244    write(*,*)trim(nf90_strerror(ierr))
245    stop
246  else
247    ierr=nf90_get_var(nid,nvarid,TESlats)
248    if (ierr.ne.NF90_NOERR) then
249      write(*,*) "Failed loading latitude data from file!"
250      write(*,*)trim(nf90_strerror(ierr))
251      stop
252    endif
253  endif
254
255  ierr=nf90_inq_varid(nid,"albedo",nvarid)
256  if (ierr.ne.NF90_NOERR) then
257    write(*,*) "Failed to find albedo in file!"
258    write(*,*)trim(nf90_strerror(ierr))
259    stop
260  else
261    ierr=nf90_get_var(nid,nvarid,TESalbs)
262    if (ierr.ne.NF90_NOERR) then
263      write(*,*) "Failed loading albedo data from file!"
264      write(*,*)trim(nf90_strerror(ierr))
265      stop
266    endif
267  endif
268
269  ierr=nf90_close(nid)
270
271! constants:
272  pi=acos(-1.)
273  radeg=180/pi
274
275  zls_old=-999 ! dummy initialization
276
277  firstcall=.false.
278endif ! of if firstcall
279
280! 1. Identify encompassing latitudes
281
282! Check that latitude is such that there is TES data to use
283! (ie: latitude 45 deg and poleward) otherwise use 'default' albedoes
284latd=lati(ig)*radeg ! latitude, in degrees
285if (icap.eq.1) then
286 ! North hemisphere
287 if (latd.lt.TESlatnmin) then
288   alb=albedice(1)
289   ! the job is done; quit this routine
290   return
291 else
292   ! find encompassing latitudes
293   if (latd.ge.TESlatn(TESlatsize)) then
294     latinf=TESlatsize
295     latsup=TESlatsize
296     rellat=0.
297   else
298     do i=1,TESlatsize-1
299       if ((latd.ge.TESlatn(i)).and.(latd.lt.TESlatn(i+1))) then
300         latinf=i
301         latsup=i+1
302         rellat=(latd-TESlatn(i))/TESdeltalat
303         exit ! found encompassing indexes; quit loop
304       endif
305     enddo
306   endif
307 endif ! of if (latd.lt.TESlatnmin)
308else ! icap=2
309 ! South hemisphere
310 if (latd.gt.TESlatsmax) then
311   alb=albedice(2)
312   ! the job is done; quit this routine
313   return
314 else
315   ! find encompassing latitudes
316   if (latd.lt.TESlats(1)) then
317     latinf=1
318     latsup=1
319     rellat=0.
320   else
321     do i=1,TESlatsize-1
322       if ((latd.ge.TESlats(i)).and.(latd.lt.TESlats(i+1))) then
323         latinf=i
324         latsup=i+1
325         rellat=(latd-TESlats(i))/TESdeltalat
326         exit ! found encompassing indexes; quit loop
327       endif
328     enddo
329   endif
330 endif ! of if (latd.gt.-45.)
331endif ! of if (icap.eq.1)
332
333! 2. Identify encompassing time indexes
334if (zls.ne.zls_old) then
335  zlsd=zls*radeg ! solar longitude, in degrees
336 
337  if (zlsd.lt.TESls(1)) then
338    tinf=TESlssize
339    tsup=1
340    reltime=0.5+zlsd/TESdeltals
341  else
342    if (zlsd.ge.TESls(TESlssize)) then
343      tinf=TESlssize
344      tsup=1
345      reltime=(360.-zlsd)/TESdeltals
346    else
347      ! look for encompassing indexes
348      do i=1,TESlssize-1
349        if ((zlsd.ge.TESls(i)).and.(zlsd.lt.TESls(i+1))) then
350          tinf=i
351          tsup=i+1
352          reltime=(zlsd-TESls(i))/TESdeltals
353          exit ! quit loop, we found the indexes
354        endif
355      enddo
356    endif
357  endif ! of if (zlsd.lt.TESls(1))
358 
359  zls_old=zls ! store current zls
360endif ! of if (zls.ne.zls_old)
361
362! 3. Identify encompassing longitudes
363lond=long(ig)*radeg ! east longitude, in degrees
364if (lond.lt.TESlon(1)) then
365  loninf=TESlonsize
366  lonsup=1
367  rellon=0.5+(180.+lond)/TESdeltalon
368else
369  if (lond.ge.TESlon(TESlonsize)) then
370    loninf=TESlonsize
371    lonsup=1
372    rellon=(180-lond)/TESdeltalon
373  else
374    do i=1,TESlonsize-1
375      if ((lond.ge.TESlon(i)).and.(lond.lt.TESlon(i+1))) then
376        loninf=i
377        lonsup=i+1
378        rellon=(lond-TESlon(i))/TESdeltalon
379        exit ! quit loop, we found the indexes
380      endif
381    enddo
382  endif ! of if (lond.ge.TESlon(TESlonsize))
383endif ! of if (lond.lt.TESlon(1))
384
385! 4. Use linear interpolation in time to build encompassing nodal values
386!    encompassing nodes are arranged as follow : 4 3
387!                                                1 2
388if (icap.eq.1) then
389  ! Northern hemisphere
390  val(1)=(1.-reltime)*TESalbn(loninf,latinf,tinf) &
391         +reltime*TESalbn(loninf,latinf,tsup)
392  val(2)=(1.-reltime)*TESalbn(lonsup,latinf,tinf) &
393         +reltime*TESalbn(lonsup,latinf,tsup)
394  val(3)=(1.-reltime)*TESalbn(lonsup,latsup,tinf) &
395         +reltime*TESalbn(lonsup,latsup,tsup)
396  val(4)=(1.-reltime)*TESalbn(loninf,latsup,tinf) &
397         +reltime*TESalbn(loninf,latsup,tsup)
398else
399  ! Southern hemisphere
400  val(1)=(1.-reltime)*TESalbs(loninf,latinf,tinf) &
401         +reltime*TESalbs(loninf,latinf,tsup)
402  val(2)=(1.-reltime)*TESalbs(lonsup,latinf,tinf) &
403         +reltime*TESalbs(lonsup,latinf,tsup)
404  val(3)=(1.-reltime)*TESalbs(lonsup,latsup,tinf) &
405         +reltime*TESalbs(lonsup,latsup,tsup)
406  val(4)=(1.-reltime)*TESalbs(loninf,latsup,tinf) &
407         +reltime*TESalbs(loninf,latsup,tsup)
408endif ! of if (icap.eq.1)
409
410! 5. Use bilinear interpolation to compute albedo
411alb=(1.-rellon)*(1.-rellat)*val(1) &
412    +rellon*(1.-rellat)*val(2) &
413    +rellon*rellat*val(3) &
414    +(1.-rellon)*rellat*val(4)
415
416! 6. Apply coefficient to interpolated TES albedo
417if (icap.eq.1) then
418  alb=alb*TESice_Ncoef
419else
420  alb=alb*TESice_Scoef
421endif ! of if (icap.eq.1)
422
423! Make sure that returned albedo is never greater than 0.90
424if (alb.gt.0.90) alb=0.90
425
426end subroutine TES_icecap_albedo
427
Note: See TracBrowser for help on using the repository browser.