source: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario_mod.F90 @ 3567

Last change on this file since 3567 was 2921, checked in by emillour, 21 months ago

Mars PCM:
Add the possibility to use MY36 dust scenario and MY36 EUV inputs
EM
File size: 12.3 KB
Line 
1module read_dust_scenario_mod
2
3implicit none
4
5contains
6
7subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,IRtoVIScoef,tau_pref_scenario)
8
9! Reading of the dust scenario file
10
11use netcdf
12use geometry_mod, only: latitude, longitude ! in radians
13use datafile_mod, only: datadir
14use dust_param_mod, only: odpref, dustscaling_mode
15use planete_h, only: year_day
16USE mod_phys_lmdz_transfert_para, ONLY: bcast
17implicit none
18
19include "callkeys.h"
20
21integer, intent(in) :: ngrid,nlayer
22real, intent(in) :: zday ! date in martian days
23real, dimension(ngrid,nlayer+1), intent(in) :: pplev
24real, dimension(ngrid), intent(in) :: IRtoVIScoef ! conversion coefficient to apply on
25                                       ! scenario absorption IR (9.3um) CDOD
26                                       ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
27real, dimension(ngrid), intent(out) :: tau_pref_scenario ! visible dust column
28                                       ! opacity at odpref from scenario
29
30! Local variables
31
32real :: realday
33integer nid,nvarid,ierr
34integer ltloop,lsloop,iloop,jloop,varloop,ig
35real, dimension(2) :: taubuf
36real tau1(4),tau
37real alt(4)
38integer latp(4),lonp(4)
39integer yinf,ysup,xinf,xsup,tinf,tsup
40real latinf,latsup,loninf,lonsup
41real latintmp,lonintmp,latdeg,londeg
42real colat,dlat,dlon,colattmp
43logical, save :: firstcall=.true.
44logical :: timeflag
45real,save :: radeg,pi
46integer :: timedim,londim,latdim
47real, dimension(:), allocatable, save :: lat,lon,time
48real, dimension(:,:,:), allocatable, save :: tautes
49integer, save :: timelen,lonlen,latlen
50character(len=33),save :: filename
51logical,save :: IRscenario ! =true if the scenarios are IR absorption opacity (variable 'cdod')
52                           ! and not VIS extinction (variable 'dustop')
53
54!$OMP THREADPRIVATE(firstcall,radeg,pi, &
55!$OMP               lat,lon,time,tautes,  &
56!$OMP               timelen,lonlen,latlen,  &
57!$OMP               filename,IRscenario)
58
59realday=mod(zday,669.)
60
61! AS: firstcall OK absolute
62if (firstcall) then
63   firstcall=.false.
64
65   pi=acos(-1.)
66   radeg=180/pi
67   
68   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
69   ! iaervar=4 means read dust_tes.nc file
70   ! iaervar=6 means read dust_cold.nc file
71   ! iaervar=7 means read dust_warm.nc file
72   ! iaervar=8 means read dust_clim.nc file
73   ! iaervar=24 means read dust_MY24.nc file
74   ! iaervar=25 means read dust_MY25.nc file
75   ! iaervar=26 means read dust_MY26.nc file, etc.
76   if (iaervar.eq.4) then
77   ! NB: 4: old TES assimilated MY24 dust scenarios (at 700Pa ref pressure!)
78     filename="dust_tes.nc"
79   else if (iaervar.eq.6) then
80     filename="dust_cold.nc"
81   else if (iaervar.eq.7) then
82     filename="dust_warm.nc"
83   else if (iaervar.eq.8) then
84     filename="dust_clim.nc"
85   else if ((iaervar.ge.24).and.(iaervar.le.36)) then
86     write(filename,'(a7,i2,a3)')"dust_MY",iaervar,".nc"
87   ! 124,125,126: old TES assimilated dust scenarios (at 700Pa ref pressure!)
88   else if  (iaervar.eq.124) then
89     filename="dust_tes_MY24.nc"
90   else if (iaervar.eq.125) then
91     filename="dust_tes_MY25.nc"
92   else if (iaervar.eq.126) then
93     filename="dust_tes_MY26.nc"
94   endif
95
96
97!$OMP MASTER
98   
99   ierr=nf90_open(trim(datadir)//"/"//trim(filename),NF90_NOWRITE,nid)
100   IF (ierr.NE.nf90_noerr) THEN
101     write(*,*)'Problem opening ',trim(filename),' (in phymars/read_dust_scenario.F90)'
102     write(*,*)'It should be in :',trim(datadir),'/'
103     write(*,*)'1) You can change this directory address in callfis.def with'
104     write(*,*)'   datadir=/path/to/datafiles'
105     write(*,*)'2) If necessary, dust*.nc (and other datafiles)'
106     write(*,*)'   can be obtained online on:'
107     write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
108     CALL abort_physic("read_dust_scenario","Cannot find file",1)
109   ENDIF
110
111   ierr=nf90_inq_dimid(nid,"Time",timedim)
112   if (ierr.ne.nf90_noerr) then
113    ! 'Time' dimension not found, look for 'time'
114    ierr=nf90_inq_dimid(nid,"time",timedim)
115    if (ierr.ne.nf90_noerr) then
116      write(*,*)"Error: read_dust_scenario <time> not found"
117    endif
118   endif
119   ierr=nf90_inquire_dimension(nid,timedim,len=timelen)
120   ierr=nf90_inq_dimid(nid,"latitude",latdim)
121   ierr=nf90_inquire_dimension(nid,latdim,len=latlen)
122   ierr=nf90_inq_dimid(nid,"longitude",londim)
123   ierr=nf90_inquire_dimension(nid,londim,len=lonlen)
124
125!$OMP END MASTER
126
127   call bcast(lonlen)
128   call bcast(latlen)
129   call bcast(timelen)
130
131   if (.not.allocated(tautes)) allocate(tautes(lonlen,latlen,timelen))
132   if (.not.allocated(lat)) allocate(lat(latlen), lon(lonlen), time(timelen))
133
134!$OMP MASTER
135
136   ! "dustop" if loading visible extinction opacity
137   ! "cdod" if loading IR absorption opacity
138   ierr=nf90_inq_varid(nid,"dustop",nvarid)
139   if (ierr.eq.nf90_noerr) then
140     IRscenario = .false.
141     ierr=nf90_get_var(nid,nvarid,tautes)
142     IF (ierr .NE. nf90_noerr) THEN
143        PRINT*, "Error: read_dust_scenario <dustop> not found"
144        write(*,*)trim(nf90_strerror(ierr))
145        call abort_physic("read_dust_scenario","dustop not found",1)
146     ENDIF
147   else
148     ! did not find "dustop" , look for "cdod"
149     IRscenario = .true.
150     ierr=nf90_inq_varid(nid,"cdod",nvarid)
151     ierr=nf90_get_var(nid,nvarid,tautes)
152     IF (ierr .NE. nf90_noerr) THEN
153        PRINT*, "Error: read_dust_scenario <cdod> not found"
154        write(*,*)trim(nf90_strerror(ierr))
155        call abort_physic("read_dust_scenario","cdod not found",1)
156     ENDIF
157   endif
158
159   ierr=nf90_inq_varid(nid,"Time",nvarid)
160   ierr=nf90_get_var(nid,nvarid,time)
161   IF (ierr .NE. nf90_noerr) THEN
162      PRINT*, "Error: read_dust_scenario <Time> not found"
163      write(*,*)trim(nf90_strerror(ierr))
164      call abort_physic("read_dust_scenario","Time not found",1)
165   ENDIF
166
167   ! Adapt time axis values if using dustinjection scheme and/or
168   ! dustscaling_mode == 2
169   IF ((dustinjection>0).or.(dustscaling_mode==2)) THEN
170     ! Sanity check: we expect to have one map per sol
171     if (timelen/=year_day) then
172       write(*,*)"read_dust_scenario error: timelen=",timelen
173       write(*,*)" whereas it was exepceted to be year_day=",year_day
174       call abort_physic("read_dust_scenario","timelen/=year_day",1)
175     endif
176     ! fill time() with intergers from 0 to timelen-1
177     write(*,*)"read_dust_scenario: adapting time axis to integer values"
178     do iloop=1,timelen
179       time(iloop)=iloop-1
180     enddo
181   ENDIF ! of IF ((dustinjection>0).or.(dustscaling_mode==2))
182
183   ierr=nf90_inq_varid(nid,"latitude",nvarid)
184   ierr=nf90_get_var(nid,nvarid,lat)
185   IF (ierr .NE. nf90_noerr) THEN
186      PRINT*, "Error: read_dust_scenario <latitude> not found"
187      write(*,*)trim(nf90_strerror(ierr))
188      call abort_physic("read_dust_scenario","latitude not found",1)
189   ENDIF
190
191   ierr=nf90_inq_varid(nid,"longitude",nvarid)
192   ierr=nf90_get_var(nid,nvarid,lon)
193   IF (ierr .NE. nf90_noerr) THEN
194      PRINT*, "Error: read_dust_scenario <longitude> not found"
195      write(*,*)trim(nf90_strerror(ierr))
196      call abort_physic("read_dust_scenario","longitude not found",1)
197   ENDIF
198
199   ierr=nf90_close(nid)
200
201!$OMP END MASTER
202
203   call bcast(tautes)
204   call bcast(time)
205   call bcast(lat)
206   call bcast(lon)
207   call bcast(IRscenario)
208
209endif ! of if (firstcall)
210
211
212do ig=1,ngrid
213
214! Find the four nearest points, arranged as follows:
215!                               1 2
216!                               3 4
217
218   latdeg=latitude(ig)*radeg  ! latitude, in degrees
219   londeg=longitude(ig)*radeg  ! longitude, in degrees east
220   colat=90-latdeg        ! colatitude, in degrees
221
222! Ehouarn: rounding effects and/or specific compiler issues
223!          sometime cause londeg to be sligthly below -180 ...
224  if (londeg.lt.-180) then
225    ! check if it is by a large amount
226    if ((londeg+180).lt.-1.e-3) then
227      write(*,*) ' ig=',ig,' londeg=',londeg
228      call abort_physic("read_dust_scenario","longitude way out of range",1)
229    else
230      londeg=-180
231    endif
232  endif
233
234 ! Find encompassing latitudes
235   if (colat<(90-lat(1))) then ! between north pole and lat(1)
236      ysup=1
237      yinf=1
238   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
239      ysup=latlen
240      yinf=latlen
241   else ! general case
242      do iloop=2,latlen
243         if(colat<(90-lat(iloop))) then
244           ysup=iloop-1
245           yinf=iloop
246           exit
247         endif
248      enddo
249   endif
250
251   latinf=lat(yinf)
252   latsup=lat(ysup)
253
254
255 ! Find encompassing longitudes
256   ! Note: in input file, lon(1)=-180.
257   if (londeg>lon(lonlen)) then
258      xsup=1
259      xinf=lonlen
260      loninf=lon(xsup)
261      lonsup=180.0 ! since lon(1)=-180.0
262   else
263      do iloop=2,lonlen
264         if(londeg<lon(iloop)) then
265              xsup=iloop
266              xinf=iloop-1
267              exit
268         endif
269      enddo
270      loninf=lon(xinf)
271      lonsup=lon(xsup)
272   endif
273
274   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
275     .OR.(ysup.lt.1)) then
276      write (*,*) "read_dust_scenario: SYSTEM ERROR on x or y in ts_gcm"
277      write (*,*) "xinf: ",xinf
278      write (*,*) "xsup: ",xsup
279      write (*,*) "yinf: ",yinf
280      write (*,*) "ysup: ",ysup
281      call abort_physic("read_dust_scenario","invalid coordinates",1)
282   endif
283
284!   loninf=lon(xinf)
285!   lonsup=lon(xsup)
286!   latinf=lat(yinf)
287!   latsup=lat(ysup)
288
289! The four neighbouring points are arranged as follows:
290!                               1 2
291!                               3 4
292
293   latp(1)=ysup
294   latp(2)=ysup
295   latp(3)=yinf
296   latp(4)=yinf
297
298   lonp(1)=xinf
299   lonp(2)=xsup
300   lonp(3)=xinf
301   lonp(4)=xsup
302
303! Linear interpolation on time, for all four neighbouring points
304
305   if ((realday<time(1)).or.(realday>time(timelen))) then
306      tinf=timelen
307      tsup=1
308      timeflag=.true.
309   else
310      timeflag=.false.
311      do iloop=2,timelen
312         if (realday<=time(iloop)) then
313            tinf=iloop-1
314            tsup=iloop
315            exit
316         endif
317      enddo
318   endif
319
320! Bilinear interpolation on the four nearest points
321
322   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
323      dlat=0
324   else
325      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
326   endif
327
328   if (lonsup==loninf) then
329      dlon=0
330   else
331      dlon=(londeg-loninf)/(lonsup-loninf)
332   endif
333
334   do iloop=1,4
335      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
336      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
337      if (timeflag) then
338         if (realday>time(timelen)) then
339            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf)))
340         else
341            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
342         endif
343      else
344         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
345      endif
346      if (tau1(iloop)<0) then
347          write (*,*) "read_dust_scenario: SYSTEM ERROR on tau"
348          write (*,*) "utime ",realday
349          write (*,*) "time(tinf) ",time(tinf)
350          write (*,*) "time(tsup) ",time(tsup)
351          write (*,*) "tau1 ",taubuf(1)
352          write (*,*) "tau2 ",taubuf(2)
353          write (*,*) "tau ",tau1(iloop)
354          call abort_physic("read_dust_scenario","negative tau!",1)
355      endif
356   enddo
357
358   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
359      write (*,*) "read_dust_scenario: SYSTEM ERROR on dlat or dlon in ts_gcm"
360      write (*,*) "dlat: ",dlat
361      write (*,*) "lat: ",latdeg
362      write (*,*) "dlon: ",dlon
363      write (*,*) "lon: ",londeg
364      call abort_physic("read_dust_scenario","negative dlat or dlon!",1)
365   endif
366
367   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
368
369   if (IRscenario) then
370     ! if the scenarios are in infrared, multiply by IRtoVIScoef(ig)
371     ! to convert from IR absorption to visible extinction opacity
372     tau_pref_scenario(ig)=tau * IRtoVIScoef(ig)
373   else
374     tau_pref_scenario(ig)=tau
375   endif
376   
377!
378enddo ! of ig=1,ngrid
379
380if (filename(1:8)=="dust_tes") then
381  ! when using old TES files, correction for:
382  ! - the reference pressure was 700Pa (unlike 610Pa now)
383  ! - the 1.3 conversion factor from IR absorbtion opacity to
384  !   IR extinction opacity
385  tau_pref_scenario(:)=tau_pref_scenario(:)*1.3*(odpref/700.)
386endif
387
388if (swrtype.eq.1) then ! Fouquart (NB: swrtype is set in callkeys.h)
389 ! when using old radiative transfer (like in MCD 4.x)
390 ! needed to decrease opacity (*0.825) to compensate overestimation of
391 ! heating rates
392  tau_pref_scenario(:)=tau_pref_scenario(:)*0.825/1.3
393endif
394
395end subroutine read_dust_scenario
396
397end module read_dust_scenario_mod
Note: See TracBrowser for help on using the repository browser.