source: trunk/MESOSCALE/LMD_MM_MARS/SRC/WPS/ungrib/src/output.F90

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

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 10.2 KB
Line 
1subroutine output(hdate, nlvl, maxlvl, plvl, interval, iflag, out_format, prefix, debug_level)
2!                                                                             !
3!*****************************************************************************!
4!  Write output to a file.
5!                                                                             !
6!    hdate:
7!    nlvl: 
8!    maxlvl:
9!    plvl:
10!    interval:
11!    iflag:  1 = output for ingest into rrpr ; 2 = output for hinterp
12!    out_format: requested output format (WPS, SI, or MM5)
13!                                                                             !
14!*****************************************************************************!
15
16  use table
17  use gridinfo
18  use storage_module
19  use filelist
20  use module_debug
21  use stringutil
22
23  implicit none
24
25  character(LEN=19) :: hdate
26  character(LEN=24) :: hdate_output
27  character(LEN=3)  :: out_format
28  character(LEN=256)  :: prefix
29  integer :: iunit = 13
30
31  real, pointer, dimension(:,:) :: scr2d
32
33  integer :: maxlvl
34  integer nlvl, debug_level
35  real , dimension(maxlvl) :: plvl
36  character (LEN=9) :: field
37  real :: level
38  integer :: sunit = 14
39  integer :: interval
40  integer :: iflag
41! Local Miscellaneous
42  integer :: k, n, m, ilev
43  integer :: ii, jj
44  real :: maxv, minv
45  real :: xplv
46  real :: xfcst = 0.
47  character (LEN=25) :: units
48  character (LEN=46) :: Desc
49  logical lopen
50
51! DATELEN:  length of date strings to use for our output file names.
52  integer :: datelen
53
54! Decide the length of date strings to use for output file names. 
55! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds.
56  if (mod(interval,3600) == 0) then
57     datelen = 13
58  elseif (mod(interval,60) == 0) then
59     datelen = 16
60  else
61     datelen = 19
62  endif
63 
64  call get_plvls(plvl, maxlvl, nlvl)
65
66  if ( debug_level .ge. 0 ) then
67  write(*,119) hdate(1:10), hdate(12:19)
68119 format(/,79('#'),//,'Inventory for date = ', A10,1x,A8,/)
69
70  write(*,advance='NO', fmt='("PRES", 2x)')
71  WRTLOOP : do n = 1, maxvar
72     do k = 1, n-1
73        if (namvar(k).eq.namvar(n)) cycle WRTLOOP
74     enddo
75     write(*,advance='NO', fmt='(1x,A8)') namvar(n)
76  enddo WRTLOOP
77  write(*,advance='YES', fmt='(1x)')
78
79  write(*,FMT='(79("-"))')
80  end if
81  KLOOP : do k = 1, nlvl
82     if ((iflag.eq.2).and.(plvl(k).gt.200100) .and. (plvl(k).lt.200200)) then
83        cycle KLOOP
84     endif
85     ilev = nint(plvl(k))
86     if ( debug_level .ge. 0 ) then
87     write(*, advance='NO', FMT='(F6.1)') plvl(k)/100.
88     end if
89     MLOOP : do m = 1, maxvar
90        do n = 1, m-1
91           if (namvar(m).eq.namvar(n)) cycle MLOOP
92        enddo
93        if ( debug_level .ge. 0 ) then
94        if (is_there(ilev,namvar(m))) then
95           write(*, advance='NO', FMT='("  X      ")')
96        else
97           if ( plvl(k).gt.200000 ) then
98             write(*, advance='NO', FMT='("  O      ")')
99           else
100             write(*, advance='NO', FMT='("         ")')
101           endif
102        endif
103        endif
104     enddo MLOOP
105     if ( debug_level .ge. 0 ) then
106     write(*,advance='YES', fmt='(1x)')
107     endif
108  enddo KLOOP
109  if ( debug_level .ge. 0 ) then
110  write(*,FMT='(79("-"))')
111  endif
112 
113  if (iflag.eq.1) then
114     if (nfiles.eq.0) then
115        open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
116             position='REWIND')
117        nfiles = nfiles + 1
118        filedates(nfiles)(1:datelen) = hdate(1:datelen)
119     else
120        DOFILES : do k = 1, nfiles
121           if (hdate(1:datelen).eq.filedates(k)(1:datelen)) then
122              open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted',&
123                   position='APPEND')
124           endif
125        enddo DOFILES
126        inquire (iunit, OPENED=LOPEN)
127        if (.not. LOPEN) then
128           open(iunit, file=trim(get_path(prefix))//'PFILE:'//HDATE(1:datelen), form='unformatted', &
129                position='REWIND')
130           nfiles = nfiles + 1
131           filedates(nfiles)(1:datelen) = hdate(1:datelen)
132        endif
133     endif
134  else if (iflag.eq.2) then
135     open(iunit, file=trim(prefix)//':'//HDATE(1:datelen), form='unformatted', &
136          position='REWIND')
137  endif
138
139!MGD  if ( debug_level .gt. 100 ) then
140!MGD     write(6,*) 'begin nloop'
141!MGD  end if
142  NLOOP : do n = 1, nlvl
143
144!MGD  if ( debug_level .gt. 100 ) then
145!MGD     write(6,*) 'begin outloop'
146!MGD  end if
147     OUTLOOP : do m = 1, maxvar
148        field = namvar(m)
149        do k = 1, m-1
150           if (field.eq.namvar(k)) cycle OUTLOOP
151        enddo
152        level = plvl(n)
153        if ((iflag.eq.2).and.(level.gt.200100) .and. (level.lt.200200)) then
154           cycle NLOOP
155        endif
156        ilev = nint(level)
157        desc = ddesc(m)
158        if (iflag.eq.2) then
159           if (desc.eq.' ') cycle OUTLOOP
160        endif
161        units = dunits(m)
162        if ((iflag.eq.1).or.(iflag.eq.2.and.desc(1:1).ne.' ')) then
163          if (is_there(ilev,field)) then
164            call get_dims(ilev, field)
165
166!MGD            if ( debug_level .gt. 100 ) then
167!MGD               write(6,*) 'call refr_storage'
168!MGD            end if
169            call refr_storage(ilev, field, scr2d, map%nx, map%ny)
170
171!MGD            if ( debug_level .gt. 100 ) then
172!MGD               write(6,*) 'back from refr'
173!MGD               write(6,*) 'out_format = ',out_format
174!MGD            end if
175
176            if (out_format(1:2) .eq. 'SI') then
177!MGD              if ( debug_level .gt. 100 ) then
178!MGD                 write(6,*) 'writing in SI format'
179!MGD              end if
180              write(iunit) 4
181              hdate_output = hdate
182              write (iunit) hdate_output, xfcst, map%source, field, units, &
183                   Desc, level, map%nx, map%ny, map%igrid
184              if (map%igrid.eq.3) then ! lamcon
185                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
186                      map%lov, map%truelat1, map%truelat2
187              elseif (map%igrid.eq.5) then ! Polar Stereographic
188                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
189                      map%lov, map%truelat1
190              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
191                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx
192              elseif (map%igrid.eq.1)then ! Mercator
193                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
194                      map%truelat1
195              else
196                 call mprintf(.true.,ERROR, &
197                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
198              endif
199              write (iunit) scr2d
200            else if (out_format(1:2) .eq. 'WP') then   
201                call mprintf(.true.,DEBUG, &
202         "writing in WPS format  iunit = %i, map%%igrid = %i",i1=iunit,i2=map%igrid)
203              write(iunit) 5
204              hdate_output = hdate
205              write (iunit) hdate_output, xfcst, map%source, field, units, &
206                   Desc, level, map%nx, map%ny, map%igrid
207              if (map%igrid.eq.3) then ! lamcon
208                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
209                      map%lov, map%truelat1, map%truelat2, map%r_earth
210              elseif (map%igrid.eq.5) then ! Polar Stereographic
211                 write (iunit) map%startloc, map%lat1, map%lon1, map%dx, map%dy, &
212                      map%lov, map%truelat1, map%r_earth
213              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
214                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
215                      map%r_earth
216              elseif (map%igrid.eq.1)then ! Mercator
217                 write (iunit) map%startloc, map%lat1, map%lon1, map%dy, map%dx, &
218                      map%truelat1, map%r_earth
219              else
220                 call mprintf(.true.,ERROR, &
221                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
222              endif
223              write (iunit) map%grid_wind
224              write (iunit) scr2d
225            else if (out_format(1:2) .eq. 'MM') then
226!MGD              if ( debug_level .gt. 100 ) then
227!MGD             write(6,*) 'writing in MM5 format'
228!MGD              end if
229              write(iunit) 3
230              hdate_output = hdate
231              write (iunit) hdate_output, xfcst, field, units, Desc, level,&
232                   map%nx, map%ny, map%igrid
233              if (map%igrid.eq.3) then ! lamcon
234                 write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
235                      map%truelat1, map%truelat2
236              elseif (map%igrid.eq.5) then ! Polar Stereographic
237                 write (iunit) map%lat1, map%lon1, map%dx, map%dy, map%lov, &
238                      map%truelat1
239              elseif (map%igrid.eq.0 .or. map%igrid.eq.4)then ! lat/lon
240                 write (iunit) map%lat1, map%lon1, map%dy, map%dx
241              elseif (map%igrid.eq.1)then ! Mercator
242                 write (iunit) map%lat1, map%lon1, map%dy, map%dx, map%truelat1
243              else
244                 call mprintf(.true.,ERROR, &
245                "Unrecognized map%%igrid: %i in subroutine output 1",i1=map%igrid)
246              endif
247              write (iunit) scr2d
248            endif
249              if ( debug_level .gt. 100 ) then
250                call mprintf(.true.,DEBUG, &
251                "hdate = %s,  xfcst = %f ",s1=hdate_output,f1=xfcst)
252                call mprintf(.true.,DEBUG, &
253           "map%%source = %s, field = %s, units = %s",s1=map%source,s2=field,s3=units)
254                call mprintf(.true.,DEBUG, &
255                    "Desc = %s, level = %f",s1=Desc,f1=level)
256                call mprintf(.true.,DEBUG, &
257                    "map%%nx = %i, map%%ny = %i",i1=map%nx,i2=map%ny)
258              else if ( debug_level .gt. 0 ) then
259                call mprintf(.true.,STDOUT, &
260              " field = %s, level = %f",s1=field,f1=level)
261                call mprintf(.true.,LOGFILE, &
262              " field = %s, level = %f",s1=field,f1=level)
263              end if
264              if ( debug_level .gt. 100 ) then
265              maxv = -99999.
266              minv = 999999.
267              do jj = 1, map%ny
268              do ii = 1, map%nx
269                if (scr2d(ii,jj) .gt. maxv) maxv = scr2d(ii,jj)
270                if (scr2d(ii,jj) .lt. minv) minv = scr2d(ii,jj)
271              enddo
272              enddo
273              call mprintf(.true.,DEBUG, &
274                 "max value = %f , min value = %f",f1=maxv,f2=minv)
275              end if
276
277              nullify(scr2d)
278
279           endif
280        endif
281     enddo OUTLOOP
282  enddo NLOOP
283
284  close(iunit)
285
286end subroutine output
287
Note: See TracBrowser for help on using the repository browser.