source: flex_extract.git/src/phgrreal.f @ f9f7b3f

ctbtodev
Last change on this file since f9f7b3f was d69b677, checked in by Anne Philipp <bscannephilipp@…>, 6 years ago

original ECMWFDATA v7.0.2 from flexpart.eu

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