source: flexpart.git/preproc/src/preconvert.f90 @ 5763793

FPv9.3.1FPv9.3.1b_testingFPv9.3.2devfp9.3.1-20161214-nc4grib2nc4_repairrelease-10univie
Last change on this file since 5763793 was 5763793, checked in by Anne Fouilloux <annefou@…>, 5 years ago

Add pre-processing programs and examples. This is still under development!

  • Property mode set to 100644
File size: 15.4 KB
Line 
1PROGRAM 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  !       FILE      PARAMETER(S)    DATA REPRESENTATION             !
28  !                                                                 !
29  !       fort.10   U,V             spherical harmonics             !
30  !       fort.11   T               regular lamda phi grid          !
31  !       fort.12   LNSP            spherical harmonics             !
32  !       fort.13   D               spherical harmonics             !
33  !       fort.14   SD,MSL,TCC,10U,                                 !
34  !                 10V,2T,2D       regular lamda phi grid          !
35  !       fort.18   Q               regular lamda phi grid          !
36  !                                                                 !
37!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
38  !                                                                 !
39  ! DESCRIPTION OF OUTPUT:                                          !
40  !                                                                 !
41  ! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
42  !                                                                 !
43  ! 15    fort.15   U,V,ETA,T,PS,                                   !
44  !                 Q,SD,MSL,TCC,                                   !
45  !                 10U,10V,2T,2D,  regular lamda phi grid          !
46  !                 LSP,CP,SSHF,                                    !
47  !                 SSR,EWSS,NSSS                                   !
48  !                                                                 !
49!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
50  !
51
52  USE PHTOGR
53  USE GRTOPH
54  USE FTRAFO
55  USE RWGRIB2
56  USE GRIB_API
57
58  IMPLICIT NONE
59
60  REAL, ALLOCATABLE, DIMENSION (:,:)   :: Z
61  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T, UV
62  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: DIV, ETA
63  REAL, ALLOCATABLE, DIMENSION (:,:)   :: DPSDL, DPSDM
64  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: PS,DPSDT
65  REAL, ALLOCATABLE, DIMENSION (:)     :: WSAVE,H
66  REAL, ALLOCATABLE, DIMENSION (:)     :: BREITE, GBREITE
67  REAL, ALLOCATABLE, DIMENSION (:)     :: AK, BK, PV
68  INTEGER                              :: NPV
69
70  ! Arrays for Gaussian grid calculations
71
72  REAL  :: X1,X2,RMS,MW,SIG,LAM
73  REAL,ALLOCATABLE :: CUA(:,:,:),CVA(:,:,:)
74
75  REAL, ALLOCATABLE, DIMENSION (:,:) :: P,PP !,P2
76  REAL, ALLOCATABLE, DIMENSION (:,:) :: XMN,HILFUV
77  REAL, ALLOCATABLE, DIMENSION (:)   :: LNPMN,LNPMN2,LNPMN3
78  REAL, ALLOCATABLE, DIMENSION (:)   :: WEIGHT
79  REAL, ALLOCATABLE, DIMENSION (:,:) :: UGVG
80  REAL, ALLOCATABLE, DIMENSION (:,:) :: DG, ETAG
81  REAL, ALLOCATABLE, DIMENSION (:,:) :: GWSAVE
82  REAL, ALLOCATABLE, DIMENSION (:)   :: PSG,HILF
83
84  ! end arrays for Gaussian grid calculations
85
86  INTEGER, ALLOCATABLE, DIMENSION (:) :: MLAT
87  INTEGER, ALLOCATABLE :: GIFAX(:,:)
88
89  REAL, PARAMETER   :: PI=ACOS(-1.D0)
90
91  REAL COSB,DAK,DBK,P00
92  REAL URLAR8,JMIN1,LLLAR8,MAXBMIN1,PIR8,DCOSB
93
94  INTEGER I,J,K,L,IERR,M,MK,NGI,NGJ
95  INTEGER LUNIT,LUNIT_OUT
96
97  INTEGER MAXL, MAXB, MLEVEL, LEVOUT,LEVMIN,LEVMAX
98  INTEGER MGAUSS,MSMOOTH, MNAUF,META
99  INTEGER MDPDETA,METAPAR
100  REAL RLO0, RLO1, RLA0, RLA1
101  CHARACTER*300 MLEVELIST
102
103  INTEGER MAUF, MANF,IFAX(10)
104
105  INTEGER igrib,iret
106
107  CHARACTER*80 FILENAME
108
109  NAMELIST /NAMGEN/ &
110       MAXL, MAXB,  &
111       MLEVEL,MLEVELIST,MNAUF,METAPAR, &
112       RLO0, RLO1, RLA0, RLA1, &
113       MGAUSS,MSMOOTH,META,&
114       MDPDETA
115
116  read (4,NAMGEN)
117
118  PRINT*, 'MAXL= ', MAXL, ' RLO0 = ', RLO0, ' RLO1 = ', RLO1
119  MAUF=INT(360.*(REAL(MAXL)-1.)/(RLO1-RLO0)+0.0001)
120  PRINT*, 'MAUF= ', MAUF
121
122  MANF=INT(REAL(MAUF)/360.*(360.+RLO0)+1.0001)
123  PRINT*, 'MANF= ', MANF
124
125  IF(MANF .gt. MAUF) MANF=MANF-MAUF
126
127
128!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129  !     GAUSS STUFF !
130!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131
132  ! computation of etadot on gaussian grid
133
134! a new grib message is loaded from an existing sample.
135! we want a regular gaussian grid
136! set the environment variable GRIB_SAMPLES_PATH
137  write (filename, "(A13,I0,A3,I0)") "EI_regular_gg",MNAUF+1,"_ml", MLEVEL
138  call grib_new_from_samples(igrib, FILENAME)
139
140  call grib_get(igrib,'numberOfPointsAlongAMeridian', NGJ)
141  ALLOCATE (MLAT(NGJ))
142
143  !     get as a integer
144  call grib_get(igrib,'pl', MLAT)
145  print*,'Number of points along a meridian  = ', NGJ
146  NGI=SUM(MLAT)
147
148  call grib_get(igrib,'numberOfVerticalCoordinateValues',mk)
149  print*, 'numberOfVerticalCoordinateValues = ', mk
150
151  IF(mk/2-1 .ne. MLEVEL) THEN
152     WRITE(*,*) 'FATAL: Number of model levels',mk, &
153          ' does not agree with', MLEVEL,' in namelist'
154     STOP
155  ENDIF
156
157  call grib_get_size(igrib,'pv',npv)
158  allocate(pv(npv))
159  call grib_get(igrib,'pv',pv)
160  ALLOCATE(AK(NPV/2))
161  ALLOCATE(BK(NPV/2))
162  AK=pv(:NPV/2)
163  BK=pv(NPV/2+1:)
164  deallocate(pv)
165! END GAUSS INFO
166
167
168
169  ! Initialization of Legendre transform on LAT/LON grid.
170
171  ALLOCATE (BREITE(MAXB))
172  ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB))
173  !$OMP PARALLEL DO
174  DO  J=1,MAXB
175     BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))* PI/180.D0)
176     CALL PLGNFA(MNAUF,BREITE(J),Z(0,J))
177  ENDDO
178  !$OMP END PARALLEL DO
179
180
181  ! Initialisation of fields for  FFT and Legendre transformation
182  !     to  Gaussian grid and back to phase space
183  ALLOCATE (GBREITE(NGJ),WEIGHT(NGJ))
184  ALLOCATE (P(0:((MNAUF+3)*(MNAUF+4))/2,NGJ/2))
185  ALLOCATE (PP(NGJ/2,0:((MNAUF+3)*(MNAUF+4))/2))
186  X1=-1.D0
187  X2=1.D0
188  CALL GAULEG(X1,X2,GBREITE,WEIGHT,NGJ)
189
190  !$OMP PARALLEL DO PRIVATE(M)
191  DO J=1,NGJ/2
192     CALL PLGNFA(MNAUF,GBREITE(J),P(:,J))
193     DO M=0,(MNAUF+3)*(MNAUF+4)/2
194        PP(J,M)=P(M,J)
195     ENDDO
196  ENDDO
197  !$OMP END PARALLEL DO
198
199
200!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201!   READ LNSP
202!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203  FILENAME='fort.12'
204  ALLOCATE (LNPMN(0:(MNAUF+1)*(MNAUF+2)-1))
205  CALL READSPECTRAL(FILENAME,LNPMN,MNAUF,1,MLEVEL,(/152/))
206  ! Call set99 an initialization routine that must be called once
207  ! before a sequence of calls to the fft routines
208  ALLOCATE (WSAVE(4*MAUF+15))
209  ALLOCATE (PS(MAXL, MAXB,1))
210  CALL SET99(WSAVE,IFAX,mauf)
211  CALL PHGCUT(LNPMN,PS,WSAVE,IFAX,Z, &
212       MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
213  CALL STATIS(MAXL,MAXB,1,EXP(PS),RMS,MW,SIG)
214  WRITE(*,'(A12,3F12.4)') 'STATISTICS : ',RMS,MW,SIG
215
216  ALLOCATE (GWSAVE(8*NGJ+15,NGJ/2))
217  ALLOCATE (GIFAX(10,NGJ))
218  DO J=1,NGJ/2
219     CALL SET99(GWSAVE(1,J),GIFAX(1,J),MLAT(J))
220  ENDDO
221  ALLOCATE (PSG(NGI),HILF(NGI))
222  ALLOCATE (LNPMN2(0:(MNAUF+1)*(MNAUF+2)-1))
223  CALL PHGR213(LNPMN,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
224  PSG=HILF
225  CALL GRPH213(LNPMN2,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
226       MNAUF,NGI,NGJ,1)
227  CALL PHGR213(LNPMN2,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
228
229
230  HILF=exp(PSG)-exp(HILF)
231
232  CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
233  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
234
235  PSG=EXP(PSG)
236  HILF=PSG
237  CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
238  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
239
240!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
241!   READ U/V
242!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
243  FILENAME='fort.10'
244  ALLOCATE (XMN(0:(MNAUF+1)*(MNAUF+2)-1, 2*MLEVEL))
245  CALL READSPECTRAL(FILENAME, &
246       XMN,MNAUF,2*MLEVEL,MLEVEL,(/131,132/))
247  ! Transform winds on gaussian grid
248  ALLOCATE (UGVG(NGI, 2*MLEVEL))
249  CALL PHGR213(XMN,UGVG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,2*MLEVEL)
250  ALLOCATE (CUA(2,4,MLEVEL))
251  ALLOCATE (CVA(2,4,MLEVEL))
252  DO K=1,MLEVEL
253     ! North Pole
254     CALL JSPPOLE(XMN(:,K),1,MNAUF,.TRUE.,CUA(:,:,K))
255     CALL JSPPOLE(XMN(:,MLEVEL+K),1,MNAUF,.TRUE.,CVA(:,:,K))
256     ! South Pole
257     CALL JSPPOLE(XMN(:,K),-1,MNAUF,.TRUE.,CUA(:,3:4,K))
258     CALL JSPPOLE(XMN(:,MLEVEL+K),-1,MNAUF,.TRUE.,CVA(:,3:4,K))
259  ENDDO
260
261  IF(MSMOOTH .ne. 0) THEN
262    DO K=1,2*MLEVEL
263       CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
264    ENDDO
265  ENDIF
266  ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
267  CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z, &
268       MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
269
270
271!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272!   READ Divergence
273!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
274  FILENAME='fort.13'
275  CALL READSPECTRAL(FILENAME,XMN,MNAUF,MLEVEL,MLEVEL,(/155/))
276  ! Tranform horizontal divergence on gaussian grid
277  ALLOCATE (DG(NGI,MLEVEL))
278  CALL PHGR213(XMN,DG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,MLEVEL)
279
280
281  ! compute gradient of LNSP (log of surface pressure) on gaussian grid
282  ALLOCATE (H(0:(MNAUF+2)*(MNAUF+3)/2))
283  ALLOCATE (DPSDL(NGI,1),DPSDM(NGI,1))
284  CALL PHGRAD(LNPMN,DPSDL,DPSDM,GWSAVE,GIFAX,P,H,MLAT,MNAUF,NGI,NGJ,1)
285
286  ! Compute vertical weed spind on gaussian grid
287  ALLOCATE (ETAG(NGI,MLEVEL))
288  CALL CONTGL(HILF,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
289       GBREITE,ETAG,MLAT,AK,BK,NGI,NGJ,MLEVEL)
290
291
292  CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
293       MNAUF,NGI,NGJ,MLEVEL)
294  IF(MSMOOTH .ne. 0) THEN
295    DO K=1,MLEVEL
296       CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
297    ENDDO
298  ENDIF
299  ALLOCATE (ETA(MAXL,MAXB,MLEVEL))
300  CALL PHGCUT(XMN,ETA,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
301
302  CALL GRPH213(XMN,HILF,GWSAVE,GIFAX,PP,WEIGHT,MLAT, MNAUF,NGI,NGJ,1)
303
304  IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,1),MNAUF,MSMOOTH)
305  ALLOCATE (DPSDT(MAXL, MAXB,1))
306  CALL PHGCUT(XMN,DPSDT,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
307
308  CALL STATIS(MAXL,MAXB,1,DPSDT,RMS,MW,SIG)
309  WRITE(*,'(A12,3F11.4)') 'STATISTICS DPSDT: ',RMS,MW,SIG
310
311  CALL GRPH213(XMN,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
312  CALL PHGCUT(XMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
313
314  CALL STATIS(MAXL,MAXB,1,PS,RMS,MW,SIG)
315  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
316
317114 DEALLOCATE(HILF,PSG,DPSDL,DPSDM,ETAG,DG,LNPMN)
318
319  DEALLOCATE(PP,P,UGVG,MLAT,GBREITE,WEIGHT,GWSAVE,XMN)
320
321  ! CREATE FILE VERTICAL.EC NEEDED BY POP MODEL
322
323  open(21,file='VERTICAL.EC')
324  write(21,'(a)')
325  write(21,'(a)') 'VERTICAL DISCRETIZATION OF POP MODEL'
326  write(21,'(a)')
327  write(21,'(i3,a)') MLEVEL,'   number of layers'
328  write(21,'(a)')
329  write(21,'(a)') '* A(NLEV+1)'
330  write(21,'(a)')
331  do  i=1,MLEVEL+1
332     write(21,'(f18.12)') AK(I)
333  enddo
334  write(21,'(a)')
335  write(21,'(a)') '* B(NLEV+1)'
336  write(21,'(a)')
337  do  i=1,MLEVEL+1
338     write(21,'(f18.12)') BK(I)
339  enddo
340  close(21)
341
342
343  !     Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi)
344  !     Now we are transforming back to the usual winds.
345  ALLOCATE (HILFUV(2*MAXL,2))
346  DO K=1,MLEVEL
347     DO J=2,MAXB-1
348        COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
349        UV(:,J,K)=UV(:,J,K)/COSB
350        UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
351     ENDDO
352     ! special treatment for poles, if necessary.
353     DO J=1,MAXB,MAXB-1
354        COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
355        if(1.0-BREITE(J)*BREITE(J) .gt. 0 .OR. MGAUSS .NE. 1) then
356           IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
357                RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
358              UV(:,J,K)=UV(:,J,K)*1.D6
359              UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6
360           else
361              UV(:,J,K)=UV(:,J,K)/COSB
362              UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
363           endif
364        else
365           HILFUV(5:MAXL,:)=0.
366           HILFUV(1:2,:)=0.
367           IF(J.EQ.MAXB) THEN
368              ! Suedpol
369              HILFUV(3:4,1)=CUA(:,4,K)
370              HILFUV(3:4,2)=CVA(:,4,K)
371           ELSE
372              ! Nordpol
373              HILFUV(3:4,1)=CUA(:,2,K)
374              HILFUV(3:4,2)=CVA(:,2,K)
375           ENDIF
376           CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL)
377           DO I=0,MAXL-1
378              IF(MANF+I.LE.MAXL) THEN
379                 UV(I+1,J,K)=HILFUV(MANF+I,1)
380              ELSE
381                 UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1)
382              ENDIF
383           ENDDO
384           CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL)
385           DO I=0,MAXL-1
386              IF(MANF+I.LE.MAXL) THEN
387                 UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2)
388              ELSE
389                 UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2)
390              ENDIF
391           ENDDO
392        endif
393     ENDDO
394  ENDDO
395
396!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397  !                      READING OF T                      !
398!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399  ALLOCATE (T(MAXL, MAXB, MLEVEL))
400  FILENAME='fort.11'
401  CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/))
402
403!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
404  !                    WRITE MODEL LEVEL DATA TO fort.15            !
405!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
406  ! open output file
407  call grib_open_file(LUNIT,'fort.15','w')
408
409  ! we use temperature on lat/lon on model levels as template for model level data
410  LUNIT_OUT=0
411  call grib_open_file(LUNIT_OUT,'fort.11','r')
412  call grib_new_from_file(LUNIT_OUT,igrib, iret)
413  call grib_close_file(LUNIT_OUT)
414
415
416  CALL WRITELATLON(LUNIT,igrib,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,131)
417
418  CALL WRITELATLON(LUNIT,igrib,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,132)
419
420  CALL WRITELATLON(LUNIT,igrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,METAPAR)
421
422  CALL WRITELATLON(LUNIT,igrib,T,MAXL,MAXB,MLEVEL,MLEVELIST,130)
423
424  CALL WRITELATLON(LUNIT,igrib,PS,MAXL,MAXB,1,'1',134)
425
426  call grib_close_file(LUNIT)
427
428
4292000 STOP 'SUCCESSFULLY FINISHED CONVERT_PRE: CONGRATULATIONS'
4303000 STOP 'ROUTINE CONVERT_PRE: ERROR'
4319999 stop 'ROUTINE CONVERT_PRE: ERROR'
432
433CONTAINS
434
435  SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG)
436    IMPLICIT NONE
437
438    integer, intent(in) :: NI, NJ, NK
439
440    REAL PHI(NI,NJ,NK),SIG,MW,RMS,P
441    INTEGER N
442
443    N=NI*NJ*NK
444
445    RMS=0.
446    MW=0.
447
448    DO  I=1,NI
449       DO  J=1,NJ
450          DO  K=1,NK
451             P=PHI(I,J,K)
452             RMS=RMS+P*P
453             MW=MW+P
454          ENDDO
455       ENDDO
456    ENDDO
457
458    RMS=SQRT(RMS/N)
459    MW=MW/N
460
461    IF(RMS*RMS-MW*MW.LT.0.) THEN         
462       SIG=0.0
463    ELSE
464       SIG=SQRT(RMS*RMS-MW*MW)
465    ENDIF
466
467    RETURN
468  END SUBROUTINE STATIS
469
470END PROGRAM PRECONVERT
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG