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