source: trunk/MESOSCALE_DEV/SRC/ARWpost/src/v5d_module.f90 @ 897

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

MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
File size: 9.8 KB
RevLine 
[207]1MODULE v5d_module
2
3   USE gridinfo_module
4   USE output_module
5   USE map_utils
6   USE misc_definitions_module
7
8
9   CONTAINS
10
11   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12   ! Name: v5d_varnames     
13   ! Purpose: Get info about the variables to send to v5d file
14   !          This routine assumes that the user only asked for variables
15   !          that are available. If variables are not available v5d will no
16   !          be working correctly.
17   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
18   SUBROUTINE v5d_varnames (v5d_nl)
19
20     IMPLICIT NONE
21
22     ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS
23     include 'v5df.h'
24
25     integer, dimension(MAXLEVELS)          :: v5d_nl
26     character (len=4000)                   :: v5d_fields
27     integer                                :: is_there, is_3d
28     character (len=1000)                   :: known_3d_fields
29     character (len=20)                     :: dummy
30
31     
32     known_3d_fields = 'P,PB,QCLOUD,QGRAUP,QICE,QNICE,QRAIN,QSNOW, &
33&QVAPOR,T,TKE,TKE_MYJ,U,V,PH,PHB,W,cape,cin,dbz,height,pressure,rh, &
34&tc,td,theta,tk,umet,vmet,wdir,wspd'
35
36
37     varnames = ' '
38     numvars = 0
39     v5d_fields = plot_these_fields(2:len_trim(plot_these_fields))
40     is_there = INDEX(v5d_fields,",")
41     DO WHILE ( is_there > 1 )
42
43       numvars = numvars  + 1
44       varnames(numvars) = v5d_fields(1:is_there-1)
45       dummy = ","//trim(varnames(numvars))//","
46
47       IF ( trim(varnames(numvars)) == 'SH2O'  .OR.  &
48            trim(varnames(numvars)) == 'SMOIS' .OR.  &
49            trim(varnames(numvars)) == 'TSLB'  ) THEN
50         print*,"  WARING:  We cannot deal with field ", trim(varnames(numvars))
51         print*,"           It will be removed from the plot list"
52         is_3d = INDEX(plot_these_fields,trim(dummy))
53         plot_these_fields = plot_these_fields(1:is_3d)//plot_these_fields(is_3d+len_trim(dummy):len_trim(plot_these_fields))
54           
55         varnames(numvars) = ' '
56         numvars = numvars - 1
57       ELSE
58         v5d_nl(numvars) = 1
59         is_3d = INDEX(known_3d_fields,trim(varnames(numvars)))
60         IF ( is_3d > 0 ) v5d_nl(numvars) = number_of_zlevs
61       END IF 
62
63       v5d_fields = v5d_fields(is_there+1:len_trim(v5d_fields))
64
65       !!! Check for doubles
66       is_there = INDEX(v5d_fields,trim(dummy))
67       DO WHILE ( is_there /= 0 )
68          v5d_fields = v5d_fields(1:is_there)//v5d_fields(is_there+len_trim(dummy):len_trim(v5d_fields))
69          is_there = INDEX(v5d_fields,trim(dummy))
70        END DO
71
72       is_there = INDEX(v5d_fields,",")
73
74     END DO
75
76 
77   END SUBROUTINE v5d_varnames
78
79   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80   ! Name: v5d_times     
81   ! Purpose: Get info about the variables to send to v5d file
82   !          This routine assumes that the user only asked for variables
83   !          that are available. If variables are not available v5d will no
84   !          be working correctly.
85   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86   SUBROUTINE v5d_times (n_times, timestamp, datestamp)
87
88     USE date_pack
89
90     IMPLICIT NONE
91
92     ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS
93     include 'v5df.h'
94
95     integer, dimension(MAXTIMES)           :: timestamp, datestamp
96     character (len=19)                     :: next_date
97     integer                                :: tt, n_times
98     integer                                :: year, month, day, hours, minutes, seconds
99
100
101     datestamp = IMISSING
102     timestamp = IMISSING
103
104     DO tt = 0,n_times
105       CALL geth_newdate(next_date, trim(start_date), tt*interval_seconds)
106
107       read(next_date(1:4),*)   year
108       read(next_date(6:7),*)   month
109       read(next_date(9:10),*)  day
110       read(next_date(12:13),*) hours
111       read(next_date(15:16),*) minutes
112       read(next_date(18:19),*) seconds
113
114
115       timestamp(tt) = hours*10000 + minutes*100 + seconds
116
117       if(month >= 2) day = day+31  ! add january
118       if(month >= 3) day = day+28  ! add february
119       if(month >= 4) day = day+31  ! add march
120       if(month >= 5) day = day+30  ! add april
121       if(month >= 6) day = day+31  ! add may
122       if(month >= 7) day = day+30  ! add june
123       if(month >= 8) day = day+31  ! add july
124       if(month >= 9) day = day+31  ! add august
125       if(month >= 10) day = day+30 ! add september
126       if(month >= 11) day = day+31 ! add october
127       if(month >= 12) day = day+30 ! add november
128       if((month > 2) .and. (mod(year,4) == 0)) day = day+1  ! get leap year day
129
130       datestamp(tt) = (year)*1000 + day
131
132     END DO
133
134   END SUBROUTINE v5d_times
135
136   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
137   ! Name: v5d_proj     
138   ! Purpose: Get map information needed to create v5d output file
139   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140   SUBROUTINE v5d_proj (date, projection, proj_args)
141
142     USE input_module
143
144     IMPLICIT NONE
145
146     character (len=19)                  :: date
147     integer                             :: projection
148     real, dimension(100)                :: proj_args
149     real                                :: xlat_s, xlat_n, xlon_w, xlon_e, xlat1, xlon1
150     real                                :: xx, yy, temp
151
152     integer                             :: istatus, wrftype
153     real, pointer, dimension(:,:,:)     :: real_array
154     character (len=3)                   :: memorder
155     character (len=128)                 :: cname, stagger
156     character (len=128), dimension(3)   :: dimnames
157     integer, dimension(3)               :: domain_start, domain_end
158     type (proj_info)                    :: projst
159     real                                :: msf
160     real                                :: hemiflag
161     real                                :: pi, degtorad
162
163
164
165     domain_start = 1
166     domain_end(1) = west_east_dim
167     domain_end(2) = south_north_dim
168     domain_end(3) = 1
169     memorder = "XY"
170     stagger = '-'
171     wrftype = 104
172
173     istatus = 0
174     IF ( iprogram .ge. 6 ) THEN
175       CALL read_spec_field(domain_start, domain_end, 'XLAT', wrftype, &
176                            memorder, stagger, dimnames, real_array, date, istatus)
177     ELSE                 
178       CALL read_spec_field(domain_start, domain_end, 'XLAT_M', wrftype, &
179                            memorder, stagger, dimnames, real_array, date, istatus)
180     END IF
181     IF ( istatus == 0) THEN
182       xlat_s = real_array(1,1,1)
183       xlat_n = real_array(1,south_north_dim,1)
184       xlat1  = real_array(1,1,1)
185     END IF
186     istatus = 0
187     IF ( iprogram .ge. 6 ) THEN
188       CALL read_spec_field(domain_start, domain_end, 'XLONG', wrftype, &
189                            memorder, stagger, dimnames, real_array, date, istatus)
190     ELSE                 
191       CALL read_spec_field(domain_start, domain_end, 'XLONG_M', wrftype, &
192                            memorder, stagger, dimnames, real_array, date, istatus)
193     END IF
194     IF ( istatus == 0) THEN
195       xlon_w = real_array(1,south_north_dim,1)
196       xlon_e = real_array(west_east_dim,south_north_dim,1)
197       xlon1  = real_array(1,1,1)
198     END IF
199
200
201
202     IF ( MAP_PROJ == 0 ) THEN        !!! ideal cases
203
204       projection = 0
205       proj_args(1) = 1.
206       proj_args(2) = 1.
207       proj_args(3) = 1.
208       proj_args(4) = 1.
209
210     ELSE IF ( MAP_PROJ == 1 ) then    !!! Lambert
211
212       projection = 2
213
214       IF ( TRUELAT1 < TRUELAT2 ) THEN   !!! truelat1 must be bigger than truelat2
215         temp = TRUELAT1
216         TRUELAT1 = TRUELAT2
217         TRUELAT2 = temp
218       END IF
219
220       CALL make_lambert(xlat_n,xlon_w,xx,yy)
221
222       proj_args(1) = TRUELAT1
223       proj_args(2) = TRUELAT2
224       proj_args(3) = -yy
225       proj_args(4) = -xx
226       proj_args(5) = -STAND_LON
227       proj_args(6) = 0.001*DX
228
229     ELSE IF ( MAP_PROJ == 2 ) THEN   !!! PS
230
231       CALL map_set(PROJ_PS,xlat1,xlon1,DX,STAND_LON,TRUELAT1,TRUELAT2, &
232                    west_east_dim, south_north_dim, projst)
233       projection = 3
234
235       IF (TRUELAT1 .ge. 0.) THEN
236         proj_args(1) = 90.
237         hemiflag = 1.
238       ELSE
239         proj_args(1) = -90.
240         hemiflag = -1.
241       END IF
242       proj_args(2) = -STAND_LON
243       proj_args(3) = float(south_north_dim) - projst%polej
244       proj_args(4) = projst%polei
245       msf  = (1.+hemiflag*sin(projst%TRUELAT1*degtorad))/ &
246              (1.+hemiflag*sin(proj_args(1)*degtorad))
247       proj_args(5) = DX* 0.001/msf
248 
249     ELSE IF ( MAP_PROJ == 3 ) THEN   !!! Mercator
250
251       projection = 1
252
253       proj_args(1) = xlat_n
254       proj_args(2) = -xlon_w
255       proj_args(3) = (xlon_w - xlon_e) / (west_east_dim - 1)
256       proj_args(4) = (xlat_n - xlat_s) / (south_north_dim - 1)
257
258     END IF
259 
260   END SUBROUTINE v5d_proj
261
262
263  SUBROUTINE make_lambert(xlatn,xlonw,xx,yy)
264
265!  imported from TOVIS5D for MM5
266
267     USE input_module
268
269     real                 :: xlatn, xlonw, xx, yy
270     real                 :: radius, degtorad, alpha
271     real                 :: theta, theta1, theta2, htheta, htheta1, htheta2, cone, dis, rad
272
273     data radius /6371./
274
275!  put map parameters in as defined by vis5d;
276!    south is positive direction, so negative ypole
277
278     pi       = 4.0*atan( 1.0 )
279     degtorad = pi/180.0
280 
281     dis      = 0.001*DX
282 
283     theta    = (90.0-xlatn) * degtorad
284     theta1   = (90.0-TRUELAT1) * degtorad
285     theta2   = (90.0-TRUELAT2) * degtorad
286     htheta   = theta /2.0
287     htheta1  = theta1/2.0
288     htheta2  = theta2/2.0
289     cone     = alog(sin(theta1)/sin(theta2))/alog(tan(htheta1)/tan(htheta2))
290
291     alpha    = cone * (xlonw - STAND_LON) * degtorad
292     rad      = radius * sin(theta1) * (tan(htheta)/tan(htheta1))**cone/cone
293
294     xx       = rad * sin(alpha) / dis
295     yy       = rad * cos(alpha) / dis
296
297  END SUBROUTINE make_lambert
298
299
300END MODULE v5d_module
Note: See TracBrowser for help on using the repository browser.