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.f90

    r36 r37  
    4040  !   Variables tth and qvh (on eta coordinates) in common block
    4141  ! 2/2015, PS: add missing paramter iret in call to grib subr
     42  ! 3/2015, PS: add some verbosity messages
     43  ! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
    4244  !
    4345  !**********************************************************************
    … …  
    113115  ! OPENING OF DATA FILE (GRIB CODE)
    114116  !
    115 5   call grib_open_file(ifile,path(3)(1:length(3)) &
     117 
     1185 continue
     119  call verboutput(verbosity,1,'      call grib_open_file '//path(3)(1:length(3))&
     120         //trim(wfname(indj)))
     121  call grib_open_file(ifile,path(3)(1:length(3)) &
    116122         //trim(wfname(indj)),'r',iret)
    117123  if (iret.ne.GRIB_SUCCESS) then
    … …  
    127133  ! GET NEXT FIELDS
    128134  !
     135  call verboutput(verbosity,2,'      call grib_new_from_file')
    129136  call grib_new_from_file(ifile,igrib,iret)
    130137  if (iret.eq.GRIB_END_OF_FILE)  then
    … …  
    135142
    136143  !first see if we read GRIB1 or GRIB2
     144  call verboutput(verbosity,2,'      call grib_get_int')
    137145  call grib_get_int(igrib,'editionNumber',gribVer,iret)
     146  call verboutput(verbosity,2,'      call grib_check')
    138147  call grib_check(iret,gribFunction,gribErrorMsg)
    139148
    140149  if (gribVer.eq.1) then ! GRIB Edition 1
    141150
    142   !print*,'GRiB Edition 1'
    143   !read the grib2 identifiers
    144   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    145   call grib_check(iret,gribFunction,gribErrorMsg)
    146   call grib_get_int(igrib,'level',isec1(8),iret)
    147   call grib_check(iret,gribFunction,gribErrorMsg)
    148 
    149   !change code for etadot to code for omega
    150   if (isec1(6).eq.77) then
    151     isec1(6)=135
    152   endif
    153 
    154   conversion_factor=1.
    155 
    156   else
    157 
    158   !print*,'GRiB Edition 2'
    159   !read the grib2 identifiers
    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)
    171   call grib_check(iret,gribFunction,gribErrorMsg)
    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 eta dot
    192     isec1(6)=135         ! indicatorOfParameter
    193   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    194     isec1(6)=135         ! indicatorOfParameter
    195   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    196     isec1(6)=151         ! indicatorOfParameter
    197   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
    198     isec1(6)=165         ! indicatorOfParameter
    199   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
    200     isec1(6)=166         ! indicatorOfParameter
    201   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
    202     isec1(6)=167         ! indicatorOfParameter
    203   elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
    204     isec1(6)=168         ! indicatorOfParameter
    205   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    206     isec1(6)=141         ! indicatorOfParameter
    207     conversion_factor=1000.
    208   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    209     isec1(6)=164         ! indicatorOfParameter
    210   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    211     isec1(6)=142         ! indicatorOfParameter
    212   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    213     isec1(6)=143         ! indicatorOfParameter
    214     conversion_factor=1000.
    215   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    216     isec1(6)=146         ! indicatorOfParameter
    217   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    218     isec1(6)=176         ! indicatorOfParameter
    219   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    220     isec1(6)=180         ! indicatorOfParameter
    221   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    222     isec1(6)=181         ! indicatorOfParameter
    223   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    224     isec1(6)=129         ! indicatorOfParameter
    225   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    226     isec1(6)=160         ! indicatorOfParameter
    227   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    228        (typSurf.eq.1)) then ! LSM
    229     isec1(6)=172         ! indicatorOfParameter
    230   else
    231     print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    232          parCat,parNum,typSurf
    233   endif
    234   if(parId .ne. isec1(6) .and. parId .ne. 77) then
    235     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
    236 !    stop
    237   endif
    238 
    239   endif
     151    !print*,'GRiB Edition 1'
     152
     153    ! read the grib2 identifiers
     154
     155    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     156    call grib_check(iret,gribFunction,gribErrorMsg)
     157    call grib_get_int(igrib,'level',isec1(8),iret)
     158    call grib_check(iret,gribFunction,gribErrorMsg)
     159
     160    !change code for etadot to code for omega
     161    if (isec1(6).eq.77) then
     162      isec1(6)=135
     163    endif
     164
     165    conversion_factor=1.
     166
     167  else ! GRIB 2
     168
     169    !print*,'GRiB Edition 2'
     170
     171    ! read the grib2 identifiers
     172
     173    call grib_get_int(igrib,'discipline',discipl,iret)
     174    call grib_check(iret,gribFunction,gribErrorMsg)
     175    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     176    call grib_check(iret,gribFunction,gribErrorMsg)
     177    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     178    call grib_check(iret,gribFunction,gribErrorMsg)
     179    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     180    call grib_check(iret,gribFunction,gribErrorMsg)
     181    call grib_get_int(igrib,'level',valSurf,iret)
     182    call grib_check(iret,gribFunction,gribErrorMsg)
     183    call grib_get_int(igrib,'paramId',parId,iret)
     184    call grib_check(iret,gribFunction,gribErrorMsg)
     185
     186    !print*,discipl,parCat,parNum,typSurf,valSurf
     187
     188    ! convert to grib1 identifiers
     189
     190    isec1(6)=-1
     191    isec1(7)=-1
     192    isec1(8)=-1
     193    isec1(8)=valSurf     ! level
     194    conversion_factor=1.
     195    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     196      isec1(6)=130         ! indicatorOfParameter
     197    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
     198      isec1(6)=131         ! indicatorOfParameter
     199    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
     200      isec1(6)=132         ! indicatorOfParameter
     201    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     202      isec1(6)=133         ! indicatorOfParameter
     203    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     204      isec1(6)=134         ! indicatorOfParameter
     205    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta
     206      isec1(6)=135         ! indicatorOfParameter
     207    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
     208      isec1(6)=135         ! indicatorOfParameter
     209    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
     210      isec1(6)=135         ! indicatorOfParameter
     211    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     212      isec1(6)=151         ! indicatorOfParameter
     213    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
     214      isec1(6)=165         ! indicatorOfParameter
     215    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
     216      isec1(6)=166         ! indicatorOfParameter
     217    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
     218      isec1(6)=167         ! indicatorOfParameter
     219    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
     220      isec1(6)=168         ! indicatorOfParameter
     221    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     222      isec1(6)=141         ! indicatorOfParameter
     223      conversion_factor=1000.
     224    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     225      isec1(6)=164         ! indicatorOfParameter
     226    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     227      isec1(6)=142         ! indicatorOfParameter
     228    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     229      isec1(6)=143         ! indicatorOfParameter
     230      conversion_factor=1000.
     231    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     232      isec1(6)=146         ! indicatorOfParameter
     233    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     234      isec1(6)=176         ! indicatorOfParameter
     235    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     236      isec1(6)=180         ! indicatorOfParameter
     237    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     238      isec1(6)=181         ! indicatorOfParameter
     239    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     240      isec1(6)=129         ! indicatorOfParameter
     241    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     242      isec1(6)=160         ! indicatorOfParameter
     243    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
     244         (typSurf.eq.1)) then ! LSM
     245      isec1(6)=172         ! indicatorOfParameter
     246    else
     247      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
     248           parCat,parNum,typSurf
     249    endif
     250    if(parId .ne. isec1(6) .and. parId .ne. 77) then
     251      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     252  !    stop
     253    endif
     254
     255  endif ! end GRIB 1 / GRIB 2 cases
    240256
    241257  !HSO  get the size and data of the values array
    … …  
    351367  end do
    352368
     369  call verboutput(verbosity,2,'      call grib_release')
    353370  call grib_release(igrib)
    354371  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    … …  
    380397
    381398  if (xglobal) then
     399    call verboutput(verbosity,2,'      call shift_field_0')
    382400    call shift_field_0(ewss,nxfield,ny)
    383401    call shift_field_0(nsss,nxfield,ny)
    … …  
    429447        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
    430448        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
     449        call verboutput(verbosity,2,'      call pbl_profile')
    431450        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
    432451             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
    … …  
    457476  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
    458477
     478  call verboutput(verbosity,1,'      leaving readwind')
    459479  return
     480
    460481888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    461482  write(*,*) ' #### ',wfname(indj),'                    #### '
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG