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

Last change on this file since 3225 was 3159, checked in by jbclement, 14 months ago

PEM:

  • A number of corrections related to r3149 for CO2 ice management:

    'co2_ice' was not transferred to 'perennial_co2ice' at the end of the PEM run;
    In the Mars PCM, 'perennial_co2ice' is now correctly handled regarding the frost in case of CO2 sublimation/condensation.

  • Addition of (CO2 and H2O) ice metamorphism: if the frost is above a given value, then the excess frost is transferred to the perennial ice.
  • Thresholds value for ice management can now be set in the "run_PEM.def".

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