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

Last change on this file since 3400 was 3343, checked in by llange, 6 months ago

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