diff --git a/src_FLEXPART91_fixes/gridcheck_gfs.f90 b/src_FLEXPART91_fixes/gridcheck_gfs.f90 new file mode 100644 index 0000000000000000000000000000000000000000..680ded16feb94a37558a098b6af678e7f355dc18 --- /dev/null +++ b/src_FLEXPART91_fixes/gridcheck_gfs.f90 @@ -0,0 +1,557 @@ +!********************************************************************** +! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * +! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * +! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * +! * +! This file is part of FLEXPART. * +! * +! FLEXPART is free software: you can redistribute it and/or modify * +! it under the terms of the GNU General Public License as published by* +! the Free Software Foundation, either version 3 of the License, or * +! (at your option) any later version. * +! * +! FLEXPART is distributed in the hope that it will be useful, * +! but WITHOUT ANY WARRANTY; without even the implied warranty of * +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +! GNU General Public License for more details. * +! * +! You should have received a copy of the GNU General Public License * +! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * +!********************************************************************** + +subroutine gridcheck + + !********************************************************************** + ! * + ! FLEXPART MODEL SUBROUTINE GRIDCHECK * + ! * + !********************************************************************** + ! * + ! AUTHOR: G. WOTAWA * + ! DATE: 1997-08-06 * + ! LAST UPDATE: 1997-10-10 * + ! * + ! Update: 1999-02-08, global fields allowed, A. Stohl* + ! CHANGE: 17/11/2005, Caroline Forster, GFS data * + ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with * + ! ECMWF grib_api * + ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * + ! ECMWF grib_api * + ! CHANGE: Jerome Brioude, sort akm * + !********************************************************************** + ! * + ! DESCRIPTION: * + ! * + ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT * + ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST- * + ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE * + ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES * + ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT * + ! ANY CALL. * + ! * + ! XLON0 geographical longitude of lower left gridpoint * + ! YLAT0 geographical latitude of lower left gridpoint * + ! NX number of grid points x-direction * + ! NY number of grid points y-direction * + ! DX grid distance x-direction * + ! DY grid distance y-direction * + ! NUVZ number of grid points for horizontal wind * + ! components in z direction * + ! NWZ number of grid points for vertical wind * + ! component in z direction * + ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid* + ! points of the polar stereographic grid): * + ! used to check the CFL criterion * + ! UVHEIGHT(1)- heights of gridpoints where u and v are * + ! UVHEIGHT(NUVZ) given * + ! WHEIGHT(1)- heights of gridpoints where w is given * + ! WHEIGHT(NWZ) * + ! * + !********************************************************************** + + use grib_api + use par_mod + use com_mod + use conv_mod + use cmapf_mod, only: stlmbr,stcm2p + + implicit none + + !HSO parameters for grib_api + integer :: ifile + integer :: iret + integer :: igrib,inc + real(kind=4) :: xaux1,xaux2,yaux1,yaux2 + real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in + integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl + !HSO end + integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip + real :: sizesouth,sizenorth,xauxa,pint + real :: akm_usort(nwzmax),akm2(nwzmax) + real,parameter :: eps=0.0001 + + ! NCEP GFS + real :: pres(nwzmax), help + + integer :: i179,i180,i181 + + ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING + + integer :: isec1(8),isec2(3) + real(kind=4) :: zsec4(jpunp) + character(len=1) :: opt + + !HSO grib api error messages + character(len=24) :: gribErrorMsg = 'Error reading grib file' + character(len=20) :: gribFunction = 'gridcheckwind_gfs' + ! + if (numbnests.ge.1) then + write(*,*) ' ###########################################' + write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:' + write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS! ' + write(*,*) ' ###########################################' + stop + endif + + iumax=0 + iwmax=0 + + if(ideltas.gt.0) then + ifn=1 + else + ifn=numbwf + endif + ! + ! OPENING OF DATA FILE (GRIB CODE) + ! +5 call grib_open_file(ifile,path(3)(1:length(3)) & + //trim(wfname(ifn)),'r',iret) + if (iret.ne.GRIB_SUCCESS) then + goto 999 ! ERROR DETECTED + endif + !turn on support for multi fields messages + call grib_multi_support_on + + ifield=0 +10 ifield=ifield+1 + ! + ! GET NEXT FIELDS + ! + call grib_new_from_file(ifile,igrib,iret) + if (iret.eq.GRIB_END_OF_FILE ) then + goto 30 ! EOF DETECTED + elseif (iret.ne.GRIB_SUCCESS) then + goto 999 ! ERROR DETECTED + endif + + !first see if we read GRIB1 or GRIB2 + call grib_get_int(igrib,'editionNumber',gribVer,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + if (gribVer.eq.1) then ! GRIB Edition 1 + + !read the grib1 identifiers + call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',isec1(8),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + !get the size and data of the values array + call grib_get_real4_array(igrib,'values',zsec4,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + else ! GRIB Edition 2 + + !read the grib2 identifiers + call grib_get_int(igrib,'discipline',discipl,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterCategory',parCat,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterNumber',parNum,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & + valSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + !convert to grib1 identifiers + isec1(6)=-1 + isec1(7)=-1 + isec1(8)=-1 + if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U + isec1(6)=33 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO + isec1(6)=7 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) & + .and.(discipl.eq.2)) then ! LSM + isec1(6)=81 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + endif + + if (isec1(6).ne.-1) then + ! get the size and data of the values array + call grib_get_real4_array(igrib,'values',zsec4,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + endif + + endif ! gribVer + + if(ifield.eq.1) then + + !get the required fields from section 2 + !store compatible to gribex input + call grib_get_int(igrib,'numberOfPointsAlongAParallel', & + isec2(2),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'numberOfPointsAlongAMeridian', & + isec2(3),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & + xaux1in,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', & + xaux2in,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', & + yaux1in,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', & + yaux2in,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + xaux1=xaux1in + xaux2=xaux2in + yaux1=yaux1in + yaux2=yaux2in + + nxfield=isec2(2) + ny=isec2(3) + if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO + xaux1=-179.0 ! 359 DEG EAST -> + xaux2=-179.0+360.-360./real(nxfield) ! TRANSFORMED TO -179 + endif ! TO 180 DEG EAST + if (xaux1.gt.180) xaux1=xaux1-360.0 + if (xaux2.gt.180) xaux2=xaux2-360.0 + if (xaux1.lt.-180) xaux1=xaux1+360.0 + if (xaux2.lt.-180) xaux2=xaux2+360.0 + if (xaux2.lt.xaux1) xaux2=xaux2+360. + xlon0=xaux1 + ylat0=yaux1 + dx=(xaux2-xaux1)/real(nxfield-1) + dy=(yaux2-yaux1)/real(ny-1) + dxconst=180./(dx*r_earth*pi) + dyconst=180./(dy*r_earth*pi) + !HSO end edits + + + ! Check whether fields are global + ! If they contain the poles, specify polar stereographic map + ! projections using the stlmbr- and stcm2p-calls + !*********************************************************** + + xauxa=abs(xaux2+dx-360.-xaux1) + if (xauxa.lt.0.001) then + nx=nxfield+1 ! field is cyclic + xglobal=.true. + if (abs(nxshift).ge.nx) & + stop 'nxshift in file par_mod is too large' + xlon0=xlon0+real(nxshift)*dx + else + nx=nxfield + xglobal=.false. + if (nxshift.ne.0) & + stop 'nxshift (par_mod) must be zero for non-global domain' + endif + nxmin1=nx-1 + nymin1=ny-1 + if (xlon0.gt.180.) xlon0=xlon0-360. + xauxa=abs(yaux1+90.) + if (xglobal.and.xauxa.lt.0.001) then + sglobal=.true. ! field contains south pole + ! Enhance the map scale by factor 3 (*2=6) compared to north-south + ! map scale + sizesouth=6.*(switchsouth+90.)/dy + call stlmbr(southpolemap,-90.,0.) + call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, & + sizesouth,switchsouth,180.) + switchsouthg=(switchsouth-ylat0)/dy + else + sglobal=.false. + switchsouthg=999999. + endif + xauxa=abs(yaux2-90.) + if (xglobal.and.xauxa.lt.0.001) then + nglobal=.true. ! field contains north pole + ! Enhance the map scale by factor 3 (*2=6) compared to north-south + ! map scale + sizenorth=6.*(90.-switchnorth)/dy + call stlmbr(northpolemap,90.,0.) + call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, & + sizenorth,switchnorth,180.) + switchnorthg=(switchnorth-ylat0)/dy + else + nglobal=.false. + switchnorthg=999999. + endif + endif ! ifield.eq.1 + + if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative' + if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large' + + ! NCEP ISOBARIC LEVELS + !********************* + + if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind + iumax=iumax+1 + pres(iumax)=real(isec1(8))*100.0 + endif + + + i179=nint(179./dx) + if (dx.lt.0.7) then + i180=nint(180./dx)+1 ! 0.5 deg data + else + i180=nint(179./dx)+1 ! 1 deg data + endif + i181=i180+1 + + + ! NCEP TERRAIN + !************* + + if((isec1(6).eq.007).and.(isec1(7).eq.001)) then + do jy=0,ny-1 + do ix=0,nxfield-1 + help=zsec4(nxfield*(ny-jy-1)+ix+1) + if(ix.le.i180) then + oro(i179+ix,jy)=help + excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED + else + oro(ix-i181,jy)=help + excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED + endif + end do + end do + endif + + ! NCEP LAND SEA MASK + !******************* + + if((isec1(6).eq.081).and.(isec1(7).eq.001)) then + do jy=0,ny-1 + do ix=0,nxfield-1 + help=zsec4(nxfield*(ny-jy-1)+ix+1) + if(ix.le.i180) then + lsm(i179+ix,jy)=help + else + lsm(ix-i181,jy)=help + endif + end do + end do + endif + + call grib_release(igrib) + + goto 10 !! READ NEXT LEVEL OR PARAMETER + ! + ! CLOSING OF INPUT DATA FILE + ! + + ! HSO +30 continue + call grib_close_file(ifile) + ! HSO end edits + + nuvz=iumax + nwz =iumax + nlev_ec=iumax + + if (nx.gt.nxmax) then + write(*,*) 'FLEXPART error: Too many grid points in x direction.' + write(*,*) 'Reduce resolution of wind fields.' + write(*,*) 'Or change parameter settings in file par_mod.' + write(*,*) nx,nxmax + stop + endif + + if (ny.gt.nymax) then + write(*,*) 'FLEXPART error: Too many grid points in y direction.' + write(*,*) 'Reduce resolution of wind fields.' + write(*,*) 'Or change parameter settings in file par_mod.' + write(*,*) ny,nymax + stop + endif + + if (nuvz.gt.nuvzmax) then + write(*,*) 'FLEXPART error: Too many u,v grid points in z '// & + 'direction.' + write(*,*) 'Reduce resolution of wind fields.' + write(*,*) 'Or change parameter settings in file par_mod.' + write(*,*) nuvz,nuvzmax + stop + endif + + if (nwz.gt.nwzmax) then + write(*,*) 'FLEXPART error: Too many w grid points in z '// & + 'direction.' + write(*,*) 'Reduce resolution of wind fields.' + write(*,*) 'Or change parameter settings in file par_mod.' + write(*,*) nwz,nwzmax + stop + endif + + ! If desired, shift all grids by nxshift grid cells + !************************************************** + + if (xglobal) then + call shift_field_0(oro,nxfield,ny) + call shift_field_0(lsm,nxfield,ny) + call shift_field_0(excessoro,nxfield,ny) + endif + + ! Output of grid info + !******************** + + write(*,*) + write(*,*) + write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', & + nuvz,nwz + write(*,*) + write(*,'(a)') 'Mother domain:' + write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', & + xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx + write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', & + ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy + write(*,*) + + + ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL + ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM + + numskip=nlev_ec-nuvz ! number of ecmwf model layers not used + ! by trajectory model + do i=1,nwz + j=numskip+i + k=nlev_ec+1+numskip+i + akm_usort(nwz-i+1)=pres(nwz-i+1) + bkm(nwz-i+1)=0.0 + end do + + !****************************** + ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted! + !****************************** + do i=1,nwz + if (akm_usort(1).gt.akm_usort(2)) then + akm(i)=akm_usort(i) + else + akm(i)=akm_usort(nwz-i+1) + endif + end do + +!*************************** +!JB sort akm which can be random + do i=1,nwzmax + akm2(i)=0. + enddo + do i=1,nwzmax + if (akm(i).gt.akm2(1)) akm2(1)=akm(i) + enddo + inc=2 + +44 do i=1,nwzmax + if(akm(i).lt.akm2(inc-1).and.akm(i).gt.akm2(inc))akm2(inc)=akm(i) + enddo + inc=inc+1 + if (inc.lt.nwzmax) goto 44 +45 continue + akm=akm2 +! end JB + + ! + ! CALCULATION OF AKZ, BKZ + ! AKZ,BKZ: model discretization parameters at the center of each model + ! layer + ! + ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0, + ! i.e. ground level + !***************************************************************************** + + do i=1,nuvz + akz(i)=akm(i) + bkz(i)=bkm(i) + end do + + ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled + ! upon the transformation to z levels. In order to save computer memory, this is + ! not done anymore in the standard version. However, this option can still be + ! switched on by replacing the following lines with those below, that are + ! currently commented out. For this, similar changes are necessary in + ! verttransform.f and verttranform_nests.f + !***************************************************************************** + + nz=nuvz + if (nz.gt.nzmax) stop 'nzmax too small' + do i=1,nuvz + aknew(i)=akz(i) + bknew(i)=bkz(i) + end do + + ! Switch on following lines to use doubled vertical resolution + !************************************************************* + !nz=nuvz+nwz-1 + !if (nz.gt.nzmax) stop 'nzmax too small' + !do 100 i=1,nwz + ! aknew(2*(i-1)+1)=akm(i) + !00 bknew(2*(i-1)+1)=bkm(i) + !do 110 i=2,nuvz + ! aknew(2*(i-1))=akz(i) + !10 bknew(2*(i-1))=bkz(i) + ! End doubled vertical resolution + + + ! Determine the uppermost level for which the convection scheme shall be applied + ! by assuming that there is no convection above 50 hPa (for standard SLP) + !***************************************************************************** + + do i=1,nuvz-2 + pint=akz(i)+bkz(i)*101325. + if (pint.lt.5000.) goto 96 + end do +96 nconvlev=i + if (nconvlev.gt.nconvlevmax-1) then + nconvlev=nconvlevmax-1 + write(*,*) 'Attention, convection only calculated up to ', & + akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa' + endif + + return + +999 write(*,*) + write(*,*) ' ###########################################'// & + '###### ' + write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:' + write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn) + write(*,*) ' ###########################################'// & + '###### ' + write(*,*) + write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!' + write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!' + write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!' + write(*,'(a)') '!!! THE "X" KEY... !!!' + write(*,*) + read(*,'(a)') opt + if(opt.eq.'X') then + stop + else + goto 5 + endif + +end subroutine gridcheck diff --git a/src_FLEXPART91_fixes/readwind_gfs.f90 b/src_FLEXPART91_fixes/readwind_gfs.f90 new file mode 100644 index 0000000000000000000000000000000000000000..64d913895e369b64683a2a22a6e228cc42956b53 --- /dev/null +++ b/src_FLEXPART91_fixes/readwind_gfs.f90 @@ -0,0 +1,719 @@ +!********************************************************************** +! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * +! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * +! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * +! * +! This file is part of FLEXPART. * +! * +! FLEXPART is free software: you can redistribute it and/or modify * +! it under the terms of the GNU General Public License as published by* +! the Free Software Foundation, either version 3 of the License, or * +! (at your option) any later version. * +! * +! FLEXPART is distributed in the hope that it will be useful, * +! but WITHOUT ANY WARRANTY; without even the implied warranty of * +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +! GNU General Public License for more details. * +! * +! You should have received a copy of the GNU General Public License * +! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * +!********************************************************************** + +subroutine readwind(indj,n,uuh,vvh,wwh) + + !*********************************************************************** + !* * + !* TRAJECTORY MODEL SUBROUTINE READWIND * + !* * + !*********************************************************************** + !* * + !* AUTHOR: G. WOTAWA * + !* DATE: 1997-08-05 * + !* LAST UPDATE: 2000-10-17, Andreas Stohl * + !* CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and * + !* qvh (on eta coordinates) in common block * + !* CHANGE: 16/11/2005, Caroline Forster, GFS data * + !* CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2 * + !* data with the ECMWF grib_api library * + !* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with * + !* ECMWF grib_api * + !* CHANGE: Jerome Brioude, update the grib1 * + !* identifier for RH2,T2,U10,V10 * + !*********************************************************************** + !* * + !* DESCRIPTION: * + !* * + !* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE * + !* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE * + !* * + !* INPUT: * + !* indj indicates number of the wind field to be read in * + !* n temporal index for meteorological fields (1 to 3)* + !* * + !* IMPORTANT VARIABLES FROM COMMON BLOCK: * + !* * + !* wfname File name of data to be read in * + !* nx,ny,nuvz,nwz expected field dimensions * + !* nlev_ec number of vertical levels ecmwf model * + !* uu,vv,ww wind fields * + !* tt,qv temperature and specific humidity * + !* ps surface pressure * + !* * + !*********************************************************************** + + use grib_api + use par_mod + use com_mod + + implicit none + + !HSO new parameters for grib_api + integer :: ifile + integer :: iret + integer :: igrib + integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl + !HSO end edits + real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) + real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax) + real :: wwh(0:nxmax-1,0:nymax-1,nwzmax) + integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax + + ! NCEP + integer :: numpt,numpu,numpv,numpw,numprh + real :: help, temp, ew + real :: elev + real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1) + real :: tlev1(0:nxmax-1,0:nymax-1) + real :: qvh2(0:nxmax-1,0:nymax-1) + + integer :: i179,i180,i181 + + ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING + !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input + + integer :: isec1(8),isec2(3) + real(kind=4) :: zsec4(jpunp) + real(kind=4) :: xaux,yaux,xaux0,yaux0 + real(kind=8) :: xauxin,yauxin + real,parameter :: eps=1.e-4 + real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1) + real :: plev1,hlev1,ff10m,fflev1 + + logical :: hflswitch,strswitch + + !HSO for grib api error messages + character(len=24) :: gribErrorMsg = 'Error reading grib file' + character(len=20) :: gribFunction = 'readwind_gfs' + + + hflswitch=.false. + strswitch=.false. + levdiff2=nlev_ec-nwz+1 + iumax=0 + iwmax=0 + + + ! OPENING OF DATA FILE (GRIB CODE) + + !HSO +5 call grib_open_file(ifile,path(3)(1:length(3)) & + //trim(wfname(indj)),'r',iret) + if (iret.ne.GRIB_SUCCESS) then + goto 888 ! ERROR DETECTED + endif + !turn on support for multi fields messages + call grib_multi_support_on + + numpt=0 + numpu=0 + numpv=0 + numpw=0 + numprh=0 + ifield=0 +10 ifield=ifield+1 + ! + ! GET NEXT FIELDS + ! + call grib_new_from_file(ifile,igrib,iret) + if (iret.eq.GRIB_END_OF_FILE) then + goto 50 ! EOF DETECTED + elseif (iret.ne.GRIB_SUCCESS) then + goto 888 ! ERROR DETECTED + endif + + !first see if we read GRIB1 or GRIB2 + call grib_get_int(igrib,'editionNumber',gribVer,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + if (gribVer.eq.1) then ! GRIB Edition 1 + + !read the grib1 identifiers + call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'level',isec1(8),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + else ! GRIB Edition 2 + + !read the grib2 identifiers + call grib_get_int(igrib,'discipline',discipl,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterCategory',parCat,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'parameterNumber',parNum,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & + valSurf,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + + !convert to grib1 identifiers + isec1(6)=-1 + isec1(7)=-1 + isec1(8)=-1 + if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T + isec1(6)=11 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U + isec1(6)=33 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V + isec1(6)=34 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W + isec1(6)=39 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH + isec1(6)=52 ! indicatorOfParameter + isec1(7)=100 ! indicatorOfTypeOfLevel + isec1(8)=valSurf/100 ! level, convert to hPa + elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103).and.(valSurf.eq.2)) then ! RH2 + isec1(6)=52 ! indicatorOfParameter + isec1(7)=105 ! indicatorOfTypeOfLevel + isec1(8)=2 + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103).and.(valSurf.eq.2)) then ! T2 + isec1(6)=11 ! indicatorOfParameter + isec1(7)=105 ! indicatorOfTypeOfLevel + isec1(8)=2 + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103).and.(valSurf.eq.10)) then ! U10 + isec1(6)=33 ! indicatorOfParameter + isec1(7)=105 ! indicatorOfTypeOfLevel + isec1(8)=10 + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103).and.(valSurf.eq.10)) then ! V10 + isec1(6)=34 ! indicatorOfParameter + isec1(7)=105 ! indicatorOfTypeOfLevel + isec1(8)=10 + elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP + isec1(6)=2 ! indicatorOfParameter + isec1(7)=102 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP + isec1(6)=1 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW + isec1(6)=66 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0 + isec1(6)=11 ! indicatorOfParameter + isec1(7)=107 ! indicatorOfTypeOfLevel + isec1(8)=0.995 ! lowest sigma level + elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0 + isec1(6)=33 ! indicatorOfParameter + isec1(7)=107 ! indicatorOfTypeOfLevel + isec1(8)=0.995 ! lowest sigma level + elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0 + isec1(6)=34 ! indicatorOfParameter + isec1(7)=107 ! indicatorOfTypeOfLevel + isec1(8)=0.995 ! lowest sigma level + elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO + isec1(6)=7 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) & + .and.(discipl.eq.2)) then ! LSM + isec1(6)=81 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH + isec1(6)=221 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP + isec1(6)=62 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP + isec1(6)=63 ! indicatorOfParameter + isec1(7)=1 ! indicatorOfTypeOfLevel + isec1(8)=0 + endif + + endif ! gribVer + + if (isec1(6).ne.-1) then + ! get the size and data of the values array + call grib_get_real4_array(igrib,'values',zsec4,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + endif + + if(ifield.eq.1) then + + !get the required fields from section 2 + !store compatible to gribex input + call grib_get_int(igrib,'numberOfPointsAlongAParallel', & + isec2(2),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_int(igrib,'numberOfPointsAlongAMeridian', & + isec2(3),iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & + xauxin,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', & + yauxin,iret) + call grib_check(iret,gribFunction,gribErrorMsg) + xaux=xauxin+real(nxshift)*dx + yaux=yauxin + + ! CHECK GRID SPECIFICATIONS + + if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT' + if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT' + if(xaux.eq.0.) xaux=-179.0 ! NCEP DATA + xaux0=xlon0 + yaux0=ylat0 + if(xaux.lt.0.) xaux=xaux+360. + if(yaux.lt.0.) yaux=yaux+360. + if(xaux0.lt.0.) xaux0=xaux0+360. + if(yaux0.lt.0.) yaux0=yaux0+360. + if(abs(xaux-xaux0).gt.eps) & + stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' + if(abs(yaux-yaux0).gt.eps) & + stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' + endif + !HSO end of edits + + i179=nint(179./dx) + if (dx.lt.0.7) then + i180=nint(180./dx)+1 ! 0.5 deg data + else + i180=nint(179./dx)+1 ! 1 deg data + endif + i181=i180+1 + + if (isec1(6).ne.-1) then + + do j=0,nymin1 + do i=0,nxfield-1 + if((isec1(6).eq.011).and.(isec1(7).eq.100)) then + ! TEMPERATURE + if((i.eq.0).and.(j.eq.0)) then + do ii=1,nuvz + if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii + end do + endif + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + tth(i179+i,j,numpt,n)=help + else + tth(i-i181,j,numpt,n)=help + endif + endif + if((isec1(6).eq.033).and.(isec1(7).eq.100)) then + ! U VELOCITY + if((i.eq.0).and.(j.eq.0)) then + do ii=1,nuvz + if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii + end do + endif + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + uuh(i179+i,j,numpu)=help + else + uuh(i-i181,j,numpu)=help + endif + endif + if((isec1(6).eq.034).and.(isec1(7).eq.100)) then + ! V VELOCITY + if((i.eq.0).and.(j.eq.0)) then + do ii=1,nuvz + if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii + end do + endif + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + vvh(i179+i,j,numpv)=help + else + vvh(i-i181,j,numpv)=help + endif + endif + if((isec1(6).eq.052).and.(isec1(7).eq.100)) then + ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER + if((i.eq.0).and.(j.eq.0)) then + do ii=1,nuvz + if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii + end do + endif + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + qvh(i179+i,j,numprh,n)=help + else + qvh(i-i181,j,numprh,n)=help + endif + endif + if((isec1(6).eq.001).and.(isec1(7).eq.001)) then + ! SURFACE PRESSURE + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + ps(i179+i,j,1,n)=help + else + ps(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.039).and.(isec1(7).eq.100)) then + ! W VELOCITY + if((i.eq.0).and.(j.eq.0)) then + do ii=1,nuvz + if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii + end do + endif + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + wwh(i179+i,j,numpw)=help + else + wwh(i-i181,j,numpw)=help + endif + endif + if((isec1(6).eq.066).and.(isec1(7).eq.001)) then + ! SNOW DEPTH + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + sd(i179+i,j,1,n)=help + else + sd(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.002).and.(isec1(7).eq.102)) then + ! MEAN SEA LEVEL PRESSURE + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + msl(i179+i,j,1,n)=help + else + msl(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.071).and.(isec1(7).eq.244)) then + ! TOTAL CLOUD COVER + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + tcc(i179+i,j,1,n)=help + else + tcc(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.033).and.(isec1(7).eq.105).and. & + (isec1(8).eq.10)) then + ! 10 M U VELOCITY + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + u10(i179+i,j,1,n)=help + else + u10(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.034).and.(isec1(7).eq.105).and. & + (isec1(8).eq.10)) then + ! 10 M V VELOCITY + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + v10(i179+i,j,1,n)=help + else + v10(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.011).and.(isec1(7).eq.105).and. & + (isec1(8).eq.02)) then + ! 2 M TEMPERATURE + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + tt2(i179+i,j,1,n)=help + else + tt2(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.017).and.(isec1(7).eq.105).and. & + (isec1(8).eq.02)) then + ! 2 M DEW POINT TEMPERATURE + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + td2(i179+i,j,1,n)=help + else + td2(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.062).and.(isec1(7).eq.001)) then + ! LARGE SCALE PREC. + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + lsprec(i179+i,j,1,n)=help + else + lsprec(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.063).and.(isec1(7).eq.001)) then + ! CONVECTIVE PREC. + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + convprec(i179+i,j,1,n)=help + else + convprec(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.007).and.(isec1(7).eq.001)) then + ! TOPOGRAPHY + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + oro(i179+i,j)=help + excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED + else + oro(i-i181,j)=help + excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED + endif + endif + if((isec1(6).eq.081).and.(isec1(7).eq.001)) then + ! LAND SEA MASK + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + lsm(i179+i,j)=help + else + lsm(i-i181,j)=help + endif + endif + if((isec1(6).eq.221).and.(isec1(7).eq.001)) then + ! MIXING HEIGHT + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + hmix(i179+i,j,1,n)=help + else + hmix(i-i181,j,1,n)=help + endif + endif + if((isec1(6).eq.052).and.(isec1(7).eq.105).and. & + (isec1(8).eq.02)) then + ! 2 M RELATIVE HUMIDITY + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + qvh2(i179+i,j)=help + else + qvh2(i-i181,j)=help + endif + endif + if((isec1(6).eq.011).and.(isec1(7).eq.107)) then + ! TEMPERATURE LOWEST SIGMA LEVEL + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + tlev1(i179+i,j)=help + else + tlev1(i-i181,j)=help + endif + endif + if((isec1(6).eq.033).and.(isec1(7).eq.107)) then + ! U VELOCITY LOWEST SIGMA LEVEL + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + ulev1(i179+i,j)=help + else + ulev1(i-i181,j)=help + endif + endif + if((isec1(6).eq.034).and.(isec1(7).eq.107)) then + ! V VELOCITY LOWEST SIGMA LEVEL + help=zsec4(nxfield*(ny-j-1)+i+1) + if(i.le.i180) then + vlev1(i179+i,j)=help + else + vlev1(i-i181,j)=help + endif + endif + + end do + end do + + endif + + if((isec1(6).eq.33).and.(isec1(7).eq.100)) then + ! NCEP ISOBARIC LEVELS + iumax=iumax+1 + endif + + call grib_release(igrib) + goto 10 !! READ NEXT LEVEL OR PARAMETER + ! + ! CLOSING OF INPUT DATA FILE + ! + + !HSO close grib file +50 continue + call grib_close_file(ifile) + + ! SENS. HEAT FLUX + sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files + hflswitch=.false. ! Heat flux not available + ! SOLAR RADIATIVE FLUXES + ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files + ! EW SURFACE STRESS + ewss=0.0 ! not available from gfs.tccz.pgrbfxx files + ! NS SURFACE STRESS + nsss=0.0 ! not available from gfs.tccz.pgrbfxx files + strswitch=.false. ! stress not available + + ! CONVERT TP TO LSP (GRIB2 only) + if (gribVer.eq.2) then + do j=0,nymin1 + do i=0,nxfield-1 + if(i.le.i180) then + if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur + lsprec(i179+i,j,1,n)= & + lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n) + else + lsprec(i179+i,j,1,n)=0 + endif + else + if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then + lsprec(i-i181,j,1,n)= & + lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n) + else + lsprec(i-i181,j,1,n)=0 + endif + endif + enddo + enddo + endif + !HSO end edits + + + ! TRANSFORM RH TO SPECIFIC HUMIDITY + + do j=0,ny-1 + do i=0,nxfield-1 + do k=1,nuvz + help=qvh(i,j,k,n) + temp=tth(i,j,k,n) + plev1=akm(k)+bkm(k)*ps(i,j,1,n) + elev=ew(temp)*help/100.0 + qvh(i,j,k,n)=xmwml*(elev/(plev1-((1.0-xmwml)*elev))) + end do + end do + end do + + ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY + ! USING BOLTON'S (1980) FORMULA + ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA + + do j=0,ny-1 + do i=0,nxfield-1 + help=qvh2(i,j) + temp=tt2(i,j,1,n) + elev=ew(temp)/100.*help/100. !vapour pressure in hPa + td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273. + if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n) + end do + end do + + if(levdiff2.eq.0) then + iwmax=nlev_ec+1 + do i=0,nxmin1 + do j=0,nymin1 + wwh(i,j,nlev_ec+1)=0. + end do + end do + endif + + + ! For global fields, assign the leftmost data column also to the rightmost + ! data column; if required, shift whole grid by nxshift grid points + !************************************************************************* + + if (xglobal) then + call shift_field_0(ewss,nxfield,ny) + call shift_field_0(nsss,nxfield,ny) + call shift_field_0(oro,nxfield,ny) + call shift_field_0(excessoro,nxfield,ny) + call shift_field_0(lsm,nxfield,ny) + call shift_field_0(ulev1,nxfield,ny) + call shift_field_0(vlev1,nxfield,ny) + call shift_field_0(tlev1,nxfield,ny) + call shift_field_0(qvh2,nxfield,ny) + call shift_field(ps,nxfield,ny,1,1,2,n) + call shift_field(sd,nxfield,ny,1,1,2,n) + call shift_field(msl,nxfield,ny,1,1,2,n) + call shift_field(tcc,nxfield,ny,1,1,2,n) + call shift_field(u10,nxfield,ny,1,1,2,n) + call shift_field(v10,nxfield,ny,1,1,2,n) + call shift_field(tt2,nxfield,ny,1,1,2,n) + call shift_field(td2,nxfield,ny,1,1,2,n) + call shift_field(lsprec,nxfield,ny,1,1,2,n) + call shift_field(convprec,nxfield,ny,1,1,2,n) + call shift_field(sshf,nxfield,ny,1,1,2,n) + call shift_field(ssr,nxfield,ny,1,1,2,n) + call shift_field(hmix,nxfield,ny,1,1,2,n) + call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n) + call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n) + call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1) + call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1) + call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1) + endif + + do i=0,nxmin1 + do j=0,nymin1 + ! Convert precip. from mm/s -> mm/hour + convprec(i,j,1,n)=convprec(i,j,1,n)*3600. + lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600. + surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2) + end do + end do + + if ((.not.hflswitch).or.(.not.strswitch)) then + ! write(*,*) 'WARNING: No flux data contained in GRIB file ', + ! + wfname(indj) + + ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD + !*************************************************************************** + + do i=0,nxmin1 + do j=0,nymin1 + hlev1=30.0 ! HEIGHT OF FIRST MODEL SIGMA LAYER + ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2) + fflev1=sqrt(ulev1(i,j)**2+vlev1(i,j)**2) + call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, & + tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, & + surfstr(i,j,1,n),sshf(i,j,1,n)) + if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200. + if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400. + end do + end do + endif + + if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT' + if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT' + + return +888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' + write(*,*) ' #### ',wfname(indj),' #### ' + write(*,*) ' #### IS NOT GRIB FORMAT !!! #### ' + stop 'Execution terminated' +999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### ' + write(*,*) ' #### ',wfname(indj),' #### ' + write(*,*) ' #### CANNOT BE OPENED !!! #### ' + stop 'Execution terminated' + +end subroutine readwind