source: trunk/mesoscale/COMMON_GCM/readtesassim.F90 @ 42

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

LMD_MM_MARS: tests avec physique completement commune avec le GCM nouvelle physique
File size: 8.1 KB
Line 
1subroutine readtesassim(ngrid,nlayer,zday,pplev,tauref)
2
3! Reading of the dust assimilation file
4
5implicit none
6
7#include "dimensions.h"
8#include "dimphys.h"
9#include "comgeomfi.h"
10#include "netcdf.inc"
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   else if ((iaervar.eq.6).or.(iaervar.eq.7)) then ! modified TES scenario
63     filename="modifieddust.nc"
64   endif
65   
66   ! Note: datafile() is defined in "datafile.h"
67   ierr=NF_OPEN(trim(datafile)//"/"//trim(filename),NF_NOWRITE,nid)
68   IF (ierr.NE.NF_NOERR) THEN
69     write(*,*)'Problem opening ',trim(filename),' (in phymars/readtesassim.F90)'
70     write(*,*)'It should be in :',trim(datafile),'/'
71     write(*,*)'1) You can change this directory address in '
72     write(*,*)'   file phymars/datafile.h'
73     write(*,*)'2) If necessary, dust.nc (and other datafiles)'
74     write(*,*)'   can be obtained online on:'
75     write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
76     CALL ABORT
77   ENDIF
78
79   ierr=NF_INQ_DIMID(nid,"Time",timedim)
80   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
81   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
82   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
83   ierr=NF_INQ_DIMID(nid,"longitude",londim)
84   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
85
86   allocate(tautes(lonlen,latlen,timelen))
87   allocate(lat(latlen), lon(lonlen), time(timelen))
88
89   ierr = NF_INQ_VARID (nid, "dustop",nvarid)
90#ifdef NC_DOUBLE
91   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tautes)
92#else
93   ierr = NF_GET_VAR_REAL(nid, nvarid, tautes)
94#endif
95   IF (ierr .NE. NF_NOERR) THEN
96      PRINT*, "Error: Readtesassim <dustop> not found"
97      stop
98   ENDIF
99   ierr = NF_INQ_VARID (nid, "Time",nvarid)
100#ifdef NC_DOUBLE
101   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
102#else
103   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
104#endif
105   IF (ierr .NE. NF_NOERR) THEN
106      PRINT*, "Error: Readtesassim <Time> not found"
107      stop
108   ENDIF
109
110   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
111#ifdef NC_DOUBLE
112   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
113#else
114   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
115#endif
116   IF (ierr .NE. NF_NOERR) THEN
117      PRINT*, "Error: Readtesassim <latitude> not found"
118      stop
119   ENDIF
120
121   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
122#ifdef NC_DOUBLE
123   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
124#else
125   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
126#endif
127   IF (ierr .NE. NF_NOERR) THEN
128      PRINT*, "Error: Readtesassim <longitude> not found"
129      stop
130   ENDIF
131
132endif ! of if (firstcall)
133
134do ig=1,ngrid
135
136! Find the four nearest points, arranged as follows:
137!                               1 2
138!                               3 4
139
140   latdeg=lati(ig)*radeg  ! latitude, in degrees
141   londeg=long(ig)*radeg  ! longitude, in degrees east
142   colat=90-latdeg        ! colatitude, in degrees
143
144 ! Find encompassing latitudes
145   if (colat<(90-lat(1))) then ! between north pole and lat(1)
146      ysup=1
147      yinf=1
148   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
149      ysup=latlen
150      yinf=latlen
151   else ! general case
152      do iloop=2,latlen
153         if(colat<(90-lat(iloop))) then
154           ysup=iloop-1
155           yinf=iloop
156           exit
157         endif
158      enddo
159   endif
160
161   latinf=lat(yinf)
162   latsup=lat(ysup)
163
164
165 ! Find encompassing longitudes
166   ! Note: in input file, lon(1)=-180.
167   if (londeg>lon(lonlen)) then
168      xsup=1
169      xinf=lonlen
170      loninf=lon(xsup)
171      lonsup=180.0 ! since lon(1)=-180.0
172   else
173      do iloop=2,lonlen
174         if(londeg<lon(iloop)) then
175              xsup=iloop
176              xinf=iloop-1
177              exit
178         endif
179      enddo
180      loninf=lon(xinf)
181      lonsup=lon(xsup)
182   endif
183
184   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
185     .OR.(ysup.lt.1)) then
186      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
187      write (*,*) "xinf: ",xinf
188      write (*,*) "xsup: ",xsup
189      write (*,*) "yinf: ",yinf
190      write (*,*) "ysup: ",ysup
191      stop
192   endif
193
194!   loninf=lon(xinf)
195!   lonsup=lon(xsup)
196!   latinf=lat(yinf)
197!   latsup=lat(ysup)
198
199! The four neighbouring points are arranged as follows:
200!                               1 2
201!                               3 4
202
203   latp(1)=ysup
204   latp(2)=ysup
205   latp(3)=yinf
206   latp(4)=yinf
207
208   lonp(1)=xinf
209   lonp(2)=xsup
210   lonp(3)=xinf
211   lonp(4)=xsup
212
213! Linear interpolation on time, for all four neighbouring points
214
215   if ((realday<time(1)).or.(realday>time(timelen))) then
216      tinf=timelen
217      tsup=1
218      timeflag=.true.
219   else
220      do iloop=2,timelen
221         if (realday<time(iloop)) then
222            tinf=iloop-1
223            tsup=iloop
224            exit
225         endif
226      enddo
227   endif
228
229! Bilinear interpolation on the four nearest points
230
231   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
232      dlat=0
233   else
234      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
235   endif
236
237   if (lonsup==loninf) then
238      dlon=0
239   else
240      dlon=(londeg-loninf)/(lonsup-loninf)
241   endif
242
243   do iloop=1,4
244      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
245      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
246      if (timeflag) then
247         if (realday>time(timelen)) then
248            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf)))
249         else
250            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
251         endif
252      else
253         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
254      endif
255      if (tau1(iloop)<0) then
256          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
257          write (*,*) "utime ",realday
258          write (*,*) "time(tinf) ",time(tinf)
259          write (*,*) "time(tsup) ",time(tsup)
260          write (*,*) "tau1 ",taubuf(1)
261          write (*,*) "tau2 ",taubuf(2)
262          write (*,*) "tau ",tau1(iloop)
263          stop
264      endif
265   enddo
266   timeflag=.false.
267
268   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
269      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
270      write (*,*) "dlat: ",dlat
271      write (*,*) "lat: ",latdeg
272      write (*,*) "dlon: ",dlon
273      write (*,*) "lon: ",londeg
274      stop
275   endif
276
277   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
278
279! Multiply by correction coefficient (value depends on SW radiative transfer)
280! NB: swrtype is set in callkeys.h
281   if (swrtype.eq.1) then ! Fouquart
282     tauref(ig)=tau*0.825
283   else
284     if (swrtype.eq.2) then ! Toon
285       tauref(ig)=tau*1.3
286        ! special case for enhanced flushing storm scenario
287        if ((iaervar.eq.7).and.(tauref(ig).gt.1.1*1.3)) then
288          tauref(ig)=2.0*tauref(ig)
289        endif
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.