Changes in src/readwind_gfs.f90 [a756649:9ca6e38] in flexpart.git


Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/readwind_gfs.f90

    • Property mode changed from 100644 to 100755
    ra756649 r9ca6e38  
    1 
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
    3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    5 !                                                                     *
    6 ! This file is part of FLEXPART.                                      *
    7 !                                                                     *
    8 ! FLEXPART is free software: you can redistribute it and/or modify    *
    9 ! it under the terms of the GNU General Public License as published by*
    10 ! the Free Software Foundation, either version 3 of the License, or   *
    11 ! (at your option) any later version.                                 *
    12 !                                                                     *
    13 ! FLEXPART is distributed in the hope that it will be useful,         *
    14 ! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
    15 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
    16 ! GNU General Public License for more details.                        *
    17 !                                                                     *
    18 ! You should have received a copy of the GNU General Public License   *
    19 ! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
    20 !**********************************************************************
     1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
     2! SPDX-License-Identifier: GPL-3.0-or-later
    213
    224subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
    235
    24   !***********************************************************************
    25   !*                                                                     *
    26   !*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
    27   !*                                                                     *
    28   !***********************************************************************
    29   !*                                                                     *
    30   !*             AUTHOR:      G. WOTAWA                                  *
    31   !*             DATE:        1997-08-05                                 *
    32   !*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
    33   !*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
    34   !*                     qvh (on eta coordinates) in common block        *
    35   !*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
    36   !*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
    37   !*                     data with the ECMWF grib_api library            *
    38   !*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    39   !*                                 ECMWF grib_api                      *
    40   !                                                                      *
    41   !   Unified ECMWF and GFS builds                                       *
    42   !   Marian Harustak, 12.5.2017                                         *
    43   !     - Renamed routine from readwind to readwind_gfs                  *
    44   !*                                                                     *
    45   !***********************************************************************
    46   !*                                                                     *
    47   !* DESCRIPTION:                                                        *
    48   !*                                                                     *
    49   !* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
    50   !* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
    51   !*                                                                     *
    52   !* INPUT:                                                              *
    53   !* indj               indicates number of the wind field to be read in *
    54   !* n                  temporal index for meteorological fields (1 to 3)*
    55   !*                                                                     *
    56   !* IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
    57   !*                                                                     *
    58   !* wfname             File name of data to be read in                  *
    59   !* nx,ny,nuvz,nwz     expected field dimensions                        *
    60   !* nlev_ec            number of vertical levels ecmwf model            *
    61   !* uu,vv,ww           wind fields                                      *
    62   !* tt,qv              temperature and specific humidity                *
    63   !* ps                 surface pressure                                 *
    64   !*                                                                     *
    65   !***********************************************************************
     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!***********************************************************************
    6648
    6749  use eccodes
    … …  
    7153  implicit none
    7254
    73   !HSO  new parameters for grib_api
     55!HSO  new parameters for grib_api
    7456  integer :: ifile
    7557  integer :: iret
    7658  integer :: igrib
    77   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
    78   !HSO end edits
     59  integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
     60!HSO end edits
    7961  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    8062  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8264  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
    8365
    84   ! NCEP
     66! NCEP
    8567  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
    8668  real :: help, temp, ew
    … …  
    9072  real :: qvh2(0:nxmax-1,0:nymax-1)
    9173
    92   integer :: i179,i180,i181
    93 
    94   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    95   !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
     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
    9678
    9779  integer :: isec1(8),isec2(3)
     80  real           :: xsec18
    9881  real(kind=4) :: zsec4(jpunp)
    9982  real(kind=4) :: xaux,yaux,xaux0,yaux0
     83  real,parameter :: eps=spacing(2.0_4*360.0_4)
    10084  real(kind=8) :: xauxin,yauxin
    101   real,parameter :: eps=1.e-4
    10285  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
    10386  real :: plev1,hlev1,ff10m,fflev1
    … …  
    10588  logical :: hflswitch,strswitch
    10689
    107   !HSO  for grib api error messages
     90!HSO  for grib api error messages
    10891  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    10992  character(len=20) :: gribFunction = 'readwind_gfs'
    … …  
    118101
    119102
    120   ! OPENING OF DATA FILE (GRIB CODE)
    121 
    122   !HSO
     103! OPENING OF DATA FILE (GRIB CODE)
     104
     105!HSO
    123106  call grib_open_file(ifile,path(3)(1:length(3)) &
    124          //trim(wfname(indj)),'r',iret)
     107       //trim(wfname(indj)),'r',iret)
    125108  if (iret.ne.GRIB_SUCCESS) then
    126109    goto 888   ! ERROR DETECTED
    127110  endif
    128   !turn on support for multi fields messages
     111!turn on support for multi fields messages
    129112  call grib_multi_support_on
    130113
    … …  
    136119  numpclwch=0
    137120  ifield=0
    138 10   ifield=ifield+1
    139   !
    140   ! GET NEXT FIELDS
    141   !
     12110 ifield=ifield+1
     122!
     123! GET NEXT FIELDS
     124!
    142125  call grib_new_from_file(ifile,igrib,iret)
    143126  if (iret.eq.GRIB_END_OF_FILE)  then
    … …  
    147130  endif
    148131
    149   !first see if we read GRIB1 or GRIB2
     132!first see if we read GRIB1 or GRIB2
    150133  call grib_get_int(igrib,'editionNumber',gribVer,iret)
    151134!  call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    153136  if (gribVer.eq.1) then ! GRIB Edition 1
    154137
    155   !read the grib1 identifiers
    156   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    157 !  call grib_check(iret,gribFunction,gribErrorMsg)
    158   call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
    159 !  call grib_check(iret,gribFunction,gribErrorMsg)
    160   call grib_get_int(igrib,'level',isec1(8),iret)
    161 !  call grib_check(iret,gribFunction,gribErrorMsg)
     138!read the grib1 identifiers
     139
     140    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     141    call grib_check(iret,gribFunction,gribErrorMsg)
     142    call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
     143    call grib_check(iret,gribFunction,gribErrorMsg)
     144    call grib_get_int(igrib,'level',isec1(8),iret)
     145    call grib_check(iret,gribFunction,gribErrorMsg)
     146!JMA / SH: isec1(8) not evaluated any more below
     147!b/c with GRIB 2 this may be a real variable
     148    xsec18 = real(isec1(8))
    162149
    163150  else ! GRIB Edition 2
    164151
    165   !read the grib2 identifiers
    166   call grib_get_string(igrib,'shortName',shortname,iret)
    167 
    168   call grib_get_int(igrib,'discipline',discipl,iret)
    169 !  call grib_check(iret,gribFunction,gribErrorMsg)
    170   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    171 !  call grib_check(iret,gribFunction,gribErrorMsg)
    172   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    173 !  call grib_check(iret,gribFunction,gribErrorMsg)
    174   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    175 !  call grib_check(iret,gribFunction,gribErrorMsg)
    176   call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
    177        valSurf,iret)
    178 !  call grib_check(iret,gribFunction,gribErrorMsg)
    179  
    180 !  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
    181   !convert to grib1 identifiers
    182   isec1(6)=-1
    183   isec1(7)=-1
    184   isec1(8)=-1
    185   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
    186     isec1(6)=11          ! indicatorOfParameter
    187     isec1(7)=100         ! indicatorOfTypeOfLevel
    188     isec1(8)=valSurf/100 ! level, convert to hPa
    189   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
    190     isec1(6)=33          ! indicatorOfParameter
    191     isec1(7)=100         ! indicatorOfTypeOfLevel
    192     isec1(8)=valSurf/100 ! level, convert to hPa
    193   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
    194     isec1(6)=34          ! indicatorOfParameter
    195     isec1(7)=100         ! indicatorOfTypeOfLevel
    196     isec1(8)=valSurf/100 ! level, convert to hPa
    197   elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
    198     isec1(6)=39          ! indicatorOfParameter
    199     isec1(7)=100         ! indicatorOfTypeOfLevel
    200     isec1(8)=valSurf/100 ! level, convert to hPa
    201   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
    202     isec1(6)=52          ! indicatorOfParameter
    203     isec1(7)=100         ! indicatorOfTypeOfLevel
    204     isec1(8)=valSurf/100 ! level, convert to hPa
    205   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
    206     isec1(6)=52          ! indicatorOfParameter
    207     isec1(7)=105         ! indicatorOfTypeOfLevel
    208     isec1(8)=2
    209   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
    210     isec1(6)=11          ! indicatorOfParameter
    211     isec1(7)=105         ! indicatorOfTypeOfLevel
    212     isec1(8)=2
    213   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
    214     isec1(6)=33          ! indicatorOfParameter
    215     isec1(7)=105         ! indicatorOfTypeOfLevel
    216     isec1(8)=10
    217   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
    218     isec1(6)=34          ! indicatorOfParameter
    219     isec1(7)=105         ! indicatorOfTypeOfLevel
    220     isec1(8)=10
    221   elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
    222     isec1(6)=153         ! indicatorOfParameter
    223     isec1(7)=100         ! indicatorOfTypeOfLevel
    224     isec1(8)=valSurf/100 ! level, convert to hPa
    225   elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
    226     isec1(6)=2           ! indicatorOfParameter
    227     isec1(7)=102         ! indicatorOfTypeOfLevel
    228     isec1(8)=0
    229   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
    230     isec1(6)=1           ! indicatorOfParameter
    231     isec1(7)=1           ! indicatorOfTypeOfLevel
    232     isec1(8)=0
    233   elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
    234     isec1(6)=66          ! indicatorOfParameter
    235     isec1(7)=1           ! indicatorOfTypeOfLevel
    236     isec1(8)=0
    237   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
    238     isec1(6)=11          ! indicatorOfParameter
    239     isec1(7)=107         ! indicatorOfTypeOfLevel
    240     isec1(8)=0.995       ! lowest sigma level
    241   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
    242     isec1(6)=33          ! indicatorOfParameter
    243     isec1(7)=107         ! indicatorOfTypeOfLevel
    244     isec1(8)=0.995       ! lowest sigma level
    245   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
    246     isec1(6)=34          ! indicatorOfParameter
    247     isec1(7)=107         ! indicatorOfTypeOfLevel
    248     isec1(8)=0.995       ! lowest sigma level
    249   elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
    250     isec1(6)=7           ! indicatorOfParameter
    251     isec1(7)=1           ! indicatorOfTypeOfLevel
    252     isec1(8)=0
    253   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
    254        .and.(discipl.eq.2)) then ! LSM
    255     isec1(6)=81          ! indicatorOfParameter
    256     isec1(7)=1           ! indicatorOfTypeOfLevel
    257     isec1(8)=0
    258   elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
    259     isec1(6)=221         ! indicatorOfParameter
    260     isec1(7)=1           ! indicatorOfTypeOfLevel
    261     isec1(8)=0
    262   elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
    263     isec1(6)=62          ! indicatorOfParameter
    264     isec1(7)=1           ! indicatorOfTypeOfLevel
    265     isec1(8)=0
    266   elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
    267     isec1(6)=63          ! indicatorOfParameter
    268     isec1(7)=1           ! indicatorOfTypeOfLevel
    269     isec1(8)=0
    270   endif
     152!read the grib2 identifiers
     153   
     154    call grib_get_string(igrib,'shortName',shortname,iret)
     155
     156    call grib_get_int(igrib,'discipline',discipl,iret)
     157    call grib_check(iret,gribFunction,gribErrorMsg)
     158    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     159    call grib_check(iret,gribFunction,gribErrorMsg)
     160    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     161    call grib_check(iret,gribFunction,gribErrorMsg)
     162    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     163    call grib_check(iret,gribFunction,gribErrorMsg)
     164!
     165    call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
     166         valSurf,iret)
     167    call grib_check(iret,gribFunction,gribErrorMsg)
     168
     169!    write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
     170!convert to grib1 identifiers
     171    isec1(6:8)=-1
     172    xsec18  =-1.0
     173!JMA / SH: isec1(8) not evaluated any more below
     174!b/c with GRIB 2 this may be a real variable
     175    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
     176      isec1(6)=11          ! indicatorOfParameter
     177      isec1(7)=100         ! indicatorOfTypeOfLevel
     178      xsec18=valSurf/100.0 ! level, convert to hPa
     179    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
     180      isec1(6)=33          ! indicatorOfParameter
     181      isec1(7)=100         ! indicatorOfTypeOfLevel
     182      xsec18=valSurf/100.0 ! level, convert to hPa
     183    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
     184      isec1(6)=34          ! indicatorOfParameter
     185      isec1(7)=100         ! indicatorOfTypeOfLevel
     186      xsec18=valSurf/100.0 ! level, convert to hPa
     187    elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
     188      isec1(6)=39          ! indicatorOfParameter
     189      isec1(7)=100         ! indicatorOfTypeOfLevel
     190      xsec18=valSurf/100.0 ! level, convert to hPa
     191    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
     192      isec1(6)=52          ! indicatorOfParameter
     193      isec1(7)=100         ! indicatorOfTypeOfLevel
     194      xsec18=valSurf/100.0 ! level, convert to hPa
     195    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
     196      isec1(6)=52          ! indicatorOfParameter
     197      isec1(7)=105         ! indicatorOfTypeOfLevel
     198      xsec18=real(2)
     199    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
     200      isec1(6)=11          ! indicatorOfParameter
     201      isec1(7)=105         ! indicatorOfTypeOfLevel
     202      xsec18=real(2)
     203    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
     204      isec1(6)=33          ! indicatorOfParameter
     205      isec1(7)=105         ! indicatorOfTypeOfLevel
     206      xsec18=real(10)
     207    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
     208      isec1(6)=34          ! indicatorOfParameter
     209      isec1(7)=105         ! indicatorOfTypeOfLevel
     210      xsec18=real(10)
     211    elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
     212      isec1(6)=153         ! indicatorOfParameter
     213      isec1(7)=100         ! indicatorOfTypeOfLevel
     214      xsec18=valSurf/100.0 ! level, convert to hPa
     215    elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
     216      isec1(6)=2           ! indicatorOfParameter
     217      isec1(7)=102         ! indicatorOfTypeOfLevel
     218      xsec18=real(0)
     219    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
     220      isec1(6)=1           ! indicatorOfParameter
     221      isec1(7)=1           ! indicatorOfTypeOfLevel
     222      xsec18=real(0)
     223    elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
     224      isec1(6)=66          ! indicatorOfParameter
     225      isec1(7)=1           ! indicatorOfTypeOfLevel
     226      xsec18=real(0)
     227    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
     228      isec1(6)=11          ! indicatorOfParameter
     229      isec1(7)=107         ! indicatorOfTypeOfLevel
     230      xsec18=0.995         ! lowest sigma level
     231    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
     232      isec1(6)=33          ! indicatorOfParameter
     233      isec1(7)=107         ! indicatorOfTypeOfLevel
     234      xsec18=0.995         ! lowest sigma level
     235    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
     236      isec1(6)=34          ! indicatorOfParameter
     237      isec1(7)=107         ! indicatorOfTypeOfLevel
     238      xsec18=0.995         ! lowest sigma level
     239    elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
     240      isec1(6)=7           ! indicatorOfParameter
     241      isec1(7)=1           ! indicatorOfTypeOfLevel
     242      xsec18=real(0)
     243    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
     244         .and.(discipl.eq.2)) then ! LSM
     245      isec1(6)=81          ! indicatorOfParameter
     246      isec1(7)=1           ! indicatorOfTypeOfLevel
     247      xsec18=real(0)
     248    elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
     249      isec1(6)=221         ! indicatorOfParameter
     250      isec1(7)=1           ! indicatorOfTypeOfLevel
     251      xsec18=real(0)
     252    elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
     253      isec1(6)=62          ! indicatorOfParameter
     254      isec1(7)=1           ! indicatorOfTypeOfLevel
     255      xsec18=real(0)
     256    elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
     257      isec1(6)=63          ! indicatorOfParameter
     258      isec1(7)=1           ! indicatorOfTypeOfLevel
     259      xsec18=real(0)
     260    endif
    271261
    272262  endif ! gribVer
    273263
    274264  if (isec1(6).ne.-1) then
    275   !  get the size and data of the values array
     265!  get the size and data of the values array
    276266    call grib_get_real4_array(igrib,'values',zsec4,iret)
    277 !    call grib_check(iret,gribFunction,gribErrorMsg)
     267    call grib_check(iret,gribFunction,gribErrorMsg)
    278268  endif
    279269
    280270  if(ifield.eq.1) then
    281271
    282   !get the required fields from section 2
    283   !store compatible to gribex input
    284   call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    285        isec2(2),iret)
    286 !  call grib_check(iret,gribFunction,gribErrorMsg)
    287   call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
    288        isec2(3),iret)
    289 !  call grib_check(iret,gribFunction,gribErrorMsg)
    290   call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    291        xauxin,iret)
    292 !  call grib_check(iret,gribFunction,gribErrorMsg)
    293   call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
    294        yauxin,iret)
    295 !  call grib_check(iret,gribFunction,gribErrorMsg)
    296   xaux=xauxin+real(nxshift)*dx
    297   yaux=yauxin
    298 
    299   ! CHECK GRID SPECIFICATIONS
     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
    300290
    301291    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    302292    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    303     if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
     293    if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
    304294    xaux0=xlon0
    305295    yaux0=ylat0
    … …  
    308298    if(xaux0.lt.0.) xaux0=xaux0+360.
    309299    if(yaux0.lt.0.) yaux0=yaux0+360.
    310     if(abs(xaux-xaux0).gt.eps) &
    311          stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    312     if(abs(yaux-yaux0).gt.eps) &
    313          stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
    314   endif
    315   !HSO end of edits
    316 
    317   i179=nint(179./dx)
    318   if (dx.lt.0.7) then
    319     i180=nint(180./dx)+1    ! 0.5 deg data
    320   else
    321     i180=nint(179./dx)+1    ! 1 deg data
    322   endif
    323   i181=i180+1
     300    if(abs(xaux-xaux0).gt.eps) then
     301      write (*, *) xaux, xaux0
     302      stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
     303    endif
     304    if(abs(yaux-yaux0).gt.eps) then
     305      write (*, *) yaux, yaux0
     306      stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     307    end if
     308  endif
     309!HSO end of edits
     310
     311  i180=nint(180./dx)
    324312
    325313  if (isec1(6).ne.-1) then
    326314
    327   do j=0,nymin1
    328     do i=0,nxfield-1
    329       if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
    330   ! TEMPERATURE
    331          if((i.eq.0).and.(j.eq.0)) then
    332             do ii=1,nuvz
    333               if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
    334             end do
    335         endif
    336         help=zsec4(nxfield*(ny-j-1)+i+1)
    337         if(i.le.i180) then
    338           tth(i179+i,j,numpt,n)=help
    339         else
    340           tth(i-i181,j,numpt,n)=help
    341         endif
    342       endif
    343       if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
    344   ! U VELOCITY
    345          if((i.eq.0).and.(j.eq.0)) then
    346             do ii=1,nuvz
    347               if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
    348             end do
    349         endif
    350         help=zsec4(nxfield*(ny-j-1)+i+1)
    351         if(i.le.i180) then
    352           uuh(i179+i,j,numpu)=help
    353         else
    354           uuh(i-i181,j,numpu)=help
    355         endif
    356       endif
    357       if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
    358   ! V VELOCITY
    359          if((i.eq.0).and.(j.eq.0)) then
    360             do ii=1,nuvz
    361               if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
    362             end do
    363         endif
    364         help=zsec4(nxfield*(ny-j-1)+i+1)
    365         if(i.le.i180) then
    366           vvh(i179+i,j,numpv)=help
    367         else
    368           vvh(i-i181,j,numpv)=help
    369         endif
    370       endif
    371       if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
    372   ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
    373          if((i.eq.0).and.(j.eq.0)) then
    374             do ii=1,nuvz
    375               if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
    376             end do
    377         endif
    378         help=zsec4(nxfield*(ny-j-1)+i+1)
    379         if(i.le.i180) then
    380           qvh(i179+i,j,numprh,n)=help
    381         else
    382           qvh(i-i181,j,numprh,n)=help
    383         endif
    384       endif
    385       if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
    386   ! SURFACE PRESSURE
    387         help=zsec4(nxfield*(ny-j-1)+i+1)
    388         if(i.le.i180) then
    389           ps(i179+i,j,1,n)=help
    390         else
    391           ps(i-i181,j,1,n)=help
    392         endif
    393       endif
    394       if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
    395   ! W VELOCITY
    396          if((i.eq.0).and.(j.eq.0)) then
    397             do ii=1,nuvz
    398               if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
    399             end do
    400         endif
    401         help=zsec4(nxfield*(ny-j-1)+i+1)
    402         if(i.le.i180) then
    403           wwh(i179+i,j,numpw)=help
    404         else
    405           wwh(i-i181,j,numpw)=help
    406         endif
    407       endif
    408       if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
    409   ! SNOW DEPTH
    410         help=zsec4(nxfield*(ny-j-1)+i+1)
    411         if(i.le.i180) then
    412           sd(i179+i,j,1,n)=help
    413         else
    414           sd(i-i181,j,1,n)=help
    415         endif
    416       endif
    417       if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
    418   ! MEAN SEA LEVEL PRESSURE
    419         help=zsec4(nxfield*(ny-j-1)+i+1)
    420         if(i.le.i180) then
    421           msl(i179+i,j,1,n)=help
    422         else
    423           msl(i-i181,j,1,n)=help
    424         endif
    425       endif
    426       if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
    427   ! TOTAL CLOUD COVER
    428         help=zsec4(nxfield*(ny-j-1)+i+1)
    429         if(i.le.i180) then
    430           tcc(i179+i,j,1,n)=help
    431         else
    432           tcc(i-i181,j,1,n)=help
    433         endif
    434       endif
    435       if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
    436            (isec1(8).eq.10)) then
    437   ! 10 M U VELOCITY
    438         help=zsec4(nxfield*(ny-j-1)+i+1)
    439         if(i.le.i180) then
    440           u10(i179+i,j,1,n)=help
    441         else
    442           u10(i-i181,j,1,n)=help
    443         endif
    444       endif
    445       if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
    446            (isec1(8).eq.10)) then
    447   ! 10 M V VELOCITY
    448         help=zsec4(nxfield*(ny-j-1)+i+1)
    449         if(i.le.i180) then
    450           v10(i179+i,j,1,n)=help
    451         else
    452           v10(i-i181,j,1,n)=help
    453         endif
    454       endif
    455       if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
    456            (isec1(8).eq.02)) then
    457   ! 2 M TEMPERATURE
    458         help=zsec4(nxfield*(ny-j-1)+i+1)
    459         if(i.le.i180) then
    460           tt2(i179+i,j,1,n)=help
    461         else
    462           tt2(i-i181,j,1,n)=help
    463         endif
    464       endif
    465       if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
    466            (isec1(8).eq.02)) then
    467   ! 2 M DEW POINT TEMPERATURE
    468         help=zsec4(nxfield*(ny-j-1)+i+1)
    469         if(i.le.i180) then
    470           td2(i179+i,j,1,n)=help
    471         else
    472           td2(i-i181,j,1,n)=help
    473         endif
    474       endif
    475       if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
    476   ! LARGE SCALE PREC.
    477         help=zsec4(nxfield*(ny-j-1)+i+1)
    478         if(i.le.i180) then
    479           lsprec(i179+i,j,1,n)=help
    480         else
    481           lsprec(i-i181,j,1,n)=help
    482         endif
    483       endif
    484       if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
    485   ! CONVECTIVE PREC.
    486         help=zsec4(nxfield*(ny-j-1)+i+1)
    487         if(i.le.i180) then
    488           convprec(i179+i,j,1,n)=help
    489         else
    490           convprec(i-i181,j,1,n)=help
    491         endif
    492       endif
    493       if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
    494   ! TOPOGRAPHY
    495         help=zsec4(nxfield*(ny-j-1)+i+1)
    496         if(i.le.i180) then
    497           oro(i179+i,j)=help
    498           excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    499         else
    500           oro(i-i181,j)=help
    501           excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    502         endif
    503       endif
    504       if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
    505   ! LAND SEA MASK
    506         help=zsec4(nxfield*(ny-j-1)+i+1)
    507         if(i.le.i180) then
    508           lsm(i179+i,j)=help
    509         else
    510           lsm(i-i181,j)=help
    511         endif
    512       endif
    513       if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
    514   ! MIXING HEIGHT
    515         help=zsec4(nxfield*(ny-j-1)+i+1)
    516         if(i.le.i180) then
    517           hmix(i179+i,j,1,n)=help
    518         else
    519           hmix(i-i181,j,1,n)=help
    520         endif
    521       endif
    522       if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
    523            (isec1(8).eq.02)) then
    524   ! 2 M RELATIVE HUMIDITY
    525         help=zsec4(nxfield*(ny-j-1)+i+1)
    526         if(i.le.i180) then
    527           qvh2(i179+i,j)=help
    528         else
    529           qvh2(i-i181,j)=help
    530         endif
    531       endif
    532       if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
    533   ! TEMPERATURE LOWEST SIGMA LEVEL
    534         help=zsec4(nxfield*(ny-j-1)+i+1)
    535         if(i.le.i180) then
    536           tlev1(i179+i,j)=help
    537         else
    538           tlev1(i-i181,j)=help
    539         endif
    540       endif
    541       if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
    542   ! U VELOCITY LOWEST SIGMA LEVEL
    543         help=zsec4(nxfield*(ny-j-1)+i+1)
    544         if(i.le.i180) then
    545           ulev1(i179+i,j)=help
    546         else
    547           ulev1(i-i181,j)=help
    548         endif
    549       endif
    550       if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
    551   ! V VELOCITY LOWEST SIGMA LEVEL
    552         help=zsec4(nxfield*(ny-j-1)+i+1)
    553         if(i.le.i180) then
    554           vlev1(i179+i,j)=help
    555         else
    556           vlev1(i-i181,j)=help
    557         endif
    558       endif
     315! write (*, *) 'nxfield: ', nxfield, i180
     316
     317    do j=0,nymin1
     318      do i=0,nxfield-1
     319        if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
     320! TEMPERATURE
     321          if((i.eq.0).and.(j.eq.0)) then
     322            numpt=minloc(abs(xsec18*100.0-akz),dim=1)
     323          endif
     324          help=zsec4(nxfield*(ny-j-1)+i+1)
     325          if (help.eq.0) then
     326            write (*, *) 'i, j: ', i, j
     327            stop 'help == 0.0'
     328          endif
     329          if(i.lt.i180) then
     330            tth(i180+i,j,numpt,n)=help
     331          else
     332            tth(i-i180,j,numpt,n)=help
     333          endif
     334        endif
     335        if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
     336! U VELOCITY
     337          if((i.eq.0).and.(j.eq.0)) then
     338            numpu=minloc(abs(xsec18*100.0-akz),dim=1)
     339          endif
     340          help=zsec4(nxfield*(ny-j-1)+i+1)
     341          if(i.lt.i180) then
     342            uuh(i180+i,j,numpu)=help
     343          else
     344            uuh(i-i180,j,numpu)=help
     345          endif
     346        endif
     347        if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
     348! V VELOCITY
     349          if((i.eq.0).and.(j.eq.0)) then
     350            numpv=minloc(abs(xsec18*100.0-akz),dim=1)
     351          endif
     352          help=zsec4(nxfield*(ny-j-1)+i+1)
     353          if(i.lt.i180) then
     354            vvh(i180+i,j,numpv)=help
     355          else
     356            vvh(i-i180,j,numpv)=help
     357          endif
     358        endif
     359        if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
     360! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
     361          if((i.eq.0).and.(j.eq.0)) then
     362            numprh=minloc(abs(xsec18*100.0-akz),dim=1)
     363          endif
     364          help=zsec4(nxfield*(ny-j-1)+i+1)
     365          if(i.lt.i180) then
     366            qvh(i180+i,j,numprh,n)=help
     367          else
     368            qvh(i-i180,j,numprh,n)=help
     369          endif
     370        endif
     371        if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
     372! SURFACE PRESSURE
     373          help=zsec4(nxfield*(ny-j-1)+i+1)
     374          if(i.lt.i180) then
     375            ps(i180+i,j,1,n)=help
     376          else
     377            ps(i-i180,j,1,n)=help
     378          endif
     379        endif
     380        if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
     381! W VELOCITY
     382          if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
     383          help=zsec4(nxfield*(ny-j-1)+i+1)
     384          if(i.lt.i180) then
     385            wwh(i180+i,j,numpw)=help
     386          else
     387            wwh(i-i180,j,numpw)=help
     388          endif
     389        endif
     390        if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
     391! SNOW DEPTH
     392          help=zsec4(nxfield*(ny-j-1)+i+1)
     393          if(i.lt.i180) then
     394            sd(i180+i,j,1,n)=help
     395          else
     396            sd(i-i180,j,1,n)=help
     397          endif
     398        endif
     399        if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
     400! MEAN SEA LEVEL PRESSURE
     401          help=zsec4(nxfield*(ny-j-1)+i+1)
     402          if(i.lt.i180) then
     403            msl(i180+i,j,1,n)=help
     404          else
     405            msl(i-i180,j,1,n)=help
     406          endif
     407        endif
     408        if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
     409! TOTAL CLOUD COVER
     410          help=zsec4(nxfield*(ny-j-1)+i+1)
     411          if(i.lt.i180) then
     412            tcc(i180+i,j,1,n)=help
     413          else
     414            tcc(i-i180,j,1,n)=help
     415          endif
     416        endif
     417        if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
     418             (nint(xsec18).eq.10)) then
     419! 10 M U VELOCITY
     420          help=zsec4(nxfield*(ny-j-1)+i+1)
     421          if(i.lt.i180) then
     422            u10(i180+i,j,1,n)=help
     423          else
     424            u10(i-i180,j,1,n)=help
     425          endif
     426        endif
     427        if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
     428             (nint(xsec18).eq.10)) then
     429! 10 M V VELOCITY
     430          help=zsec4(nxfield*(ny-j-1)+i+1)
     431          if(i.lt.i180) then
     432            v10(i180+i,j,1,n)=help
     433          else
     434            v10(i-i180,j,1,n)=help
     435          endif
     436        endif
     437        if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
     438             (nint(xsec18).eq.2)) then
     439! 2 M TEMPERATURE
     440          help=zsec4(nxfield*(ny-j-1)+i+1)
     441          if(i.lt.i180) then
     442            tt2(i180+i,j,1,n)=help
     443          else
     444            tt2(i-i180,j,1,n)=help
     445          endif
     446        endif
     447        if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
     448             (nint(xsec18).eq.2)) then
     449! 2 M DEW POINT TEMPERATURE
     450          help=zsec4(nxfield*(ny-j-1)+i+1)
     451          if(i.lt.i180) then
     452            td2(i180+i,j,1,n)=help
     453          else
     454            td2(i-i180,j,1,n)=help
     455          endif
     456        endif
     457        if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
     458! LARGE SCALE PREC.
     459          help=zsec4(nxfield*(ny-j-1)+i+1)
     460          if(i.lt.i180) then
     461            lsprec(i180+i,j,1,n)=help
     462          else
     463            lsprec(i-i180,j,1,n)=help
     464          endif
     465        endif
     466        if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
     467! CONVECTIVE PREC.
     468          help=zsec4(nxfield*(ny-j-1)+i+1)
     469          if(i.lt.i180) then
     470            convprec(i180+i,j,1,n)=help
     471          else
     472            convprec(i-i180,j,1,n)=help
     473          endif
     474        endif
     475        if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
     476! TOPOGRAPHY
     477          help=zsec4(nxfield*(ny-j-1)+i+1)
     478          if(i.lt.i180) then
     479            oro(i180+i,j)=help
     480            excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     481          else
     482            oro(i-i180,j)=help
     483            excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     484          endif
     485        endif
     486        if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
     487! LAND SEA MASK
     488          help=zsec4(nxfield*(ny-j-1)+i+1)
     489          if(i.lt.i180) then
     490            lsm(i180+i,j)=help
     491          else
     492            lsm(i-i180,j)=help
     493          endif
     494        endif
     495        if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
     496! MIXING HEIGHT
     497          help=zsec4(nxfield*(ny-j-1)+i+1)
     498          if(i.lt.i180) then
     499            hmix(i180+i,j,1,n)=help
     500          else
     501            hmix(i-i180,j,1,n)=help
     502          endif
     503        endif
     504        if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
     505             (nint(xsec18).eq.02)) then
     506! 2 M RELATIVE HUMIDITY
     507          help=zsec4(nxfield*(ny-j-1)+i+1)
     508          if(i.lt.i180) then
     509            qvh2(i180+i,j)=help
     510          else
     511            qvh2(i-i180,j)=help
     512          endif
     513        endif
     514        if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
     515! TEMPERATURE LOWEST SIGMA LEVEL
     516          help=zsec4(nxfield*(ny-j-1)+i+1)
     517          if(i.lt.i180) then
     518            tlev1(i180+i,j)=help
     519          else
     520            tlev1(i-i180,j)=help
     521          endif
     522        endif
     523        if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
     524! U VELOCITY LOWEST SIGMA LEVEL
     525          help=zsec4(nxfield*(ny-j-1)+i+1)
     526          if(i.lt.i180) then
     527            ulev1(i180+i,j)=help
     528          else
     529            ulev1(i-i180,j)=help
     530          endif
     531        endif
     532        if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
     533! V VELOCITY LOWEST SIGMA LEVEL
     534          help=zsec4(nxfield*(ny-j-1)+i+1)
     535          if(i.lt.i180) then
     536            vlev1(i180+i,j)=help
     537          else
     538            vlev1(i-i180,j)=help
     539          endif
     540        endif
    559541! SEC & IP 12/2018 read GFS clouds
    560       if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg]
    561          if((i.eq.0).and.(j.eq.0)) then
    562             do ii=1,nuvz
    563               if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
    564             end do
    565         endif
    566         help=zsec4(nxfield*(ny-j-1)+i+1)
    567         if(i.le.i180) then
    568           clwch(i179+i,j,numpclwch,n)=help
    569         else
    570           clwch(i-i181,j,numpclwch,n)=help
    571         endif
    572         readclouds=.true.
    573         sumclouds=.true.
    574 !        readclouds=.false.
    575 !       sumclouds=.false.
    576       endif
    577 
    578 
     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
    579558    end do
    580   end do
    581559
    582560  endif
    583561
    584562  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
    585   ! NCEP ISOBARIC LEVELS
     563! NCEP ISOBARIC LEVELS
    586564    iumax=iumax+1
    587565  endif
    … …  
    589567  call grib_release(igrib)
    590568  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    591   !
    592   ! CLOSING OF INPUT DATA FILE
    593   !
    594 
    595   !HSO close grib file
    596 50   continue
     569!
     570! CLOSING OF INPUT DATA FILE
     571!
     572
     573!HSO close grib file
     57450 continue
    597575  call grib_close_file(ifile)
    598576
    599   ! SENS. HEAT FLUX
     577! SENS. HEAT FLUX
    600578  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
    601579  hflswitch=.false.    ! Heat flux not available
    602   ! SOLAR RADIATIVE FLUXES
     580! SOLAR RADIATIVE FLUXES
    603581  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
    604   ! EW SURFACE STRESS
     582! EW SURFACE STRESS
    605583  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
    606   ! NS SURFACE STRESS
     584! NS SURFACE STRESS
    607585  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
    608586  strswitch=.false.    ! stress not available
    609587
    610   ! CONVERT TP TO LSP (GRIB2 only)
     588! CONVERT TP TO LSP (GRIB2 only)
    611589  if (gribVer.eq.2) then
    612     do j=0,nymin1
    613     do i=0,nxfield-1
    614      if(i.le.i180) then
    615      if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
    616          lsprec(i179+i,j,1,n)= &
    617               lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
    618      else
    619          lsprec(i179+i,j,1,n)=0
    620      endif
    621      else
    622      if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
    623           lsprec(i-i181,j,1,n)= &
    624                lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
    625      else
    626           lsprec(i-i181,j,1,n)=0
    627      endif
    628      endif
    629     enddo
    630     enddo
    631   endif
    632   !HSO end edits
    633 
    634 
    635   ! TRANSFORM RH TO SPECIFIC HUMIDITY
     590    lsprec(0:nxfield-1, 0:nymin1, 1, n) = max( 0.0 , lsprec(0:nxfield-1,0:nymin1,1,n)-convprec(0:nxfield-1,0:nymin1,1,n) )
     591  endif
     592!HSO end edits
     593
     594
     595! TRANSFORM RH TO SPECIFIC HUMIDITY
    636596
    637597  do j=0,ny-1
    … …  
    640600        help=qvh(i,j,k,n)
    641601        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
    642606        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
    643607        elev=ew(temp)*help/100.0
    … …  
    647611  end do
    648612
    649   ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
    650   ! USING BOLTON'S (1980) FORMULA
    651   ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
     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
    652616
    653617  do j=0,ny-1
    654618    do i=0,nxfield-1
    655         help=qvh2(i,j)
    656         temp=tt2(i,j,1,n)
    657         elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
    658         td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
    659         if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
     619      help=qvh2(i,j)
     620      temp=tt2(i,j,1,n)
     621      if (temp .eq. 0.0) then
     622        write (*, *) i, j, n
     623        temp = 273.0
     624      endif
     625      elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
     626      td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
     627      if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
    660628    end do
    661629  end do
    … …  
    671639
    672640
    673   ! For global fields, assign the leftmost data column also to the rightmost
    674   ! data column; if required, shift whole grid by nxshift grid points
    675   !*************************************************************************
     641! For global fields, assign the leftmost data column also to the rightmost
     642! data column; if required, shift whole grid by nxshift grid points
     643!*************************************************************************
    676644
    677645  if (xglobal) then
    … …  
    709677  do i=0,nxmin1
    710678    do j=0,nymin1
    711   ! Convert precip. from mm/s -> mm/hour
     679! Convert precip. from mm/s -> mm/hour
    712680      convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
    713681      lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
    … …  
    717685
    718686  if ((.not.hflswitch).or.(.not.strswitch)) then
    719   !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
    720   !    +  wfname(indj)
    721 
    722   ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
    723   !***************************************************************************
     687!  write(*,*) 'WARNING: No flux data contained in GRIB file ',
     688!    +  wfname(indj)
     689
     690! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
     691!***************************************************************************
    724692
    725693    do i=0,nxmin1
    … …  
    741709
    742710  return
    743 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     711888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    744712  write(*,*) ' #### ',wfname(indj),'                    #### '
    745713  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
    746714  stop 'Execution terminated'
    747 999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     715999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    748716  write(*,*) ' #### ',wfname(indj),'                    #### '
    749717  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG