diff --git a/Source/Fortran/calc_etadot.f90 b/Source/Fortran/calc_etadot.f90
new file mode 100644
index 0000000000000000000000000000000000000000..542997d20f47c0e4b7dec9968d30dacf70fd72f7
--- /dev/null
+++ b/Source/Fortran/calc_etadot.f90
@@ -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
diff --git a/Source/Fortran/calc_etadot.out b/Source/Fortran/calc_etadot.out
new file mode 120000
index 0000000000000000000000000000000000000000..799fa15cd10f62ddc15acfb5a6037f9da29d84c7
--- /dev/null
+++ b/Source/Fortran/calc_etadot.out
@@ -0,0 +1 @@
+calc_etadot_fast.out
\ No newline at end of file
diff --git a/Source/Fortran/changelog.txt b/Source/Fortran/changelog.txt
new file mode 100644
index 0000000000000000000000000000000000000000..18de918707c431461a2e93aadc56c0cb19be8e88
--- /dev/null
+++ b/Source/Fortran/changelog.txt
@@ -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
+
+
diff --git a/Source/Fortran/ftrafo.f b/Source/Fortran/ftrafo.f
deleted file mode 100644
index affdccdcf8b2a439fbd5fd35435ac076eefe7ee5..0000000000000000000000000000000000000000
--- a/Source/Fortran/ftrafo.f
+++ /dev/null
@@ -1,504 +0,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
diff --git a/Source/Fortran/ftrafo.f90 b/Source/Fortran/ftrafo.f90
new file mode 100644
index 0000000000000000000000000000000000000000..7995967ac40b0002d0858dd22494cb9b9fadd4ae
--- /dev/null
+++ b/Source/Fortran/ftrafo.f90
@@ -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
diff --git a/Source/Fortran/grphreal.f b/Source/Fortran/grphreal.f
deleted file mode 100644
index dae342bf336d149d97ffa55234c0b1375eff2661..0000000000000000000000000000000000000000
--- a/Source/Fortran/grphreal.f
+++ /dev/null
@@ -1,188 +0,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
diff --git a/Source/Fortran/grphreal.f90 b/Source/Fortran/grphreal.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b9e61f20736dfffe4c8db64bb65b74581c59f9f7
--- /dev/null
+++ b/Source/Fortran/grphreal.f90
@@ -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
diff --git a/Source/Fortran/makefile.local.gfortran b/Source/Fortran/makefile.local.gfortran
deleted file mode 100644
index 0e03b8d7c1335cce4f9a8baa1035c904d85df176..0000000000000000000000000000000000000000
--- a/Source/Fortran/makefile.local.gfortran
+++ /dev/null
@@ -1,63 +0,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
-#
-###############################################################################
diff --git a/Source/Fortran/makefile.CRAY b/Source/Fortran/makefile_cray
similarity index 88%
rename from Source/Fortran/makefile.CRAY
rename to Source/Fortran/makefile_cray
index cf851930dff3ad62abf6f30f16e3926d43b208b2..329fad55612498c574aa964e988acfc3bcafedc5 100644
--- a/Source/Fortran/makefile.CRAY
+++ b/Source/Fortran/makefile_cray
@@ -24,7 +24,7 @@ LDFLAGS = $(OPT)
BINDIR = .
-EXE = CONVERT2
+EXE = calc_etadot
.f.o:
@@ -46,8 +46,8 @@ grphreal.o: 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}
+$(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}
###############################################################################
diff --git a/Source/Fortran/makefile_debug b/Source/Fortran/makefile_debug
index 187234eaf5f7e1bd3769f9c2030a260e757f3623..e87cf2f190ebb1787184dc7c53bc1434e8eee8f3 100644
--- a/Source/Fortran/makefile_debug
+++ b/Source/Fortran/makefile_debug
@@ -1,62 +1,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
+# 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 = -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
-
-$(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}
-
+GRIB_API_LIB= -Bstatic -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper
+EMOSLIB=-lemosR64
+LIB = $(GRIB_API_LIB) $(EMOSLIB)
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
+GRIB_API_INCLUDE_DIR=/usr/include
+INC = -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
+
+.DEFAULT:
+ -echo $@ does not exist.
+
+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
+
+clean: -rm -f $(OBJ) ${EXE} $(MOD)
+
+${EXE}: $(OBJ)
+ $(FC) $(OBJ) -o ${EXE} $(LDFLAGS)
diff --git a/Source/Fortran/makefile_debug2 b/Source/Fortran/makefile_debug2
deleted file mode 100644
index c950176147a957f046cbec65881ecb0f03690609..0000000000000000000000000000000000000000
--- a/Source/Fortran/makefile_debug2
+++ /dev/null
@@ -1,66 +0,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
-#
-###############################################################################
-
diff --git a/Source/Fortran/makefile.ecgate b/Source/Fortran/makefile_ecgate
similarity index 89%
rename from Source/Fortran/makefile.ecgate
rename to Source/Fortran/makefile_ecgate
index 0c46c5b3c9e0255339c9a0d28e83648997ceab0e..dace8ad8983469b8549da1c07665dd3eee9cfb7c 100644
--- a/Source/Fortran/makefile.ecgate
+++ b/Source/Fortran/makefile_ecgate
@@ -25,7 +25,7 @@ LDFLAGS = $(OPT)
BINDIR = .
-EXE = CONVERT2
+EXE = calc_etadot
.f.o:
@@ -47,8 +47,8 @@ grphreal.o: 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}
+$(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}
###############################################################################
diff --git a/Source/Fortran/makefile_fast b/Source/Fortran/makefile_fast
index c5a4d75258e3c52fdbdfc5ad96db9b288a9928ff..fbb60d9d210b754a366f729ae593412733abbc09 100644
--- a/Source/Fortran/makefile_fast
+++ b/Source/Fortran/makefile_fast
@@ -1,62 +1,53 @@
###############################################################################
#
# 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
-
-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}
-
+GRIB_API_LIB= -Bstatic -lgrib_api_f90 -lgrib_api -Bdynamic -lm -ljasper
+EMOSLIB=-lemosR64
+LIB = $(GRIB_API_LIB) $(EMOSLIB)
-###############################################################################
-#
-# End of the Makefile
-#
-###############################################################################
+GRIB_API_INCLUDE_DIR=/usr/include
+INC = -I. -I$(GRIB_API_INCLUDE_DIR)
+
+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
+
+.DEFAULT:
+ -echo $@ does not exist.
+
+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
+
+clean:
+ -rm -f $(OBJ) ${EXE} $(MOD)
+
+${EXE}: $(OBJ)
+ $(FC) $(OBJ) -o ${EXE} $(LDFLAGS)
diff --git a/Source/Fortran/makefile_local_gfortran b/Source/Fortran/makefile_local_gfortran
new file mode 120000
index 0000000000000000000000000000000000000000..052b83a4d088620dabe01455d0094404382e9acd
--- /dev/null
+++ b/Source/Fortran/makefile_local_gfortran
@@ -0,0 +1 @@
+makefile_fast
\ No newline at end of file
diff --git a/Source/Fortran/makefile.local.ifort b/Source/Fortran/makefile_local_ifort
similarity index 89%
rename from Source/Fortran/makefile.local.ifort
rename to Source/Fortran/makefile_local_ifort
index 2744a9e38af4069a9dc4f9bc3fad90abda4e4ed3..310dc9b728d9621b44ca99a5f687f7c59741e3a2 100644
--- a/Source/Fortran/makefile.local.ifort
+++ b/Source/Fortran/makefile_local_ifort
@@ -28,7 +28,7 @@ LDFLAGS = $(OPT)
BINDIR = .
-EXE = CONVERT2.s8.ifort
+EXE = calc_etadot.s8.ifort
.f.o:
@@ -50,8 +50,8 @@ grphreal.o: 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}
+$(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}
###############################################################################
diff --git a/Source/Fortran/phgrreal.f b/Source/Fortran/phgrreal.f
deleted file mode 100644
index aa3658c917fd35d15dd4c2b9959f4fcf074923a1..0000000000000000000000000000000000000000
--- a/Source/Fortran/phgrreal.f
+++ /dev/null
@@ -1,553 +0,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
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/Source/Fortran/phgrreal.f90 b/Source/Fortran/phgrreal.f90
new file mode 100644
index 0000000000000000000000000000000000000000..2cf9ca38483dae75df047a885be7aa53b4dbd976
--- /dev/null
+++ b/Source/Fortran/phgrreal.f90
@@ -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
diff --git a/Source/Fortran/posnam.f b/Source/Fortran/posnam.f
deleted file mode 100644
index c5d12d2b9928e581f67ef0c6388dd3e641693aed..0000000000000000000000000000000000000000
--- a/Source/Fortran/posnam.f
+++ /dev/null
@@ -1,25 +0,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
-
diff --git a/Source/Fortran/posnam.f90 b/Source/Fortran/posnam.f90
new file mode 100644
index 0000000000000000000000000000000000000000..614bf7f843c6be4c6ae599d8e14a0b5d630174cd
--- /dev/null
+++ b/Source/Fortran/posnam.f90
@@ -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
diff --git a/Source/Fortran/preconvert.f90 b/Source/Fortran/preconvert.f90
deleted file mode 100644
index c28610f2aecb31a3f4c3860d7a482ae0fa85d23f..0000000000000000000000000000000000000000
--- a/Source/Fortran/preconvert.f90
+++ /dev/null
@@ -1,807 +0,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
-
diff --git a/Source/Fortran/rwGRIB2.f90 b/Source/Fortran/rwgrib2.f90
similarity index 60%
rename from Source/Fortran/rwGRIB2.f90
rename to Source/Fortran/rwgrib2.f90
index 734d06ca185685df6a584d4c91fbe617eddc26b6..240547ee7222f83796283c21886e056a856af6e1 100644
--- a/Source/Fortran/rwGRIB2.f90
+++ b/Source/Fortran/rwgrib2.f90
@@ -1,15 +1,17 @@
- MODULE RWGRIB2
+ MODULE RWGRIB2
- CONTAINS
+ CONTAINS
- SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
+ SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
- USE GRIB_API
+!! Read a field from GRIB file on lat-lon grid
- IMPLICIT NONE
+ USE GRIB_API
- integer :: ifile
- integer :: iret
+ IMPLICIT NONE
+
+ integer :: ifile
+ integer :: iret
integer :: n,mk,parid,nm
integer :: i,k
integer,dimension(:),allocatable :: igrib
@@ -67,112 +69,112 @@
call grib_get(igrib(i),'level',level)
kloop: do k=1,nm
- if(parid .eq. mpar(k)) then
+ 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
+ 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
- endif
+ end if
! print*,i
END DO iloop
- write(*,*) 'readlatlon: ',i-1,' records read'
+!! write(*,*) 'readlatlon: ',i-1,' records read'
DO i=1,n
call grib_release(igrib(i))
END DO
- deallocate(values)
+ if (allocated(values)) 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
+ 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")
-! endif
-! if(MPAR(j) .eq. 77) then
+! end if
+! 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 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
+ END SUBROUTINE WRITELATLON
- SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,&
- MAXLEV,MPAR,A,B)
+ SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,MAXLEV,MPAR,A,B)
- USE GRIB_API
+!! read a GRIB file in spherical harmonics
- IMPLICIT NONE
+ USE GRIB_API
+ IMPLICIT NONE
- integer :: ifile
+ integer :: ifile
integer :: iret
integer :: n,mk,div,nm,k
integer :: i,j,parid
@@ -181,8 +183,8 @@
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
+ integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
+ character*(*):: filename
call grib_open_file(ifile, TRIM(FILENAME),'r')
@@ -225,28 +227,28 @@ character*(*):: filename
! 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
+ 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
+ end if
- enddo kloop
- if(k .gt. nm .and. parid .ne. mpar(nm)) then
+ 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
- endif
+ end if
! print*,i
END DO iloop
- write(*,*) 'readspectral: ',i-1,' records read'
+!! write(*,*) 'readspectral: ',i-1,' records read'
DO i=1,n
call grib_release(igrib(i))
@@ -255,11 +257,9 @@ character*(*):: filename
deallocate(values)
deallocate(igrib)
+ A=pv(1:1+MAXLEV)
+ B=pv(2+MAXLEV:2*MAXLEV+2)
+ END SUBROUTINE READSPECTRAL
- A=pv(1:1+MAXLEV)
- B=pv(2+MAXLEV:2*MAXLEV+2)
-
- END SUBROUTINE READSPECTRAL
-
- END MODULE RWGRIB2
+END MODULE RWGRIB2
diff --git a/Testing/Regression/FortranEtadot/readme.txt b/Testing/Regression/FortranEtadot/readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bf4f0331f5db545cc096cf36139c5f79ec28acc0
--- /dev/null
+++ b/Testing/Regression/FortranEtadot/readme.txt
@@ -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
diff --git a/Testing/Regression/FortranEtadot/run_regrtest.sh b/Testing/Regression/FortranEtadot/run_regrtest.sh
index 3d2fac853b7d9190d08006d52084d22c5ed2c668..e00a47e26bc0f6541e6347eafef17af49a30b6f7 100755
--- a/Testing/Regression/FortranEtadot/run_regrtest.sh
+++ b/Testing/Regression/FortranEtadot/run_regrtest.sh
@@ -8,7 +8,8 @@
# SPDX-License-Identifier: MIT-0
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/
path=../${path1}
diff --git a/Testing/Regression/FortranEtadot/runtimes.csv b/Testing/Regression/FortranEtadot/runtimes.csv
index 811aeb2dd95c374df3b8d19201408edea5385488..f3f6d5a3298c7ee487468baa9a9b14525a63104c 100644
--- a/Testing/Regression/FortranEtadot/runtimes.csv
+++ b/Testing/Regression/FortranEtadot/runtimes.csv
@@ -14,3 +14,28 @@
0007a71, 'reference', 'latlon_fast', 0m0.369s, 0m0.568s, 0m0.012s
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
+