W2GCONV.inc
References to this file elsewhere.
1 SUBROUTINE W2GCONV(IGRD, JGRD, KGRD, JCAP, KMAX, INTVL ,
2 C DPHI, LMAX, KLMAX, MEND1, NEND1, JEND1, IMAXG, JMAXG,
3 C IMAX, IOUT, JMAX, JOUT, IMX, JOUTHF, JMXGHF, KMX2, LMX2, MNWAV,
4 C ht, psfc, t, u, v, q, A, B,
5 C ims , jms, kms, ime , jme, kme,
6 C ids , jds, kds, ide , jde, kde,
7 C its , jts, kts, ite , jte, kte )
8
9 INTEGER ims , jms, kms, ime , jme, kme,
10 C ids , jds, kds, ide , jde, kde,
11 C its , jts, kts, ite , jte, kte
12
13 INTEGER kst
14 DIMENSION A(kms:kme+1), B(kms:kme+1)
15 DIMENSION ht(ims:ime,jms:jme), psfc(ims:ime,jms:jme)
16 DIMENSION t(ims:ime,jms:jme,kms:kme),q(ims:ime,jms:jme,kms:kme)
17 DIMENSION u(ims:ime,jms:jme,kms:kme),v(ims:ime,jms:jme,kms:kme)
18 C
19 DIMENSION ID(5)
20 DIMENSION LAG(MEND1,NEND1)
21 DIMENSION IFAX(10),TRIGS(IMAX)
22 DIMENSION PNM (MNWAV,JOUTHF),DPNM(MNWAV,JOUTHF),
23 1 ESINCL(JOUT),ECOSCL(JOUT),GCOSCL(JMAXG)
24 C
25 DIMENSION QDAT1(KMX2,MNWAV),QDAT2(KMX2,MNWAV),
26 1 QPHIS( 2,MNWAV)
27 C
28 CrizviDIMENSION PS (IOUT*JOUT),PSX (IOUT*JOUT),PSY (IOUT*JOUT)
29 DIMENSION IDA (IOUT*JOUT)
30 DIMENSION EDAT1(IMX *JOUT*KLMAX),EDAT2(IMX *JOUT*KLMAX)
31 DIMENSION TMP(IOUT*JOUT)
32 DIMENSION TMP1(IOUT*JOUT*KLMAX),TMP2(IOUT*JOUT*KLMAX)
33 C
34 Crizvi REAL*8 GAUL(JMAXG),GAUW(JMAXG)
35 REAL*8 PPNM(MNWAV),HHNM(MNWAV)
36 C
37 C
38 DIMENSION KTOUT(82)
39 C
40 DATA IT,JT,BS,FM/4,4,0.0,1.0E10/
41 DATA ER,GASR,GRAV,GAMMA/6371.E03,287.04,9.8,5.0E-3/
42 DATA ITSAV,IPSL/1,1/
43 DATA RHMIN,TDMAX/1.0E-6,31.0/
44 DATA KTOUT/82*999/
45 DATA IWVHST,ITPGFL
46 1 / 1, 3/
47 DATA NGH,NWU,NWV,NTM,NTD,NRR,NRV,NPV
48 1 / 1, 2, 3, 4, 5, 6, 7, 8/
49 DATA KTPART /-1/
50 DATA KTOHK,IHKFL/42,99/
51
52 NAMELIST /NAMOUT/ KTPART,KTOUT,ITSAV,KTOHK,ID_PROC
53 NAMELIST /NAMFIL/ IWVHST,ITPGFL
54 NAMELIST /NAMDBG/ NOMTN,IDEBUG,IT,JT,IPSL
55 NAMELIST /NAMSLP/ MINKTI
56 C------------------------------------------------------------------
57 C Check parameters
58 C------------------------------------------------------------------
59
60 IF( IOUT /= (ide-ids+1) .or. JOUT /= (jde-jds+1) .or.
61 C KLMAX /= (kte-kts+1) ) then
62 print*,' In W2GCONV IOUT,JOUT, KLMAX = ',
63 C IOUT,JOUT, KLMAX
64 print*,' Dims mismatch '
65 stop
66 end if
67
68 IF (IGRD.NE.IOUT.OR.JGRD.NE.JOUT) THEN
69 WRITE(*,*)'OUT GRID DOES NOT MATCH'
70 WRITE(*,*)'IOUT=',IOUT,'JOUT=',JOUT
71 WRITE(*,*)'IGRD=',IGRD,'JGRD=',JGRD
72 STOP 9998
73 END IF
74 C#######################################################################
75 C +++ COMPUTE AND READ IN CONSTANTS +++
76 C#######################################################################
77 READ (95,NAMOUT)
78 READ (95,NAMFIL)
79 READ (95,NAMDBG)
80 READ (95,NAMSLP)
81 WRITE(96,NAMOUT)
82 WRITE(96,NAMFIL)
83 WRITE(96,NAMDBG)
84 WRITE(96,NAMSLP)
85 C------------------------------------------------------------------
86 C Init variables
87 C------------------------------------------------------------------
88 NNMAX=12/MINKTI+1
89 C------------------------------------------------------------------
90 C CHECK DATE SURFACE DATA FILE
91 C------------------------------------------------------------------
92 IO=IOUT
93 JO=JOUT
94 IJOUT=IOUT*JOUT
95 IJMAX=IMAX*JOUT
96 IJKOUT=IOUT*JOUT*KLMAX
97 GINV =1./GRAV
98 C-----------------------------------------------------------------------
99 CALL SETARY(LAG,MEND1,NEND1,JEND1)
100 C-----------------------------------------------------------------------
101 C +++ SET UP FFT COEFFICIENTS +++
102 C-----------------------------------------------------------------------
103 TRIGS(:) = 0.
104 CALL RFFTIM(IMAX,TRIGS,IFAX)
105 C-----------------------------------------------------------------------
106 C +++ SETUP GAUSSIAN LATITUDES +++
107 C-----------------------------------------------------------------------
108 CLSW CALL GAUSS(GAUL,GAUW,JMAXG)
109 CLSW COSCL(1:JMAXG)=GAUL(1:JMAXG)
110 C-----------------------------------------------------------------------
111 C +++ SETUP ZNM ETC +++
112 C-----------------------------------------------------------------------
113 CALL ZNME2PXX
114 I(MEND1,MNWAV,JOUTHF,
115 O PNM ,DPNM ,
116 W PPNM ,HHNM)
117 PAI=4.0*ATAN(1.0)
118 DPAI=PAI/FLOAT(JOUT-1)
119 DO 130 J=1,JOUT
120 ESINCL(J)=SIN(DPAI*FLOAT(J-1))
121 ECOSCL(J)=COS(DPAI*FLOAT(J-1))
122 130 CONTINUE
123 C-----------------------------------------------------------------------
124 C +++ READ IN TOPOGRAPHY FILE +++
125 C-----------------------------------------------------------------------
126 IF(ITPGFL.NE.0) THEN
127 READ(ITPGFL) MDIM,DPHIX,IDIM,JDIM,JHMSPH
128 DPHIX=DPHIX
129 IDIM=IDIM
130 JDIM=JDIM
131 IF(MDIM.NE.MEND1.OR.JHMSPH.NE.0) THEN
132 WRITE(96,901) MDIM,MEND1,JHMSPH
133 901 FORMAT(1H ,'*** ERROR IN TOPOGRPHY FILE ***'/
134 1 ,1X,'(MDIM,MEND1)=(',I3,',',I3,')'/
135 2 ,1X,'JHMSPH = ',I3/)
136 STOP 3002
137 END IF
138 READ(ITPGFL) ((QDAT1(KKK,NNN),KKK=1,2),NNN=1,MNWAV)
139 CLOSE (ITPGFL)
140 C
141 CALL REOWV
142 1 (QDAT1,QPHIS,MEND1,NEND1,JEND1,MNWAV,KMX2,2,0,0,2,LAG)
143 !shc-rizvi start
144 print*,'1. Doing Wave to grid transform for Topo'
145 CALL W2G
146 I (MEND1,NEND1,JEND1,MNWAV,IMAX,JOUT,IMX,JOUTHF, 1,
147 I IFAX ,TRIGS,PNM ,QPHIS,
148 O EDAT1,
149 W EDAT2)
150 !shc-rizvi end
151 C
152 ELSE
153 READ(IWVHST,END=1033) ID,KT,ISTP,FSEC,NREC,
154 1 KM,(A(K),K=1,KM),(B(K),K=1,KM)
155 if( IOUT /= (ide-ids+1) .or. JOUT /= (jde-jds+1) .or.
156 C KLMAX /= (kte-kts+1) ) then
157 print*,' In W2GCONV IOUT,JOUT, KLMAX = ',
158 C IOUT,JOUT, KLMAX
159 print*,' Dims mismatch '
160 stop
161 endif
162 IF(KT.EQ.-1.OR.KT.EQ.0) THEN
163 READ(IWVHST) (QDAT2(KKK,1),KKK=1,2*MNWAV),
164 1 (QDAT1(KKK,1),KKK=1,2*MNWAV)
165 Crizvi CALL REOWV
166 Crizvi 1 (QDAT1,QPHIS,MEND1,NEND1,JEND1,MNWAV,KMX2,2,0,0,2,LAG)
167 !shc-rizvi start
168 print*,'2. Doing Wave to grid transform for Topo'
169 CALL W2G
170 I (MEND1,NEND1,JEND1,MNWAV,IMAX,JOUT,IMX,JOUTHF, 1,
171 I IFAX ,TRIGS,PNM ,QDAT1,
172 O EDAT1,
173 W EDAT2)
174 !shc-rizvi end
175 READ(IWVHST)
176 READ(IWVHST)
177 READ(IWVHST)
178 READ(IWVHST)
179 GO TO 1033
180 ELSE
181 READ(IWVHST)
182 READ(IWVHST)
183 READ(IWVHST)
184 READ(IWVHST)
185 READ(IWVHST)
186 END IF
187 END IF
188 1033 CONTINUE
189
190 REWIND(IWVHST)
191 C
192 CALL CUT(TMP, IOUT,JOUT, 1,EDAT1,IMAX,INTVL)
193 DO 160 IJ=1,IJOUT
194 TMP(IJ) =GINV*TMP(IJ)
195 160 CONTINUE
196 IF(IDEBUG.GE.1) THEN
197 CALL GOUT(TMP,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'PHIS',1)
198 END IF
199 C write(989,'(10f10.3)') (ETPG(IJ),IJ=1,IJOUT)
200 call reorder_for_wrf(TMP,ht,IOUT,JOUT, 1,
201 C ims , jms, 1, ime , jme, 1,
202 C ids , jds, 1, ide , jde, 1,
203 C its , jts, 1, ite , jte, 1)
204 C
205 C#######################################################################
206 C +++ READ IN HISTORY FILE +++
207 C +++ WAVE TO GRID CONVERSION START +++
208 C#######################################################################
209 C...IWVHST; WAVE HISTORY FILE
210 C...ISFCFL; GAUSS GRID SURFACE FILE
211 C...IGHIST; EQUAL LAT-LON GRID P-HISTORY FILE
212 C IF IGHIST=0, NO FILE OUTPUT
213 C...IPGUES; EQUAL LAT-LON GRID P-GUESS FILE
214 C IF IPGUES=0, NO FILE OUTPUT
215 C...IXXXFL; GRID ETA-TEMPORARY FILE FOR XXX PASSED TO ETA-P CONVERSION
216 C IF IXXXFL=0, NO XXX TEMPORARY FILE OUTPUT
217 C XXX IS TMP,SPH,U,V,DIV,OMG
218 C
219 ICOUNT=1
220 C
221 2000 CONTINUE
222 C
223 IF(KTOUT(ICOUNT).EQ.999.AND.KTOUT(1).NE.-999) GOTO 1001
224 READ(IWVHST) ID,KT,ISTP,FSEC,NREC,
225 1 KM,(A(K),K=1,KM),(B(K),K=1,KM)
226 WRITE(96,666) ID,KT
227 666 FORMAT(1H ,'HISTORY FILE READ IN AT',6I5)
228 C IF(KT.LE.0) THEN
229 C WRITE(96,*)'KT=',KT
230 C END IF
231 C
232 IF(KM.NE.KMAX) THEN
233 WRITE(96,902)
234 902 FORMAT(1X ,'*** ERROR IN NO. OF ETALEV ***')
235 STOP 3001
236 END IF
237 !
238 IF(KT.EQ.KTOUT(ICOUNT).OR.KTOUT(1).EQ.-999) THEN
239 ICOUNT=ICOUNT+1
240 WRITE(96,905) ID,KT
241 905 FORMAT(1H ,2X,'YEAR=',I4,4X,'MONTH=',I2,4X,'DAY=',I2,4X,'HOUR=',
242 * I2,4X,'WEEK=',I1,4X,'FCST=',I3)
243 C---------------------------------------------------------------------
244 C...SURFACE PRESSURE AND ITS GRADIENTS
245 C---------------------------------------------------------------------
246 READ(IWVHST) (QDAT1(KMN,1),KMN=1,2*MNWAV)
247 WRITE(96,*) ID,KT,'qdat1'
248 C...IPXY=1;OUTPUT IS PS,D(PS)/DX,-D(PS)/DY
249 IPXY=1
250 print*,' Doing Wave to grid transform for Psfc'
251 CALL W2GPXY
252 I (MEND1,NEND1,JEND1,MNWAV,IMAX,JOUT,IMX,JOUTHF,
253 I IFAX ,TRIGS,PNM ,DPNM ,IPXY,IPSL,QDAT1,
254 O EDAT1,
255 W EDAT2)
256 CALL CUT(TMP ,IOUT,JOUT,1,EDAT1( 1),IMAX,INTVL)
257 Crizvi CALL CUT(PSX ,IOUT,JOUT,1,EDAT1( IJMAX+1),IMAX,INTVL)
258 Crizvi CALL CUT(PSY ,IOUT,JOUT,1,EDAT1(2*IJMAX+1),IMAX,INTVL)
259 IF(IDEBUG.GE.1) THEN
260 CALL GOUT(TMP,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'PSFC',1)
261 ! CALL GOUT(PSX,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'DPSX',1)
262 ! CALL GOUT(PSY,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'DPSY',1)
263 END IF
264 call reorder_for_wrf(TMP, psfc,IOUT,JOUT, 1,
265 C ims , jms, 1, ime , jme, 1,
266 C ids , jds, 1, ide , jde, 1,
267 C its , jts, 1, ite , jte, 1)
268 C---------------------------------------------------------------------
269 C...TEMPERATURE
270 C---------------------------------------------------------------------
271 READ (IWVHST) QDAT1
272 print*,' Doing Wave to grid transform for Temperature'
273 CALL W2G
274 I (MEND1,NEND1,JEND1,MNWAV,IMAX,JOUT,IMX,JOUTHF,KMAX,
275 I IFAX ,TRIGS,PNM ,QDAT1,
276 O EDAT1,
277 W EDAT2)
278 CALL CUT(TMP1,IOUT,JOUT,KMAX,EDAT1,IMAX,INTVL)
279 IF(IDEBUG.GE.1) THEN
280 K=KMAX
281 CALL GOUT(TMP1,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'ETMP',K)
282 END IF
283 call reorder_for_wrf(TMP1, t, IOUT,JOUT, KLMAX,
284 C ims , jms, kms, ime , jme, kme,
285 C ids , jds, kds, ide , jde, kde,
286 C its , jts, kts, ite , jte, kte )
287 C----------------------------------------------------------------------
288 C... VORTICITY & DIVERGENCE
289 C----------------------------------------------------------------------
290 C...VORTICITY
291 READ(IWVHST) QDAT1
292 C...DIVERGENCE
293 READ(IWVHST) QDAT2
294 print*,' Doing Wave to grid transform for Wind'
295 CALL W2GUV
296 I (MEND1,NEND1,JEND1 ,MNWAV,IMAX,JOUT,IMX ,JOUTHF,KMAX,
297 I IFAX ,TRIGS,ESINCL,ER ,PNM ,DPNM,QDAT1,QDAT2 ,
298 O EDAT1,EDAT2,
299 W TMP2)
300 CALL CUT(TMP1,IOUT,JOUT,KMAX,EDAT1,IMAX,INTVL)
301 CALL CUT(TMP2,IOUT,JOUT,KMAX,EDAT2,IMAX,INTVL)
302 IF(IDEBUG.GE.1) THEN
303 K=KMAX
304 CALL GOUT(TMP1,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'EU ',K)
305 CALL GOUT(TMP2,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'EV ',K)
306 END IF
307 C
308 call reorder_for_wrf(TMP1, U ,IOUT,JOUT, KLMAX,
309 C ims , jms, kms, ime , jme, kme,
310 C ids , jds, kds, ide , jde, kde,
311 C its , jts, kts, ite , jte, kte )
312 call reorder_for_wrf(TMP2, V ,IOUT,JOUT, KLMAX,
313 C ims , jms, kms, ime , jme, kme,
314 C ids , jds, kds, ide , jde, kde,
315 C its , jts, kts, ite , jte, kte )
316 C----------------------------------------------------------------------
317 C...SPECIFIC HUMIDITY
318 C----------------------------------------------------------------------
319 READ(IWVHST) QDAT1
320 print*,' Doing Wave to grid transform for Moisture'
321 CALL W2G
322 I (MEND1,NEND1,JEND1,MNWAV,IMAX,JOUT,IMX,JOUTHF,KMAX,
323 I IFAX ,TRIGS,PNM ,QDAT1,
324 O EDAT1,
325 W EDAT2)
326 print*,' Done Wave to grid transform for Moisture calling cut'
327 CALL CUT(TMP1,IOUT,JOUT,KMAX,EDAT1,IMAX,INTVL)
328 C rizvi fix negative humidity to 0.000001
329 CALL FIX_NEG_MOIST(TMP1,IOUT,JOUT,KMAX)
330 print*,' Done cut calling reorder_for_wrf'
331 c---------------------------------------------------------------
332 call reorder_for_wrf(TMP1, Q,IOUT,JOUT, KLMAX,
333 C ims , jms, kms, ime , jme, kme,
334 C ids , jds, kds, ide , jde, kde,
335 C its , jts, kts, ite , jte, kte )
336 print*,' Done calling reorder_for_wrf'
337 c---------------------------------------------------------------
338 IF(IDEBUG.GE.1) THEN
339 K=KMAX
340 CALL GOUT(TMP1,IDA,IO,JO,1,1,IO,JO,IT,JT,BS,FM,'ESPH',K)
341 END IF
342 ELSE
343 READ(IWVHST)
344 READ(IWVHST)
345 READ(IWVHST)
346 READ(IWVHST)
347 READ(IWVHST)
348 GO TO 2000
349 END IF
350 C
351 C
352 GO TO 2000
353 C
354 1001 CONTINUE
355 RETURN
356 END SUBROUTINE W2GCONV
357
358 SUBROUTINE FIX_NEG_MOIST(X,N1,N2,N3)
359 C---------------------------------------------------------
360 C This fixes less than 0.00000 input array to 0.000001
361 C Syed RH Rizvi 08/24/2004
362 C---------------------------------------------------------
363 DIMENSION X(N1,N2,N3)
364 DO K = 1,N3
365 DO J= 1,N2
366 DO I = 1, N1
367 C IF(X(I,J,K) < 0) X(I,J,K) = 0.
368 IF(X(I,J,K) < 0.000001) X(I,J,K) = 0.000001
369 END DO
370 END DO
371 END DO
372 RETURN
373 END SUBROUTINE FIX_NEG_MOIST
374
375 SUBROUTINE reorder_for_wrf_2d(kma, X, n1,n2 ,
376 C ims , jms, ime , jme,
377 C ids , jds, ide , jde,
378 C its , jts, ite , jte )
379
380 IMPLICIT NONE
381 INTEGER,intent(in) ::
382 C ims , jms, ime , jme,
383 C ids , jds, ide , jde,
384 C its , jts, ite , jte
385 integer, intent(in) :: n1,n2
386 real, intent(in) :: kma(n1*n2)
387 real, intent( out) :: X(ims:ime,jms:jme)
388 real, dimension(1:n1,1:n2) :: wrf
389 real :: tmp(n1,n2)
390 integer :: i,j,ist
391
392 C
393 ist = 1
394 do j=1,n2
395 do i=1,n1
396 tmp(i,j) = kma(ist)
397 ist = ist + 1
398 end do
399 end do
400 C
401 do j= 1,n2
402 wrf(1:n1,j) = tmp(1:n1,n2-j+1)
403 end do
404 C
405 do j = jts, jte
406 X(its:ite,j) = wrf(its:ite,j)
407 end do
408
409 RETURN
410 END SUBROUTINE reorder_for_wrf_2d
411 SUBROUTINE reorder_for_wrf_3d(kma, X, n1,n2,n3,
412 C ims , jms, kms, ime , jme, kme,
413 C ids , jds, kds, ide , jde, kde,
414 C its , jts, kts, ite , jte, kte )
415
416 IMPLICIT NONE
417 INTEGER,intent(in) ::
418 C ims , jms, kms, ime , jme, kme,
419 C ids , jds, kds, ide , jde, kde,
420 C its , jts, kts, ite , jte, kte
421 integer, intent(in) :: n1,n2,n3
422 real, intent(in) :: kma(n1*n2*n3)
423 real, intent(OUT) :: X(ims:ime,jms:jme,kms:kme)
424 real, dimension(1:n1,1:n2,1:n3) :: wrf
425 real :: tmp(n1,n2,n3)
426 integer :: i,j,k,ist
427
428 C
429 ist = 1
430 do k=1,n3
431 do j=1,n2
432 do i=1,n1
433 tmp(i,j,k) = kma(ist)
434 ist = ist + 1
435 end do
436 end do
437 end do
438 C
439 do k =1,n3
440 do j= 1,n2
441 wrf(1:n1,j,k) = tmp(1:n1,n2-j+1,k)
442 end do
443 end do
444 C
445 do k=kts,kte
446 do j = jts, jte
447 X(its:ite,j,k) = wrf(its:ite,j,k)
448 end do
449 end do
450 C
451 RETURN
452 END SUBROUTINE reorder_for_wrf_3d
453
454 SUBROUTINE reorder_for_wrf(kma, wrf, n1,n2,n3,
455 C ims , jms, kms, ime , jme, kme,
456 C ids , jds, kds, ide , jde, kde,
457 C its , jts, kts, ite , jte, kte )
458 IMPLICIT NONE
459 INTEGER,intent(in) ::
460 C ims , jms, kms, ime , jme, kme,
461 C ids , jds, kds, ide , jde, kde,
462 C its , jts, kts, ite , jte, kte
463 integer, intent(in) :: n1,n2,n3
464 real, intent(in) :: kma(n1*n2*n3)
465 real, intent(OUT) :: wrf(ims:ime,jms:jme,kms:kme)
466 real :: XWRF(ims:ime,jms:jme,kms:kme)
467 real :: tmp(n1,n2,n3), xtmp(n1,n2,n3)
468 integer :: i,j,k, half
469 C
470 tmp(1:n1,1:n2,1:n3) =
471 & RESHAPE( kma(1:n1*n2*n3),(/n1,n2,n3/) )
472 C
473 do k =1,n3
474 do j= 1,n2
475 xtmp(1:n1,j,k) = tmp(1:n1,n2-j+1,k)
476 end do
477 end do
478 C
479 do k=kts,kte
480 do j = jts, jte
481 XWRF(its:ite,j,k) = xtmp(its:ite,j,k)
482 end do
483 end do
484 C
485 half = (ite-its+1)/2
486 do k=kts,kte
487 do j= jts,jte
488 do i=its,ite
489 if( i <= half)then
490 WRF(half+i,j,k) = XWRF(i,j,k)
491 else
492 WRF(i-half,j,k) = XWRF(i,j,k)
493 end if
494 end do
495 end do
496 end do
497 C
498 RETURN
499 END SUBROUTINE reorder_for_wrf