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

Last change on this file since 4086 was 4034, checked in by emillour, 3 months ago

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