source: trunk/WRF.COMMON/WRFV2/mars_lmd/libf/phymars/meso_readtesassim.F90

Last change on this file was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 8.0 KB
Line 
1subroutine meso_readtesassim(ngrid,nlayer,zday,pplev,tauref,day_translate)
2
3! Reading of the dust assimilation file
4
5implicit none
6
7include "dimensions.h"
8include "dimphys.h"
9include "comgeomfi.h"
10include "netcdf.inc"
11include "datafile.h"
12
13integer, intent(in) :: ngrid,nlayer
14real, intent(in) :: zday ! date in martian days
15real, dimension(ngrid,nlayer+1), intent(in) :: pplev
16real, dimension(ngrid), intent(out) :: tauref
17
18! Local variables
19
20real :: realday
21integer nid,nvarid,ierr
22integer ltloop,lsloop,iloop,jloop,varloop,ig
23real, dimension(2) :: taubuf
24real tau1(4),tau
25real alt(4)
26integer latp(4),lonp(4)
27integer yinf,ysup,xinf,xsup,tinf,tsup
28real latinf,latsup,loninf,lonsup
29real latintmp,lonintmp,latdeg,londeg
30real colat,dlat,dlon,colattmp
31logical, save :: firstcall=.true.
32logical :: timeflag
33real radeg,pi
34integer :: timedim,londim,latdim
35real, dimension(:), allocatable, save :: lat,lon,time
36real, dimension(:,:,:), allocatable, save :: tautes
37integer, save :: timelen,lonlen,latlen
38INTEGER, external:: lnblnk
39
40!!WRF added
41INTEGER :: day_translate
42
43pi=acos(-1.)
44radeg=180/pi
45realday=mod(zday,669.)
46
47if (firstcall) then
48   firstcall=.false.
49
50   !! Note: datafile() is defined in "datafile.h"
51   !print *,datafile(1:lnblnk(datafile))
52   !ierr=NF_OPEN (datafile(1:lnblnk(datafile))//"/input_totoptdep.nc",NF_NOWRITE,nid)
53
54!!-------------------------------------
55!! iaervar is used here to communicate the day starting the diagfi
56!! -- in addition, inpu_totoptdep.nc must be prepared
57day_translate=day_translate-1000
58ierr=NF_OPEN ("./input_totoptdep.nc",NF_NOWRITE,nid)
59!!-------------------------------------
60
61
62   IF (ierr.NE.NF_NOERR) THEN
63     !write(*,*)'Problem opening dust.nc (in phymars/readtesassim.F90)'
64     !write(*,*)'It should be in :',datafile(1:lnblnk(datafile)),'/'
65     !write(*,*)'1) You can change this directory address in '
66     !write(*,*)'   file phymars/datafile.h'
67     !write(*,*)'2) If necessary, dust.nc (and other datafiles)'
68     !write(*,*)'   can be obtained online on:'
69     !write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
70     write(*,*)'Problem opening input_totoptdep.nc ...'
71     CALL ABORT
72   ENDIF
73
74   ierr=NF_INQ_DIMID(nid,"time",timedim)
75   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
76   ierr=NF_INQ_DIMID(nid,"lat",latdim)
77   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
78   ierr=NF_INQ_DIMID(nid,"lon",londim)
79   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
80
81   allocate(tautes(lonlen,latlen,timelen))
82   allocate(lat(latlen), lon(lonlen), time(timelen))
83
84   ierr = NF_INQ_VARID (nid, "totoptdep",nvarid)
85#ifdef NC_DOUBLE
86   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tautes)
87#else
88   ierr = NF_GET_VAR_REAL(nid, nvarid, tautes)
89#endif
90   IF (ierr .NE. NF_NOERR) THEN
91      PRINT*, "Error: Readtesassim <totoptdep> not found"
92      stop
93   ENDIF
94
95!!****WRF
96!! diagfi files feature the dust opacity divided by 700 Pa
97        tautes=tautes*700.
98!!****WRF
99
100
101   ierr = NF_INQ_VARID (nid, "time",nvarid)
102#ifdef NC_DOUBLE
103   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
104#else
105   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
106#endif
107   IF (ierr .NE. NF_NOERR) THEN
108      PRINT*, "Error: Readtesassim <Time> not found"
109      stop
110   ENDIF
111
112!!****WRF
113PRINT *, "**** read assimilation dust"
114PRINT *, "prescribed starting day ...", day_translate
115time=time+day_translate
116!!****WRF
117
118   ierr = NF_INQ_VARID (nid, "lat",nvarid)
119#ifdef NC_DOUBLE
120   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
121#else
122   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
123#endif
124   IF (ierr .NE. NF_NOERR) THEN
125      PRINT*, "Error: Readtesassim <latitude> not found"
126      stop
127   ENDIF
128
129   ierr = NF_INQ_VARID (nid, "lon",nvarid)
130#ifdef NC_DOUBLE
131   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
132#else
133   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
134#endif
135   IF (ierr .NE. NF_NOERR) THEN
136      PRINT*, "Error: Readtesassim <longitude> not found"
137      stop
138   ENDIF
139
140endif ! of if (firstcall)
141
142do ig=1,ngrid
143
144! Find the four nearest points, arranged as follows:
145!                               1 2
146!                               3 4
147
148   latdeg=lati(ig)*radeg  ! latitude, in degrees
149   londeg=long(ig)*radeg  ! longitude, in degrees east
150   colat=90-latdeg        ! colatitude, in degrees
151
152 ! Find encompassing latitudes
153   if (colat<(90-lat(1))) then ! between north pole and lat(1)
154      ysup=1
155      yinf=1
156   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
157      ysup=latlen
158      yinf=latlen
159   else ! general case
160      do iloop=2,latlen
161         if(colat<(90-lat(iloop))) then
162           ysup=iloop-1
163           yinf=iloop
164           exit
165         endif
166      enddo
167   endif
168
169   latinf=lat(yinf)
170   latsup=lat(ysup)
171
172
173 ! Find encompassing longitudes
174   ! Note: in input file, lon(1)=-180.
175   if (londeg>lon(lonlen)) then
176      xsup=1
177      xinf=lonlen
178      loninf=lon(xsup)
179      lonsup=180.0 ! since lon(1)=-180.0
180   else
181      do iloop=2,lonlen
182         if(londeg<lon(iloop)) then
183              xsup=iloop
184              xinf=iloop-1
185              exit
186         endif
187      enddo
188      loninf=lon(xinf)
189      lonsup=lon(xsup)
190   endif
191
192   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
193     .OR.(ysup.lt.1)) then
194      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
195      write (*,*) "xinf: ",xinf
196      write (*,*) "xsup: ",xsup
197      write (*,*) "yinf: ",yinf
198      write (*,*) "ysup: ",ysup
199      stop
200   endif
201
202!   loninf=lon(xinf)
203!   lonsup=lon(xsup)
204!   latinf=lat(yinf)
205!   latsup=lat(ysup)
206
207! The four neighbouring points are arranged as follows:
208!                               1 2
209!                               3 4
210
211   latp(1)=ysup
212   latp(2)=ysup
213   latp(3)=yinf
214   latp(4)=yinf
215
216   lonp(1)=xinf
217   lonp(2)=xsup
218   lonp(3)=xinf
219   lonp(4)=xsup
220
221! Linear interpolation on time, for all four neighbouring points
222
223   if ((realday<time(1)).or.(realday>time(timelen))) then
224      tinf=timelen
225      tsup=1
226      timeflag=.true.
227   else
228      do iloop=2,timelen
229         if (realday<time(iloop)) then
230            tinf=iloop-1
231            tsup=iloop
232            exit
233         endif
234      enddo
235   endif
236
237! Bilinear interpolation on the four nearest points
238
239   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
240      dlat=0
241   else
242      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
243   endif
244
245   if (lonsup==loninf) then
246      dlon=0
247   else
248      dlon=(londeg-loninf)/(lonsup-loninf)
249   endif
250
251   do iloop=1,4
252      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
253      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
254      if (timeflag) then
255         if (realday>time(timelen)) then
256            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf)))
257         else
258            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
259         endif
260      else
261         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
262      endif
263      if (tau1(iloop)<0) then
264          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
265          write (*,*) "utime ",realday
266          write (*,*) "time(tinf) ",time(tinf)
267          write (*,*) "time(tsup) ",time(tsup)
268          write (*,*) "tau1 ",taubuf(1)
269          write (*,*) "tau2 ",taubuf(2)
270          write (*,*) "tau ",tau1(iloop)
271          stop
272      endif
273   enddo
274   timeflag=.false.
275
276   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
277      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
278      write (*,*) "dlat: ",dlat
279      write (*,*) "lat: ",latdeg
280      write (*,*) "dlon: ",dlon
281      write (*,*) "lon: ",londeg
282      stop
283   endif
284
285   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
286
287!!   tauref(ig)=tau*700/pplev(ig,1)*0.8
288    tauref(ig)=tau*0.825
289
290!-------------
291! from Oxford
292!-------------
293! + tuning ... TES is little bit too hot compared to radio-occultations
294! 0.825 is the MY24 value
295!tauref(ig)=tau*700.*0.825
296!tauref(ig)=tau*700.
297!tauref(ig)=tau*700.*1.2
298!tauref(ig)=tau*700.*0.6
299!-------------------------------------------------
300! NB: multiply by 700. is now done elsewhere
301!-------------------------------------------------
302
303
304enddo
305
306end
Note: See TracBrowser for help on using the repository browser.