Changeset 467460a in flexpart.git


Ignore:
Timestamp:
Mar 13, 2021, 10:00:37 AM (3 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
GFS_025, dev
Children:
9ca6e38
Parents:
759df5f
Message:

First commit of files from Hachinger

Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/calcmatrix.f90

    r92fab65 r467460a  
    7171    k = kuvz-1
    7272    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    73     pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
    74     phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
     73      pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
     74      phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
    7575    else
    7676      phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))

  • src/ew.f90

    r92fab65 r467460a  
    1414
    1515  ew=0.
    16   if(x.le.0.) stop 'sorry: t not in [k]'
     16  if(x.le.0.) then
     17    write(*,*) 'in ew.f90: x=', x
     18    stop 'sorry: t not in [k]'
     19  end if
     20
    1721  y=373.16/x
    1822  a=-7.90298*(y-1.)

  • src/gridcheck_gfs.f90

    • Property mode changed from 100644 to 100755
    ra803521 r467460a  
    7575  real :: sizesouth,sizenorth,xauxa,pint
    7676  real :: akm_usort(nwzmax)
    77   real,parameter :: eps=0.0001
     77  real,parameter :: eps=spacing(2.0_4*360.0_4)
    7878
    7979  ! NCEP GFS
    8080  real :: pres(nwzmax), help
    8181
    82   integer :: i179,i180,i181
     82  integer :: i180
    8383
    8484  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    224224  nxfield=isec2(2)
    225225  ny=isec2(3)
    226   if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
    227     xaux1=-179.0                             ! 359 DEG EAST ->
    228     xaux2=-179.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
     226  if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO
     227    xaux1=-180.0                             ! 359 DEG EAST ->
     228    xaux2=-180.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
    229229  endif                                      ! TO 180 DEG EAST
    230230  if (xaux1.gt.180) xaux1=xaux1-360.0
    … …  
    305305
    306306
    307   i179=nint(179./dx)
    308   if (dx.lt.0.7) then
    309     i180=nint(180./dx)+1    ! 0.5 deg data
    310   else
    311     i180=nint(179./dx)+1    ! 1 deg data
    312   endif
    313   i181=i180+1
     307  i180=nint(180./dx)    ! 0.5 deg data
    314308
    315309
    … …  
    321315      do ix=0,nxfield-1
    322316        help=zsec4(nxfield*(ny-jy-1)+ix+1)
    323         if(ix.le.i180) then
    324           oro(i179+ix,jy)=help
    325           excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     317        if(ix.lt.i180) then
     318          oro(i180+ix,jy)=help
     319          excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    326320        else
    327           oro(ix-i181,jy)=help
    328           excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     321          oro(ix-i180,jy)=help
     322          excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    329323        endif
    330324      end do
    … …  
    339333      do ix=0,nxfield-1
    340334        help=zsec4(nxfield*(ny-jy-1)+ix+1)
    341         if(ix.le.i180) then
    342           lsm(i179+ix,jy)=help
     335        if(ix.lt.i180) then
     336          lsm(i180+ix,jy)=help
    343337        else
    344           lsm(ix-i181,jy)=help
     338          lsm(ix-i180,jy)=help
    345339        endif
    346340      end do
    … …  
    413407    write(*,*)
    414408    write(*,*)
    415     write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
     409  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
    416410         nuvz,nwz
    417411    write(*,*)

  • src/par_mod.f90

    r759df5f r467460a  
    146146
    147147! GFS
    148    integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
     148  !   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
     149  ! GFS 0.25
     150  integer,parameter :: nxmax=1441,nymax=721,nuvzmax=138,nwzmax=138,nzmax=138
    149151   integer :: nxshift=0 ! shift not fixed for the executable
    150152
    … …  
    212214  !**************************************************
    213215
    214   integer,parameter :: maxpart=100000
     216  integer,parameter :: maxpart=10000000
    215217  integer,parameter :: maxspec=6
    216218
    … …  
    256258  !*********************************
    257259
    258   integer,parameter :: maxrand=1000000
     260  integer,parameter :: maxrand=10000000
    259261
    260262  ! maxrand                 number of random numbers used

  • src/readreleases.f90

    • Property mode changed from 100644 to 100755
    r92fab65 r467460a  
    546546      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
    547547        write(*,*) 'FLEXPART MODEL ERROR'
    548         write(*,*) 'Release starts before simulation begins or ends'
     548        write(*,*) 'Release starts before simulation begins or ends (1)'
    549549        write(*,*) 'after simulation stops.'
    550         write(*,*) 'Make files COMMAND and RELEASES consistent.'
     550        write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
    551551        stop
    552552      endif
    … …  
    561561      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
    562562        write(*,*) 'FLEXPART MODEL ERROR'
    563         write(*,*) 'Release starts before simulation begins or ends'
     563        write(*,*) 'Release starts before simulation begins or ends (2)'
    564564        write(*,*) 'after simulation stops.'
    565         write(*,*) 'Make files COMMAND and RELEASES consistent.'
     565        write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
    566566        stop
    567567      endif

  • src/readwind_gfs.f90

    • Property mode changed from 100644 to 100755
    ra803521 r467460a  
    7272  real :: qvh2(0:nxmax-1,0:nymax-1)
    7373
    74   integer :: i179,i180,i181
     74  integer :: i180
    7575
    7676  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    7878
    7979  integer :: isec1(8),isec2(3)
     80  real           :: xsec18
    8081  real(kind=4) :: zsec4(jpunp)
    8182  real(kind=4) :: xaux,yaux,xaux0,yaux0
     83  real,parameter :: eps=spacing(2.0_4*360.0_4)
    8284  real(kind=8) :: xauxin,yauxin
    83   real,parameter :: eps=1.e-4
    8485  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
    8586  real :: plev1,hlev1,ff10m,fflev1
    … …  
    137138  !read the grib1 identifiers
    138139  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    139 !  call grib_check(iret,gribFunction,gribErrorMsg)
     140  call grib_check(iret,gribFunction,gribErrorMsg)
    140141  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
    141 !  call grib_check(iret,gribFunction,gribErrorMsg)
     142  call grib_check(iret,gribFunction,gribErrorMsg)
    142143  call grib_get_int(igrib,'level',isec1(8),iret)
    143 !  call grib_check(iret,gribFunction,gribErrorMsg)
     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))
    144148
    145149  else ! GRIB Edition 2
    … …  
    149153
    150154  call grib_get_int(igrib,'discipline',discipl,iret)
    151 !  call grib_check(iret,gribFunction,gribErrorMsg)
     155  call grib_check(iret,gribFunction,gribErrorMsg)
    152156  call grib_get_int(igrib,'parameterCategory',parCat,iret)
    153 !  call grib_check(iret,gribFunction,gribErrorMsg)
     157  call grib_check(iret,gribFunction,gribErrorMsg)
    154158  call grib_get_int(igrib,'parameterNumber',parNum,iret)
    155 !  call grib_check(iret,gribFunction,gribErrorMsg)
     159  call grib_check(iret,gribFunction,gribErrorMsg)
    156160  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    157 !  call grib_check(iret,gribFunction,gribErrorMsg)
     161  call grib_check(iret,gribFunction,gribErrorMsg)
    158162  call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
    159163       valSurf,iret)
    160 !  call grib_check(iret,gribFunction,gribErrorMsg)
     164  call grib_check(iret,gribFunction,gribErrorMsg)
    161165 
    162166!  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
    163167  !convert to grib1 identifiers
    164   isec1(6)=-1
    165   isec1(7)=-1
    166   isec1(8)=-1
     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
    167172  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
    168173    isec1(6)=11          ! indicatorOfParameter
    169174    isec1(7)=100         ! indicatorOfTypeOfLevel
    170     isec1(8)=valSurf/100 ! level, convert to hPa
     175    xsec18=valSurf/100.0 ! level, convert to hPa
    171176  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
    172177    isec1(6)=33          ! indicatorOfParameter
    173178    isec1(7)=100         ! indicatorOfTypeOfLevel
    174     isec1(8)=valSurf/100 ! level, convert to hPa
     179    xsec18=valSurf/100.0 ! level, convert to hPa
    175180  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
    176181    isec1(6)=34          ! indicatorOfParameter
    177182    isec1(7)=100         ! indicatorOfTypeOfLevel
    178     isec1(8)=valSurf/100 ! level, convert to hPa
     183    xsec18=valSurf/100.0 ! level, convert to hPa
    179184  elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
    180185    isec1(6)=39          ! indicatorOfParameter
    181186    isec1(7)=100         ! indicatorOfTypeOfLevel
    182     isec1(8)=valSurf/100 ! level, convert to hPa
     187    xsec18=valSurf/100.0 ! level, convert to hPa
    183188  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
    184189    isec1(6)=52          ! indicatorOfParameter
    185190    isec1(7)=100         ! indicatorOfTypeOfLevel
    186     isec1(8)=valSurf/100 ! level, convert to hPa
     191    xsec18=valSurf/100.0 ! level, convert to hPa
    187192  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
    188193    isec1(6)=52          ! indicatorOfParameter
    189194    isec1(7)=105         ! indicatorOfTypeOfLevel
    190     isec1(8)=2
     195    xsec18=real(2)
    191196  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
    192197    isec1(6)=11          ! indicatorOfParameter
    193198    isec1(7)=105         ! indicatorOfTypeOfLevel
    194     isec1(8)=2
     199    xsec18=real(2)
    195200  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
    196201    isec1(6)=33          ! indicatorOfParameter
    197202    isec1(7)=105         ! indicatorOfTypeOfLevel
    198     isec1(8)=10
     203    xsec18=real(10)
    199204  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
    200205    isec1(6)=34          ! indicatorOfParameter
    201206    isec1(7)=105         ! indicatorOfTypeOfLevel
    202     isec1(8)=10
     207    xsec18=real(10)
    203208  elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
    204209    isec1(6)=153         ! indicatorOfParameter
    205210    isec1(7)=100         ! indicatorOfTypeOfLevel
    206     isec1(8)=valSurf/100 ! level, convert to hPa
     211    xsec18=valSurf/100.0 ! level, convert to hPa
    207212  elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
    208213    isec1(6)=2           ! indicatorOfParameter
    209214    isec1(7)=102         ! indicatorOfTypeOfLevel
    210     isec1(8)=0
     215    xsec18=real(0)
    211216  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
    212217    isec1(6)=1           ! indicatorOfParameter
    213218    isec1(7)=1           ! indicatorOfTypeOfLevel
    214     isec1(8)=0
     219    xsec18=real(0)
    215220  elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
    216221    isec1(6)=66          ! indicatorOfParameter
    217222    isec1(7)=1           ! indicatorOfTypeOfLevel
    218     isec1(8)=0
     223    xsec18=real(0)
    219224  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
    220225    isec1(6)=11          ! indicatorOfParameter
    221226    isec1(7)=107         ! indicatorOfTypeOfLevel
    222     isec1(8)=0.995       ! lowest sigma level
     227    xsec18=0.995         ! lowest sigma level
    223228  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
    224229    isec1(6)=33          ! indicatorOfParameter
    225230    isec1(7)=107         ! indicatorOfTypeOfLevel
    226     isec1(8)=0.995       ! lowest sigma level
     231    xsec18=0.995         ! lowest sigma level
    227232  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
    228233    isec1(6)=34          ! indicatorOfParameter
    229234    isec1(7)=107         ! indicatorOfTypeOfLevel
    230     isec1(8)=0.995       ! lowest sigma level
     235    xsec18=0.995         ! lowest sigma level
    231236  elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
    232237    isec1(6)=7           ! indicatorOfParameter
    233238    isec1(7)=1           ! indicatorOfTypeOfLevel
    234     isec1(8)=0
     239    xsec18=real(0)
    235240  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
    236241       .and.(discipl.eq.2)) then ! LSM
    237242    isec1(6)=81          ! indicatorOfParameter
    238243    isec1(7)=1           ! indicatorOfTypeOfLevel
    239     isec1(8)=0
     244    xsec18=real(0)
    240245  elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
    241246    isec1(6)=221         ! indicatorOfParameter
    242247    isec1(7)=1           ! indicatorOfTypeOfLevel
    243     isec1(8)=0
     248    xsec18=real(0)
    244249  elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
    245250    isec1(6)=62          ! indicatorOfParameter
    246251    isec1(7)=1           ! indicatorOfTypeOfLevel
    247     isec1(8)=0
     252    xsec18=real(0)
    248253  elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
    249254    isec1(6)=63          ! indicatorOfParameter
    250255    isec1(7)=1           ! indicatorOfTypeOfLevel
    251     isec1(8)=0
     256    xsec18=real(0)
    252257  endif
    253258
    … …  
    257262  !  get the size and data of the values array
    258263    call grib_get_real4_array(igrib,'values',zsec4,iret)
    259 !    call grib_check(iret,gribFunction,gribErrorMsg)
     264    call grib_check(iret,gribFunction,gribErrorMsg)
    260265  endif
    261266
    … …  
    266271  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    267272       isec2(2),iret)
    268 !  call grib_check(iret,gribFunction,gribErrorMsg)
     273  call grib_check(iret,gribFunction,gribErrorMsg)
    269274  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
    270275       isec2(3),iret)
    271 !  call grib_check(iret,gribFunction,gribErrorMsg)
     276  call grib_check(iret,gribFunction,gribErrorMsg)
    272277  call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    273278       xauxin,iret)
    274 !  call grib_check(iret,gribFunction,gribErrorMsg)
     279  call grib_check(iret,gribFunction,gribErrorMsg)
    275280  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
    276281       yauxin,iret)
    277 !  call grib_check(iret,gribFunction,gribErrorMsg)
     282  call grib_check(iret,gribFunction,gribErrorMsg)
    278283  xaux=xauxin+real(nxshift)*dx
    279284  yaux=yauxin
    … …  
    283288    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    284289    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    285     if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
     290    if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
    286291    xaux0=xlon0
    287292    yaux0=ylat0
    … …  
    290295    if(xaux0.lt.0.) xaux0=xaux0+360.
    291296    if(yaux0.lt.0.) yaux0=yaux0+360.
    292     if(abs(xaux-xaux0).gt.eps) &
     297    if(abs(xaux-xaux0).gt.eps) then
     298         write (*, *) xaux, xaux0
    293299         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    294     if(abs(yaux-yaux0).gt.eps) &
     300    endif
     301    if(abs(yaux-yaux0).gt.eps) then
     302         write (*, *) yaux, yaux0
    295303         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     304    end if
    296305  endif
    297306  !HSO end of edits
    298307
    299   i179=nint(179./dx)
    300   if (dx.lt.0.7) then
    301     i180=nint(180./dx)+1    ! 0.5 deg data
    302   else
    303     i180=nint(179./dx)+1    ! 1 deg data
    304   endif
    305   i181=i180+1
     308  i180=nint(180./dx)
    306309
    307310  if (isec1(6).ne.-1) then
     311
     312! write (*, *) 'nxfield: ', nxfield, i180
    308313
    309314  do j=0,nymin1
    310315    do i=0,nxfield-1
    311316      if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
    312   ! TEMPERATURE
    313          if((i.eq.0).and.(j.eq.0)) then
    314             do ii=1,nuvz
    315               if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
    316             end do
    317         endif
    318         help=zsec4(nxfield*(ny-j-1)+i+1)
    319         if(i.le.i180) then
    320           tth(i179+i,j,numpt,n)=help
    321         else
    322           tth(i-i181,j,numpt,n)=help
     317! 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
    323343        endif
    324344      endif
    … …  
    327347         if((i.eq.0).and.(j.eq.0)) then
    328348            do ii=1,nuvz
    329               if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
     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
    330353            end do
    331354        endif
    332355        help=zsec4(nxfield*(ny-j-1)+i+1)
    333         if(i.le.i180) then
    334           uuh(i179+i,j,numpu)=help
    335         else
    336           uuh(i-i181,j,numpu)=help
     356        if(i.lt.i180) then
     357          uuh(i180+i,j,numpu)=help
     358        else
     359          uuh(i-i180,j,numpu)=help
    337360        endif
    338361      endif
    … …  
    341364         if((i.eq.0).and.(j.eq.0)) then
    342365            do ii=1,nuvz
    343               if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
     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
    344370            end do
    345371        endif
    346372        help=zsec4(nxfield*(ny-j-1)+i+1)
    347         if(i.le.i180) then
    348           vvh(i179+i,j,numpv)=help
    349         else
    350           vvh(i-i181,j,numpv)=help
     373        if(i.lt.i180) then
     374          vvh(i180+i,j,numpv)=help
     375        else
     376          vvh(i-i180,j,numpv)=help
    351377        endif
    352378      endif
    … …  
    355381         if((i.eq.0).and.(j.eq.0)) then
    356382            do ii=1,nuvz
    357               if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
     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
    358387            end do
    359388        endif
    360389        help=zsec4(nxfield*(ny-j-1)+i+1)
    361         if(i.le.i180) then
    362           qvh(i179+i,j,numprh,n)=help
    363         else
    364           qvh(i-i181,j,numprh,n)=help
     390        if(i.lt.i180) then
     391          qvh(i180+i,j,numprh,n)=help
     392        else
     393          qvh(i-i180,j,numprh,n)=help
    365394        endif
    366395      endif
    … …  
    368397  ! SURFACE PRESSURE
    369398        help=zsec4(nxfield*(ny-j-1)+i+1)
    370         if(i.le.i180) then
    371           ps(i179+i,j,1,n)=help
    372         else
    373           ps(i-i181,j,1,n)=help
     399        if(i.lt.i180) then
     400          ps(i180+i,j,1,n)=help
     401        else
     402          ps(i-i180,j,1,n)=help
    374403        endif
    375404      endif
    … …  
    378407         if((i.eq.0).and.(j.eq.0)) then
    379408            do ii=1,nuvz
    380               if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
     409              if ((xsec18*100.0).eq.akz(ii)) numpw=ii
    381410            end do
    382411        endif
    383412        help=zsec4(nxfield*(ny-j-1)+i+1)
    384         if(i.le.i180) then
    385           wwh(i179+i,j,numpw)=help
    386         else
    387           wwh(i-i181,j,numpw)=help
     413        if(i.lt.i180) then
     414          wwh(i180+i,j,numpw)=help
     415        else
     416          wwh(i-i180,j,numpw)=help
    388417        endif
    389418      endif
    … …  
    391420  ! SNOW DEPTH
    392421        help=zsec4(nxfield*(ny-j-1)+i+1)
    393         if(i.le.i180) then
    394           sd(i179+i,j,1,n)=help
    395         else
    396           sd(i-i181,j,1,n)=help
     422        if(i.lt.i180) then
     423          sd(i180+i,j,1,n)=help
     424        else
     425          sd(i-i180,j,1,n)=help
    397426        endif
    398427      endif
    … …  
    400429  ! MEAN SEA LEVEL PRESSURE
    401430        help=zsec4(nxfield*(ny-j-1)+i+1)
    402         if(i.le.i180) then
    403           msl(i179+i,j,1,n)=help
    404         else
    405           msl(i-i181,j,1,n)=help
     431        if(i.lt.i180) then
     432          msl(i180+i,j,1,n)=help
     433        else
     434          msl(i-i180,j,1,n)=help
    406435        endif
    407436      endif
    … …  
    409438  ! TOTAL CLOUD COVER
    410439        help=zsec4(nxfield*(ny-j-1)+i+1)
    411         if(i.le.i180) then
    412           tcc(i179+i,j,1,n)=help
    413         else
    414           tcc(i-i181,j,1,n)=help
     440        if(i.lt.i180) then
     441          tcc(i180+i,j,1,n)=help
     442        else
     443          tcc(i-i180,j,1,n)=help
    415444        endif
    416445      endif
    417446      if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
    418            (isec1(8).eq.10)) then
     447         (nint(xsec18).eq.10)) then
    419448  ! 10 M U VELOCITY
    420449        help=zsec4(nxfield*(ny-j-1)+i+1)
    421         if(i.le.i180) then
    422           u10(i179+i,j,1,n)=help
    423         else
    424           u10(i-i181,j,1,n)=help
     450        if(i.lt.i180) then
     451          u10(i180+i,j,1,n)=help
     452        else
     453          u10(i-i180,j,1,n)=help
    425454        endif
    426455      endif
    427456      if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
    428            (isec1(8).eq.10)) then
     457         (nint(xsec18).eq.10)) then
    429458  ! 10 M V VELOCITY
    430459        help=zsec4(nxfield*(ny-j-1)+i+1)
    431         if(i.le.i180) then
    432           v10(i179+i,j,1,n)=help
    433         else
    434           v10(i-i181,j,1,n)=help
     460        if(i.lt.i180) then
     461          v10(i180+i,j,1,n)=help
     462        else
     463          v10(i-i180,j,1,n)=help
    435464        endif
    436465      endif
    437466      if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
    438            (isec1(8).eq.02)) then
     467         (nint(xsec18).eq.2)) then
    439468  ! 2 M TEMPERATURE
    440469        help=zsec4(nxfield*(ny-j-1)+i+1)
    441         if(i.le.i180) then
    442           tt2(i179+i,j,1,n)=help
    443         else
    444           tt2(i-i181,j,1,n)=help
     470        if(i.lt.i180) then
     471          tt2(i180+i,j,1,n)=help
     472        else
     473          tt2(i-i180,j,1,n)=help
    445474        endif
    446475      endif
    447476      if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
    448            (isec1(8).eq.02)) then
     477         (nint(xsec18).eq.2)) then
    449478  ! 2 M DEW POINT TEMPERATURE
    450479        help=zsec4(nxfield*(ny-j-1)+i+1)
    451         if(i.le.i180) then
    452           td2(i179+i,j,1,n)=help
    453         else
    454           td2(i-i181,j,1,n)=help
     480        if(i.lt.i180) then
     481          td2(i180+i,j,1,n)=help
     482        else
     483          td2(i-i180,j,1,n)=help
    455484        endif
    456485      endif
    … …  
    458487  ! LARGE SCALE PREC.
    459488        help=zsec4(nxfield*(ny-j-1)+i+1)
    460         if(i.le.i180) then
    461           lsprec(i179+i,j,1,n)=help
    462         else
    463           lsprec(i-i181,j,1,n)=help
     489        if(i.lt.i180) then
     490          lsprec(i180+i,j,1,n)=help
     491        else
     492          lsprec(i-i180,j,1,n)=help
    464493        endif
    465494      endif
    … …  
    467496  ! CONVECTIVE PREC.
    468497        help=zsec4(nxfield*(ny-j-1)+i+1)
    469         if(i.le.i180) then
    470           convprec(i179+i,j,1,n)=help
    471         else
    472           convprec(i-i181,j,1,n)=help
     498        if(i.lt.i180) then
     499          convprec(i180+i,j,1,n)=help
     500        else
     501          convprec(i-i180,j,1,n)=help
    473502        endif
    474503      endif
    … …  
    476505  ! TOPOGRAPHY
    477506        help=zsec4(nxfield*(ny-j-1)+i+1)
    478         if(i.le.i180) then
    479           oro(i179+i,j)=help
    480           excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
    481         else
    482           oro(i-i181,j)=help
    483           excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     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
    484513        endif
    485514      endif
    … …  
    487516  ! LAND SEA MASK
    488517        help=zsec4(nxfield*(ny-j-1)+i+1)
    489         if(i.le.i180) then
    490           lsm(i179+i,j)=help
    491         else
    492           lsm(i-i181,j)=help
     518        if(i.lt.i180) then
     519          lsm(i180+i,j)=help
     520        else
     521          lsm(i-i180,j)=help
    493522        endif
    494523      endif
    … …  
    496525  ! MIXING HEIGHT
    497526        help=zsec4(nxfield*(ny-j-1)+i+1)
    498         if(i.le.i180) then
    499           hmix(i179+i,j,1,n)=help
    500         else
    501           hmix(i-i181,j,1,n)=help
     527        if(i.lt.i180) then
     528          hmix(i180+i,j,1,n)=help
     529        else
     530          hmix(i-i180,j,1,n)=help
    502531        endif
    503532      endif
    504533      if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
    505            (isec1(8).eq.02)) then
     534         (nint(xsec18).eq.02)) then
    506535  ! 2 M RELATIVE HUMIDITY
    507536        help=zsec4(nxfield*(ny-j-1)+i+1)
    508         if(i.le.i180) then
    509           qvh2(i179+i,j)=help
    510         else
    511           qvh2(i-i181,j)=help
     537        if(i.lt.i180) then
     538          qvh2(i180+i,j)=help
     539        else
     540          qvh2(i-i180,j)=help
    512541        endif
    513542      endif
    … …  
    515544  ! TEMPERATURE LOWEST SIGMA LEVEL
    516545        help=zsec4(nxfield*(ny-j-1)+i+1)
    517         if(i.le.i180) then
    518           tlev1(i179+i,j)=help
    519         else
    520           tlev1(i-i181,j)=help
     546        if(i.lt.i180) then
     547          tlev1(i180+i,j)=help
     548        else
     549          tlev1(i-i180,j)=help
    521550        endif
    522551      endif
    … …  
    524553  ! U VELOCITY LOWEST SIGMA LEVEL
    525554        help=zsec4(nxfield*(ny-j-1)+i+1)
    526         if(i.le.i180) then
    527           ulev1(i179+i,j)=help
    528         else
    529           ulev1(i-i181,j)=help
     555        if(i.lt.i180) then
     556          ulev1(i180+i,j)=help
     557        else
     558          ulev1(i-i180,j)=help
    530559        endif
    531560      endif
    … …  
    533562  ! V VELOCITY LOWEST SIGMA LEVEL
    534563        help=zsec4(nxfield*(ny-j-1)+i+1)
    535         if(i.le.i180) then
    536           vlev1(i179+i,j)=help
    537         else
    538           vlev1(i-i181,j)=help
     564        if(i.lt.i180) then
     565          vlev1(i180+i,j)=help
     566        else
     567          vlev1(i-i180,j)=help
    539568        endif
    540569      endif
    541570! SEC & IP 12/2018 read GFS clouds
    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             do ii=1,nuvz
    545               if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
    546             end do
    547         endif
    548         help=zsec4(nxfield*(ny-j-1)+i+1)
    549         if(i.le.i180) then
    550           clwch(i179+i,j,numpclwch,n)=help
    551         else
    552           clwch(i-i181,j,numpclwch,n)=help
    553         endif
    554         readclouds=.true.
    555         sumclouds=.true.
    556 !        readclouds=.false.
    557 !       sumclouds=.false.
    558       endif
     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
    559587
    560588
    … …  
    592620  ! CONVERT TP TO LSP (GRIB2 only)
    593621  if (gribVer.eq.2) then
    594     do j=0,nymin1
    595     do i=0,nxfield-1
    596      if(i.le.i180) then
    597      if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
    598          lsprec(i179+i,j,1,n)= &
    599               lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
    600      else
    601          lsprec(i179+i,j,1,n)=0
    602      endif
    603      else
    604      if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
    605           lsprec(i-i181,j,1,n)= &
    606                lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
    607      else
    608           lsprec(i-i181,j,1,n)=0
    609      endif
    610      endif
    611     enddo
    612     enddo
     622    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) )
    613623  endif
    614624  !HSO end edits
    … …  
    622632        help=qvh(i,j,k,n)
    623633        temp=tth(i,j,k,n)
     634        if (temp .eq. 0.0) then
     635           write (*, *) i, j, k, n
     636           temp = 273.0
     637        endif
    624638        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
    625639        elev=ew(temp)*help/100.0
    … …  
    637651        help=qvh2(i,j)
    638652        temp=tt2(i,j,1,n)
     653        if (temp .eq. 0.0) then
     654           write (*, *) i, j, n
     655           temp = 273.0
     656        endif
    639657        elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
    640658        td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.

  • src/richardson.f90

    • Property mode changed from 100644 to 100755
    r92fab65 r467460a  
    141141    zold=z
    142142  end do
    143   k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
     143  ! k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
    14414420   continue
    145145
    146146  ! Determine Richardson number between the critical levels
    147147  !********************************************************
     148  ! JMA: It may happen that k >= nuvz:
     149  ! JMA: In that case: k is set to nuvz
     150
     151  if (k .gt. nuvz) k = nuvz ! JMA
    148152
    149153  zl1=zold
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG