Changeset dfa7dbd in flex_extract.git for Source/Fortran/rwgrib2.f90


Ignore:
Timestamp:
Aug 21, 2019, 7:19:43 PM (5 years ago)
Author:
pesei <petra seibert -@…>
Branches:
master, ctbto, dev
Children:
5b5589b
Parents:
8c55c02
Message:

changes in the Fortran part and associated regression test

2019-08-21 PS
introduce the "new" versions of source files:

all .f90 free format
code beautification
regression test is OK

make new local gfortran makefiles, remove parameters not needed
anymore

change filenames rwgrib.f90 (all lower), preconvert to calc_etadot,
adapt messages and comments in calc_etadot.f90
adapt all makefiles to new filenames
adapt success message of logfiles in regression test references
redo regression test OK

provide softlinks for standards:

calc_etadot.out -> calc_etadot_fast.out
makefile_local_gfortran -> makefile_fast

provide changelog.txt in Fortran
provide readme.txt in Regression/FortranEtadot?

File:
1 moved

Legend:

Unmodified
Added
Removed

  • Source/Fortran/rwgrib2.f90

    rba99230 rdfa7dbd  
    1       MODULE RWGRIB2
    2 
    3       CONTAINS
    4 
    5       SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
    6 
    7       USE GRIB_API
    8 
    9       IMPLICIT NONE
    10 
    11         integer                            ::  ifile
    12         integer                            ::  iret
     1 MODULE RWGRIB2
     2
     3 CONTAINS
     4
     5 SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
     6
     7!! Read a field from GRIB file on lat-lon grid
     8
     9 USE GRIB_API
     10
     11 IMPLICIT NONE
     12
     13   integer                            ::  ifile
     14   integer                            ::  iret
    1315   integer                            ::  n,mk,parid,nm
    1416   integer                            ::  i,k
    … …  
    6870
    6971 kloop:  do k=1,nm
    70         if(parid .eq. mpar(k)) then
     72        if (parid .eq. mpar(k)) then
    7173!         l(k)=l(k)+1
    7274         feld(:,:,(k-1)*div+level)=reshape(values,(/maxl,maxb/))
    7375!         print*,(k-1)*div+l(k),parid
    7476         exit kloop
    75         endif
    76       enddo kloop
    77       if(k .gt. nm .and. parid .ne. mpar(nm)) then
     77        end if
     78      end do kloop
     79      if (k .gt. nm .and. parid .ne. mpar(nm)) then
    7880        write(*,*) k,nm,parid,mpar(nm)
    7981        write(*,*) 'ERROR readlatlon: parameter ',parid,'is not',mpar
    8082        stop
    81       endif
     83      end if
    8284
    8385!      print*,i
    8486   END DO iloop
    85    write(*,*) 'readlatlon: ',i-1,' records read'
     87!!   write(*,*) 'readlatlon: ',i-1,' records read'
    8688 
    8789   DO i=1,n
    … …  
    8991   END DO
    9092 
    91    deallocate(values)
     93   if (allocated(values)) deallocate(values)
    9294   deallocate(igrib)
    9395
    94       END SUBROUTINE READLATLON
    95 
    96       SUBROUTINE WRITELATLON(iunit,igrib,ogrib,FELD,MAXL,MAXB,MLEVEL,&
    97       MLEVELIST,MSTRIDE,MPAR)
    98 
    99       USE GRIB_API
    100 
    101       IMPLICIT NONE
    102 
    103       INTEGER IFIELD,MLEVEL,MNAUF,I,J,K,L,MSTRIDE,IERR,JOUT
    104       INTEGER MPAR(MSTRIDE),MAXL,MAXB,LEVMIN,LEVMAX
    105       INTEGER IUNIT,igrib,ogrib
    106       REAL ZSEC4(MAXL*MAXB)
    107       REAL    FELD(MAXL,MAXB,MLEVEL)
    108       CHARACTER*(*) MLEVELIST
    109       INTEGER ILEVEL(MLEVEL),MLINDEX(MLEVEL+1),LLEN
    110 
    111     ! parse MLEVELIST
    112    
    113       LLEN=len(trim(MLEVELIST))
    114       if(index(MLEVELIST,'to') .ne. 0 .or. index(MLEVELIST,'TO') .ne. 0) THEN
    115         i=index(MLEVELIST,'/')
    116         read(MLEVELIST(1:i-1),*) LEVMIN
    117         i=index(MLEVELIST,'/',.true.)
    118         read(MLEVELIST(i+1:LLEN),*) LEVMAX
    119         l=0
    120         do i=LEVMIN,LEVMAX
    121           l=l+1
    122           ILEVEL(l)=i
    123         enddo
    124       else
    125         l=1
    126         MLINDEX(1)=0
    127         do i=1,LLEN
    128           if(MLEVELIST(i:i) .eq. '/') THEN
    129             l=l+1
    130             MLINDEX(l)=i
    131           endif
    132         enddo
    133         MLINDEX(l+1)=LLEN+1
    134         do i=1,l
    135           read(MLEVELIST(MLINDEX(i)+1:MLINDEX(i+1)-1),*) ILEVEL(i)
    136         enddo
    137       endif
    138 
    139       DO k=1,l
    140         call grib_set(igrib,"level",ILEVEL(k))
    141         DO j=1,MSTRIDE
    142          call grib_set(igrib,"paramId",MPAR(j))
    143 !         if(MPAR(j) .eq. 87) then
     96   END SUBROUTINE READLATLON
     97
     98   SUBROUTINE WRITELATLON(iunit,igrib,ogrib,FELD,MAXL,MAXB,MLEVEL,&
     99   MLEVELIST,MSTRIDE,MPAR)
     100
     101!! write a field on lat-lon grid to GRIB file
     102
     103   USE GRIB_API
     104
     105   IMPLICIT NONE
     106
     107   INTEGER IFIELD,MLEVEL,MNAUF,I,J,K,L,MSTRIDE,IERR,JOUT
     108   INTEGER MPAR(MSTRIDE),MAXL,MAXB,LEVMIN,LEVMAX
     109   INTEGER IUNIT,igrib,ogrib
     110   REAL ZSEC4(MAXL*MAXB)
     111   REAL    FELD(MAXL,MAXB,MLEVEL)
     112   CHARACTER*(*) MLEVELIST
     113   INTEGER ILEVEL(MLEVEL),MLINDEX(MLEVEL+1),LLEN
     114
     115 ! parse MLEVELIST
     116
     117   LLEN=len(trim(MLEVELIST))
     118   if (index(MLEVELIST,'to') .ne. 0 .or. index(MLEVELIST,'TO') .ne. 0) THEN
     119     i=index(MLEVELIST,'/')
     120     read(MLEVELIST(1:i-1),*) LEVMIN
     121     i=index(MLEVELIST,'/',.true.)
     122     read(MLEVELIST(i+1:LLEN),*) LEVMAX
     123     l=0
     124     do i=LEVMIN,LEVMAX
     125       l=l+1
     126       ILEVEL(l)=i
     127     end do
     128   else
     129     l=1
     130     MLINDEX(1)=0
     131     do i=1,LLEN
     132       if (MLEVELIST(i:i) .eq. '/') THEN
     133                l=l+1
     134                MLINDEX(l)=i
     135              end if
     136            end do
     137            MLINDEX(l+1)=LLEN+1
     138            do i=1,l
     139              read(MLEVELIST(MLINDEX(i)+1:MLINDEX(i+1)-1),*) ILEVEL(i)
     140            end do
     141   end if
     142
     143   DO k=1,l
     144     call grib_set(igrib,"level",ILEVEL(k))
     145     DO j=1,MSTRIDE
     146       call grib_set(igrib,"paramId",MPAR(j))
     147!         if (MPAR(j) .eq. 87) then
    144148!           call grib_set(igrib,"shortName","etadot")
    145149!           call grib_set(igrib,"units","Pa,s**-1")
    146 !         endif
    147 !         if(MPAR(j) .eq. 77) then
     150!         end if
     151!         if (MPAR(j) .eq. 77) then
    148152!           call grib_set(igrib,"shortName","etadot")
    149153!           call grib_set(igrib,"units","s**-1")
    150 !         endif
    151          if(l .ne. mlevel) then
    152            zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,ILEVEL(k)),(/maxl*maxb/))
    153          else
    154            zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,k),(/maxl*maxb/))
    155          endif
    156          call grib_set(igrib,"values",zsec4)
    157 
    158          call grib_write(igrib,iunit)
    159 
    160         ENDDO
    161       ENDDO
    162 
    163 
    164 
    165       END SUBROUTINE WRITELATLON
    166 
    167       SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,&
    168         MAXLEV,MPAR,A,B)
    169 
    170       USE GRIB_API
    171 
    172       IMPLICIT NONE
    173 
    174 
    175         integer                            ::  ifile
     154!         end if
     155              if (l .ne. mlevel) then
     156         zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,ILEVEL(k)),(/maxl*maxb/))
     157              else
     158                zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,k),(/maxl*maxb/))
     159              end if
     160       call grib_set(igrib,"values",zsec4)
     161
     162       call grib_write(igrib,iunit)
     163
     164     END DO
     165   END DO
     166
     167   END SUBROUTINE WRITELATLON
     168
     169   SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,MAXLEV,MPAR,A,B)
     170
     171!!  read a GRIB file in spherical harmonics
     172
     173   USE GRIB_API
     174
     175   IMPLICIT NONE
     176
     177   integer                            ::  ifile
    176178   integer                            ::  iret
    177179   integer                            ::  n,mk,div,nm,k
    … …  
    182184   REAL :: A(MAXLEV+1),B(MAXLEV+1),pv(2*MAXLEV+2)
    183185   REAL:: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
    184 integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
    185 character*(*):: filename                             
     186  integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
     187   character*(*):: filename                             
    186188 
    187189   call grib_open_file(ifile, TRIM(FILENAME),'r')
    … …  
    226228
    227229      call grib_get(igrib(i),'paramId',parid)
    228     nm=size(mpar)
    229     div=mlevel/nm
    230     kloop:  do k=1,nm
    231         if(parid .eq. mpar(k)) then
     230      nm=size(mpar)
     231      div=mlevel/nm
     232      kloop:  do k=1,nm
     233        if (parid .eq. mpar(k)) then
    232234         l(k)=l(k)+1
    233235         cxmn(:,(k-1)*div+l(k))=values(1:(MNAUF+1)*(MNAUF+2))
    234236!         print*,(k-1)*div+l(k),parid
    235237         exit kloop
    236         endif
     238        end if
    237239       
    238       enddo kloop
    239       if(k .gt. nm .and. parid .ne. mpar(nm)) then
     240      end do kloop
     241      if (k .gt. nm .and. parid .ne. mpar(nm)) then
    240242        write(*,*) k,nm,parid,mpar(nm)
    241243        write(*,*) 'ERROR readspectral: parameter ',parid,'is not',mpar
    242244        stop
    243       endif
     245      end if
    244246
    245247!      print*,i
    … …  
    247249   END DO iloop
    248250
    249    write(*,*) 'readspectral: ',i-1,' records read'
     251!!   write(*,*) 'readspectral: ',i-1,' records read'
    250252 
    251253   DO i=1,n
    … …  
    256258   deallocate(igrib)
    257259
    258 
    259 
    260         A=pv(1:1+MAXLEV)
    261         B=pv(2+MAXLEV:2*MAXLEV+2)
    262 
    263       END SUBROUTINE READSPECTRAL
    264 
    265       END MODULE RWGRIB2
     260   A=pv(1:1+MAXLEV)
     261   B=pv(2+MAXLEV:2*MAXLEV+2)
     262
     263   END SUBROUTINE READSPECTRAL
     264
     265END MODULE RWGRIB2
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG