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