Changes in / [4c52d0b:4138764] in flexpart.git


Ignore:
Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/calcmatrix.f90

    r9ca6e38 r92fab65  
    6767
    6868  phconv(1) = psconv
    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) )
     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) )
     80
     81  ! initialize mass fractions
     82    do kk=1,nconvlev
     83      fmassfrac(k,kk)=0.
    7784    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)
     85  end do
    8386
    84     do k = 1,nuvz-1
    85       qsconv(k) = f_qvsat( pconv(k), tconv(k) )
    86     end do
    87   end if
    88    
    89 ! initialize mass fractions
    90   fmassfrac(1:nuvz-1,1:nconvlev)=0.0
    9187
    9288  !note that Emanuel says it is important
    … …  
    9793  !******************
    9894
    99   cbmfold = cbmf
    100 ! Convert pressures to hPa, as required by Emanuel scheme
    101 !********************************************************
     95    cbmfold = cbmf
     96  ! Convert pressures to hPa, as required by Emanuel scheme
     97  !********************************************************
    10298!!$    do k=1,nconvlev     !old
    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)
     99    do k=1,nconvlev+1      !bugfix
     100      pconv_hpa(k)=pconv(k)/100.
     101      phconv_hpa(k)=phconv(k)/100.
     102    end do
     103    phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
     104    call convect(nconvlevmax, nconvlev, delt, iflag, &
     105         precip, wd, tprime, qprime, cbmf)
    110106
    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
     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
    118114
    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
     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
    126122
    127 ! Update fmassfrac
    128 ! account for mass displaced from level k to level k
     123  ! Update fmassfrac
     124  ! account for mass displaced from level k to level k
    129125
    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)
     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
    137135    end do
    138     fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
    139   end do
    140136
    141 200 continue
     137200   continue
    142138
    143139end subroutine calcmatrix

  • src/ew.f90

    r467460a r92fab65  
    1414
    1515  ew=0.
    16   if(x.le.0.) then
    17     write(*,*) 'in ew.f90: x=', x
    18     stop 'sorry: t not in [k]'
    19   end if
    20 
     16  if(x.le.0.) stop 'sorry: t not in [k]'
    2117  y=373.16/x
    2218  a=-7.90298*(y-1.)

  • src/gridcheck_gfs.f90

    • Property mode changed from 100755 to 100644
    r467460a ra756649  
    7575  real :: sizesouth,sizenorth,xauxa,pint
    7676  real :: akm_usort(nwzmax)
    77   real,parameter :: eps=spacing(2.0_4*360.0_4)
     77  real,parameter :: eps=0.0001
    7878
    7979  ! NCEP GFS
    8080  real :: pres(nwzmax), help
    8181
    82   integer :: i180
     82  integer :: i179,i180,i181
    8383
    8484  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    224224  nxfield=isec2(2)
    225225  ny=isec2(3)
    226   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
     226  if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
     227    xaux1=-179.0                             ! 359 DEG EAST ->
     228    xaux2=-179.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
    229229  endif                                      ! TO 180 DEG EAST
    230230  if (xaux1.gt.180) xaux1=xaux1-360.0
    … …  
    305305
    306306
    307   i180=nint(180./dx)    ! 0.5 deg data
     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
    308314
    309315
    … …  
    315321      do ix=0,nxfield-1
    316322        help=zsec4(nxfield*(ny-jy-1)+ix+1)
    317         if(ix.lt.i180) then
    318           oro(i180+ix,jy)=help
    319           excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     323        if(ix.le.i180) then
     324          oro(i179+ix,jy)=help
     325          excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    320326        else
    321           oro(ix-i180,jy)=help
    322           excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     327          oro(ix-i181,jy)=help
     328          excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    323329        endif
    324330      end do
    … …  
    333339      do ix=0,nxfield-1
    334340        help=zsec4(nxfield*(ny-jy-1)+ix+1)
    335         if(ix.lt.i180) then
    336           lsm(i180+ix,jy)=help
     341        if(ix.le.i180) then
     342          lsm(i179+ix,jy)=help
    337343        else
    338           lsm(ix-i180,jy)=help
     344          lsm(ix-i181,jy)=help
    339345        endif
    340346      end do
    … …  
    407413    write(*,*)
    408414    write(*,*)
    409   write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
     415    write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
    410416         nuvz,nwz
    411417    write(*,*)

  • src/par_mod.f90

    raa939a9 r759df5f  
    147147! GFS
    148148   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
    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
     149   integer :: nxshift=0 ! shift not fixed for the executable
    152150
    153151

  • src/readreleases.f90

    • Property mode changed from 100755 to 100644
    r467460a r92fab65  
    546546      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
    547547        write(*,*) 'FLEXPART MODEL ERROR'
    548         write(*,*) 'Release starts before simulation begins or ends (1)'
     548        write(*,*) 'Release starts before simulation begins or ends'
    549549        write(*,*) 'after simulation stops.'
    550         write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
     550        write(*,*) 'Make files COMMAND and RELEASES consistent.'
    551551        stop
    552552      endif
    … …  
    561561      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
    562562        write(*,*) 'FLEXPART MODEL ERROR'
    563         write(*,*) 'Release starts before simulation begins or ends (2)'
     563        write(*,*) 'Release starts before simulation begins or ends'
    564564        write(*,*) 'after simulation stops.'
    565         write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
     565        write(*,*) 'Make files COMMAND and RELEASES consistent.'
    566566        stop
    567567      endif

  • src/readwind_gfs.f90

    • Property mode changed from 100755 to 100644
    r9ca6e38 ra756649  
    44subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
    55
    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 !***********************************************************************
     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  !***********************************************************************
    4848
    4949  use eccodes
    … …  
    5353  implicit none
    5454
    55 !HSO  new parameters for grib_api
     55  !HSO  new parameters for grib_api
    5656  integer :: ifile
    5757  integer :: iret
    5858  integer :: igrib
    59   integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
    60 !HSO end edits
     59  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     60  !HSO end edits
    6161  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    6262  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    6464  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
    6565
    66 ! NCEP
     66  ! NCEP
    6767  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
    6868  real :: help, temp, ew
    … …  
    7272  real :: qvh2(0:nxmax-1,0:nymax-1)
    7373
    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
     74  integer :: i179,i180,i181
     75
     76  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
     77  !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
    7878
    7979  integer :: isec1(8),isec2(3)
    80   real           :: xsec18
    8180  real(kind=4) :: zsec4(jpunp)
    8281  real(kind=4) :: xaux,yaux,xaux0,yaux0
    83   real,parameter :: eps=spacing(2.0_4*360.0_4)
    8482  real(kind=8) :: xauxin,yauxin
     83  real,parameter :: eps=1.e-4
    8584  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
    8685  real :: plev1,hlev1,ff10m,fflev1
    … …  
    8887  logical :: hflswitch,strswitch
    8988
    90 !HSO  for grib api error messages
     89  !HSO  for grib api error messages
    9190  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    9291  character(len=20) :: gribFunction = 'readwind_gfs'
    … …  
    101100
    102101
    103 ! OPENING OF DATA FILE (GRIB CODE)
    104 
    105 !HSO
     102  ! OPENING OF DATA FILE (GRIB CODE)
     103
     104  !HSO
    106105  call grib_open_file(ifile,path(3)(1:length(3)) &
    107        //trim(wfname(indj)),'r',iret)
     106         //trim(wfname(indj)),'r',iret)
    108107  if (iret.ne.GRIB_SUCCESS) then
    109108    goto 888   ! ERROR DETECTED
    110109  endif
    111 !turn on support for multi fields messages
     110  !turn on support for multi fields messages
    112111  call grib_multi_support_on
    113112
    … …  
    119118  numpclwch=0
    120119  ifield=0
    121 10 ifield=ifield+1
    122 !
    123 ! GET NEXT FIELDS
    124 !
     12010   ifield=ifield+1
     121  !
     122  ! GET NEXT FIELDS
     123  !
    125124  call grib_new_from_file(ifile,igrib,iret)
    126125  if (iret.eq.GRIB_END_OF_FILE)  then
    … …  
    130129  endif
    131130
    132 !first see if we read GRIB1 or GRIB2
     131  !first see if we read GRIB1 or GRIB2
    133132  call grib_get_int(igrib,'editionNumber',gribVer,iret)
    134133!  call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    136135  if (gribVer.eq.1) then ! GRIB Edition 1
    137136
    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))
     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)
    149144
    150145  else ! GRIB Edition 2
    151146
    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
     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
    261253
    262254  endif ! gribVer
    263255
    264256  if (isec1(6).ne.-1) then
    265 !  get the size and data of the values array
     257  !  get the size and data of the values array
    266258    call grib_get_real4_array(igrib,'values',zsec4,iret)
    267     call grib_check(iret,gribFunction,gribErrorMsg)
     259!    call grib_check(iret,gribFunction,gribErrorMsg)
    268260  endif
    269261
    270262  if(ifield.eq.1) then
    271263
    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
     264  !get the required fields from section 2
     265  !store compatible to gribex input
     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)*dx
     279  yaux=yauxin
     280
     281  ! CHECK GRID SPECIFICATIONS
    290282
    291283    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    292284    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    293     if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
     285    if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
    294286    xaux0=xlon0
    295287    yaux0=ylat0
    … …  
    298290    if(xaux0.lt.0.) xaux0=xaux0+360.
    299291    if(yaux0.lt.0.) yaux0=yaux0+360.
    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)
     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
    312306
    313307  if (isec1(6).ne.-1) then
    314308
    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
     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
    541541! 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             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
     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
    558561    end do
     562  end do
    559563
    560564  endif
    561565
    562566  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
    563 ! NCEP ISOBARIC LEVELS
     567  ! NCEP ISOBARIC LEVELS
    564568    iumax=iumax+1
    565569  endif
    … …  
    567571  call grib_release(igrib)
    568572  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    569 !
    570 ! CLOSING OF INPUT DATA FILE
    571 !
    572 
    573 !HSO close grib file
    574 50 continue
     573  !
     574  ! CLOSING OF INPUT DATA FILE
     575  !
     576
     577  !HSO close grib file
     57850   continue
    575579  call grib_close_file(ifile)
    576580
    577 ! SENS. HEAT FLUX
     581  ! SENS. HEAT FLUX
    578582  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
    579583  hflswitch=.false.    ! Heat flux not available
    580 ! SOLAR RADIATIVE FLUXES
     584  ! SOLAR RADIATIVE FLUXES
    581585  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
    582 ! EW SURFACE STRESS
     586  ! EW SURFACE STRESS
    583587  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
    584 ! NS SURFACE STRESS
     588  ! NS SURFACE STRESS
    585589  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
    586590  strswitch=.false.    ! stress not available
    587591
    588 ! CONVERT TP TO LSP (GRIB2 only)
     592  ! CONVERT TP TO LSP (GRIB2 only)
    589593  if (gribVer.eq.2) then
    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
     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
    596618
    597619  do j=0,ny-1
    … …  
    600622        help=qvh(i,j,k,n)
    601623        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
    606624        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
    607625        elev=ew(temp)*help/100.0
    … …  
    611629  end do
    612630
    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
     631  ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
     632  ! USING BOLTON'S (1980) FORMULA
     633  ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
    616634
    617635  do j=0,ny-1
    618636    do i=0,nxfield-1
    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)
     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)
    628642    end do
    629643  end do
    … …  
    639653
    640654
    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 !*************************************************************************
     655  ! For global fields, assign the leftmost data column also to the rightmost
     656  ! data column; if required, shift whole grid by nxshift grid points
     657  !*************************************************************************
    644658
    645659  if (xglobal) then
    … …  
    677691  do i=0,nxmin1
    678692    do j=0,nymin1
    679 ! Convert precip. from mm/s -> mm/hour
     693  ! Convert precip. from mm/s -> mm/hour
    680694      convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
    681695      lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
    … …  
    685699
    686700  if ((.not.hflswitch).or.(.not.strswitch)) then
    687 !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
    688 !    +  wfname(indj)
    689 
    690 ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
    691 !***************************************************************************
     701  !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
     702  !    +  wfname(indj)
     703
     704  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
     705  !***************************************************************************
    692706
    693707    do i=0,nxmin1
    … …  
    709723
    710724  return
    711 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     725888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    712726  write(*,*) ' #### ',wfname(indj),'                    #### '
    713727  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
    714728  stop 'Execution terminated'
    715 999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     729999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    716730  write(*,*) ' #### ',wfname(indj),'                    #### '
    717731  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '

  • src/richardson.f90

    • Property mode changed from 100755 to 100644
    r467460a r92fab65  
    141141    zold=z
    142142  end do
    143   ! k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
     143  k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
    14414420   continue
    145145
    146146  ! Determine Richardson number between the critical levels
    147147  !********************************************************
    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
    152148
    153149  zl1=zold
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG