Changeset 9ca6e38 in flexpart.git


Ignore:
Timestamp:
Mar 15, 2021, 9:48:30 AM (3 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
GFS_025, dev
Children:
aa939a9
Parents:
467460a
Message:

minor edits

Location:
src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • src/calcmatrix.f90

    r467460a r9ca6e38  
    6767
    6868  phconv(1) = psconv
    69   ! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
    70   do kuvz = 2,nuvz
    71     k = kuvz-1
    72     if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     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
    7373      pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
    7474      phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
    75     else
    76       phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
    77     endif
    78     dpr(k) = phconv(k) - phconv(kuvz)
    79     qsconv(k) = f_qvsat( pconv(k), tconv(k) )
     75      dpr(k) = phconv(k) - phconv(kuvz)
     76      qsconv(k) = f_qvsat( pconv(k), tconv(k) )
     77    end do
     78  else
     79! JMA / SH Bugfix phconv was set in loop with access to undefined pconv(nuvz)
     80    phconv(2:nuvz-1) = 0.5*(pconv(2:nuvz-1)+pconv(1:nuvz-2))
     81    phconv(nuvz) = pconv(nuvz-1)
     82    dpr(1:nuvz-1) = phconv(1:nuvz-1)-phconv(2:nuvz)
    8083
    81   ! initialize mass fractions
    82     do kk=1,nconvlev
    83       fmassfrac(k,kk)=0.
     84    do k = 1,nuvz-1
     85      qsconv(k) = f_qvsat( pconv(k), tconv(k) )
    8486    end do
    85   end do
    86 
     87  end if
     88   
     89! initialize mass fractions
     90  fmassfrac(1:nuvz-1,1:nconvlev)=0.0
    8791
    8892  !note that Emanuel says it is important
     
    9397  !******************
    9498
    95     cbmfold = cbmf
    96   ! Convert pressures to hPa, as required by Emanuel scheme
    97   !********************************************************
     99  cbmfold = cbmf
     100! Convert pressures to hPa, as required by Emanuel scheme
     101!********************************************************
    98102!!$    do k=1,nconvlev     !old
    99     do k=1,nconvlev+1      !bugfix
    100       pconv_hpa(k)=pconv(k)/100.
    101       phconv_hpa(k)=phconv(k)/100.
     103  do k=1,nconvlev+1      !bugfix
     104    pconv_hpa(k)=pconv(k)/100.
     105    phconv_hpa(k)=phconv(k)/100.
     106  end do
     107  phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
     108  call convect(nconvlevmax, nconvlev, delt, iflag, &
     109       precip, wd, tprime, qprime, cbmf)
     110
     111! do not update fmassfrac and cloudbase massflux
     112! if no convection takes place or
     113! if a CFL criterion is violated in convect43c.f
     114  if (iflag .ne. 1 .and. iflag .ne. 4) then
     115    cbmf=cbmfold
     116    goto 200
     117  endif
     118
     119! do not update fmassfrac and cloudbase massflux
     120! if the old and the new cloud base mass
     121! fluxes are zero
     122  if (cbmf.le.0..and.cbmfold.le.0.) then
     123    cbmf=cbmfold
     124    goto 200
     125  endif
     126
     127! Update fmassfrac
     128! account for mass displaced from level k to level k
     129
     130  lconv = .true.
     131  do k=1,nconvtop
     132    rlevmass = dpr(k)/ga
     133    summe = 0.
     134    do kk=1,nconvtop
     135      fmassfrac(k,kk) = delt*fmass(k,kk)
     136      summe = summe + fmassfrac(k,kk)
    102137    end do
    103     phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
    104     call convect(nconvlevmax, nconvlev, delt, iflag, &
    105          precip, wd, tprime, qprime, cbmf)
     138    fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
     139  end do
    106140
    107   ! do not update fmassfrac and cloudbase massflux
    108   ! if no convection takes place or
    109   ! if a CFL criterion is violated in convect43c.f
    110    if (iflag .ne. 1 .and. iflag .ne. 4) then
    111      cbmf=cbmfold
    112      goto 200
    113    endif
    114 
    115   ! do not update fmassfrac and cloudbase massflux
    116   ! if the old and the new cloud base mass
    117   ! fluxes are zero
    118    if (cbmf.le.0..and.cbmfold.le.0.) then
    119      cbmf=cbmfold
    120      goto 200
    121    endif
    122 
    123   ! Update fmassfrac
    124   ! account for mass displaced from level k to level k
    125 
    126    lconv = .true.
    127     do k=1,nconvtop
    128       rlevmass = dpr(k)/ga
    129       summe = 0.
    130       do kk=1,nconvtop
    131         fmassfrac(k,kk) = delt*fmass(k,kk)
    132         summe = summe + fmassfrac(k,kk)
    133       end do
    134       fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
    135     end do
    136 
    137 200   continue
     141200 continue
    138142
    139143end subroutine calcmatrix
  • src/readwind_gfs.f90

    r467460a r9ca6e38  
    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,valSurf,discipl
    60   !HSO end edits
     59  integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
     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
     
    7474  integer :: i180
    7575
    76   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    77   !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
     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)
     
    8888  logical :: hflswitch,strswitch
    8989
    90   !HSO  for grib api error messages
     90!HSO  for grib api error messages
    9191  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    9292  character(len=20) :: gribFunction = 'readwind_gfs'
     
    101101
    102102
    103   ! OPENING OF DATA FILE (GRIB CODE)
    104 
    105   !HSO
     103! OPENING OF DATA FILE (GRIB CODE)
     104
     105!HSO
    106106  call grib_open_file(ifile,path(3)(1:length(3)) &
    107          //trim(wfname(indj)),'r',iret)
     107       //trim(wfname(indj)),'r',iret)
    108108  if (iret.ne.GRIB_SUCCESS) then
    109109    goto 888   ! ERROR DETECTED
    110110  endif
    111   !turn on support for multi fields messages
     111!turn on support for multi fields messages
    112112  call grib_multi_support_on
    113113
     
    119119  numpclwch=0
    120120  ifield=0
    121 10   ifield=ifield+1
    122   !
    123   ! GET NEXT FIELDS
    124   !
     12110 ifield=ifield+1
     122!
     123! GET NEXT FIELDS
     124!
    125125  call grib_new_from_file(ifile,igrib,iret)
    126126  if (iret.eq.GRIB_END_OF_FILE)  then
     
    130130  endif
    131131
    132   !first see if we read GRIB1 or GRIB2
     132!first see if we read GRIB1 or GRIB2
    133133  call grib_get_int(igrib,'editionNumber',gribVer,iret)
    134134!  call grib_check(iret,gribFunction,gribErrorMsg)
     
    136136  if (gribVer.eq.1) then ! GRIB Edition 1
    137137
    138   !read the grib1 identifiers
    139   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    140   call grib_check(iret,gribFunction,gribErrorMsg)
    141   call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
    142   call grib_check(iret,gribFunction,gribErrorMsg)
    143   call grib_get_int(igrib,'level',isec1(8),iret)
    144   call grib_check(iret,gribFunction,gribErrorMsg)
    145   !JMA / SH: isec1(8) not evaluated any more below
    146   !b/c with GRIB 2 this may be a real variable
    147   xsec18 = real(isec1(8))
     138!read the grib1 identifiers
     139
     140    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     141    call grib_check(iret,gribFunction,gribErrorMsg)
     142    call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
     143    call grib_check(iret,gribFunction,gribErrorMsg)
     144    call grib_get_int(igrib,'level',isec1(8),iret)
     145    call grib_check(iret,gribFunction,gribErrorMsg)
     146!JMA / SH: isec1(8) not evaluated any more below
     147!b/c with GRIB 2 this may be a real variable
     148    xsec18 = real(isec1(8))
    148149
    149150  else ! GRIB Edition 2
    150151
    151   !read the grib2 identifiers
    152   call grib_get_string(igrib,'shortName',shortname,iret)
    153 
    154   call grib_get_int(igrib,'discipline',discipl,iret)
    155   call grib_check(iret,gribFunction,gribErrorMsg)
    156   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    157   call grib_check(iret,gribFunction,gribErrorMsg)
    158   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    159   call grib_check(iret,gribFunction,gribErrorMsg)
    160   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    161   call grib_check(iret,gribFunction,gribErrorMsg)
    162   call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
    163        valSurf,iret)
    164   call grib_check(iret,gribFunction,gribErrorMsg)
    165  
    166 !  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
    167   !convert to grib1 identifiers
    168   isec1(6:8)=-1
    169   xsec18  =-1.0
    170   !JMA / SH: isec1(8) not evaluated any more below
    171   !b/c with GRIB 2 this may be a real variable
    172   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
    173     isec1(6)=11          ! indicatorOfParameter
    174     isec1(7)=100         ! indicatorOfTypeOfLevel
    175     xsec18=valSurf/100.0 ! level, convert to hPa
    176   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
    177     isec1(6)=33          ! indicatorOfParameter
    178     isec1(7)=100         ! indicatorOfTypeOfLevel
    179     xsec18=valSurf/100.0 ! level, convert to hPa
    180   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
    181     isec1(6)=34          ! indicatorOfParameter
    182     isec1(7)=100         ! indicatorOfTypeOfLevel
    183     xsec18=valSurf/100.0 ! level, convert to hPa
    184   elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
    185     isec1(6)=39          ! indicatorOfParameter
    186     isec1(7)=100         ! indicatorOfTypeOfLevel
    187     xsec18=valSurf/100.0 ! level, convert to hPa
    188   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
    189     isec1(6)=52          ! indicatorOfParameter
    190     isec1(7)=100         ! indicatorOfTypeOfLevel
    191     xsec18=valSurf/100.0 ! level, convert to hPa
    192   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
    193     isec1(6)=52          ! indicatorOfParameter
    194     isec1(7)=105         ! indicatorOfTypeOfLevel
    195     xsec18=real(2)
    196   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
    197     isec1(6)=11          ! indicatorOfParameter
    198     isec1(7)=105         ! indicatorOfTypeOfLevel
    199     xsec18=real(2)
    200   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
    201     isec1(6)=33          ! indicatorOfParameter
    202     isec1(7)=105         ! indicatorOfTypeOfLevel
    203     xsec18=real(10)
    204   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
    205     isec1(6)=34          ! indicatorOfParameter
    206     isec1(7)=105         ! indicatorOfTypeOfLevel
    207     xsec18=real(10)
    208   elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
    209     isec1(6)=153         ! indicatorOfParameter
    210     isec1(7)=100         ! indicatorOfTypeOfLevel
    211     xsec18=valSurf/100.0 ! level, convert to hPa
    212   elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
    213     isec1(6)=2           ! indicatorOfParameter
    214     isec1(7)=102         ! indicatorOfTypeOfLevel
    215     xsec18=real(0)
    216   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
    217     isec1(6)=1           ! indicatorOfParameter
    218     isec1(7)=1           ! indicatorOfTypeOfLevel
    219     xsec18=real(0)
    220   elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
    221     isec1(6)=66          ! indicatorOfParameter
    222     isec1(7)=1           ! indicatorOfTypeOfLevel
    223     xsec18=real(0)
    224   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
    225     isec1(6)=11          ! indicatorOfParameter
    226     isec1(7)=107         ! indicatorOfTypeOfLevel
    227     xsec18=0.995         ! lowest sigma level
    228   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
    229     isec1(6)=33          ! indicatorOfParameter
    230     isec1(7)=107         ! indicatorOfTypeOfLevel
    231     xsec18=0.995         ! lowest sigma level
    232   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
    233     isec1(6)=34          ! indicatorOfParameter
    234     isec1(7)=107         ! indicatorOfTypeOfLevel
    235     xsec18=0.995         ! lowest sigma level
    236   elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
    237     isec1(6)=7           ! indicatorOfParameter
    238     isec1(7)=1           ! indicatorOfTypeOfLevel
    239     xsec18=real(0)
    240   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
    241        .and.(discipl.eq.2)) then ! LSM
    242     isec1(6)=81          ! indicatorOfParameter
    243     isec1(7)=1           ! indicatorOfTypeOfLevel
    244     xsec18=real(0)
    245   elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
    246     isec1(6)=221         ! indicatorOfParameter
    247     isec1(7)=1           ! indicatorOfTypeOfLevel
    248     xsec18=real(0)
    249   elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
    250     isec1(6)=62          ! indicatorOfParameter
    251     isec1(7)=1           ! indicatorOfTypeOfLevel
    252     xsec18=real(0)
    253   elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
    254     isec1(6)=63          ! indicatorOfParameter
    255     isec1(7)=1           ! indicatorOfTypeOfLevel
    256     xsec18=real(0)
    257   endif
     152!read the grib2 identifiers
     153   
     154    call grib_get_string(igrib,'shortName',shortname,iret)
     155
     156    call grib_get_int(igrib,'discipline',discipl,iret)
     157    call grib_check(iret,gribFunction,gribErrorMsg)
     158    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     159    call grib_check(iret,gribFunction,gribErrorMsg)
     160    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     161    call grib_check(iret,gribFunction,gribErrorMsg)
     162    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     163    call grib_check(iret,gribFunction,gribErrorMsg)
     164!
     165    call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
     166         valSurf,iret)
     167    call grib_check(iret,gribFunction,gribErrorMsg)
     168
     169!    write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
     170!convert to grib1 identifiers
     171    isec1(6:8)=-1
     172    xsec18  =-1.0
     173!JMA / SH: isec1(8) not evaluated any more below
     174!b/c with GRIB 2 this may be a real variable
     175    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
     176      isec1(6)=11          ! indicatorOfParameter
     177      isec1(7)=100         ! indicatorOfTypeOfLevel
     178      xsec18=valSurf/100.0 ! level, convert to hPa
     179    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
     180      isec1(6)=33          ! indicatorOfParameter
     181      isec1(7)=100         ! indicatorOfTypeOfLevel
     182      xsec18=valSurf/100.0 ! level, convert to hPa
     183    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
     184      isec1(6)=34          ! indicatorOfParameter
     185      isec1(7)=100         ! indicatorOfTypeOfLevel
     186      xsec18=valSurf/100.0 ! level, convert to hPa
     187    elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
     188      isec1(6)=39          ! indicatorOfParameter
     189      isec1(7)=100         ! indicatorOfTypeOfLevel
     190      xsec18=valSurf/100.0 ! level, convert to hPa
     191    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
     192      isec1(6)=52          ! indicatorOfParameter
     193      isec1(7)=100         ! indicatorOfTypeOfLevel
     194      xsec18=valSurf/100.0 ! level, convert to hPa
     195    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
     196      isec1(6)=52          ! indicatorOfParameter
     197      isec1(7)=105         ! indicatorOfTypeOfLevel
     198      xsec18=real(2)
     199    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
     200      isec1(6)=11          ! indicatorOfParameter
     201      isec1(7)=105         ! indicatorOfTypeOfLevel
     202      xsec18=real(2)
     203    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
     204      isec1(6)=33          ! indicatorOfParameter
     205      isec1(7)=105         ! indicatorOfTypeOfLevel
     206      xsec18=real(10)
     207    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
     208      isec1(6)=34          ! indicatorOfParameter
     209      isec1(7)=105         ! indicatorOfTypeOfLevel
     210      xsec18=real(10)
     211    elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
     212      isec1(6)=153         ! indicatorOfParameter
     213      isec1(7)=100         ! indicatorOfTypeOfLevel
     214      xsec18=valSurf/100.0 ! level, convert to hPa
     215    elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
     216      isec1(6)=2           ! indicatorOfParameter
     217      isec1(7)=102         ! indicatorOfTypeOfLevel
     218      xsec18=real(0)
     219    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
     220      isec1(6)=1           ! indicatorOfParameter
     221      isec1(7)=1           ! indicatorOfTypeOfLevel
     222      xsec18=real(0)
     223    elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
     224      isec1(6)=66          ! indicatorOfParameter
     225      isec1(7)=1           ! indicatorOfTypeOfLevel
     226      xsec18=real(0)
     227    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
     228      isec1(6)=11          ! indicatorOfParameter
     229      isec1(7)=107         ! indicatorOfTypeOfLevel
     230      xsec18=0.995         ! lowest sigma level
     231    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
     232      isec1(6)=33          ! indicatorOfParameter
     233      isec1(7)=107         ! indicatorOfTypeOfLevel
     234      xsec18=0.995         ! lowest sigma level
     235    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
     236      isec1(6)=34          ! indicatorOfParameter
     237      isec1(7)=107         ! indicatorOfTypeOfLevel
     238      xsec18=0.995         ! lowest sigma level
     239    elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
     240      isec1(6)=7           ! indicatorOfParameter
     241      isec1(7)=1           ! indicatorOfTypeOfLevel
     242      xsec18=real(0)
     243    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
     244         .and.(discipl.eq.2)) then ! LSM
     245      isec1(6)=81          ! indicatorOfParameter
     246      isec1(7)=1           ! indicatorOfTypeOfLevel
     247      xsec18=real(0)
     248    elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
     249      isec1(6)=221         ! indicatorOfParameter
     250      isec1(7)=1           ! indicatorOfTypeOfLevel
     251      xsec18=real(0)
     252    elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
     253      isec1(6)=62          ! indicatorOfParameter
     254      isec1(7)=1           ! indicatorOfTypeOfLevel
     255      xsec18=real(0)
     256    elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
     257      isec1(6)=63          ! indicatorOfParameter
     258      isec1(7)=1           ! indicatorOfTypeOfLevel
     259      xsec18=real(0)
     260    endif
    258261
    259262  endif ! gribVer
    260263
    261264  if (isec1(6).ne.-1) then
    262   !  get the size and data of the values array
     265!  get the size and data of the values array
    263266    call grib_get_real4_array(igrib,'values',zsec4,iret)
    264267    call grib_check(iret,gribFunction,gribErrorMsg)
     
    267270  if(ifield.eq.1) then
    268271
    269   !get the required fields from section 2
    270   !store compatible to gribex input
    271   call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    272        isec2(2),iret)
    273   call grib_check(iret,gribFunction,gribErrorMsg)
    274   call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
    275        isec2(3),iret)
    276   call grib_check(iret,gribFunction,gribErrorMsg)
    277   call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    278        xauxin,iret)
    279   call grib_check(iret,gribFunction,gribErrorMsg)
    280   call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
    281        yauxin,iret)
    282   call grib_check(iret,gribFunction,gribErrorMsg)
    283   xaux=xauxin+real(nxshift)*dx
    284   yaux=yauxin
    285 
    286   ! 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
    287290
    288291    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
     
    296299    if(yaux0.lt.0.) yaux0=yaux0+360.
    297300    if(abs(xaux-xaux0).gt.eps) then
    298          write (*, *) xaux, xaux0
    299          stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
     301      write (*, *) xaux, xaux0
     302      stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    300303    endif
    301304    if(abs(yaux-yaux0).gt.eps) then
    302          write (*, *) yaux, yaux0
    303          stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     305      write (*, *) yaux, yaux0
     306      stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
    304307    end if
    305308  endif
    306   !HSO end of edits
     309!HSO end of edits
    307310
    308311  i180=nint(180./dx)
     
    312315! write (*, *) 'nxfield: ', nxfield, i180
    313316
    314   do j=0,nymin1
    315     do i=0,nxfield-1
    316       if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
     317    do j=0,nymin1
     318      do i=0,nxfield-1
     319        if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
    317320! TEMPERATURE
    318         if((i.eq.0).and.(j.eq.0)) then
    319           do ii=1,nuvz+1
    320              write(*,*) 'xsec18, akz(ii), xsec18*100.0-akz(ii), spacing(akz(ii), spacing(xsec18*100.0)',&
    321                   & xsec18, akz(ii), xsec18*100.0-akz(ii), spacing(akz(ii)), spacing(xsec18*100.0)
    322              
    323             if (abs(xsec18*100.0-akz(ii)) < &
    324                  10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then
    325               numpt=ii
    326             end if
    327           end do
    328         endif
    329         help=zsec4(nxfield*(ny-j-1)+i+1)
    330         if (help.eq.0) then
    331           write (*, *) 'i, j: ', i, j
    332           stop 'help == 0.0'
    333         endif
    334         if(i.lt.i180) then
    335 ! ESO: dbg out-of-bounds issue
    336           if (numpt < 1) then
    337             write(*,*) 'nuvzmax, nuvz', nuvzmax, nuvz
    338             write(*,*) 'i180+i,j,numpt,n', i180+i,j,numpt,n
    339           end if
    340           tth(i180+i,j,numpt,n)=help
    341         else
    342           tth(i-i180,j,numpt,n)=help
    343         endif
    344       endif
    345       if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
    346   ! U VELOCITY
    347          if((i.eq.0).and.(j.eq.0)) then
    348             do ii=1,nuvz
    349              if (abs(xsec18*100.0-akz(ii)) < &
    350                 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then
    351                numpu=ii
    352              end if
    353             end do
    354         endif
    355         help=zsec4(nxfield*(ny-j-1)+i+1)
    356         if(i.lt.i180) then
    357           uuh(i180+i,j,numpu)=help
    358         else
    359           uuh(i-i180,j,numpu)=help
    360         endif
    361       endif
    362       if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
    363   ! V VELOCITY
    364          if((i.eq.0).and.(j.eq.0)) then
    365             do ii=1,nuvz
    366              if (abs(xsec18*100.0-akz(ii)) < &
    367                 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then
    368                 numpv=ii
    369              end if
    370             end do
    371         endif
    372         help=zsec4(nxfield*(ny-j-1)+i+1)
    373         if(i.lt.i180) then
    374           vvh(i180+i,j,numpv)=help
    375         else
    376           vvh(i-i180,j,numpv)=help
    377         endif
    378       endif
    379       if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
    380   ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
    381          if((i.eq.0).and.(j.eq.0)) then
    382             do ii=1,nuvz
    383              if (abs(xsec18*100.0-akz(ii)) < &
    384                 10.0*max(spacing(akz(ii)),spacing(xsec18*100.0))) then
    385                numprh=ii
    386              end if
    387             end do
    388         endif
    389         help=zsec4(nxfield*(ny-j-1)+i+1)
    390         if(i.lt.i180) then
    391           qvh(i180+i,j,numprh,n)=help
    392         else
    393           qvh(i-i180,j,numprh,n)=help
    394         endif
    395       endif
    396       if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
    397   ! SURFACE PRESSURE
    398         help=zsec4(nxfield*(ny-j-1)+i+1)
    399         if(i.lt.i180) then
    400           ps(i180+i,j,1,n)=help
    401         else
    402           ps(i-i180,j,1,n)=help
    403         endif
    404       endif
    405       if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
    406   ! W VELOCITY
    407          if((i.eq.0).and.(j.eq.0)) then
    408             do ii=1,nuvz
    409               if ((xsec18*100.0).eq.akz(ii)) numpw=ii
    410             end do
    411         endif
    412         help=zsec4(nxfield*(ny-j-1)+i+1)
    413         if(i.lt.i180) then
    414           wwh(i180+i,j,numpw)=help
    415         else
    416           wwh(i-i180,j,numpw)=help
    417         endif
    418       endif
    419       if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
    420   ! SNOW DEPTH
    421         help=zsec4(nxfield*(ny-j-1)+i+1)
    422         if(i.lt.i180) then
    423           sd(i180+i,j,1,n)=help
    424         else
    425           sd(i-i180,j,1,n)=help
    426         endif
    427       endif
    428       if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
    429   ! MEAN SEA LEVEL PRESSURE
    430         help=zsec4(nxfield*(ny-j-1)+i+1)
    431         if(i.lt.i180) then
    432           msl(i180+i,j,1,n)=help
    433         else
    434           msl(i-i180,j,1,n)=help
    435         endif
    436       endif
    437       if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
    438   ! TOTAL CLOUD COVER
    439         help=zsec4(nxfield*(ny-j-1)+i+1)
    440         if(i.lt.i180) then
    441           tcc(i180+i,j,1,n)=help
    442         else
    443           tcc(i-i180,j,1,n)=help
    444         endif
    445       endif
    446       if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
    447          (nint(xsec18).eq.10)) then
    448   ! 10 M U VELOCITY
    449         help=zsec4(nxfield*(ny-j-1)+i+1)
    450         if(i.lt.i180) then
    451           u10(i180+i,j,1,n)=help
    452         else
    453           u10(i-i180,j,1,n)=help
    454         endif
    455       endif
    456       if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
    457          (nint(xsec18).eq.10)) then
    458   ! 10 M V VELOCITY
    459         help=zsec4(nxfield*(ny-j-1)+i+1)
    460         if(i.lt.i180) then
    461           v10(i180+i,j,1,n)=help
    462         else
    463           v10(i-i180,j,1,n)=help
    464         endif
    465       endif
    466       if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
    467          (nint(xsec18).eq.2)) then
    468   ! 2 M TEMPERATURE
    469         help=zsec4(nxfield*(ny-j-1)+i+1)
    470         if(i.lt.i180) then
    471           tt2(i180+i,j,1,n)=help
    472         else
    473           tt2(i-i180,j,1,n)=help
    474         endif
    475       endif
    476       if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
    477          (nint(xsec18).eq.2)) then
    478   ! 2 M DEW POINT TEMPERATURE
    479         help=zsec4(nxfield*(ny-j-1)+i+1)
    480         if(i.lt.i180) then
    481           td2(i180+i,j,1,n)=help
    482         else
    483           td2(i-i180,j,1,n)=help
    484         endif
    485       endif
    486       if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
    487   ! LARGE SCALE PREC.
    488         help=zsec4(nxfield*(ny-j-1)+i+1)
    489         if(i.lt.i180) then
    490           lsprec(i180+i,j,1,n)=help
    491         else
    492           lsprec(i-i180,j,1,n)=help
    493         endif
    494       endif
    495       if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
    496   ! CONVECTIVE PREC.
    497         help=zsec4(nxfield*(ny-j-1)+i+1)
    498         if(i.lt.i180) then
    499           convprec(i180+i,j,1,n)=help
    500         else
    501           convprec(i-i180,j,1,n)=help
    502         endif
    503       endif
    504       if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
    505   ! TOPOGRAPHY
    506         help=zsec4(nxfield*(ny-j-1)+i+1)
    507         if(i.lt.i180) then
    508           oro(i180+i,j)=help
    509           excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    510         else
    511           oro(i-i180,j)=help
    512           excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    513         endif
    514       endif
    515       if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
    516   ! LAND SEA MASK
    517         help=zsec4(nxfield*(ny-j-1)+i+1)
    518         if(i.lt.i180) then
    519           lsm(i180+i,j)=help
    520         else
    521           lsm(i-i180,j)=help
    522         endif
    523       endif
    524       if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
    525   ! MIXING HEIGHT
    526         help=zsec4(nxfield*(ny-j-1)+i+1)
    527         if(i.lt.i180) then
    528           hmix(i180+i,j,1,n)=help
    529         else
    530           hmix(i-i180,j,1,n)=help
    531         endif
    532       endif
    533       if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
    534          (nint(xsec18).eq.02)) then
    535   ! 2 M RELATIVE HUMIDITY
    536         help=zsec4(nxfield*(ny-j-1)+i+1)
    537         if(i.lt.i180) then
    538           qvh2(i180+i,j)=help
    539         else
    540           qvh2(i-i180,j)=help
    541         endif
    542       endif
    543       if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
    544   ! TEMPERATURE LOWEST SIGMA LEVEL
    545         help=zsec4(nxfield*(ny-j-1)+i+1)
    546         if(i.lt.i180) then
    547           tlev1(i180+i,j)=help
    548         else
    549           tlev1(i-i180,j)=help
    550         endif
    551       endif
    552       if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
    553   ! U VELOCITY LOWEST SIGMA LEVEL
    554         help=zsec4(nxfield*(ny-j-1)+i+1)
    555         if(i.lt.i180) then
    556           ulev1(i180+i,j)=help
    557         else
    558           ulev1(i-i180,j)=help
    559         endif
    560       endif
    561       if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
    562   ! V VELOCITY LOWEST SIGMA LEVEL
    563         help=zsec4(nxfield*(ny-j-1)+i+1)
    564         if(i.lt.i180) then
    565           vlev1(i180+i,j)=help
    566         else
    567           vlev1(i-i180,j)=help
    568         endif
    569       endif
     321          if((i.eq.0).and.(j.eq.0)) then
     322            numpt=minloc(abs(xsec18*100.0-akz),dim=1)
     323          endif
     324          help=zsec4(nxfield*(ny-j-1)+i+1)
     325          if (help.eq.0) then
     326            write (*, *) 'i, j: ', i, j
     327            stop 'help == 0.0'
     328          endif
     329          if(i.lt.i180) then
     330            tth(i180+i,j,numpt,n)=help
     331          else
     332            tth(i-i180,j,numpt,n)=help
     333          endif
     334        endif
     335        if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
     336! U VELOCITY
     337          if((i.eq.0).and.(j.eq.0)) then
     338            numpu=minloc(abs(xsec18*100.0-akz),dim=1)
     339          endif
     340          help=zsec4(nxfield*(ny-j-1)+i+1)
     341          if(i.lt.i180) then
     342            uuh(i180+i,j,numpu)=help
     343          else
     344            uuh(i-i180,j,numpu)=help
     345          endif
     346        endif
     347        if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
     348! V VELOCITY
     349          if((i.eq.0).and.(j.eq.0)) then
     350            numpv=minloc(abs(xsec18*100.0-akz),dim=1)
     351          endif
     352          help=zsec4(nxfield*(ny-j-1)+i+1)
     353          if(i.lt.i180) then
     354            vvh(i180+i,j,numpv)=help
     355          else
     356            vvh(i-i180,j,numpv)=help
     357          endif
     358        endif
     359        if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
     360! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
     361          if((i.eq.0).and.(j.eq.0)) then
     362            numprh=minloc(abs(xsec18*100.0-akz),dim=1)
     363          endif
     364          help=zsec4(nxfield*(ny-j-1)+i+1)
     365          if(i.lt.i180) then
     366            qvh(i180+i,j,numprh,n)=help
     367          else
     368            qvh(i-i180,j,numprh,n)=help
     369          endif
     370        endif
     371        if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
     372! SURFACE PRESSURE
     373          help=zsec4(nxfield*(ny-j-1)+i+1)
     374          if(i.lt.i180) then
     375            ps(i180+i,j,1,n)=help
     376          else
     377            ps(i-i180,j,1,n)=help
     378          endif
     379        endif
     380        if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
     381! W VELOCITY
     382          if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
     383          help=zsec4(nxfield*(ny-j-1)+i+1)
     384          if(i.lt.i180) then
     385            wwh(i180+i,j,numpw)=help
     386          else
     387            wwh(i-i180,j,numpw)=help
     388          endif
     389        endif
     390        if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
     391! SNOW DEPTH
     392          help=zsec4(nxfield*(ny-j-1)+i+1)
     393          if(i.lt.i180) then
     394            sd(i180+i,j,1,n)=help
     395          else
     396            sd(i-i180,j,1,n)=help
     397          endif
     398        endif
     399        if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
     400! MEAN SEA LEVEL PRESSURE
     401          help=zsec4(nxfield*(ny-j-1)+i+1)
     402          if(i.lt.i180) then
     403            msl(i180+i,j,1,n)=help
     404          else
     405            msl(i-i180,j,1,n)=help
     406          endif
     407        endif
     408        if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
     409! TOTAL CLOUD COVER
     410          help=zsec4(nxfield*(ny-j-1)+i+1)
     411          if(i.lt.i180) then
     412            tcc(i180+i,j,1,n)=help
     413          else
     414            tcc(i-i180,j,1,n)=help
     415          endif
     416        endif
     417        if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
     418             (nint(xsec18).eq.10)) then
     419! 10 M U VELOCITY
     420          help=zsec4(nxfield*(ny-j-1)+i+1)
     421          if(i.lt.i180) then
     422            u10(i180+i,j,1,n)=help
     423          else
     424            u10(i-i180,j,1,n)=help
     425          endif
     426        endif
     427        if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
     428             (nint(xsec18).eq.10)) then
     429! 10 M V VELOCITY
     430          help=zsec4(nxfield*(ny-j-1)+i+1)
     431          if(i.lt.i180) then
     432            v10(i180+i,j,1,n)=help
     433          else
     434            v10(i-i180,j,1,n)=help
     435          endif
     436        endif
     437        if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
     438             (nint(xsec18).eq.2)) then
     439! 2 M TEMPERATURE
     440          help=zsec4(nxfield*(ny-j-1)+i+1)
     441          if(i.lt.i180) then
     442            tt2(i180+i,j,1,n)=help
     443          else
     444            tt2(i-i180,j,1,n)=help
     445          endif
     446        endif
     447        if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
     448             (nint(xsec18).eq.2)) then
     449! 2 M DEW POINT TEMPERATURE
     450          help=zsec4(nxfield*(ny-j-1)+i+1)
     451          if(i.lt.i180) then
     452            td2(i180+i,j,1,n)=help
     453          else
     454            td2(i-i180,j,1,n)=help
     455          endif
     456        endif
     457        if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
     458! LARGE SCALE PREC.
     459          help=zsec4(nxfield*(ny-j-1)+i+1)
     460          if(i.lt.i180) then
     461            lsprec(i180+i,j,1,n)=help
     462          else
     463            lsprec(i-i180,j,1,n)=help
     464          endif
     465        endif
     466        if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
     467! CONVECTIVE PREC.
     468          help=zsec4(nxfield*(ny-j-1)+i+1)
     469          if(i.lt.i180) then
     470            convprec(i180+i,j,1,n)=help
     471          else
     472            convprec(i-i180,j,1,n)=help
     473          endif
     474        endif
     475        if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
     476! TOPOGRAPHY
     477          help=zsec4(nxfield*(ny-j-1)+i+1)
     478          if(i.lt.i180) then
     479            oro(i180+i,j)=help
     480            excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     481          else
     482            oro(i-i180,j)=help
     483            excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     484          endif
     485        endif
     486        if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
     487! LAND SEA MASK
     488          help=zsec4(nxfield*(ny-j-1)+i+1)
     489          if(i.lt.i180) then
     490            lsm(i180+i,j)=help
     491          else
     492            lsm(i-i180,j)=help
     493          endif
     494        endif
     495        if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
     496! MIXING HEIGHT
     497          help=zsec4(nxfield*(ny-j-1)+i+1)
     498          if(i.lt.i180) then
     499            hmix(i180+i,j,1,n)=help
     500          else
     501            hmix(i-i180,j,1,n)=help
     502          endif
     503        endif
     504        if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
     505             (nint(xsec18).eq.02)) then
     506! 2 M RELATIVE HUMIDITY
     507          help=zsec4(nxfield*(ny-j-1)+i+1)
     508          if(i.lt.i180) then
     509            qvh2(i180+i,j)=help
     510          else
     511            qvh2(i-i180,j)=help
     512          endif
     513        endif
     514        if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
     515! TEMPERATURE LOWEST SIGMA LEVEL
     516          help=zsec4(nxfield*(ny-j-1)+i+1)
     517          if(i.lt.i180) then
     518            tlev1(i180+i,j)=help
     519          else
     520            tlev1(i-i180,j)=help
     521          endif
     522        endif
     523        if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
     524! U VELOCITY LOWEST SIGMA LEVEL
     525          help=zsec4(nxfield*(ny-j-1)+i+1)
     526          if(i.lt.i180) then
     527            ulev1(i180+i,j)=help
     528          else
     529            ulev1(i-i180,j)=help
     530          endif
     531        endif
     532        if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
     533! V VELOCITY LOWEST SIGMA LEVEL
     534          help=zsec4(nxfield*(ny-j-1)+i+1)
     535          if(i.lt.i180) then
     536            vlev1(i180+i,j)=help
     537          else
     538            vlev1(i-i180,j)=help
     539          endif
     540        endif
    570541! SEC & IP 12/2018 read GFS clouds
    571 ! ESO TODO: do as for hachinger version
    572       ! if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg]
    573       !    if((i.eq.0).and.(j.eq.0)) then
    574       !       do ii=1,nuvz
    575       !         if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
    576       !       end do
    577       !   endif
    578       !   help=zsec4(nxfield*(ny-j-1)+i+1)
    579       !   if(i.le.i180) then
    580       !     clwch(i179+i,j,numpclwch,n)=help
    581       !   else
    582       !     clwch(i-i181,j,numpclwch,n)=help
    583       !   endif
    584       !   readclouds=.true.
    585       !   sumclouds=.true.
    586       ! endif
    587 
    588 
     542        if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg]
     543          if((i.eq.0).and.(j.eq.0)) then
     544            numpclwch=minloc(abs(xsec18*100.0-akz),dim=1)
     545          endif
     546          help=zsec4(nxfield*(ny-j-1)+i+1)
     547          if(i.lt.i180) then
     548            clwch(i180+i,j,numpclwch,n)=help
     549          else
     550            clwch(i-i180,j,numpclwch,n)=help
     551          endif
     552          readclouds=.true.
     553          sumclouds=.true.
     554        endif
     555
     556
     557      end do
    589558    end do
    590   end do
    591559
    592560  endif
    593561
    594562  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
    595   ! NCEP ISOBARIC LEVELS
     563! NCEP ISOBARIC LEVELS
    596564    iumax=iumax+1
    597565  endif
     
    599567  call grib_release(igrib)
    600568  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    601   !
    602   ! CLOSING OF INPUT DATA FILE
    603   !
    604 
    605   !HSO close grib file
    606 50   continue
     569!
     570! CLOSING OF INPUT DATA FILE
     571!
     572
     573!HSO close grib file
     57450 continue
    607575  call grib_close_file(ifile)
    608576
    609   ! SENS. HEAT FLUX
     577! SENS. HEAT FLUX
    610578  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
    611579  hflswitch=.false.    ! Heat flux not available
    612   ! SOLAR RADIATIVE FLUXES
     580! SOLAR RADIATIVE FLUXES
    613581  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
    614   ! EW SURFACE STRESS
     582! EW SURFACE STRESS
    615583  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
    616   ! NS SURFACE STRESS
     584! NS SURFACE STRESS
    617585  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
    618586  strswitch=.false.    ! stress not available
    619587
    620   ! CONVERT TP TO LSP (GRIB2 only)
     588! CONVERT TP TO LSP (GRIB2 only)
    621589  if (gribVer.eq.2) then
    622590    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) )
    623591  endif
    624   !HSO end edits
    625 
    626 
    627   ! TRANSFORM RH TO SPECIFIC HUMIDITY
     592!HSO end edits
     593
     594
     595! TRANSFORM RH TO SPECIFIC HUMIDITY
    628596
    629597  do j=0,ny-1
     
    633601        temp=tth(i,j,k,n)
    634602        if (temp .eq. 0.0) then
    635            write (*, *) i, j, k, n
    636            temp = 273.0
     603          write (*, *) i, j, k, n
     604          temp = 273.0
    637605        endif
    638606        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
     
    643611  end do
    644612
    645   ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
    646   ! USING BOLTON'S (1980) FORMULA
    647   ! 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
    648616
    649617  do j=0,ny-1
    650618    do i=0,nxfield-1
    651         help=qvh2(i,j)
    652         temp=tt2(i,j,1,n)
    653         if (temp .eq. 0.0) then
    654            write (*, *) i, j, n
    655            temp = 273.0
    656         endif
    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