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

Last change on this file since 3599 was 3586, checked in by jbclement, 8 days ago

Mars PCM:

  • Albedo is now set properly in every situation given the presence of CO2/H2O ice and frost + water albedo is done solely in "physiq".
  • Small update to make the display of code version status/diff clearer.

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