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

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

Mars GCM:

  • Update of the read_dust_scenario routine: when input dust scenarios file contain variable "dustop", it is assumed to be visible extinction opacity, and if it is "cdod" (recent change, for the MCDv5.1 dust scenarios), then it is IR absorption opacity (and is multiplied by 2.6 to be converted to visible extinction opacity).

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