Context Navigation

calc_etadot.f90 @ dfa7dbd

ctbtodev

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

changes in the Fortran part and associated regression test

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

all .f90 free format
code beautification
regression test is OK

make new local gfortran makefiles, remove parameters not needed
anymore

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

provide softlinks for standards:

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

provide changelog.txt in Fortran
provide readme.txt in Regression/FortranEtadot?

Property mode set to 100644

File size: 26.0 KB

Line
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(4MAUF+15),WSAVE2(4MAUF+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, 2MLEVEL),HILFUV(2MAXL,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, 2MLEVEL))
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,2MLEVEL)
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,2MLEVEL)
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)=2ETAR(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=NINJNK
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 (RMSRMS-MWMW .LT. 0.) THEN
807	SIG=0.0
808	ELSE
809	SIG=SQRT(RMSRMS-MWMW)
810	END IF
811
812	RETURN
813
814	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format