| 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 | ! 2020-06-25 Commented out by PS - not needed anymore |
|---|
| 490 | |
|---|
| 491 | ! OPEN(21,FILE='VERTICAL.EC') |
|---|
| 492 | ! WRITE(21,'(A)') |
|---|
| 493 | ! WRITE(21,'(A)') 'VERTICAL DISCRETIZATION OF POP MODEL' |
|---|
| 494 | ! WRITE(21,'(A)') |
|---|
| 495 | ! write(21,'(i3,a)') MLEVEL,' number of layers' |
|---|
| 496 | ! WRITE(21,'(A)') |
|---|
| 497 | ! WRITE(21,'(A)') '* A(NLEV+1)' |
|---|
| 498 | ! WRITE(21,'(A)') |
|---|
| 499 | ! DO 205 I=1,MLEVEL+1 |
|---|
| 500 | !205 WRITE(21,'(F18.12)') AK(I) |
|---|
| 501 | ! WRITE(21,'(A)') |
|---|
| 502 | ! WRITE(21,'(A)') '* B(NLEV+1)' |
|---|
| 503 | ! WRITE(21,'(A)') |
|---|
| 504 | ! DO 210 I=1,MLEVEL+1 |
|---|
| 505 | !210 WRITE(21,'(F18.12)') BK(I) |
|---|
| 506 | ! CLOSE(21) |
|---|
| 507 | |
|---|
| 508 | !------------------------------------------------------------------ |
|---|
| 509 | ! READING OF OMEGA |
|---|
| 510 | !------------------------------------------------------------------ |
|---|
| 511 | |
|---|
| 512 | IF (MOMEGA .NE. 0 ) THEN |
|---|
| 513 | |
|---|
| 514 | ALLOCATE (OMR(MAXL, MAXB, MLEVEL)) |
|---|
| 515 | |
|---|
| 516 | FILENAME='fort.19' |
|---|
| 517 | CALL READLATLON(FILENAME,OMR,MAXL,MAXB,MLEVEL,(/135/)) |
|---|
| 518 | |
|---|
| 519 | IF (MOMEGADIFF .NE. 0) THEN |
|---|
| 520 | |
|---|
| 521 | DO K=1,MLEVEL |
|---|
| 522 | CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG) |
|---|
| 523 | WRITE(*,'(A12,I3,3F12.4)') ' ETA: ',K,RMS,MW,SIG |
|---|
| 524 | CALL STATIS(MAXL,MAXB,1,OMR(:,:,K),RMS,MW,SIG) |
|---|
| 525 | WRITE(*,'(A12,I3,3F12.4)') ' OMEGA: ',K,RMS,MW,SIG |
|---|
| 526 | CALL STATIS(MAXL,MAXB,1,OM(:,:,K)-OMR(:,:,K),RMS,MW,SIG) |
|---|
| 527 | WRITE(*,'(A12,I3,3F12.4)') 'OMEGA DIFF: ',K,RMS,MW,SIG |
|---|
| 528 | END DO |
|---|
| 529 | |
|---|
| 530 | END IF |
|---|
| 531 | END IF |
|---|
| 532 | |
|---|
| 533 | !------------------------------------------------------------------ |
|---|
| 534 | ! READING OF ETA |
|---|
| 535 | !------------------------------------------------------------------ |
|---|
| 536 | |
|---|
| 537 | IF (META .NE. 0 ) THEN |
|---|
| 538 | |
|---|
| 539 | ALLOCATE (ETAR(MAXL, MAXB, MLEVEL)) |
|---|
| 540 | |
|---|
| 541 | P00=101325. |
|---|
| 542 | FILENAME='fort.21' |
|---|
| 543 | CALL READLATLON(FILENAME,ETAR,MAXL,MAXB,MLEVEL,(/77/)) |
|---|
| 544 | |
|---|
| 545 | IF(MDPDETA .EQ. 1) THEN |
|---|
| 546 | DO K=1,MLEVEL |
|---|
| 547 | DAK=AK(K+1)-AK(K) |
|---|
| 548 | DBK=BK(K+1)-BK(K) |
|---|
| 549 | DO J=1,MAXB |
|---|
| 550 | DO I=1,MAXL |
|---|
| 551 | ETAR(I,J,K)=2*ETAR(I,J,K)*PS(I,J,1)*(DAK/PS(I,J,1)+DBK)/ & |
|---|
| 552 | (DAK/P00+DBK) |
|---|
| 553 | IF (K .GT. 1) ETAR(I,J,K)=ETAR(I,J,K)-ETAR(I,J,K-1) |
|---|
| 554 | END DO |
|---|
| 555 | END DO |
|---|
| 556 | END DO |
|---|
| 557 | END IF |
|---|
| 558 | |
|---|
| 559 | IF (METADIFF .NE. 0 ) THEN |
|---|
| 560 | |
|---|
| 561 | DO K=1,MLEVEL |
|---|
| 562 | CALL STATIS(MAXL,MAXB,1,ETA(:,:,K),RMS,MW,SIG) |
|---|
| 563 | WRITE(*,'(A12,I3,3F12.4)') ' ETA: ',K,RMS,MW,SIG |
|---|
| 564 | CALL STATIS(MAXL,MAXB,1,ETAR(:,:,K),RMS,MW,SIG) |
|---|
| 565 | WRITE(*,'(A12,I3,3F12.4)') ' ETAR: ',K,RMS,MW,SIG |
|---|
| 566 | CALL STATIS(MAXL,MAXB,1,ETA(:,:,K)-ETAR(:,:,K),RMS,MW,SIG) |
|---|
| 567 | WRITE(*,'(A12,I3,3F12.4)') 'ETA DIFF: ',K,RMS,MW,SIG |
|---|
| 568 | END DO |
|---|
| 569 | DO K=1,MLEVEL |
|---|
| 570 | WRITE(*,'(I3,2F12.4)') K,ETA(1,MAXB/2,K),ETAR(1,MAXB/2,K) |
|---|
| 571 | END DO |
|---|
| 572 | ELSE |
|---|
| 573 | ETA=ETAR |
|---|
| 574 | END IF |
|---|
| 575 | END IF |
|---|
| 576 | |
|---|
| 577 | ALLOCATE (T(MAXL, MAXB, MLEVEL)) |
|---|
| 578 | ALLOCATE (QA(MAXL, MAXB, MLEVEL)) |
|---|
| 579 | |
|---|
| 580 | !------------------------------------------------------------------ |
|---|
| 581 | !! READING OF T |
|---|
| 582 | !------------------------------------------------------------------ |
|---|
| 583 | |
|---|
| 584 | ! OPENING OF UNBLOCKED GRIB FILE |
|---|
| 585 | |
|---|
| 586 | FILENAME='fort.11' |
|---|
| 587 | CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/)) |
|---|
| 588 | |
|---|
| 589 | !------------------------------------------------------------------ |
|---|
| 590 | !! READING OF SPECIFIC HUMIDITY |
|---|
| 591 | !------------------------------------------------------------------ |
|---|
| 592 | |
|---|
| 593 | FILENAME='fort.17' |
|---|
| 594 | CALL READLATLON(FILENAME,QA,MAXL,MAXB,MLEVEL,(/133/)) |
|---|
| 595 | |
|---|
| 596 | !------------------------------------------------------------------ |
|---|
| 597 | ! TEST READING OF UV from MARS (debug only) |
|---|
| 598 | !------------------------------------------------------------------ |
|---|
| 599 | ! FILENAME='fort.22' |
|---|
| 600 | ! CALL READLATLON(FILENAME,UV2,MAXL,MAXB,2*MLEVEL,2,(/131,132/)) |
|---|
| 601 | |
|---|
| 602 | !------------------------------------------------------------------ |
|---|
| 603 | !! WRITE MODEL LEVEL DATA TO fort.15 |
|---|
| 604 | !------------------------------------------------------------------ |
|---|
| 605 | |
|---|
| 606 | !! Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi) |
|---|
| 607 | !! Now we are transforming back to the usual winds. |
|---|
| 608 | |
|---|
| 609 | DO K=1,MLEVEL |
|---|
| 610 | DO J=2,MAXB-1 |
|---|
| 611 | COSB=SQRT(1.0-(BREITE(J))*(BREITE(J))) |
|---|
| 612 | UV(:,J,K)=UV(:,J,K)/COSB |
|---|
| 613 | UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB |
|---|
| 614 | END DO |
|---|
| 615 | |
|---|
| 616 | ! special treatment for poles, if necessary. |
|---|
| 617 | DO J=1,MAXB,MAXB-1 |
|---|
| 618 | COSB=SQRT(1.0-(BREITE(J))*(BREITE(J))) |
|---|
| 619 | IF (1.0-BREITE(J)*BREITE(J) .GT. 0 .OR. MGAUSS .NE. 1) THEN |
|---|
| 620 | IF (RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. & |
|---|
| 621 | RLA1 .EQ. 90.0 .AND. J .EQ. 1) THEN |
|---|
| 622 | UV(:,J,K)=UV(:,J,K)*1.D6 |
|---|
| 623 | UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6 |
|---|
| 624 | ELSE |
|---|
| 625 | UV(:,J,K)=UV(:,J,K)/COSB |
|---|
| 626 | UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB |
|---|
| 627 | END IF |
|---|
| 628 | ELSE |
|---|
| 629 | HILFUV(5:MAXL,:)=0. |
|---|
| 630 | HILFUV(1:2,:)=0. |
|---|
| 631 | IF (J.EQ.MAXB) THEN |
|---|
| 632 | ! Suedpol |
|---|
| 633 | HILFUV(3:4,1)=CUA(:,4,K) |
|---|
| 634 | HILFUV(3:4,2)=CVA(:,4,K) |
|---|
| 635 | ELSE |
|---|
| 636 | ! Nordpol |
|---|
| 637 | HILFUV(3:4,1)=CUA(:,2,K) |
|---|
| 638 | HILFUV(3:4,2)=CVA(:,2,K) |
|---|
| 639 | END IF |
|---|
| 640 | CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL,-1) |
|---|
| 641 | DO I=0,MAXL-1 |
|---|
| 642 | IF (MANF+I .LE. MAXL) THEN |
|---|
| 643 | UV(I+1,J,K)=HILFUV(MANF+I,1) |
|---|
| 644 | ELSE |
|---|
| 645 | UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1) |
|---|
| 646 | END IF |
|---|
| 647 | END DO |
|---|
| 648 | CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL,-1) |
|---|
| 649 | DO I=0,MAXL-1 |
|---|
| 650 | IF (MANF+I .LE. MAXL) THEN |
|---|
| 651 | UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2) |
|---|
| 652 | ELSE |
|---|
| 653 | UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2) |
|---|
| 654 | END IF |
|---|
| 655 | END DO |
|---|
| 656 | end if |
|---|
| 657 | END DO |
|---|
| 658 | END DO |
|---|
| 659 | |
|---|
| 660 | ! open output file |
|---|
| 661 | call grib_open_file(LUNIT,'fort.15','w') |
|---|
| 662 | |
|---|
| 663 | ! we use temperature on lat/lon on model levels as template for model level data |
|---|
| 664 | LUNIT2=0 |
|---|
| 665 | CALL GRIB_OPEN_FILE(LUNIT2,'fort.11','R') |
|---|
| 666 | CALL GRIB_NEW_FROM_FILE(LUNIT2,IGRIB(1), IRET) |
|---|
| 667 | CALL GRIB_CLOSE_FILE(LUNIT2) |
|---|
| 668 | |
|---|
| 669 | |
|---|
| 670 | CALL WRITELATLON & |
|---|
| 671 | (LUNIT,IGRIB(1),OGRIB,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/131/)) |
|---|
| 672 | |
|---|
| 673 | CALL WRITELATLON & |
|---|
| 674 | (LUNIT,IGRIB(1),OGRIB,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,1,(/132/)) |
|---|
| 675 | |
|---|
| 676 | IF (MDPDETA .ne. 1 .AND. MGAUSS .EQ. 0 .and. META .eq. 1) THEN |
|---|
| 677 | CALL WRITELATLON & |
|---|
| 678 | (LUNIT,IGRIB(1),OGRIB,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/77/)) |
|---|
| 679 | ELSE |
|---|
| 680 | CALL WRITELATLON & |
|---|
| 681 | (LUNIT,IGRIB(1),OGRIB,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/METAPAR/)) |
|---|
| 682 | END IF |
|---|
| 683 | |
|---|
| 684 | CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,T,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/130/)) |
|---|
| 685 | |
|---|
| 686 | CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,PS,MAXL,MAXB,1,'1',1,(/134/)) |
|---|
| 687 | |
|---|
| 688 | CALL GRIB_SET(IGRIB(1),"levelType","ml") |
|---|
| 689 | CALL GRIB_SET(IGRIB(1),"typeOfLevel","hybrid") |
|---|
| 690 | CALL WRITELATLON(LUNIT,IGRIB(1),OGRIB,QA,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/133/)) |
|---|
| 691 | |
|---|
| 692 | |
|---|
| 693 | IF (MOMEGA .EQ. 1) THEN |
|---|
| 694 | CALL GRIB_OPEN_FILE(LUNIT2,'fort.25','w') |
|---|
| 695 | CALL WRITELATLON & |
|---|
| 696 | (LUNIT2,IGRIB(1),OGRIB,OMR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/)) |
|---|
| 697 | |
|---|
| 698 | IF (MOMEGADIFF .EQ. 1) THEN |
|---|
| 699 | CALL WRITELATLON(LUNIT2,IGRIB(1),OGRIB,DPSDT,MAXL,MAXB,1,'1',1,(/158/)) |
|---|
| 700 | OM=OM-OMR |
|---|
| 701 | CALL WRITELATLON & |
|---|
| 702 | (LUNIT2,IGRIB(1),OGRIB,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/)) |
|---|
| 703 | CALL GRIB_CLOSE_FILE(LUNIT2) |
|---|
| 704 | END IF |
|---|
| 705 | END IF |
|---|
| 706 | |
|---|
| 707 | IF (META .EQ. 1 .AND. METADIFF .EQ. 1) THEN |
|---|
| 708 | CALL GRIB_OPEN_FILE(LUNIT2,'fort.26','w') |
|---|
| 709 | CALL WRITELATLON & |
|---|
| 710 | (LUNIT2,IGRIB(1),OGRIB,ETAR,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/135/)) |
|---|
| 711 | ! IF (MOMEGADIFF .EQ. 1) THEN |
|---|
| 712 | CALL WRITELATLON(LUNIT2,IGRIB(1),OGRIB,DPSDT,MAXL,MAXB,1,'1',1,(/158/)) |
|---|
| 713 | OM=ETA-ETAR |
|---|
| 714 | CALL WRITELATLON & |
|---|
| 715 | (LUNIT2,IGRIB(1),OGRIB,OM,MAXL,MAXB,MLEVEL,MLEVELIST,1,(/001/)) |
|---|
| 716 | CALL GRIB_CLOSE_FILE(LUNIT2) |
|---|
| 717 | ! END IF |
|---|
| 718 | END IF |
|---|
| 719 | |
|---|
| 720 | CALL GRIB_CLOSE_FILE(LUNIT) |
|---|
| 721 | |
|---|
| 722 | 2000 STOP 'SUCCESSFULLY FINISHED calc_etadot: CONGRATULATIONS' |
|---|
| 723 | 3000 STOP 'ROUTINE calc_etadot: ERROR' |
|---|
| 724 | 9999 stop 'ROUTINE calc_etadot: ERROR' |
|---|
| 725 | |
|---|
| 726 | END |
|---|
| 727 | |
|---|
| 728 | !------------------------------------------------------------------ |
|---|
| 729 | INTEGER FUNCTION IA (FIELD1,NI,NJ,NK,G) |
|---|
| 730 | |
|---|
| 731 | |
|---|
| 732 | !------------------------------------------------------------------ |
|---|
| 733 | !! Calculate something that is roughly log10( maxval(field1)/g ) [PS] |
|---|
| 734 | !------------------------------------------------------------------ |
|---|
| 735 | |
|---|
| 736 | |
|---|
| 737 | IMPLICIT NONE |
|---|
| 738 | |
|---|
| 739 | INTEGER :: I,J,K |
|---|
| 740 | INTEGER, INTENT(IN) :: NI,NJ,NK |
|---|
| 741 | REAL, INTENT(IN) :: FIELD1(NI,NJ,NK) |
|---|
| 742 | REAL, INTENT(IN) :: G |
|---|
| 743 | REAL :: RMIN,RMAX,XMAX,A,A1,A2 |
|---|
| 744 | |
|---|
| 745 | RMAX=FIELD1(1,1,1) |
|---|
| 746 | RMIN=FIELD1(1,1,1) |
|---|
| 747 | |
|---|
| 748 | DO 100 K=1,NK |
|---|
| 749 | DO 100 J=1,NJ |
|---|
| 750 | DO 100 I=1,NI |
|---|
| 751 | IF (FIELD1(I,J,K) .GT. RMAX) RMAX=FIELD1(I,J,K) |
|---|
| 752 | IF (FIELD1(I,J,K) .LT. RMIN) RMIN=FIELD1(I,J,K) |
|---|
| 753 | 100 CONTINUE |
|---|
| 754 | |
|---|
| 755 | IF (ABS(RMIN) .GT. RMAX .OR. ABS(RMIN) .EQ. RMAX) THEN |
|---|
| 756 | XMAX=ABS(RMIN) |
|---|
| 757 | ELSE |
|---|
| 758 | XMAX=RMAX |
|---|
| 759 | END IF |
|---|
| 760 | |
|---|
| 761 | IF (XMAX .EQ. 0) THEN |
|---|
| 762 | IA = 0 |
|---|
| 763 | RETURN |
|---|
| 764 | END IF |
|---|
| 765 | |
|---|
| 766 | A1=LOG10( (G/10.d0)/XMAX ) |
|---|
| 767 | A2=LOG10( G/XMAX ) |
|---|
| 768 | IF (A1 .gt. A2) THEN |
|---|
| 769 | A=A2 |
|---|
| 770 | ELSE |
|---|
| 771 | A=A1 |
|---|
| 772 | END IF |
|---|
| 773 | |
|---|
| 774 | IF (A .GT. 0) IA=INT(A) |
|---|
| 775 | IF (A .LT. 0) IA=INT(A-1.0) |
|---|
| 776 | |
|---|
| 777 | RETURN |
|---|
| 778 | |
|---|
| 779 | END |
|---|
| 780 | |
|---|
| 781 | SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG) |
|---|
| 782 | |
|---|
| 783 | !------------------------------------------------------------------ |
|---|
| 784 | !! calculate mean, rms, stdev |
|---|
| 785 | !------------------------------------------------------------------ |
|---|
| 786 | |
|---|
| 787 | IMPLICIT REAL (A-H,O-Z) |
|---|
| 788 | |
|---|
| 789 | REAL PHI(NI,NJ,NK),SIG,MW,RMS,P |
|---|
| 790 | |
|---|
| 791 | N=NI*NJ*NK |
|---|
| 792 | |
|---|
| 793 | RMS=0. |
|---|
| 794 | MW=0. |
|---|
| 795 | ! 10.86 sinstead of 11.04 sec |
|---|
| 796 | DO 10 K=1,NK |
|---|
| 797 | DO 10 J=1,NJ |
|---|
| 798 | DO 10 I=1,NI |
|---|
| 799 | P=PHI(I,J,K) |
|---|
| 800 | RMS=RMS+P*P |
|---|
| 801 | MW=MW+P |
|---|
| 802 | 10 CONTINUE |
|---|
| 803 | |
|---|
| 804 | RMS=SQRT(RMS/N) |
|---|
| 805 | MW=MW/N |
|---|
| 806 | |
|---|
| 807 | IF (RMS*RMS-MW*MW .LT. 0.) THEN |
|---|
| 808 | SIG=0.0 |
|---|
| 809 | ELSE |
|---|
| 810 | SIG=SQRT(RMS*RMS-MW*MW) |
|---|
| 811 | END IF |
|---|
| 812 | |
|---|
| 813 | RETURN |
|---|
| 814 | |
|---|
| 815 | END |
|---|
