Changeset 8a65cb0 in flexpart.git for src/readwind.f90


Ignore:
Timestamp:
Mar 2, 2015, 3:11:55 PM (9 years ago)
Author:
Espen Sollum ATMOS <espen@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
1d207bb
Parents:
60403cd
Message:

Added code, makefile for dev branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/readwind.f90

    r4fbe7a5 r8a65cb0  
    2222subroutine readwind(indj,n,uuh,vvh,wwh)
    2323
    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: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
    34   !                                 ECMWF grib_api                      *
    35   !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    36   !                                 ECMWF grib_api                      *
    37   !                                                                     *
    38   !**********************************************************************
    39   !  Changes, Bernd C. Krueger, Feb. 2001:
    40   !   Variables tth and qvh (on eta coordinates) in common block
    41   !**********************************************************************
    42   !                                                                     *
    43   ! DESCRIPTION:                                                        *
    44   !                                                                     *
    45   ! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
    46   ! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
    47   !                                                                     *
    48   ! INPUT:                                                              *
    49   ! indj               indicates number of the wind field to be read in *
    50   ! n                  temporal index for meteorological fields (1 to 3)*
    51   !                                                                     *
    52   ! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
    53   !                                                                     *
    54   ! wfname             File name of data to be read in                  *
    55   ! nx,ny,nuvz,nwz     expected field dimensions                        *
    56   ! nlev_ec            number of vertical levels ecmwf model            *
    57   ! uu,vv,ww           wind fields                                      *
    58   ! tt,qv              temperature and specific humidity                *
    59   ! ps                 surface pressure                                 *
    60   !                                                                     *
    61   !**********************************************************************
    62 
    63   use GRIB_API
     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: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
     34!                                 ECMWF grib_api                      *
     35!             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
     36!                                 ECMWF grib_api                      *
     37!                                                                     *
     38!**********************************************************************
     39!  Changes, Bernd C. Krueger, Feb. 2001:
     40!   Variables tth and qvh (on eta coordinates) in common block
     41!**********************************************************************
     42!                                                                     *
     43! DESCRIPTION:                                                        *
     44!                                                                     *
     45! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
     46! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
     47!                                                                     *
     48! INPUT:                                                              *
     49! indj               indicates number of the wind field to be read in *
     50! n                  temporal index for meteorological fields (1 to 3)*
     51!                                                                     *
     52! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
     53!                                                                     *
     54! wfname             File name of data to be read in                  *
     55! nx,ny,nuvz,nwz     expected field dimensions                        *
     56! nlev_ec            number of vertical levels ecmwf model            *
     57! uu,vv,ww           wind fields                                      *
     58! tt,qv              temperature and specific humidity                *
     59! ps                 surface pressure                                 *
     60!                                                                     *
     61!**********************************************************************
     62
     63  use grib_api
    6464  use par_mod
    6565  use com_mod
     
    6767  implicit none
    6868
    69   !HSO  parameters for grib_api
     69!  include 'grib_api.h'
     70
     71!HSO  parameters for grib_api
    7072  integer :: ifile
    7173  integer :: iret
     
    7375  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
    7476  integer :: gotGrid
    75   !HSO  end
     77!HSO  end
    7678
    7779  real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     
    8082  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
    8183
    82   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    83 
    84   ! dimension of isec2 at least (22+n), where n is the number of parallels or
    85   ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
    86 
    87   ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
    88   ! coordinate parameters
     84! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
     85
     86! dimension of isec2 at least (22+n), where n is the number of parallels or
     87! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
     88
     89! dimension of zsec2 at least (10+nn), where nn is the number of vertical
     90! coordinate parameters
    8991
    9092  integer :: isec1(56),isec2(22+nxmax+nymax)
    9193  real(kind=4) :: zsec4(jpunp)
    92   real(kind=4) :: xaux,yaux
     94  real(kind=4) :: xaux,yaux,xaux0,yaux0
    9395  real(kind=8) :: xauxin,yauxin
    9496  real,parameter :: eps=1.e-4
     
    9698  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    9799
    98   logical :: hflswitch,strswitch
    99 
    100   !HSO  grib api error messages
     100  logical :: hflswitch,strswitch, readcloud
     101
     102!HSO  grib api error messages
    101103  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    102104  character(len=20) :: gribFunction = 'readwind'
     
    104106  hflswitch=.false.
    105107  strswitch=.false.
     108!hg
     109  readcloud=.false.
     110!hg end
    106111  levdiff2=nlev_ec-nwz+1
    107112  iumax=0
    108113  iwmax=0
    109114
    110   !
    111   ! OPENING OF DATA FILE (GRIB CODE)
    112   !
    113 5   call grib_open_file(ifile,path(3)(1:length(3)) &
    114          //trim(wfname(indj)),'r',iret)
     115!
     116! OPENING OF DATA FILE (GRIB CODE)
     117!
     1185 call grib_open_file(ifile,path(3)(1:length(3)) &
     119       //trim(wfname(indj)),'r',iret)
    115120  if (iret.ne.GRIB_SUCCESS) then
    116121    goto 888   ! ERROR DETECTED
    117122  endif
    118   !turn on support for multi fields messages */
    119   !call grib_multi_support_on
     123!turn on support for multi fields messages */
     124!call grib_multi_support_on
    120125
    121126  gotGrid=0
    122127  ifield=0
    123 10   ifield=ifield+1
    124   !
    125   ! GET NEXT FIELDS
    126   !
     12810 ifield=ifield+1
     129!
     130! GET NEXT FIELDS
     131!
    127132  call grib_new_from_file(ifile,igrib,iret)
    128133  if (iret.eq.GRIB_END_OF_FILE)  then
     
    132137  endif
    133138
    134   !first see if we read GRIB1 or GRIB2
     139!first see if we read GRIB1 or GRIB2
    135140  call grib_get_int(igrib,'editionNumber',gribVer,iret)
    136141  call grib_check(iret,gribFunction,gribErrorMsg)
     
    138143  if (gribVer.eq.1) then ! GRIB Edition 1
    139144
    140   !print*,'GRiB Edition 1'
    141   !read the grib2 identifiers
    142   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),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 
    147   !change code for etadot to code for omega
    148   if (isec1(6).eq.77) then
    149     isec1(6)=135
    150   endif
    151 
    152   conversion_factor=1.
     145!print*,'GRiB Edition 1'
     146!read the grib2 identifiers
     147    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     148    call grib_check(iret,gribFunction,gribErrorMsg)
     149    call grib_get_int(igrib,'level',isec1(8),iret)
     150    call grib_check(iret,gribFunction,gribErrorMsg)
     151
     152!change code for etadot to code for omega
     153    if (isec1(6).eq.77) then
     154      isec1(6)=135
     155    endif
     156
     157    conversion_factor=1.
    153158
    154159  else
    155160
    156   !print*,'GRiB Edition 2'
    157   !read the grib2 identifiers
    158   call grib_get_int(igrib,'discipline',discipl,iret)
    159   call grib_check(iret,gribFunction,gribErrorMsg)
    160   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    161   call grib_check(iret,gribFunction,gribErrorMsg)
    162   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    163   call grib_check(iret,gribFunction,gribErrorMsg)
    164   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    165   call grib_check(iret,gribFunction,gribErrorMsg)
    166   call grib_get_int(igrib,'level',valSurf,iret)
    167   call grib_check(iret,gribFunction,gribErrorMsg)
    168   call grib_get_int(igrib,'paramId',parId,iret)
    169   call grib_check(iret,gribFunction,gribErrorMsg)
    170 
    171   !print*,discipl,parCat,parNum,typSurf,valSurf
    172 
    173   !convert to grib1 identifiers
    174   isec1(6)=-1
    175   isec1(7)=-1
    176   isec1(8)=-1
    177   isec1(8)=valSurf     ! level
    178   conversion_factor=1.
    179   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    180     isec1(6)=130         ! indicatorOfParameter
    181   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
    182     isec1(6)=131         ! indicatorOfParameter
    183   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
    184     isec1(6)=132         ! indicatorOfParameter
    185   elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    186     isec1(6)=133         ! indicatorOfParameter
    187   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    188     isec1(6)=134         ! indicatorOfParameter
    189   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    190     isec1(6)=135         ! indicatorOfParameter
    191   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    192     isec1(6)=135         ! indicatorOfParameter
    193   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    194     isec1(6)=151         ! indicatorOfParameter
    195   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
    196     isec1(6)=165         ! indicatorOfParameter
    197   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
    198     isec1(6)=166         ! indicatorOfParameter
    199   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
    200     isec1(6)=167         ! indicatorOfParameter
    201   elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
    202     isec1(6)=168         ! indicatorOfParameter
    203   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    204     isec1(6)=141         ! indicatorOfParameter
    205     conversion_factor=1000.
    206   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    207     isec1(6)=164         ! indicatorOfParameter
    208   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    209     isec1(6)=142         ! indicatorOfParameter
    210   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    211     isec1(6)=143         ! indicatorOfParameter
    212     conversion_factor=1000.
    213   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    214     isec1(6)=146         ! indicatorOfParameter
    215   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    216     isec1(6)=176         ! indicatorOfParameter
    217   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    218     isec1(6)=180         ! indicatorOfParameter
    219   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    220     isec1(6)=181         ! indicatorOfParameter
    221   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    222     isec1(6)=129         ! indicatorOfParameter
    223   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    224     isec1(6)=160         ! indicatorOfParameter
    225   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    226        (typSurf.eq.1)) then ! LSM
    227     isec1(6)=172         ! indicatorOfParameter
    228   else
    229     print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    230          parCat,parNum,typSurf
    231   endif
    232   if(parId .ne. isec1(6) .and. parId .ne. 77) then
    233     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     161!print*,'GRiB Edition 2'
     162!read the grib2 identifiers
     163    call grib_get_int(igrib,'discipline',discipl,iret)
     164    call grib_check(iret,gribFunction,gribErrorMsg)
     165    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     166    call grib_check(iret,gribFunction,gribErrorMsg)
     167    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     168    call grib_check(iret,gribFunction,gribErrorMsg)
     169    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     170    call grib_check(iret,gribFunction,gribErrorMsg)
     171    call grib_get_int(igrib,'level',valSurf,iret)
     172    call grib_check(iret,gribFunction,gribErrorMsg)
     173    call grib_get_int(igrib,'paramId',parId,iret)
     174    call grib_check(iret,gribFunction,gribErrorMsg)
     175
     176!print*,discipl,parCat,parNum,typSurf,valSurf
     177
     178!convert to grib1 identifiers
     179    isec1(6)=-1
     180    isec1(7)=-1
     181    isec1(8)=-1
     182    isec1(8)=valSurf     ! level
     183    conversion_factor=1.
     184    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     185      isec1(6)=130         ! indicatorOfParameter
     186    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
     187      isec1(6)=131         ! indicatorOfParameter
     188    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
     189      isec1(6)=132         ! indicatorOfParameter
     190    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     191      isec1(6)=133         ! indicatorOfParameter
     192!hg
     193    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     194      isec1(6)=246         ! indicatorOfParameter
     195    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     196      isec1(6)=247         ! indicatorOfParameter
     197 !hg end
     198    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     199      isec1(6)=134         ! indicatorOfParameter
     200    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     201      isec1(6)=135         ! indicatorOfParameter
     202    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
     203      isec1(6)=135         ! indicatorOfParameter
     204    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     205      isec1(6)=151         ! indicatorOfParameter
     206    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
     207      isec1(6)=165         ! indicatorOfParameter
     208    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
     209      isec1(6)=166         ! indicatorOfParameter
     210    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
     211      isec1(6)=167         ! indicatorOfParameter
     212    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
     213      isec1(6)=168         ! indicatorOfParameter
     214    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     215      isec1(6)=141         ! indicatorOfParameter
     216      conversion_factor=1000.
     217    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     218      isec1(6)=164         ! indicatorOfParameter
     219    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     220      isec1(6)=142         ! indicatorOfParameter
     221    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     222      isec1(6)=143         ! indicatorOfParameter
     223      conversion_factor=1000.
     224    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     225      isec1(6)=146         ! indicatorOfParameter
     226    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     227      isec1(6)=176         ! indicatorOfParameter
     228    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     229      isec1(6)=180         ! indicatorOfParameter
     230    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     231      isec1(6)=181         ! indicatorOfParameter
     232    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     233      isec1(6)=129         ! indicatorOfParameter
     234    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     235      isec1(6)=160         ! indicatorOfParameter
     236    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
     237         (typSurf.eq.1)) then ! LSM
     238      isec1(6)=172         ! indicatorOfParameter
     239    else
     240      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
     241           parCat,parNum,typSurf
     242    endif
     243    if(parId .ne. isec1(6) .and. parId .ne. 77) then
     244      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
    234245!    stop
    235   endif
    236 
    237   endif
    238 
    239   !HSO  get the size and data of the values array
     246    endif
     247
     248  endif
     249
     250!HSO  get the size and data of the values array
    240251  if (isec1(6).ne.-1) then
    241252    call grib_get_real4_array(igrib,'values',zsec4,iret)
     
    243254  endif
    244255
    245   !HSO  get the required fields from section 2 in a gribex compatible manner
     256!HSO  get the required fields from section 2 in a gribex compatible manner
    246257  if (ifield.eq.1) then
    247   call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
    248   call grib_check(iret,gribFunction,gribErrorMsg)
    249   call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
    250   call grib_check(iret,gribFunction,gribErrorMsg)
    251   call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
    252   call grib_check(iret,gribFunction,gribErrorMsg)
    253   ! CHECK GRID SPECIFICATIONS
    254   if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    255   if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    256   if(isec2(12)/2-1.ne.nlev_ec) &
    257   stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
     258    call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
     259    call grib_check(iret,gribFunction,gribErrorMsg)
     260    call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
     261    call grib_check(iret,gribFunction,gribErrorMsg)
     262    call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
     263    call grib_check(iret,gribFunction,gribErrorMsg)
     264! CHECK GRID SPECIFICATIONS
     265    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
     266    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
     267    if(isec2(12)/2-1.ne.nlev_ec) &
     268         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
    258269  endif ! ifield
    259270
    260   !HSO  get the second part of the grid dimensions only from GRiB1 messages
     271!HSO  get the second part of the grid dimensions only from GRiB1 messages
    261272  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
    262273    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
     
    273284    if (xaux.gt.180.) xaux=xaux-360.0
    274285    if(abs(xaux-xlon0).gt.eps) &
    275     stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
     286         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    276287    if(abs(yaux-ylat0).gt.eps) &
    277     stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     288         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
    278289    gotGrid=1
    279290  endif ! gotGrid
     
    292303        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
    293304             qvh(i,j,nlev_ec-k+2,n) = 0.
    294   !        this is necessary because the gridded data may contain
    295   !        spurious negative values
     305!        this is necessary because the gridded data may contain
     306!        spurious negative values
    296307      endif
    297308      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
     
    336347      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
    337348           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
    338   !sec        strswitch=.true.
     349!sec        strswitch=.true.
    339350      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
    340351           zsec4(nxfield*(ny-j-1)+i+1)/ga
     
    345356      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
    346357      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
     358!hg READING CLOUD FIELDS ASWELL
     359      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content
     360        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     361        readclouds = .true.
     362        !write(*,*) 'found water!'
     363      endif
     364
     365      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
     366        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     367        !write(*,*) 'found ice!'
     368      endif
     369!hg end
    347370
    348371    end do
     
    351374  call grib_release(igrib)
    352375  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    353   !
    354   ! CLOSING OF INPUT DATA FILE
    355   !
    356 
    357 50   call grib_close_file(ifile)
    358 
    359   !error message if no fields found with correct first longitude in it
     376!
     377! CLOSING OF INPUT DATA FILE
     378!
     379
     38050 call grib_close_file(ifile)
     381
     382
     383!error message if no fields found with correct first longitude in it
    360384  if (gotGrid.eq.0) then
    361385    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
     
    373397  endif
    374398
    375   ! For global fields, assign the leftmost data column also to the rightmost
    376   ! data column; if required, shift whole grid by nxshift grid points
    377   !*************************************************************************
     399! For global fields, assign the leftmost data column also to the rightmost
     400! data column; if required, shift whole grid by nxshift grid points
     401!*************************************************************************
    378402
    379403  if (xglobal) then
     
    400424    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
    401425    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
     426!hg
     427    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
     428    call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
     429!hg end
     430
    402431  endif
    403432
     
    412441         wfname(indj)
    413442
    414   ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
    415   ! As ECMWF has increased the model resolution, such that now the first model
    416   ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
    417   ! (3rd model level in FLEXPART) for the profile method
    418   !***************************************************************************
     443! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
     444! As ECMWF has increased the model resolution, such that now the first model
     445! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
     446! (3rd model level in FLEXPART) for the profile method
     447!***************************************************************************
    419448
    420449    do i=0,nxmin1
     
    437466
    438467
    439   ! Assign 10 m wind to model level at eta=1.0 to have one additional model
    440   ! level at the ground
    441   ! Specific humidity is taken the same as at one level above
    442   ! Temperature is taken as 2 m temperature
    443   !**************************************************************************
    444 
    445      do i=0,nxmin1
    446         do j=0,nymin1
    447            uuh(i,j,1)=u10(i,j,1,n)
    448            vvh(i,j,1)=v10(i,j,1,n)
    449            qvh(i,j,1,n)=qvh(i,j,2,n)
    450            tth(i,j,1,n)=tt2(i,j,1,n)
    451         end do
    452      end do
     468! Assign 10 m wind to model level at eta=1.0 to have one additional model
     469! level at the ground
     470! Specific humidity is taken the same as at one level above
     471! Temperature is taken as 2 m temperature
     472!**************************************************************************
     473
     474  do i=0,nxmin1
     475    do j=0,nymin1
     476      uuh(i,j,1)=u10(i,j,1,n)
     477      vvh(i,j,1)=v10(i,j,1,n)
     478      qvh(i,j,1,n)=qvh(i,j,2,n)
     479      tth(i,j,1,n)=tt2(i,j,1,n)
     480    end do
     481  end do
    453482
    454483  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
     
    456485
    457486  return
    458 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     487
     488888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    459489  write(*,*) ' #### ',wfname(indj),'                    #### '
    460490  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
    461491  stop 'Execution terminated'
    462 999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    463   write(*,*) ' #### ',wfname(indj),'                    #### '
    464   write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
    465   stop 'Execution terminated'
    466492
    467493end subroutine readwind
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG