source: flex_extract.git/src/preconvert.f90 @ d69b677

ctbtodev
Last change on this file since d69b677 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: 26.3 KB
Line 
1      PROGRAM PRECONVERT
2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3!                                                                 !
4! PROGRAM PRECONVERT - PREPARES INPUT DATA FOR POP MODEL METEOR-  !
5!                      OLOGICAL PREPROCESSOR                      !
6!                                                                 !
7!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
8!                                                                 !
9! CALCULATION OF ETAPOINT ON A REGULAR LAMDA/PHI GRID AND WRITING !
10! U,V,ETAPOINT,T,PS,Q,SD,MSL,TCC,10U, 10V, 2T,2D,LSP,CP,SSHF,SSR, !
11! EWSS,NSSS TO AN OUTPUT FILE (GRIB 1 or 2 FORMAT).               !
12!                                                                 !
13! AUTHORS: L. HAIMBERGER, G. WOTAWA, 1994-04                      !
14!                     adapted: A. BECK                            !
15!                     2003-05-11                                  !
16!          L. Haimberger 2006-12    V2.0                          !
17!                    modified to handle arbitrary regular grids   !
18!                    and T799 resolution data                     !
19!          L. Haimberger 2010-03    V4.0                          !
20!                    modified to grib edition 2 fields            !
21!                    and T1279 resolution data                    !
22!                                                                 !
23!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
24!                                                                 !
25! DESCRIPTION OF NEEDED INPUT:                                    !
26!                                                                 !
27! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
28!                                                                 !
29! 11    fort.11   T,U,V           regular lamda phi grid          !
30! 12    fort.12   D               regular lamda phi grid          !
31! 13    fort.13   LNSP          fort.13  spherical harmonics             !
32! 14    fort.14   SD,MSL,TCC,10U,                                 !
33!                 10V,2T,2D       regular lamda phi grid          !
34! 16    fort.16   LSP,CP,SSHF,                                    !
35!                 SSR,EWSS,NSSS   regular lamda phi grid          !
36! 17    fort.17   Q               regular lamda phi grid          !
37!                                                                 !
38!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
39!                                                                 !
40! DESCRIPTION OF OUTPUT:                                          !
41!                                                                 !
42! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
43!                                                                 !
44! 15    fort.15   U,V,ETA,T,PS,                                   !
45!                 Q,SD,MSL,TCC,                                   !
46!                 10U,10V,2T,2D,  regular lamda phi grid          !
47!                 LSP,CP,SSHF,                                    !
48!                 SSR,EWSS,NSSS                                   !
49!                                                                 !
50!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51!
52
53      USE PHTOGR
54      USE GRTOPH
55      USE FTRAFO
56      USE RWGRIB2
57      USE GRIB_API
58
59      IMPLICIT NONE
60
61      REAL, ALLOCATABLE, DIMENSION (:,:) :: LNPS
62      REAL, ALLOCATABLE, DIMENSION (:,:) :: Z
63      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T, UV , UV2
64      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: QA,OM,OMR
65      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: DIV, ETA,ETAR
66      REAL, ALLOCATABLE, DIMENSION (:,:) :: DPSDL, DPSDM
67      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: PS,DPSDT
68      REAL, ALLOCATABLE, DIMENSION (:,:,:) :: SURF,FLUX,OROLSM
69      REAL, ALLOCATABLE, DIMENSION (:) :: WSAVE,H,SINL,COSL,WSAVE2
70      REAL, ALLOCATABLE, DIMENSION (:) :: BREITE, GBREITE,AK, BK,pv
71
72! Arrays for Gaussian grid calculations
73
74      REAL  :: X1,X2,RMS,MW,SIG,LAM
75      REAL,ALLOCATABLE :: CUA(:,:,:),CVA(:,:,:)
76
77      REAL, ALLOCATABLE, DIMENSION (:,:) :: P,PP !,P2
78      REAL, ALLOCATABLE, DIMENSION (:,:) :: XMN,HILFUV
79      REAL, ALLOCATABLE, DIMENSION (:) :: LNPMN,LNPMN2,LNPMN3
80      REAL, ALLOCATABLE, DIMENSION (:) :: WEIGHT
81      REAL, ALLOCATABLE, DIMENSION (:,:) :: UGVG
82      REAL, ALLOCATABLE, DIMENSION (:,:) :: DG, ETAG
83      REAL, ALLOCATABLE, DIMENSION (:,:) :: GWSAVE
84      REAL, ALLOCATABLE, DIMENSION (:) :: PSG,HILF
85
86! end arrays for Gaussian grid calculations
87
88      INTEGER, ALLOCATABLE, DIMENSION (:) :: MLAT,MPSURF,MPFLUX,MPORO,MPAR
89      INTEGER, ALLOCATABLE :: GIFAX(:,:)
90
91      REAL PI,COSB,DAK,DBK,P00
92      REAL URLAR8,JMIN1,LLLAR8,MAXBMIN1,PIR8,DCOSB
93
94      INTEGER I,J,K,L,IERR,M,LTEST,MK,NGI,NGJ
95      INTEGER MFLUX,MSURF,MORO
96      INTEGER LUNIT,LUNIT2
97
98      INTEGER MAXL, MAXB, MLEVEL, LEVOUT,LEVMIN,LEVMAX
99      INTEGER MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH, MNAUF,META,METADIFF
100      INTEGER MDPDETA,METAPAR
101      REAL RLO0, RLO1, RLA0, RLA1
102      CHARACTER*300 MLEVELIST
103
104      INTEGER MAUF, MANF,IFAX(10)
105
106      INTEGER IGRIB(1),iret,ogrib
107
108      CHARACTER*80 FILENAME
109
110      NAMELIST /NAMGEN/ &
111                 MAXL, MAXB, &
112                 MLEVEL,MLEVELIST,MNAUF,METAPAR, &
113                 RLO0, RLO1, RLA0, RLA1, &
114                 MOMEGA,MOMEGADIFF,MGAUSS,MSMOOTH,META,METADIFF,&
115                 MDPDETA
116
117      LTEST=1
118
119      call posnam (4,'NAMGEN')
120      read (4,NAMGEN)
121
122      MAUF=INT(360.*(REAL(MAXL)-1.)/(RLO1-RLO0)+0.0001)
123!      PRINT*, MAUF
124
125      MANF=INT(REAL(MAUF)/360.*(360.+RLO0)+1.0001)
126      IF(MANF .gt. MAUF) MANF=MANF-MAUF
127
128
129!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130!                       ALLOCATE VARIABLES                        !
131!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
132
133      ALLOCATE (LNPS(0:(MNAUF+1)*(MNAUF+2)-1,1))
134
135      ALLOCATE (H(0:(MNAUF+2)*(MNAUF+3)/2))
136
137     
138      ALLOCATE (OM(MAXL, MAXB, MLEVEL))
139
140      ALLOCATE (ETA(MAXL,MAXB,MLEVEL))
141
142      ALLOCATE (PS(MAXL, MAXB,1),DPSDT(MAXL, MAXB,1))
143
144     
145      ALLOCATE (WSAVE(4*MAUF+15),WSAVE2(4*MAUF+15))
146
147      ALLOCATE (BREITE(MAXB),AK(MLEVEL+1),BK(MLEVEL+1),pv(2*mlevel+2))
148     
149      ALLOCATE (MPAR(2))
150
151      ALLOCATE (COSL(MAXL),SINL(MAXL))
152               
153      ALLOCATE (CUA(2,4,MLEVEL),CVA(2,4,MLEVEL))
154       
155!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156!     GAUSS STUFF !
157!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158
159
160      IF(MGAUSS .EQ. 1) THEN
161      LUNIT=0
162      FILENAME='fort.18'
163
164      call grib_open_file(LUNIT, TRIM(FILENAME),'r')
165 
166      call grib_new_from_file(LUNIT,igrib(1), iret)
167 
168! we can close the file
169      call grib_close_file(LUNIT)
170 
171!      call grib_get(igrib(1),'gridType', j)
172
173      NGJ=MNAUF+1
174
175      ALLOCATE (GWSAVE(8*NGJ+15,NGJ/2))
176      ALLOCATE(GIFAX(10,NGJ))
177      ALLOCATE (GBREITE(NGJ),WEIGHT(NGJ))
178      ALLOCATE (MLAT(NGJ))
179      ALLOCATE (P(0:((MNAUF+3)*(MNAUF+4))/2,NGJ/2))
180      ALLOCATE (PP(NGJ/2,0:((MNAUF+3)*(MNAUF+4))/2))
181      ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB))
182
183      call grib_get(igrib(1),'numberOfPointsAlongAMeridian', NGJ)
184 
185      !     get as a integer
186      call grib_get(igrib(1),'pl', MLAT)
187
188      NGI=SUM(MLAT)
189
190      call grib_get(igrib(1),'numberOfVerticalCoordinateValues',mk)
191
192      IF(mk/2-1 .ne. MLEVEL) THEN
193        WRITE(*,*) 'FATAL: Number of model levels',mk, &
194          ' does not agree with', MLEVEL,' in namelist'
195        STOP
196      ENDIF
197      call grib_get(igrib(1),'pv',pv)
198        AK=pv(1:1+MLEVEL)
199        BK=pv(2+MLEVEL:2*MLEVEL+2)
200
201      ALLOCATE (LNPMN(0:(MNAUF+1)*(MNAUF+2)-1))
202      ALLOCATE (LNPMN2(0:(MNAUF+1)*(MNAUF+2)-1))
203      ALLOCATE (UGVG(NGI, 2*MLEVEL),HILFUV(2*MAXL,2))
204
205
206      ALLOCATE (DPSDL(NGI,1),DPSDM(NGI,1))
207
208      ALLOCATE (PSG(NGI),HILF(NGI))
209      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
210!      ALLOCATE (UV2(MAXL, MAXB, 2*MLEVEL))
211
212      ALLOCATE (XMN(0:(MNAUF+1)*(MNAUF+2)-1, 2*MLEVEL))
213      ALLOCATE (DG(NGI,MLEVEL),ETAG(NGI,MLEVEL))
214
215! Initialisieren  Legendretransformation
216!       auf das LaT/LON Gitter
217
218      PI=ACOS(-1.D0)
219!$OMP PARALLEL DO
220      DO 20 J=1,MAXB
221
222      BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))* PI/180.D0)
223
224      CALL PLGNFA(MNAUF,BREITE(J),Z(0,J))
225
22620    CONTINUE
227!$OMP END PARALLEL DO
228
229! Avoid possible Pole problem
230!      IF(RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
231!      IF(RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
232
233! Initialisation of fields for  FFT and Legendre transformation
234!       to  Gaussian grid and back to phase space
235        X1=-1.D0
236        X2=1.D0
237        CALL GAULEG(X1,X2,GBREITE,WEIGHT,NGJ)
238
239!$OMP PARALLEL DO PRIVATE(M)
240        DO J=1,NGJ/2
241               CALL PLGNFA(MNAUF,GBREITE(J),P(:,J))
242                DO M=0,(MNAUF+3)*(MNAUF+4)/2
243                  PP(J,M)=P(M,J)
244        ENDDO
245        ENDDO
246!$OMP END PARALLEL DO
247
248
249!       MPAR(1)=152
250        FILENAME='fort.12'
251        CALL READSPECTRAL(FILENAME,LNPMN,MNAUF,1,MLEVEL,(/152/),AK,BK)
252!      goto 111
253        CALL SET99(WSAVE,IFAX,mauf)
254        CALL PHGCUT(LNPMN,PS,WSAVE,IFAX,Z, &
255      MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
256      CALL STATIS(MAXL,MAXB,1,EXP(PS),RMS,MW,SIG)
257      WRITE(*,'(A12,3F12.4)') 'STATISTICS: ',RMS,MW,SIG
258
259        DO J=1,NGJ/2
260          CALL SET99(GWSAVE(1,J),GIFAX(1,J),MLAT(J))
261        ENDDO
262        CALL PHGR213(LNPMN,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
263        PSG=HILF
264        CALL GRPH213(LNPMN2,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
265      MNAUF,NGI,NGJ,1)
266        CALL PHGR213(LNPMN2,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
267
268
269        HILF=exp(PSG)-exp(HILF)
270
271      CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
272      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
273
274        PSG=EXP(PSG)
275        HILF=PSG
276      CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
277      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
278
279  111         FILENAME='fort.10'
280       CALL READSPECTRAL(FILENAME, &
281      XMN,MNAUF,2*MLEVEL,MLEVEL,(/131,132/),AK,BK)
282!       Transformieren des Windes auf das Gaussgitter   
283        CALL PHGR213(XMN,UGVG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,2*MLEVEL)
284       DO K=1,MLEVEL
285! North Pole
286            CALL JSPPOLE(XMN(:,K),1,MNAUF,.TRUE.,CUA(:,:,K))
287            CALL JSPPOLE(XMN(:,MLEVEL+K),1,MNAUF,.TRUE.,CVA(:,:,K))
288! South Pole
289            CALL JSPPOLE(XMN(:,K),-1,MNAUF,.TRUE.,CUA(:,3:4,K))
290            CALL JSPPOLE(XMN(:,MLEVEL+K),-1,MNAUF,.TRUE.,CVA(:,3:4,K))
291       ENDDO
292   
293        DO K=1,2*MLEVEL
294          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
295        ENDDO
296        CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z, &
297      MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
298
299
300 112        FILENAME='fort.13'
301      CALL READSPECTRAL(FILENAME,XMN,MNAUF,MLEVEL,MLEVEL,(/155/),AK,BK)
302!       Transformieren der horizontalen Divergenz auf das Gaussgitter                   
303        CALL PHGR213(XMN,DG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,MLEVEL)
304
305
306!       Berechnung des Gradienten des Logarithmus des Bodendrucks
307!       auf dem Gaussgitter     
308        CALL PHGRAD(LNPMN,DPSDL,DPSDM,GWSAVE,GIFAX,P,H,MLAT,MNAUF,NGI,NGJ,1)
309
310!       Berechnung der Vertikalgeschwindigkeit auf dem Gaussgitter
311      CALL CONTGL(HILF,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
312      GBREITE,ETAG,MLAT,AK,BK,NGI,NGJ,MLEVEL)
313
314 
315        CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
316      MNAUF,NGI,NGJ,MLEVEL)
317        DO K=1,MLEVEL
318          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
319        ENDDO
320        CALL PHGCUT(XMN,ETA,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
321
322        CALL GRPH213(XMN,HILF,GWSAVE,GIFAX,PP,WEIGHT,MLAT, MNAUF,NGI,NGJ,1)
323
324        IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,1),MNAUF,MSMOOTH)
325        CALL PHGCUT(XMN,DPSDT,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
326!       GOTO 114
327
328      CALL STATIS(MAXL,MAXB,1,DPSDT,RMS,MW,SIG)
329      WRITE(*,'(A12,3F11.4)') 'STATISTICS DPSDT: ',RMS,MW,SIG
330
331       IF(MOMEGADIFF .ne. 0) THEN
332!       Berechnung von Omega auf dem Gaussgitter
333        CALL OMEGA(PSG,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
334      GBREITE,ETAG,MLAT,AK,BK,NGI ,NGJ,MLEVEL)
335
336        CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,&
337      MNAUF,NGI,NGJ,MLEVEL)
338        DO K=1,MLEVEL
339          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
340        ENDDO
341        CALL PHGCUT(XMN,OM,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
342
343       ENDIF !MOMEGA
344
345        CALL GRPH213(XMN,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
346        CALL PHGCUT(XMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
347
348      CALL STATIS(MAXL,MAXB,1,PS,RMS,MW,SIG)
349      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
350
351 114  DEALLOCATE(HILF,PSG,DPSDL,DPSDM,ETAG,DG,LNPMN)
352
353!      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
354!       CALL GRPH213(XMN,UGVG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,
355!     *MNAUF,NGI,NGJ,2*MLEVEL)
356!        DO K=1,2*MLEVEL
357!          IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
358!        ENDDO
359!        CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z,
360!     *MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
361        DEALLOCATE(PP,P,UGVG,MLAT,GBREITE,WEIGHT,GWSAVE,XMN)
362
363!        CALL ETAGAUSS(Z,WSAVE
364!     *,BREITE,UV,ETA,OM,PS,
365!     *MAUF,MAXB,MAXL,MANF,MNAUF,MLEVEL,MSMOOTH)
366
367      ELSE
368
369!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
370!          READING OF PREPARED METEOROLOGICAL FIELDS              !
371!                                                                 !
372!          THE FOLLOWING FIELDS ARE EXPECTED:                     !
373!                                                                 !
374!          UNIT 11: T,U,V        (REGULAR GRID)                   !
375!          UNIT 17: Q            (REGULAR GRID)                   !
376!          UNIT 13: D            (REGULAR GRID)                   !
377!          UNIT 12: LNSP         (SPHERICAL HARMONICS)            !
378!          UNIT 14: SURFACE DATA (REGULAR GRID)                   !
379!          UNIT 16: FLUX DATA    (REGULAR GRID)                   !
380!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
381!
382      ALLOCATE (MLAT(MAXB))
383      MLAT=MAXL
384      ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,1))
385      ALLOCATE (DPSDL(MAXL,MAXB),DPSDM(MAXL,MAXB))
386      ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL),DIV(MAXL,MAXB,MLEVEL))
387
388!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
389!                  READING OF SURFACE PRESSURE                    !
390!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
391
392      FILENAME='fort.12'
393        CALL READSPECTRAL(FILENAME,LNPS,MNAUF,1,MLEVEL,(/152/),AK,BK)
394
395!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
396!                      READING OF U,V                      !
397!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
398!
399! OPENING OF UNBLOCKED GRIB FILE
400!
401      FILENAME='fort.10'
402      CALL READLATLON(FILENAME,UV,MAXL,MAXB,2*MLEVEL,(/131,132/))
403
404
405      PI=ACOS(-1.D0)
406      DO J=1,MAXB
407
408        BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))*PI/180.D0)
409
410      ENDDO
411! Avoid possible Pole problem
412!      IF(RLA0 .EQ. -90.0) BREITE(MAXB)=sin(-89.99*PI/180.d0)
413!      IF(RLA1 .EQ. 90.0)  BREITE(1)=sin(89.99*PI/180.d0)
414
415      DO K=1,2*MLEVEL
416        DO J=1,MAXB
417          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
418          IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
419             RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
420            UV(:,J,K)=UV(:,J,K)/1.D6
421          else
422            UV(:,J,K)=UV(:,J,K)*COSB
423          endif
424      ENDDO
425      ENDDO
426
427!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428!                     READING OF LNSP on grid                !
429!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
430
431! For debugging only
432!      FILENAME='LNSPG_G.20060330.600'
433!      INQUIRE(FILE=FILENAME,EXIST=EX)
434!      CALL READLATLON(FILENAME,QA,
435!     *MAXL,MAXB,1,1,(/152/))
436
437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
438!                     READING OF DIVERGENCE                       !
439!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
440
441      IF(META .EQ. 0 .OR. METADIFF .EQ. 1) THEN
442        FILENAME='fort.13'
443        CALL READLATLON(FILENAME,DIV,MAXL,MAXB,MLEVEL,(/155/))
444      ENDIF
445
446
447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448!       CALCULATION OF ETAPOINT --> TOTAL TIME DERIVATIVE OF       !
449!      ECMWF VERTICAL COORDINATE ETA MULTIPLIED BY DERIVATIVE     !
450!      OF PRESSURE IN ETA DIRECTION                               !
451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452
453! Initialisieren  Legendretransformation
454!       auf das LaT/LON Gitter
455! Without Gaussian grid calculation Legendre Polynomials are calculated
456! only for one latitude to save space
457
458      DO J=1,MAXB
459
460        CALL PLGNFA(MNAUF,BREITE(J),Z(0,1))
461
462        CALL PHGCUT(LNPS,PS(:,J,1),WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,1,1)
463
464
465        IF(META .EQ. 0 .or. METADIFF .EQ. 1 ) THEN
466          CALL PHGRACUT(LNPS,DPSDL(:,J),DPSDM(:,J),WSAVE,IFAX,Z,H,MAUF, &
467      MNAUF,MAXL,1,MANF,1)
468        ENDIF
469      ENDDO
470     
471      PS=EXP(PS)
472
473! For debugging only
474      CALL STATIS(MAXL,MAXB,1,PS(:,:,1),RMS,MW,SIG)
475      WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
476     
477     
478       IF(MOMEGADIFF .ne. 0) THEN
479
480         CALL OMEGA(PS,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
481      BREITE,OM,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL)
482       ENDIF
483     
484       IF(META .EQ. 0 .OR. METADIFF .ne. 0) THEN
485         DPSDT=PS
486         CALL CONTGL(DPSDT,DPSDL,DPSDM,DIV,UV(:,:,1),UV(:,:,MLEVEL+1), &
487      BREITE,ETA,MLAT,AK,BK,MAXL*MAXB,MAXB,MLEVEL)
488       ENDIF
489
490      ENDIF ! MGAUSS
491
492! CREATE FILE VERTICAL.EC NEEDED BY POP MODEL
493
494      open(21,file='VERTICAL.EC')
495      write(21,'(a)')
496      write(21,'(a)') 'VERTICAL DISCRETIZATION OF POP MODEL'
497      write(21,'(a)')
498      write(21,'(i3,a)') MLEVEL,'   number of layers'
499      write(21,'(a)')
500      write(21,'(a)') '* A(NLEV+1)'
501      write(21,'(a)')
502      do 205 i=1,MLEVEL+1
503205      write(21,'(f18.12)') AK(I)
504      write(21,'(a)')
505      write(21,'(a)') '* B(NLEV+1)'
506      write(21,'(a)')
507      do 210 i=1,MLEVEL+1
508210      write(21,'(f18.12)') BK(I)
509      close(21)
510
511!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
512!                     READING OF OMEGA                       !
513!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
514
515      IF(MOMEGA .NE. 0 ) THEN
516
517
518
519         ALLOCATE (OMR(MAXL, MAXB, MLEVEL))
520 
521         FILENAME='fort.19'
522         CALL READLATLON(FILENAME,OMR,MAXL,MAXB,MLEVEL,(/135/))
523
524      IF(MOMEGADIFF .NE. 0 ) THEN
525
526      DO K=1,MLEVEL
527        CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)       
528        WRITE(*,'(A12,I3,3F11.4)') '       ETA: ',K,RMS,MW,SIG
529        CALL STATIS(MAXL,MAXB,1,OMR(:,:,K),RMS,MW,SIG)       
530        WRITE(*,'(A12,I3,3F11.4)') '     OMEGA: ',K,RMS,MW,SIG
531        CALL STATIS(MAXL,MAXB,1,OM(:,:,K)-OMR(:,:,K),RMS,MW,SIG)       
532        WRITE(*,'(A12,I3,3F11.4)') 'OMEGA DIFF: ',K,RMS,MW,SIG
533      ENDDO
534
535      ENDIF
536      ENDIF
537
538!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
539!                     READING OF ETA                       !
540!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
541
542      IF(META .NE. 0 ) THEN
543
544         ALLOCATE (ETAR(MAXL, MAXB, MLEVEL))
545 
546         P00=101325.
547         FILENAME='fort.21'
548         CALL READLATLON(FILENAME,ETAR,MAXL,MAXB,MLEVEL,(/77/))
549
550         if(MDPDETA .EQ. 1) THEN
551          DO K=1,MLEVEL
552           DAK=AK(K+1)-AK(K)
553           DBK=BK(K+1)-BK(K)
554           DO J=1,MAXB
555             DO I=1,MAXL
556               ETAR(I,J,K)=2*ETAR(I,J,K)*PS(I,J,1)*(DAK/PS(I,J,1)+DBK)/ &
557                          (DAK/P00+DBK)
558               IF(K .GT. 1) ETAR(I,J,K)=ETAR(I,J,K)-ETAR(I,J,K-1)
559             ENDDO
560           ENDDO
561          ENDDO
562         ENDIF
563
564        IF(METADIFF .NE. 0 ) THEN
565
566         DO K=1,MLEVEL
567          CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG)       
568          WRITE(*,'(A12,I3,3F11.4)') '       ETA: ',K,RMS,MW,SIG
569          CALL STATIS(MAXL,MAXB,1,ETAR(:,:,K),RMS,MW,SIG)       
570          WRITE(*,'(A12,I3,3F11.4)') '     ETAR: ',K,RMS,MW,SIG
571          CALL STATIS(MAXL,MAXB,1,ETA(:,:,K)-ETAR(:,:,K),RMS,MW,SIG)       
572          WRITE(*,'(A12,I3,3F11.4)') 'ETA DIFF: ',K,RMS,MW,SIG
573         ENDDO
574         DO K=1,MLEVEL
575          WRITE(*,'(I3,2F11.4)') K,ETA(1,MAXB/2,K),ETAR(1,MAXB/2,K)
576         ENDDO
577        ELSE
578          ETA=ETAR
579        ENDIF
580      ENDIF
581
582      ALLOCATE (T(MAXL, MAXB, MLEVEL))
583      ALLOCATE (QA(MAXL, MAXB, MLEVEL))
584!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585!                      READING OF T                      !
586!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
587!
588! OPENING OF UNBLOCKED GRIB FILE
589!
590      FILENAME='fort.11'
591      CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/))
592
593!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
594!                     READING OF SPECIFIC HUMIDITY                !
595!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596
597      FILENAME='fort.17'
598      CALL READLATLON(FILENAME,QA,MAXL,MAXB,MLEVEL,(/133/))
599
600!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
601!                     TEST READING OF UV from MARS (debug only)   !
602!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
603
604!      FILENAME='fort.22'
605!      CALL READLATLON(FILENAME,UV2,MAXL,MAXB,2*MLEVEL,2,(/131,132/))
606
607
608
609!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
610!                    WRITE MODEL LEVEL DATA TO fort.15            !
611!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
612
613!     Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi)
614!     Now we are transforming back to the usual winds.
615      DO K=1,MLEVEL
616        DO J=2,MAXB-1
617          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
618          UV(:,J,K)=UV(:,J,K)/COSB
619          UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
620        ENDDO
621! special treatment for poles, if necessary.
622        DO J=1,MAXB,MAXB-1
623          COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
624          if(1.0-BREITE(J)*BREITE(J) .gt. 0 .OR. MGAUSS .NE. 1) then
625            IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
626             RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
627             UV(:,J,K)=UV(:,J,K)*1.D6
628             UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6
629            else
630             UV(:,J,K)=UV(:,J,K)/COSB
631             UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
632            endif
633          else
634              HILFUV(5:MAXL,:)=0.
635              HILFUV(1:2,:)=0.
636            IF(J.EQ.MAXB) THEN
637! Suedpol
638              HILFUV(3:4,1)=CUA(:,4,K)
639              HILFUV(3:4,2)=CVA(:,4,K)
640            ELSE
641! Nordpol
642              HILFUV(3:4,1)=CUA(:,2,K)
643              HILFUV(3:4,2)=CVA(:,2,K)
644            ENDIF
645              CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
646              DO I=0,MAXL-1
647                IF(MANF+I.LE.MAXL) THEN
648                  UV(I+1,J,K)=HILFUV(MANF+I,1)
649                ELSE
650                  UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1)
651                ENDIF
652              ENDDO
653              CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL,-1)
654              DO I=0,MAXL-1
655                IF(MANF+I.LE.MAXL) THEN
656                  UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2)
657                ELSE
658                  UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2)
659                ENDIF
660              ENDDO
661          endif
662      ENDDO
663      ENDDO
664
665! open output file
666      call grib_open_file(LUNIT,'fort.15','w')
667
668! we use temperature on lat/lon on model levels as template for model level data
669      LUNIT2=0
670      call grib_open_file(LUNIT2,'fort.11','r')
671      call grib_new_from_file(LUNIT2,igrib(1), iret)
672      call grib_close_file(LUNIT2)
673
674
675      CALL WRITELATLON(LUNIT,igrib(1),ogrib,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/131/))
676
677      CALL WRITELATLON(LUNIT,igrib(1),ogrib,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/132/))
678
679      IF(MDPDETA .ne. 1 .AND. MGAUSS .EQ. 0 .and. META .eq. 1) THEN
680        CALL WRITELATLON(LUNIT,igrib(1),ogrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/77/))
681      ELSE
682        CALL WRITELATLON(LUNIT,igrib(1),ogrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/METAPAR/))
683      ENDIF
684
685      CALL WRITELATLON(LUNIT,igrib(1),ogrib,T,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/130/))
686
687      CALL WRITELATLON(LUNIT,igrib(1),ogrib,PS,MAXL,MAXB,1,'1',1,(/134/))
688 
689      call grib_set(igrib(1),"levelType","ml")
690      call grib_set(igrib(1),"typeOfLevel","hybrid")
691      CALL WRITELATLON(LUNIT,igrib(1),ogrib,QA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/133/))
692
693 
694      IF(MOMEGA .EQ. 1) THEN
695        call grib_open_file(LUNIT2,'fort.25','w')
696        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OMR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
697
698        IF(MOMEGADIFF .EQ. 1) THEN
699
700          CALL WRITELATLON(LUNIT2,igrib(1),ogrib,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
701
702        OM=OM-OMR
703        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
704      call grib_close_file(LUNIT2)
705        ENDIF
706      ENDIF
707
708      IF(META .EQ. 1 .and. METADIFF .EQ. 1) THEN
709        call grib_open_file(LUNIT2,'fort.26','w')
710        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,ETAR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/))
711
712!        IF(MOMEGADIFF .EQ. 1) THEN
713
714          CALL WRITELATLON(LUNIT2,igrib(1),ogrib,DPSDT,MAXL,MAXB,1,'1',1,(/158/))
715
716        OM=ETA-ETAR
717        CALL WRITELATLON(LUNIT2,igrib(1),ogrib,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/))
718      call grib_close_file(LUNIT2)
719!        ENDIF
720      ENDIF
721
722
723      call grib_close_file(LUNIT)
724
725
726
727 2000 STOP 'SUCCESSFULLY FINISHED CONVERT_PRE: CONGRATULATIONS'
728 3000 STOP 'ROUTINE CONVERT_PRE: ERROR'
729 9999 stop 'ROUTINE CONVERT_PRE: ERROR'
730      END
731
732     
733
734!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
735      INTEGER FUNCTION IA (FIELD1,NI,NJ,NK,G)
736
737      IMPLICIT NONE
738      INTEGER NI,NJ,NK,I,J,K
739      REAL FIELD1(NI,NJ,NK)
740      REAL G
741      REAL RMIN,RMAX,XMAX,A,A1,A2
742
743      RMAX=FIELD1(1,1,1)
744      RMIN=FIELD1(1,1,1)
745     
746      DO 100 K=1,NK
747      DO 100 J=1,NJ
748        DO 100 I=1,NI
749          IF (FIELD1(I,J,K).GT.RMAX)RMAX=FIELD1(I,J,K)
750          IF (FIELD1(I,J,K).LT.RMIN)RMIN=FIELD1(I,J,K)
751100   CONTINUE
752
753      IF (ABS(RMIN).GT.RMAX.OR.ABS(RMIN).EQ.RMAX) THEN
754      XMAX=ABS(RMIN)
755      ELSE
756      XMAX=RMAX
757      ENDIF
758     
759      IF (XMAX.EQ.0) THEN
760      IA = 0
761      RETURN
762      ENDIF
763     
764      A1=LOG10 ((G/10.d0)/XMAX)
765      A2=LOG10 ( G/XMAX )
766      IF(A1 .gt. A2) THEN
767        A=A2
768      ELSE
769        A=A1
770      ENDIF
771     
772      IF (A.GT.0) IA=INT(A)
773      IF (A.LT.0) IA=INT(A-1.0)
774     
775      RETURN
776      END
777     
778      SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG)
779        IMPLICIT REAL (A-H,O-Z)
780
781      REAL PHI(NI,NJ,NK),SIG,MW,RMS,P
782 
783      N=NI*NJ*NK
784 
785      RMS=0.
786        MW=0.
787 
788      DO 10 I=1,NI
789      DO 10 J=1,NJ
790        DO 10 K=1,NK
791              P=PHI(I,J,K)
792          RMS=RMS+P*P
793                MW=MW+P
79410    CONTINUE
795
796      RMS=SQRT(RMS/N)
797        MW=MW/N
798       
799        IF(RMS*RMS-MW*MW.LT.0.) THEN     
800          SIG=0.0
801        ELSE
802          SIG=SQRT(RMS*RMS-MW*MW)
803        ENDIF
804 
805      RETURN
806      END
807     
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG