source: trunk/LMDZ.MARS/libf/phymars/readtesassim.F90 @ 292

Last change on this file since 292 was 249, checked in by emillour, 13 years ago

Mars GCM:

  • Fixed bug in readtesassim (timeflag was not always set before being used)

EM
File size: 8.2 KB
Line 
1subroutine readtesassim(ngrid,nlayer,zday,pplev,tauref)
2
3! Reading of the dust assimilation file
4
5use netcdf
6implicit none
7
8#include "dimensions.h"
9#include "dimphys.h"
10#include "comgeomfi.h"
11#include "datafile.h"
12#include "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) :: tauref
18
19! Local variables
20
21real :: realday
22integer nid,nvarid,ierr
23integer ltloop,lsloop,iloop,jloop,varloop,ig
24real, dimension(2) :: taubuf
25real tau1(4),tau
26real alt(4)
27integer latp(4),lonp(4)
28integer yinf,ysup,xinf,xsup,tinf,tsup
29real latinf,latsup,loninf,lonsup
30real latintmp,lonintmp,latdeg,londeg
31real colat,dlat,dlon,colattmp
32logical, save :: firstcall=.true.
33logical :: timeflag
34real,save :: radeg,pi
35integer :: timedim,londim,latdim
36real, dimension(:), allocatable, save :: lat,lon,time
37real, dimension(:,:,:), allocatable, save :: tautes
38integer, save :: timelen,lonlen,latlen
39character(len=33) :: filename
40
41realday=mod(zday,669.)
42
43if (firstcall) then
44   firstcall=.false.
45
46   pi=acos(-1.)
47   radeg=180/pi
48   
49   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
50   ! iaervar=4 means read dust_tes.nc file
51   ! iaervar=24 means read dust_tes_MY24.nc file
52   ! iaervar=25 means read dust_tes_MY25.nc file
53   ! iaervar=26 means read dust_tes_MY26.nc file
54   if (iaervar.eq.4) then
55     filename="dust_tes.nc"
56   else if (iaervar.eq.24) then
57     filename="dust_tes_MY24.nc" ! which currently is identical to dust_tes.nc
58   else if (iaervar.eq.25) then
59     filename="dust_tes_MY25.nc"
60   else if (iaervar.eq.26) then
61     filename="dust_tes_MY26.nc"
62   endif
63   
64   ! Note: datafile() is defined in "datafile.h"
65   !ierr=NF_OPEN(trim(datafile)//"/"//trim(filename),NF_NOWRITE,nid)
66   ierr=nf90_open(trim(datafile)//"/"//trim(filename),NF90_NOWRITE,nid)
67   IF (ierr.NE.nf90_noerr) THEN
68     write(*,*)'Problem opening ',trim(filename),' (in phymars/readtesassim.F90)'
69     write(*,*)'It should be in :',trim(datafile),'/'
70     write(*,*)'1) You can change this directory address in '
71     write(*,*)'   file phymars/datafile.h'
72     write(*,*)'2) If necessary, dust.nc (and other datafiles)'
73     write(*,*)'   can be obtained online on:'
74     write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
75     CALL ABORT
76   ENDIF
77
78   ierr=nf90_inq_dimid(nid,"Time",timedim)
79   ierr=nf90_inquire_dimension(nid,timedim,len=timelen)
80   ierr=nf90_inq_dimid(nid,"latitude",latdim)
81   ierr=nf90_inquire_dimension(nid,latdim,len=latlen)
82   ierr=nf90_inq_dimid(nid,"longitude",londim)
83   ierr=nf90_inquire_dimension(nid,londim,len=lonlen)
84
85
86   allocate(tautes(lonlen,latlen,timelen))
87   allocate(lat(latlen), lon(lonlen), time(timelen))
88
89   ierr=nf90_inq_varid(nid,"dustop",nvarid)
90   ierr=nf90_get_var(nid,nvarid,tautes)
91   IF (ierr .NE. nf90_noerr) THEN
92      PRINT*, "Error: Readtesassim <dustop> not found"
93      write(*,*)trim(nf90_strerror(ierr))
94      stop
95   ENDIF
96
97   ierr=nf90_inq_varid(nid,"Time",nvarid)
98   ierr=nf90_get_var(nid,nvarid,time)
99   IF (ierr .NE. nf90_noerr) THEN
100      PRINT*, "Error: Readtesassim <Time> not found"
101      write(*,*)trim(nf90_strerror(ierr))
102      stop
103   ENDIF
104
105   ierr=nf90_inq_varid(nid,"latitude",nvarid)
106   ierr=nf90_get_var(nid,nvarid,lat)
107   IF (ierr .NE. nf90_noerr) THEN
108      PRINT*, "Error: Readtesassim <latitude> not found"
109      write(*,*)trim(nf90_strerror(ierr))
110      stop
111   ENDIF
112
113   ierr=nf90_inq_varid(nid,"longitude",nvarid)
114   ierr=nf90_get_var(nid,nvarid,lon)
115   IF (ierr .NE. nf90_noerr) THEN
116      PRINT*, "Error: Readtesassim <longitude> not found"
117      write(*,*)trim(nf90_strerror(ierr))
118      stop
119   ENDIF
120
121   ierr=nf90_close(nid)
122
123endif ! of if (firstcall)
124
125do ig=1,ngrid
126
127! Find the four nearest points, arranged as follows:
128!                               1 2
129!                               3 4
130
131   latdeg=lati(ig)*radeg  ! latitude, in degrees
132   londeg=long(ig)*radeg  ! longitude, in degrees east
133   colat=90-latdeg        ! colatitude, in degrees
134
135! Ehouarn: rounding effects and/or specific compiler issues
136!          sometime cause londeg to be sligthly below -180 ...
137  if (londeg.lt.-180) then
138    ! check if it is by a large amount
139    if ((londeg+180).lt.-1.e-3) then
140      write(*,*) 'reattesassim: error!!'
141      write(*,*) ' ig=',ig,' londeg=',londeg
142      stop
143    else
144      londeg=-180
145    endif
146  endif
147
148 ! Find encompassing latitudes
149   if (colat<(90-lat(1))) then ! between north pole and lat(1)
150      ysup=1
151      yinf=1
152   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
153      ysup=latlen
154      yinf=latlen
155   else ! general case
156      do iloop=2,latlen
157         if(colat<(90-lat(iloop))) then
158           ysup=iloop-1
159           yinf=iloop
160           exit
161         endif
162      enddo
163   endif
164
165   latinf=lat(yinf)
166   latsup=lat(ysup)
167
168
169 ! Find encompassing longitudes
170   ! Note: in input file, lon(1)=-180.
171   if (londeg>lon(lonlen)) then
172      xsup=1
173      xinf=lonlen
174      loninf=lon(xsup)
175      lonsup=180.0 ! since lon(1)=-180.0
176   else
177      do iloop=2,lonlen
178         if(londeg<lon(iloop)) then
179              xsup=iloop
180              xinf=iloop-1
181              exit
182         endif
183      enddo
184      loninf=lon(xinf)
185      lonsup=lon(xsup)
186   endif
187
188   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
189     .OR.(ysup.lt.1)) then
190      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
191      write (*,*) "xinf: ",xinf
192      write (*,*) "xsup: ",xsup
193      write (*,*) "yinf: ",yinf
194      write (*,*) "ysup: ",ysup
195      stop
196   endif
197
198!   loninf=lon(xinf)
199!   lonsup=lon(xsup)
200!   latinf=lat(yinf)
201!   latsup=lat(ysup)
202
203! The four neighbouring points are arranged as follows:
204!                               1 2
205!                               3 4
206
207   latp(1)=ysup
208   latp(2)=ysup
209   latp(3)=yinf
210   latp(4)=yinf
211
212   lonp(1)=xinf
213   lonp(2)=xsup
214   lonp(3)=xinf
215   lonp(4)=xsup
216
217! Linear interpolation on time, for all four neighbouring points
218
219   if ((realday<time(1)).or.(realday>time(timelen))) then
220      tinf=timelen
221      tsup=1
222      timeflag=.true.
223   else
224      timeflag=.false.
225      do iloop=2,timelen
226         if (realday<time(iloop)) then
227            tinf=iloop-1
228            tsup=iloop
229            exit
230         endif
231      enddo
232   endif
233
234! Bilinear interpolation on the four nearest points
235
236   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
237      dlat=0
238   else
239      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
240   endif
241
242   if (lonsup==loninf) then
243      dlon=0
244   else
245      dlon=(londeg-loninf)/(lonsup-loninf)
246   endif
247
248   do iloop=1,4
249      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
250      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
251      if (timeflag) then
252         if (realday>time(timelen)) then
253            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf)))
254         else
255            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
256         endif
257      else
258         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
259      endif
260      if (tau1(iloop)<0) then
261          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
262          write (*,*) "utime ",realday
263          write (*,*) "time(tinf) ",time(tinf)
264          write (*,*) "time(tsup) ",time(tsup)
265          write (*,*) "tau1 ",taubuf(1)
266          write (*,*) "tau2 ",taubuf(2)
267          write (*,*) "tau ",tau1(iloop)
268          stop
269      endif
270   enddo
271
272   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
273      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
274      write (*,*) "dlat: ",dlat
275      write (*,*) "lat: ",latdeg
276      write (*,*) "dlon: ",dlon
277      write (*,*) "lon: ",londeg
278      stop
279   endif
280
281   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
282
283! Multiply by correction coefficient (value depends on SW radiative transfer)
284! NB: swrtype is set in callkeys.h
285   if (swrtype.eq.1) then ! Fouquart
286     tauref(ig)=tau*0.825
287   else
288     if (swrtype.eq.2) then ! Toon
289       tauref(ig)=tau*1.3
290     else
291       write(*,*) "readtesassim: wrong value for flag swrtype !"
292       stop
293     endif
294   endif ! of if (swrtype.eq.1)
295enddo
296
297end
Note: See TracBrowser for help on using the repository browser.