source: trunk/MESOSCALE/LMD_MM_MARS/SRC/WPS/ungrib/src/datint.F90 @ 3992

Last change on this file since 3992 was 11, checked in by aslmd, 15 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 13.7 KB
Line 
1subroutine datint(fuldates, nful, hstart, ntimes, interval, out_format, prefix)
2!                                                                             !
3!*****************************************************************************!
4!                                                                             !
5!   interpolate missing data in time
6!    out_format: requested output format
7!                                                                             !
8!*****************************************************************************!
9  use gridinfo
10  use storage_module
11  use module_debug
12  implicit none
13  integer :: nful
14  integer :: interval
15  character(len=*), dimension(nful) :: fuldates
16  character(len=*) :: hstart
17  integer :: ntimes
18
19  character(len=24) :: hdate = "0000-00-00_00:00:00.0000"
20  character(len=24) :: hdate_output, jdate
21  character(len=9) :: field
22  character(len=25) :: units
23  character(len=46) :: desc
24  character(LEN=3)  :: out_format
25  character(LEN=256)  :: prefix
26  real :: xfcst
27
28  real :: level
29  real, allocatable, dimension(:,:) :: scr2d, bfr2d
30  integer :: iful, intervala, intervalb, ifv
31  real :: awt
32  integer :: itime
33
34! DATELEN:  length of date strings to use for our output file names.
35  integer :: datelen
36
37! Decide the length of date strings to use for output file names. 
38! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
39  if (mod(interval,3600) == 0) then
40     datelen = 13
41  else if (mod(interval, 60) == 0) then
42     datelen = 16
43  else
44     datelen = 19
45  end if
46
47  call mprintf(.true.,STDOUT,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...")
48  call mprintf(.true.,LOGFILE,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...")
49
50  TIMELOOP : do itime = 1, ntimes
51     call geth_newdate(hdate(1:19), hstart(1:19), (itime-1)*interval)
52     call mprintf(.true.,STDOUT,"Looking for data at time %s",s1=hdate(1:datelen))
53     call mprintf(.true.,LOGFILE,"Looking for data at time %s",s1=hdate(1:datelen))
54     do iful = 1, nful
55        if (fuldates(iful).eq.hdate) then
56           call mprintf(.true.,STDOUT,"Found file:      %s:%s", &
57                 s1=trim(prefix),s2=hdate(1:datelen))
58           call mprintf(.true.,LOGFILE,"Found file:      %s:%s", &
59                 s1=trim(prefix),s2=hdate(1:datelen))
60           cycle TIMELOOP
61        else if ((fuldates(iful).lt.hdate) .and. &
62             (fuldates(iful+1).gt.hdate) )then
63
64           call mprintf(.true.,STDOUT,"Found surrounding files:      %s: %s  and %s: %s", &
65               s1=trim(prefix),s2=fuldates(iful)(1:datelen), &
66               s3=trim(prefix),s4=fuldates(iful+1)(1:datelen))
67           call mprintf(.true.,LOGFILE,"Found surrounding files:      %s: %s  and %s: %s", &
68               s1=trim(prefix),s2=fuldates(iful)(1:datelen), &
69               s3=trim(prefix),s4=fuldates(iful+1)(1:datelen))
70           call mprintf(.true.,STDOUT,"Interpolating to create file:      %s: %s", &
71                 s1=trim(prefix),s2=hdate(1:datelen))
72           call mprintf(.true.,LOGFILE,"Interpolating to create file:      %s: %s", &
73                 s1=trim(prefix),s2=hdate(1:datelen))
74           call geth_idts(hdate(1:19), fuldates(iful)(1:19), intervalA)
75           call mprintf(.true.,STDOUT,"A Time Difference = %f",f1=float(intervalA) / 3600.)
76           call mprintf(.true.,LOGFILE,"A Time Difference = %f",f1=float(intervalA) / 3600.)
77           call geth_idts(fuldates(iful+1)(1:19), hdate(1:19), intervalB)
78           call mprintf(.true.,STDOUT,"B Time Difference = %f",f1=float(intervalB) / 3600.)
79           call mprintf(.true.,LOGFILE,"B Time Difference = %f",f1=float(intervalB) / 3600.)
80           AWT = 1. - (float(intervalA)/float(intervalA+intervalB))
81
82           open(10, file=trim(prefix)//':'//fuldates(iful)(1:datelen), form='unformatted', &
83                status='old')
84           call clear_storage
85           READLOOP1 : do
86              read(10, end=44) ifv
87              if ( ifv .eq. 5) then     ! WPS
88                read (10) jdate, xfcst, map%source, field, units, desc, level, &
89                     map%nx, map%ny, map%igrid
90                select case (map%igrid)
91                case (0, 4)
92                   read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
93                case (3)
94                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
95                        map%lov, map%truelat1, map%truelat2, map%r_earth
96                case (5)
97                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
98                        map%lov, map%truelat1, map%r_earth
99                case default
100                  call mprintf(.true.,ERROR, &
101                    "Unrecognized map%%igrid: %i in DATINT 1",i1=map%igrid)
102                end select
103                read (10) map%grid_wind
104              else if ( ifv .eq. 4 ) then          ! SI
105                read (10) jdate, xfcst, map%source, field, units, desc, level, &
106                      map%nx, map%ny, map%igrid
107                select case (map%igrid)
108                case (0, 4)
109                   read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx
110                case (3)
111                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
112                        map%lov, map%truelat1, map%truelat2
113                case (5)
114                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
115                        map%lov, map%truelat1
116                case default
117                  call mprintf(.true.,ERROR, &
118                    "Unrecognized map%%igrid: %i in DATINT 2",i1=map%igrid)
119                end select
120              else if ( ifv .eq. 3 ) then          ! MM5
121                read(10) jdate, xfcst, field, units, desc, level,&
122                     map%nx, map%ny, map%igrid
123                select case (map%igrid)
124                case (3)      ! lamcon
125                   read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
126                        map%truelat1, map%truelat2
127                case (5)      ! Polar Stereographic
128                   read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
129                        map%truelat1
130                case (0, 4)      ! lat/lon
131                   read (10) map%lat1, map%lon1, map%dy, map%dx
132                case (1)      ! Mercator
133                   read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
134                case default
135                  call mprintf(.true.,ERROR, &
136                    "Unrecognized map%%igrid: %i in DATINT 3",i1=map%igrid)
137                end select
138              else
139                call mprintf(.true.,ERROR, &
140                  "Unknown out_format: %i in DATINT ",i1=ifv)
141              endif
142              allocate(scr2d(map%nx, map%ny))
143              read (10) scr2d
144              call put_storage(nint(level), field, scr2d, map%nx, map%ny)
145              deallocate(scr2d)
146           enddo READLOOP1
14744         close(10)
148
149           open(10, file=trim(prefix)//':'//fuldates(iful+1)(1:datelen), status='old', &
150                form = 'unformatted')
151           open(11, file=trim(prefix)//':'//hdate(1:datelen), status='new', form='unformatted')
152           READLOOP2 : do
153              read (10,END=45) ifv
154              if ( ifv .eq. 5) then     ! WPS
155                read (10) jdate, xfcst, map%source, field, units, desc, level, &
156                      map%nx, map%ny, map%igrid
157                select case (map%igrid)
158                case (0, 4)
159                   read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth
160                case (3)
161                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
162                        map%lov, map%truelat1, map%truelat2, map%r_earth
163                case (5)
164                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
165                      map%lov, map%truelat1, map%r_earth
166                case default
167                   call mprintf(.true.,ERROR, &
168                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
169                end select
170                read (10) map%grid_wind
171              else if ( ifv .eq. 4 ) then          ! SI
172                read (10) jdate, xfcst, map%source, field, units, desc, level, &
173                      map%nx, map%ny, map%igrid
174                select case (map%igrid)
175                case (0, 4)
176                   read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx
177                case (1)
178                   read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1
179                case (3)
180                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
181                      map%lov, map%truelat1, map%truelat2
182                case (5)
183                   read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
184                      map%lov, map%truelat1
185                case default
186                   call mprintf(.true.,ERROR, &
187                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
188                end select
189
190              else if ( ifv .eq. 3 ) then          ! MM5
191                read(10) jdate, xfcst, field, units, desc, level,&
192                     map%nx, map%ny, map%igrid
193                select case (map%igrid)
194                case (3)    ! lamcon
195                   read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
196                        map%truelat1, map%truelat2
197                case (5)    ! Polar Stereographic
198                   read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
199                      map%truelat1
200                case (0, 4)    ! lat/lon
201                   read (10) map%lat1, map%lon1, map%dy, map%dx
202                case (1)    ! Mercator
203                   read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
204                case default
205                   call mprintf(.true.,ERROR, &
206                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
207                end select
208
209              else
210                call mprintf(.true.,ERROR, &
211                  "Unknown out_format: %i in DATINT ",i1=ifv)
212              endif
213
214              allocate(scr2d(map%nx, map%ny))
215              read (10) scr2d
216              if (is_there(nint(level), field)) then
217                 allocate(bfr2d(map%nx,map%ny))
218                 call get_storage(nint(level), field, bfr2d, map%nx, map%ny)
219                 scr2d = bfr2d * (AWT) + scr2d * (1.-AWT)
220                 hdate_output = hdate
221
222                 if (out_format(1:2) .eq. 'SI') then
223                 write(11) ifv
224                 write(11) hdate_output, xfcst, map%source, field, units, desc, &
225                      level, map%nx, map%ny, map%igrid
226                 if (map%igrid == 0 .or. map%igrid == 4) then
227                    write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx
228                 elseif (map%igrid == 1) then
229                    write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1
230                 elseif (map%igrid == 3) then
231                    write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
232                         map%lov, map%truelat1, map%truelat2
233                 elseif (map%igrid == 5) then
234                    write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
235                         map%lov, map%truelat1
236                 else
237                   call mprintf(.true.,ERROR, &
238                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
239                 endif
240
241                 else if (out_format(1:2) .eq. 'WP') then
242                 write(11) ifv
243                 write(11) hdate_output, xfcst, map%source, field, units, desc, &
244                      level, map%nx, map%ny, map%igrid
245                 if (map%igrid == 0 .or. map%igrid == 4) then
246                    write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
247                              map%r_earth
248                 elseif (map%igrid == 1) then
249                    write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
250                              map%truelat1, map%r_earth
251                 elseif (map%igrid == 3) then
252                    write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
253                         map%lov, map%truelat1, map%truelat2, map%r_earth
254                 elseif (map%igrid == 5) then
255                    write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
256                         map%lov, map%truelat1, map%r_earth
257                 else
258                   call mprintf(.true.,ERROR, &
259                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
260                 endif
261                 write(11) map%grid_wind
262
263                 else if (out_format(1:2) .eq. 'MM') then
264                 write (11) ifv
265                 write (11) hdate_output, xfcst, field, units, Desc, level,&
266                   map%nx, map%ny, map%igrid
267                 if (map%igrid .eq. 3) then ! lamcon
268                   write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
269                      map%truelat1, map%truelat2
270                 elseif (map%igrid .eq. 5) then ! Polar Stereographic
271                   write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
272                      map%truelat1
273                 elseif (map%igrid .eq. 0 .or. map%igrid .eq. 4)then ! lat/lon
274                   write (11) map%lat1, map%lon1, map%dy, map%dx
275                 elseif (map%igrid.eq.1)then ! Mercator
276                   write (11) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
277                 else
278                   call mprintf(.true.,ERROR, &
279                     "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid)
280                 endif
281                 endif
282                 write(11) scr2d
283              else
284                 call mprintf(.true.,ERROR, &
285                  "hdate = %s , fuldates = %s %s, Field = %s",s1=hdate,s2=fuldates(iful),s3=fuldates(iful+1),s4=field)
286              endif
287              deallocate(scr2d, bfr2d)
288           enddo READLOOP2
28945         close(10)
290           close(11)
291           cycle TIMELOOP
292        endif
293     enddo
294
295     call mprintf(.true.,ERROR, &
296        "Data not found: %s",s1=hdate)
297     
298  enddo TIMELOOP
299
300  call mprintf(.true.,STDOUT, &
301   "End Subroutine DATINT.")
302  call mprintf(.true.,LOGFILE, &
303   "End Subroutine DATINT.")
304
305end subroutine datint
Note: See TracBrowser for help on using the repository browser.