Changes in / [4138764:4c52d0b] in flexpart.git
- Location:
- src
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
src/calcmatrix.f90
r92fab65 r9ca6e38 67 67 68 68 phconv(1) = psconv 69 ! Emanuel subroutine needs pressure in hPa, therefore convert all pressures 70 do kuvz = 2,nuvz 71 k = kuvz-1 72 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 73 pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv) 74 phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv) 75 else 76 phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k)) 77 endif 78 dpr(k) = phconv(k) - phconv(kuvz) 79 qsconv(k) = f_qvsat( pconv(k), tconv(k) ) 69 ! Emanuel subroutine needs pressure in hPa, therefore convert all pressures 70 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 71 do kuvz = 2,nuvz 72 k = kuvz-1 73 pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv) 74 phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv) 75 dpr(k) = phconv(k) - phconv(kuvz) 76 qsconv(k) = f_qvsat( pconv(k), tconv(k) ) 77 end do 78 else 79 ! JMA / SH Bugfix phconv was set in loop with access to undefined pconv(nuvz) 80 phconv(2:nuvz-1) = 0.5*(pconv(2:nuvz-1)+pconv(1:nuvz-2)) 81 phconv(nuvz) = pconv(nuvz-1) 82 dpr(1:nuvz-1) = phconv(1:nuvz-1)-phconv(2:nuvz) 80 83 81 ! initialize mass fractions 82 do kk=1,nconvlev 83 fmassfrac(k,kk)=0. 84 do k = 1,nuvz-1 85 qsconv(k) = f_qvsat( pconv(k), tconv(k) ) 84 86 end do 85 end do 86 87 end if 88 89 ! initialize mass fractions 90 fmassfrac(1:nuvz-1,1:nconvlev)=0.0 87 91 88 92 !note that Emanuel says it is important … … 93 97 !****************** 94 98 95 96 97 99 cbmfold = cbmf 100 ! Convert pressures to hPa, as required by Emanuel scheme 101 !******************************************************** 98 102 !!$ do k=1,nconvlev !old 99 do k=1,nconvlev+1 !bugfix 100 pconv_hpa(k)=pconv(k)/100. 101 phconv_hpa(k)=phconv(k)/100. 103 do k=1,nconvlev+1 !bugfix 104 pconv_hpa(k)=pconv(k)/100. 105 phconv_hpa(k)=phconv(k)/100. 106 end do 107 phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100. 108 call convect(nconvlevmax, nconvlev, delt, iflag, & 109 precip, wd, tprime, qprime, cbmf) 110 111 ! do not update fmassfrac and cloudbase massflux 112 ! if no convection takes place or 113 ! if a CFL criterion is violated in convect43c.f 114 if (iflag .ne. 1 .and. iflag .ne. 4) then 115 cbmf=cbmfold 116 goto 200 117 endif 118 119 ! do not update fmassfrac and cloudbase massflux 120 ! if the old and the new cloud base mass 121 ! fluxes are zero 122 if (cbmf.le.0..and.cbmfold.le.0.) then 123 cbmf=cbmfold 124 goto 200 125 endif 126 127 ! Update fmassfrac 128 ! account for mass displaced from level k to level k 129 130 lconv = .true. 131 do k=1,nconvtop 132 rlevmass = dpr(k)/ga 133 summe = 0. 134 do kk=1,nconvtop 135 fmassfrac(k,kk) = delt*fmass(k,kk) 136 summe = summe + fmassfrac(k,kk) 102 137 end do 103 phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100. 104 call convect(nconvlevmax, nconvlev, delt, iflag, & 105 precip, wd, tprime, qprime, cbmf) 138 fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe 139 end do 106 140 107 ! do not update fmassfrac and cloudbase massflux 108 ! if no convection takes place or 109 ! if a CFL criterion is violated in convect43c.f 110 if (iflag .ne. 1 .and. iflag .ne. 4) then 111 cbmf=cbmfold 112 goto 200 113 endif 114 115 ! do not update fmassfrac and cloudbase massflux 116 ! if the old and the new cloud base mass 117 ! fluxes are zero 118 if (cbmf.le.0..and.cbmfold.le.0.) then 119 cbmf=cbmfold 120 goto 200 121 endif 122 123 ! Update fmassfrac 124 ! account for mass displaced from level k to level k 125 126 lconv = .true. 127 do k=1,nconvtop 128 rlevmass = dpr(k)/ga 129 summe = 0. 130 do kk=1,nconvtop 131 fmassfrac(k,kk) = delt*fmass(k,kk) 132 summe = summe + fmassfrac(k,kk) 133 end do 134 fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe 135 end do 136 137 200 continue 141 200 continue 138 142 139 143 end subroutine calcmatrix -
src/ew.f90
r92fab65 r467460a 14 14 15 15 ew=0. 16 if(x.le.0.) stop 'sorry: t not in [k]' 16 if(x.le.0.) then 17 write(*,*) 'in ew.f90: x=', x 18 stop 'sorry: t not in [k]' 19 end if 20 17 21 y=373.16/x 18 22 a=-7.90298*(y-1.) -
src/gridcheck_gfs.f90
- Property mode changed from 100644 to 100755
ra756649 r467460a 75 75 real :: sizesouth,sizenorth,xauxa,pint 76 76 real :: akm_usort(nwzmax) 77 real,parameter :: eps= 0.000177 real,parameter :: eps=spacing(2.0_4*360.0_4) 78 78 79 79 ! NCEP GFS 80 80 real :: pres(nwzmax), help 81 81 82 integer :: i1 79,i180,i18182 integer :: i180 83 83 84 84 ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING … … 224 224 nxfield=isec2(2) 225 225 ny=isec2(3) 226 if((abs(xaux1).lt.eps).and.(xaux2.ge.359 )) then ! NCEP DATA FROM 0 TO227 xaux1=-1 79.0 ! 359 DEG EAST ->228 xaux2=-1 79.0+360.-360./real(nxfield) ! TRANSFORMED TO -179226 if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO 227 xaux1=-180.0 ! 359 DEG EAST -> 228 xaux2=-180.0+360.-360./real(nxfield) ! TRANSFORMED TO -179 229 229 endif ! TO 180 DEG EAST 230 230 if (xaux1.gt.180) xaux1=xaux1-360.0 … … 305 305 306 306 307 i179=nint(179./dx) 308 if (dx.lt.0.7) then 309 i180=nint(180./dx)+1 ! 0.5 deg data 310 else 311 i180=nint(179./dx)+1 ! 1 deg data 312 endif 313 i181=i180+1 307 i180=nint(180./dx) ! 0.5 deg data 314 308 315 309 … … 321 315 do ix=0,nxfield-1 322 316 help=zsec4(nxfield*(ny-jy-1)+ix+1) 323 if(ix.l e.i180) then324 oro(i1 79+ix,jy)=help325 excessoro(i1 79+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED317 if(ix.lt.i180) then 318 oro(i180+ix,jy)=help 319 excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 326 320 else 327 oro(ix-i18 1,jy)=help328 excessoro(ix-i18 1,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED321 oro(ix-i180,jy)=help 322 excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 329 323 endif 330 324 end do … … 339 333 do ix=0,nxfield-1 340 334 help=zsec4(nxfield*(ny-jy-1)+ix+1) 341 if(ix.l e.i180) then342 lsm(i1 79+ix,jy)=help335 if(ix.lt.i180) then 336 lsm(i180+ix,jy)=help 343 337 else 344 lsm(ix-i18 1,jy)=help338 lsm(ix-i180,jy)=help 345 339 endif 346 340 end do … … 413 407 write(*,*) 414 408 write(*,*) 415 write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &409 write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', & 416 410 nuvz,nwz 417 411 write(*,*) -
src/par_mod.f90
r759df5f raa939a9 147 147 ! GFS 148 148 integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 149 integer :: nxshift=0 ! shift not fixed for the executable 149 ! GFS 0.25 150 ! integer,parameter :: nxmax=1441,nymax=721,nuvzmax=138,nwzmax=138,nzmax=138 151 integer :: nxshift=0 ! shift not fixed for the executable 150 152 151 153 -
src/readreleases.f90
- Property mode changed from 100644 to 100755
r92fab65 r467460a 546 546 if ((jul1.lt.bdate).or.(jul2.gt.edate)) then 547 547 write(*,*) 'FLEXPART MODEL ERROR' 548 write(*,*) 'Release starts before simulation begins or ends '548 write(*,*) 'Release starts before simulation begins or ends (1)' 549 549 write(*,*) 'after simulation stops.' 550 write(*,*) 'Make files COMMAND and RELEASES consistent .'550 write(*,*) 'Make files COMMAND and RELEASES consistent (1).' 551 551 stop 552 552 endif … … 561 561 if ((jul1.lt.edate).or.(jul2.gt.bdate)) then 562 562 write(*,*) 'FLEXPART MODEL ERROR' 563 write(*,*) 'Release starts before simulation begins or ends '563 write(*,*) 'Release starts before simulation begins or ends (2)' 564 564 write(*,*) 'after simulation stops.' 565 write(*,*) 'Make files COMMAND and RELEASES consistent .'565 write(*,*) 'Make files COMMAND and RELEASES consistent (2).' 566 566 stop 567 567 endif -
src/readwind_gfs.f90
- Property mode changed from 100644 to 100755
ra756649 r9ca6e38 4 4 subroutine readwind_gfs(indj,n,uuh,vvh,wwh) 5 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 6 !*********************************************************************** 7 !* * 8 !* TRAJECTORY MODEL SUBROUTINE READWIND * 9 !* * 10 !*********************************************************************** 11 !* * 12 !* AUTHOR: G. WOTAWA * 13 !* DATE: 1997-08-05 * 14 !* LAST UPDATE: 2000-10-17, Andreas Stohl * 15 !* CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and * 16 !* qvh (on eta coordinates) in common block * 17 !* CHANGE: 16/11/2005, Caroline Forster, GFS data * 18 !* CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2 * 19 !* data with the ECMWF grib_api library * 20 !* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * 21 !* ECMWF grib_api * 22 ! * 23 ! Unified ECMWF and GFS builds * 24 ! Marian Harustak, 12.5.2017 * 25 ! - Renamed routine from readwind to readwind_gfs * 26 !* * 27 !*********************************************************************** 28 !* * 29 !* DESCRIPTION: * 30 !* * 31 !* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE * 32 !* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE * 33 !* * 34 !* INPUT: * 35 !* indj indicates number of the wind field to be read in * 36 !* n temporal index for meteorological fields (1 to 3)* 37 !* * 38 !* IMPORTANT VARIABLES FROM COMMON BLOCK: * 39 !* * 40 !* wfname File name of data to be read in * 41 !* nx,ny,nuvz,nwz expected field dimensions * 42 !* nlev_ec number of vertical levels ecmwf model * 43 !* uu,vv,ww wind fields * 44 !* tt,qv temperature and specific humidity * 45 !* ps surface pressure * 46 !* * 47 !*********************************************************************** 48 48 49 49 use eccodes … … 53 53 implicit none 54 54 55 55 !HSO new parameters for grib_api 56 56 integer :: ifile 57 57 integer :: iret 58 58 integer :: igrib 59 integer :: gribVer,parCat,parNum,typSurf, valSurf,discipl60 59 integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf 60 !HSO end edits 61 61 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) 62 62 real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax) … … 64 64 integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax 65 65 66 66 ! NCEP 67 67 integer :: numpt,numpu,numpv,numpw,numprh,numpclwch 68 68 real :: help, temp, ew … … 72 72 real :: qvh2(0:nxmax-1,0:nymax-1) 73 73 74 integer :: i1 79,i180,i18175 76 77 74 integer :: i180 75 76 ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING 77 !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input 78 78 79 79 integer :: isec1(8),isec2(3) 80 real :: xsec18 80 81 real(kind=4) :: zsec4(jpunp) 81 82 real(kind=4) :: xaux,yaux,xaux0,yaux0 83 real,parameter :: eps=spacing(2.0_4*360.0_4) 82 84 real(kind=8) :: xauxin,yauxin 83 real,parameter :: eps=1.e-484 85 real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1) 85 86 real :: plev1,hlev1,ff10m,fflev1 … … 87 88 logical :: hflswitch,strswitch 88 89 89 90 !HSO for grib api error messages 90 91 character(len=24) :: gribErrorMsg = 'Error reading grib file' 91 92 character(len=20) :: gribFunction = 'readwind_gfs' … … 100 101 101 102 102 103 104 103 ! OPENING OF DATA FILE (GRIB CODE) 104 105 !HSO 105 106 call grib_open_file(ifile,path(3)(1:length(3)) & 106 107 //trim(wfname(indj)),'r',iret) 107 108 if (iret.ne.GRIB_SUCCESS) then 108 109 goto 888 ! ERROR DETECTED 109 110 endif 110 111 !turn on support for multi fields messages 111 112 call grib_multi_support_on 112 113 … … 118 119 numpclwch=0 119 120 ifield=0 120 10 121 122 123 121 10 ifield=ifield+1 122 ! 123 ! GET NEXT FIELDS 124 ! 124 125 call grib_new_from_file(ifile,igrib,iret) 125 126 if (iret.eq.GRIB_END_OF_FILE) then … … 129 130 endif 130 131 131 132 !first see if we read GRIB1 or GRIB2 132 133 call grib_get_int(igrib,'editionNumber',gribVer,iret) 133 134 ! call grib_check(iret,gribFunction,gribErrorMsg) … … 135 136 if (gribVer.eq.1) then ! GRIB Edition 1 136 137 137 !read the grib1 identifiers 138 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 139 ! call grib_check(iret,gribFunction,gribErrorMsg) 140 call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret) 141 ! call grib_check(iret,gribFunction,gribErrorMsg) 142 call grib_get_int(igrib,'level',isec1(8),iret) 143 ! call grib_check(iret,gribFunction,gribErrorMsg) 138 !read the grib1 identifiers 139 140 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 141 call grib_check(iret,gribFunction,gribErrorMsg) 142 call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret) 143 call grib_check(iret,gribFunction,gribErrorMsg) 144 call grib_get_int(igrib,'level',isec1(8),iret) 145 call grib_check(iret,gribFunction,gribErrorMsg) 146 !JMA / SH: isec1(8) not evaluated any more below 147 !b/c with GRIB 2 this may be a real variable 148 xsec18 = real(isec1(8)) 144 149 145 150 else ! GRIB Edition 2 146 151 147 !read the grib2 identifiers 148 call grib_get_string(igrib,'shortName',shortname,iret) 149 150 call grib_get_int(igrib,'discipline',discipl,iret) 151 ! call grib_check(iret,gribFunction,gribErrorMsg) 152 call grib_get_int(igrib,'parameterCategory',parCat,iret) 153 ! call grib_check(iret,gribFunction,gribErrorMsg) 154 call grib_get_int(igrib,'parameterNumber',parNum,iret) 155 ! call grib_check(iret,gribFunction,gribErrorMsg) 156 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 157 ! call grib_check(iret,gribFunction,gribErrorMsg) 158 call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & 159 valSurf,iret) 160 ! call grib_check(iret,gribFunction,gribErrorMsg) 161 162 ! write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname 163 !convert to grib1 identifiers 164 isec1(6)=-1 165 isec1(7)=-1 166 isec1(8)=-1 167 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T 168 isec1(6)=11 ! indicatorOfParameter 169 isec1(7)=100 ! indicatorOfTypeOfLevel 170 isec1(8)=valSurf/100 ! level, convert to hPa 171 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U 172 isec1(6)=33 ! indicatorOfParameter 173 isec1(7)=100 ! indicatorOfTypeOfLevel 174 isec1(8)=valSurf/100 ! level, convert to hPa 175 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V 176 isec1(6)=34 ! indicatorOfParameter 177 isec1(7)=100 ! indicatorOfTypeOfLevel 178 isec1(8)=valSurf/100 ! level, convert to hPa 179 elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W 180 isec1(6)=39 ! indicatorOfParameter 181 isec1(7)=100 ! indicatorOfTypeOfLevel 182 isec1(8)=valSurf/100 ! level, convert to hPa 183 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH 184 isec1(6)=52 ! indicatorOfParameter 185 isec1(7)=100 ! indicatorOfTypeOfLevel 186 isec1(8)=valSurf/100 ! level, convert to hPa 187 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2 188 isec1(6)=52 ! indicatorOfParameter 189 isec1(7)=105 ! indicatorOfTypeOfLevel 190 isec1(8)=2 191 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2 192 isec1(6)=11 ! indicatorOfParameter 193 isec1(7)=105 ! indicatorOfTypeOfLevel 194 isec1(8)=2 195 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10 196 isec1(6)=33 ! indicatorOfParameter 197 isec1(7)=105 ! indicatorOfTypeOfLevel 198 isec1(8)=10 199 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10 200 isec1(6)=34 ! indicatorOfParameter 201 isec1(7)=105 ! indicatorOfTypeOfLevel 202 isec1(8)=10 203 elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]: 204 isec1(6)=153 ! indicatorOfParameter 205 isec1(7)=100 ! indicatorOfTypeOfLevel 206 isec1(8)=valSurf/100 ! level, convert to hPa 207 elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP 208 isec1(6)=2 ! indicatorOfParameter 209 isec1(7)=102 ! indicatorOfTypeOfLevel 210 isec1(8)=0 211 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP 212 isec1(6)=1 ! indicatorOfParameter 213 isec1(7)=1 ! indicatorOfTypeOfLevel 214 isec1(8)=0 215 elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW 216 isec1(6)=66 ! indicatorOfParameter 217 isec1(7)=1 ! indicatorOfTypeOfLevel 218 isec1(8)=0 219 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0 220 isec1(6)=11 ! indicatorOfParameter 221 isec1(7)=107 ! indicatorOfTypeOfLevel 222 isec1(8)=0.995 ! lowest sigma level 223 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0 224 isec1(6)=33 ! indicatorOfParameter 225 isec1(7)=107 ! indicatorOfTypeOfLevel 226 isec1(8)=0.995 ! lowest sigma level 227 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0 228 isec1(6)=34 ! indicatorOfParameter 229 isec1(7)=107 ! indicatorOfTypeOfLevel 230 isec1(8)=0.995 ! lowest sigma level 231 elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO 232 isec1(6)=7 ! indicatorOfParameter 233 isec1(7)=1 ! indicatorOfTypeOfLevel 234 isec1(8)=0 235 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) & 236 .and.(discipl.eq.2)) then ! LSM 237 isec1(6)=81 ! indicatorOfParameter 238 isec1(7)=1 ! indicatorOfTypeOfLevel 239 isec1(8)=0 240 elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH 241 isec1(6)=221 ! indicatorOfParameter 242 isec1(7)=1 ! indicatorOfTypeOfLevel 243 isec1(8)=0 244 elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP 245 isec1(6)=62 ! indicatorOfParameter 246 isec1(7)=1 ! indicatorOfTypeOfLevel 247 isec1(8)=0 248 elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP 249 isec1(6)=63 ! indicatorOfParameter 250 isec1(7)=1 ! indicatorOfTypeOfLevel 251 isec1(8)=0 252 endif 152 !read the grib2 identifiers 153 154 call grib_get_string(igrib,'shortName',shortname,iret) 155 156 call grib_get_int(igrib,'discipline',discipl,iret) 157 call grib_check(iret,gribFunction,gribErrorMsg) 158 call grib_get_int(igrib,'parameterCategory',parCat,iret) 159 call grib_check(iret,gribFunction,gribErrorMsg) 160 call grib_get_int(igrib,'parameterNumber',parNum,iret) 161 call grib_check(iret,gribFunction,gribErrorMsg) 162 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 163 call grib_check(iret,gribFunction,gribErrorMsg) 164 ! 165 call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & 166 valSurf,iret) 167 call grib_check(iret,gribFunction,gribErrorMsg) 168 169 ! write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname 170 !convert to grib1 identifiers 171 isec1(6:8)=-1 172 xsec18 =-1.0 173 !JMA / SH: isec1(8) not evaluated any more below 174 !b/c with GRIB 2 this may be a real variable 175 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T 176 isec1(6)=11 ! indicatorOfParameter 177 isec1(7)=100 ! indicatorOfTypeOfLevel 178 xsec18=valSurf/100.0 ! level, convert to hPa 179 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U 180 isec1(6)=33 ! indicatorOfParameter 181 isec1(7)=100 ! indicatorOfTypeOfLevel 182 xsec18=valSurf/100.0 ! level, convert to hPa 183 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V 184 isec1(6)=34 ! indicatorOfParameter 185 isec1(7)=100 ! indicatorOfTypeOfLevel 186 xsec18=valSurf/100.0 ! level, convert to hPa 187 elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W 188 isec1(6)=39 ! indicatorOfParameter 189 isec1(7)=100 ! indicatorOfTypeOfLevel 190 xsec18=valSurf/100.0 ! level, convert to hPa 191 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH 192 isec1(6)=52 ! indicatorOfParameter 193 isec1(7)=100 ! indicatorOfTypeOfLevel 194 xsec18=valSurf/100.0 ! level, convert to hPa 195 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2 196 isec1(6)=52 ! indicatorOfParameter 197 isec1(7)=105 ! indicatorOfTypeOfLevel 198 xsec18=real(2) 199 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2 200 isec1(6)=11 ! indicatorOfParameter 201 isec1(7)=105 ! indicatorOfTypeOfLevel 202 xsec18=real(2) 203 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10 204 isec1(6)=33 ! indicatorOfParameter 205 isec1(7)=105 ! indicatorOfTypeOfLevel 206 xsec18=real(10) 207 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10 208 isec1(6)=34 ! indicatorOfParameter 209 isec1(7)=105 ! indicatorOfTypeOfLevel 210 xsec18=real(10) 211 elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]: 212 isec1(6)=153 ! indicatorOfParameter 213 isec1(7)=100 ! indicatorOfTypeOfLevel 214 xsec18=valSurf/100.0 ! level, convert to hPa 215 elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP 216 isec1(6)=2 ! indicatorOfParameter 217 isec1(7)=102 ! indicatorOfTypeOfLevel 218 xsec18=real(0) 219 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP 220 isec1(6)=1 ! indicatorOfParameter 221 isec1(7)=1 ! indicatorOfTypeOfLevel 222 xsec18=real(0) 223 elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW 224 isec1(6)=66 ! indicatorOfParameter 225 isec1(7)=1 ! indicatorOfTypeOfLevel 226 xsec18=real(0) 227 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0 228 isec1(6)=11 ! indicatorOfParameter 229 isec1(7)=107 ! indicatorOfTypeOfLevel 230 xsec18=0.995 ! lowest sigma level 231 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0 232 isec1(6)=33 ! indicatorOfParameter 233 isec1(7)=107 ! indicatorOfTypeOfLevel 234 xsec18=0.995 ! lowest sigma level 235 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0 236 isec1(6)=34 ! indicatorOfParameter 237 isec1(7)=107 ! indicatorOfTypeOfLevel 238 xsec18=0.995 ! lowest sigma level 239 elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO 240 isec1(6)=7 ! indicatorOfParameter 241 isec1(7)=1 ! indicatorOfTypeOfLevel 242 xsec18=real(0) 243 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) & 244 .and.(discipl.eq.2)) then ! LSM 245 isec1(6)=81 ! indicatorOfParameter 246 isec1(7)=1 ! indicatorOfTypeOfLevel 247 xsec18=real(0) 248 elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH 249 isec1(6)=221 ! indicatorOfParameter 250 isec1(7)=1 ! indicatorOfTypeOfLevel 251 xsec18=real(0) 252 elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP 253 isec1(6)=62 ! indicatorOfParameter 254 isec1(7)=1 ! indicatorOfTypeOfLevel 255 xsec18=real(0) 256 elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP 257 isec1(6)=63 ! indicatorOfParameter 258 isec1(7)=1 ! indicatorOfTypeOfLevel 259 xsec18=real(0) 260 endif 253 261 254 262 endif ! gribVer 255 263 256 264 if (isec1(6).ne.-1) then 257 265 ! get the size and data of the values array 258 266 call grib_get_real4_array(igrib,'values',zsec4,iret) 259 !call grib_check(iret,gribFunction,gribErrorMsg)267 call grib_check(iret,gribFunction,gribErrorMsg) 260 268 endif 261 269 262 270 if(ifield.eq.1) then 263 271 264 265 266 call grib_get_int(igrib,'numberOfPointsAlongAParallel', &267 isec2(2),iret)268 !call grib_check(iret,gribFunction,gribErrorMsg)269 call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &270 isec2(3),iret)271 !call grib_check(iret,gribFunction,gribErrorMsg)272 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &273 xauxin,iret)274 !call grib_check(iret,gribFunction,gribErrorMsg)275 call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &276 yauxin,iret)277 !call grib_check(iret,gribFunction,gribErrorMsg)278 xaux=xauxin+real(nxshift)*dx279 yaux=yauxin280 281 272 !get the required fields from section 2 273 !store compatible to gribex input 274 call grib_get_int(igrib,'numberOfPointsAlongAParallel', & 275 isec2(2),iret) 276 call grib_check(iret,gribFunction,gribErrorMsg) 277 call grib_get_int(igrib,'numberOfPointsAlongAMeridian', & 278 isec2(3),iret) 279 call grib_check(iret,gribFunction,gribErrorMsg) 280 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & 281 xauxin,iret) 282 call grib_check(iret,gribFunction,gribErrorMsg) 283 call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', & 284 yauxin,iret) 285 call grib_check(iret,gribFunction,gribErrorMsg) 286 xaux=xauxin+real(nxshift)*dx 287 yaux=yauxin 288 289 ! CHECK GRID SPECIFICATIONS 282 290 283 291 if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT' 284 292 if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT' 285 if(xaux.eq.0.) xaux=-1 79.0 ! NCEP DATA293 if(xaux.eq.0.) xaux=-180.0 ! NCEP DATA 286 294 xaux0=xlon0 287 295 yaux0=ylat0 … … 290 298 if(xaux0.lt.0.) xaux0=xaux0+360. 291 299 if(yaux0.lt.0.) yaux0=yaux0+360. 292 if(abs(xaux-xaux0).gt.eps) & 293 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' 294 if(abs(yaux-yaux0).gt.eps) & 295 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' 296 endif 297 !HSO end of edits 298 299 i179=nint(179./dx) 300 if (dx.lt.0.7) then 301 i180=nint(180./dx)+1 ! 0.5 deg data 302 else 303 i180=nint(179./dx)+1 ! 1 deg data 304 endif 305 i181=i180+1 300 if(abs(xaux-xaux0).gt.eps) then 301 write (*, *) xaux, xaux0 302 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' 303 endif 304 if(abs(yaux-yaux0).gt.eps) then 305 write (*, *) yaux, yaux0 306 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' 307 end if 308 endif 309 !HSO end of edits 310 311 i180=nint(180./dx) 306 312 307 313 if (isec1(6).ne.-1) then 308 314 309 do j=0,nymin1 310 do i=0,nxfield-1 311 if((isec1(6).eq.011).and.(isec1(7).eq.100)) then 312 ! TEMPERATURE 313 if((i.eq.0).and.(j.eq.0)) then 314 do ii=1,nuvz 315 if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii 316 end do 317 endif 318 help=zsec4(nxfield*(ny-j-1)+i+1) 319 if(i.le.i180) then 320 tth(i179+i,j,numpt,n)=help 321 else 322 tth(i-i181,j,numpt,n)=help 323 endif 324 endif 325 if((isec1(6).eq.033).and.(isec1(7).eq.100)) then 326 ! U VELOCITY 327 if((i.eq.0).and.(j.eq.0)) then 328 do ii=1,nuvz 329 if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii 330 end do 331 endif 332 help=zsec4(nxfield*(ny-j-1)+i+1) 333 if(i.le.i180) then 334 uuh(i179+i,j,numpu)=help 335 else 336 uuh(i-i181,j,numpu)=help 337 endif 338 endif 339 if((isec1(6).eq.034).and.(isec1(7).eq.100)) then 340 ! V VELOCITY 341 if((i.eq.0).and.(j.eq.0)) then 342 do ii=1,nuvz 343 if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii 344 end do 345 endif 346 help=zsec4(nxfield*(ny-j-1)+i+1) 347 if(i.le.i180) then 348 vvh(i179+i,j,numpv)=help 349 else 350 vvh(i-i181,j,numpv)=help 351 endif 352 endif 353 if((isec1(6).eq.052).and.(isec1(7).eq.100)) then 354 ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER 355 if((i.eq.0).and.(j.eq.0)) then 356 do ii=1,nuvz 357 if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii 358 end do 359 endif 360 help=zsec4(nxfield*(ny-j-1)+i+1) 361 if(i.le.i180) then 362 qvh(i179+i,j,numprh,n)=help 363 else 364 qvh(i-i181,j,numprh,n)=help 365 endif 366 endif 367 if((isec1(6).eq.001).and.(isec1(7).eq.001)) then 368 ! SURFACE PRESSURE 369 help=zsec4(nxfield*(ny-j-1)+i+1) 370 if(i.le.i180) then 371 ps(i179+i,j,1,n)=help 372 else 373 ps(i-i181,j,1,n)=help 374 endif 375 endif 376 if((isec1(6).eq.039).and.(isec1(7).eq.100)) then 377 ! W VELOCITY 378 if((i.eq.0).and.(j.eq.0)) then 379 do ii=1,nuvz 380 if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii 381 end do 382 endif 383 help=zsec4(nxfield*(ny-j-1)+i+1) 384 if(i.le.i180) then 385 wwh(i179+i,j,numpw)=help 386 else 387 wwh(i-i181,j,numpw)=help 388 endif 389 endif 390 if((isec1(6).eq.066).and.(isec1(7).eq.001)) then 391 ! SNOW DEPTH 392 help=zsec4(nxfield*(ny-j-1)+i+1) 393 if(i.le.i180) then 394 sd(i179+i,j,1,n)=help 395 else 396 sd(i-i181,j,1,n)=help 397 endif 398 endif 399 if((isec1(6).eq.002).and.(isec1(7).eq.102)) then 400 ! MEAN SEA LEVEL PRESSURE 401 help=zsec4(nxfield*(ny-j-1)+i+1) 402 if(i.le.i180) then 403 msl(i179+i,j,1,n)=help 404 else 405 msl(i-i181,j,1,n)=help 406 endif 407 endif 408 if((isec1(6).eq.071).and.(isec1(7).eq.244)) then 409 ! TOTAL CLOUD COVER 410 help=zsec4(nxfield*(ny-j-1)+i+1) 411 if(i.le.i180) then 412 tcc(i179+i,j,1,n)=help 413 else 414 tcc(i-i181,j,1,n)=help 415 endif 416 endif 417 if((isec1(6).eq.033).and.(isec1(7).eq.105).and. & 418 (isec1(8).eq.10)) then 419 ! 10 M U VELOCITY 420 help=zsec4(nxfield*(ny-j-1)+i+1) 421 if(i.le.i180) then 422 u10(i179+i,j,1,n)=help 423 else 424 u10(i-i181,j,1,n)=help 425 endif 426 endif 427 if((isec1(6).eq.034).and.(isec1(7).eq.105).and. & 428 (isec1(8).eq.10)) then 429 ! 10 M V VELOCITY 430 help=zsec4(nxfield*(ny-j-1)+i+1) 431 if(i.le.i180) then 432 v10(i179+i,j,1,n)=help 433 else 434 v10(i-i181,j,1,n)=help 435 endif 436 endif 437 if((isec1(6).eq.011).and.(isec1(7).eq.105).and. & 438 (isec1(8).eq.02)) then 439 ! 2 M TEMPERATURE 440 help=zsec4(nxfield*(ny-j-1)+i+1) 441 if(i.le.i180) then 442 tt2(i179+i,j,1,n)=help 443 else 444 tt2(i-i181,j,1,n)=help 445 endif 446 endif 447 if((isec1(6).eq.017).and.(isec1(7).eq.105).and. & 448 (isec1(8).eq.02)) then 449 ! 2 M DEW POINT TEMPERATURE 450 help=zsec4(nxfield*(ny-j-1)+i+1) 451 if(i.le.i180) then 452 td2(i179+i,j,1,n)=help 453 else 454 td2(i-i181,j,1,n)=help 455 endif 456 endif 457 if((isec1(6).eq.062).and.(isec1(7).eq.001)) then 458 ! LARGE SCALE PREC. 459 help=zsec4(nxfield*(ny-j-1)+i+1) 460 if(i.le.i180) then 461 lsprec(i179+i,j,1,n)=help 462 else 463 lsprec(i-i181,j,1,n)=help 464 endif 465 endif 466 if((isec1(6).eq.063).and.(isec1(7).eq.001)) then 467 ! CONVECTIVE PREC. 468 help=zsec4(nxfield*(ny-j-1)+i+1) 469 if(i.le.i180) then 470 convprec(i179+i,j,1,n)=help 471 else 472 convprec(i-i181,j,1,n)=help 473 endif 474 endif 475 if((isec1(6).eq.007).and.(isec1(7).eq.001)) then 476 ! TOPOGRAPHY 477 help=zsec4(nxfield*(ny-j-1)+i+1) 478 if(i.le.i180) then 479 oro(i179+i,j)=help 480 excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 481 else 482 oro(i-i181,j)=help 483 excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 484 endif 485 endif 486 if((isec1(6).eq.081).and.(isec1(7).eq.001)) then 487 ! LAND SEA MASK 488 help=zsec4(nxfield*(ny-j-1)+i+1) 489 if(i.le.i180) then 490 lsm(i179+i,j)=help 491 else 492 lsm(i-i181,j)=help 493 endif 494 endif 495 if((isec1(6).eq.221).and.(isec1(7).eq.001)) then 496 ! MIXING HEIGHT 497 help=zsec4(nxfield*(ny-j-1)+i+1) 498 if(i.le.i180) then 499 hmix(i179+i,j,1,n)=help 500 else 501 hmix(i-i181,j,1,n)=help 502 endif 503 endif 504 if((isec1(6).eq.052).and.(isec1(7).eq.105).and. & 505 (isec1(8).eq.02)) then 506 ! 2 M RELATIVE HUMIDITY 507 help=zsec4(nxfield*(ny-j-1)+i+1) 508 if(i.le.i180) then 509 qvh2(i179+i,j)=help 510 else 511 qvh2(i-i181,j)=help 512 endif 513 endif 514 if((isec1(6).eq.011).and.(isec1(7).eq.107)) then 515 ! TEMPERATURE LOWEST SIGMA LEVEL 516 help=zsec4(nxfield*(ny-j-1)+i+1) 517 if(i.le.i180) then 518 tlev1(i179+i,j)=help 519 else 520 tlev1(i-i181,j)=help 521 endif 522 endif 523 if((isec1(6).eq.033).and.(isec1(7).eq.107)) then 524 ! U VELOCITY LOWEST SIGMA LEVEL 525 help=zsec4(nxfield*(ny-j-1)+i+1) 526 if(i.le.i180) then 527 ulev1(i179+i,j)=help 528 else 529 ulev1(i-i181,j)=help 530 endif 531 endif 532 if((isec1(6).eq.034).and.(isec1(7).eq.107)) then 533 ! V VELOCITY LOWEST SIGMA LEVEL 534 help=zsec4(nxfield*(ny-j-1)+i+1) 535 if(i.le.i180) then 536 vlev1(i179+i,j)=help 537 else 538 vlev1(i-i181,j)=help 539 endif 540 endif 315 ! write (*, *) 'nxfield: ', nxfield, i180 316 317 do j=0,nymin1 318 do i=0,nxfield-1 319 if((isec1(6).eq.011).and.(isec1(7).eq.100)) then 320 ! TEMPERATURE 321 if((i.eq.0).and.(j.eq.0)) then 322 numpt=minloc(abs(xsec18*100.0-akz),dim=1) 323 endif 324 help=zsec4(nxfield*(ny-j-1)+i+1) 325 if (help.eq.0) then 326 write (*, *) 'i, j: ', i, j 327 stop 'help == 0.0' 328 endif 329 if(i.lt.i180) then 330 tth(i180+i,j,numpt,n)=help 331 else 332 tth(i-i180,j,numpt,n)=help 333 endif 334 endif 335 if((isec1(6).eq.033).and.(isec1(7).eq.100)) then 336 ! U VELOCITY 337 if((i.eq.0).and.(j.eq.0)) then 338 numpu=minloc(abs(xsec18*100.0-akz),dim=1) 339 endif 340 help=zsec4(nxfield*(ny-j-1)+i+1) 341 if(i.lt.i180) then 342 uuh(i180+i,j,numpu)=help 343 else 344 uuh(i-i180,j,numpu)=help 345 endif 346 endif 347 if((isec1(6).eq.034).and.(isec1(7).eq.100)) then 348 ! V VELOCITY 349 if((i.eq.0).and.(j.eq.0)) then 350 numpv=minloc(abs(xsec18*100.0-akz),dim=1) 351 endif 352 help=zsec4(nxfield*(ny-j-1)+i+1) 353 if(i.lt.i180) then 354 vvh(i180+i,j,numpv)=help 355 else 356 vvh(i-i180,j,numpv)=help 357 endif 358 endif 359 if((isec1(6).eq.052).and.(isec1(7).eq.100)) then 360 ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER 361 if((i.eq.0).and.(j.eq.0)) then 362 numprh=minloc(abs(xsec18*100.0-akz),dim=1) 363 endif 364 help=zsec4(nxfield*(ny-j-1)+i+1) 365 if(i.lt.i180) then 366 qvh(i180+i,j,numprh,n)=help 367 else 368 qvh(i-i180,j,numprh,n)=help 369 endif 370 endif 371 if((isec1(6).eq.001).and.(isec1(7).eq.001)) then 372 ! SURFACE PRESSURE 373 help=zsec4(nxfield*(ny-j-1)+i+1) 374 if(i.lt.i180) then 375 ps(i180+i,j,1,n)=help 376 else 377 ps(i-i180,j,1,n)=help 378 endif 379 endif 380 if((isec1(6).eq.039).and.(isec1(7).eq.100)) then 381 ! W VELOCITY 382 if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1) 383 help=zsec4(nxfield*(ny-j-1)+i+1) 384 if(i.lt.i180) then 385 wwh(i180+i,j,numpw)=help 386 else 387 wwh(i-i180,j,numpw)=help 388 endif 389 endif 390 if((isec1(6).eq.066).and.(isec1(7).eq.001)) then 391 ! SNOW DEPTH 392 help=zsec4(nxfield*(ny-j-1)+i+1) 393 if(i.lt.i180) then 394 sd(i180+i,j,1,n)=help 395 else 396 sd(i-i180,j,1,n)=help 397 endif 398 endif 399 if((isec1(6).eq.002).and.(isec1(7).eq.102)) then 400 ! MEAN SEA LEVEL PRESSURE 401 help=zsec4(nxfield*(ny-j-1)+i+1) 402 if(i.lt.i180) then 403 msl(i180+i,j,1,n)=help 404 else 405 msl(i-i180,j,1,n)=help 406 endif 407 endif 408 if((isec1(6).eq.071).and.(isec1(7).eq.244)) then 409 ! TOTAL CLOUD COVER 410 help=zsec4(nxfield*(ny-j-1)+i+1) 411 if(i.lt.i180) then 412 tcc(i180+i,j,1,n)=help 413 else 414 tcc(i-i180,j,1,n)=help 415 endif 416 endif 417 if((isec1(6).eq.033).and.(isec1(7).eq.105).and. & 418 (nint(xsec18).eq.10)) then 419 ! 10 M U VELOCITY 420 help=zsec4(nxfield*(ny-j-1)+i+1) 421 if(i.lt.i180) then 422 u10(i180+i,j,1,n)=help 423 else 424 u10(i-i180,j,1,n)=help 425 endif 426 endif 427 if((isec1(6).eq.034).and.(isec1(7).eq.105).and. & 428 (nint(xsec18).eq.10)) then 429 ! 10 M V VELOCITY 430 help=zsec4(nxfield*(ny-j-1)+i+1) 431 if(i.lt.i180) then 432 v10(i180+i,j,1,n)=help 433 else 434 v10(i-i180,j,1,n)=help 435 endif 436 endif 437 if((isec1(6).eq.011).and.(isec1(7).eq.105).and. & 438 (nint(xsec18).eq.2)) then 439 ! 2 M TEMPERATURE 440 help=zsec4(nxfield*(ny-j-1)+i+1) 441 if(i.lt.i180) then 442 tt2(i180+i,j,1,n)=help 443 else 444 tt2(i-i180,j,1,n)=help 445 endif 446 endif 447 if((isec1(6).eq.017).and.(isec1(7).eq.105).and. & 448 (nint(xsec18).eq.2)) then 449 ! 2 M DEW POINT TEMPERATURE 450 help=zsec4(nxfield*(ny-j-1)+i+1) 451 if(i.lt.i180) then 452 td2(i180+i,j,1,n)=help 453 else 454 td2(i-i180,j,1,n)=help 455 endif 456 endif 457 if((isec1(6).eq.062).and.(isec1(7).eq.001)) then 458 ! LARGE SCALE PREC. 459 help=zsec4(nxfield*(ny-j-1)+i+1) 460 if(i.lt.i180) then 461 lsprec(i180+i,j,1,n)=help 462 else 463 lsprec(i-i180,j,1,n)=help 464 endif 465 endif 466 if((isec1(6).eq.063).and.(isec1(7).eq.001)) then 467 ! CONVECTIVE PREC. 468 help=zsec4(nxfield*(ny-j-1)+i+1) 469 if(i.lt.i180) then 470 convprec(i180+i,j,1,n)=help 471 else 472 convprec(i-i180,j,1,n)=help 473 endif 474 endif 475 if((isec1(6).eq.007).and.(isec1(7).eq.001)) then 476 ! TOPOGRAPHY 477 help=zsec4(nxfield*(ny-j-1)+i+1) 478 if(i.lt.i180) then 479 oro(i180+i,j)=help 480 excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 481 else 482 oro(i-i180,j)=help 483 excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 484 endif 485 endif 486 if((isec1(6).eq.081).and.(isec1(7).eq.001)) then 487 ! LAND SEA MASK 488 help=zsec4(nxfield*(ny-j-1)+i+1) 489 if(i.lt.i180) then 490 lsm(i180+i,j)=help 491 else 492 lsm(i-i180,j)=help 493 endif 494 endif 495 if((isec1(6).eq.221).and.(isec1(7).eq.001)) then 496 ! MIXING HEIGHT 497 help=zsec4(nxfield*(ny-j-1)+i+1) 498 if(i.lt.i180) then 499 hmix(i180+i,j,1,n)=help 500 else 501 hmix(i-i180,j,1,n)=help 502 endif 503 endif 504 if((isec1(6).eq.052).and.(isec1(7).eq.105).and. & 505 (nint(xsec18).eq.02)) then 506 ! 2 M RELATIVE HUMIDITY 507 help=zsec4(nxfield*(ny-j-1)+i+1) 508 if(i.lt.i180) then 509 qvh2(i180+i,j)=help 510 else 511 qvh2(i-i180,j)=help 512 endif 513 endif 514 if((isec1(6).eq.011).and.(isec1(7).eq.107)) then 515 ! TEMPERATURE LOWEST SIGMA LEVEL 516 help=zsec4(nxfield*(ny-j-1)+i+1) 517 if(i.lt.i180) then 518 tlev1(i180+i,j)=help 519 else 520 tlev1(i-i180,j)=help 521 endif 522 endif 523 if((isec1(6).eq.033).and.(isec1(7).eq.107)) then 524 ! U VELOCITY LOWEST SIGMA LEVEL 525 help=zsec4(nxfield*(ny-j-1)+i+1) 526 if(i.lt.i180) then 527 ulev1(i180+i,j)=help 528 else 529 ulev1(i-i180,j)=help 530 endif 531 endif 532 if((isec1(6).eq.034).and.(isec1(7).eq.107)) then 533 ! V VELOCITY LOWEST SIGMA LEVEL 534 help=zsec4(nxfield*(ny-j-1)+i+1) 535 if(i.lt.i180) then 536 vlev1(i180+i,j)=help 537 else 538 vlev1(i-i180,j)=help 539 endif 540 endif 541 541 ! SEC & IP 12/2018 read GFS clouds 542 if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg] 543 if((i.eq.0).and.(j.eq.0)) then 544 do ii=1,nuvz 545 if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii 546 end do 547 endif 548 help=zsec4(nxfield*(ny-j-1)+i+1) 549 if(i.le.i180) then 550 clwch(i179+i,j,numpclwch,n)=help 551 else 552 clwch(i-i181,j,numpclwch,n)=help 553 endif 554 readclouds=.true. 555 sumclouds=.true. 556 ! readclouds=.false. 557 ! sumclouds=.false. 558 endif 559 560 542 if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg] 543 if((i.eq.0).and.(j.eq.0)) then 544 numpclwch=minloc(abs(xsec18*100.0-akz),dim=1) 545 endif 546 help=zsec4(nxfield*(ny-j-1)+i+1) 547 if(i.lt.i180) then 548 clwch(i180+i,j,numpclwch,n)=help 549 else 550 clwch(i-i180,j,numpclwch,n)=help 551 endif 552 readclouds=.true. 553 sumclouds=.true. 554 endif 555 556 557 end do 561 558 end do 562 end do563 559 564 560 endif 565 561 566 562 if((isec1(6).eq.33).and.(isec1(7).eq.100)) then 567 563 ! NCEP ISOBARIC LEVELS 568 564 iumax=iumax+1 569 565 endif … … 571 567 call grib_release(igrib) 572 568 goto 10 !! READ NEXT LEVEL OR PARAMETER 573 574 575 576 577 578 50 569 ! 570 ! CLOSING OF INPUT DATA FILE 571 ! 572 573 !HSO close grib file 574 50 continue 579 575 call grib_close_file(ifile) 580 576 581 577 ! SENS. HEAT FLUX 582 578 sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files 583 579 hflswitch=.false. ! Heat flux not available 584 580 ! SOLAR RADIATIVE FLUXES 585 581 ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files 586 582 ! EW SURFACE STRESS 587 583 ewss=0.0 ! not available from gfs.tccz.pgrbfxx files 588 584 ! NS SURFACE STRESS 589 585 nsss=0.0 ! not available from gfs.tccz.pgrbfxx files 590 586 strswitch=.false. ! stress not available 591 587 592 588 ! CONVERT TP TO LSP (GRIB2 only) 593 589 if (gribVer.eq.2) then 594 do j=0,nymin1 595 do i=0,nxfield-1 596 if(i.le.i180) then 597 if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur 598 lsprec(i179+i,j,1,n)= & 599 lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n) 600 else 601 lsprec(i179+i,j,1,n)=0 602 endif 603 else 604 if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then 605 lsprec(i-i181,j,1,n)= & 606 lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n) 607 else 608 lsprec(i-i181,j,1,n)=0 609 endif 610 endif 611 enddo 612 enddo 613 endif 614 !HSO end edits 615 616 617 ! TRANSFORM RH TO SPECIFIC HUMIDITY 590 lsprec(0:nxfield-1, 0:nymin1, 1, n) = max( 0.0 , lsprec(0:nxfield-1,0:nymin1,1,n)-convprec(0:nxfield-1,0:nymin1,1,n) ) 591 endif 592 !HSO end edits 593 594 595 ! TRANSFORM RH TO SPECIFIC HUMIDITY 618 596 619 597 do j=0,ny-1 … … 622 600 help=qvh(i,j,k,n) 623 601 temp=tth(i,j,k,n) 602 if (temp .eq. 0.0) then 603 write (*, *) i, j, k, n 604 temp = 273.0 605 endif 624 606 plev1=akm(k)+bkm(k)*ps(i,j,1,n) 625 607 elev=ew(temp)*help/100.0 … … 629 611 end do 630 612 631 632 633 613 ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY 614 ! USING BOLTON'S (1980) FORMULA 615 ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA 634 616 635 617 do j=0,ny-1 636 618 do i=0,nxfield-1 637 help=qvh2(i,j) 638 temp=tt2(i,j,1,n) 639 elev=ew(temp)/100.*help/100. !vapour pressure in hPa 640 td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273. 641 if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n) 619 help=qvh2(i,j) 620 temp=tt2(i,j,1,n) 621 if (temp .eq. 0.0) then 622 write (*, *) i, j, n 623 temp = 273.0 624 endif 625 elev=ew(temp)/100.*help/100. !vapour pressure in hPa 626 td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273. 627 if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n) 642 628 end do 643 629 end do … … 653 639 654 640 655 656 657 641 ! For global fields, assign the leftmost data column also to the rightmost 642 ! data column; if required, shift whole grid by nxshift grid points 643 !************************************************************************* 658 644 659 645 if (xglobal) then … … 691 677 do i=0,nxmin1 692 678 do j=0,nymin1 693 679 ! Convert precip. from mm/s -> mm/hour 694 680 convprec(i,j,1,n)=convprec(i,j,1,n)*3600. 695 681 lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600. … … 699 685 700 686 if ((.not.hflswitch).or.(.not.strswitch)) then 701 702 703 704 705 687 ! write(*,*) 'WARNING: No flux data contained in GRIB file ', 688 ! + wfname(indj) 689 690 ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD 691 !*************************************************************************** 706 692 707 693 do i=0,nxmin1 … … 723 709 724 710 return 725 888 711 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 726 712 write(*,*) ' #### ',wfname(indj),' #### ' 727 713 write(*,*) ' #### IS NOT GRIB FORMAT !!! #### ' 728 714 stop 'Execution terminated' 729 999 715 999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 730 716 write(*,*) ' #### ',wfname(indj),' #### ' 731 717 write(*,*) ' #### CANNOT BE OPENED !!! #### ' -
src/richardson.f90
- Property mode changed from 100644 to 100755
r92fab65 r467460a 141 141 zold=z 142 142 end do 143 k=k-1 ! ESO: make sure k <= nuvz (ticket #139)143 ! k=k-1 ! ESO: make sure k <= nuvz (ticket #139) 144 144 20 continue 145 145 146 146 ! Determine Richardson number between the critical levels 147 147 !******************************************************** 148 ! JMA: It may happen that k >= nuvz: 149 ! JMA: In that case: k is set to nuvz 150 151 if (k .gt. nuvz) k = nuvz ! JMA 148 152 149 153 zl1=zold
Note: See TracChangeset
for help on using the changeset viewer.