source: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90 @ 2437

Last change on this file since 2437 was 2425, checked in by emillour, 5 years ago

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