merge standard release WRF/WPS V3.0.1.1 into wrffire
[wrffire.git] / wrfv2_fire / external / io_grib_share / io_grib_share.F
blob70b9ed0ca4641551ff31e47305ae820b2eae7bc4
2 ! Todd Hutchinson
3 ! WSI
4 ! August 17, 2005
6 ! Routines in this file are shared by io_grib1 and io_grib2
9 !*****************************************************************************
11 SUBROUTINE get_dims(MemoryOrder, Start, End, ndim, x_start, x_end, y_start, &
12      y_end, z_start, z_end)
13   IMPLICIT NONE
14   CHARACTER (LEN=*)    ,INTENT(IN)    :: MemoryOrder
15   INTEGER              ,INTENT(OUT)   :: ndim,x_start,x_end,y_start
16   INTEGER              ,INTENT(OUT)   :: y_end,z_start,z_end
17   integer ,dimension(*),intent(in)    :: Start, End
18   CHARACTER (LEN=1)                   :: char
19   INTEGER                             :: idx
20   CHARACTER (LEN=3)                   :: MemoryOrderLcl
22   x_start = 1
23   x_end   = 1
24   y_start = 1
25   y_end   = 1
26   z_start = 1
27   z_end   = 1
29   !
30   ! Note: Need to add "char == 'S'" for boundary conditions
31   !
33   ndim = 0
35   ! Fix for out-of-bounds references
36   MemoryOrderLcl = '   '
37   do idx=1,len_trim(MemoryOrder)
38      MemoryOrderLcl(idx:idx) = MemoryOrder(idx:idx)
39   enddo
40   !
41   ! First, do the special boundary cases.  These do not seem to 
42   !    
43   if ((MemoryOrderLcl(1:3) .eq. 'XSZ') &
44        .or. (MemoryOrderLcl(1:3) .eq. 'XEZ')) then
45      x_start = Start(3)
46      x_end = End(3)
47      y_start = Start(1)
48      y_end = End(1)
49      z_start = Start(2)
50      z_end = End(2)
51      ndim = 3
52   else if ((MemoryOrderLcl(1:3) .eq. 'YSZ') .or. &
53        (MemoryOrderLcl(1:3) .eq. 'YEZ')) then
54      x_start = Start(1)
55      x_end = End(1)
56      y_start = Start(3)
57      y_end = End(3)
58      z_start = Start(2)
59      z_end = End(2)
60      ndim = 3
61   else if ((MemoryOrderLcl(1:2) .eq. 'YS') .or. &
62        (MemoryOrderLcl(1:2) .eq. 'YE')) then
63      x_start = Start(1)
64      x_end = End(1)
65      y_start = Start(2)
66      y_end = End(2)
67      ndim = 2
68   else if ((MemoryOrderLcl(1:2) .eq. 'XS') .or. &
69        (MemoryOrderLcl(1:2) .eq. 'XE')) then
70      x_start = Start(2)
71      x_end = End(2)
72      y_start = Start(1)
73      y_end = End(1)
74      ndim = 2
75   else if ((MemoryOrderLcl(1:1) .eq. 'C') .or. (MemoryOrderLcl(1:1) .eq. 'c')) then 
76      ! This is for "non-decomposed" fields
77      x_start = Start(1)
78      x_end = End(1)
79 !     y_start = Start(2)
80 !     y_end = End(2)
81 !     z_start = Start(3)
82 !     z_end = End(3)
83      ndim = 3
84   else
85      do idx=1,len_trim(MemoryOrderLcl)
86         char = MemoryOrderLcl(idx:idx)
87         if ((char == 'X') .or. (char == 'x')) then
88            x_start = Start(idx)
89            x_end   = End(idx)
90            ndim = ndim + 1
91         else if ((char == 'Y') .or. (char == 'y')) then
92            y_start = Start(idx)
93            y_end   = End(idx)
94            ndim = ndim + 1
95         else if ((char == 'Z') .or. (char == 'z')) then
96            z_start = Start(idx)
97            z_end   = End(idx)
98            ndim = ndim + 1
99         else if (char == '0') then
100            ! Do nothing, this indicates field is a scalar.
101            ndim = 0
102         else
103            call wrf_message('Invalid Dimension in get_dims: '//char)
104         endif
105      enddo
106   endif
108 END SUBROUTINE get_dims
110 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
112 SUBROUTINE geth_idts (ndate, odate, idts)
114   IMPLICIT NONE
116   !  From 2 input mdates ('YYYY-MM-DD HH:MM:SS.ffff'), 
117   !  compute the time difference.
119   !  on entry     -  ndate  -  the new hdate.
120   !                  odate  -  the old hdate.
122   !  on exit      -  idts    -  the change in time in seconds.
124   CHARACTER (LEN=*) , INTENT(INOUT) :: ndate, odate
125   REAL              , INTENT(OUT)   :: idts
127   !  Local Variables
129   !  yrnew    -  indicates the year associated with "ndate"
130   !  yrold    -  indicates the year associated with "odate"
131   !  monew    -  indicates the month associated with "ndate"
132   !  moold    -  indicates the month associated with "odate"
133   !  dynew    -  indicates the day associated with "ndate"
134   !  dyold    -  indicates the day associated with "odate"
135   !  hrnew    -  indicates the hour associated with "ndate"
136   !  hrold    -  indicates the hour associated with "odate"
137   !  minew    -  indicates the minute associated with "ndate"
138   !  miold    -  indicates the minute associated with "odate"
139   !  scnew    -  indicates the second associated with "ndate"
140   !  scold    -  indicates the second associated with "odate"
141   !  i        -  loop counter
142   !  mday     -  a list assigning the number of days in each month
144   CHARACTER (LEN=24) :: tdate
145   INTEGER :: olen, nlen
146   INTEGER :: yrnew, monew, dynew, hrnew, minew, scnew
147   INTEGER :: yrold, moold, dyold, hrold, miold, scold
148   INTEGER :: mday(12), i, newdys, olddys
149   LOGICAL :: npass, opass
150   INTEGER :: isign
151   CHARACTER (LEN=300) :: wrf_err_message
152   INTEGER :: ndfeb
154   IF (odate.GT.ndate) THEN
155      isign = -1
156      tdate=ndate
157      ndate=odate
158      odate=tdate
159   ELSE
160      isign = 1
161   END IF
163   !  Assign the number of days in a months
165   mday( 1) = 31
166   mday( 2) = 28
167   mday( 3) = 31
168   mday( 4) = 30
169   mday( 5) = 31
170   mday( 6) = 30
171   mday( 7) = 31
172   mday( 8) = 31
173   mday( 9) = 30
174   mday(10) = 31
175   mday(11) = 30
176   mday(12) = 31
178   !  Break down old hdate into parts
180   hrold = 0
181   miold = 0
182   scold = 0
183   olen = LEN(odate)
185   READ(odate(1:4),  '(I4)') yrold
186   READ(odate(6:7),  '(I2)') moold
187   READ(odate(9:10), '(I2)') dyold
188   IF (olen.GE.13) THEN
189      READ(odate(12:13),'(I2)') hrold
190      IF (olen.GE.16) THEN
191         READ(odate(15:16),'(I2)') miold
192         IF (olen.GE.19) THEN
193            READ(odate(18:19),'(I2)') scold
194         END IF
195      END IF
196   END IF
198   !  Break down new hdate into parts
200   hrnew = 0
201   minew = 0
202   scnew = 0
203   nlen = LEN(ndate)
205   READ(ndate(1:4),  '(I4)') yrnew
206   READ(ndate(6:7),  '(I2)') monew
207   READ(ndate(9:10), '(I2)') dynew
208   IF (nlen.GE.13) THEN
209      READ(ndate(12:13),'(I2)') hrnew
210      IF (nlen.GE.16) THEN
211         READ(ndate(15:16),'(I2)') minew
212         IF (nlen.GE.19) THEN
213            READ(ndate(18:19),'(I2)') scnew
214         END IF
215      END IF
216   END IF
218   !  Check that the dates make sense.
220   npass = .true.
221   opass = .true.
223   !  Check that the month of NDATE makes sense.
225   IF ((monew.GT.12).or.(monew.LT.1)) THEN
226      PRINT*, 'GETH_IDTS:  Month of NDATE = ', monew
227      npass = .false.
228   END IF
230   !  Check that the month of ODATE makes sense.
232   IF ((moold.GT.12).or.(moold.LT.1)) THEN
233      PRINT*, 'GETH_IDTS:  Month of ODATE = ', moold
234      opass = .false.
235   END IF
237   !  Check that the day of NDATE makes sense.
239   IF (monew.ne.2) THEN
240      ! ...... For all months but February
241      IF ((dynew.GT.mday(monew)).or.(dynew.LT.1)) THEN
242         PRINT*, 'GETH_IDTS:  Day of NDATE = ', dynew
243         npass = .false.
244      END IF
245   ELSE IF (monew.eq.2) THEN
246      ! ...... For February
247      IF ((dynew.GT.ndfeb(yrnew)).OR.(dynew.LT.1)) THEN
248         PRINT*, 'GETH_IDTS:  Day of NDATE = ', dynew
249         npass = .false.
250      END IF
251   END IF
253   !  Check that the day of ODATE makes sense.
255   IF (moold.ne.2) THEN
256      ! ...... For all months but February
257      IF ((dyold.GT.mday(moold)).or.(dyold.LT.1)) THEN
258         PRINT*, 'GETH_IDTS:  Day of ODATE = ', dyold
259         opass = .false.
260      END IF
261   ELSE IF (moold.eq.2) THEN
262      ! ....... For February
263      IF ((dyold.GT.ndfeb(yrold)).or.(dyold.LT.1)) THEN
264         PRINT*, 'GETH_IDTS:  Day of ODATE = ', dyold
265         opass = .false.
266      END IF
267   END IF
269   !  Check that the hour of NDATE makes sense.
271   IF ((hrnew.GT.23).or.(hrnew.LT.0)) THEN
272      PRINT*, 'GETH_IDTS:  Hour of NDATE = ', hrnew
273      npass = .false.
274   END IF
276   !  Check that the hour of ODATE makes sense.
278   IF ((hrold.GT.23).or.(hrold.LT.0)) THEN
279      PRINT*, 'GETH_IDTS:  Hour of ODATE = ', hrold
280      opass = .false.
281   END IF
283   !  Check that the minute of NDATE makes sense.
285   IF ((minew.GT.59).or.(minew.LT.0)) THEN
286      PRINT*, 'GETH_IDTS:  Minute of NDATE = ', minew
287      npass = .false.
288   END IF
290   !  Check that the minute of ODATE makes sense.
292   IF ((miold.GT.59).or.(miold.LT.0)) THEN
293      PRINT*, 'GETH_IDTS:  Minute of ODATE = ', miold
294      opass = .false.
295   END IF
297   !  Check that the second of NDATE makes sense.
299   IF ((scnew.GT.59).or.(scnew.LT.0)) THEN
300      PRINT*, 'GETH_IDTS:  SECOND of NDATE = ', scnew
301      npass = .false.
302   END IF
304   !  Check that the second of ODATE makes sense.
306   IF ((scold.GT.59).or.(scold.LT.0)) THEN
307      PRINT*, 'GETH_IDTS:  Second of ODATE = ', scold
308      opass = .false.
309   END IF
311   IF (.not. npass) THEN
312      WRITE( wrf_err_message , * ) &
313           'module_date_time: geth_idts: Bad NDATE: ', ndate(1:nlen)
314      CALL wrf_error_fatal ( TRIM ( wrf_err_message ) )
315   END IF
317   IF (.not. opass) THEN
318      WRITE( wrf_err_message , * ) &
319           'module_date_time: geth_idts: Bad ODATE: ', odate(1:olen)
320      CALL wrf_error_fatal ( TRIM ( wrf_err_message ) )
321   END IF
323   !  Date Checks are completed.  Continue.
325   !  Compute number of days from 1 January ODATE, 00:00:00 until ndate
326   !  Compute number of hours from 1 January ODATE, 00:00:00 until ndate
327   !  Compute number of minutes from 1 January ODATE, 00:00:00 until ndate
329   newdys = 0
330   DO i = yrold, yrnew - 1
331      newdys = newdys + (365 + (ndfeb(i)-28))
332   END DO
334   IF (monew .GT. 1) THEN
335      mday(2) = ndfeb(yrnew)
336      DO i = 1, monew - 1
337         newdys = newdys + mday(i)
338      END DO
339      mday(2) = 28
340   END IF
342   newdys = newdys + dynew-1
344   !  Compute number of hours from 1 January ODATE, 00:00:00 until odate
345   !  Compute number of minutes from 1 January ODATE, 00:00:00 until odate
347   olddys = 0
349   IF (moold .GT. 1) THEN
350      mday(2) = ndfeb(yrold)
351      DO i = 1, moold - 1
352         olddys = olddys + mday(i)
353      END DO
354      mday(2) = 28
355   END IF
357   olddys = olddys + dyold-1
359   !  Determine the time difference in seconds
361   idts = (newdys - olddys) * 86400
362   idts = idts + (hrnew - hrold) * 3600
363   idts = idts + (minew - miold) * 60
364   idts = idts + (scnew - scold)
366   IF (isign .eq. -1) THEN
367      tdate=ndate
368      ndate=odate
369      odate=tdate
370      idts = idts * isign
371   END IF
373 END SUBROUTINE geth_idts
375 !*****************************************************************************
377 SUBROUTINE get_vert_stag(VarName,Stagger,vert_stag)
378   
379   character (LEN=*) :: VarName
380   character (LEN=*) :: Stagger
381   logical           :: vert_stag
383   if ((index(Stagger,'Z') > 0) .or. (VarName .eq. 'DNW') &
384        .or.(VarName .eq. 'RDNW')) then
385      vert_stag = .true.
386   else
387      vert_stag = .false.
388   endif
389 end SUBROUTINE
391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393 FUNCTION ndfeb ( year ) RESULT (num_days)
394   
395   ! Compute the number of days in February for the given year
396   
397   IMPLICIT NONE
398   
399   INTEGER :: year
400   INTEGER :: num_days
401   
402   num_days = 28 ! By default, February has 28 days ...
403   IF (MOD(year,4).eq.0) THEN  
404      num_days = 29  ! But every four years, it has 29 days ...
405      IF (MOD(year,100).eq.0) THEN
406         num_days = 28  ! Except every 100 years, when it has 28 days ...
407         IF (MOD(year,400).eq.0) THEN
408            num_days = 29  ! Except every 400 years, when it has 29 days.
409         END IF
410      END IF
411   END IF
412   
413 END FUNCTION ndfeb
415 !*****************************************************************************
417 SUBROUTINE get_dimvals(MemoryOrder,x,y,z,dims)
419   IMPLICIT NONE
420   CHARACTER (LEN=*)    ,INTENT(IN)    :: MemoryOrder
421   INTEGER              ,INTENT(IN)    :: x,y,z
422   INTEGER, DIMENSION(*),INTENT(OUT)   :: dims
423   INTEGER                             :: idx
424   CHARACTER (LEN=1) :: char
425   CHARACTER (LEN=3) :: MemoryOrderLcl
427   dims(1) = 1
428   dims(2) = 1
429   dims(3) = 1
431   ! Fix for out-of-bounds references
432   MemoryOrderLcl = '   '
433   do idx=1,len_trim(MemoryOrder)
434      MemoryOrderLcl(idx:idx) = MemoryOrder(idx:idx)
435   enddo
437   !
438   ! Note: Need to add "char == 'S'" for boundary conditions
439   !
441   if ((MemoryOrderLcl(1:3) .eq. 'XSZ') &
442        .or. (MemoryOrderLcl(1:3) .eq. 'XEZ')) then
443      dims(1) = y
444      dims(2) = z
445      dims(3) = x
446   else if ((MemoryOrderLcl(1:3) .eq. 'YSZ') .or. &
447        (MemoryOrderLcl(1:3) .eq. 'YEZ')) then
448      dims(1) = x
449      dims(2) = z
450      dims(3) = y
451   else if ((MemoryOrderLcl(1:2) .eq. 'YS') .or. &
452        (MemoryOrderLcl(1:2) .eq. 'YE')) then
453      dims(1) = x
454      dims(2) = y
455      dims(3) = z
456   else if ((MemoryOrderLcl(1:2) .eq. 'XS') .or. &
457        (MemoryOrderLcl(1:2) .eq. 'XE')) then
458      dims(1) = y
459      dims(2) = x
460      dims(3) = z
461   else if ((MemoryOrderLcl(1:1) .eq. 'C') .or. &
462        (MemoryOrderLcl(1:1) .eq. 'c')) then
463      ! Non-decomposed field
464      dims(1) = x
465      dims(2) = y
466      dims(3) = z
467   else 
468      do idx=1,len_trim(MemoryOrderLcl)
469         char = MemoryOrderLcl(idx:idx)
470         if ((char == 'X') .or. (char == 'x')) then
471            dims(idx) = x
472         else if ((char == 'Y') .or. (char == 'y')) then
473            dims(idx) = y
474         else if ((char == 'Z') .or. (char == 'z')) then
475            dims(idx) = z
476         else if (char == '0') then
477            ! This is a scalar, do nothing.
478         else
479            call wrf_message ('Invalid Dimension in get_dimvals: '//char)
480         endif
481      enddo
482   endif
484 END SUBROUTINE get_dimvals
486 !*****************************************************************************
488 SUBROUTINE get_soil_layers(VarName,soil_layers)
489   
490   character (LEN=*) :: VarName
491   logical           :: soil_layers
493   if ((VarName .eq. 'ZS') .or. (VarName .eq. 'DZS') &
494        .or.(VarName .eq. 'TSLB') .or. (VarName .eq. 'SMOIS') &
495        .or. (VarName .eq. 'SH2O') .or. (VarName .eq. 'KEEPFR3DFLAG') &
496        .or. (VarName .eq. 'SMFR3D')) then
497      soil_layers = .true.
498   else
499      soil_layers = .false.
500   endif
501 end SUBROUTINE
503 !*****************************************************************************
505 SUBROUTINE Transpose_grib(MemoryOrder, di, FieldType, Field, &
506      Start1, End1, Start2, End2, Start3, End3, data, zidx, numrows, numcols)
508   IMPLICIT NONE
510 #include "wrf_io_flags.h"
512   CHARACTER (LEN=*),INTENT(IN)    :: MemoryOrder
513   INTEGER          ,INTENT(IN)    :: Start1,End1,Start2,End2,Start3,End3
514   INTEGER          ,INTENT(IN)    :: di
515   integer          ,intent(inout) :: &
516        Field(di,Start1:End1,Start2:End2,Start3:End3)
517   INTEGER          ,intent(in)    :: FieldType
518   real             ,intent(in)    :: data(*)
519   INTEGER          ,INTENT(IN)    :: zidx, numcols, numrows
520   INTEGER, DIMENSION(3)           :: dims
521   INTEGER                         :: col, row
522   LOGICAL                         :: logicaltype
523   CHARACTER (LEN=1000)            :: msg
524      
525   if ((FieldType == WRF_REAL) .or. (FieldType == WRF_DOUBLE)) then
526      do col=1,numcols
527         do row=1,numrows
528            call get_dimvals(MemoryOrder,col,row,zidx,dims)
529            Field(1:di,dims(1),dims(2),dims(3)) = &
530                 TRANSFER(data((row-1)*numcols+col),Field,1)
531         enddo
532      enddo
533   else if (FieldType == WRF_INTEGER) then
534      do col=1,numcols
535         do row=1,numrows
536            call get_dimvals(MemoryOrder,col,row,zidx,dims)
537            Field(1:di,dims(1),dims(2),dims(3)) = data((row-1)*numcols+col)
538         enddo
539      enddo
540   else
541      write (msg,*)'Reading of type ',FieldType,'from grib data not supported'
542      call wrf_message(msg)
543   endif
546 ! This following seciton is for the logical type.  This caused some problems
547 !   on certain platforms.
549 !  else if (FieldType == WRF_LOGICAL) then
550 !     do col=1,numcols
551 !        do row=1,numrows
552 !           call get_dimvals(MemoryOrder,col,row,zidx,dims)
553 !           Field(1:di,dims(1),dims(2),dims(3)) = &
554 !                TRANSFER(data((row-1)*numcols+col),logicaltype,1)
555 !        enddo
556 !     enddo
557   
558   
559 end SUBROUTINE
561 !*****************************************************************************
563 SUBROUTINE Transpose1D_grib(MemoryOrder, di, FieldType, Field, &
564      Start1, End1, Start2, End2, Start3, End3, data, nelems)
566   IMPLICIT NONE
568 #include "wrf_io_flags.h"
570   CHARACTER (LEN=*),INTENT(IN)    :: MemoryOrder
571   INTEGER          ,INTENT(IN)    :: Start1,End1,Start2,End2,Start3,End3
572   INTEGER          ,INTENT(IN)    :: di
573   integer          ,intent(inout) :: &
574        Field(di,Start1:End1,Start2:End2,Start3:End3)
575   INTEGER          ,intent(in)    :: FieldType
576   real             ,intent(in)    :: data(*)
577   LOGICAL                         :: logicaltype
578   CHARACTER (LEN=1000)            :: msg
579   integer                         :: elemnum,nelems
581   if ((FieldType == WRF_REAL) .or. (FieldType == WRF_DOUBLE)) then
582      do elemnum=1,nelems
583         Field(1:di,elemnum,1,1) = TRANSFER(data(elemnum),Field,1)
584      enddo
585   else if (FieldType == WRF_INTEGER) then
586      do elemnum=1,nelems
587         Field(1:di,elemnum,1,1) = TRANSFER(data(elemnum),Field,1)
588      enddo
589   else
590      write (msg,*)'Reading of type ',FieldType,'from grib1 data not supported'
591      call wrf_message(msg)
592   endif
595 ! This following seciton is for the logical type.  This caused some problems
596 !   on certain platforms.
598 !  else if (FieldType == WRF_LOGICAL) then
599 !     do col=1,numcols
600 !        do row=1,numrows
601 !           call get_dimvals(MemoryOrder,col,row,zidx,dims)
602 !           Field(1:di,dims(1),dims(2),dims(3)) = &
603 !                TRANSFER(data((row-1)*numcols+col),logicaltype,1)
604 !        enddo
605 !     enddo
606   
607   
608 end SUBROUTINE Transpose1D_grib
610 !*****************************************************************************
612 ! Takes a starting date (startTime) in WRF format (yyyy-mm-dd_hh:mm:ss), 
613 !   adds an input number of seconds to the time, and outputs a new date 
614 !   (endTime) in WRF format.
616 !*****************************************************************************
618 subroutine advance_wrf_time(startTime, addsecs, endTime)
619   implicit none
621   integer          , intent(in)  :: addsecs
622   character (len=*), intent(in)  :: startTime
623   character (len=*), intent(out) :: endTime
624   integer :: syear,smonth,sday,shour,smin,ssec
625   integer :: days_in_month(12)
627   read(startTime,'(I4.4,1X,I2.2,1X,I2.2,1X,I2.2,1X,I2.2,1X,I2.2)') &
628        syear,smonth,sday,shour,smin,ssec
629   
630   ssec = ssec + addsecs
632   do while (ssec .ge. 60) 
633      smin = smin + 1
634      ssec = ssec - 60
635   enddo
637   do while (smin .ge. 60)
638      shour = shour + 1
639      smin = smin - 60
640   enddo
642   do while (shour .ge. 24)
643      sday = sday + 1
644      shour = shour - 24
645   enddo
648   days_in_month(1) = 31
649   if (((mod(syear,4) .eq. 0) .and. (mod(syear,100) .ne. 0)) &
650        .or. (mod(syear,400) .eq. 0)) then
651      days_in_month(2) = 29
652   else
653      days_in_month(2) = 28
654   endif
655   days_in_month(3) = 31
656   days_in_month(4) = 30
657   days_in_month(5) = 31
658   days_in_month(6) = 30
659   days_in_month(7) = 31
660   days_in_month(8) = 31
661   days_in_month(9) = 30
662   days_in_month(10) = 31
663   days_in_month(11) = 30
664   days_in_month(12) = 31
667   do while (sday .gt. days_in_month(smonth))
668      sday = sday - days_in_month(smonth)
669      smonth = smonth + 1
670      if (smonth .gt. 12) then
671         smonth = 1
672         syear = syear + 1
673      endif
674   enddo
675   
677   write(endTime,'(I4.4,A,I2.2,A,I2.2,A,I2.2,A,I2.2,A,I2.2)') &
678        syear,'-',smonth,'-',sday,'_',shour,':',smin,':',ssec
680   return
682 end subroutine