source: trunk/MESOSCALE/LMD_MM_MARS/SRC/DEV/readtesassim.F90.old @ 3567

Last change on this file since 3567 was 1243, checked in by aslmd, 11 years ago

LMDZ.MARS & MESOSCALE. Bye Bye meso_inc in libf/phymars
File size: 6.9 KB
Line 
1subroutine readtesassim(ngrid,nlayer,zday,pplev,tauref)
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
40pi=acos(-1.)
41radeg=180/pi
42realday=mod(zday,669.)
43
44if (firstcall) then
45   firstcall=.false.
46
47   ! Note: datafile() is defined in "datafile.h"
48   ierr=NF_OPEN (datafile(1:lnblnk(datafile))//"/dust_tes.nc",NF_NOWRITE,nid)
49   IF (ierr.NE.NF_NOERR) THEN
50     write(*,*)'Problem opening dust.nc (in phymars/readtesassim.F90)'
51     write(*,*)'It should be in :',datafile(1:lnblnk(datafile)),'/'
52     write(*,*)'1) You can change this directory address in '
53     write(*,*)'   file phymars/datafile.h'
54     write(*,*)'2) If necessary, dust.nc (and other datafiles)'
55     write(*,*)'   can be obtained online on:'
56     write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
57     CALL ABORT
58   ENDIF
59
60   ierr=NF_INQ_DIMID(nid,"Time",timedim)
61   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
62   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
63   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
64   ierr=NF_INQ_DIMID(nid,"longitude",londim)
65   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
66
67   allocate(tautes(lonlen,latlen,timelen))
68   allocate(lat(latlen), lon(lonlen), time(timelen))
69
70   ierr = NF_INQ_VARID (nid, "dustop",nvarid)
71#ifdef NC_DOUBLE
72   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tautes)
73#else
74   ierr = NF_GET_VAR_REAL(nid, nvarid, tautes)
75#endif
76   IF (ierr .NE. NF_NOERR) THEN
77      PRINT*, "Error: Readtesassim <dustop> not found"
78      stop
79   ENDIF
80   ierr = NF_INQ_VARID (nid, "Time",nvarid)
81#ifdef NC_DOUBLE
82   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
83#else
84   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
85#endif
86   IF (ierr .NE. NF_NOERR) THEN
87      PRINT*, "Error: Readtesassim <Time> not found"
88      stop
89   ENDIF
90
91   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
92#ifdef NC_DOUBLE
93   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
94#else
95   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
96#endif
97   IF (ierr .NE. NF_NOERR) THEN
98      PRINT*, "Error: Readtesassim <latitude> not found"
99      stop
100   ENDIF
101
102   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
103#ifdef NC_DOUBLE
104   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
105#else
106   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
107#endif
108   IF (ierr .NE. NF_NOERR) THEN
109      PRINT*, "Error: Readtesassim <longitude> not found"
110      stop
111   ENDIF
112
113endif ! of if (firstcall)
114
115do ig=1,ngrid
116
117! Find the four nearest points, arranged as follows:
118!                               1 2
119!                               3 4
120
121   latdeg=lati(ig)*radeg  ! latitude, in degrees
122   londeg=long(ig)*radeg  ! longitude, in degrees east
123   colat=90-latdeg        ! colatitude, in degrees
124
125 ! Find encompassing latitudes
126   if (colat<(90-lat(1))) then ! between north pole and lat(1)
127      ysup=1
128      yinf=1
129   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
130      ysup=latlen
131      yinf=latlen
132   else ! general case
133      do iloop=2,latlen
134         if(colat<(90-lat(iloop))) then
135           ysup=iloop-1
136           yinf=iloop
137           exit
138         endif
139      enddo
140   endif
141
142   latinf=lat(yinf)
143   latsup=lat(ysup)
144
145
146 ! Find encompassing longitudes
147   ! Note: in input file, lon(1)=-180.
148   if (londeg>lon(lonlen)) then
149      xsup=1
150      xinf=lonlen
151      loninf=lon(xsup)
152      lonsup=180.0 ! since lon(1)=-180.0
153   else
154      do iloop=2,lonlen
155         if(londeg<lon(iloop)) then
156              xsup=iloop
157              xinf=iloop-1
158              exit
159         endif
160      enddo
161      loninf=lon(xinf)
162      lonsup=lon(xsup)
163   endif
164
165   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
166     .OR.(ysup.lt.1)) then
167      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
168      write (*,*) "xinf: ",xinf
169      write (*,*) "xsup: ",xsup
170      write (*,*) "yinf: ",yinf
171      write (*,*) "ysup: ",ysup
172      stop
173   endif
174
175!   loninf=lon(xinf)
176!   lonsup=lon(xsup)
177!   latinf=lat(yinf)
178!   latsup=lat(ysup)
179
180! The four neighbouring points are arranged as follows:
181!                               1 2
182!                               3 4
183
184   latp(1)=ysup
185   latp(2)=ysup
186   latp(3)=yinf
187   latp(4)=yinf
188
189   lonp(1)=xinf
190   lonp(2)=xsup
191   lonp(3)=xinf
192   lonp(4)=xsup
193
194! Linear interpolation on time, for all four neighbouring points
195
196   if ((realday<time(1)).or.(realday>time(timelen))) then
197      tinf=timelen
198      tsup=1
199      timeflag=.true.
200   else
201      do iloop=2,timelen
202         if (realday<time(iloop)) then
203            tinf=iloop-1
204            tsup=iloop
205            exit
206         endif
207      enddo
208   endif
209
210! Bilinear interpolation on the four nearest points
211
212   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
213      dlat=0
214   else
215      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
216   endif
217
218   if (lonsup==loninf) then
219      dlon=0
220   else
221      dlon=(londeg-loninf)/(lonsup-loninf)
222   endif
223
224   do iloop=1,4
225      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
226      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
227      if (timeflag) then
228         if (realday>time(timelen)) then
229            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf)))
230         else
231            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
232         endif
233      else
234         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
235      endif
236      if (tau1(iloop)<0) then
237          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
238          write (*,*) "utime ",realday
239          write (*,*) "time(tinf) ",time(tinf)
240          write (*,*) "time(tsup) ",time(tsup)
241          write (*,*) "tau1 ",taubuf(1)
242          write (*,*) "tau2 ",taubuf(2)
243          write (*,*) "tau ",tau1(iloop)
244          stop
245      endif
246   enddo
247   timeflag=.false.
248
249   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
250      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
251      write (*,*) "dlat: ",dlat
252      write (*,*) "lat: ",latdeg
253      write (*,*) "dlon: ",dlon
254      write (*,*) "lon: ",londeg
255      stop
256   endif
257
258   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
259
260!   tauref(ig)=tau*700/pplev(ig,1)*0.8
261    tauref(ig)=tau*0.825
262enddo
263
264end
Note: See TracBrowser for help on using the repository browser.