Index: Source/Fortran/calc_etadot.f90
===================================================================
--- Source/Fortran/calc_etadot.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/calc_etadot.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,814 @@
+PROGRAM calc_etadot
+
+!! Prepare input data for FLEXPART, esp. vertical velocity as 
+!! etadot or etadot * dp/deta
+
+!*----------------------------------------------------------------
+! author: L. Haimberger
+! date:   03/2010
+! version: V4.0                       
+!                                                                 
+!## Program calc_etadot 
+!  
+! **Prepares input data for POP model meteorological preprocessor**
+!                                                                 
+!-----------------------------------------------------------------
+!                                                                
+! Calculation of etapoint on a regular \(\lambda-\phi\) grid and writing 
+! `U,V,ETAPOINT,T,PS,Q,SD,MSL,TCC,10U, 10V, 2T,2D,LSP,CP,SSHF,SSR, 
+! EWSS,NSSS`
+! to an output file (input and output in GRIB 1 or 2 format).               
+!  etapoint is defined as the total time derivative of 
+!  ECMWF vertical coordinate eta multiplied by the derivative
+!  of pressure with respect to eta:
+!  \[\frac{\mathrm{d}\eta}{\mathrm{d}t}\frac{\partial p}{\partial \eta}\]
+!                                                                 
+!### Version history and authors:
+! - 04/1994: Leopold Haimberger, Gerhard Wotawa
+!
+! - 2003-05-11: Alexander Beck                            
+!
+! - 12/2006: L. Haimberger V2.0,
+!           handle arbitrary regular grids and T799 resolution data                    
+!
+! - 03/2010: L. Haimberger V4.0,
+!           handle GRIB edition 2 fields and T1279 resolution data     
+! - 04-06/2019: Petra Seibert, 
+!            beautify code and add FORD documentation
+!                                                                
+!-----------------------------------------------------------------
+!  #                                                              
+!## Input required:
+!                                                                
+!     UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION            
+!                                                                    
+!     11    fort.11   T,U,V           regular lambda phi grid
+!     12    fort.12   D               regular lambda phi grid   
+!     13    fort.13   LNSP            spherical harmonics
+!     14    fort.14   SD,MSL,TCC,10U,                                
+!                     10V,2T,2D       regular lambda phi grid     
+!     16    fort.16   LSP,CP,SSHF,         
+!                     SSR,EWSS,NSSS   regular lambda phi grid
+!     17    fort.17   Q               regular lambda phi grid
+!                                                                
+!------------------------------------------------------------------
+!                                                                
+!### Output produced:
+!                                                                
+!    UNIT  FILE      PARAMETER(S)      DATA REPRESENTATION            
+!                                                                   
+!    15    fort.15   `U,V,ETA,T,PS,                                  
+!                    `Q,SD,MSL,TCC,`                                  
+!                    `10U,10V,2T,2D,`  regular lambda phi grid         
+!                    `LSP,CP,SSHF,`                                   
+!                    `SSR,EWSS,NSSS`                                  
+!                                                                
+!------------------------------------------------------------------
+
+  USE PHTOGR
+  USE GRTOPH
+  USE FTRAFO
+  USE RWGRIB2
+  USE GRIB_API
+
+  IMPLICIT NONE
+
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: LNPS
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: Z
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T, UV , UV2
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: QA,OM,OMR
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: DIV, ETA,ETAR
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: DPSDL, DPSDM
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: PS,DPSDT
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: SURF,FLUX,OROLSM
+  REAL, ALLOCATABLE, DIMENSION (:) :: WSAVE,H,SINL,COSL,WSAVE2
+  REAL, ALLOCATABLE, DIMENSION (:) :: BREITE, GBREITE,AK, BK,pv
+
+! Arrays for Gaussian grid calculations
+  REAL  :: X1,X2,RMS,MW,SIG,LAM
+  REAL,ALLOCATABLE :: CUA(:,:,:),CVA(:,:,:)
+
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: P,PP,P2
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: XMN,HILFUV
+  REAL, ALLOCATABLE, DIMENSION (:) :: LNPMN,LNPMN2,LNPMN3
+  REAL, ALLOCATABLE, DIMENSION (:) :: WEIGHT
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: UGVG
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: DG, ETAG
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: GWSAVE
+  REAL, ALLOCATABLE, DIMENSION (:) :: PSG,HILF
+
+! end arrays for Gaussian grid calculations
+
+  INTEGER, ALLOCATABLE, DIMENSION (:) :: MLAT,MPSURF,MPFLUX,MPORO,MPAR
+  INTEGER, ALLOCATABLE :: GIFAX(:,:)
+
+  REAL PI,COSB,DAK,DBK,P00
+  REAL URLAR8,JMIN1,LLLAR8,MAXBMIN1,PIR8,DCOSB
+
+  INTEGER I,J,K,L,IERR,M,LTEST,MK,NGI,NGJ
+  INTEGER MFLUX,MSURF,MORO
+  INTEGER LUNIT,LUNIT2
+
+  INTEGER MAXL, MAXB, MLEVEL, LEVOUT,LEVMIN,LEVMAX
+  INTEGER MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH, MNAUF,META,METADIFF
+  INTEGER MDPDETA,METAPAR
+  REAL RLO0, RLO1, RLA0, RLA1
+  CHARACTER*300 MLEVELIST
+
+  INTEGER MAUF, MANF,IFAX(10)
+
+  INTEGER IGRIB(1),iret,ogrib
+
+  CHARACTER*80 FILENAME
+
+  NAMELIST /NAMGEN/ &
+    MAXL, MAXB, &
+    MLEVEL,MLEVELIST,MNAUF,METAPAR, &
+    RLO0, RLO1, RLA0, RLA1, &
+    MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH,META,METADIFF,&
+    MDPDETA
+
+  LTEST=1
+
+  CALL POSNAM (4,'NAMGEN')
+  READ (4,NAMGEN)
+
+  MAUF=INT(360.*(REAL(MAXL)-1.)/(RLO1-RLO0)+0.0001)
+!      PRINT*, MAUF
+
+  MANF=INT(REAL(MAUF)/360.*(360.+RLO0)+1.0001)
+  IF (MANF .gt. MAUF) MANF=MANF-MAUF
+
+
+!------------------------------------------------------------------
+!! ALLOCATE VARIABLES                       
+!------------------------------------------------------------------
+
+  ALLOCATE (LNPS(0:(MNAUF+1)*(MNAUF+2)-1,1))
+  ALLOCATE (H(0:(MNAUF+2)*(MNAUF+3)/2))
+  ALLOCATE (OM(MAXL, MAXB, MLEVEL))
+  ALLOCATE (ETA(MAXL,MAXB,MLEVEL))
+  ALLOCATE (PS(MAXL, MAXB,1),DPSDT(MAXL, MAXB,1))
+  ALLOCATE (WSAVE(4*MAUF+15),WSAVE2(4*MAUF+15))
+  ALLOCATE (BREITE(MAXB),AK(MLEVEL+1),BK(MLEVEL+1),PV(2*MLEVEL+2))
+  ALLOCATE (MPAR(2))
+  ALLOCATE (COSL(MAXL),SINL(MAXL))
+  ALLOCATE (CUA(2,4,MLEVEL),CVA(2,4,MLEVEL))
+
+!------------------------------------------------------------------
+! GAUSS STUFF                                                
+!------------------------------------------------------------------
+
+  IF (MGAUSS .EQ. 1) THEN
+    LUNIT=0
+    FILENAME='fort.18'
+
+    CALL GRIB_OPEN_FILE(LUNIT, TRIM(FILENAME),'R')
+
+    CALL GRIB_NEW_FROM_FILE(LUNIT,IGRIB(1), IRET)
+
+! we can close the file
+    CALL GRIB_CLOSE_FILE(LUNIT)
+!      call grib_get(igrib(1),'gridType', j)
+
+    NGJ=MNAUF+1
+
+    ALLOCATE (GWSAVE(8*NGJ+15,NGJ/2))
+    ALLOCATE(GIFAX(10,NGJ))
+    ALLOCATE (GBREITE(NGJ),WEIGHT(NGJ))
+    ALLOCATE (MLAT(NGJ))
+    ALLOCATE (P(0:((MNAUF+3)*(MNAUF+4))/2,NGJ/2))
+    ALLOCATE (PP(NGJ/2,0:((MNAUF+3)*(MNAUF+4))/2))
+    ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB))
+
+    CALL GRIB_GET(IGRIB(1),'numberOfPointsAlongAMeridian', NGJ)
+
+!   get as a integer
+    call grib_get(igrib(1),'pl', MLAT)
+
+    NGI=SUM(MLAT)
+
+    CALL GRIB_GET(IGRIB(1),'numberOfVerticalCoordinateValues',MK)
+
+    IF (MK/2-1 .NE. MLEVEL) THEN
+      WRITE(*,*) 'FATAL: Number of model levels',mk, &
+        ' does not agree with', MLEVEL,' in namelist'
+      STOP
+    END IF
+    call grib_get(igrib(1),'pv',pv)
+    AK=PV(1:1+MLEVEL)
+    BK=PV(2+MLEVEL:2*MLEVEL+2)
+
+    ALLOCATE (LNPMN(0:(MNAUF+1)*(MNAUF+2)-1))
+    ALLOCATE (LNPMN2(0:(MNAUF+1)*(MNAUF+2)-1))
+    ALLOCATE (UGVG(NGI, 2*MLEVEL),HILFUV(2*MAXL,2))
+    ALLOCATE (DPSDL(NGI,1),DPSDM(NGI,1))
+    ALLOCATE (PSG(NGI),HILF(NGI))
+    ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
+!      ALLOCATE (UV2(MAXL, MAXB, 2*MLEVEL))
+    ALLOCATE (XMN(0:(MNAUF+1)*(MNAUF+2)-1, 2*MLEVEL))
+    ALLOCATE (DG(NGI,MLEVEL),ETAG(NGI,MLEVEL))
+
+!! Initialisieren  Legendretransformation auf das LaT/LON Gitter  
+
+    PI=ACOS(-1.D0)
+
+!$OMP PARALLEL DO
+    DO 20 J=1,MAXB
+      BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))* PI/180.D0)
+      CALL PLGNFA(MNAUF,BREITE(J),Z(0,J))
+20  CONTINUE
+!$OMP END PARALLEL DO
+
+! Avoid possible Pole problem
+!      IF (RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
+!      IF (RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
+
+!* Initialisation of fields for FFT and Legendre transformation
+! to Gaussian grid and back to phase space
+    X1=-1.D0
+    X2=1.D0
+    CALL GAULEG(X1,X2,GBREITE,WEIGHT,NGJ)
+
+!$OMP PARALLEL DO PRIVATE(M)
+    DO J=1,NGJ/2
+      CALL PLGNFA(MNAUF,GBREITE(J),P(:,J))
+      DO M=0,(MNAUF+3)*(MNAUF+4)/2
+        PP(J,M)=P(M,J)
+      END DO
+    END DO
+!$OMP END PARALLEL DO
+
+!       MPAR(1)=152
+    FILENAME='fort.12' 
+!!  read LNSP in SH
+    CALL READSPECTRAL(FILENAME,LNPMN,MNAUF,1,MLEVEL,(/152/),AK,BK)
+    CALL SET99(WSAVE,IFAX,mauf)
+    CALL PHGCUT(LNPMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+    CALL STATIS(MAXL,MAXB,1,EXP(PS),RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS PS: ',RMS,MW,SIG
+
+    DO J=1,NGJ/2
+      CALL SET99(GWSAVE(1,J),GIFAX(1,J),MLAT(J))
+    END DO
+    CALL PHGR213(LNPMN,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
+    PSG=HILF
+    CALL GRPH213(LNPMN2,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
+    CALL PHGR213(LNPMN2,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
+
+
+    HILF=exp(PSG)-exp(HILF)
+
+    CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS ratio of PS: ',RMS,MW,SIG
+
+    PSG=EXP(PSG)
+    HILF=PSG
+    CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS PSG: ',RMS,MW,SIG
+
+111 FILENAME='fort.10'
+!!  read u,v in SH
+    CALL READSPECTRAL(FILENAME,XMN,MNAUF,2*MLEVEL,MLEVEL,(/131,132/),AK,BK)
+
+!!  Transformieren des Windes auf das Gaussgitter  
+    CALL PHGR213(XMN,UGVG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,2*MLEVEL)
+    DO K=1,MLEVEL
+! North Pole
+      CALL JSPPOLE(XMN(:,K),1,MNAUF,.TRUE.,CUA(:,:,K))
+      CALL JSPPOLE(XMN(:,MLEVEL+K),1,MNAUF,.TRUE.,CVA(:,:,K))
+! South Pole
+      CALL JSPPOLE(XMN(:,K),-1,MNAUF,.TRUE.,CUA(:,3:4,K))
+      CALL JSPPOLE(XMN(:,MLEVEL+K),-1,MNAUF,.TRUE.,CVA(:,3:4,K))
+    END DO
+
+    DO K=1,2*MLEVEL
+      IF (MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+    END DO
+    CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
+
+
+112 FILENAME='fort.13'
+!!  read DIV in SH
+    CALL READSPECTRAL(FILENAME,XMN,MNAUF,MLEVEL,MLEVEL,(/155/),AK,BK)
+!! Transformieren der horizontalen Divergenz auf das Gaussgitter  
+    CALL PHGR213(XMN,DG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,MLEVEL)
+    CALL STATIS(MAXL,MAXB,1,DG,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3p,3F12.4)') 'STATISTICS DG-PS: ',RMS,MW,SIG
+
+!! Berechnung des Gradienten des Logarithmus des Bodendrucks auf dem Gaussgitter  
+    CALL PHGRAD(LNPMN,DPSDL,DPSDM,GWSAVE,GIFAX,P,H,MLAT,MNAUF,NGI,NGJ,1)
+
+!! Berechnung der Vertikalgeschwindigkeit auf dem Gaussgitter  
+    CALL CONTGL(HILF,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
+      GBREITE,ETAG,MLAT,AK,BK,NGI,NGJ,MLEVEL)
+! note that HILF is ps on input and  dpsdt*ps on output
+      
+     CALL STATIS(MAXL,MAXB,1,ETAG,RMS,MW,SIG)
+     WRITE(*,'(A,T20,3p,3F12.4)') 'STATISTICS ETAG-PS: ',RMS,MW,SIG
+    CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,MLEVEL)
+     CALL STATIS(MAXL,MAXB,1,ETAG,RMS,MW,SIG)
+     WRITE(*,'(A,T20,3p,3F12.4)') 'STATISTICS ETAG-PS: ',RMS,MW,SIG
+    DO K=1,MLEVEL
+      IF (MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+    END DO
+
+    CALL PHGCUT(XMN,ETA,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+    CALL STATIS(MAXL,MAXB,1,ETA,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3p,3F12.4)') 'STATISTICS ETA-PS: ',RMS,MW,SIG
+
+    CALL GRPH213(XMN,HILF,GWSAVE,GIFAX,PP,WEIGHT,MLAT, MNAUF,NGI,NGJ,1)
+    CALL STATIS(MAXL,MAXB,1,HILF,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3p,3F12.4)') 'STATISTICS HILF-PS: ',RMS,MW,SIG
+
+    IF (MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,1),MNAUF,MSMOOTH)
+    CALL PHGCUT(XMN,DPSDT,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+
+    CALL STATIS(MAXL,MAXB,1,DPSDT,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS DPSDT: ',RMS,MW,SIG
+
+    IF (MOMEGADIFF .ne. 0) THEN
+!! Berechnung von Omega auf dem Gaussgitter  
+
+      CALL OMEGA(PSG,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
+        GBREITE,ETAG,MLAT,AK,BK,NGI ,NGJ,MLEVEL)
+      CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,MLEVEL)
+      DO K=1,MLEVEL
+        IF (MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+      END DO
+      CALL PHGCUT(XMN,OM,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+    END IF ! MOMEGA
+
+    CALL GRPH213(XMN,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
+
+    CALL STATIS(MAXL,MAXB,1,PSG,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS PSG-PS: ',RMS,MW,SIG
+
+    CALL PHGCUT(XMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+
+    CALL STATIS(MAXL,MAXB,1,PS,RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS PS: ',RMS,MW,SIG
+
+114 DEALLOCATE(HILF,PSG,DPSDL,DPSDM,ETAG,DG,LNPMN)
+
+!      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
+! CALL GRPH213(XMN,UGVG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,
+!     *MNAUF,NGI,NGJ,2*MLEVEL)
+!        DO K=1,2*MLEVEL
+!          IF (MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+!        END DO
+!        CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z,
+!     *MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
+    DEALLOCATE(PP,P,UGVG,MLAT,GBREITE,WEIGHT,GWSAVE,XMN)
+!        CALL ETAGAUSS(Z,WSAVE
+!     *,BREITE,UV,ETA,OM,PS,
+!     *MAUF,MAXB,MAXL,MANF,MNAUF,MLEVEL,MSMOOTH)
+
+  ELSE
+
+!-----------------------------------------------------------------
+!     READING OF PREPARED METEOROLOGICAL FIELDS             
+!                                                                
+!     THE FOLLOWING FIELDS ARE EXPECTED:                    
+!                                                                
+!     UNIT 11: T,U,V        (REGULAR GRID)                  
+!     UNIT 17: Q            (REGULAR GRID)                  
+!     UNIT 13: D            (REGULAR GRID)                  
+!     UNIT 12: LNSP         (SPHERICAL HARMONICS)           
+!     UNIT 14: SURFACE DATA (REGULAR GRID)                  
+!     UNIT 16: FLUX DATA    (REGULAR GRID)                  
+!------------------------------------------------------------------
+
+    ALLOCATE (MLAT(MAXB))
+    MLAT=MAXL
+    ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,1))
+    ALLOCATE (DPSDL(MAXL,MAXB),DPSDM(MAXL,MAXB))
+    ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL),DIV(MAXL,MAXB,MLEVEL))
+
+!------------------------------------------------------------------
+!! READING OF SURFACE PRESSURE                   
+!------------------------------------------------------------------
+
+    FILENAME='fort.12'
+    CALL READSPECTRAL(FILENAME,LNPS,MNAUF,1,MLEVEL,(/152/),AK,BK)
+
+!------------------------------------------------------------------
+!! READING OF U,V                     
+!------------------------------------------------------------------
+
+! OPENING OF UNBLOCKED GRIB FILE
+
+    FILENAME='fort.10'
+    CALL READLATLON(FILENAME,UV,MAXL,MAXB,2*MLEVEL,(/131,132/))
+
+    PI=ACOS(-1.D0)
+    DO J=1,MAXB
+      BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))*PI/180.D0)
+    END DO
+
+! Avoid possible Pole problem
+!      IF (RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
+!      IF (RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
+
+    DO K=1,2*MLEVEL
+      DO J=1,MAXB
+        COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
+        IF (RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
+            RLA1 .EQ.  90.0 .AND. J .EQ. 1) THEN
+          UV(:,J,K)=UV(:,J,K)/1.D6
+        ELSE
+          UV(:,J,K)=UV(:,J,K)*COSB
+        END IF
+      END DO
+    END DO
+
+!------------------------------------------------------------------
+!! READING OF LNSP on grid               
+!------------------------------------------------------------------
+
+! For debugging only
+!      FILENAME='LNSPG_G.20060330.600'
+!      INQUIRE(FILE=FILENAME,EXIST=EX)
+!      CALL READLATLON(FILENAME,QA,
+!     *MAXL,MAXB,1,1,(/152/))
+
+!------------------------------------------------------------------
+!! READING OF DIVERGENCE                      
+!------------------------------------------------------------------
+
+    IF (META .EQ. 0 .OR. METADIFF .EQ. 1) THEN
+      FILENAME='fort.13'
+      CALL READLATLON(FILENAME,DIV,MAXL,MAXB,MLEVEL,(/155/))
+    END IF
+
+
+!------------------------------------------------------------------
+!
+!  Calculation of etapoint --> total time derivative of    
+!  ECMWF vertical coordinate eta multiplied by the derivative
+!  of pressure with respect to eta:
+!  \[\frac{\mathrm{d}\eta}{\mathrm{d}t}\frac{\partial p}{\partial \eta}\]
+!------------------------------------------------------------------
+
+!* Initialisieren  Legendretransformation auf das LaT/LON Gitter
+!! Without Gaussian grid calculation Legendre Polynomials are calculated
+!! only for one latitude to save space
+
+
+
+    DO J=1,MAXB
+      CALL PLGNFA(MNAUF,BREITE(J),Z(0,1))
+      CALL PHGCUT(LNPS,PS(:,J,1),WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,1,1)
+      IF (META .EQ. 0 .OR. METADIFF .EQ. 1 ) THEN
+        CALL PHGRACUT(LNPS,DPSDL(:,J),DPSDM(:,J),WSAVE,IFAX,Z,H,MAUF, &
+          MNAUF,MAXL,1,MANF,1)
+      END IF
+    END DO
+
+    PS=EXP(PS)
+
+! For debugging only
+    CALL STATIS(MAXL,MAXB,1,PS(:,:,1),RMS,MW,SIG)
+    WRITE(*,'(A,T20,3F12.4)') 'STATISTICS: ',RMS,MW,SIG
+
+    IF (MOMEGADIFF .ne. 0) THEN
+      CALL OMEGA(PS,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
+        BREITE,OM,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL)
+    END IF
+
+    IF (META .EQ. 0 .OR. METADIFF .ne. 0) THEN
+      DPSDT=PS
+      CALL CONTGL(DPSDT,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
+        BREITE,ETA,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL)
+    END IF
+
+  END IF ! MGAUSS
+
+!! CREATE FILE VERTICAL.EC NEEDED BY POP MODEL  
+
+  OPEN(21,FILE='VERTICAL.EC')
+  WRITE(21,'(A)')
+  WRITE(21,'(A)') 'VERTICAL DISCRETIZATION OF POP MODEL'
+  WRITE(21,'(A)')
+  write(21,'(i3,a)') MLEVEL,'   number of layers'
+  WRITE(21,'(A)')
+  WRITE(21,'(A)') '* A(NLEV+1)'
+  WRITE(21,'(A)')
+  DO 205 I=1,MLEVEL+1
+205 WRITE(21,'(F18.12)') AK(I)
+  WRITE(21,'(A)')
+  WRITE(21,'(A)') '* B(NLEV+1)'
+  WRITE(21,'(A)')
+  DO 210 I=1,MLEVEL+1
+210 WRITE(21,'(F18.12)') BK(I)
+  CLOSE(21)
+
+!------------------------------------------------------------------
+! READING OF OMEGA                           
+!------------------------------------------------------------------
+
+  IF (MOMEGA .NE. 0 ) THEN
+
+    ALLOCATE (OMR(MAXL, MAXB, MLEVEL))
+
+    FILENAME='fort.19'
+    CALL READLATLON(FILENAME,OMR,MAXL,MAXB,MLEVEL,(/135/))
+
+    IF (MOMEGADIFF .NE. 0) THEN
+
+      DO K=1,MLEVEL
+        CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') '       ETA: ',K,RMS,MW,SIG
+        CALL STATIS(MAXL,MAXB,1,OMR(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') '     OMEGA: ',K,RMS,MW,SIG
+        CALL STATIS(MAXL,MAXB,1,OM(:,:,K)-OMR(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') 'OMEGA DIFF: ',K,RMS,MW,SIG
+      END DO
+
+    END IF
+  END IF
+
+!------------------------------------------------------------------
+! READING OF ETA                             
+!------------------------------------------------------------------
+
+  IF (META .NE. 0 ) THEN
+
+    ALLOCATE (ETAR(MAXL, MAXB, MLEVEL))
+
+    P00=101325.
+    FILENAME='fort.21'
+    CALL READLATLON(FILENAME,ETAR,MAXL,MAXB,MLEVEL,(/77/))
+
+    IF(MDPDETA .EQ. 1) THEN
+      DO K=1,MLEVEL
+        DAK=AK(K+1)-AK(K)
+        DBK=BK(K+1)-BK(K)
+        DO J=1,MAXB
+          DO I=1,MAXL
+            ETAR(I,J,K)=2*ETAR(I,J,K)*PS(I,J,1)*(DAK/PS(I,J,1)+DBK)/ &
+              (DAK/P00+DBK)
+            IF (K .GT. 1) ETAR(I,J,K)=ETAR(I,J,K)-ETAR(I,J,K-1)
+          END DO
+        END DO
+      END DO
+    END IF
+
+    IF (METADIFF .NE. 0 ) THEN
+
+      DO K=1,MLEVEL
+        CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') '       ETA: ',K,RMS,MW,SIG
+        CALL STATIS(MAXL,MAXB,1,ETAR(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') '     ETAR: ',K,RMS,MW,SIG
+        CALL STATIS(MAXL,MAXB,1,ETA(:,:,K)-ETAR(:,:,K),RMS,MW,SIG)
+        WRITE(*,'(A12,I3,3F12.4)') 'ETA DIFF: ',K,RMS,MW,SIG
+      END DO
+      DO K=1,MLEVEL
+        WRITE(*,'(I3,2F12.4)') K,ETA(1,MAXB/2,K),ETAR(1,MAXB/2,K)
+      END DO
+    ELSE
+      ETA=ETAR
+    END IF
+  END IF
+
+  ALLOCATE (T(MAXL, MAXB, MLEVEL))
+  ALLOCATE (QA(MAXL, MAXB, MLEVEL))
+
+!------------------------------------------------------------------
+!! READING OF T                     
+!------------------------------------------------------------------
+
+! OPENING OF UNBLOCKED GRIB FILE
+
+  FILENAME='fort.11'
+  CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/))
+
+!------------------------------------------------------------------
+!! READING OF SPECIFIC HUMIDITY               
+!------------------------------------------------------------------
+
+  FILENAME='fort.17'
+  CALL READLATLON(FILENAME,QA,MAXL,MAXB,MLEVEL,(/133/))
+
+!------------------------------------------------------------------
+!                     TEST READING OF UV from MARS (debug only)  
+!------------------------------------------------------------------
+!      FILENAME='fort.22'
+!      CALL READLATLON(FILENAME,UV2,MAXL,MAXB,2*MLEVEL,2,(/131,132/))
+
+!------------------------------------------------------------------
+!! WRITE MODEL LEVEL DATA TO fort.15           
+!------------------------------------------------------------------
+
+!!     Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi)
+!!     Now we are transforming back to the usual winds.
+
+  DO K=1,MLEVEL
+    DO J=2,MAXB-1
+      COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
+      UV(:,J,K)=UV(:,J,K)/COSB
+      UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
+    END DO
+
+! special treatment for poles, if necessary.
+    DO J=1,MAXB,MAXB-1
+      COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
+      IF (1.0-BREITE(J)*BREITE(J) .GT. 0 .OR. MGAUSS .NE. 1) THEN
+        IF (RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
+            RLA1 .EQ.  90.0 .AND. J .EQ. 1) THEN
+          UV(:,J,K)=UV(:,J,K)*1.D6
+          UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6
+        ELSE
+          UV(:,J,K)=UV(:,J,K)/COSB
+          UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
+        END IF
+      ELSE
+        HILFUV(5:MAXL,:)=0.
+        HILFUV(1:2,:)=0.
+        IF (J.EQ.MAXB) THEN
+! Suedpol
+          HILFUV(3:4,1)=CUA(:,4,K)
+          HILFUV(3:4,2)=CVA(:,4,K)
+        ELSE
+! Nordpol
+          HILFUV(3:4,1)=CUA(:,2,K)
+          HILFUV(3:4,2)=CVA(:,2,K)
+        END IF
+        CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
+        DO I=0,MAXL-1
+          IF (MANF+I .LE. MAXL) THEN
+            UV(I+1,J,K)=HILFUV(MANF+I,1)
+          ELSE
+            UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1)
+          END IF
+        END DO
+        CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
+        DO I=0,MAXL-1
+          IF (MANF+I .LE. MAXL) THEN
+            UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2)
+          ELSE
+            UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2)
+          END IF
+        END DO
+      end if
+    END DO
+  END DO
+
+! open output file
+  call grib_open_file(LUNIT,'fort.15','w')
+
+! we use temperature on lat/lon on model levels as template for model level data
+  LUNIT2=0
+  CALL GRIB_OPEN_FILE(LUNIT2,'fort.11','R')
+  CALL GRIB_NEW_FROM_FILE(LUNIT2,IGRIB(1), IRET)
+  CALL GRIB_CLOSE_FILE(LUNIT2)
+
+
+  CALL WRITELATLON &
+    (LUNIT,IGRIB(1),OGRIB,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/131/))
+
+  CALL WRITELATLON &
+    (LUNIT,IGRIB(1),OGRIB,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/132/))
+
+  IF (MDPDETA .ne. 1 .AND. MGAUSS .EQ. 0 .and. META .eq. 1) THEN
+    CALL WRITELATLON &
+      (LUNIT,IGRIB(1),OGRIB,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/77/))
+  ELSE
+    CALL WRITELATLON &
+      (LUNIT,IGRIB(1),OGRIB,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/METAPAR/))
+  END IF
+
+  CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,T,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/130/))
+
+  CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,PS,MAXL,MAXB,1,'1',1,(/134/))
+
+  CALL GRIB_SET(IGRIB(1),"levelType","ml")
+  CALL GRIB_SET(IGRIB(1),"typeOfLevel","hybrid")
+  CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,QA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/133/))
+
+
+  IF (MOMEGA .EQ. 1) THEN
+    CALL GRIB_OPEN_FILE(LUNIT2,'fort.25','w')
+    CALL WRITELATLON  & 
+      (LUNIT2,IGRIB(1),OGRIB,OMR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
+
+    IF (MOMEGADIFF .EQ. 1) THEN
+      CALL WRITELATLON(LUNIT2,IGRIB(1),OGRIB,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
+      OM=OM-OMR
+      CALL WRITELATLON &
+        (LUNIT2,IGRIB(1),OGRIB,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
+      CALL GRIB_CLOSE_FILE(LUNIT2)
+    END IF
+  END IF
+
+  IF (META .EQ. 1 .AND. METADIFF .EQ. 1) THEN
+    CALL GRIB_OPEN_FILE(LUNIT2,'fort.26','w')
+    CALL WRITELATLON &
+      (LUNIT2,IGRIB(1),OGRIB,ETAR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
+!        IF (MOMEGADIFF .EQ. 1) THEN
+    CALL WRITELATLON(LUNIT2,IGRIB(1),OGRIB,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
+    OM=ETA-ETAR
+    CALL WRITELATLON &
+      (LUNIT2,IGRIB(1),OGRIB,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
+    CALL GRIB_CLOSE_FILE(LUNIT2)
+!        END IF
+  END IF
+
+  CALL GRIB_CLOSE_FILE(LUNIT)
+
+2000 STOP 'SUCCESSFULLY FINISHED calc_etadot: CONGRATULATIONS'
+3000 STOP 'ROUTINE calc_etadot: ERROR'
+9999 stop 'ROUTINE calc_etadot: ERROR'
+
+END
+
+!------------------------------------------------------------------
+INTEGER FUNCTION IA (FIELD1,NI,NJ,NK,G)
+
+
+!------------------------------------------------------------------
+!! Calculate something that is roughly log10( maxval(field1)/g ) [PS]           
+!------------------------------------------------------------------
+
+
+  IMPLICIT NONE
+
+  INTEGER :: I,J,K
+  INTEGER, INTENT(IN) :: NI,NJ,NK
+  REAL, INTENT(IN) :: FIELD1(NI,NJ,NK)
+  REAL, INTENT(IN) :: G
+  REAL  :: RMIN,RMAX,XMAX,A,A1,A2
+
+  RMAX=FIELD1(1,1,1)
+  RMIN=FIELD1(1,1,1)
+
+  DO 100 K=1,NK
+    DO 100 J=1,NJ
+      DO 100 I=1,NI
+        IF (FIELD1(I,J,K) .GT. RMAX) RMAX=FIELD1(I,J,K)
+        IF (FIELD1(I,J,K) .LT. RMIN) RMIN=FIELD1(I,J,K)
+100 CONTINUE
+
+  IF (ABS(RMIN) .GT. RMAX .OR. ABS(RMIN) .EQ. RMAX) THEN
+    XMAX=ABS(RMIN)
+  ELSE
+    XMAX=RMAX
+  END IF
+
+  IF (XMAX .EQ. 0) THEN
+    IA = 0
+    RETURN
+  END IF
+
+  A1=LOG10( (G/10.d0)/XMAX )
+  A2=LOG10( G/XMAX )
+  IF (A1 .gt. A2) THEN
+    A=A2
+  ELSE
+    A=A1
+  END IF
+
+  IF (A .GT. 0) IA=INT(A)
+  IF (A .LT. 0) IA=INT(A-1.0)
+
+  RETURN
+
+END
+
+SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG)
+
+!------------------------------------------------------------------
+!! calculate mean, rms, stdev
+!------------------------------------------------------------------
+
+  IMPLICIT REAL (A-H,O-Z)
+
+  REAL PHI(NI,NJ,NK),SIG,MW,RMS,P
+
+  N=NI*NJ*NK
+
+  RMS=0.
+  MW=0.
+! 10.86 sinstead of 11.04 sec
+      DO 10 K=1,NK
+    DO 10 J=1,NJ
+  DO 10 I=1,NI
+        P=PHI(I,J,K)
+        RMS=RMS+P*P
+        MW=MW+P
+10 CONTINUE
+
+  RMS=SQRT(RMS/N)
+  MW=MW/N
+
+  IF (RMS*RMS-MW*MW .LT. 0.) THEN
+    SIG=0.0
+  ELSE
+    SIG=SQRT(RMS*RMS-MW*MW)
+  END IF
+
+  RETURN
+
+END
Index: Source/Fortran/calc_etadot.out
===================================================================
--- Source/Fortran/calc_etadot.out	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/calc_etadot.out	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,1 @@
+calc_etadot_fast.out
Index: Source/Fortran/changelog.txt
===================================================================
--- Source/Fortran/changelog.txt	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/changelog.txt	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,25 @@
+2019-08-19 PS
+rwGRIB2.f90 subroutine readlatlon: check alloc status deallocate values
+  - in case of empty input file 
+  NOTE: check why we run into that at all!
+  Improve by realloc only if dim has changed
+
+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
+   
+   
Index: urce/Fortran/ftrafo.f
===================================================================
--- Source/Fortran/ftrafo.f	(revision None)
+++ 	(revision )
@@ -1,504 +1,0 @@
-	MODULE FTRAFO
-
-	CONTAINS
-
-
-
-C
-C	Implementierung der spektralen Transformationsmethode unter Verwendung
-C	des reduzierten Gauss'schen Gitters
-C	
-C Berechnung der scale winds aus Vorticity und Divergenz
-C uebergibt man in XMN die Divergenz, so wird der divergente Anteil des
-C Windes (XPHI=Ud,XPHI=Vd) zurueckgegeben, uebergibt man die Vorticity, so 
-C erhaelt man den rotationellen Wind (XLAM=Vrot,XPHI=-Urot).
-C Summiert man beide, erhaelt man den gesamten Scale wind
-C	GWSAVE ist ein Hilfsfeld fuer die FFT
-C	P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
-C	MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
-C	MNAUF gibt die spektrale Aufloesung an,
-C	NI = Anzahl der Gauss'schen Gitterpunkte pro Flaeche
-C	NJ = Anzahl der Gauss'schen Breiten,
-C	NK = Anzahl der Niveaus
- 
-	SUBROUTINE VDTOUV(XMN,XLAM,XPHI,GWSAVE,IFAX,P,MLAT,MNAUF,NI,NJ,NK)
-
-
-	USE PHTOGR
-
-	IMPLICIT NONE
-	INTEGER   J,N,NI,NJ,NK,MNAUF,GGIND(NJ/2)
-	INTEGER		MLAT(NJ),IFAX(10,NJ)
-	REAL			XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
-	REAL		P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
-	REAL		H(0:(MNAUF+2)*(MNAUF+3)/2)
-	REAL			XLAM(NI,NK),XPHI(NI,NK)
-	REAL			GWSAVE(8*NJ+15,NJ/2)
-	REAL 		SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
-	REAL	  RT,IT
-
-	GGIND(1)=0
-	DO 4 J = 2,NJ/2	   	   
-	  GGIND(J)=GGIND(J-1)+MLAT(J-1)
-4	CONTINUE
-!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
-	DO 5 J = 1,NJ/2
-	   CALL VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,GGIND(J),
-     *MLAT,MNAUF,NI,NJ,NK)   
- 5	CONTINUE
-!$OMP END PARALLEL DO
-	RETURN
-	END SUBROUTINE VDTOUV
-
-	SUBROUTINE VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,
-     *GGIND,MLAT,MNAUF,NI,NJ,NK) 
-
-	USE PHTOGR
-
-	IMPLICIT NONE
-	INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
-	INTEGER		MLAT(NJ),IFAX(10,NJ)
-	REAL     	UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
-	REAL 			VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
-	REAL			XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
-	REAL		P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
-	REAL		H(0:(MNAUF+2)*(MNAUF+3)/2)
-	REAL			XLAM(NI,NK),XPHI(NI,NK)
-	REAL			GWSAVE(8*NJ+15,NJ/2)
-	REAL 		ERAD,SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
-	REAL	  FAC(0:MNAUF),RT,IT
-  
-	
-	ERAD = 6367470.D0
-	
-	FAC(0)=0.D0
-	DO 12 N=1,MNAUF
-	  FAC(N)=-ERAD/DBLE(N)/DBLE(N+1)
-12	CONTINUE
-
-	 CALL DPLGND(MNAUF,P(0,J),H)
-	 DO 3 K = 1,NK
-	  LL=0
-	  LLP=0
-	  LLH=0
-	  DO 2 M = 0,MNAUF
-	   SCR=0.D0
-	   SCI=0.D0
-	   ACR=0.D0
-	   ACI=0.D0
-	   MUSCR=0.D0
-	   MUSCI=0.D0
-	   MUACR=0.D0
-	   MUACI=0.D0
-	   LLS=LL
-	   LLPS=LLP
-	   LLHS=LLH
-	   IF(2*M+1.LT.MLAT(J)) THEN
-	      DO 1 N = M,MNAUF,2
-	       RT=XMN(2*LL,K)*FAC(N)
-	       IT=XMN(2*LL+1,K)*FAC(N)
-	       SCR =SCR+ RT*P(LLP,J)
-	       SCI =SCI+ IT*P(LLP,J)
-	       MUACR =MUACR+ RT*H(LLH)
-	       MUACI =MUACI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 1	      CONTINUE
-	      LL=LLS+1
-	      LLP=LLPS+1
-	      LLH=LLHS+1
-	      DO 11 N = M+1,MNAUF,2
-	       RT=XMN(2*LL,K)*FAC(N)
-	       IT=XMN(2*LL+1,K)*FAC(N)
-	       ACR =ACR+ RT*P(LLP,J)
-	       ACI =ACI+ IT*P(LLP,J)
-	       MUSCR =MUSCR+ RT*H(LLH)
-	       MUSCI =MUSCI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 11	      CONTINUE
-	   ENDIF
-	   LL=LLS+(MNAUF-M+1)
-	   LLP=LLPS+(MNAUF-M+3)
-	   LLH=LLHS+(MNAUF-M+2)
-
-	   UFOUC(2*M)=-M*(SCI-ACI)
-	   UFOUC(2*M+1)=M*(SCR-ACR)
-	   VFOUC(2*M)=-M*(SCI+ACI)
-	   VFOUC(2*M+1)=M*(SCR+ACR)
-	   
-	   MUFOUC(2*M)=-(MUSCR-MUACR)
-	   MUFOUC(2*M+1)=-(MUSCI-MUACI)
-	   MVFOUC(2*M)=-(MUSCR+MUACR)
-	   MVFOUC(2*M+1)=-(MUSCI+MUACI)
- 2	  CONTINUE
-			
-	 CALL RFOURTR(VFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
-	 CALL RFOURTR(UFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
-			
-	 CALL RFOURTR(MVFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
-	 CALL RFOURTR(MUFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
-
-3		CONTINUE
-
-      RETURN
-      END SUBROUTINE VDUVSUB
-
-C Berechnung des Gradienten eines Skalars aus dem Feld des
-C	Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
-C	Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
-C	GWSAVE ist ein Hilfsfeld fuer die FFT
-C	P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
-C	MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
-C	MNAUF gibt die spektrale Aufloesung an, 
-C	NI = Anzahl der Gauss'schen Gitterpunkte,
-C	NJ = Anzahl der Gauss'schen Breiten,
-C	NK = Anzahl der Niveaus
- 
-	SUBROUTINE PHGRAD(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MLAT,
-     *MNAUF,NI,NJ,NK)
-
-	USE PHTOGR
-	IMPLICIT NONE
-	INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
-	INTEGER		MLAT(NJ),IFAX(10,NJ)
-	REAL     	UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
-	REAL 			VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
-	REAL			XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
-	REAL		P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
-	REAL		H(0:(MNAUF+2)*(MNAUF+3)/2)
-	REAL			XLAM(NI,NK),XPHI(NI,NK)
-	REAL			GWSAVE(8*NJ+15,NJ/2)
-	REAL      ERAD
-	REAL 		SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
-	
-	ERAD = 6367470.0
-	
-	GGIND=0
-	DO 4 J = 1,NJ/2
-	 CALL DPLGND(MNAUF,P(0,J),H)
-	 DO 3 K = 1,NK
-	  LL=0
-	  LLP=0
-	  LLH=0
-	  DO 2 M = 0,MNAUF
-	   SCR=0.D0
-	   SCI=0.D0
-	   ACR=0.D0
-	   ACI=0.D0
-	   MUSCR=0.D0
-	   MUSCI=0.D0
-	   MUACR=0.D0
-	   MUACI=0.D0
-	   LLS=LL
-	   LLPS=LLP
-	   LLHS=LLH
-	   IF(2*M+1.LT.MLAT(J)) THEN
-	      DO 1 N = M,MNAUF,2
-	       RT=XMN(2*LL,K)
-	       IT=XMN(2*LL+1,K)
-	       SCR =SCR+ RT*P(LLP,J)
-	       SCI =SCI+ IT*P(LLP,J)
-	       MUACR =MUACR+RT*H(LLH)
-	       MUACI =MUACI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 1	      CONTINUE
-	      LL=LLS+1
-	      LLP=LLPS+1
-	      LLH=LLHS+1
-	      DO 11 N = M+1,MNAUF,2
-	       RT=XMN(2*LL,K)
-	       IT=XMN(2*LL+1,K)
-	       ACR =ACR+ RT*P(LLP,J)
-	       ACI =ACI+ IT*P(LLP,J)
-	       MUSCR =MUSCR+ RT*H(LLH)
-	       MUSCI =MUSCI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 11	      CONTINUE
-	   ENDIF
-          LL=LLS+(MNAUF-M+1)
-          LLP=LLPS+(MNAUF-M+3)
-          LLH=LLHS+(MNAUF-M+2)
-
-		UFOUC(2*M)=-M*(SCI-ACI)/ERAD
-		UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
-		VFOUC(2*M)=-M*(SCI+ACI)/ERAD
-		VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
-
-		MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
-		MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
-		MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
-		MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
-2	CONTINUE
-
-	 CALL RFOURTR(VFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
-	 CALL RFOURTR(UFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
-			
-	 CALL RFOURTR(MVFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
-	 CALL RFOURTR(MUFOUC,
-     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
-         XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
-			
-3		CONTINUE
-		GGIND=GGIND+MLAT(J)
-4	CONTINUE
-
-
-	RETURN
-	END SUBROUTINE PHGRAD
-
-C Berechnung des Gradienten eines Skalars aus dem Feld des
-C	Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
-C	Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
-C	GWSAVE ist ein Hilfsfeld fuer die FFT
-C	P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
-C	MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
-C	MNAUF gibt die spektrale Aufloesung an, 
-C	NI = Anzahl der Gauss'schen Gitterpunkte,
-C	NJ = Anzahl der Gauss'schen Breiten,
-C	NK = Anzahl der Niveaus
- 
-	SUBROUTINE PHGRACUT(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MAUF,
-     *MNAUF,NI,NJ,MANF,NK)
-
-	USE PHTOGR
-	IMPLICIT NONE
-	INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
-	INTEGER		MAUF,MANF,I,IFAX(10)
-	REAL     	UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
-	REAL 			VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
-	REAL			XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
-	REAL		P(0:(MNAUF+3)*(MNAUF+4)/2,NJ)
-	REAL		H(0:(MNAUF+2)*(MNAUF+3)/2)
-	REAL			XLAM(NI,NJ,NK),XPHI(NI,NJ,NK)
-	REAL			HLAM(MAXAUF,2),HPHI(MAXAUF,2)
-	REAL			GWSAVE(4*MAUF+15)
-	REAL      ERAD
-	REAL 		SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
-	
-	ERAD = 6367470.0
-	
-	GGIND=0
-	DO 4 J = 1,NJ
-	 CALL DPLGND(MNAUF,P(0,J),H)
-	 DO 3 K = 1,NK
-	  LL=0
-	  LLP=0
-	  LLH=0
-	  DO 2 M = 0,MNAUF
-	   SCR=0.D0
-	   SCI=0.D0
-	   ACR=0.D0
-	   ACI=0.D0
-	   MUSCR=0.D0
-	   MUSCI=0.D0
-	   MUACR=0.D0
-	   MUACI=0.D0
-	   LLS=LL
-	   LLPS=LLP
-	   LLHS=LLH
-	   IF(2*M+1.LT.MAUF) THEN
-	      DO 1 N = M,MNAUF,2
-	       RT=XMN(2*LL,K)
-	       IT=XMN(2*LL+1,K)
-	       SCR =SCR+ RT*P(LLP,J)
-	       SCI =SCI+ IT*P(LLP,J)
-	       MUACR =MUACR+RT*H(LLH)
-	       MUACI =MUACI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 1	      CONTINUE
-	      LL=LLS+1
-	      LLP=LLPS+1
-	      LLH=LLHS+1
-	      DO 11 N = M+1,MNAUF,2
-	       RT=XMN(2*LL,K)
-	       IT=XMN(2*LL+1,K)
-	       ACR =ACR+ RT*P(LLP,J)
-	       ACI =ACI+ IT*P(LLP,J)
-	       MUSCR =MUSCR+ RT*H(LLH)
-	       MUSCI =MUSCI+ IT*H(LLH)
-	       LL=LL+2
-	       LLP=LLP+2
-	       LLH=LLH+2
- 11	      CONTINUE
-	   ENDIF
-          LL=LLS+(MNAUF-M+1)
-          LLP=LLPS+(MNAUF-M+3)
-          LLH=LLHS+(MNAUF-M+2)
-
-		UFOUC(2*M)=-M*(SCI-ACI)/ERAD
-		UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
-		VFOUC(2*M)=-M*(SCI+ACI)/ERAD
-		VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
-
-		MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
-		MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
-		MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
-		MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
-2			CONTINUE
-			
-		CALL RFOURTR(VFOUC,
-     *GWSAVE,IFAX,MNAUF,MAUF,1)
-			
-		CALL RFOURTR(MVFOUC,
-     *GWSAVE,IFAX,MNAUF,MAUF,1)
-			
-            DO 6 I=0,NI-1
-                  IF(MANF+I.LE. MAUF) THEN
-                    XLAM(I+1,J,K)=VFOUC(MANF+I-1)
-                    XPHI(I+1,J,K)=MVFOUC(MANF+I-1)
-                  ELSE
-                    XLAM(I+1,J,K)=VFOUC(MANF-MAUF+I-1)
-                    XPHI(I+1,J,K)=MVFOUC(MANF-MAUF+I-1)
-                  ENDIF
-    6   CONTINUE
-3		CONTINUE
-		GGIND=GGIND+MAUF
-4	CONTINUE
-
-	RETURN
-	END SUBROUTINE PHGRACUT
-
-C Berechnung der Divergenz aus dem Windfeld (U,V)
-C	im Phasenraum. Zurueckgegeben werden die Felder der
-C	Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
-C	GWSAVE ist ein Hilfsfeld fuer die FFT
-C	P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
-C	MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
-C	MNAUF gibt die spektrale Aufloesung an, 
-C	NI = Anzahl der Gauss'schen Gitterpunkte,
-C	NJ = Anzahl der Gauss'schen Breiten,
-C	NK = Anzahl der Niveaus
-C Beachte, dass das Windfeld eine um 1 erhoehte Aufloesung in mu-Richtung hat.
- 
-	SUBROUTINE CONTGL(PS,DPSDL,DPSDM,DIV,U,V,BREITE,ETA,
-     *MLAT,A,B,NI,NJ,NK)
-				
-	IMPLICIT NONE
-	
-	INTEGER NI,NJ,NK,I,J,K,MLAT(NJ),L
-	
-	REAL A(NK+1),B(NK+1)
-	REAL PS(NI),DPSDL(NI),DPSDM(NI)
-	REAL DIV(NI,NK),U(NI,NK),V(NI,NK),ETA(NI,NK)
-	REAL BREITE(NJ)
-	
-	REAL DIVT1,DIVT2,POB,PUN,DPSDT,COSB
-	
-        L=0
-	DO 4 J=1,NJ
-		COSB=(1.0-BREITE(J)*BREITE(J))
-		DO 3 I=1,MLAT(J)
-                        L=L+1
-			DIVT1=0.0
-			DIVT2=0.0
-			DO 1 K=1,NK
-				POB=A(K)+B(K)*PS(L)
-				PUN=A(K+1)+B(K+1)*PS(L)
-				
-				DIVT1=DIVT1+DIV(L,K)*(PUN-POB)
-				if(cosb .gt. 0.) then
-                                  DIVT2=DIVT2+(B(K+1)-B(K))*PS(L)*
-     *(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
-                                endif
-     
-     		ETA(L,K)=-DIVT1-DIVT2
-1			CONTINUE
-
-			DPSDT=(-DIVT1-DIVT2)/PS(L)
-			
-			DO 2 K=1,NK
-				ETA(L,K)=ETA(L,K)-DPSDT*B(K+1)*PS(L)
-2			CONTINUE
-			PS(L)=DPSDT*PS(L)
-3		CONTINUE
-4	CONTINUE
-	RETURN
-	END SUBROUTINE CONTGL
-	
-C	OMEGA berechnet omega im Hybridkoordinatensystem
-C	PS ist der Bodendruck,
-C	DPSDL,DPSDM sind die Komponenten des Gradienten des Logarithmus des
-C	Bodendrucks
-C	DIV,U,V sind die horizontale Divergenz und das horizontale Windfeld
-C	BREITE ist das Feld der Gauss'schen Breiten
-C	E ist omega,
-
-	SUBROUTINE OMEGA(PS,DPSDL,DPSDM,DIV,U,V,BREITE,E,MLAT,A,B,NGI
-     *	 ,NGJ,MKK) 
-				
-	IMPLICIT NONE
-		
-	INTEGER I,J,K,L,NGI,NGJ,MKK,MLAT(NGJ)
-	
-	REAL PS(NGI),DPSDL(NGI),DPSDM(NGI),A(MKK+1),B(MKK+1)
-	REAL DIV(NGI,MKK),U(NGI,MKK),V(NGI,MKK),E(NGI,MKK)
-	REAL BREITE(NGJ)
-	
-	REAL DIVT1,DIVT2,POB,PUN,DP,X,Y,COSB
-	REAL DIVT3(MKK+2)
-	
-	L=0
-	DO 4 J=1,NGJ
-		COSB=(1.0-BREITE(J)*BREITE(J))
-		DO 3 I=1,MLAT(J)
-			L=L+1
-			DIVT1=0.0
-			DIVT2=0.0
-			DIVT3(1)=0.0
-			DO 1 K=1,MKK
-				POB=A(K)+B(K)*PS(L)
-				PUN=A(K+1)+B(K+1)*PS(L)
-				DP=PUN-POB
-				
-		Y=PS(L)*(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
-     		IF(K.LT.3) THEN
-     		X=0.0
-     		ELSE
-				X=(B(K+1)-B(K))*Y
-				ENDIF
-
-				DIVT1=DIVT1+DIV(L,K)*DP			
-				DIVT2=DIVT2+X
-     
-     		DIVT3(K+1)=-DIVT1-DIVT2
-     		
-     		IF(K.GT.1) THEN
-				E(L,K) = 0.5*(POB+PUN)/DP*Y*
-     *((B(K+1)-B(K))+(A(K+1)*B(K)-A(K)*B(K+1))/
-     *DP*LOG(PUN/POB))
-				ELSE
-					E(L,K) = 0.0
-				ENDIF
-								
-			E(L,K) = E(L,K)+0.5*(DIVT3(K)+DIVT3(K+1))
-
-1			CONTINUE
-3		CONTINUE
-4	CONTINUE
-	RETURN
-	END SUBROUTINE OMEGA
-	
-	END MODULE FTRAFO
Index: Source/Fortran/ftrafo.f90
===================================================================
--- Source/Fortran/ftrafo.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/ftrafo.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,504 @@
+MODULE FTRAFO
+
+!! Implementation of the spectral transformation using reduced the Gaussian grid
+
+CONTAINS
+
+! Implementierung der spektralen Transformationsmethode unter Verwendung
+! des reduzierten Gauss'schen Gitters
+
+
+  SUBROUTINE VDTOUV(XMN,XLAM,XPHI,GWSAVE,IFAX,P,MLAT,MNAUF,NI,NJ,NK)
+
+!! Berechnung der scale winds aus Vorticity und Divergenz
+!! uebergibt man in XMN die Divergenz, so wird der divergente Anteil des
+!! Windes (XPHI=Ud,XPHI=Vd) zurueckgegeben, uebergibt man die Vorticity, so
+!! erhaelt man den rotationellen Wind (XLAM=Vrot,XPHI=-Urot).
+!! Summiert man beide, erhaelt man den gesamten Scale wind
+! GWSAVE ist ein Hilfsfeld fuer die FFT
+! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
+! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
+! MNAUF gibt die spektrale Aufloesung an,
+! NI = Anzahl der Gauss'schen Gitterpunkte pro Flaeche
+! NJ = Anzahl der Gauss'schen Breiten,
+! NK = Anzahl der Niveaus
+
+    USE PHTOGR
+
+    IMPLICIT NONE
+    INTEGER J,N,NI,NJ,NK,MNAUF,GGIND(NJ/2)
+    INTEGER MLAT(NJ),IFAX(10,NJ)
+    REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
+    REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
+    REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
+    REAL XLAM(NI,NK),XPHI(NI,NK)
+    REAL GWSAVE(8*NJ+15,NJ/2)
+    REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
+    REAL RT,IT
+
+    GGIND(1)=0
+    DO 4 J = 2,NJ/2
+      GGIND(J)=GGIND(J-1)+MLAT(J-1)
+4   CONTINUE
+
+!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
+    DO 5 J = 1,NJ/2
+      CALL VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,GGIND(J),MLAT,MNAUF,NI,NJ,NK)
+5   CONTINUE
+!$OMP END PARALLEL DO
+
+    RETURN
+
+  END SUBROUTINE VDTOUV
+
+  SUBROUTINE VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,GGIND,MLAT,MNAUF,NI,NJ,NK)
+
+    USE PHTOGR
+
+    IMPLICIT NONE
+
+    INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
+    INTEGER MLAT(NJ),IFAX(10,NJ)
+    REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
+    REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
+    REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
+    REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
+    REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
+    REAL XLAM(NI,NK),XPHI(NI,NK)
+    REAL GWSAVE(8*NJ+15,NJ/2)
+    REAL ERAD,SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
+    REAL FAC(0:MNAUF),RT,IT
+
+    ERAD = 6367470.D0
+
+    FAC(0)=0.D0
+    DO 12 N=1,MNAUF
+      FAC(N)=-ERAD/DBLE(N)/DBLE(N+1)
+12  CONTINUE
+
+    CALL DPLGND(MNAUF,P(0,J),H)
+    DO 3 K = 1,NK
+      LL=0
+      LLP=0
+      LLH=0
+      DO 2 M = 0,MNAUF
+        SCR=0.D0
+        SCI=0.D0
+        ACR=0.D0
+        ACI=0.D0
+        MUSCR=0.D0
+        MUSCI=0.D0
+        MUACR=0.D0
+        MUACI=0.D0
+        LLS=LL
+        LLPS=LLP
+        LLHS=LLH
+        IF (2*M+1 .LT. MLAT(J)) THEN
+          DO 1 N = M,MNAUF,2
+            RT=XMN(2*LL,K)*FAC(N)
+            IT=XMN(2*LL+1,K)*FAC(N)
+            SCR =SCR+ RT*P(LLP,J)
+            SCI =SCI+ IT*P(LLP,J)
+            MUACR =MUACR+ RT*H(LLH)
+            MUACI =MUACI+ IT*H(LLH)
+            LL=LL+2
+            LLP=LLP+2
+            LLH=LLH+2
+1         CONTINUE
+          LL=LLS+1
+          LLP=LLPS+1
+          LLH=LLHS+1
+          DO 11 N = M+1,MNAUF,2
+            RT=XMN(2*LL,K)*FAC(N)
+            IT=XMN(2*LL+1,K)*FAC(N)
+            ACR =ACR+ RT*P(LLP,J)
+            ACI =ACI+ IT*P(LLP,J)
+            MUSCR =MUSCR+ RT*H(LLH)
+            MUSCI =MUSCI+ IT*H(LLH)
+            LL=LL+2
+            LLP=LLP+2
+            LLH=LLH+2
+11        CONTINUE
+        END IF
+        LL=LLS+(MNAUF-M+1)
+        LLP=LLPS+(MNAUF-M+3)
+        LLH=LLHS+(MNAUF-M+2)
+
+        UFOUC(2*M)=-M*(SCI-ACI)
+        UFOUC(2*M+1)=M*(SCR-ACR)
+        VFOUC(2*M)=-M*(SCI+ACI)
+        VFOUC(2*M+1)=M*(SCR+ACR)
+
+        MUFOUC(2*M)=-(MUSCR-MUACR)
+        MUFOUC(2*M+1)=-(MUSCI-MUACI)
+        MVFOUC(2*M)=-(MUSCR+MUACR)
+        MVFOUC(2*M+1)=-(MUSCI+MUACI)
+2     CONTINUE
+
+      CALL RFOURTR(VFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
+      CALL RFOURTR(UFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
+
+      CALL RFOURTR(MVFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
+      CALL RFOURTR(MUFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
+
+3   CONTINUE
+
+    RETURN
+
+  END SUBROUTINE VDUVSUB
+
+
+  SUBROUTINE PHGRAD(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MLAT,MNAUF,NI,NJ,NK)
+
+!! Berechnung des Gradienten eines Skalars aus dem Feld des
+!! Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
+!! Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
+! GWSAVE ist ein Hilfsfeld fuer die FFT
+! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
+! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
+! MNAUF gibt die spektrale Aufloesung an,
+! NI = Anzahl der Gauss'schen Gitterpunkte,
+! NJ = Anzahl der Gauss'schen Breiten,
+! NK = Anzahl der Niveaus
+
+    USE PHTOGR
+
+    IMPLICIT NONE
+
+    INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
+    INTEGER MLAT(NJ),IFAX(10,NJ)
+    REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
+    REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
+    REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
+    REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
+    REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
+    REAL XLAM(NI,NK),XPHI(NI,NK)
+    REAL GWSAVE(8*NJ+15,NJ/2)
+    REAL ERAD
+    REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
+
+    ERAD = 6367470.0
+
+    GGIND=0
+    DO 4 J = 1,NJ/2
+      CALL DPLGND(MNAUF,P(0,J),H)
+      DO 3 K = 1,NK
+        LL=0
+        LLP=0
+        LLH=0
+        DO 2 M = 0,MNAUF
+          SCR=0.D0
+          SCI=0.D0
+          ACR=0.D0
+          ACI=0.D0
+          MUSCR=0.D0
+          MUSCI=0.D0
+          MUACR=0.D0
+          MUACI=0.D0
+          LLS=LL
+          LLPS=LLP
+          LLHS=LLH
+          IF (2*M+1 .LT. MLAT(J)) THEN
+            DO 1 N = M,MNAUF,2
+              RT=XMN(2*LL,K)
+              IT=XMN(2*LL+1,K)
+              SCR =SCR+ RT*P(LLP,J)
+              SCI =SCI+ IT*P(LLP,J)
+              MUACR =MUACR+RT*H(LLH)
+              MUACI =MUACI+ IT*H(LLH)
+              LL=LL+2
+              LLP=LLP+2
+              LLH=LLH+2
+1           CONTINUE
+            LL=LLS+1
+            LLP=LLPS+1
+            LLH=LLHS+1
+            DO 11 N = M+1,MNAUF,2
+              RT=XMN(2*LL,K)
+              IT=XMN(2*LL+1,K)
+              ACR =ACR+ RT*P(LLP,J)
+              ACI =ACI+ IT*P(LLP,J)
+              MUSCR =MUSCR+ RT*H(LLH)
+              MUSCI =MUSCI+ IT*H(LLH)
+              LL=LL+2
+              LLP=LLP+2
+              LLH=LLH+2
+11          CONTINUE
+          END IF
+          LL=LLS+(MNAUF-M+1)
+          LLP=LLPS+(MNAUF-M+3)
+          LLH=LLHS+(MNAUF-M+2)
+
+          UFOUC(2*M)=-M*(SCI-ACI)/ERAD
+          UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
+          VFOUC(2*M)=-M*(SCI+ACI)/ERAD
+          VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
+
+          MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
+          MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
+          MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
+          MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
+2       CONTINUE
+
+        CALL RFOURTR(VFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+        XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
+        CALL RFOURTR(UFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+        XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
+
+        CALL RFOURTR(MVFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+        XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
+        CALL RFOURTR(MUFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
+        XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
+
+3     CONTINUE
+      GGIND=GGIND+MLAT(J)
+4   CONTINUE
+
+
+    RETURN
+
+  END SUBROUTINE PHGRAD
+
+
+  SUBROUTINE PHGRACUT(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MAUF,MNAUF,NI,NJ,MANF,NK)
+
+!! Berechnung des Gradienten eines Skalars aus dem Feld des
+!! Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
+!! Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
+! GWSAVE ist ein Hilfsfeld fuer die FFT
+! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
+! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
+! MNAUF gibt die spektrale Aufloesung an,
+! NI = Anzahl der Gauss'schen Gitterpunkte,
+! NJ = Anzahl der Gauss'schen Breiten,
+! NK = Anzahl der Niveaus
+
+    USE PHTOGR
+
+    IMPLICIT NONE
+
+    INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
+    INTEGER MAUF,MANF,I,IFAX(10)
+    REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
+    REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
+    REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
+    REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ)
+    REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
+    REAL XLAM(NI,NJ,NK),XPHI(NI,NJ,NK)
+    REAL HLAM(MAXAUF,2),HPHI(MAXAUF,2)
+    REAL GWSAVE(4*MAUF+15)
+    REAL ERAD
+    REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
+
+    ERAD = 6367470.0
+
+    GGIND=0
+    DO 4 J = 1,NJ
+      CALL DPLGND(MNAUF,P(0,J),H)
+      DO 3 K = 1,NK
+        LL=0
+        LLP=0
+        LLH=0
+        DO 2 M = 0,MNAUF
+          SCR=0.D0
+          SCI=0.D0
+          ACR=0.D0
+          ACI=0.D0
+          MUSCR=0.D0
+          MUSCI=0.D0
+          MUACR=0.D0
+          MUACI=0.D0
+          LLS=LL
+          LLPS=LLP
+          LLHS=LLH
+          IF (2*M+1 .LT. MAUF) THEN
+            DO 1 N = M,MNAUF,2
+              RT=XMN(2*LL,K)
+              IT=XMN(2*LL+1,K)
+              SCR =SCR+ RT*P(LLP,J)
+              SCI =SCI+ IT*P(LLP,J)
+              MUACR =MUACR+RT*H(LLH)
+              MUACI =MUACI+ IT*H(LLH)
+              LL=LL+2
+              LLP=LLP+2
+              LLH=LLH+2
+1           CONTINUE
+            LL=LLS+1
+            LLP=LLPS+1
+            LLH=LLHS+1
+            DO 11 N = M+1,MNAUF,2
+              RT=XMN(2*LL,K)
+              IT=XMN(2*LL+1,K)
+              ACR =ACR+ RT*P(LLP,J)
+              ACI =ACI+ IT*P(LLP,J)
+              MUSCR =MUSCR+ RT*H(LLH)
+              MUSCI =MUSCI+ IT*H(LLH)
+              LL=LL+2
+              LLP=LLP+2
+              LLH=LLH+2
+11          CONTINUE
+          END IF
+          LL=LLS+(MNAUF-M+1)
+          LLP=LLPS+(MNAUF-M+3)
+          LLH=LLHS+(MNAUF-M+2)
+
+          UFOUC(2*M)=-M*(SCI-ACI)/ERAD
+          UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
+          VFOUC(2*M)=-M*(SCI+ACI)/ERAD
+          VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
+
+          MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
+          MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
+          MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
+          MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
+2       CONTINUE
+
+        CALL RFOURTR(VFOUC,GWSAVE,IFAX,MNAUF,MAUF,1)
+
+        CALL RFOURTR(MVFOUC,GWSAVE,IFAX,MNAUF,MAUF,1)
+
+        DO 6 I=0,NI-1
+          IF (MANF+I .LE.  MAUF) THEN
+            XLAM(I+1,J,K)=VFOUC(MANF+I-1)
+            XPHI(I+1,J,K)=MVFOUC(MANF+I-1)
+          ELSE
+            XLAM(I+1,J,K)=VFOUC(MANF-MAUF+I-1)
+            XPHI(I+1,J,K)=MVFOUC(MANF-MAUF+I-1)
+          END IF
+6       CONTINUE
+3     CONTINUE
+      GGIND=GGIND+MAUF
+4   CONTINUE
+
+    RETURN
+
+  END SUBROUTINE PHGRACUT
+
+  SUBROUTINE CONTGL(PS,DPSDL,DPSDM,DIV,U,V,BREITE,ETA,MLAT,A,B,NI,NJ,NK)
+
+!! Berechnung der Divergenz aus dem Windfeld (U,V)
+!! im Phasenraum. Zurueckgegeben werden die Felder der
+!! Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
+! GWSAVE ist ein Hilfsfeld fuer die FFT
+! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
+! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
+! MNAUF gibt die spektrale Aufloesung an,
+! NI = Anzahl der Gauss'schen Gitterpunkte,
+! NJ = Anzahl der Gauss'schen Breiten,
+! NK = Anzahl der Niveaus
+! Beachte, dass das Windfeld eine um 1 erhoehte Aufloesung in mu-Richtung hat.
+
+    IMPLICIT NONE
+
+    INTEGER NI,NJ,NK,I,J,K,MLAT(NJ),L
+
+    REAL A(NK+1),B(NK+1)
+    REAL PS(NI),DPSDL(NI),DPSDM(NI)
+    REAL DIV(NI,NK),U(NI,NK),V(NI,NK),ETA(NI,NK)
+    REAL BREITE(NJ)
+
+    REAL DIVT1,DIVT2,POB,PUN,DPSDT,COSB
+
+    L=0
+    DO 4 J=1,NJ
+      COSB=(1.0-BREITE(J)*BREITE(J))
+      DO 3 I=1,MLAT(J)
+        L=L+1
+        DIVT1=0.0
+        DIVT2=0.0
+        DO 1 K=1,NK
+          POB=A(K)+B(K)*PS(L)
+          PUN=A(K+1)+B(K+1)*PS(L)
+
+          DIVT1=DIVT1+DIV(L,K)*(PUN-POB)
+          IF (COSB .GT. 0.) THEN
+            DIVT2=DIVT2+(B(K+1)-B(K))*PS(L)* &
+              (U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
+          END IF
+
+          ETA(L,K)=-DIVT1-DIVT2
+1       CONTINUE
+
+        DPSDT=(-DIVT1-DIVT2)/PS(L)
+
+        DO 2 K=1,NK
+          ETA(L,K)=ETA(L,K)-DPSDT*B(K+1)*PS(L)
+2       CONTINUE
+        PS(L)=DPSDT*PS(L)
+3     CONTINUE
+4   CONTINUE
+
+    RETURN
+
+  END SUBROUTINE CONTGL
+
+  SUBROUTINE OMEGA(PS,DPSDL,DPSDM,DIV,U,V,BREITE,E,MLAT,A,B,NGI,NGJ,MKK)
+
+!! calculates $\omega$ in the hybrid ($\eta$-) coordinate system
+
+! OMEGA berechnet omega im Hybridkoordinatensystem
+! PS ist der Bodendruck,
+! DPSDL,DPSDM sind die Komponenten des Gradienten des Logarithmus des
+! Bodendrucks
+! DIV,U,V sind die horizontale Divergenz und das horizontale Windfeld
+! BREITE ist das Feld der Gauss'schen Breiten
+! E ist omega,
+
+    IMPLICIT NONE
+
+    INTEGER I,J,K,L,NGI,NGJ,MKK,MLAT(NGJ)
+
+    REAL PS(NGI),DPSDL(NGI),DPSDM(NGI),A(MKK+1),B(MKK+1)
+    REAL DIV(NGI,MKK),U(NGI,MKK),V(NGI,MKK),E(NGI,MKK)
+    REAL BREITE(NGJ)
+
+    REAL DIVT1,DIVT2,POB,PUN,DP,X,Y,COSB
+    REAL DIVT3(MKK+2)
+
+    L=0
+    DO 4 J=1,NGJ
+      COSB=(1.0-BREITE(J)*BREITE(J))
+      DO 3 I=1,MLAT(J)
+        L=L+1
+        DIVT1=0.0
+        DIVT2=0.0
+        DIVT3(1)=0.0
+        DO 1 K=1,MKK
+          POB=A(K)+B(K)*PS(L)
+          PUN=A(K+1)+B(K+1)*PS(L)
+          DP=PUN-POB
+
+          Y=PS(L)*(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
+          IF (K .LT. 3) THEN
+            X=0.0
+          ELSE
+            X=(B(K+1)-B(K))*Y
+          END IF
+
+          DIVT1=DIVT1+DIV(L,K)*DP
+          DIVT2=DIVT2+X
+
+          DIVT3(K+1)=-DIVT1-DIVT2
+
+          IF (K .GT. 1) THEN
+            E(L,K) = 0.5*(POB+PUN)/ &
+              DP*Y*((B(K+1)-B(K))+(A(K+1)*B(K)-A(K)*B(K+1))/DP*LOG(PUN/POB))
+          ELSE
+            E(L,K) = 0.0
+          END IF
+
+          E(L,K) = E(L,K)+0.5*(DIVT3(K)+DIVT3(K+1))
+
+1       CONTINUE
+3     CONTINUE
+4   CONTINUE
+
+    RETURN
+
+  END SUBROUTINE OMEGA
+
+END MODULE FTRAFO
Index: urce/Fortran/grphreal.f
===================================================================
--- Source/Fortran/grphreal.f	(revision None)
+++ 	(revision )
@@ -1,188 +1,0 @@
-      MODULE GRTOPH
-
-      USE PHTOGR
-
-      CONTAINS
-C
-      SUBROUTINE GRPH213(CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,
-     *MNAUF,MAXL,MAXB,MLEVEL)
-
-C     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
-C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
-C     RAUM AUF KUGELKOORDINATEN DURCH
-C
-C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
-C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
-C			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
-C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
-C			WSAVE  = Working Array fuer Fouriertransformation
-C     Z 	 = LEGENDREFUNKTIONSWERTE
-C
-C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
-C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
-C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
-C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
-C
-      IMPLICIT REAL (A-H,O-Z)
-
-
-C			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
-C			Gauss'schen Gitters
-	INTEGER MLAT(MAXB),ISIZE,IFAX(10,MAXB)
-			
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-      REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
-      
-C      LOGICAL*1 USED(((216*217)/2+1)*160)
-
-      DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-      REAL FELD(MAXL,MLEVEL)
-      DIMENSION WSAVE(8*MAXB+15,MAXB/2)
-      REAL*8 W(MAXB) 
-      DIMENSION IND(MAXB)
-			
-			
-      	IND(1)=0
-      DO 6 J=2,MAXB/2
-	 IND(j)=IND(J-1)+MLAT(J-1)
- 6    CONTINUE
-!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
-      DO 16 L=1,MLEVEL
-        CALL GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,
-     *MNAUF,MAXL,MAXB,MLEVEL)
-16    CONTINUE
-!$omp end parallel do
-
-
-      RETURN
-      END SUBROUTINE GRPH213
-C
-      SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,
-     *MNAUF,MAXL,MAXB,MLEVEL)
-
-C     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
-C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
-C     RAUM AUF KUGELKOORDINATEN DURCH
-C
-C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
-C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
-C			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
-C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
-C			WSAVE  = Working Array fuer Fouriertransformation
-C     Z 	 = LEGENDREFUNKTIONSWERTE
-C
-C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
-C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
-C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
-C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
-C
-      IMPLICIT REAL (A-H,O-Z)
-
-C    FELD DER FOURIERKOEFFIZIENTEN
-      REAL CXMS(4*(MNAUF+1))
-      REAL CXMA(4*(MNAUF+1))
-      REAL,ALLOCATABLE :: CXM(:,:)
-
-C			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
-C			Gauss'schen Gitters
-			INTEGER MLAT(MAXB),ISIZE
-			
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-      REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
-      
-C      LOGICAL*1 USED(((216*217)/2+1)*160)
-
-      REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-      REAL FELD(MAXL,MLEVEL)
-      REAL WSAVE(8*MAXB+15,MAXB/2)
-      INTEGER IFAX(10,MAXB)
-      REAL W(MAXB)
-      INTEGER IND(MAXB)
-			
-      ALLOCATE(CXM( 4*MAXB,MAXB))
-	  DO 5 J=1,MAXB/2
-            CXMS(1:MLAT(J))=FELD(IND(J)+1:IND(J)+MLAT(J),L)
-	    CALL RFOUFTR(CXMS,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
-            CXMA(1:MLAT(J))=FELD(MAXL-IND(J)-MLAT(J)+1:MAXL-IND(J),L)
-	    CALL RFOUFTR(CXMA,
-     *WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
-	   	DO 4 I=1,2*(MNAUF+1)
-	    	CXM(I,J)=CXMS(I)+CXMA(I)
-	    	CXM(I,MAXB+1-J)=CXMS(I)-CXMA(I)
-4			CONTINUE
-    5 	CONTINUE
-        CALL LGTR213(CXMN(0,L),CXM,Z,W,MLAT,MNAUF,MAXB)
-
-      DEALLOCATE(CXM)
-
-      RETURN
-      END SUBROUTINE GRPHSUB
-C
-        SUBROUTINE LGTR213(CXMN,CXM,Z,W,MLAT,MNAUF,MAXB)
-        IMPLICIT REAL (A-H,O-Z)
-        INTEGER MLAT(MAXB)
-        DIMENSION CXM(0:4*MAXB-1,MAXB)
-        DIMENSION CXMN(0:2*(((MNAUF+1)*MNAUF)/2+MNAUF)+1)
-        REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
-        REAL*8 W(MAXB),CR,CI,HILF
-        LOGICAL EVEN
-C
-C     DIESE ROUTINE BERECHNET DIE KFFKs CXMN
-C
-			LL=0
-			LLP=0
-      DO 1 I=0,MNAUF
-        KM=0
-    9   KM=KM+1
-        IF(MLAT(KM).LE.2*I) THEN
-           GOTO 9
-        ENDIF
-        DO 2 J=I,MNAUF
-          CR=0
-          CI=0
-          EVEN=MOD(I+J,2).EQ.0
-          IF(EVEN) THEN
-            DO 3 K=KM,MAXB/2
-        			  HILF=W(K)*Z(K,LLP)
-        	      CR=CR+CXM(2*I,K)*HILF
-        	      CI=CI+CXM(2*I+1,K)*HILF
-    3     	CONTINUE
-          ELSE
-            DO 4 K=KM,MAXB/2
-        			  HILF=W(K)*Z(K,LLP)
-        	      CR=CR+CXM(2*I,MAXB+1-K)*HILF
-        	      CI=CI+CXM(2*I+1,MAXB+1-K)*HILF
-    4     	CONTINUE
-          ENDIF
-    5     CXMN(2*LL)=CR
-          CXMN(2*LL+1)=CI
-          LL=LL+1
-          LLP=LLP+1
-    2   CONTINUE
-        LLP=LLP+2
-    1 CONTINUE
-        RETURN
-        END SUBROUTINE LGTR213
-C
-
-C
-      SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
-C     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION CXM(0:2*MAXL-1)
-      DIMENSION FELD(MAXL),TRIGS(2*MAXL)
-      DIMENSION WSAVE(MAXAUF)
-      INTEGER IFAX(10)
-
-
-C NORMIERUNG...
-      WSAVE(1)=CXM(MAXL-1)
-
-      CXM(1:MAXL)=CXM(0:MAXL-1)/2
-      CXM(0)=WSAVE(1)/2
-!      CALL CFFTF(MAXL,CXM,WSAVE)
-      CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,-1)
-      RETURN
-      END SUBROUTINE RFOUFTR
-
-      END MODULE GRTOPH
Index: Source/Fortran/grphreal.f90
===================================================================
--- Source/Fortran/grphreal.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/grphreal.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,185 @@
+MODULE GRTOPH
+
+  USE PHTOGR
+
+CONTAINS
+
+  SUBROUTINE GRPH213(CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+
+!! WRONG>>> DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
+!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+!! RAUM AUF KUGELKOORDINATEN DURCH
+
+! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+! CXM   = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
+! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
+! WSAVE  = Working Array fuer Fouriertransformation
+! Z   = LEGENDREFUNKTIONSWERTE
+!
+! MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
+! MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+! MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+! MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+
+    IMPLICIT REAL (A-H,O-Z)
+
+!   Anzahl der Gitterpunkte pro Breitenkreis des reduzierten
+!   Gauss'schen Gitters
+    INTEGER MLAT(MAXB),ISIZE,IFAX(10,MAXB)
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
+
+!   LOGICAL*1 USED(((216*217)/2+1)*160)
+
+    DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    DIMENSION WSAVE(8*MAXB+15,MAXB/2)
+    REAL W(MAXB)
+    DIMENSION IND(MAXB)
+
+    IND(1)=0
+    DO 6 J=2,MAXB/2
+      IND(j)=IND(J-1)+MLAT(J-1)
+6   CONTINUE
+!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
+    DO 16 L=1,MLEVEL
+      CALL GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+16  CONTINUE
+!$omp end parallel do
+
+    RETURN
+
+  END SUBROUTINE GRPH213
+
+  SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+
+!! DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
+!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+!! RAUM AUF KUGELKOORDINATEN DURCH
+!
+! CXMN  = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+!         CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+! CXM   = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
+! FELD  = FELD DER METEOROLOGISCHEN VARIABLEN
+! WSAVE = Working Array fuer Fouriertransformation
+! Z     = LEGENDREFUNKTIONSWERTE
+!
+! MNAUF ANZAHL DER FOURIERKOEFFIZIENTEN
+! MAXL  ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+! MAXB  ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+
+    IMPLICIT REAL (A-H,O-Z)
+
+!   FELD DER FOURIERKOEFFIZIENTEN
+    REAL CXMS(4*(MNAUF+1))
+    REAL CXMA(4*(MNAUF+1))
+    REAL,ALLOCATABLE :: CXM(:,:)
+
+!   Anzahl der Gitterpunkte pro Breitenkreis des reduzierten
+!   Gauss'schen Gitters
+    INTEGER MLAT(MAXB),ISIZE
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
+
+!   LOGICAL*1 USED(((216*217)/2+1)*160)
+
+    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    REAL WSAVE(8*MAXB+15,MAXB/2)
+    INTEGER IFAX(10,MAXB)
+    REAL W(MAXB)
+    INTEGER IND(MAXB)
+
+    ALLOCATE(CXM( 4*MAXB,MAXB))
+    DO 5 J=1,MAXB/2
+      CXMS(1:MLAT(J))=FELD(IND(J)+1:IND(J)+MLAT(J),L)
+      CALL RFOUFTR(CXMS,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      CXMA(1:MLAT(J))=FELD(MAXL-IND(J)-MLAT(J)+1:MAXL-IND(J),L)
+      CALL RFOUFTR(CXMA,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
+      DO 4 I=1,2*(MNAUF+1)
+        CXM(I,J)=CXMS(I)+CXMA(I)
+        CXM(I,MAXB+1-J)=CXMS(I)-CXMA(I)
+4     CONTINUE
+5   CONTINUE
+    CALL LGTR213(CXMN(0,L),CXM,Z,W,MLAT,MNAUF,MAXB)
+
+    DEALLOCATE(CXM)
+
+    RETURN
+    
+  END SUBROUTINE GRPHSUB
+!
+  SUBROUTINE LGTR213(CXMN,CXM,Z,W,MLAT,MNAUF,MAXB)
+
+!!     DIESE ROUTINE BERECHNET DIE KFFKs CXMN
+
+    IMPLICIT REAL (A-H,O-Z)
+    INTEGER MLAT(MAXB)
+    DIMENSION CXM(0:4*MAXB-1,MAXB)
+    DIMENSION CXMN(0:2*(((MNAUF+1)*MNAUF)/2+MNAUF)+1)
+    REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
+    REAL*8 W(MAXB),CR,CI,HILF
+    LOGICAL EVEN
+
+
+    LL=0
+    LLP=0
+    DO 1 I=0,MNAUF
+      KM=0
+9     KM=KM+1
+      IF (MLAT(KM) .LE. 2*I) THEN
+        GOTO 9
+      END IF
+      DO 2 J=I,MNAUF
+        CR=0
+        CI=0
+        EVEN=MOD(I+J,2) .EQ. 0
+        IF (EVEN) THEN
+          DO 3 K=KM,MAXB/2
+            HILF=W(K)*Z(K,LLP)
+            CR=CR+CXM(2*I,K)*HILF
+            CI=CI+CXM(2*I+1,K)*HILF
+3         CONTINUE
+        ELSE
+          DO 4 K=KM,MAXB/2
+            HILF=W(K)*Z(K,LLP)
+            CR=CR+CXM(2*I,MAXB+1-K)*HILF
+            CI=CI+CXM(2*I+1,MAXB+1-K)*HILF
+4         CONTINUE
+        END IF
+5       CXMN(2*LL)=CR
+        CXMN(2*LL+1)=CI
+        LL=LL+1
+        LLP=LLP+1
+2     CONTINUE
+      LLP=LLP+2
+1   CONTINUE
+    RETURN
+    
+  END SUBROUTINE LGTR213
+
+  SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
+!
+! BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
+
+    IMPLICIT REAL (A-H,O-Z)
+    DIMENSION CXM(0:2*MAXL-1)
+    DIMENSION FELD(MAXL),TRIGS(2*MAXL)
+    DIMENSION WSAVE(MAXAUF)
+    INTEGER IFAX(10)
+
+! NORMIERUNG...
+    WSAVE(1)=CXM(MAXL-1)
+
+    CXM(1:MAXL)=CXM(0:MAXL-1)/2
+    CXM(0)=WSAVE(1)/2
+!    CALL CFFTF(MAXL,CXM,WSAVE)
+    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,-1)
+    RETURN
+  END SUBROUTINE RFOUFTR
+
+END MODULE GRTOPH
Index: urce/Fortran/makefile.CRAY
===================================================================
--- Source/Fortran/makefile.CRAY	(revision None)
+++ 	(revision )
@@ -1,57 +1,0 @@
-###############################################################################
-#
-# Top level Makefile for ECMWFDATA7.0 software
-#
-# Last modified:  December 1, 2015
-#
-###############################################################################
-
-
-.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
-           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
-
-OPT	= 
-DEBUG	= -g
-LIB     =  $(GRIB_API_LIBS) $(EMOSLIB)
-
-FC=ftn  $(F90FLAGS)
-F90C=ftn    $(F90FLAGS)
-
-FFLAGS =  $(OPT) -I. -r8 -I$(GRIB_API_INCLUDE_DIR)
-F90FLAGS =  $(OPT) -I. -r8 -I$(GRIB_API_INCLUDE_DIR)
-
-LDFLAGS =  $(OPT)
-
-BINDIR  =      .
-
-EXE	=      CONVERT2
-
-
-.f.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f
-.f90.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f90
-
-all:	${EXE}
-
-clean:
-	rm *.o *.mod ${EXE}
-
-phgrreal.o: phgrreal.f
-	$(F90C) -c -g -O3 phgrreal.f
-
-grphreal.o: grphreal.f
-	$(F90C) -c -g -O3 grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90C) -c -g -O3 ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
Index: urce/Fortran/makefile.ecgate
===================================================================
--- Source/Fortran/makefile.ecgate	(revision None)
+++ 	(revision )
@@ -1,58 +1,0 @@
-###############################################################################
-#
-# Top level Makefile for ECMWFDATA7.0 software
-#
-# Last modified:  December 1, 2015
-#
-###############################################################################
-
-
-.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
-           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
-
-OPT	= -O3
-DEBUG	= -g
-LIB     =  $(ECCODES_LIB) $(EMOSLIB)   
-
-FC=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian  
-F90C=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian 
-
-FFLAGS =  $(OPT) -I. -I$(ECCODES_INCLUDE_DIR) 
-
-F90FLAGS =  $(OPT) -I. -I$(ECCODES_INCLUDE_DIR)
-
-LDFLAGS =  $(OPT)
-
-BINDIR  =      .
-
-EXE	=      CONVERT2
-
-
-.f.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f
-.f90.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f90
-
-all:	${EXE}
-
-clean:
-	rm *.o *.mod ${EXE}
-
-phgrreal.o: phgrreal.f
-	$(F90C) -c -g -O3 -fopenmp phgrreal.f
-
-grphreal.o: grphreal.f
-	$(F90C) -c -g -O3 -fopenmp grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90C) -c -g -O3 -fopenmp ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
Index: urce/Fortran/makefile.local.gfortran
===================================================================
--- Source/Fortran/makefile.local.gfortran	(revision None)
+++ 	(revision )
@@ -1,63 +1,0 @@
-###############################################################################
-#
-# Top level Makefile for ECMWFDATA7.0 software
-#
-# Last modified:  December 1, 2015
-#
-###############################################################################
-
-
-.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
-           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
-
-
-GRIB_API_INCLUDE_DIR=/usr/local/gcc-4.9.3/grib_api-1.14.3/include
-GRIB_API_LIB= -L/usr/local/gcc-4.9.3/grib_api-1.14.3/lib -Bstatic  -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic  -lm  -ljasper
-
-EMOSLIB=-lemosR64
-
-OPT	= -O3 -fopenmp
-DEBUG	= -O3
-LIB     =  $(GRIB_API_LIB) $(EMOSLIB)
-
-FC=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian  
-F90C=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian 
-
-FFLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
-F90FLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
-
-LDFLAGS =  $(OPT)
-
-BINDIR  =      .
-
-EXE	=      CONVERT2
-
-
-.f.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f
-.f90.o:
-	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f90
-
-all:	${EXE}
-
-clean:
-	rm *.o *.mod ${EXE}
-
-phgrreal.o: phgrreal.f
-	$(F90C) -c -g -O3 -fopenmp phgrreal.f
-
-grphreal.o: grphreal.f
-	$(F90C) -c -g -O3 -fopenmp grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90C) -c -g -O3 -fopenmp ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
Index: urce/Fortran/makefile.local.ifort
===================================================================
--- Source/Fortran/makefile.local.ifort	(revision None)
+++ 	(revision )
@@ -1,61 +1,0 @@
-###############################################################################
-#
-# Top level Makefile for ECMWFDATA7.0 software
-#
-# Last modified:  December 1, 2015
-#
-###############################################################################
-
-
-.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
-           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
-
-GRIB_API_INCLUDE_DIR=/usr/local/include
-GRIB_API_LIB= -L/usr/local/lib -leccodes_f90  -leccodes -lm  -ljasper
-
-
-OPT	= -g -O3 -mcmodel=medium -unroll -inline -heap-arrays 32 
-DEBUG	= ''
-LIB     =  $(GRIB_API_LIB) -lemosR64 -lgfortran
-
-FC=/opt/intel/composerxe/bin/ifort -132  -traceback -r8
-F90C=/opt/intel/composerxe/bin/ifort -132  -traceback -r8
-
-FFLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
-F90FLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
-
-LDFLAGS =  $(OPT)
-
-BINDIR  =      .
-
-EXE	=      CONVERT2.s8.ifort
-
-
-.f.o:
-	$(F90C) -c $(F90FLAGS) -132 $(DEBUG) $*.f
-.f90.o:
-	$(F90C) -c $(F90FLAGS) -132 $(DEBUG) $*.f90
-
-all:	${EXE}
-
-clean:
-	rm *.o *.mod ${EXE}
-
-phgrreal.o: phgrreal.f
-	$(F90C) -c -g -O3 phgrreal.f
-
-grphreal.o: grphreal.f
-	$(F90C) -c -g -O3 grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90C) -c -g -O3 ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
Index: Source/Fortran/makefile_cray
===================================================================
--- Source/Fortran/makefile_cray	(revision None)
+++ Source/Fortran/makefile_cray	(revision None)
@@ -0,0 +1,57 @@
+###############################################################################
+#
+# Top level Makefile for ECMWFDATA7.0 software
+#
+# Last modified:  December 1, 2015
+#
+###############################################################################
+
+
+.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
+           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
+
+OPT	= 
+DEBUG	= -g
+LIB     =  $(GRIB_API_LIBS) $(EMOSLIB)
+
+FC=ftn  $(F90FLAGS)
+F90C=ftn    $(F90FLAGS)
+
+FFLAGS =  $(OPT) -I. -r8 -I$(GRIB_API_INCLUDE_DIR)
+F90FLAGS =  $(OPT) -I. -r8 -I$(GRIB_API_INCLUDE_DIR)
+
+LDFLAGS =  $(OPT)
+
+BINDIR  =      .
+
+EXE	=      calc_etadot
+
+
+.f.o:
+	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f
+.f90.o:
+	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f90
+
+all:	${EXE}
+
+clean:
+	rm *.o *.mod ${EXE}
+
+phgrreal.o: phgrreal.f
+	$(F90C) -c -g -O3 phgrreal.f
+
+grphreal.o: grphreal.f
+	$(F90C) -c -g -O3 grphreal.f
+
+ftrafo.o: ftrafo.f
+	$(F90C) -c -g -O3 ftrafo.f
+
+$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwgrib2.o  posnam.o calc_etadot.o
+	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwgrib2.o posnam.o calc_etadot.o ${LIB}
+
+
+###############################################################################
+#
+# End of the Makefile
+#
+###############################################################################
Index: Source/Fortran/makefile_debug
===================================================================
--- Source/Fortran/makefile_debug	(revision None)
+++ Source/Fortran/makefile_debug	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -2,9 +2,10 @@
 #
 # Makefile for flex_extract, Fortran code to calculate etadot
+# Makefile created using by mkmf 19.3.0
 #
 # Copyright: Leopold Haimberger, Petra Seibert
 # SPDX-License-Identifier: GPL-2.0 
 # 
-# Version for a machine with grib_api installed on standard paths
+# Version for a machine with grib_api and emoslib installed on standard paths
 # full debugging
 #
@@ -12,51 +13,40 @@
 
 
-GRIB_API_INCLUDE_DIR=/usr/include
-GRIB_API_LIB=  -Bstatic  -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic  -lm  -ljasper
+EXE      =  calc_etadot_debug.out
+
+GRIB_API_LIB=  -Bstatic -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper
 EMOSLIB=-lemosR64
-
 LIB =  $(GRIB_API_LIB) $(EMOSLIB)
 
-F90 = gfortran -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore -ffixed-line-length-132 -fopenmp -fconvert=big-endian
+GRIB_API_INCLUDE_DIR=/usr/include
+INC = -I. -I$(GRIB_API_INCLUDE_DIR)
 
-F90FLAGS =  -I. -I$(GRIB_API_INCLUDE_DIR)
+FC = gfortran 
+OPT = -g -Og -fbacktrace -fcheck=all
+FFLAGS =   $(OPT) $(LIB) $(INC) -fdefault-real-8 -fopenmp -fconvert=big-endian 
+LDFLAGS =  $(OPT) $(LIB) -fopenmp
+SRC = ./rwgrib2.f90 ./calc_etadot.f90 ./ftrafo.f90 ./grphreal.f90 ./posnam.f90 ./phgrreal.f90
+OBJ = rwgrib2.o calc_etadot.o ftrafo.o grphreal.o posnam.o phgrreal.o
+MOD = ftrafo.mod  grtoph.mod	phtogr.mod  rwgrib2.mod
 
-OPT = -g -fbacktrace -fbounds-check
+.DEFAULT:
+	-echo $@ does not exist.
 
-BINDIR   =  .
-EXE      =  calc_etadot_debug.out
+all: ${EXE}
+ftrafo.o: ./ftrafo.f90 phgrreal.o
+	$(FC)  $(FFLAGS)  -c	./ftrafo.f90
+grphreal.o: ./grphreal.f90 phgrreal.o
+	$(FC)  $(FFLAGS)  -c	./grphreal.f90
+phgrreal.o: ./phgrreal.f90
+	$(FC)  $(FFLAGS)  -c	./phgrreal.f90
+posnam.o: ./posnam.f90
+	$(FC)  $(FFLAGS)  -c	./posnam.f90
+calc_etadot.o: ./calc_etadot.f90 phgrreal.o grphreal.o ftrafo.o rwgrib2.o
+	$(FC)  $(FFLAGS)  -c	./calc_etadot.f90
+rwgrib2.o: ./rwgrib2.f90
+	$(FC)  $(FFLAGS)  -c	./rwgrib2.f90
 
-.f.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f
-.f90.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f90
+clean: -rm -f $(OBJ) ${EXE} $(MOD)
 
-all:	${EXE}
-
-clean:
-	rm *.o *.mod #${EXE}
-
-preconvert.o: preconvert.f90
-	$(F90) -c $(OPT) $(F90FLAGS) preconvert.f90
-        
-phgrreal.o: phgrreal.f
-	$(F90) -c $(OPT)   phgrreal.f
-
-rwGRIB2.o: rwGRIB2.f90
-	$(F90) -c $(OPT) $(F90FLAGS)  rwGRIB2.f90
-
-grphreal.o: grphreal.f
-	$(F90) -c $(OPT)   grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90) -c $(OPT)   ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
+${EXE}: $(OBJ)
+	$(FC) $(OBJ) -o ${EXE}  $(LDFLAGS)
Index: urce/Fortran/makefile_debug2
===================================================================
--- Source/Fortran/makefile_debug2	(revision 2d1c338814a949a88dd9f5b55ddec4efbd7b60e6)
+++ 	(revision )
@@ -1,66 +1,0 @@
-###############################################################################
-#
-# Makefile for flex_extract, Fortran code to calculate etadot
-#
-# Copyright: Leopold Haimberger, Petra Seibert
-# SPDX-License-Identifier: GPL-2.0 
-# 
-# Version for a machine with grib_api installed on standard paths
-# full debugging
-#
-###############################################################################
-
-
-GRIB_API_INCLUDE_DIR=/usr/include
-GRIB_API_LIB=  -Bstatic  -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic  -lm  -ljasper
-EMOSLIB=-lemosR64
-
-LIB =  $(GRIB_API_LIB) $(EMOSLIB)
-
-F90 = gfortran -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore -ffixed-line-length-132 -fopenmp -fconvert=big-endian
-
-F90FLAGS =  -I. -I$(GRIB_API_INCLUDE_DIR)
-
-OPT = -g -fbacktrace -fbounds-check
-
-BINDIR   =  .
-EXE      =  calc_etadot_debug.out
-
-.f.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f
-.f90.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f90
-
-all:	${EXE}
-
-clean:
-	rm *.o *.mod #${EXE}
-
-preconvert.o: preconvert.f90
-	$(F90) -c $(OPT) $(F90FLAGS) preconvert.f90
-        
-phgrreal.o: phgrreal.f
-	$(F90) -c $(OPT)   phgrreal.f
-
-rwGRIB2.o: rwGRIB2.f90
-	$(F90) -c $(OPT) $(F90FLAGS)  rwGRIB2.f90
-
-grphreal.o: grphreal.f
-	$(F90) -c $(OPT)   grphreal.f
-
-ftrafo.o: ftrafo.f
-	$(F90) -c $(OPT)   ftrafo.f
-rpassm.o: rpassm.f
-	$(F90) -c $(OPT)   rpassm.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o rpassm.o
-
-	$(F90) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o rpassm.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
-
Index: Source/Fortran/makefile_ecgate
===================================================================
--- Source/Fortran/makefile_ecgate	(revision None)
+++ Source/Fortran/makefile_ecgate	(revision None)
@@ -0,0 +1,58 @@
+###############################################################################
+#
+# Top level Makefile for ECMWFDATA7.0 software
+#
+# Last modified:  December 1, 2015
+#
+###############################################################################
+
+
+.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
+           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
+
+OPT	= -O3
+DEBUG	= -g
+LIB     =  $(ECCODES_LIB) $(EMOSLIB)   
+
+FC=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian  
+F90C=gfortran   -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore  -ffixed-line-length-132 -fopenmp  -fconvert=big-endian 
+
+FFLAGS =  $(OPT) -I. -I$(ECCODES_INCLUDE_DIR) 
+
+F90FLAGS =  $(OPT) -I. -I$(ECCODES_INCLUDE_DIR)
+
+LDFLAGS =  $(OPT)
+
+BINDIR  =      .
+
+EXE	=      calc_etadot
+
+
+.f.o:
+	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f
+.f90.o:
+	$(F90C) -c $(F90FLAGS) $(DEBUG) $*.f90
+
+all:	${EXE}
+
+clean:
+	rm *.o *.mod ${EXE}
+
+phgrreal.o: phgrreal.f
+	$(F90C) -c -g -O3 -fopenmp phgrreal.f
+
+grphreal.o: grphreal.f
+	$(F90C) -c -g -O3 -fopenmp grphreal.f
+
+ftrafo.o: ftrafo.f
+	$(F90C) -c -g -O3 -fopenmp ftrafo.f
+
+$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwgrib2.o  posnam.o calc_etadot.o
+	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwgrib2.o posnam.o calc_etadot.o ${LIB}
+
+
+###############################################################################
+#
+# End of the Makefile
+#
+###############################################################################
Index: Source/Fortran/makefile_fast
===================================================================
--- Source/Fortran/makefile_fast	(revision None)
+++ Source/Fortran/makefile_fast	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -2,61 +2,52 @@
 #
 # Makefile for flex_extract, Fortran code to calculate etadot
+# Makefile created using by mkmf 19.3.0
 #
 # Copyright: Leopold Haimberger, Petra Seibert
 # SPDX-License-Identifier: GPL-2.0 
 # 
-# Version for a machine with grib_api installed on standard paths
-# compiled for fast runs
+# Version for a machine with grib_api and emoslib installed on standard paths
+# full debugging
 #
 ###############################################################################
 
 
-GRIB_API_INCLUDE_DIR=/usr/include
-GRIB_API_LIB=  -Bstatic  -lgrib_api_f77 -lgrib_api_f90 -lgrib_api -Bdynamic  -lm  -ljasper
-EMOSLIB=-lemosR64
-
-LIB     =  $(GRIB_API_LIB) $(EMOSLIB)
-
-F90 = gfortran -m64 -fdefault-real-8 -fcray-pointer -fno-second-underscore -ffixed-line-length-132 -fopenmp -fconvert=big-endian
-
-F90FLAGS =  -I. -I$(GRIB_API_INCLUDE_DIR)
-
-OPT = -O3
-
-BINDIR   =  .
 EXE      =  calc_etadot_fast.out
 
-.f.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f
-.f90.o:
-	$(F90) -c $(F90FLAGS) $(OPT) $*.f90
+GRIB_API_LIB=  -Bstatic -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper
+EMOSLIB=-lemosR64
+LIB =  $(GRIB_API_LIB) $(EMOSLIB)
 
-all:	${EXE}
+GRIB_API_INCLUDE_DIR=/usr/include
+INC = -I. -I$(GRIB_API_INCLUDE_DIR)
 
-clean:
-	rm *.o *.mod #${EXE}
+FC = gfortran 
+OPT = -O3 -march=native
+FFLAGS =   $(OPT) $(LIB) $(INC) -fdefault-real-8 -fopenmp -fconvert=big-endian 
+LDFLAGS =  $(OPT) $(LIB) -fopenmp
+SRC = ./rwgrib2.f90 ./calc_etadot.f90 ./ftrafo.f90 ./grphreal.f90 ./posnam.f90 ./phgrreal.f90
+OBJ = rwgrib2.o calc_etadot.o ftrafo.o grphreal.o posnam.o phgrreal.o
+MOD = ftrafo.mod  grtoph.mod	phtogr.mod  rwgrib2.mod
 
-preconvert.o: preconvert.f90
-	$(F90) -c $(OPT) $(F90FLAGS) preconvert.f90
-        
-phgrreal.o: phgrreal.f
-	$(F90) -c $(OPT)   phgrreal.f
+.DEFAULT:
+	-echo $@ does not exist.
 
-rwGRIB2.o: rwGRIB2.f90
-	$(F90) -c $(OPT) $(F90FLAGS)  rwGRIB2.f90
+all: ${EXE}
+ftrafo.o: ./ftrafo.f90 phgrreal.o
+	$(FC)  $(FFLAGS)  -c	./ftrafo.f90
+grphreal.o: ./grphreal.f90 phgrreal.o
+	$(FC)  $(FFLAGS)  -c	./grphreal.f90
+phgrreal.o: ./phgrreal.f90
+	$(FC)  $(FFLAGS)  -c	./phgrreal.f90
+posnam.o: ./posnam.f90
+	$(FC)  $(FFLAGS)  -c	./posnam.f90
+calc_etadot.o: ./calc_etadot.f90 phgrreal.o grphreal.o ftrafo.o rwgrib2.o
+	$(FC)  $(FFLAGS)  -c	./calc_etadot.f90
+rwgrib2.o: ./rwgrib2.f90
+	$(FC)  $(FFLAGS)  -c	./rwgrib2.f90
 
-grphreal.o: grphreal.f
-	$(F90) -c $(OPT)   grphreal.f
+clean: 
+	-rm -f $(OBJ) ${EXE} $(MOD)
 
-ftrafo.o: ftrafo.f
-	$(F90) -c $(OPT)   ftrafo.f
-
-$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwGRIB2.o  posnam.o preconvert.o
-	$(F90) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwGRIB2.o posnam.o preconvert.o ${LIB}
-
-
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
+${EXE}: $(OBJ)
+	$(FC) $(OBJ) -o ${EXE}  $(LDFLAGS)
Index: Source/Fortran/makefile_local_gfortran
===================================================================
--- Source/Fortran/makefile_local_gfortran	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/makefile_local_gfortran	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,1 @@
+makefile_fast
Index: Source/Fortran/makefile_local_ifort
===================================================================
--- Source/Fortran/makefile_local_ifort	(revision None)
+++ Source/Fortran/makefile_local_ifort	(revision None)
@@ -0,0 +1,61 @@
+###############################################################################
+#
+# Top level Makefile for ECMWFDATA7.0 software
+#
+# Last modified:  December 1, 2015
+#
+###############################################################################
+
+
+.SUFFIXES: .o .c .c~ .f .f~ .F90 .f90 .f90~ .f95 .f95~ .F .F~ .y .y~ .l .l~ \
+           .s .s~ .sh .sh~ .h .h~ .C .C~ .a
+
+GRIB_API_INCLUDE_DIR=/usr/local/include
+GRIB_API_LIB= -L/usr/local/lib -leccodes_f90  -leccodes -lm  -ljasper
+
+
+OPT	= -g -O3 -mcmodel=medium -unroll -inline -heap-arrays 32 
+DEBUG	= ''
+LIB     =  $(GRIB_API_LIB) -lemosR64 -lgfortran
+
+FC=/opt/intel/composerxe/bin/ifort -132  -traceback -r8
+F90C=/opt/intel/composerxe/bin/ifort -132  -traceback -r8
+
+FFLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
+F90FLAGS =  $(OPT) -I. -I$(GRIB_API_INCLUDE_DIR)
+
+LDFLAGS =  $(OPT)
+
+BINDIR  =      .
+
+EXE	=      calc_etadot.s8.ifort
+
+
+.f.o:
+	$(F90C) -c $(F90FLAGS) -132 $(DEBUG) $*.f
+.f90.o:
+	$(F90C) -c $(F90FLAGS) -132 $(DEBUG) $*.f90
+
+all:	${EXE}
+
+clean:
+	rm *.o *.mod ${EXE}
+
+phgrreal.o: phgrreal.f
+	$(F90C) -c -g -O3 phgrreal.f
+
+grphreal.o: grphreal.f
+	$(F90C) -c -g -O3 grphreal.f
+
+ftrafo.o: ftrafo.f
+	$(F90C) -c -g -O3 ftrafo.f
+
+$(BINDIR)/${EXE}:	phgrreal.o grphreal.o ftrafo.o rwgrib2.o  posnam.o calc_etadot.o
+	$(F90C) $(DEBUG) $(OPT) -o $(BINDIR)/${EXE} ftrafo.o phgrreal.o grphreal.o rwgrib2.o posnam.o calc_etadot.o ${LIB}
+
+
+###############################################################################
+#
+# End of the Makefile
+#
+###############################################################################
Index: urce/Fortran/phgrreal.f
===================================================================
--- Source/Fortran/phgrreal.f	(revision None)
+++ 	(revision )
@@ -1,553 +1,0 @@
-      MODULE PHTOGR
-
-      INTEGER, PARAMETER :: MAXAUF=36000
-
-      CONTAINS
-
-      SUBROUTINE PHGR213(CXMN,FELD,WSAVE,IFAX,Z,MLAT,MNAUF,
-     *MAXL,MAXB,MLEVEL)
-
-C     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
-C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
-C     RAUM AUF DAS REDUZIERTE GAUSS'SCHE GITTER DURCH
-C
-C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
-C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
-C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
-C	WSAVE  = Working Array fuer Fouriertransformation
-C     Z 	 = LEGENDREFUNKTIONSWERTE
-C
-C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
-C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
-C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
-C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
-C
-      IMPLICIT NONE
-
-C			Anzahl der Gitterpunkte auf jedem Breitenkreis
-      INTEGER MLAT(MAXB/2)
-      INTEGER K,MAXL,MAXB,MLEVEL,MNAUF
-      INTEGER IND(MAXB)
-      
-			
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-      REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
-
-      REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-      REAL FELD(MAXL,MLEVEL)
-      REAL WSAVE(8*MAXB+15,MAXB/2)
-      INTEGER :: IFAX(10,MAXB)
-      
-      IND(1)=0
-     	DO 7 K=2,MAXB/2
-       	IND(K)=IND(K-1)+MLAT(K-1)
-7 		CONTINUE
-
-!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
-     	DO 17 K=1,MAXB/2
-        CALL PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,
-     *MNAUF,MAXL,MAXB,MLEVEL)
-
-17 		CONTINUE
-!$OMP END PARALLEL DO
-
-        RETURN
-        END SUBROUTINE PHGR213
-C
-C
-      SUBROUTINE PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,
-     *MNAUF,MAXL,MAXB,MLEVEL)
-
-      IMPLICIT NONE
-
-      INTEGER MLAT(MAXB/2)
-      INTEGER K,L,I,J,LLS,LLPS,LL,LLP,MAXL,MAXB,MLEVEL,MNAUF
-      INTEGER IND(MAXB)
-      INTEGER :: IFAX(10,MAXB)
-      
-			
-C    FELD DER FOURIERKOEFFIZIENTEN
-      REAL :: CXMS(0:MAXAUF-1),CXMA(0:MAXAUF-1)
-
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-      REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
-      REAL ACR,ACI,SCR,SCI
-
-      REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-      REAL FELD(MAXL,MLEVEL)
-      REAL WSAVE(8*MAXB+15,MAXB/2)
-
-      	DO 6 L=1,MLEVEL
-      	 LL=0
-      	 LLP=0
-         DO 1 I=0,MNAUF
-          SCR=0.D0
-          SCI=0.D0
-          ACR=0.D0
-          ACI=0.D0
-          LLS=LL
-          LLPS=LLP
-          IF(2*I+1.LT.MLAT(K)) THEN
-C	Innerste Schleife aufgespalten um if-Abfrage zu sparen
-            DO 18 J=I,MNAUF,2
-              SCR=SCR+Z(LLP,K)*CXMN(2*LL,L)
-              SCI=SCI+Z(LLP,K)*CXMN(2*LL+1,L)
-              LL=LL+2            
-              LLP=LLP+2            
-18          CONTINUE
-            LL=LLS+1
-            LLP=LLPS+1
-            DO 19 J=I+1,MNAUF,2
-              ACR=ACR+Z(LLP,K)*CXMN(2*LL,L)
-              ACI=ACI+Z(LLP,K)*CXMN(2*LL+1,L)
-              LL=LL+2            
-              LLP=LLP+2            
-19          CONTINUE
-          ENDIF
-          LL=LLS+(MNAUF-I+1)
-          LLP=LLPS+(MNAUF-I+3)
-          CXMS(2*I)=SCR+ACR
-          CXMS(2*I+1)=SCI+ACI
-          CXMA(2*I)=SCR-ACR
-          CXMA(2*I+1)=SCI-ACI
-    1    CONTINUE
-C         CALL FOURTR(CXMS,FELD(IND(k)+1,L),WSAVE(:,K),MNAUF,
-C     *MLAT(K),1)
-C         CALL FOURTR(CXMA,FELD(MAXL-IND(k)-MLAT(K)+1,L),
-C     *WSAVE(:,K),MNAUF,MLAT(K),1)
-         CALL RFOURTR(CXMS,WSAVE(:,K),IFAX(:,K),MNAUF,
-     *MLAT(K),1)
-          FELD(IND(k)+1:IND(K)+MLAT(K),L)=CXMS(0:MLAT(K)-1)
-         CALL RFOURTR(CXMA,
-     *WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K),1)
-         FELD(MAXL-IND(k)-MLAT(K)+1:MAXL-IND(k),L)=CXMA(0:MLAT(K)-1)
-C         WRITE(*,*) IND+1,FELD(IND+1,L)
-6 			CONTINUE
-
-      END SUBROUTINE PHSYM
-
-      SUBROUTINE PHGCUT(CXMN,FELD,WSAVE,IFAX,Z,
-     *                  MNAUF,MMAX,MAUF,MANF,MAXL,MAXB,MLEVEL)
-
-C     DIE ROUTINE FUEHRT EINE TRANSFORMATION EINER
-C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
-C     RAUM AUF KUGELKOORDINATEN DURCH. Es kann ein Teilausschnitt
-C			Der Erde angegeben werden. Diese Routine ist langsamer als
-C			phgrph
-C
-C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
-C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
-C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
-C     BREITE = SINUS DER GEOGRAFISCHEN BREITEN
-C
-C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
-C     MAUF     ANZAHL DER LAENGEN UND DER FOURIERKOEFFIZIENTEN
-C     MANF     ANFANG DES LAENGENBEREICHS FUER DAS GITTER,
-C              AUF DAS INTERPOLIERT WERDEN SOLL
-C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
-C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
-C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
-C
-        IMPLICIT REAL (A-H,O-Z)
-
-C    FELD DER FOURIERKOEFFIZIENTEN
-
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-        REAL Z(0:((MMAX+3)*(MMAX+4))/2,MAXB)
-
-        DIMENSION CXMN(0:(MMAX+1)*(MMAX+2)-1,MLEVEL)
-        REAL FELD(MAXL,MAXB,MLEVEL)
-      DIMENSION WSAVE(4*MAUF+15)
-      INTEGER:: IFAX(10)
-      
-      LOGICAL SYM
-
-C
-C      write(*,*)mauf,mnauf,manf,maxl
-
-     
-        IF(MAUF.LE.MNAUF) WRITE(*,*) 'TOO COARSE LONGITUDE RESOLUTION'
-          IF((MANF.LT.1).OR.(MAXL.LT.1).OR.
-     *       (MANF.GT.MAUF).OR.(MAXL.GT.MAUF)) THEN
-            WRITE(*,*) 'WRONG LONGITUDE RANGE',MANF,MAXL
-            STOP
-          ENDIF
-
-C Pruefe, ob Ausgabegitter symmetrisch zum Aequator ist
-C Wenn ja soll Symmetrie der Legendrepolynome ausgenutzt werden
-      IF(MAXB .GT. 4) THEN
-          SYM=.TRUE.
-	    DO 11 J=5,5
-	      IF(ABS(ABS(Z(100,J))-ABS(Z(100,MAXB+1-J))).GT.1E-11)
-     *    SYM=.FALSE.
-C	      WRITE(*,*) ABS(Z(100,J)),ABS(Z(100,MAXB+1-J))
-11	  CONTINUE
-      WRITE(*,*) 'Symmetrisch: ',SYM
-      ELSE
-        SYM=.FALSE.
-      ENDIF
-
-	
-		IF(SYM) THEN
-!$OMP PARALLEL DO 
-      	DO J=1,(MAXB+1)/2
-          CALL PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,
-     *MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
-
-        ENDDO
-!$OMP END PARALLEL DO
-		ELSE
-!$OMP PARALLEL DO 
-        DO J=1,MAXB
-            CALL PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,
-     *J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
-        ENDDO
-!$OMP END PARALLEL DO
-
-      ENDIF
-
-
-        RETURN
-        END SUBROUTINE PHGCUT
-
-        SUBROUTINE PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,
-     *MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
-
-        IMPLICIT REAL (A-H,O-Z)
-
-C    FELD DER FOURIERKOEFFIZIENTEN
-
-        REAL :: CXM(0:MAXAUF-1),CXMA(0:MAXAUF-1)
-
-
-C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
-        REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
-        REAL SCR,SCI,ACR,ACI
-
-        DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-        REAL FELD(MAXL,MAXB,MLEVEL)
-      DIMENSION WSAVE(4*MAUF+15)
-        INTEGER :: IFAX(10)
-
-        DO 16 L=1,MLEVEL
-      	 LL=0
-      	 LLP=0
-         DO 17 I=0,MNAUF
-          SCR=0.D0
-          SCI=0.D0
-          ACR=0.D0
-          ACI=0.D0
-          LLS=LL
-          LLPS=LLP
-C	Innerste Schleife aufgespalten um if-Abfrage zu sparen
-          DO 18 K=I,MNAUF,2
-            SCR=SCR+Z(LLP,J)*CXMN(2*LL,L)
-            SCI=SCI+Z(LLP,J)*CXMN(2*LL+1,L)
-            LL=LL+2            
-            LLP=LLP+2            
-18        CONTINUE
-          LL=LLS+1
-          LLP=LLPS+1
-          DO 19 K=I+1,MNAUF,2
-            ACR=ACR+Z(LLP,J)*CXMN(2*LL,L)
-            ACI=ACI +Z(LLP,J)*CXMN(2*LL+1,L)
-            LL=LL+2            
-            LLP=LLP+2            
-19        CONTINUE
-          LL=LLS+MNAUF-I+1
-          LLP=LLPS+MNAUF-I+3
-          CXM(2*I)=SCR+ACR
-          CXM(2*I+1)=SCI+ACI
-          CXMA(2*I)=SCR-ACR
-          CXMA(2*I+1)=SCI-ACI
-17       CONTINUE
-         
-         CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
-         DO 26 I=0,MAXL-1
-           IF(MANF+I.LE.MAUF) THEN
-             FELD(I+1,J,L)=CXM(MANF+I-1)
-           ELSE
-             FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
-           ENDIF
-26       CONTINUE
-         CALL RFOURTR(CXMA,WSAVE,IFAX,MNAUF,MAUF,1)
-         DO 36 I=0,MAXL-1
-           IF(MANF+I.LE.MAUF) THEN
-             FELD(I+1,MAXB+1-J,L)=CXMA(MANF+I-1)
-           ELSE
-             FELD(I+1,MAXB+1-J,L)=CXMA(MANF-MAUF+I-1)
-           ENDIF
-36       CONTINUE
-16 			CONTINUE
-
-      END SUBROUTINE PHSYMCUT
-
-      SUBROUTINE PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,
-     *J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
-
-        IMPLICIT NONE
-        INTEGER,intent(in) :: MNAUF,MAUF,MANF,J,MAXL,MAXB,MLEVEL
-        REAL :: CXM(0:MAXAUF-1)
-        REAL,intent(in) :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
-
-        REAL,intent(in) :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-
-        REAL,intent(in) :: WSAVE(4*MAUF+15)
-
-        REAL :: FELD(MAXL,MAXB,MLEVEL)
-        INTEGER :: IFAX(10)
-
-        INTEGER I,L
-
-          DO L=1,MLEVEL
-            CALL LEGTR(CXMN(:,L),CXM,Z(:,J),MNAUF,MAUF)
-            CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
-           
-            DO I=0,MAXL-1
-                  IF(MANF+I.LE.MAUF) THEN
-                    FELD(I+1,J,L)=CXM(MANF+I-1)
-                  ELSE
-                    FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
-                  ENDIF
-            ENDDO
-          ENDDO
-      END SUBROUTINE PHGPNS
-C
-        SUBROUTINE LEGTR(CXMN,CXM,Z,MNAUF,MAUF)
-        IMPLICIT NONE
-        INTEGER MNAUF,MAUF,LL,LLP,I,J
-        REAL CXM(0:MAXAUF-1)
-        REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1)
-        REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
-        REAL CI,CR
-C
-C     DIESE ROUTINE BERECHNET DIE FOURIERKOEFFIZIENTEN CXM
-C
-        LL=0
-        LLP=0
-        DO 1 I=0,MNAUF
-            CR=0.D0
-            CI=0.D0
-            DO 2 J=I,MNAUF
-              CR=CR+Z(LLP)*CXMN(2*LL)
-              CI=CI+Z(LLP)*CXMN(2*LL+1)
-              LL=LL+1
-              LLP=LLP+1
-    2     CONTINUE
-          LLP=LLP+2
-          CXM(2*I)=CR
-          CXM(2*I+1)=CI
-    1 CONTINUE
-        RETURN
-        END SUBROUTINE LEGTR
-C
-C     
-C
-      SUBROUTINE RFOURTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
-C     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
-      IMPLICIT REAL (A-H,O-Z)
-      DIMENSION CXM(0:MAXAUF-1)
-      REAL :: WSAVE(2*MAXL),TRIGS(2*MAXL)
-      INTEGER IFAX(10)
-
-      DO I=MNAUF+1,MAXL-1
-          CXM(2*I)=0.0
-          CXM(2*I+1)=0.0
-      ENDDO
-      CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,1)
-      DO I=0,MAXL-1
-        CXM(I)=CXM(I+1)
-      ENDDO
-
-      RETURN
-      END SUBROUTINE RFOURTR
-C     
-C     
-      SUBROUTINE GAULEG(X1,X2,X,W,N)
-C     BERECHNET DIE GAUSS+SCHEN BREITEN
-        IMPLICIT REAL (A-H,O-Z)
-        DIMENSION X(N),W(N)
-      PARAMETER (EPS=3.D-14)
-      M=(N+1)/2
-      XM=0.5D0*(X2+X1)
-      XL=0.5D0*(X2-X1)
-      DO 12 I=1,M
-        Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
-1       CONTINUE
-          P1=1.D0
-          P2=0.D0
-          DO 11 J=1,N
-            P3=P2
-            P2=P1
-            P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
-11        CONTINUE
-          PP=N*(Z*P1-P2)/(Z*Z-1.D0)
-          Z1=Z
-          Z=Z1-P1/PP
-        IF(ABS(Z-Z1).GT.EPS)GO TO 1
-        X(I)=XM-XL*Z
-        X(N+1-I)=XM+XL*Z
-        W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
-        W(N+1-I)=W(I)
-12    CONTINUE
-      RETURN
-      END SUBROUTINE GAULEG
-C
-C
-        SUBROUTINE PLGNFA(LL,X,Z)
-C
-C PLGNDN BERECHNET ALLE NORMIERTEN ASSOZIIERTEN
-C LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
-C UND SCHREIBT SIE IN DAS FELD Z
-C Die Polynome sind wie im ECMWF indiziert, d.h.
-C P00,P10,P11,P20,P21,P22,...
-C	Ansonsten ist die Routine analog zu PLGNDN
-C X IST DER COSINUS DES ZENITWINKELS ODER
-C       DER SINUS DER GEOGRAFISCHEN BREITE
-C
-        IMPLICIT REAL (A-H,O-Z)
-        DIMENSION Z(0:((LL+3)*(LL+4))/2)
-C
-      L=LL+2
-      I=1
-      Z(0)=1.D0
-      FACT=1.D0
-      POT=1.D0
-      SOMX2=DSQRT(1.D0-X*X)
-        DO 14 J=0,L
-          DJ=DBLE(J)
-        IF(J.GT.0) THEN
-          FACT=FACT*(2.D0*DJ-1.D0)/(2.D0*DJ)
-            POT=POT*SOMX2
-          Z(I)=DSQRT((2.D0*DJ+1.D0)*FACT)*POT
-          I=I+1
-        ENDIF
-          IF(J.LT.L) THEN
-            Z(I)=X*
-     *DSQRT((4.D0*DJ*DJ+8.D0*DJ+3.D0)/(2.D0*DJ+1.D0))*Z(I-1)
-            I=I+1
-          ENDIF
-          DK=DJ+2.D0
-        DO 14 K=J+2,L
-          DDK=(DK*DK-DJ*DJ)
-          Z(I)=X*DSQRT((4.D0*DK*DK-1.D0)/DDK)*Z(I-1)-
-     *    DSQRT(((2.D0*DK+1.D0)*(DK-DJ-1.D0)*(DK+DJ-1.D0))/
-     *    ((2.D0*DK-3.D0)*DDK))*Z(I-2)
-          DK=DK+1.D0
-          I=I+1
-14    CONTINUE
-      RETURN
-      END SUBROUTINE PLGNFA
-
-
-      SUBROUTINE DPLGND(MNAUF,Z,DZ)
-C
-C DPLGND BERECHNET DIE ABLEITUNG DER NORMIERTEN ASSOZIIERTEN
-C LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
-C UND SCHREIBT SIE IN DAS FELD DZ
-C DIE REIHENFOLGE IST
-C P00(X),P01(X),P11(X),P02(X),P12(X),P22(X),..PLL(X)
-C
-        IMPLICIT REAL (A-H,O-Z)
-        DIMENSION Z(0:((MNAUF+3)*(MNAUF+4))/2)
-        DIMENSION DZ(0:((MNAUF+2)*(MNAUF+3))/2)
-C
-				IF(Z(0).NE.1.D0) THEN
-			WRITE(*,*) 'DPLGND: Z(0) must be 1.0'
-					STOP
-				ENDIF
-				
-				LLP=0
-				LLH=0
-        DO 1 I=0,MNAUF+1
-            DO 2 J=I,MNAUF+1
-              IF(I.EQ.J) THEN
-              	WURZELA=
-     *DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
-                DZ(LLH)=DBLE(J)*WURZELA*Z(LLP+1)
-              ELSE
-              		WURZELB=
-     *DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
-                DZ(LLH)=
-     *DBLE(J)*WURZELB*Z(LLP+1)-DBLE(J+1)*WURZELA*Z(LLP-1)
-     						WURZELA=WURZELB
-	        ENDIF
-	        LLH=LLH+1              
-          LLP=LLP+1
-2	      CONTINUE
-        LLP=LLP+1
-1	CONTINUE
-	RETURN
-	END SUBROUTINE DPLGND
-	
-
-* Spectral Filter of Sardeshmukh and Hoskins (1984, MWR)
-* MM=Spectral truncation of field
-* MMAX= Spectral truncation of filter
-*
-      SUBROUTINE SPFILTER(FELDMN,MM,MMAX)
-
-      IMPLICIT NONE
-
-      INTEGER MM,MMAX,I,J,K,L
-      REAL FELDMN(0:(MM+1)*(MM+2)-1)
-      REAL KMAX,SMAX,FAK
-
-      SMAX=0.1
-      KMAX=-ALOG(SMAX)
-      KMAX=KMAX/(float(MMAX)*float(MMAX+1))**2
-c      WRITE(*,*)'alogsmax',alog(smax),'KMAX:',KMAX
-        l=0
-        do i=0,MM
-          do j=i,MM
-c          write(*,*) i,j,feld(k),feld(k)*exp(-KMAX*(j*(j+1))**2)
-            if(j .le. MMAX) then
-c               fak=exp(-KMAX*(j*(j+1))**2)
-               fak=1.0
-               feldmn(2*l)=feldmn(2*l)*fak
-               feldmn(2*l+1)=feldmn(2*l+1)*fak
-            else
-               feldmn(2*l)=0.
-               feldmn(2*l+1)=0.
-            endif
-            l=l+1
-          enddo
-        enddo
-      END SUBROUTINE SPFILTER
-
-      END MODULE PHTOGR
-    
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
Index: Source/Fortran/phgrreal.f90
===================================================================
--- Source/Fortran/phgrreal.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/phgrreal.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,516 @@
+MODULE PHTOGR
+
+  INTEGER, PARAMETER :: MAXAUF=36000
+
+CONTAINS
+
+  SUBROUTINE PHGR213(CXMN,FELD,WSAVE,IFAX,Z,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+
+!! DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
+!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+!! RAUM AUF DAS REDUZIERTE GAUSS'SCHE GITTER DURCH
+!
+! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
+!	WSAVE  = Working Array fuer Fouriertransformation
+! Z 	    = LEGENDREFUNKTIONSWERTE
+!
+! MNAUF  ANZAHL DER FOURIERKOEFFIZIENTEN
+! MAXL   ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+! MAXB   ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+
+    IMPLICIT NONE
+
+!			Anzahl der Gitterpunkte auf jedem Breitenkreis
+    INTEGER MLAT(MAXB/2)
+    INTEGER K,MAXL,MAXB,MLEVEL,MNAUF
+    INTEGER IND(MAXB)
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
+
+    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    REAL WSAVE(8*MAXB+15,MAXB/2)
+    INTEGER :: IFAX(10,MAXB)
+
+    IND(1)=0
+    DO 7 K=2,MAXB/2
+      IND(K)=IND(K-1)+MLAT(K-1)
+7   CONTINUE
+
+!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
+    DO 17 K=1,MAXB/2
+      CALL PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+17  CONTINUE
+!$OMP END PARALLEL DO
+
+    RETURN
+
+  END SUBROUTINE PHGR213
+
+  SUBROUTINE PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
+
+    IMPLICIT NONE
+
+    INTEGER MLAT(MAXB/2)
+    INTEGER K,L,I,J,LLS,LLPS,LL,LLP,MAXL,MAXB,MLEVEL,MNAUF
+    INTEGER IND(MAXB)
+    INTEGER :: IFAX(10,MAXB)
+
+
+!   FELD DER FOURIERKOEFFIZIENTEN
+    REAL :: CXMS(0:MAXAUF-1),CXMA(0:MAXAUF-1)
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
+    REAL ACR,ACI,SCR,SCI
+
+    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    REAL WSAVE(8*MAXB+15,MAXB/2)
+
+    DO 6 L=1,MLEVEL
+      LL=0
+      LLP=0
+      DO 1 I=0,MNAUF
+        SCR=0.D0
+        SCI=0.D0
+        ACR=0.D0
+        ACI=0.D0
+        LLS=LL
+        LLPS=LLP
+        IF (2*I+1 .LT. MLAT(K)) THEN
+!	        Innerste Schleife aufgespalten um if-Abfrage zu sparen
+          DO 18 J=I,MNAUF,2
+            SCR=SCR+Z(LLP,K)*CXMN(2*LL,L)
+            SCI=SCI+Z(LLP,K)*CXMN(2*LL+1,L)
+            LL=LL+2
+            LLP=LLP+2
+18        CONTINUE
+          LL=LLS+1
+          LLP=LLPS+1
+          DO 19 J=I+1,MNAUF,2
+            ACR=ACR+Z(LLP,K)*CXMN(2*LL,L)
+            ACI=ACI+Z(LLP,K)*CXMN(2*LL+1,L)
+            LL=LL+2
+            LLP=LLP+2
+19        CONTINUE
+        END IF
+        LL=LLS+(MNAUF-I+1)
+        LLP=LLPS+(MNAUF-I+3)
+        CXMS(2*I)=SCR+ACR
+        CXMS(2*I+1)=SCI+ACI
+        CXMA(2*I)=SCR-ACR
+        CXMA(2*I+1)=SCI-ACI
+1     CONTINUE
+!     CALL FOURTR(CXMS,FELD(IND(k)+1,L),WSAVE(:,K),MNAUF,*MLAT(K),1)
+!     CALL FOURTR(CXMA,FELD(MAXL-IND(k)-MLAT(K)+1,L),WSAVE(:,K),MNAUF,MLAT(K),1)
+      CALL RFOURTR(CXMS,WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K),1)
+      FELD(IND(k)+1:IND(K)+MLAT(K),L)=CXMS(0:MLAT(K)-1)
+      CALL RFOURTR(CXMA,WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K),1)
+      FELD(MAXL-IND(k)-MLAT(K)+1:MAXL-IND(k),L)=CXMA(0:MLAT(K)-1)
+!      WRITE(*,*) IND+1,FELD(IND+1,L)
+6   CONTINUE
+
+  END SUBROUTINE PHSYM
+
+  SUBROUTINE PHGCUT(CXMN,FELD,WSAVE,IFAX,Z, &
+    MNAUF,MMAX,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+!! DIE ROUTINE FUEHRT EINE TRANSFORMATION EINER
+!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+!! RAUM AUF KUGELKOORDINATEN DURCH. Es kann ein Teilausschnitt
+!!	Der Erde angegeben werden. Diese Routine ist langsamer als phgrph
+
+! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
+! BREITE = SINUS DER GEOGRAFISCHEN BREITEN
+!
+! MNAUF  ANZAHL DER FOURIERKOEFFIZIENTEN
+! MAUF   ANZAHL DER LAENGEN UND DER FOURIERKOEFFIZIENTEN
+! MANF   ANFANG DES LAENGENBEREICHS FUER DAS GITTER,
+!        AUF DAS INTERPOLIERT WERDEN SOLL
+! MAXL   ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+! MAXB   ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+
+    IMPLICIT REAL (A-H,O-Z)
+
+!   FELD DER FOURIERKOEFFIZIENTEN
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MMAX+3)*(MMAX+4))/2,MAXB)
+
+    DIMENSION CXMN(0:(MMAX+1)*(MMAX+2)-1,MLEVEL)
+    REAL FELD(MAXL,MAXB,MLEVEL)
+    DIMENSION WSAVE(4*MAUF+15)
+    INTEGER:: IFAX(10)
+
+    LOGICAL SYM
+
+!    write(*,*)mauf,mnauf,manf,maxl
+
+    IF (MAUF .LE. MNAUF) WRITE(*,*) 'TOO COARSE LONGITUDE RESOLUTION'
+    IF (MANF .LT. 1    .OR. MAXL .LT. 1 .OR. &
+        MANF .GT. MAUF .OR. MAXL .GT. MAUF) THEN
+      WRITE(*,*) 'WRONG LONGITUDE RANGE',MANF,MAXL
+      STOP
+    END IF
+
+! Pruefe, ob Ausgabegitter symmetrisch zum Aequator ist
+! Wenn ja soll Symmetrie der Legendrepolynome ausgenutzt werden
+    IF (MAXB .GT. 4) THEN
+      SYM=.TRUE.
+      DO 11 J=5,5
+        IF (ABS(ABS(Z(100,J))-ABS(Z(100,MAXB+1-J))) .GT. 1E-11) SYM=.FALSE.
+!	      WRITE(*,*) ABS(Z(100,J)),ABS(Z(100,MAXB+1-J))
+11    CONTINUE
+!!      WRITE(*,*) 'Symmetrisch: ',SYM
+    ELSE
+      SYM=.FALSE.
+    END IF
+
+
+    IF (SYM) THEN
+
+!$OMP PARALLEL DO
+      DO J=1,(MAXB+1)/2
+        CALL PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
+      END DO
+!$OMP END PARALLEL DO
+
+    ELSE
+
+!$OMP PARALLEL DO
+      DO J=1,MAXB
+        CALL PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+      END DO
+!$OMP END PARALLEL DO
+
+    END IF
+
+    RETURN
+
+  END SUBROUTINE PHGCUT
+
+  SUBROUTINE PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
+
+    IMPLICIT REAL (A-H,O-Z)
+
+!   FELD DER FOURIERKOEFFIZIENTEN
+
+    REAL :: CXM(0:MAXAUF-1),CXMA(0:MAXAUF-1)
+
+!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
+    REAL SCR,SCI,ACR,ACI
+
+    DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MAXB,MLEVEL)
+    DIMENSION WSAVE(4*MAUF+15)
+    INTEGER :: IFAX(10)
+
+    DO 16 L=1,MLEVEL
+      LL=0
+      LLP=0
+      DO 17 I=0,MNAUF
+        SCR=0.D0
+        SCI=0.D0
+        ACR=0.D0
+        ACI=0.D0
+        LLS=LL
+        LLPS=LLP
+!	      Innerste Schleife aufgespalten um if-Abfrage zu sparen
+        DO 18 K=I,MNAUF,2
+          SCR=SCR+Z(LLP,J)*CXMN(2*LL,L)
+          SCI=SCI+Z(LLP,J)*CXMN(2*LL+1,L)
+          LL=LL+2
+          LLP=LLP+2
+18      CONTINUE
+        LL=LLS+1
+        LLP=LLPS+1
+        DO 19 K=I+1,MNAUF,2
+          ACR=ACR+Z(LLP,J)*CXMN(2*LL,L)
+          ACI=ACI +Z(LLP,J)*CXMN(2*LL+1,L)
+          LL=LL+2
+          LLP=LLP+2
+19      CONTINUE
+        LL=LLS+MNAUF-I+1
+        LLP=LLPS+MNAUF-I+3
+        CXM(2*I)=SCR+ACR
+        CXM(2*I+1)=SCI+ACI
+        CXMA(2*I)=SCR-ACR
+        CXMA(2*I+1)=SCI-ACI
+17    CONTINUE
+
+      CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
+      DO 26 I=0,MAXL-1
+        IF (MANF+I .LE. MAUF) THEN
+          FELD(I+1,J,L)=CXM(MANF+I-1)
+        ELSE
+          FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
+        END IF
+26    CONTINUE
+      CALL RFOURTR(CXMA,WSAVE,IFAX,MNAUF,MAUF,1)
+      DO 36 I=0,MAXL-1
+        IF (MANF+I .LE. MAUF) THEN
+          FELD(I+1,MAXB+1-J,L)=CXMA(MANF+I-1)
+        ELSE
+          FELD(I+1,MAXB+1-J,L)=CXMA(MANF-MAUF+I-1)
+        END IF
+36    CONTINUE
+16  CONTINUE
+
+  END SUBROUTINE PHSYMCUT
+
+  SUBROUTINE PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+    IMPLICIT NONE
+    
+    INTEGER,INTENT(IN) :: MNAUF,MAUF,MANF,J,MAXL,MAXB,MLEVEL
+
+    REAL :: CXM(0:MAXAUF-1)
+    REAL,INTENT(IN) :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
+    REAL,INTENT(IN) :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL,INTENT(IN) :: WSAVE(4*MAUF+15)
+    REAL :: FELD(MAXL,MAXB,MLEVEL)
+
+    INTEGER :: IFAX(10)
+    INTEGER I,L
+
+    DO L=1,MLEVEL
+      CALL LEGTR(CXMN(:,L),CXM,Z(:,J),MNAUF,MAUF)
+      CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
+
+      DO I=0,MAXL-1
+        IF (MANF+I .LE. MAUF) THEN
+          FELD(I+1,J,L)=CXM(MANF+I-1)
+        ELSE
+          FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
+        END IF
+      END DO
+    END DO
+  END SUBROUTINE PHGPNS
+
+  SUBROUTINE LEGTR(CXMN,CXM,Z,MNAUF,MAUF)
+
+!!   DIESE ROUTINE BERECHNET DIE FOURIERKOEFFIZIENTEN CXM
+
+
+    IMPLICIT NONE
+
+    INTEGER MNAUF,MAUF,LL,LLP,I,J
+    REAL CXM(0:MAXAUF-1)
+    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1)
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
+    REAL CI,CR
+
+    LL=0
+    LLP=0
+    DO 1 I=0,MNAUF
+      CR=0.D0
+      CI=0.D0
+      DO 2 J=I,MNAUF
+        CR=CR+Z(LLP)*CXMN(2*LL)
+        CI=CI+Z(LLP)*CXMN(2*LL+1)
+        LL=LL+1
+        LLP=LLP+1
+2     CONTINUE
+      LLP=LLP+2
+      CXM(2*I)=CR
+      CXM(2*I+1)=CI
+1   CONTINUE
+    RETURN
+    
+  END SUBROUTINE LEGTR
+
+  SUBROUTINE RFOURTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
+
+!!     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
+
+    IMPLICIT REAL (A-H,O-Z)
+
+    DIMENSION CXM(0:MAXAUF-1)
+    REAL :: WSAVE(2*MAXL),TRIGS(2*MAXL)
+    INTEGER IFAX(10)
+
+    DO I=MNAUF+1,MAXL-1
+      CXM(2*I)=0.0
+      CXM(2*I+1)=0.0
+    END DO
+
+    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,1)
+
+    DO I=0,MAXL-1
+      CXM(I)=CXM(I+1)
+    END DO
+
+    RETURN
+
+  END SUBROUTINE RFOURTR
+
+  SUBROUTINE GAULEG(X1,X2,X,W,N)
+
+!! BERECHNET DIE GAUSS+SCHEN BREITEN
+
+    IMPLICIT REAL (A-H,O-Z)
+
+    DIMENSION X(N),W(N)
+    PARAMETER (EPS=3.D-14)
+
+    M=(N+1)/2
+    XM=0.5D0*(X2+X1)
+    XL=0.5D0*(X2-X1)
+    DO 12 I=1,M
+      Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
+1     CONTINUE
+      P1=1.D0
+      P2=0.D0
+      DO 11 J=1,N
+        P3=P2
+        P2=P1
+        P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
+11    CONTINUE
+      PP=N*(Z*P1-P2)/(Z*Z-1.D0)
+      Z1=Z
+      Z=Z1-P1/PP
+      IF (ABS(Z-Z1) .GT. EPS)GO TO 1
+      X(I)=XM-XL*Z
+      X(N+1-I)=XM+XL*Z
+      W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
+      W(N+1-I)=W(I)
+12  CONTINUE
+
+    RETURN
+
+  END SUBROUTINE GAULEG
+
+
+  SUBROUTINE PLGNFA(LL,X,Z)
+
+!! PLGNDN BERECHNET ALLE NORMIERTEN ASSOZIIERTEN
+!! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
+!! UND SCHREIBT SIE IN DAS FELD Z
+! Die Polynome sind wie im ECMWF indiziert, d.h.
+! P00,P10,P11,P20,P21,P22,...
+!	Ansonsten ist die Routine analog zu PLGNDN
+! X IST DER COSINUS DES ZENITWINKELS ODER
+!       DER SINUS DER GEOGRAFISCHEN BREITE
+
+    IMPLICIT REAL (A-H,O-Z)
+
+    DIMENSION Z(0:((LL+3)*(LL+4))/2)
+
+    L=LL+2
+    I=1
+    Z(0)=1.D0
+    FACT=1.D0
+    POT=1.D0
+    SOMX2=DSQRT(1.D0-X*X)
+    DO 14 J=0,L
+      DJ=DBLE(J)
+      IF (J .GT. 0) THEN
+        FACT=FACT*(2.D0*DJ-1.D0)/(2.D0*DJ)
+        POT=POT*SOMX2
+        Z(I)=DSQRT((2.D0*DJ+1.D0)*FACT)*POT
+        I=I+1
+      END IF
+      IF (J .LT. L) THEN
+        Z(I)=X*DSQRT((4.D0*DJ*DJ+8.D0*DJ+3.D0)/(2.D0*DJ+1.D0))*Z(I-1)
+        I=I+1
+      END IF
+      DK=DJ+2.D0
+      DO 14 K=J+2,L
+        DDK=(DK*DK-DJ*DJ)
+        Z(I)=X*DSQRT((4.D0*DK*DK-1.D0)/DDK)*Z(I-1)- &
+          DSQRT(((2.D0*DK+1.D0)*(DK-DJ-1.D0)*(DK+DJ-1.D0))/ &
+          ((2.D0*DK-3.D0)*DDK))*Z(I-2)
+        DK=DK+1.D0
+        I=I+1
+14  CONTINUE
+
+    RETURN
+
+  END SUBROUTINE PLGNFA
+
+  SUBROUTINE DPLGND(MNAUF,Z,DZ)
+
+!! DPLGND BERECHNET DIE ABLEITUNG DER NORMIERTEN ASSOZIIERTEN
+!! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
+!! UND SCHREIBT SIE IN DAS FELD DZ
+! DIE REIHENFOLGE IST
+! P00(X),P01(X),P11(X),P02(X),P12(X),P22(X),..PLL(X)
+
+    IMPLICIT REAL (A-H,O-Z)
+
+    DIMENSION Z(0:((MNAUF+3)*(MNAUF+4))/2)
+    DIMENSION DZ(0:((MNAUF+2)*(MNAUF+3))/2)
+
+    IF (Z(0) .NE. 1.D0) THEN
+      WRITE(*,*) 'DPLGND: Z(0) must be 1.0'
+      STOP
+    END IF
+
+    LLP=0
+    LLH=0
+    DO 1 I=0,MNAUF+1
+      DO 2 J=I,MNAUF+1
+        IF (I .EQ. J) THEN
+          WURZELA=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
+          DZ(LLH)=DBLE(J)*WURZELA*Z(LLP+1)
+        ELSE
+          WURZELB=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
+          DZ(LLH)=DBLE(J)*WURZELB*Z(LLP+1)-DBLE(J+1)*WURZELA*Z(LLP-1)
+          WURZELA=WURZELB
+        END IF
+        LLH=LLH+1
+        LLP=LLP+1
+2     CONTINUE
+      LLP=LLP+1
+1   CONTINUE
+
+    RETURN
+
+  END SUBROUTINE DPLGND
+
+
+  SUBROUTINE SPFILTER(FELDMN,MM,MMAX)
+
+!! Spectral Filter of Sardeshmukh and Hoskins (1984, MWR)
+! MM=Spectral truncation of field
+! MMAX= Spectral truncation of filter
+
+    IMPLICIT NONE
+
+    INTEGER MM,MMAX,I,J,K,L
+    REAL FELDMN(0:(MM+1)*(MM+2)-1)
+    REAL KMAX,SMAX,FAK
+
+    SMAX=0.1
+    KMAX=-ALOG(SMAX)
+    KMAX=KMAX/(float(MMAX)*float(MMAX+1))**2
+!    WRITE(*,*)'alogsmax',alog(smax),'KMAX:',KMAX
+    L=0
+    DO I=0,MM
+      DO J=I,MM
+!        WRITE(*,*) I,J,FELD(K),FELD(K)*EXP(-KMAX*(J*(J+1))**2)
+        IF(J .LE. MMAX) THEN
+!          FAK=EXP(-KMAX*(J*(J+1))**2)
+          FAK=1.0
+          FELDMN(2*L)=FELDMN(2*L)*FAK
+          FELDMN(2*L+1)=FELDMN(2*L+1)*FAK
+        ELSE
+          FELDMN(2*L)=0.
+          FELDMN(2*L+1)=0.
+        END IF
+        L=L+1
+      END DO
+    END DO
+    
+  END SUBROUTINE SPFILTER
+
+END MODULE PHTOGR
Index: urce/Fortran/posnam.f
===================================================================
--- Source/Fortran/posnam.f	(revision None)
+++ 	(revision )
@@ -1,25 +1,0 @@
-      SUBROUTINE POSNAM(KULNAM,CDNAML)
-!-------------------------------------
-
-!--- position in namelist file. 
-!    author:  Mats Hamrud, ECMWF
-
-      INTEGER, INTENT(IN)       :: KULNAM
-      CHARACTER*(*), INTENT(IN) :: CDNAML
-      CHARACTER*120 CLINE
-      CHARACTER*1 CLTEST
-      REWIND(KULNAM)
-      ILEN=LEN(CDNAML)
- 102  CONTINUE
-      CLINE=' '
-      READ(KULNAM,'(A)') CLINE
-      IND1=INDEX(CLINE,'&'//CDNAML)
-      IF(IND1.EQ.0) GO TO 102
-      CLTEST=CLINE(IND1+ILEN+1:IND1+ILEN+1)
-      IF((LGE(CLTEST,'0').AND.LLE(CLTEST,'9')).OR.
-     &   (LGE(CLTEST,'A').AND.LLE(CLTEST,'Z'))) GO TO 102
-      BACKSPACE(KULNAM)
-
-      RETURN
-      END SUBROUTINE POSNAM
-
Index: Source/Fortran/posnam.f90
===================================================================
--- Source/Fortran/posnam.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Source/Fortran/posnam.f90	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,25 @@
+  SUBROUTINE POSNAM(KULNAM,CDNAML)
+
+ !! position in namelist file.
+ ! author:  Mats Hamrud, ECMWF
+
+    INTEGER, INTENT(IN)       :: KULNAM
+    CHARACTER*(*), INTENT(IN) :: CDNAML
+    CHARACTER*120 CLINE
+    CHARACTER*1 CLTEST
+    
+    REWIND(KULNAM)
+    ILEN = LEN(CDNAML)
+102 CONTINUE
+      CLINE = ' '
+      READ(KULNAM,'(A)') CLINE
+      IND1 = INDEX(CLINE,'&'//CDNAML)
+      IF (IND1 .EQ. 0) GO TO 102
+      CLTEST = CLINE(IND1+ILEN+1:IND1+ILEN+1)
+    IF (LGE(CLTEST,'0') .AND. LLE(CLTEST,'9') .OR. &
+        LGE(CLTEST,'A') .AND. LLE(CLTEST,'Z')) GOTO 102
+    BACKSPACE(KULNAM)
+
+    RETURN
+    
+  END SUBROUTINE POSNAM
Index: urce/Fortran/preconvert.f90
===================================================================
--- Source/Fortran/preconvert.f90	(revision None)
+++ 	(revision )
@@ -1,807 +1,0 @@
-      PROGRAM PRECONVERT
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                                                                 !
-! PROGRAM PRECONVERT - PREPARES INPUT DATA FOR POP MODEL METEOR-  !
-!                      OLOGICAL PREPROCESSOR                      !
-!                                                                 !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                                                                 !
-! CALCULATION OF ETAPOINT ON A REGULAR LAMDA/PHI GRID AND WRITING !
-! U,V,ETAPOINT,T,PS,Q,SD,MSL,TCC,10U, 10V, 2T,2D,LSP,CP,SSHF,SSR, !
-! EWSS,NSSS TO AN OUTPUT FILE (GRIB 1 or 2 FORMAT).               ! 
-!                                                                 !
-! AUTHORS: L. HAIMBERGER, G. WOTAWA, 1994-04                      !
-!                     adapted: A. BECK                            !
-!                     2003-05-11                                  !
-!          L. Haimberger 2006-12    V2.0                          !
-!                    modified to handle arbitrary regular grids   !
-!                    and T799 resolution data                     !
-!          L. Haimberger 2010-03    V4.0                          !
-!                    modified to grib edition 2 fields            !
-!                    and T1279 resolution data                    !
-!                                                                 !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                                                                 !
-! DESCRIPTION OF NEEDED INPUT:                                    !
-!                                                                 !
-! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
-!                                                                 !
-! 11    fort.11   T,U,V           regular lamda phi grid          !
-! 12    fort.12   D               regular lamda phi grid          !
-! 13    fort.13   LNSP          fort.13  spherical harmonics             !
-! 14    fort.14   SD,MSL,TCC,10U,                                 !
-!                 10V,2T,2D       regular lamda phi grid          !
-! 16    fort.16   LSP,CP,SSHF,                                    !
-!                 SSR,EWSS,NSSS   regular lamda phi grid          !
-! 17    fort.17   Q               regular lamda phi grid          !
-!                                                                 !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                                                                 !
-! DESCRIPTION OF OUTPUT:                                          !
-!                                                                 !
-! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
-!                                                                 !
-! 15    fort.15   U,V,ETA,T,PS,                                   !
-!                 Q,SD,MSL,TCC,                                   !
-!                 10U,10V,2T,2D,  regular lamda phi grid          !
-!                 LSP,CP,SSHF,                                    !
-!                 SSR,EWSS,NSSS                                   !
-!                                                                 !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-
-      USE PHTOGR
-      USE GRTOPH
-      USE FTRAFO
-      USE RWGRIB2
-      USE GRIB_API
-
-      IMPLICIT NONE
-
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: LNPS
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: Z
-      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T, UV , UV2
-      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: QA,OM,OMR
-      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: DIV, ETA,ETAR
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: DPSDL, DPSDM
-      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: PS,DPSDT
-      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: SURF,FLUX,OROLSM
-      REAL, ALLOCATABLE, DIMENSION (:) :: WSAVE,H,SINL,COSL,WSAVE2
-      REAL, ALLOCATABLE, DIMENSION (:) :: BREITE, GBREITE,AK, BK,pv
-
-! Arrays for Gaussian grid calculations
-
-      REAL  :: X1,X2,RMS,MW,SIG,LAM
-      REAL,ALLOCATABLE :: CUA(:,:,:),CVA(:,:,:)
-
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: P,PP !,P2
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: XMN,HILFUV
-      REAL, ALLOCATABLE, DIMENSION (:) :: LNPMN,LNPMN2,LNPMN3
-      REAL, ALLOCATABLE, DIMENSION (:) :: WEIGHT
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: UGVG
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: DG, ETAG
-      REAL, ALLOCATABLE, DIMENSION (:,:) :: GWSAVE
-      REAL, ALLOCATABLE, DIMENSION (:) :: PSG,HILF
-
-! end arrays for Gaussian grid calculations
-
-      INTEGER, ALLOCATABLE, DIMENSION (:) :: MLAT,MPSURF,MPFLUX,MPORO,MPAR
-      INTEGER, ALLOCATABLE :: GIFAX(:,:)
-
-      REAL PI,COSB,DAK,DBK,P00
-      REAL URLAR8,JMIN1,LLLAR8,MAXBMIN1,PIR8,DCOSB
-
-      INTEGER I,J,K,L,IERR,M,LTEST,MK,NGI,NGJ
-      INTEGER MFLUX,MSURF,MORO
-      INTEGER LUNIT,LUNIT2
-
-      INTEGER MAXL, MAXB, MLEVEL, LEVOUT,LEVMIN,LEVMAX
-      INTEGER MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH, MNAUF,META,METADIFF
-      INTEGER MDPDETA,METAPAR
-      REAL RLO0, RLO1, RLA0, RLA1
-      CHARACTER*300 MLEVELIST
-
-      INTEGER MAUF, MANF,IFAX(10)
-
-      INTEGER IGRIB(1),iret,ogrib
-
-      CHARACTER*80 FILENAME
-
-      NAMELIST /NAMGEN/ &
-                 MAXL, MAXB, &
-                 MLEVEL,MLEVELIST,MNAUF,METAPAR, &
-                 RLO0, RLO1, RLA0, RLA1, &
-                 MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH,META,METADIFF,&
-                 MDPDETA
-
-      LTEST=1
-
-      call posnam (4,'NAMGEN')
-      read (4,NAMGEN)
-
-      MAUF=INT(360.*(REAL(MAXL)-1.)/(RLO1-RLO0)+0.0001)
-!      PRINT*, MAUF
-
-      MANF=INT(REAL(MAUF)/360.*(360.+RLO0)+1.0001)
-      IF(MANF .gt. MAUF) MANF=MANF-MAUF
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                       ALLOCATE VARIABLES                        !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      ALLOCATE (LNPS(0:(MNAUF+1)*(MNAUF+2)-1,1))
-
-      ALLOCATE (H(0:(MNAUF+2)*(MNAUF+3)/2))
-
-     
-      ALLOCATE (OM(MAXL, MAXB, MLEVEL))
-
-      ALLOCATE (ETA(MAXL,MAXB,MLEVEL))
-
-      ALLOCATE (PS(MAXL, MAXB,1),DPSDT(MAXL, MAXB,1))
-
-      
-      ALLOCATE (WSAVE(4*MAUF+15),WSAVE2(4*MAUF+15))
-
-      ALLOCATE (BREITE(MAXB),AK(MLEVEL+1),BK(MLEVEL+1),pv(2*mlevel+2))
-      
-      ALLOCATE (MPAR(2))
-
-      ALLOCATE (COSL(MAXL),SINL(MAXL))
-		
-      ALLOCATE (CUA(2,4,MLEVEL),CVA(2,4,MLEVEL))
-    	
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!     GAUSS STUFF !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-      IF(MGAUSS .EQ. 1) THEN
-      LUNIT=0
-      FILENAME='fort.18'
-
-      call grib_open_file(LUNIT, TRIM(FILENAME),'r')
- 
-      call grib_new_from_file(LUNIT,igrib(1), iret)
- 
-! we can close the file
-      call grib_close_file(LUNIT)
- 
-!      call grib_get(igrib(1),'gridType', j)
-
-      NGJ=MNAUF+1
-
-      ALLOCATE (GWSAVE(8*NGJ+15,NGJ/2))
-      ALLOCATE(GIFAX(10,NGJ))
-      ALLOCATE (GBREITE(NGJ),WEIGHT(NGJ))
-      ALLOCATE (MLAT(NGJ))
-      ALLOCATE (P(0:((MNAUF+3)*(MNAUF+4))/2,NGJ/2))
-      ALLOCATE (PP(NGJ/2,0:((MNAUF+3)*(MNAUF+4))/2))
-      ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB))
-
-      call grib_get(igrib(1),'numberOfPointsAlongAMeridian', NGJ)
- 
-      !     get as a integer
-      call grib_get(igrib(1),'pl', MLAT)
-
-      NGI=SUM(MLAT)
-
-      call grib_get(igrib(1),'numberOfVerticalCoordinateValues',mk)
-
-      IF(mk/2-1 .ne. MLEVEL) THEN 
-        WRITE(*,*) 'FATAL: Number of model levels',mk, &
-          ' does not agree with', MLEVEL,' in namelist'
-        STOP
-      ENDIF
-      call grib_get(igrib(1),'pv',pv)
-        AK=pv(1:1+MLEVEL)
-        BK=pv(2+MLEVEL:2*MLEVEL+2)
-
-      ALLOCATE (LNPMN(0:(MNAUF+1)*(MNAUF+2)-1))
-      ALLOCATE (LNPMN2(0:(MNAUF+1)*(MNAUF+2)-1))
-      ALLOCATE (UGVG(NGI, 2*MLEVEL),HILFUV(2*MAXL,2))
-
-
-      ALLOCATE (DPSDL(NGI,1),DPSDM(NGI,1))
-
-      ALLOCATE (PSG(NGI),HILF(NGI))
-      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
-!      ALLOCATE (UV2(MAXL, MAXB, 2*MLEVEL))
-
-      ALLOCATE (XMN(0:(MNAUF+1)*(MNAUF+2)-1, 2*MLEVEL))
-      ALLOCATE (DG(NGI,MLEVEL),ETAG(NGI,MLEVEL))
-
-! Initialisieren  Legendretransformation
-!	auf das LaT/LON Gitter
-
-      PI=ACOS(-1.D0)
-!$OMP PARALLEL DO
-      DO 20 J=1,MAXB
-
-      BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))* PI/180.D0)
-
-      CALL PLGNFA(MNAUF,BREITE(J),Z(0,J))
-
-20    CONTINUE
-!$OMP END PARALLEL DO
-
-! Avoid possible Pole problem
-!      IF(RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
-!      IF(RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
-
-! Initialisation of fields for  FFT and Legendre transformation
-!	to  Gaussian grid and back to phase space
-	X1=-1.D0
-	X2=1.D0
-	CALL GAULEG(X1,X2,GBREITE,WEIGHT,NGJ)
-
-!$OMP PARALLEL DO PRIVATE(M)
-	DO J=1,NGJ/2
-               CALL PLGNFA(MNAUF,GBREITE(J),P(:,J))
-		DO M=0,(MNAUF+3)*(MNAUF+4)/2
-		  PP(J,M)=P(M,J)
-	ENDDO
-	ENDDO
-!$OMP END PARALLEL DO
-
-
-!       MPAR(1)=152
-        FILENAME='fort.12'
-        CALL READSPECTRAL(FILENAME,LNPMN,MNAUF,1,MLEVEL,(/152/),AK,BK)
-!      goto 111
-        CALL SET99(WSAVE,IFAX,mauf)
-        CALL PHGCUT(LNPMN,PS,WSAVE,IFAX,Z, &
-      MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
-      CALL STATIS(MAXL,MAXB,1,EXP(PS),RMS,MW,SIG)
-      WRITE(*,'(A12,3F12.4)') 'STATISTICS: ',RMS,MW,SIG
-
-        DO J=1,NGJ/2
-          CALL SET99(GWSAVE(1,J),GIFAX(1,J),MLAT(J))
-        ENDDO
-	CALL PHGR213(LNPMN,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
-        PSG=HILF
-	CALL GRPH213(LNPMN2,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
-      MNAUF,NGI,NGJ,1)
-	CALL PHGR213(LNPMN2,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
-
-
-        HILF=exp(PSG)-exp(HILF)
-
-      CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
-      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
-
-        PSG=EXP(PSG)
-        HILF=PSG
-      CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
-      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
-
-  111         FILENAME='fort.10'
-       CALL READSPECTRAL(FILENAME, &
-      XMN,MNAUF,2*MLEVEL,MLEVEL,(/131,132/),AK,BK) 
-!	Transformieren des Windes auf das Gaussgitter	
-	CALL PHGR213(XMN,UGVG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,2*MLEVEL)
-       DO K=1,MLEVEL
-! North Pole
-            CALL JSPPOLE(XMN(:,K),1,MNAUF,.TRUE.,CUA(:,:,K))
-            CALL JSPPOLE(XMN(:,MLEVEL+K),1,MNAUF,.TRUE.,CVA(:,:,K))
-! South Pole
-            CALL JSPPOLE(XMN(:,K),-1,MNAUF,.TRUE.,CUA(:,3:4,K))
-            CALL JSPPOLE(XMN(:,MLEVEL+K),-1,MNAUF,.TRUE.,CVA(:,3:4,K))
-       ENDDO
-    
-        DO K=1,2*MLEVEL
-          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
-        ENDDO
-        CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z, &
-      MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
-
-
- 112        FILENAME='fort.13'
-      CALL READSPECTRAL(FILENAME,XMN,MNAUF,MLEVEL,MLEVEL,(/155/),AK,BK)
-!	Transformieren der horizontalen Divergenz auf das Gaussgitter			
-	CALL PHGR213(XMN,DG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,MLEVEL)
-
-
-!	Berechnung des Gradienten des Logarithmus des Bodendrucks 
-!       auf dem Gaussgitter	
-	CALL PHGRAD(LNPMN,DPSDL,DPSDM,GWSAVE,GIFAX,P,H,MLAT,MNAUF,NGI,NGJ,1)
-
-!	Berechnung der Vertikalgeschwindigkeit auf dem Gaussgitter
-      CALL CONTGL(HILF,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
-      GBREITE,ETAG,MLAT,AK,BK,NGI,NGJ,MLEVEL)
-
- 
- 	CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
-      MNAUF,NGI,NGJ,MLEVEL)
-        DO K=1,MLEVEL
-          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
-        ENDDO
-        CALL PHGCUT(XMN,ETA,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
-
-	CALL GRPH213(XMN,HILF,GWSAVE,GIFAX,PP,WEIGHT,MLAT, MNAUF,NGI,NGJ,1)
-
-        IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,1),MNAUF,MSMOOTH)
-        CALL PHGCUT(XMN,DPSDT,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
-!       GOTO 114
-
-      CALL STATIS(MAXL,MAXB,1,DPSDT,RMS,MW,SIG)
-      WRITE(*,'(A12,3F11.4)') 'STATISTICS DPSDT: ',RMS,MW,SIG
-
-       IF(MOMEGADIFF .ne. 0) THEN
-!	Berechnung von Omega auf dem Gaussgitter
-	CALL OMEGA(PSG,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
-      GBREITE,ETAG,MLAT,AK,BK,NGI ,NGJ,MLEVEL) 
-
-	CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,&
-      MNAUF,NGI,NGJ,MLEVEL)
-        DO K=1,MLEVEL
-          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
-        ENDDO
-        CALL PHGCUT(XMN,OM,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
-
-       ENDIF !MOMEGA
-
-	CALL GRPH213(XMN,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
-        CALL PHGCUT(XMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
-
-      CALL STATIS(MAXL,MAXB,1,PS,RMS,MW,SIG)
-      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
-
- 114  DEALLOCATE(HILF,PSG,DPSDL,DPSDM,ETAG,DG,LNPMN)
-
-!      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
-!	CALL GRPH213(XMN,UGVG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,
-!     *MNAUF,NGI,NGJ,2*MLEVEL)
-!        DO K=1,2*MLEVEL
-!          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
-!        ENDDO
-!        CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z,
-!     *MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
-        DEALLOCATE(PP,P,UGVG,MLAT,GBREITE,WEIGHT,GWSAVE,XMN)
-
-!        CALL ETAGAUSS(Z,WSAVE
-!     *,BREITE,UV,ETA,OM,PS,
-!     *MAUF,MAXB,MAXL,MANF,MNAUF,MLEVEL,MSMOOTH)
-
-      ELSE
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!          READING OF PREPARED METEOROLOGICAL FIELDS              !
-!                                                                 !
-!          THE FOLLOWING FIELDS ARE EXPECTED:                     !
-!                                                                 !
-!          UNIT 11: T,U,V        (REGULAR GRID)                   !
-!          UNIT 17: Q            (REGULAR GRID)                   !
-!          UNIT 13: D            (REGULAR GRID)                   !
-!          UNIT 12: LNSP         (SPHERICAL HARMONICS)            !
-!          UNIT 14: SURFACE DATA (REGULAR GRID)                   !
-!          UNIT 16: FLUX DATA    (REGULAR GRID)                   !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-      ALLOCATE (MLAT(MAXB))
-      MLAT=MAXL
-      ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,1))
-      ALLOCATE (DPSDL(MAXL,MAXB),DPSDM(MAXL,MAXB))
-      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL),DIV(MAXL,MAXB,MLEVEL))
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                  READING OF SURFACE PRESSURE                    !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      FILENAME='fort.12'
-        CALL READSPECTRAL(FILENAME,LNPS,MNAUF,1,MLEVEL,(/152/),AK,BK)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                      READING OF U,V                      !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-! OPENING OF UNBLOCKED GRIB FILE
-!
-      FILENAME='fort.10'
-      CALL READLATLON(FILENAME,UV,MAXL,MAXB,2*MLEVEL,(/131,132/))
-
-
-      PI=ACOS(-1.D0)
-      DO J=1,MAXB
-
-        BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))*PI/180.D0)
-
-      ENDDO
-! Avoid possible Pole problem
-!      IF(RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
-!      IF(RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
-
-      DO K=1,2*MLEVEL
-        DO J=1,MAXB
-          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
-          IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
-             RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
-            UV(:,J,K)=UV(:,J,K)/1.D6
-          else
-            UV(:,J,K)=UV(:,J,K)*COSB
-          endif
-      ENDDO
-      ENDDO
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     READING OF LNSP on grid                !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! For debugging only
-!      FILENAME='LNSPG_G.20060330.600'
-!      INQUIRE(FILE=FILENAME,EXIST=EX)
-!      CALL READLATLON(FILENAME,QA,
-!     *MAXL,MAXB,1,1,(/152/))
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     READING OF DIVERGENCE                       !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      IF(META .EQ. 0 .OR. METADIFF .EQ. 1) THEN
-        FILENAME='fort.13'
-        CALL READLATLON(FILENAME,DIV,MAXL,MAXB,MLEVEL,(/155/))
-      ENDIF
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!       CALCULATION OF ETAPOINT --> TOTAL TIME DERIVATIVE OF       !
-!      ECMWF VERTICAL COORDINATE ETA MULTIPLIED BY DERIVATIVE     !
-!      OF PRESSURE IN ETA DIRECTION                               !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-! Initialisieren  Legendretransformation
-!	auf das LaT/LON Gitter
-! Without Gaussian grid calculation Legendre Polynomials are calculated
-! only for one latitude to save space
-
-      DO J=1,MAXB
-
-        CALL PLGNFA(MNAUF,BREITE(J),Z(0,1))
-
-        CALL PHGCUT(LNPS,PS(:,J,1),WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,1,1)
-
-
-        IF(META .EQ. 0 .or. METADIFF .EQ. 1 ) THEN
-          CALL PHGRACUT(LNPS,DPSDL(:,J),DPSDM(:,J),WSAVE,IFAX,Z,H,MAUF, &
-      MNAUF,MAXL,1,MANF,1)
-        ENDIF
-      ENDDO
-     
-      PS=EXP(PS)
-
-! For debugging only
-      CALL STATIS(MAXL,MAXB,1,PS(:,:,1),RMS,MW,SIG)
-      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
-      
-      
-       IF(MOMEGADIFF .ne. 0) THEN
-
-         CALL OMEGA(PS,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
-      BREITE,OM,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL) 
-       ENDIF
-     
-       IF(META .EQ. 0 .OR. METADIFF .ne. 0) THEN
-         DPSDT=PS
-         CALL CONTGL(DPSDT,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
-      BREITE,ETA,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL)
-       ENDIF
-
-      ENDIF ! MGAUSS
-
-! CREATE FILE VERTICAL.EC NEEDED BY POP MODEL
-
-      open(21,file='VERTICAL.EC')
-      write(21,'(a)')
-      write(21,'(a)') 'VERTICAL DISCRETIZATION OF POP MODEL'
-      write(21,'(a)')
-      write(21,'(i3,a)') MLEVEL,'   number of layers'
-      write(21,'(a)')
-      write(21,'(a)') '* A(NLEV+1)'
-      write(21,'(a)')
-      do 205 i=1,MLEVEL+1
-205      write(21,'(f18.12)') AK(I)
-      write(21,'(a)')
-      write(21,'(a)') '* B(NLEV+1)'
-      write(21,'(a)')
-      do 210 i=1,MLEVEL+1
-210      write(21,'(f18.12)') BK(I)
-      close(21)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     READING OF OMEGA                       !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      IF(MOMEGA .NE. 0 ) THEN
-
-
-
-         ALLOCATE (OMR(MAXL, MAXB, MLEVEL))
- 
-         FILENAME='fort.19'
-         CALL READLATLON(FILENAME,OMR,MAXL,MAXB,MLEVEL,(/135/))
-
-      IF(MOMEGADIFF .NE. 0 ) THEN
-
-      DO K=1,MLEVEL
-        CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)        
-        WRITE(*,'(A12,I3,3F11.4)') '       ETA: ',K,RMS,MW,SIG
-        CALL STATIS(MAXL,MAXB,1,OMR(:,:,K),RMS,MW,SIG)        
-        WRITE(*,'(A12,I3,3F11.4)') '     OMEGA: ',K,RMS,MW,SIG
-        CALL STATIS(MAXL,MAXB,1,OM(:,:,K)-OMR(:,:,K),RMS,MW,SIG)        
-        WRITE(*,'(A12,I3,3F11.4)') 'OMEGA DIFF: ',K,RMS,MW,SIG
-      ENDDO
-
-      ENDIF
-      ENDIF
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     READING OF ETA                       !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      IF(META .NE. 0 ) THEN
-
-         ALLOCATE (ETAR(MAXL, MAXB, MLEVEL))
- 
-         P00=101325.
-         FILENAME='fort.21'
-         CALL READLATLON(FILENAME,ETAR,MAXL,MAXB,MLEVEL,(/77/))
-
-         if(MDPDETA .EQ. 1) THEN
-          DO K=1,MLEVEL
-           DAK=AK(K+1)-AK(K)
-           DBK=BK(K+1)-BK(K)
-           DO J=1,MAXB
-             DO I=1,MAXL
-               ETAR(I,J,K)=2*ETAR(I,J,K)*PS(I,J,1)*(DAK/PS(I,J,1)+DBK)/ &
-                          (DAK/P00+DBK)
-               IF(K .GT. 1) ETAR(I,J,K)=ETAR(I,J,K)-ETAR(I,J,K-1)
-             ENDDO
-           ENDDO
-          ENDDO
-         ENDIF
-
-        IF(METADIFF .NE. 0 ) THEN
-
-         DO K=1,MLEVEL
-          CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)        
-          WRITE(*,'(A12,I3,3F11.4)') '       ETA: ',K,RMS,MW,SIG
-          CALL STATIS(MAXL,MAXB,1,ETAR(:,:,K),RMS,MW,SIG)        
-          WRITE(*,'(A12,I3,3F11.4)') '     ETAR: ',K,RMS,MW,SIG
-          CALL STATIS(MAXL,MAXB,1,ETA(:,:,K)-ETAR(:,:,K),RMS,MW,SIG)        
-          WRITE(*,'(A12,I3,3F11.4)') 'ETA DIFF: ',K,RMS,MW,SIG
-         ENDDO
-         DO K=1,MLEVEL
-          WRITE(*,'(I3,2F11.4)') K,ETA(1,MAXB/2,K),ETAR(1,MAXB/2,K)
-         ENDDO
-        ELSE
-          ETA=ETAR
-        ENDIF
-      ENDIF
-
-      ALLOCATE (T(MAXL, MAXB, MLEVEL))
-      ALLOCATE (QA(MAXL, MAXB, MLEVEL))
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                      READING OF T                      !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-! OPENING OF UNBLOCKED GRIB FILE
-!
-      FILENAME='fort.11'
-      CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/))
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     READING OF SPECIFIC HUMIDITY                !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      FILENAME='fort.17'
-      CALL READLATLON(FILENAME,QA,MAXL,MAXB,MLEVEL,(/133/))
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                     TEST READING OF UV from MARS (debug only)   !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!      FILENAME='fort.22'
-!      CALL READLATLON(FILENAME,UV2,MAXL,MAXB,2*MLEVEL,2,(/131,132/))
-
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                    WRITE MODEL LEVEL DATA TO fort.15            !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!     Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi)
-!     Now we are transforming back to the usual winds. 
-      DO K=1,MLEVEL
-        DO J=2,MAXB-1
-          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
-          UV(:,J,K)=UV(:,J,K)/COSB
-          UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
-        ENDDO
-! special treatment for poles, if necessary. 
-        DO J=1,MAXB,MAXB-1
-          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
-          if(1.0-BREITE(J)*BREITE(J) .gt. 0 .OR. MGAUSS .NE. 1) then
-            IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
-             RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
-             UV(:,J,K)=UV(:,J,K)*1.D6
-             UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6
-            else
-             UV(:,J,K)=UV(:,J,K)/COSB
-             UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
-            endif
-          else
-              HILFUV(5:MAXL,:)=0.
-              HILFUV(1:2,:)=0.
-            IF(J.EQ.MAXB) THEN
-! Suedpol
-              HILFUV(3:4,1)=CUA(:,4,K)
-              HILFUV(3:4,2)=CVA(:,4,K)
-            ELSE
-! Nordpol
-              HILFUV(3:4,1)=CUA(:,2,K)
-              HILFUV(3:4,2)=CVA(:,2,K)
-            ENDIF
-              CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
-              DO I=0,MAXL-1
-                IF(MANF+I.LE.MAXL) THEN
-                  UV(I+1,J,K)=HILFUV(MANF+I,1)
-                ELSE
-                  UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1)
-                ENDIF
-              ENDDO
-              CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
-              DO I=0,MAXL-1
-                IF(MANF+I.LE.MAXL) THEN
-                  UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2)
-                ELSE
-                  UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2)
-                ENDIF
-              ENDDO
-          endif
-      ENDDO
-      ENDDO
-
-! open output file
-      call grib_open_file(LUNIT,'fort.15','w')
-
-! we use temperature on lat/lon on model levels as template for model level data
-      LUNIT2=0
-      call grib_open_file(LUNIT2,'fort.11','r')
-      call grib_new_from_file(LUNIT2,igrib(1), iret)
-      call grib_close_file(LUNIT2)
-
-
-      CALL WRITELATLON(LUNIT,igrib(1),ogrib,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/131/))
-
-      CALL WRITELATLON(LUNIT,igrib(1),ogrib,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/132/))
-
-      IF(MDPDETA .ne. 1 .AND. MGAUSS .EQ. 0 .and. META .eq. 1) THEN
-        CALL WRITELATLON(LUNIT,igrib(1),ogrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/77/))
-      ELSE
-        CALL WRITELATLON(LUNIT,igrib(1),ogrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/METAPAR/))
-      ENDIF
-
-      CALL WRITELATLON(LUNIT,igrib(1),ogrib,T,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/130/))
-
-      CALL WRITELATLON(LUNIT,igrib(1),ogrib,PS,MAXL,MAXB,1,'1',1,(/134/))
- 
-      call grib_set(igrib(1),"levelType","ml")
-      call grib_set(igrib(1),"typeOfLevel","hybrid")
-      CALL WRITELATLON(LUNIT,igrib(1),ogrib,QA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/133/))
-
-  
-      IF(MOMEGA .EQ. 1) THEN
-        call grib_open_file(LUNIT2,'fort.25','w')
-        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OMR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
-
-        IF(MOMEGADIFF .EQ. 1) THEN
-
-          CALL WRITELATLON(LUNIT2,igrib(1),ogrib,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
-
-        OM=OM-OMR
-        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
-      call grib_close_file(LUNIT2)
-        ENDIF
-      ENDIF
-
-      IF(META .EQ. 1 .and. METADIFF .EQ. 1) THEN
-        call grib_open_file(LUNIT2,'fort.26','w')
-        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,ETAR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
-
-!        IF(MOMEGADIFF .EQ. 1) THEN
-
-          CALL WRITELATLON(LUNIT2,igrib(1),ogrib,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
-
-        OM=ETA-ETAR
-        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
-      call grib_close_file(LUNIT2)
-!        ENDIF
-      ENDIF
-
-
-      call grib_close_file(LUNIT)
-
-
-
- 2000 STOP 'SUCCESSFULLY FINISHED CONVERT_PRE: CONGRATULATIONS'
- 3000 STOP 'ROUTINE CONVERT_PRE: ERROR'
- 9999 stop 'ROUTINE CONVERT_PRE: ERROR'
-      END
-
-      
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      INTEGER FUNCTION IA (FIELD1,NI,NJ,NK,G)
-
-      IMPLICIT NONE
-      INTEGER NI,NJ,NK,I,J,K
-      REAL FIELD1(NI,NJ,NK)
-      REAL G
-      REAL RMIN,RMAX,XMAX,A,A1,A2
-
-      RMAX=FIELD1(1,1,1)
-      RMIN=FIELD1(1,1,1)
-      
-      DO 100 K=1,NK
-      DO 100 J=1,NJ
-        DO 100 I=1,NI
-          IF (FIELD1(I,J,K).GT.RMAX)RMAX=FIELD1(I,J,K)
-          IF (FIELD1(I,J,K).LT.RMIN)RMIN=FIELD1(I,J,K)
-100   CONTINUE
-
-      IF (ABS(RMIN).GT.RMAX.OR.ABS(RMIN).EQ.RMAX) THEN
-      XMAX=ABS(RMIN)
-      ELSE
-      XMAX=RMAX
-      ENDIF
-      
-      IF (XMAX.EQ.0) THEN
-      IA = 0
-      RETURN
-      ENDIF
-      
-      A1=LOG10 ((G/10.d0)/XMAX)
-      A2=LOG10 ( G/XMAX )
-      IF(A1 .gt. A2) THEN
-        A=A2
-      ELSE 
-        A=A1
-      ENDIF
-      
-      IF (A.GT.0) IA=INT(A)
-      IF (A.LT.0) IA=INT(A-1.0)
-      
-      RETURN
-      END
-      
-      SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG)
-	IMPLICIT REAL (A-H,O-Z)
-
-      REAL PHI(NI,NJ,NK),SIG,MW,RMS,P
- 
-      N=NI*NJ*NK
- 
-      RMS=0.
-	MW=0.
- 
-      DO 10 I=1,NI
-      DO 10 J=1,NJ
-        DO 10 K=1,NK
-	      P=PHI(I,J,K)
-          RMS=RMS+P*P
- 		MW=MW+P
-10    CONTINUE
-
-      RMS=SQRT(RMS/N)
-	MW=MW/N
-	
-	IF(RMS*RMS-MW*MW.LT.0.) THEN	  
-	  SIG=0.0
-	ELSE
-	  SIG=SQRT(RMS*RMS-MW*MW)
-	ENDIF
- 
-      RETURN
-      END
-      
Index: urce/Fortran/rwGRIB2.f90
===================================================================
--- Source/Fortran/rwGRIB2.f90	(revision None)
+++ 	(revision )
@@ -1,265 +1,0 @@
-      MODULE RWGRIB2
-
-      CONTAINS
-
-      SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
-
-      USE GRIB_API
-
-      IMPLICIT NONE
-
-        integer                            ::  ifile
-        integer                            ::  iret
-   integer                            ::  n,mk,parid,nm
-   integer                            ::  i,k
-   integer,dimension(:),allocatable   ::  igrib
-   integer                            ::  numberOfPointsAlongAParallel
-   integer                            ::  numberOfPointsAlongAMeridian
-   real, dimension(:), allocatable    ::  values
-   integer                            ::  numberOfValues
-   real,dimension(maxl,maxb,mlevel)   ::  feld  
-   integer:: maxl,maxb,mlevel,mstride,mpar(:),irest,div,level
-   integer :: l(size(mpar))
-   character*(*):: filename   
-   
-   feld=0.                          
-
-   call grib_open_file(ifile, TRIM(FILENAME),'r')
- 
-   ! count the messages in the file
-   call grib_count_in_file(ifile,n)
-   allocate(igrib(n))
-   igrib=-1
- 
-   ! Load the messages from the file.
-   DO i=1,n
-      call grib_new_from_file(ifile,igrib(i), iret)
-   END DO
- 
-   ! we can close the file
-   call grib_close_file(ifile)
- 
-      nm=size(mpar)
-      div=mlevel/nm
-      l=0
-      
-   ! Loop on all the messages in memory
-  iloop:  DO i=1,n
-!      write(*,*) 'processing message number ',i
-      !     get as a integer
-      call grib_get(igrib(i),'numberOfPointsAlongAParallel', &
-           numberOfPointsAlongAParallel)
- 
-      !     get as a integer
-      call grib_get(igrib(i),'numberOfPointsAlongAMeridian', &
-           numberOfPointsAlongAMeridian)
-
-      call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
-
-      call grib_get_size(igrib(i),'values',numberOfValues)
-!      write(*,*) 'numberOfValues=',numberOfValues
- 
-      allocate(values(numberOfValues), stat=iret)
-      !     get data values
-      call grib_get(igrib(i),'values',values)
-
-      call grib_get(igrib(i),'paramId',parid)
-      call grib_get(igrib(i),'level',level)
-
- kloop:  do k=1,nm
-        if(parid .eq. mpar(k)) then
-!         l(k)=l(k)+1
-         feld(:,:,(k-1)*div+level)=reshape(values,(/maxl,maxb/))
-!         print*,(k-1)*div+l(k),parid
-         exit kloop
-        endif
-      enddo kloop
-      if(k .gt. nm .and. parid .ne. mpar(nm)) then
-        write(*,*) k,nm,parid,mpar(nm)
-        write(*,*) 'ERROR readlatlon: parameter ',parid,'is not',mpar
-        stop
-      endif
-
-!      print*,i
-   END DO iloop
-   write(*,*) 'readlatlon: ',i-1,' records read'
- 
-   DO i=1,n
-     call grib_release(igrib(i))
-   END DO
- 
-   deallocate(values)
-   deallocate(igrib)
-
-      END SUBROUTINE READLATLON
-
-      SUBROUTINE WRITELATLON(iunit,igrib,ogrib,FELD,MAXL,MAXB,MLEVEL,&
-      MLEVELIST,MSTRIDE,MPAR)
-
-      USE GRIB_API
-
-      IMPLICIT NONE
-
-      INTEGER IFIELD,MLEVEL,MNAUF,I,J,K,L,MSTRIDE,IERR,JOUT
-      INTEGER MPAR(MSTRIDE),MAXL,MAXB,LEVMIN,LEVMAX
-      INTEGER IUNIT,igrib,ogrib
-      REAL ZSEC4(MAXL*MAXB)
-      REAL    FELD(MAXL,MAXB,MLEVEL)
-      CHARACTER*(*) MLEVELIST
-      INTEGER ILEVEL(MLEVEL),MLINDEX(MLEVEL+1),LLEN
-
-    ! parse MLEVELIST
-   
-      LLEN=len(trim(MLEVELIST))
-      if(index(MLEVELIST,'to') .ne. 0 .or. index(MLEVELIST,'TO') .ne. 0) THEN
-        i=index(MLEVELIST,'/')
-        read(MLEVELIST(1:i-1),*) LEVMIN
-        i=index(MLEVELIST,'/',.true.)
-        read(MLEVELIST(i+1:LLEN),*) LEVMAX
-        l=0
-        do i=LEVMIN,LEVMAX
-          l=l+1
-          ILEVEL(l)=i
-        enddo
-      else
-        l=1
-	MLINDEX(1)=0
-        do i=1,LLEN
-          if(MLEVELIST(i:i) .eq. '/') THEN
-	    l=l+1
-	    MLINDEX(l)=i
-	  endif
-	enddo
-	MLINDEX(l+1)=LLEN+1
-	do i=1,l
-	  read(MLEVELIST(MLINDEX(i)+1:MLINDEX(i+1)-1),*) ILEVEL(i)
-	enddo
-      endif 
-
-      DO k=1,l
-        call grib_set(igrib,"level",ILEVEL(k))
-        DO j=1,MSTRIDE
-         call grib_set(igrib,"paramId",MPAR(j))
-!         if(MPAR(j) .eq. 87) then
-!           call grib_set(igrib,"shortName","etadot")
-!           call grib_set(igrib,"units","Pa,s**-1")
-!         endif
-!         if(MPAR(j) .eq. 77) then
-!           call grib_set(igrib,"shortName","etadot")
-!           call grib_set(igrib,"units","s**-1")
-!         endif
-	 if(l .ne. mlevel) then
-           zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,ILEVEL(k)),(/maxl*maxb/))
-	 else
-	   zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,k),(/maxl*maxb/))
-	 endif
-         call grib_set(igrib,"values",zsec4)
-
-         call grib_write(igrib,iunit)
-
-        ENDDO
-      ENDDO
-
-
-
-      END SUBROUTINE WRITELATLON
-
-      SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,&
-        MAXLEV,MPAR,A,B)
-
-      USE GRIB_API
-
-      IMPLICIT NONE
-
-
-        integer                            ::  ifile
-   integer                            ::  iret
-   integer                            ::  n,mk,div,nm,k
-   integer                            ::  i,j,parid
-   integer,dimension(:),allocatable   ::  igrib
-   real, dimension(:), allocatable    ::  values
-   integer                            ::  numberOfValues,maxlev
-   REAL :: A(MAXLEV+1),B(MAXLEV+1),pv(2*MAXLEV+2)
-   REAL:: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
-integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
-character*(*):: filename                             
- 
-   call grib_open_file(ifile, TRIM(FILENAME),'r')
- 
-   ! count the messages in the file
-   call grib_count_in_file(ifile,n)
-   allocate(igrib(n))
-   igrib=-1
- 
-   ! Load the messages from the file.
-   DO i=1,n
-      call grib_new_from_file(ifile,igrib(i), iret)
-   END DO
- 
-   ! we can close the file
-   call grib_close_file(ifile)
- 
-    l=0
-   ! Loop on all the messages in memory
-   iloop: DO i=1,n
-   ! write(*,*) 'processing message number ',i
-      !     get as a integer
-      call grib_get(igrib(i),'pentagonalResolutionParameterJ', j)
-
-      call grib_get_size(igrib(i),'values',numberOfValues)
-   !   write(*,*) 'numberOfValues=',numberOfValues
- 
-      call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
-
-      call grib_get(igrib(i),'level',ilev)
-
-      
-
-      call grib_get(igrib(i),'pv',pv)
-
-      allocate(values(numberOfValues), stat=iret)
-      !     get data values
-      call grib_get(igrib(i),'values',values)
-
-!      IOFFSET=mod(i-1,MSTRIDE)*(mk/2-1)
-!           CXMN(:,IOFFSET+ilev)=values(1:(MNAUF+1)*(MNAUF+2))
-
-      call grib_get(igrib(i),'paramId',parid)
-    nm=size(mpar)
-    div=mlevel/nm
-    kloop:  do k=1,nm
-        if(parid .eq. mpar(k)) then
-         l(k)=l(k)+1
-         cxmn(:,(k-1)*div+l(k))=values(1:(MNAUF+1)*(MNAUF+2))
-!         print*,(k-1)*div+l(k),parid
-         exit kloop
-        endif
-        
-      enddo kloop
-      if(k .gt. nm .and. parid .ne. mpar(nm)) then
-        write(*,*) k,nm,parid,mpar(nm)
-        write(*,*) 'ERROR readspectral: parameter ',parid,'is not',mpar
-        stop
-      endif
-
-!      print*,i
-
-   END DO iloop
-
-   write(*,*) 'readspectral: ',i-1,' records read'
- 
-   DO i=1,n
-     call grib_release(igrib(i))
-   END DO
- 
-   deallocate(values)
-   deallocate(igrib)
-
-
-
-        A=pv(1:1+MAXLEV)
-        B=pv(2+MAXLEV:2*MAXLEV+2)
-
-      END SUBROUTINE READSPECTRAL
-
-      END MODULE RWGRIB2
Index: Source/Fortran/rwgrib2.f90
===================================================================
--- Source/Fortran/rwgrib2.f90	(revision None)
+++ Source/Fortran/rwgrib2.f90	(revision None)
@@ -0,0 +1,265 @@
+ MODULE RWGRIB2
+
+ CONTAINS
+
+ SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
+
+!! Read a field from GRIB file on lat-lon grid
+
+ USE GRIB_API
+
+ IMPLICIT NONE
+
+   integer                            ::  ifile
+   integer                            ::  iret
+   integer                            ::  n,mk,parid,nm
+   integer                            ::  i,k
+   integer,dimension(:),allocatable   ::  igrib
+   integer                            ::  numberOfPointsAlongAParallel
+   integer                            ::  numberOfPointsAlongAMeridian
+   real, dimension(:), allocatable    ::  values
+   integer                            ::  numberOfValues
+   real,dimension(maxl,maxb,mlevel)   ::  feld  
+   integer:: maxl,maxb,mlevel,mstride,mpar(:),irest,div,level
+   integer :: l(size(mpar))
+   character*(*):: filename   
+   
+   feld=0.                          
+
+   call grib_open_file(ifile, TRIM(FILENAME),'r')
+ 
+   ! count the messages in the file
+   call grib_count_in_file(ifile,n)
+   allocate(igrib(n))
+   igrib=-1
+ 
+   ! Load the messages from the file.
+   DO i=1,n
+      call grib_new_from_file(ifile,igrib(i), iret)
+   END DO
+ 
+   ! we can close the file
+   call grib_close_file(ifile)
+ 
+      nm=size(mpar)
+      div=mlevel/nm
+      l=0
+      
+   ! Loop on all the messages in memory
+  iloop:  DO i=1,n
+!      write(*,*) 'processing message number ',i
+      !     get as a integer
+      call grib_get(igrib(i),'numberOfPointsAlongAParallel', &
+           numberOfPointsAlongAParallel)
+ 
+      !     get as a integer
+      call grib_get(igrib(i),'numberOfPointsAlongAMeridian', &
+           numberOfPointsAlongAMeridian)
+
+      call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
+
+      call grib_get_size(igrib(i),'values',numberOfValues)
+!      write(*,*) 'numberOfValues=',numberOfValues
+ 
+      allocate(values(numberOfValues), stat=iret)
+      !     get data values
+      call grib_get(igrib(i),'values',values)
+
+      call grib_get(igrib(i),'paramId',parid)
+      call grib_get(igrib(i),'level',level)
+
+ kloop:  do k=1,nm
+        if (parid .eq. mpar(k)) then
+!         l(k)=l(k)+1
+         feld(:,:,(k-1)*div+level)=reshape(values,(/maxl,maxb/))
+!         print*,(k-1)*div+l(k),parid
+         exit kloop
+        end if
+      end do kloop
+      if (k .gt. nm .and. parid .ne. mpar(nm)) then
+        write(*,*) k,nm,parid,mpar(nm)
+        write(*,*) 'ERROR readlatlon: parameter ',parid,'is not',mpar
+        stop
+      end if
+
+!      print*,i
+   END DO iloop
+!!   write(*,*) 'readlatlon: ',i-1,' records read'
+ 
+   DO i=1,n
+     call grib_release(igrib(i))
+   END DO
+ 
+   if (allocated(values)) deallocate(values)
+   deallocate(igrib)
+
+   END SUBROUTINE READLATLON
+
+   SUBROUTINE WRITELATLON(iunit,igrib,ogrib,FELD,MAXL,MAXB,MLEVEL,&
+   MLEVELIST,MSTRIDE,MPAR)
+
+!! write a field on lat-lon grid to GRIB file
+
+   USE GRIB_API
+
+   IMPLICIT NONE
+
+   INTEGER IFIELD,MLEVEL,MNAUF,I,J,K,L,MSTRIDE,IERR,JOUT
+   INTEGER MPAR(MSTRIDE),MAXL,MAXB,LEVMIN,LEVMAX
+   INTEGER IUNIT,igrib,ogrib
+   REAL ZSEC4(MAXL*MAXB)
+   REAL    FELD(MAXL,MAXB,MLEVEL)
+   CHARACTER*(*) MLEVELIST
+   INTEGER ILEVEL(MLEVEL),MLINDEX(MLEVEL+1),LLEN
+
+ ! parse MLEVELIST
+
+   LLEN=len(trim(MLEVELIST))
+   if (index(MLEVELIST,'to') .ne. 0 .or. index(MLEVELIST,'TO') .ne. 0) THEN
+     i=index(MLEVELIST,'/')
+     read(MLEVELIST(1:i-1),*) LEVMIN
+     i=index(MLEVELIST,'/',.true.)
+     read(MLEVELIST(i+1:LLEN),*) LEVMAX
+     l=0
+     do i=LEVMIN,LEVMAX
+       l=l+1
+       ILEVEL(l)=i
+     end do
+   else
+     l=1
+     MLINDEX(1)=0
+     do i=1,LLEN
+       if (MLEVELIST(i:i) .eq. '/') THEN
+	        l=l+1
+	        MLINDEX(l)=i
+	      end if
+	    end do
+	    MLINDEX(l+1)=LLEN+1
+	    do i=1,l
+	      read(MLEVELIST(MLINDEX(i)+1:MLINDEX(i+1)-1),*) ILEVEL(i)
+	    end do
+   end if 
+
+   DO k=1,l
+     call grib_set(igrib,"level",ILEVEL(k))
+     DO j=1,MSTRIDE
+       call grib_set(igrib,"paramId",MPAR(j))
+!         if (MPAR(j) .eq. 87) then
+!           call grib_set(igrib,"shortName","etadot")
+!           call grib_set(igrib,"units","Pa,s**-1")
+!         end if
+!         if (MPAR(j) .eq. 77) then
+!           call grib_set(igrib,"shortName","etadot")
+!           call grib_set(igrib,"units","s**-1")
+!         end if
+	      if (l .ne. mlevel) then
+         zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,ILEVEL(k)),(/maxl*maxb/))
+	      else
+	        zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,k),(/maxl*maxb/))
+	      end if
+       call grib_set(igrib,"values",zsec4)
+
+       call grib_write(igrib,iunit)
+
+     END DO
+   END DO
+
+   END SUBROUTINE WRITELATLON
+
+   SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,MAXLEV,MPAR,A,B)
+
+!!  read a GRIB file in spherical harmonics
+
+   USE GRIB_API
+
+   IMPLICIT NONE
+
+   integer                            ::  ifile
+   integer                            ::  iret
+   integer                            ::  n,mk,div,nm,k
+   integer                            ::  i,j,parid
+   integer,dimension(:),allocatable   ::  igrib
+   real, dimension(:), allocatable    ::  values
+   integer                            ::  numberOfValues,maxlev
+   REAL :: A(MAXLEV+1),B(MAXLEV+1),pv(2*MAXLEV+2)
+   REAL:: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+  integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
+   character*(*):: filename                             
+ 
+   call grib_open_file(ifile, TRIM(FILENAME),'r')
+ 
+   ! count the messages in the file
+   call grib_count_in_file(ifile,n)
+   allocate(igrib(n))
+   igrib=-1
+ 
+   ! Load the messages from the file.
+   DO i=1,n
+      call grib_new_from_file(ifile,igrib(i), iret)
+   END DO
+ 
+   ! we can close the file
+   call grib_close_file(ifile)
+ 
+    l=0
+   ! Loop on all the messages in memory
+   iloop: DO i=1,n
+   ! write(*,*) 'processing message number ',i
+      !     get as a integer
+      call grib_get(igrib(i),'pentagonalResolutionParameterJ', j)
+
+      call grib_get_size(igrib(i),'values',numberOfValues)
+   !   write(*,*) 'numberOfValues=',numberOfValues
+ 
+      call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
+
+      call grib_get(igrib(i),'level',ilev)
+
+      
+
+      call grib_get(igrib(i),'pv',pv)
+
+      allocate(values(numberOfValues), stat=iret)
+      !     get data values
+      call grib_get(igrib(i),'values',values)
+
+!      IOFFSET=mod(i-1,MSTRIDE)*(mk/2-1)
+!           CXMN(:,IOFFSET+ilev)=values(1:(MNAUF+1)*(MNAUF+2))
+
+      call grib_get(igrib(i),'paramId',parid)
+      nm=size(mpar)
+      div=mlevel/nm
+      kloop:  do k=1,nm
+        if (parid .eq. mpar(k)) then
+         l(k)=l(k)+1
+         cxmn(:,(k-1)*div+l(k))=values(1:(MNAUF+1)*(MNAUF+2))
+!         print*,(k-1)*div+l(k),parid
+         exit kloop
+        end if
+        
+      end do kloop
+      if (k .gt. nm .and. parid .ne. mpar(nm)) then
+        write(*,*) k,nm,parid,mpar(nm)
+        write(*,*) 'ERROR readspectral: parameter ',parid,'is not',mpar
+        stop
+      end if
+
+!      print*,i
+
+   END DO iloop
+
+!!   write(*,*) 'readspectral: ',i-1,' records read'
+ 
+   DO i=1,n
+     call grib_release(igrib(i))
+   END DO
+ 
+   deallocate(values)
+   deallocate(igrib)
+
+   A=pv(1:1+MAXLEV)
+   B=pv(2+MAXLEV:2*MAXLEV+2)
+
+   END SUBROUTINE READSPECTRAL
+
+END MODULE RWGRIB2
Index: Testing/Regression/FortranEtadot/readme.txt
===================================================================
--- Testing/Regression/FortranEtadot/readme.txt	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
+++ Testing/Regression/FortranEtadot/readme.txt	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -0,0 +1,33 @@
+HOW TO DO REGRESSION TESTS OF THE FORTRAN CODE
+
+1. Go to flex_extract/Testing/Regression/FortranEtadot
+2. Download the tarball containing the input files and reference outputs
+2. Untar the tarball
+3. Create a wirking directory: mkdir Work
+4. Compile the unmodified Fortran code with makefile_fast and makefile_debug
+   (is in flex_extract/Source/Fortran)
+5. Run a regression test to see whether the current Fortran code gives
+   output consistent with the reference output.
+   If not, carefully check why (machine-dependent small deviation?)
+   The output from the regression run is in 'Outputs' (automatically created).
+   If you need a new reference, you could remove or rename 'Outputs_ref', 
+   and then run
+   ./mk_outputdirs.sh
+   ./run_ref.sh
+   to create a new reference version. 
+6. Work on the code and use the 'run_ref.sh' script to test your results.
+
+Note 1: The regression tests except those with "high" in their name will
+      altogether run in about 1 minute. The "high" tests (hemispherical data)
+      can take many minutes and also require up to ca. 20 GB of memory.
+      Therefore, the script can be invoked as
+      ./run_ref.sh omithigh
+      to omit the "high" tests. For single development steps this should be
+      sufficient. When you are satisfied, run the "high" tests at the end.
+Note 2: The test scripts contain
+      export OMP_NUM_THREADS=4 # you may want to change this
+      export OMP_PLACES=cores
+      You should set OMP_NUM_THREADS to the number of physical cores of your
+      test machine or less.
+      OMP environment variables are explained on
+      https://gcc.gnu.org/onlinedocs/libgomp/#toc-OpenMP-Environment-Variables
Index: Testing/Regression/FortranEtadot/run_regrtest.sh
===================================================================
--- Testing/Regression/FortranEtadot/run_regrtest.sh	(revision 6bc4b425ba942c4c9d48c0c37fa5667b00c457ad)
+++ Testing/Regression/FortranEtadot/run_regrtest.sh	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -9,5 +9,6 @@
 
 export OMP_NUM_THREADS=4 # you may want to change this
-export OMP_PLACES=sockets
+export OMP_PLACES=cores
+export OMP_DISPLAY_ENV=verbose
 testhome=`pwd`
 path1=../../../Source/Fortran/
Index: Testing/Regression/FortranEtadot/runtimes.csv
===================================================================
--- Testing/Regression/FortranEtadot/runtimes.csv	(revision 2d1c338814a949a88dd9f5b55ddec4efbd7b60e6)
+++ Testing/Regression/FortranEtadot/runtimes.csv	(revision dfa7dbd2f70a8aa071dbb9f8b4928984af159da3)
@@ -15,2 +15,27 @@
 0007a71, 'reference', 'latlonall_debug', 0m0.405s, 0m0.536s, 0m0.008s
 0007a71, 'reference', 'latlonall_fast', 0m0.350s, 0m0.496s, 0m0.004s
+6bc4b42, 'reference', 'etadot_debug', 0m0.322s, 0m0.360s, 0m0.012s
+6bc4b42, 'reference', 'etadot_fast', 0m0.311s, 0m0.348s, 0m0.008s
+6bc4b42, 'reference', 'etadotall_debug', 0m0.387s, 0m0.428s, 0m0.008s
+6bc4b42, 'reference', 'etadotall_fast', 0m0.382s, 0m0.424s, 0m0.004s
+6bc4b42, 'reference', 'gauss_debug', 0m3.346s, 0m8.776s, 0m0.152s
+6bc4b42, 'reference', 'gauss_fast', 0m2.060s, 0m4.036s, 0m0.108s
+6bc4b42, 'reference', 'gaussall_debug', 0m13.690s, 0m38.840s, 0m0.292s
+6bc4b42, 'reference', 'gaussall_fast', 0m7.507s, 0m16.524s, 0m0.308s
+6bc4b42, 'reference', 'latlon_debug', 0m0.387s, 0m0.420s, 0m0.016s
+6bc4b42, 'reference', 'latlon_fast', 0m0.335s, 0m0.376s, 0m0.004s
+6bc4b42, 'reference', 'latlonall_debug', 0m0.351s, 0m0.392s, 0m0.004s
+6bc4b42, 'reference', 'latlonall_fast', 0m0.338s, 0m0.380s, 0m0.004s
+8c55c02, 'reference', 'etadot_debug', 0m0.263s, 0m0.304s, 0m0.004s
+8c55c02, 'reference', 'etadot_fast', 0m0.261s, 0m0.296s, 0m0.008s
+8c55c02, 'reference', 'etadotall_debug', 0m0.348s, 0m0.388s, 0m0.004s
+8c55c02, 'reference', 'etadotall_fast', 0m0.391s, 0m0.428s, 0m0.008s
+8c55c02, 'reference', 'gauss_debug', 0m1.985s, 0m3.940s, 0m0.104s
+8c55c02, 'reference', 'gauss_fast', 0m1.976s, 0m3.912s, 0m0.108s
+8c55c02, 'reference', 'gaussall_debug', 0m7.653s, 0m16.852s, 0m0.368s
+8c55c02, 'reference', 'gaussall_fast', 0m7.533s, 0m16.456s, 0m0.408s
+8c55c02, 'reference', 'latlon_debug', 0m0.334s, 0m0.376s, 0m0.004s
+8c55c02, 'reference', 'latlon_fast', 0m0.332s, 0m0.368s, 0m0.008s
+8c55c02, 'reference', 'latlonall_debug', 0m0.394s, 0m0.424s, 0m0.016s
+8c55c02, 'reference', 'latlonall_fast', 0m0.344s, 0m0.376s, 0m0.012s
+
