Changeset 9ca6e38 in flexpart.git
- Timestamp:
- Mar 15, 2021, 9:48:30 AM (3 years ago)
- Branches:
- GFS_025, dev
- Children:
- aa939a9
- Parents:
- 467460a
- Location:
- src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
src/calcmatrix.f90
r467460a r9ca6e38 67 67 68 68 phconv(1) = psconv 69 70 do kuvz = 2,nuvz71 k = kuvz-172 if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then69 ! 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 73 pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv) 74 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) ) 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/readwind_gfs.f90
r467460a 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 … … 74 74 integer :: i180 75 75 76 77 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) … … 88 88 logical :: hflswitch,strswitch 89 89 90 90 !HSO for grib api error messages 91 91 character(len=24) :: gribErrorMsg = 'Error reading grib file' 92 92 character(len=20) :: gribFunction = 'readwind_gfs' … … 101 101 102 102 103 104 105 103 ! OPENING OF DATA FILE (GRIB CODE) 104 105 !HSO 106 106 call grib_open_file(ifile,path(3)(1:length(3)) & 107 107 //trim(wfname(indj)),'r',iret) 108 108 if (iret.ne.GRIB_SUCCESS) then 109 109 goto 888 ! ERROR DETECTED 110 110 endif 111 111 !turn on support for multi fields messages 112 112 call grib_multi_support_on 113 113 … … 119 119 numpclwch=0 120 120 ifield=0 121 10 122 123 124 121 10 ifield=ifield+1 122 ! 123 ! GET NEXT FIELDS 124 ! 125 125 call grib_new_from_file(ifile,igrib,iret) 126 126 if (iret.eq.GRIB_END_OF_FILE) then … … 130 130 endif 131 131 132 132 !first see if we read GRIB1 or GRIB2 133 133 call grib_get_int(igrib,'editionNumber',gribVer,iret) 134 134 ! call grib_check(iret,gribFunction,gribErrorMsg) … … 136 136 if (gribVer.eq.1) then ! GRIB Edition 1 137 137 138 !read the grib1 identifiers 139 call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) 140 call grib_check(iret,gribFunction,gribErrorMsg) 141 call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret) 142 call grib_check(iret,gribFunction,gribErrorMsg) 143 call grib_get_int(igrib,'level',isec1(8),iret) 144 call grib_check(iret,gribFunction,gribErrorMsg) 145 !JMA / SH: isec1(8) not evaluated any more below 146 !b/c with GRIB 2 this may be a real variable 147 xsec18 = real(isec1(8)) 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)) 148 149 149 150 else ! GRIB Edition 2 150 151 151 !read the grib2 identifiers 152 call grib_get_string(igrib,'shortName',shortname,iret) 153 154 call grib_get_int(igrib,'discipline',discipl,iret) 155 call grib_check(iret,gribFunction,gribErrorMsg) 156 call grib_get_int(igrib,'parameterCategory',parCat,iret) 157 call grib_check(iret,gribFunction,gribErrorMsg) 158 call grib_get_int(igrib,'parameterNumber',parNum,iret) 159 call grib_check(iret,gribFunction,gribErrorMsg) 160 call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) 161 call grib_check(iret,gribFunction,gribErrorMsg) 162 call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & 163 valSurf,iret) 164 call grib_check(iret,gribFunction,gribErrorMsg) 165 166 ! write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname 167 !convert to grib1 identifiers 168 isec1(6:8)=-1 169 xsec18 =-1.0 170 !JMA / SH: isec1(8) not evaluated any more below 171 !b/c with GRIB 2 this may be a real variable 172 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T 173 isec1(6)=11 ! indicatorOfParameter 174 isec1(7)=100 ! indicatorOfTypeOfLevel 175 xsec18=valSurf/100.0 ! level, convert to hPa 176 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U 177 isec1(6)=33 ! indicatorOfParameter 178 isec1(7)=100 ! indicatorOfTypeOfLevel 179 xsec18=valSurf/100.0 ! level, convert to hPa 180 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V 181 isec1(6)=34 ! indicatorOfParameter 182 isec1(7)=100 ! indicatorOfTypeOfLevel 183 xsec18=valSurf/100.0 ! level, convert to hPa 184 elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W 185 isec1(6)=39 ! indicatorOfParameter 186 isec1(7)=100 ! indicatorOfTypeOfLevel 187 xsec18=valSurf/100.0 ! level, convert to hPa 188 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH 189 isec1(6)=52 ! indicatorOfParameter 190 isec1(7)=100 ! indicatorOfTypeOfLevel 191 xsec18=valSurf/100.0 ! level, convert to hPa 192 elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2 193 isec1(6)=52 ! indicatorOfParameter 194 isec1(7)=105 ! indicatorOfTypeOfLevel 195 xsec18=real(2) 196 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2 197 isec1(6)=11 ! indicatorOfParameter 198 isec1(7)=105 ! indicatorOfTypeOfLevel 199 xsec18=real(2) 200 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10 201 isec1(6)=33 ! indicatorOfParameter 202 isec1(7)=105 ! indicatorOfTypeOfLevel 203 xsec18=real(10) 204 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10 205 isec1(6)=34 ! indicatorOfParameter 206 isec1(7)=105 ! indicatorOfTypeOfLevel 207 xsec18=real(10) 208 elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]: 209 isec1(6)=153 ! indicatorOfParameter 210 isec1(7)=100 ! indicatorOfTypeOfLevel 211 xsec18=valSurf/100.0 ! level, convert to hPa 212 elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP 213 isec1(6)=2 ! indicatorOfParameter 214 isec1(7)=102 ! indicatorOfTypeOfLevel 215 xsec18=real(0) 216 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP 217 isec1(6)=1 ! indicatorOfParameter 218 isec1(7)=1 ! indicatorOfTypeOfLevel 219 xsec18=real(0) 220 elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW 221 isec1(6)=66 ! indicatorOfParameter 222 isec1(7)=1 ! indicatorOfTypeOfLevel 223 xsec18=real(0) 224 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0 225 isec1(6)=11 ! indicatorOfParameter 226 isec1(7)=107 ! indicatorOfTypeOfLevel 227 xsec18=0.995 ! lowest sigma level 228 elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0 229 isec1(6)=33 ! indicatorOfParameter 230 isec1(7)=107 ! indicatorOfTypeOfLevel 231 xsec18=0.995 ! lowest sigma level 232 elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0 233 isec1(6)=34 ! indicatorOfParameter 234 isec1(7)=107 ! indicatorOfTypeOfLevel 235 xsec18=0.995 ! lowest sigma level 236 elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO 237 isec1(6)=7 ! indicatorOfParameter 238 isec1(7)=1 ! indicatorOfTypeOfLevel 239 xsec18=real(0) 240 elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) & 241 .and.(discipl.eq.2)) then ! LSM 242 isec1(6)=81 ! indicatorOfParameter 243 isec1(7)=1 ! indicatorOfTypeOfLevel 244 xsec18=real(0) 245 elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH 246 isec1(6)=221 ! indicatorOfParameter 247 isec1(7)=1 ! indicatorOfTypeOfLevel 248 xsec18=real(0) 249 elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP 250 isec1(6)=62 ! indicatorOfParameter 251 isec1(7)=1 ! indicatorOfTypeOfLevel 252 xsec18=real(0) 253 elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP 254 isec1(6)=63 ! indicatorOfParameter 255 isec1(7)=1 ! indicatorOfTypeOfLevel 256 xsec18=real(0) 257 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 258 261 259 262 endif ! gribVer 260 263 261 264 if (isec1(6).ne.-1) then 262 265 ! get the size and data of the values array 263 266 call grib_get_real4_array(igrib,'values',zsec4,iret) 264 267 call grib_check(iret,gribFunction,gribErrorMsg) … … 267 270 if(ifield.eq.1) then 268 271 269 270 271 call grib_get_int(igrib,'numberOfPointsAlongAParallel', &272 isec2(2),iret)273 call grib_check(iret,gribFunction,gribErrorMsg)274 call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &275 isec2(3),iret)276 call grib_check(iret,gribFunction,gribErrorMsg)277 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &278 xauxin,iret)279 call grib_check(iret,gribFunction,gribErrorMsg)280 call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &281 yauxin,iret)282 call grib_check(iret,gribFunction,gribErrorMsg)283 xaux=xauxin+real(nxshift)*dx284 yaux=yauxin285 286 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 287 290 288 291 if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT' … … 296 299 if(yaux0.lt.0.) yaux0=yaux0+360. 297 300 if(abs(xaux-xaux0).gt.eps) then 298 299 301 write (*, *) xaux, xaux0 302 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' 300 303 endif 301 304 if(abs(yaux-yaux0).gt.eps) then 302 303 305 write (*, *) yaux, yaux0 306 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' 304 307 end if 305 308 endif 306 309 !HSO end of edits 307 310 308 311 i180=nint(180./dx) … … 312 315 ! write (*, *) 'nxfield: ', nxfield, i180 313 316 314 do j=0,nymin1315 do i=0,nxfield-1316 if((isec1(6).eq.011).and.(isec1(7).eq.100)) then317 do j=0,nymin1 318 do i=0,nxfield-1 319 if((isec1(6).eq.011).and.(isec1(7).eq.100)) then 317 320 ! TEMPERATURE 318 if((i.eq.0).and.(j.eq.0)) then 319 do ii=1,nuvz+1 320 write(*,*) 'xsec18, akz(ii), xsec18*100.0-akz(ii), spacing(akz(ii), spacing(xsec18*100.0)',& 321 & xsec18, akz(ii), xsec18*100.0-akz(ii), spacing(akz(ii)), spacing(xsec18*100.0) 322 323 if (abs(xsec18*100.0-akz(ii)) < & 324 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then 325 numpt=ii 326 end if 327 end do 328 endif 329 help=zsec4(nxfield*(ny-j-1)+i+1) 330 if (help.eq.0) then 331 write (*, *) 'i, j: ', i, j 332 stop 'help == 0.0' 333 endif 334 if(i.lt.i180) then 335 ! ESO: dbg out-of-bounds issue 336 if (numpt < 1) then 337 write(*,*) 'nuvzmax, nuvz', nuvzmax, nuvz 338 write(*,*) 'i180+i,j,numpt,n', i180+i,j,numpt,n 339 end if 340 tth(i180+i,j,numpt,n)=help 341 else 342 tth(i-i180,j,numpt,n)=help 343 endif 344 endif 345 if((isec1(6).eq.033).and.(isec1(7).eq.100)) then 346 ! U VELOCITY 347 if((i.eq.0).and.(j.eq.0)) then 348 do ii=1,nuvz 349 if (abs(xsec18*100.0-akz(ii)) < & 350 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then 351 numpu=ii 352 end if 353 end do 354 endif 355 help=zsec4(nxfield*(ny-j-1)+i+1) 356 if(i.lt.i180) then 357 uuh(i180+i,j,numpu)=help 358 else 359 uuh(i-i180,j,numpu)=help 360 endif 361 endif 362 if((isec1(6).eq.034).and.(isec1(7).eq.100)) then 363 ! V VELOCITY 364 if((i.eq.0).and.(j.eq.0)) then 365 do ii=1,nuvz 366 if (abs(xsec18*100.0-akz(ii)) < & 367 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then 368 numpv=ii 369 end if 370 end do 371 endif 372 help=zsec4(nxfield*(ny-j-1)+i+1) 373 if(i.lt.i180) then 374 vvh(i180+i,j,numpv)=help 375 else 376 vvh(i-i180,j,numpv)=help 377 endif 378 endif 379 if((isec1(6).eq.052).and.(isec1(7).eq.100)) then 380 ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER 381 if((i.eq.0).and.(j.eq.0)) then 382 do ii=1,nuvz 383 if (abs(xsec18*100.0-akz(ii)) < & 384 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then 385 numprh=ii 386 end if 387 end do 388 endif 389 help=zsec4(nxfield*(ny-j-1)+i+1) 390 if(i.lt.i180) then 391 qvh(i180+i,j,numprh,n)=help 392 else 393 qvh(i-i180,j,numprh,n)=help 394 endif 395 endif 396 if((isec1(6).eq.001).and.(isec1(7).eq.001)) then 397 ! SURFACE PRESSURE 398 help=zsec4(nxfield*(ny-j-1)+i+1) 399 if(i.lt.i180) then 400 ps(i180+i,j,1,n)=help 401 else 402 ps(i-i180,j,1,n)=help 403 endif 404 endif 405 if((isec1(6).eq.039).and.(isec1(7).eq.100)) then 406 ! W VELOCITY 407 if((i.eq.0).and.(j.eq.0)) then 408 do ii=1,nuvz 409 if ((xsec18*100.0).eq.akz(ii)) numpw=ii 410 end do 411 endif 412 help=zsec4(nxfield*(ny-j-1)+i+1) 413 if(i.lt.i180) then 414 wwh(i180+i,j,numpw)=help 415 else 416 wwh(i-i180,j,numpw)=help 417 endif 418 endif 419 if((isec1(6).eq.066).and.(isec1(7).eq.001)) then 420 ! SNOW DEPTH 421 help=zsec4(nxfield*(ny-j-1)+i+1) 422 if(i.lt.i180) then 423 sd(i180+i,j,1,n)=help 424 else 425 sd(i-i180,j,1,n)=help 426 endif 427 endif 428 if((isec1(6).eq.002).and.(isec1(7).eq.102)) then 429 ! MEAN SEA LEVEL PRESSURE 430 help=zsec4(nxfield*(ny-j-1)+i+1) 431 if(i.lt.i180) then 432 msl(i180+i,j,1,n)=help 433 else 434 msl(i-i180,j,1,n)=help 435 endif 436 endif 437 if((isec1(6).eq.071).and.(isec1(7).eq.244)) then 438 ! TOTAL CLOUD COVER 439 help=zsec4(nxfield*(ny-j-1)+i+1) 440 if(i.lt.i180) then 441 tcc(i180+i,j,1,n)=help 442 else 443 tcc(i-i180,j,1,n)=help 444 endif 445 endif 446 if((isec1(6).eq.033).and.(isec1(7).eq.105).and. & 447 (nint(xsec18).eq.10)) then 448 ! 10 M U VELOCITY 449 help=zsec4(nxfield*(ny-j-1)+i+1) 450 if(i.lt.i180) then 451 u10(i180+i,j,1,n)=help 452 else 453 u10(i-i180,j,1,n)=help 454 endif 455 endif 456 if((isec1(6).eq.034).and.(isec1(7).eq.105).and. & 457 (nint(xsec18).eq.10)) then 458 ! 10 M V VELOCITY 459 help=zsec4(nxfield*(ny-j-1)+i+1) 460 if(i.lt.i180) then 461 v10(i180+i,j,1,n)=help 462 else 463 v10(i-i180,j,1,n)=help 464 endif 465 endif 466 if((isec1(6).eq.011).and.(isec1(7).eq.105).and. & 467 (nint(xsec18).eq.2)) then 468 ! 2 M TEMPERATURE 469 help=zsec4(nxfield*(ny-j-1)+i+1) 470 if(i.lt.i180) then 471 tt2(i180+i,j,1,n)=help 472 else 473 tt2(i-i180,j,1,n)=help 474 endif 475 endif 476 if((isec1(6).eq.017).and.(isec1(7).eq.105).and. & 477 (nint(xsec18).eq.2)) then 478 ! 2 M DEW POINT TEMPERATURE 479 help=zsec4(nxfield*(ny-j-1)+i+1) 480 if(i.lt.i180) then 481 td2(i180+i,j,1,n)=help 482 else 483 td2(i-i180,j,1,n)=help 484 endif 485 endif 486 if((isec1(6).eq.062).and.(isec1(7).eq.001)) then 487 ! LARGE SCALE PREC. 488 help=zsec4(nxfield*(ny-j-1)+i+1) 489 if(i.lt.i180) then 490 lsprec(i180+i,j,1,n)=help 491 else 492 lsprec(i-i180,j,1,n)=help 493 endif 494 endif 495 if((isec1(6).eq.063).and.(isec1(7).eq.001)) then 496 ! CONVECTIVE PREC. 497 help=zsec4(nxfield*(ny-j-1)+i+1) 498 if(i.lt.i180) then 499 convprec(i180+i,j,1,n)=help 500 else 501 convprec(i-i180,j,1,n)=help 502 endif 503 endif 504 if((isec1(6).eq.007).and.(isec1(7).eq.001)) then 505 ! TOPOGRAPHY 506 help=zsec4(nxfield*(ny-j-1)+i+1) 507 if(i.lt.i180) then 508 oro(i180+i,j)=help 509 excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 510 else 511 oro(i-i180,j)=help 512 excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED 513 endif 514 endif 515 if((isec1(6).eq.081).and.(isec1(7).eq.001)) then 516 ! LAND SEA MASK 517 help=zsec4(nxfield*(ny-j-1)+i+1) 518 if(i.lt.i180) then 519 lsm(i180+i,j)=help 520 else 521 lsm(i-i180,j)=help 522 endif 523 endif 524 if((isec1(6).eq.221).and.(isec1(7).eq.001)) then 525 ! MIXING HEIGHT 526 help=zsec4(nxfield*(ny-j-1)+i+1) 527 if(i.lt.i180) then 528 hmix(i180+i,j,1,n)=help 529 else 530 hmix(i-i180,j,1,n)=help 531 endif 532 endif 533 if((isec1(6).eq.052).and.(isec1(7).eq.105).and. & 534 (nint(xsec18).eq.02)) then 535 ! 2 M RELATIVE HUMIDITY 536 help=zsec4(nxfield*(ny-j-1)+i+1) 537 if(i.lt.i180) then 538 qvh2(i180+i,j)=help 539 else 540 qvh2(i-i180,j)=help 541 endif 542 endif 543 if((isec1(6).eq.011).and.(isec1(7).eq.107)) then 544 ! TEMPERATURE LOWEST SIGMA LEVEL 545 help=zsec4(nxfield*(ny-j-1)+i+1) 546 if(i.lt.i180) then 547 tlev1(i180+i,j)=help 548 else 549 tlev1(i-i180,j)=help 550 endif 551 endif 552 if((isec1(6).eq.033).and.(isec1(7).eq.107)) then 553 ! U VELOCITY LOWEST SIGMA LEVEL 554 help=zsec4(nxfield*(ny-j-1)+i+1) 555 if(i.lt.i180) then 556 ulev1(i180+i,j)=help 557 else 558 ulev1(i-i180,j)=help 559 endif 560 endif 561 if((isec1(6).eq.034).and.(isec1(7).eq.107)) then 562 ! V VELOCITY LOWEST SIGMA LEVEL 563 help=zsec4(nxfield*(ny-j-1)+i+1) 564 if(i.lt.i180) then 565 vlev1(i180+i,j)=help 566 else 567 vlev1(i-i180,j)=help 568 endif 569 endif 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 570 541 ! SEC & IP 12/2018 read GFS clouds 571 ! ESO TODO: do as for hachinger version 572 ! if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg] 573 ! if((i.eq.0).and.(j.eq.0)) then 574 ! do ii=1,nuvz 575 ! if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii 576 ! end do 577 ! endif 578 ! help=zsec4(nxfield*(ny-j-1)+i+1) 579 ! if(i.le.i180) then 580 ! clwch(i179+i,j,numpclwch,n)=help 581 ! else 582 ! clwch(i-i181,j,numpclwch,n)=help 583 ! endif 584 ! readclouds=.true. 585 ! sumclouds=.true. 586 ! endif 587 588 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 589 558 end do 590 end do591 559 592 560 endif 593 561 594 562 if((isec1(6).eq.33).and.(isec1(7).eq.100)) then 595 563 ! NCEP ISOBARIC LEVELS 596 564 iumax=iumax+1 597 565 endif … … 599 567 call grib_release(igrib) 600 568 goto 10 !! READ NEXT LEVEL OR PARAMETER 601 602 603 604 605 606 50 569 ! 570 ! CLOSING OF INPUT DATA FILE 571 ! 572 573 !HSO close grib file 574 50 continue 607 575 call grib_close_file(ifile) 608 576 609 577 ! SENS. HEAT FLUX 610 578 sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files 611 579 hflswitch=.false. ! Heat flux not available 612 580 ! SOLAR RADIATIVE FLUXES 613 581 ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files 614 582 ! EW SURFACE STRESS 615 583 ewss=0.0 ! not available from gfs.tccz.pgrbfxx files 616 584 ! NS SURFACE STRESS 617 585 nsss=0.0 ! not available from gfs.tccz.pgrbfxx files 618 586 strswitch=.false. ! stress not available 619 587 620 588 ! CONVERT TP TO LSP (GRIB2 only) 621 589 if (gribVer.eq.2) then 622 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) ) 623 591 endif 624 625 626 627 592 !HSO end edits 593 594 595 ! TRANSFORM RH TO SPECIFIC HUMIDITY 628 596 629 597 do j=0,ny-1 … … 633 601 temp=tth(i,j,k,n) 634 602 if (temp .eq. 0.0) then 635 636 603 write (*, *) i, j, k, n 604 temp = 273.0 637 605 endif 638 606 plev1=akm(k)+bkm(k)*ps(i,j,1,n) … … 643 611 end do 644 612 645 646 647 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 648 616 649 617 do j=0,ny-1 650 618 do i=0,nxfield-1 651 652 653 654 655 656 657 658 659 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) 660 628 end do 661 629 end do … … 671 639 672 640 673 674 675 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 !************************************************************************* 676 644 677 645 if (xglobal) then … … 709 677 do i=0,nxmin1 710 678 do j=0,nymin1 711 679 ! Convert precip. from mm/s -> mm/hour 712 680 convprec(i,j,1,n)=convprec(i,j,1,n)*3600. 713 681 lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600. … … 717 685 718 686 if ((.not.hflswitch).or.(.not.strswitch)) then 719 720 721 722 723 687 ! write(*,*) 'WARNING: No flux data contained in GRIB file ', 688 ! + wfname(indj) 689 690 ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD 691 !*************************************************************************** 724 692 725 693 do i=0,nxmin1 … … 741 709 742 710 return 743 888 711 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 744 712 write(*,*) ' #### ',wfname(indj),' #### ' 745 713 write(*,*) ' #### IS NOT GRIB FORMAT !!! #### ' 746 714 stop 'Execution terminated' 747 999 715 999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' 748 716 write(*,*) ' #### ',wfname(indj),' #### ' 749 717 write(*,*) ' #### CANNOT BE OPENED !!! #### '
Note: See TracChangeset
for help on using the changeset viewer.