source: flex_extract.git/Source/Fortran/phgrreal.f90 @ dfa7dbd

ctbtodev
Last change on this file since dfa7dbd was dfa7dbd, checked in by pesei <petra seibert -@…>, 5 years ago

changes in the Fortran part and associated regression test

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

all .f90 free format
code beautification
regression test is OK

make new local gfortran makefiles, remove parameters not needed
anymore

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

provide softlinks for standards:

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

provide changelog.txt in Fortran
provide readme.txt in Regression/FortranEtadot?
  • Property mode set to 100644
File size: 12.4 KB
Line 
1MODULE PHTOGR
2
3  INTEGER, PARAMETER :: MAXAUF=36000
4
5CONTAINS
6
7  SUBROUTINE PHGR213(CXMN,FELD,WSAVE,IFAX,Z,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
8
9!! DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
10!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
11!! RAUM AUF DAS REDUZIERTE GAUSS'SCHE GITTER DURCH
12!
13! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
14!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
15! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
16!       WSAVE  = Working Array fuer Fouriertransformation
17! Z         = LEGENDREFUNKTIONSWERTE
18!
19! MNAUF  ANZAHL DER FOURIERKOEFFIZIENTEN
20! MAXL   ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
21! MAXB   ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
22! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
23
24    IMPLICIT NONE
25
26!                       Anzahl der Gitterpunkte auf jedem Breitenkreis
27    INTEGER MLAT(MAXB/2)
28    INTEGER K,MAXL,MAXB,MLEVEL,MNAUF
29    INTEGER IND(MAXB)
30
31!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
32    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
33
34    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
35    REAL FELD(MAXL,MLEVEL)
36    REAL WSAVE(8*MAXB+15,MAXB/2)
37    INTEGER :: IFAX(10,MAXB)
38
39    IND(1)=0
40    DO 7 K=2,MAXB/2
41      IND(K)=IND(K-1)+MLAT(K-1)
427   CONTINUE
43
44!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
45    DO 17 K=1,MAXB/2
46      CALL PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
4717  CONTINUE
48!$OMP END PARALLEL DO
49
50    RETURN
51
52  END SUBROUTINE PHGR213
53
54  SUBROUTINE PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
55
56    IMPLICIT NONE
57
58    INTEGER MLAT(MAXB/2)
59    INTEGER K,L,I,J,LLS,LLPS,LL,LLP,MAXL,MAXB,MLEVEL,MNAUF
60    INTEGER IND(MAXB)
61    INTEGER :: IFAX(10,MAXB)
62
63
64!   FELD DER FOURIERKOEFFIZIENTEN
65    REAL :: CXMS(0:MAXAUF-1),CXMA(0:MAXAUF-1)
66
67!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
68    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
69    REAL ACR,ACI,SCR,SCI
70
71    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
72    REAL FELD(MAXL,MLEVEL)
73    REAL WSAVE(8*MAXB+15,MAXB/2)
74
75    DO 6 L=1,MLEVEL
76      LL=0
77      LLP=0
78      DO 1 I=0,MNAUF
79        SCR=0.D0
80        SCI=0.D0
81        ACR=0.D0
82        ACI=0.D0
83        LLS=LL
84        LLPS=LLP
85        IF (2*I+1 .LT. MLAT(K)) THEN
86!               Innerste Schleife aufgespalten um if-Abfrage zu sparen
87          DO 18 J=I,MNAUF,2
88            SCR=SCR+Z(LLP,K)*CXMN(2*LL,L)
89            SCI=SCI+Z(LLP,K)*CXMN(2*LL+1,L)
90            LL=LL+2
91            LLP=LLP+2
9218        CONTINUE
93          LL=LLS+1
94          LLP=LLPS+1
95          DO 19 J=I+1,MNAUF,2
96            ACR=ACR+Z(LLP,K)*CXMN(2*LL,L)
97            ACI=ACI+Z(LLP,K)*CXMN(2*LL+1,L)
98            LL=LL+2
99            LLP=LLP+2
10019        CONTINUE
101        END IF
102        LL=LLS+(MNAUF-I+1)
103        LLP=LLPS+(MNAUF-I+3)
104        CXMS(2*I)=SCR+ACR
105        CXMS(2*I+1)=SCI+ACI
106        CXMA(2*I)=SCR-ACR
107        CXMA(2*I+1)=SCI-ACI
1081     CONTINUE
109!     CALL FOURTR(CXMS,FELD(IND(k)+1,L),WSAVE(:,K),MNAUF,*MLAT(K),1)
110!     CALL FOURTR(CXMA,FELD(MAXL-IND(k)-MLAT(K)+1,L),WSAVE(:,K),MNAUF,MLAT(K),1)
111      CALL RFOURTR(CXMS,WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K),1)
112      FELD(IND(k)+1:IND(K)+MLAT(K),L)=CXMS(0:MLAT(K)-1)
113      CALL RFOURTR(CXMA,WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K),1)
114      FELD(MAXL-IND(k)-MLAT(K)+1:MAXL-IND(k),L)=CXMA(0:MLAT(K)-1)
115!      WRITE(*,*) IND+1,FELD(IND+1,L)
1166   CONTINUE
117
118  END SUBROUTINE PHSYM
119
120  SUBROUTINE PHGCUT(CXMN,FELD,WSAVE,IFAX,Z, &
121    MNAUF,MMAX,MAUF,MANF,MAXL,MAXB,MLEVEL)
122
123!! DIE ROUTINE FUEHRT EINE TRANSFORMATION EINER
124!! FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
125!! RAUM AUF KUGELKOORDINATEN DURCH. Es kann ein Teilausschnitt
126!!      Der Erde angegeben werden. Diese Routine ist langsamer als phgrph
127
128! CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
129!          CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
130! FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
131! BREITE = SINUS DER GEOGRAFISCHEN BREITEN
132!
133! MNAUF  ANZAHL DER FOURIERKOEFFIZIENTEN
134! MAUF   ANZAHL DER LAENGEN UND DER FOURIERKOEFFIZIENTEN
135! MANF   ANFANG DES LAENGENBEREICHS FUER DAS GITTER,
136!        AUF DAS INTERPOLIERT WERDEN SOLL
137! MAXL   ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
138! MAXB   ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
139! MLEVEL ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
140
141    IMPLICIT REAL (A-H,O-Z)
142
143!   FELD DER FOURIERKOEFFIZIENTEN
144
145!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
146    REAL Z(0:((MMAX+3)*(MMAX+4))/2,MAXB)
147
148    DIMENSION CXMN(0:(MMAX+1)*(MMAX+2)-1,MLEVEL)
149    REAL FELD(MAXL,MAXB,MLEVEL)
150    DIMENSION WSAVE(4*MAUF+15)
151    INTEGER:: IFAX(10)
152
153    LOGICAL SYM
154
155!    write(*,*)mauf,mnauf,manf,maxl
156
157    IF (MAUF .LE. MNAUF) WRITE(*,*) 'TOO COARSE LONGITUDE RESOLUTION'
158    IF (MANF .LT. 1    .OR. MAXL .LT. 1 .OR. &
159        MANF .GT. MAUF .OR. MAXL .GT. MAUF) THEN
160      WRITE(*,*) 'WRONG LONGITUDE RANGE',MANF,MAXL
161      STOP
162    END IF
163
164! Pruefe, ob Ausgabegitter symmetrisch zum Aequator ist
165! Wenn ja soll Symmetrie der Legendrepolynome ausgenutzt werden
166    IF (MAXB .GT. 4) THEN
167      SYM=.TRUE.
168      DO 11 J=5,5
169        IF (ABS(ABS(Z(100,J))-ABS(Z(100,MAXB+1-J))) .GT. 1E-11) SYM=.FALSE.
170!             WRITE(*,*) ABS(Z(100,J)),ABS(Z(100,MAXB+1-J))
17111    CONTINUE
172!!      WRITE(*,*) 'Symmetrisch: ',SYM
173    ELSE
174      SYM=.FALSE.
175    END IF
176
177
178    IF (SYM) THEN
179
180!$OMP PARALLEL DO
181      DO J=1,(MAXB+1)/2
182        CALL PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
183      END DO
184!$OMP END PARALLEL DO
185
186    ELSE
187
188!$OMP PARALLEL DO
189      DO J=1,MAXB
190        CALL PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
191      END DO
192!$OMP END PARALLEL DO
193
194    END IF
195
196    RETURN
197
198  END SUBROUTINE PHGCUT
199
200  SUBROUTINE PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
201
202    IMPLICIT REAL (A-H,O-Z)
203
204!   FELD DER FOURIERKOEFFIZIENTEN
205
206    REAL :: CXM(0:MAXAUF-1),CXMA(0:MAXAUF-1)
207
208!   FELD DER LEGENDREPOLYNOME FUER EINE BREITE
209    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
210    REAL SCR,SCI,ACR,ACI
211
212    DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
213    REAL FELD(MAXL,MAXB,MLEVEL)
214    DIMENSION WSAVE(4*MAUF+15)
215    INTEGER :: IFAX(10)
216
217    DO 16 L=1,MLEVEL
218      LL=0
219      LLP=0
220      DO 17 I=0,MNAUF
221        SCR=0.D0
222        SCI=0.D0
223        ACR=0.D0
224        ACI=0.D0
225        LLS=LL
226        LLPS=LLP
227!             Innerste Schleife aufgespalten um if-Abfrage zu sparen
228        DO 18 K=I,MNAUF,2
229          SCR=SCR+Z(LLP,J)*CXMN(2*LL,L)
230          SCI=SCI+Z(LLP,J)*CXMN(2*LL+1,L)
231          LL=LL+2
232          LLP=LLP+2
23318      CONTINUE
234        LL=LLS+1
235        LLP=LLPS+1
236        DO 19 K=I+1,MNAUF,2
237          ACR=ACR+Z(LLP,J)*CXMN(2*LL,L)
238          ACI=ACI +Z(LLP,J)*CXMN(2*LL+1,L)
239          LL=LL+2
240          LLP=LLP+2
24119      CONTINUE
242        LL=LLS+MNAUF-I+1
243        LLP=LLPS+MNAUF-I+3
244        CXM(2*I)=SCR+ACR
245        CXM(2*I+1)=SCI+ACI
246        CXMA(2*I)=SCR-ACR
247        CXMA(2*I+1)=SCI-ACI
24817    CONTINUE
249
250      CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
251      DO 26 I=0,MAXL-1
252        IF (MANF+I .LE. MAUF) THEN
253          FELD(I+1,J,L)=CXM(MANF+I-1)
254        ELSE
255          FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
256        END IF
25726    CONTINUE
258      CALL RFOURTR(CXMA,WSAVE,IFAX,MNAUF,MAUF,1)
259      DO 36 I=0,MAXL-1
260        IF (MANF+I .LE. MAUF) THEN
261          FELD(I+1,MAXB+1-J,L)=CXMA(MANF+I-1)
262        ELSE
263          FELD(I+1,MAXB+1-J,L)=CXMA(MANF-MAUF+I-1)
264        END IF
26536    CONTINUE
26616  CONTINUE
267
268  END SUBROUTINE PHSYMCUT
269
270  SUBROUTINE PHGPNS(CXMN,FELD,Z,WSAVE,IFAX,J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
271
272    IMPLICIT NONE
273   
274    INTEGER,INTENT(IN) :: MNAUF,MAUF,MANF,J,MAXL,MAXB,MLEVEL
275
276    REAL :: CXM(0:MAXAUF-1)
277    REAL,INTENT(IN) :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
278    REAL,INTENT(IN) :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
279    REAL,INTENT(IN) :: WSAVE(4*MAUF+15)
280    REAL :: FELD(MAXL,MAXB,MLEVEL)
281
282    INTEGER :: IFAX(10)
283    INTEGER I,L
284
285    DO L=1,MLEVEL
286      CALL LEGTR(CXMN(:,L),CXM,Z(:,J),MNAUF,MAUF)
287      CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF,1)
288
289      DO I=0,MAXL-1
290        IF (MANF+I .LE. MAUF) THEN
291          FELD(I+1,J,L)=CXM(MANF+I-1)
292        ELSE
293          FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
294        END IF
295      END DO
296    END DO
297  END SUBROUTINE PHGPNS
298
299  SUBROUTINE LEGTR(CXMN,CXM,Z,MNAUF,MAUF)
300
301!!   DIESE ROUTINE BERECHNET DIE FOURIERKOEFFIZIENTEN CXM
302
303
304    IMPLICIT NONE
305
306    INTEGER MNAUF,MAUF,LL,LLP,I,J
307    REAL CXM(0:MAXAUF-1)
308    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1)
309    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
310    REAL CI,CR
311
312    LL=0
313    LLP=0
314    DO 1 I=0,MNAUF
315      CR=0.D0
316      CI=0.D0
317      DO 2 J=I,MNAUF
318        CR=CR+Z(LLP)*CXMN(2*LL)
319        CI=CI+Z(LLP)*CXMN(2*LL+1)
320        LL=LL+1
321        LLP=LLP+1
3222     CONTINUE
323      LLP=LLP+2
324      CXM(2*I)=CR
325      CXM(2*I+1)=CI
3261   CONTINUE
327    RETURN
328   
329  END SUBROUTINE LEGTR
330
331  SUBROUTINE RFOURTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
332
333!!     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
334
335    IMPLICIT REAL (A-H,O-Z)
336
337    DIMENSION CXM(0:MAXAUF-1)
338    REAL :: WSAVE(2*MAXL),TRIGS(2*MAXL)
339    INTEGER IFAX(10)
340
341    DO I=MNAUF+1,MAXL-1
342      CXM(2*I)=0.0
343      CXM(2*I+1)=0.0
344    END DO
345
346    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,1)
347
348    DO I=0,MAXL-1
349      CXM(I)=CXM(I+1)
350    END DO
351
352    RETURN
353
354  END SUBROUTINE RFOURTR
355
356  SUBROUTINE GAULEG(X1,X2,X,W,N)
357
358!! BERECHNET DIE GAUSS+SCHEN BREITEN
359
360    IMPLICIT REAL (A-H,O-Z)
361
362    DIMENSION X(N),W(N)
363    PARAMETER (EPS=3.D-14)
364
365    M=(N+1)/2
366    XM=0.5D0*(X2+X1)
367    XL=0.5D0*(X2-X1)
368    DO 12 I=1,M
369      Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
3701     CONTINUE
371      P1=1.D0
372      P2=0.D0
373      DO 11 J=1,N
374        P3=P2
375        P2=P1
376        P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
37711    CONTINUE
378      PP=N*(Z*P1-P2)/(Z*Z-1.D0)
379      Z1=Z
380      Z=Z1-P1/PP
381      IF (ABS(Z-Z1) .GT. EPS)GO TO 1
382      X(I)=XM-XL*Z
383      X(N+1-I)=XM+XL*Z
384      W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
385      W(N+1-I)=W(I)
38612  CONTINUE
387
388    RETURN
389
390  END SUBROUTINE GAULEG
391
392
393  SUBROUTINE PLGNFA(LL,X,Z)
394
395!! PLGNDN BERECHNET ALLE NORMIERTEN ASSOZIIERTEN
396!! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
397!! UND SCHREIBT SIE IN DAS FELD Z
398! Die Polynome sind wie im ECMWF indiziert, d.h.
399! P00,P10,P11,P20,P21,P22,...
400!       Ansonsten ist die Routine analog zu PLGNDN
401! X IST DER COSINUS DES ZENITWINKELS ODER
402!       DER SINUS DER GEOGRAFISCHEN BREITE
403
404    IMPLICIT REAL (A-H,O-Z)
405
406    DIMENSION Z(0:((LL+3)*(LL+4))/2)
407
408    L=LL+2
409    I=1
410    Z(0)=1.D0
411    FACT=1.D0
412    POT=1.D0
413    SOMX2=DSQRT(1.D0-X*X)
414    DO 14 J=0,L
415      DJ=DBLE(J)
416      IF (J .GT. 0) THEN
417        FACT=FACT*(2.D0*DJ-1.D0)/(2.D0*DJ)
418        POT=POT*SOMX2
419        Z(I)=DSQRT((2.D0*DJ+1.D0)*FACT)*POT
420        I=I+1
421      END IF
422      IF (J .LT. L) THEN
423        Z(I)=X*DSQRT((4.D0*DJ*DJ+8.D0*DJ+3.D0)/(2.D0*DJ+1.D0))*Z(I-1)
424        I=I+1
425      END IF
426      DK=DJ+2.D0
427      DO 14 K=J+2,L
428        DDK=(DK*DK-DJ*DJ)
429        Z(I)=X*DSQRT((4.D0*DK*DK-1.D0)/DDK)*Z(I-1)- &
430          DSQRT(((2.D0*DK+1.D0)*(DK-DJ-1.D0)*(DK+DJ-1.D0))/ &
431          ((2.D0*DK-3.D0)*DDK))*Z(I-2)
432        DK=DK+1.D0
433        I=I+1
43414  CONTINUE
435
436    RETURN
437
438  END SUBROUTINE PLGNFA
439
440  SUBROUTINE DPLGND(MNAUF,Z,DZ)
441
442!! DPLGND BERECHNET DIE ABLEITUNG DER NORMIERTEN ASSOZIIERTEN
443!! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
444!! UND SCHREIBT SIE IN DAS FELD DZ
445! DIE REIHENFOLGE IST
446! P00(X),P01(X),P11(X),P02(X),P12(X),P22(X),..PLL(X)
447
448    IMPLICIT REAL (A-H,O-Z)
449
450    DIMENSION Z(0:((MNAUF+3)*(MNAUF+4))/2)
451    DIMENSION DZ(0:((MNAUF+2)*(MNAUF+3))/2)
452
453    IF (Z(0) .NE. 1.D0) THEN
454      WRITE(*,*) 'DPLGND: Z(0) must be 1.0'
455      STOP
456    END IF
457
458    LLP=0
459    LLH=0
460    DO 1 I=0,MNAUF+1
461      DO 2 J=I,MNAUF+1
462        IF (I .EQ. J) THEN
463          WURZELA=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
464          DZ(LLH)=DBLE(J)*WURZELA*Z(LLP+1)
465        ELSE
466          WURZELB=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
467          DZ(LLH)=DBLE(J)*WURZELB*Z(LLP+1)-DBLE(J+1)*WURZELA*Z(LLP-1)
468          WURZELA=WURZELB
469        END IF
470        LLH=LLH+1
471        LLP=LLP+1
4722     CONTINUE
473      LLP=LLP+1
4741   CONTINUE
475
476    RETURN
477
478  END SUBROUTINE DPLGND
479
480
481  SUBROUTINE SPFILTER(FELDMN,MM,MMAX)
482
483!! Spectral Filter of Sardeshmukh and Hoskins (1984, MWR)
484! MM=Spectral truncation of field
485! MMAX= Spectral truncation of filter
486
487    IMPLICIT NONE
488
489    INTEGER MM,MMAX,I,J,K,L
490    REAL FELDMN(0:(MM+1)*(MM+2)-1)
491    REAL KMAX,SMAX,FAK
492
493    SMAX=0.1
494    KMAX=-ALOG(SMAX)
495    KMAX=KMAX/(float(MMAX)*float(MMAX+1))**2
496!    WRITE(*,*)'alogsmax',alog(smax),'KMAX:',KMAX
497    L=0
498    DO I=0,MM
499      DO J=I,MM
500!        WRITE(*,*) I,J,FELD(K),FELD(K)*EXP(-KMAX*(J*(J+1))**2)
501        IF(J .LE. MMAX) THEN
502!          FAK=EXP(-KMAX*(J*(J+1))**2)
503          FAK=1.0
504          FELDMN(2*L)=FELDMN(2*L)*FAK
505          FELDMN(2*L+1)=FELDMN(2*L+1)*FAK
506        ELSE
507          FELDMN(2*L)=0.
508          FELDMN(2*L+1)=0.
509        END IF
510        L=L+1
511      END DO
512    END DO
513   
514  END SUBROUTINE SPFILTER
515
516END MODULE PHTOGR
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG