Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / var / convertor / kma_wave2grid / W2GCONV.inc
blobea612f8a83e614b6e4fce9f74ed2736cba5b3bd9
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)
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)
25       DIMENSION QDAT1(KMX2,MNWAV),QDAT2(KMX2,MNWAV),
26      1          QPHIS(   2,MNWAV)
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)
34 Crizvi      REAL*8 GAUL(JMAXG),GAUW(JMAXG)
35       REAL*8 PPNM(MNWAV),HHNM(MNWAV)
38       DIMENSION KTOUT(82)
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/
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------------------------------------------------------------------
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 
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)
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
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
189         REWIND(IWVHST)
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)
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
218       ICOUNT=1
220  2000 CONTINUE
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                                                           
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
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
351       GO TO 2000
353  1001 CONTINUE
354       RETURN
355       END SUBROUTINE W2GCONV
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     
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   )
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
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 
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  )
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
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 
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
450       RETURN
451       END SUBROUTINE reorder_for_wrf_3d
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
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 
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
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
497        RETURN
498       END SUBROUTINE reorder_for_wrf