Ignore:
Timestamp:
Apr 22, 2015, 3:42:35 PM (9 years ago)
Author:
pesei
Message:

Wet dep quick fix and other small changes. Wet depo quick fix not final yet.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/petra/src/readwind_nests.f90

    r36 r37  
    3434  !                                                                            *
    3535  !*****************************************************************************
    36   !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
     36  !  CHANGES:
     37  !
     38  !  2/2001, Bernc C. Krueger:                                              *
    3739  !        Variables tthn and qvhn (on eta coordinates) in common block        *
    38   !  CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api    *
    39   !  CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api    *
    40   ! 2/2015, PS: add missing paramter iret in call to grib subr
     40  !  11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api            *
     41  !  03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api            *
     42  !  2/2015, PS: add missing paramter iret in call to grib subr
     43  !  3/2015, PS: add verbosity output
     44  !  3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
    4145  !
    4246  !*****************************************************************************
    … …  
    8488  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    8589  character(len=20) :: gribFunction = 'readwind_nests'
     90 
     91  call verboutput(verbosity,1,'      entered readwind_nests')
    8692
    8793  do l=1,numbnests
    … …  
    100106  !
    101107
    102 5   call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
     1085 continue
     109  call verboutput(verbosity,1,'        call grib_open_file '//&
     110    path(numpath+2*(l-1)+1)(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)))
     111  call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
    103112         (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r')
    104113  if (iret.ne.GRIB_SUCCESS) then
    … …  
    127136  if (gribVer.eq.1) then ! GRIB Edition 1
    128137
    129   !print*,'GRiB Edition 1'
    130   !read the grib2 identifiers
    131   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    132   call grib_check(iret,gribFunction,gribErrorMsg)
    133   call grib_get_int(igrib,'level',isec1(8),iret)
    134   call grib_check(iret,gribFunction,gribErrorMsg)
    135 
    136   !change code for etadot to code for omega
    137   if (isec1(6).eq.77) then
    138     isec1(6)=135
    139   endif
    140 
    141   conversion_factor=1.
    142 
     138    !print*,'GRiB Edition 1'
     139   
     140    ! read the grib2 identifiers
     141   
     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.
    143153
    144154  else
    145155
    146   !print*,'GRiB Edition 2'
    147   !read the grib2 identifiers
    148   call grib_get_int(igrib,'discipline',discipl,iret)
    149   call grib_check(iret,gribFunction,gribErrorMsg)
    150   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    151   call grib_check(iret,gribFunction,gribErrorMsg)
    152   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    153   call grib_check(iret,gribFunction,gribErrorMsg)
    154   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    155   call grib_check(iret,gribFunction,gribErrorMsg)
    156   call grib_get_int(igrib,'level',valSurf,iret)
    157   call grib_check(iret,gribFunction,gribErrorMsg)
    158   call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
    159   call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
    160 
    161   !print*,discipl,parCat,parNum,typSurf,valSurf
    162 
    163   !convert to grib1 identifiers
    164   isec1(6)=-1
    165   isec1(7)=-1
    166   isec1(8)=-1
    167   isec1(8)=valSurf     ! level
    168    conversion_factor=1.
    169   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    170     isec1(6)=130         ! indicatorOfParameter
    171   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
    172     isec1(6)=131         ! indicatorOfParameter
    173   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
    174     isec1(6)=132         ! indicatorOfParameter
    175   elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    176     isec1(6)=133         ! indicatorOfParameter
    177   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    178     isec1(6)=134         ! indicatorOfParameter
    179   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
    180     isec1(6)=135         ! indicatorOfParameter
    181   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
    182     isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
    183   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    184     isec1(6)=151         ! indicatorOfParameter
    185   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
    186     isec1(6)=165         ! indicatorOfParameter
    187   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
    188     isec1(6)=166         ! indicatorOfParameter
    189   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
    190     isec1(6)=167         ! indicatorOfParameter
    191   elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
    192     isec1(6)=168         ! indicatorOfParameter
    193   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    194     isec1(6)=141         ! indicatorOfParameter
    195     conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
    196   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
    197     isec1(6)=164         ! indicatorOfParameter
    198   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
    199     isec1(6)=142         ! indicatorOfParameter
    200   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    201     isec1(6)=143         ! indicatorOfParameter
    202     conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
    203   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    204     isec1(6)=146         ! indicatorOfParameter
    205   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    206     isec1(6)=176         ! indicatorOfParameter
    207   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
    208     isec1(6)=180         ! indicatorOfParameter
    209   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
    210     isec1(6)=181         ! indicatorOfParameter
    211   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    212     isec1(6)=129         ! indicatorOfParameter
    213    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
    214     isec1(6)=160         ! indicatorOfParameter
    215   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    216        (typSurf.eq.1)) then ! LSM
    217     isec1(6)=172         ! indicatorOfParameter
    218   else
    219     print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    220          parCat,parNum,typSurf
    221   endif
    222   if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
    223     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
    224 !    stop
    225   endif
     156    !print*,'GRiB Edition 2'
     157
     158    ! read the grib2 identifiers
     159
     160    call grib_get_int(igrib,'discipline',discipl,iret)
     161    call grib_check(iret,gribFunction,gribErrorMsg)
     162    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     163    call grib_check(iret,gribFunction,gribErrorMsg)
     164    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     165    call grib_check(iret,gribFunction,gribErrorMsg)
     166    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     167    call grib_check(iret,gribFunction,gribErrorMsg)
     168    call grib_get_int(igrib,'level',valSurf,iret)
     169    call grib_check(iret,gribFunction,gribErrorMsg)
     170    call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consistent with new readwind.f90
     171    call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consistent with new readwind.f90
     172
     173    !print*,discipl,parCat,parNum,typSurf,valSurf
     174
     175    !convert to grib1 identifiers
     176    isec1(6)=-1
     177    isec1(7)=-1
     178    isec1(8)=-1
     179    isec1(8)=valSurf     ! level
     180     conversion_factor=1.
     181    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     182      isec1(6)=130         ! indicatorOfParameter
     183    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
     184      isec1(6)=131         ! indicatorOfParameter
     185    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
     186      isec1(6)=132         ! indicatorOfParameter
     187    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     188      isec1(6)=133         ! indicatorOfParameter
     189    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     190      isec1(6)=134         ! indicatorOfParameter
     191    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta
     192      isec1(6)=135         ! indicatorOfParameter
     193    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
     194      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
     195    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
     196      isec1(6)=135         ! indicatorOfParameter !added by mc to make it consistent with new readwind.f90
     197    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     198      isec1(6)=151         ! indicatorOfParameter
     199    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
     200      isec1(6)=165         ! indicatorOfParameter
     201    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
     202      isec1(6)=166         ! indicatorOfParameter
     203    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
     204      isec1(6)=167         ! indicatorOfParameter
     205    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
     206      isec1(6)=168         ! indicatorOfParameter
     207    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     208      isec1(6)=141         ! indicatorOfParameter
     209      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
     210    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new readwind.f90
     211      isec1(6)=164         ! indicatorOfParameter
     212    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new readwind.f90
     213      isec1(6)=142         ! indicatorOfParameter
     214    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     215      isec1(6)=143         ! indicatorOfParameter
     216      conversion_factor=1000. !added by mc to make it consistent with new readwind.f90
     217    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     218      isec1(6)=146         ! indicatorOfParameter
     219    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     220      isec1(6)=176         ! indicatorOfParameter
     221    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new readwind.f90
     222      isec1(6)=180         ! indicatorOfParameter
     223    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new readwind.f90
     224      isec1(6)=181         ! indicatorOfParameter
     225    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     226      isec1(6)=129         ! indicatorOfParameter
     227     elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new readwind.f90
     228      isec1(6)=160         ! indicatorOfParameter
     229    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
     230         (typSurf.eq.1)) then ! LSM
     231      isec1(6)=172         ! indicatorOfParameter
     232    else
     233      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
     234           parCat,parNum,typSurf
     235    endif
     236    if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new readwind.f90
     237      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
     238  !    stop
     239    endif
    226240
    227241  endif
    … …  
    246260  ! CHECK GRID SPECIFICATIONS
    247261  if(isec2(2).ne.nxn(l)) stop &
    248   'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
     262    'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
    249263  if(isec2(3).ne.nyn(l)) stop &
    250   'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
    251   if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
    252   IZATION NOT CONSISTENT FOR A NESTING LEVEL'
     264    'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
     265  if(isec2(12)/2-1.ne.nlev_ec) stop &
     266    'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT FOR A NESTING LEVEL'
    253267  endif ! ifield
    254268
    255269  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    256  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
     270 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consistent with new readwind.f90
    257271    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    258272         xauxin,iret)
    … …  
    294308             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    295309        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
    296              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
     310             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90!
    297311        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
    298312             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    … …  
    312326        endif
    313327        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
    314           convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
     328          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consistent with new readwind.f90
    315329          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
    316330        endif
    … …  
    421435  end do
    422436
     437  call verboutput(verbosity,1,'      leaving readwind_nests')
     438
    423439  return
     440
    424441888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    425442  write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
    … …  
    431448  write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
    432449  write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
     450  call verboutput(verbosity,1,'      leaving readwind')
    433451
    434452end subroutine readwind_nests
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG