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

Last change on this file since 2616 was 2609, checked in by romain.vande, 3 years ago

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