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

Last change on this file since 2448 was 2304, checked in by emillour, 6 years ago

Mars GCM:
Some code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort"
EM
File size: 14.1 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
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_ice, &
10                     emissiv, albedodat
11implicit none
12
13include"callkeys.h"
14
15! arguments:
16real,intent(in) :: zls ! solar longitude (rad)
17integer,intent(in) :: ngrid
18real,intent(in) :: piceco2(ngrid) ! amount of CO2 ice on the surface (kg/m2)
19real,intent(out) :: psolaralb(ngrid,2) ! albedo of the surface
20real,intent(out) :: emisref(ngrid) ! emissivity of the surface
21
22
23! local variables:
24logical,save :: firstcall=.true.
25integer :: ig,icap
26
27! 1. Initializations
28! AS: OK firstcall absolute
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_p("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_p("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_p("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 (latitude(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 geometry_mod, only: latitude, longitude ! in radians
90use surfdat_h, only: albedice, TESice_Ncoef, TESice_Scoef
91use datafile_mod, only: datadir
92use netcdf, only: nf90_open, NF90_NOWRITE, NF90_NOERR, &
93                  nf90_strerror, nf90_inq_varid, nf90_get_var, nf90_close
94                 
95implicit none
96
97! arguments:
98real,intent(in) :: zls ! solar longitude (rad)
99integer,intent(in) :: ig ! grid point index
100real,intent(out) :: alb ! (interpolated) TES ice albedo at that grid point
101integer :: icap ! =1: Northern hemisphere =2: Southern hemisphere
102
103! local variables:
104logical,save :: firstcall=.true.
105real,save :: zls_old ! value of zls from a previous call
106integer,save :: tinf,tsup ! encompassing time indexes of TES data
107real,save :: reltime ! relative position in-between time indexes (in [0;1])
108integer :: latinf,latsup ! encompassing latitude indexes of TES data
109real :: rellat ! relative position in-between latitude indexes (in[0;1])
110integer :: loninf,lonsup ! encompassing longitude indexes of TES data
111real :: rellon !relative position in-between longitude indexes (in[0;1])
112real,save :: pi,radeg ! to convert radians to degrees
113real :: zlsd ! solar longitude, in degrees
114real :: latd ! latitude, in degrees
115real :: lond ! longitude, in degrees
116integer :: i
117character(len=20),parameter :: modname="TES_icecap_albedo"
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  ierr=nf90_open(trim(datadir)//"/npsc_albedo.nc",NF90_NOWRITE,nid)
154  IF (ierr.NE.NF90_NOERR) THEN
155    write(*,*)'Problem opening npsc_albedo.nc (phymars/albedocaps.F90)'
156    write(*,*)'It should be in :',trim(datadir),'/'
157    write(*,*)'1) You can change this directory address in callfis.def with'
158    write(*,*)'   datadir=/path/to/datafiles'
159    write(*,*)'2) If necessary, npsc_albedo.nc (and other datafiles)'
160    write(*,*)'   can be obtained online on:'
161    write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
162    CALL abort_physic(modname,"missing input file",1)
163  ELSE
164    write(*,*) "albedocaps: using file ",trim(datadir)//"/npsc_albedo.nc"
165  ENDIF
166 
167  ierr=nf90_inq_varid(nid,"longitude",nvarid)
168  if (ierr.ne.NF90_NOERR) then
169    write(*,*) "Failed to find longitude in file!"
170    write(*,*)trim(nf90_strerror(ierr))
171    call abort_physic(modname,"failed finding longitude",1)
172  else
173    ierr=nf90_get_var(nid,nvarid,TESlon)
174    if (ierr.ne.NF90_NOERR) then
175      write(*,*) "Failed loading longitude data from file!"
176      write(*,*)trim(nf90_strerror(ierr))
177      call abort_physic(modname,"failed loading longitude",1)
178    endif
179  endif
180 
181  ierr=nf90_inq_varid(nid,"latitude",nvarid)
182  if (ierr.ne.NF90_NOERR) then
183    write(*,*) "Failed to find latitude in file!"
184    write(*,*)trim(nf90_strerror(ierr))
185    call abort_physic(modname,"failed finding latitude",1)
186  else
187    ierr=nf90_get_var(nid,nvarid,TESlatn)
188    if (ierr.ne.NF90_NOERR) then
189      write(*,*) "Failed loading latitude data from file!"
190      write(*,*)trim(nf90_strerror(ierr))
191      call abort_physic(modname,"failed loading latitude",1)
192    endif
193  endif
194 
195  ierr=nf90_inq_varid(nid,"time",nvarid)
196  if (ierr.ne.NF90_NOERR) then
197    write(*,*) "Failed to find time in file!"
198    write(*,*)trim(nf90_strerror(ierr))
199    call abort_physic(modname,"failed finding time",1)
200  else
201    ierr=nf90_get_var(nid,nvarid,TESls)
202    if (ierr.ne.NF90_NOERR) then
203      write(*,*) "Failed loading time data from file!"
204      write(*,*)trim(nf90_strerror(ierr))
205      call abort_physic(modname,"failed loading time",1)
206    endif
207  endif
208 
209  ierr=nf90_inq_varid(nid,"albedo",nvarid)
210  if (ierr.ne.NF90_NOERR) then
211    write(*,*) "Failed to find albedo in file!"
212    write(*,*)trim(nf90_strerror(ierr))
213    call abort_physic(modname,"failed finding albedo",1)
214  else
215    ierr=nf90_get_var(nid,nvarid,TESalbn)
216    if (ierr.ne.NF90_NOERR) then
217      write(*,*) "Failed loading albedo data from file!"
218      write(*,*)trim(nf90_strerror(ierr))
219      call abort_physic(modname,"failed loading albedo",1)
220    endif
221  endif
222
223  ierr=nf90_close(nid)
224
225! Load albedoes for Southern Hemisphere
226  ierr=nf90_open(trim(datadir)//"/spsc_albedo.nc",NF90_NOWRITE,nid)
227  IF (ierr.NE.NF90_NOERR) THEN
228    write(*,*)'Problem opening spsc_albedo.nc (phymars/albedocaps.F90)'
229    write(*,*)'It should be in :',trim(datadir),'/'
230    write(*,*)'1) You can change this directory address in callfis.def with'
231    write(*,*)'   datadir=/path/to/datafiles'
232    write(*,*)'2) If necessary, spsc_albedo.nc (and other datafiles)'
233    write(*,*)'   can be obtained online on:'
234    write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
235    CALL abort_physic(modname,"missing input file",1)
236  ELSE
237    write(*,*) "albedocaps: using file ",trim(datadir)//"/spsc_albedo.nc"
238  ENDIF
239
240  ierr=nf90_inq_varid(nid,"latitude",nvarid)
241  if (ierr.ne.NF90_NOERR) then
242    write(*,*) "Failed to find latitude in file!"
243    write(*,*)trim(nf90_strerror(ierr))
244    call abort_physic(modname,"failed finding latitude",1)
245  else
246    ierr=nf90_get_var(nid,nvarid,TESlats)
247    if (ierr.ne.NF90_NOERR) then
248      write(*,*) "Failed loading latitude data from file!"
249      write(*,*)trim(nf90_strerror(ierr))
250      call abort_physic(modname,"failed loading latitude",1)
251    endif
252  endif
253
254  ierr=nf90_inq_varid(nid,"albedo",nvarid)
255  if (ierr.ne.NF90_NOERR) then
256    write(*,*) "Failed to find albedo in file!"
257    write(*,*)trim(nf90_strerror(ierr))
258    call abort_physic(modname,"failed finding albedo",1)
259  else
260    ierr=nf90_get_var(nid,nvarid,TESalbs)
261    if (ierr.ne.NF90_NOERR) then
262      write(*,*) "Failed loading albedo data from file!"
263      write(*,*)trim(nf90_strerror(ierr))
264      call abort_physic(modname,"failed loading albedo",1)
265    endif
266  endif
267
268  ierr=nf90_close(nid)
269
270! constants:
271  pi=acos(-1.)
272  radeg=180/pi
273
274  zls_old=-999 ! dummy initialization
275
276  firstcall=.false.
277endif ! of if firstcall
278
279! 1. Identify encompassing latitudes
280
281! Check that latitude is such that there is TES data to use
282! (ie: latitude 45 deg and poleward) otherwise use 'default' albedoes
283latd=latitude(ig)*radeg ! latitude, in degrees
284if (icap.eq.1) then
285 ! North hemisphere
286 if (latd.lt.TESlatnmin) then
287   alb=albedice(1)
288   ! the job is done; quit this routine
289   return
290 else
291   ! find encompassing latitudes
292   if (latd.ge.TESlatn(TESlatsize)) then
293     latinf=TESlatsize
294     latsup=TESlatsize
295     rellat=0.
296   else
297     do i=1,TESlatsize-1
298       if ((latd.ge.TESlatn(i)).and.(latd.lt.TESlatn(i+1))) then
299         latinf=i
300         latsup=i+1
301         rellat=(latd-TESlatn(i))/TESdeltalat
302         exit ! found encompassing indexes; quit loop
303       endif
304     enddo
305   endif
306 endif ! of if (latd.lt.TESlatnmin)
307else ! icap=2
308 ! South hemisphere
309 if (latd.gt.TESlatsmax) then
310   alb=albedice(2)
311   ! the job is done; quit this routine
312   return
313 else
314   ! find encompassing latitudes
315   if (latd.lt.TESlats(1)) then
316     latinf=1
317     latsup=1
318     rellat=0.
319   else
320     do i=1,TESlatsize-1
321       if ((latd.ge.TESlats(i)).and.(latd.lt.TESlats(i+1))) then
322         latinf=i
323         latsup=i+1
324         rellat=(latd-TESlats(i))/TESdeltalat
325         exit ! found encompassing indexes; quit loop
326       endif
327     enddo
328   endif
329 endif ! of if (latd.gt.-45.)
330endif ! of if (icap.eq.1)
331
332! 2. Identify encompassing time indexes
333if (zls.ne.zls_old) then
334  zlsd=zls*radeg ! solar longitude, in degrees
335 
336  if (zlsd.lt.TESls(1)) then
337    tinf=TESlssize
338    tsup=1
339    reltime=0.5+zlsd/TESdeltals
340  else
341    if (zlsd.ge.TESls(TESlssize)) then
342      tinf=TESlssize
343      tsup=1
344      reltime=(360.-zlsd)/TESdeltals
345    else
346      ! look for encompassing indexes
347      do i=1,TESlssize-1
348        if ((zlsd.ge.TESls(i)).and.(zlsd.lt.TESls(i+1))) then
349          tinf=i
350          tsup=i+1
351          reltime=(zlsd-TESls(i))/TESdeltals
352          exit ! quit loop, we found the indexes
353        endif
354      enddo
355    endif
356  endif ! of if (zlsd.lt.TESls(1))
357 
358  zls_old=zls ! store current zls
359endif ! of if (zls.ne.zls_old)
360
361! 3. Identify encompassing longitudes
362lond=longitude(ig)*radeg ! east longitude, in degrees
363if (lond.lt.TESlon(1)) then
364  loninf=TESlonsize
365  lonsup=1
366  rellon=0.5+(180.+lond)/TESdeltalon
367else
368  if (lond.ge.TESlon(TESlonsize)) then
369    loninf=TESlonsize
370    lonsup=1
371    rellon=(180-lond)/TESdeltalon
372  else
373    do i=1,TESlonsize-1
374      if ((lond.ge.TESlon(i)).and.(lond.lt.TESlon(i+1))) then
375        loninf=i
376        lonsup=i+1
377        rellon=(lond-TESlon(i))/TESdeltalon
378        exit ! quit loop, we found the indexes
379      endif
380    enddo
381  endif ! of if (lond.ge.TESlon(TESlonsize))
382endif ! of if (lond.lt.TESlon(1))
383
384! 4. Use linear interpolation in time to build encompassing nodal values
385!    encompassing nodes are arranged as follow : 4 3
386!                                                1 2
387if (icap.eq.1) then
388  ! Northern hemisphere
389  val(1)=(1.-reltime)*TESalbn(loninf,latinf,tinf) &
390         +reltime*TESalbn(loninf,latinf,tsup)
391  val(2)=(1.-reltime)*TESalbn(lonsup,latinf,tinf) &
392         +reltime*TESalbn(lonsup,latinf,tsup)
393  val(3)=(1.-reltime)*TESalbn(lonsup,latsup,tinf) &
394         +reltime*TESalbn(lonsup,latsup,tsup)
395  val(4)=(1.-reltime)*TESalbn(loninf,latsup,tinf) &
396         +reltime*TESalbn(loninf,latsup,tsup)
397else
398  ! Southern hemisphere
399  val(1)=(1.-reltime)*TESalbs(loninf,latinf,tinf) &
400         +reltime*TESalbs(loninf,latinf,tsup)
401  val(2)=(1.-reltime)*TESalbs(lonsup,latinf,tinf) &
402         +reltime*TESalbs(lonsup,latinf,tsup)
403  val(3)=(1.-reltime)*TESalbs(lonsup,latsup,tinf) &
404         +reltime*TESalbs(lonsup,latsup,tsup)
405  val(4)=(1.-reltime)*TESalbs(loninf,latsup,tinf) &
406         +reltime*TESalbs(loninf,latsup,tsup)
407endif ! of if (icap.eq.1)
408
409! 5. Use bilinear interpolation to compute albedo
410alb=(1.-rellon)*(1.-rellat)*val(1) &
411    +rellon*(1.-rellat)*val(2) &
412    +rellon*rellat*val(3) &
413    +(1.-rellon)*rellat*val(4)
414
415! 6. Apply coefficient to interpolated TES albedo
416if (icap.eq.1) then
417  alb=alb*TESice_Ncoef
418else
419  alb=alb*TESice_Scoef
420endif ! of if (icap.eq.1)
421
422! Make sure that returned albedo is never greater than 0.90
423if (alb.gt.0.90) alb=0.90
424
425end subroutine TES_icecap_albedo
426
Note: See TracBrowser for help on using the repository browser.