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

Last change on this file since 1292 was 1278, checked in by emillour, 11 years ago

Mars GCM:

  • Add possibility to use clim or MY31 dust scenarios (and realistic EUV for MY31).

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