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

Last change on this file since 1944 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

  • remove "comorbit.h" since it is no longer used.
  • turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

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