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

Last change on this file since 3581 was 3581, checked in by jbclement, 2 weeks ago

Mars PCM:

  • Deletion of the variable 'albedo_co2_cap' which was redundant with 'albedo_perennialco2'.
  • 'albedo_perennialco2' is now a vector of dimension 2 to account for different albedos of CO2 ice between the northern and southern hemispheres (by default 0.6 and 0.85 respectively).

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