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

Last change on this file since 3567 was 3418, checked in by llange, 4 months ago

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