Changeset 4c52d0b in flexpart.git


Ignore:
Timestamp:
Apr 8, 2021, 9:27:43 AM (3 years ago)
Author:
Espen Sollum <eso@…>
Branches:
dev
Children:
b1f28c3
Parents:
aa939a9 (diff), 4138764 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'dev' into GFS_025

Files:
14 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    rf976c5d r969c1fe  
     1log*
     2noh*
     3diff*
    14FP_ecmwf_gfortran*
    25*.o
    … …  
    710FLEXPART
    811FLEXPART_MPI
     12FLEX*

  • src/calcfluxes.f90

    r92fab65 r4138764  
    6464       (jyave.le.numygrid-1)) then
    6565    do kz=1,numzgrid                ! determine height of cell
    66       if (outheighthalf(kz).gt.zold) goto 11
     66!      if (outheighthalf(kz).gt.zold) goto 11
     67      if (outheight(kz).gt.zold) goto 11 !sec, use upper layer instead of mid layer
    6768    end do
    686911   k1=min(numzgrid,kz)
    6970    do kz=1,numzgrid                ! determine height of cell
    70       if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
     71!      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
     72      if (outheight(kz).gt.ztra1(jpart)) goto 21
    7173    end do
    727421   k2=min(numzgrid,kz)
    … …  
    9698
    9799    if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
    98       ix1=int((xold*dx+xoutshift)/dxout+0.5)
    99       ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
     100      ix1=int((xold*dx+xoutshift)/dxout)              ! flux throught the western boundary (sec)
     101      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout)      ! flux throught the western boundary (sec)
     102!      ix1=int((xold*dx+xoutshift)/dxout+0.5)
     103!      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
    100104      do k=1,nspec
    101105        do ix=ix1,ix2-1
    … …  
    145149  if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
    146150       (ixave.le.numxgrid-1)) then
    147     jy1=int((yold*dy+youtshift)/dyout+0.5)
    148     jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
     151    jy1=int((yold*dy+youtshift)/dyout)         ! flux throught the southern boundary (sec)
     152    jy2=int((ytra1(jpart)*dy+youtshift)/dyout) ! flux throught the southern boundary (sec)
     153!    jy1=int((yold*dy+youtshift)/dyout+0.5)
     154!    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
    149155
    150156    do k=1,nspec

  • src/init_domainfill.f90

    r92fab65 r4138764  
    209209                   ztra1(numpart+jj)=height(nz)-0.5
    210210
    211 
    212211! Interpolate PV to the particle position
    213212!****************************************
    … …  
    216215              ixp=ixm+1
    217216              jyp=jym+1
     217              if (jyp.gt.180) then
     218                 write (*,*) 'init_domainfill, over: ',jyp,jym,ytra1(numpart+jj),jy,ran1(idummy),ny
     219                 jyp=jym
     220              endif
    218221              ddx=xtra1(numpart+jj)-real(ixm)
    219222              ddy=ytra1(numpart+jj)-real(jym)

  • src/netcdf_output_mod.f90

    • Property mode changed from 100644 to 100755
    r759df5f r4138764  
    3636  use outg_mod,  only: outheight,oroout,densityoutgrid,factor3d,volume,&
    3737                       wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
    38                        area,arean,volumen, orooutn, p0out, t0out
     38                       area,arean,volumen, orooutn, areaeast, areanorth, p0out, t0out
    3939  use par_mod,   only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
    4040                       unitoutrecept,unitoutreceptppt, nxmax,unittmp, &
    4141                       write_p0t0
    42   use com_mod,   only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
     42  use com_mod,   only: path,length,ldirect,bdate,ibdate,ibtime,iedate,ietime, &
    4343                       loutstep,loutaver,loutsample,outlon0,outlat0,&
    4444                       numxgrid,numygrid,dxout,dyout,numzgrid, height, &
    … …  
    6767
    6868  public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
    69        &concoutput_nest_netcdf,concoutput_surf_netcdf
     69       &concoutput_nest_netcdf,concoutput_surf_netcdf,fluxoutput_netcdf
    7070
    7171!  include 'netcdf.inc'
    … …  
    15201520end subroutine concoutput_surf_nest_netcdf
    15211521
     1522subroutine fluxoutput_netcdf(itime)
     1523     
     1524  !                            i
     1525  !*****************************************************************************
     1526  !                                                                            *
     1527  !     Output of the gridded fluxes.                                          *
     1528  !     Eastward, westward, northward, southward, upward and downward gross    *
     1529  !     fluxes are written to output file in either sparse matrix or grid dump *
     1530  !     format, whichever is more efficient.                                   *
     1531  !                                                                            *
     1532  !     Author: A. Stohl                                                       *
     1533  !                                                                            *
     1534  !     04 April 2000                                                          *
     1535  !     netcdfoutput S. Eckhardt, 2020                                         *
     1536  !                                                                            *
     1537  !*****************************************************************************
     1538
     1539   use flux_mod
     1540
     1541  implicit none
     1542
     1543  real(kind=dp) :: jul
     1544  integer :: itime,ix,jy,kz,ks,nage,jjjjmmdd,ihmmss,kp,i
     1545 
     1546  character(len=255) :: ncfname
     1547  character :: adate*8,atime*6,timeunit*32,anspec*3
     1548
     1549  integer :: ncid
     1550  integer :: timeDimID, latDimID, lonDimID, levDimID
     1551  integer :: nspecDimID, npointDimID, nageclassDimID, ncharDimID, pointspecDimID
     1552  integer :: tID, lonID, latID, levID, lageID, fluxID
     1553  integer, dimension(6)       :: dIDs
     1554  integer :: cache_size
     1555  real, allocatable, dimension(:) :: coord
     1556
     1557
     1558  ! Determine current calendar date, needed for the file name
     1559  !**********************************************************
     1560
     1561  jul=bdate+real(itime,kind=dp)/86400._dp
     1562  call caldate(jul,jjjjmmdd,ihmmss)
     1563  write(adate,'(i8.8)') jjjjmmdd
     1564  write(atime,'(i6.6)') ihmmss
     1565
     1566  ncfname=path(2)(1:length(2))//'grid_flux_'//adate// &
     1567       atime//'.nc'
     1568
     1569  ! setting cache size in bytes. It is set to 4 times the largest data block that is written
     1570  !   size_type x nx x ny x nz
     1571  ! create file
     1572
     1573  cache_size = 16 * numxgrid * numygrid * numzgrid
     1574  call nf90_err(nf90_create(trim(ncfname), cmode = nf90_hdf5, ncid = ncid, &
     1575    cache_size = cache_size)) 
     1576
     1577  ! create dimensions:
     1578  !*************************
     1579  ! time
     1580  call nf90_err(nf90_def_dim(ncid, 'time', nf90_unlimited, timeDimID))
     1581  timeunit = 'seconds since '//adate(1:4)//'-'//adate(5:6)// &
     1582     '-'//adate(7:8)//' '//atime(1:2)//':'//atime(3:4)
     1583
     1584  call nf90_err(nf90_def_dim(ncid, 'longitude', numxgrid, lonDimID))
     1585  call nf90_err(nf90_def_dim(ncid, 'latitude', numygrid, latDimID))
     1586  call nf90_err(nf90_def_dim(ncid, 'height', numzgrid, levDimID))
     1587  call nf90_err(nf90_def_dim(ncid, 'numspec', nspec, nspecDimID))
     1588  call nf90_err(nf90_def_dim(ncid, 'pointspec', maxpointspec_act, pointspecDimID))
     1589  call nf90_err(nf90_def_dim(ncid, 'nageclass', nageclass, nageclassDimID))
     1590  call nf90_err(nf90_def_dim(ncid, 'nchar', 45, ncharDimID))
     1591  call nf90_err(nf90_def_dim(ncid, 'numpoint', numpoint, npointDimID))
     1592
     1593
     1594  ! create variables
     1595  !*************************
     1596
     1597  ! time
     1598  call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ timeDimID /), tID))
     1599  call nf90_err(nf90_put_att(ncid, tID, 'units', timeunit))
     1600  call nf90_err(nf90_put_att(ncid, tID, 'calendar', 'proleptic_gregorian'))
     1601  timeID = tID
     1602
     1603  ! lon
     1604  call nf90_err(nf90_def_var(ncid, 'longitude', nf90_float, (/ lonDimID /), lonID))
     1605  call nf90_err(nf90_put_att(ncid, lonID, 'long_name', 'longitude in degree east'))
     1606  call nf90_err(nf90_put_att(ncid, lonID, 'axis', 'Lon'))
     1607  call nf90_err(nf90_put_att(ncid, lonID, 'units', 'degrees_east'))
     1608  call nf90_err(nf90_put_att(ncid, lonID, 'standard_name', 'grid_longitude'))
     1609  call nf90_err(nf90_put_att(ncid, lonID, 'description', 'grid cell centers'))
     1610
     1611  ! lat
     1612  call nf90_err(nf90_def_var(ncid, 'latitude', nf90_float, (/ latDimID /), latID))
     1613  call nf90_err(nf90_put_att(ncid, latID, 'long_name', 'latitude in degree north'))
     1614  call nf90_err(nf90_put_att(ncid, latID, 'axis', 'Lat'))
     1615  call nf90_err(nf90_put_att(ncid, latID, 'units', 'degrees_north'))
     1616  call nf90_err(nf90_put_att(ncid, latID, 'standard_name', 'grid_latitude'))
     1617  call nf90_err(nf90_put_att(ncid, latID, 'description', 'grid cell centers'))
     1618
     1619  ! height
     1620  call nf90_err(nf90_def_var(ncid, 'height', nf90_float, (/ levDimID /), levID))
     1621  call nf90_err(nf90_put_att(ncid, levID, 'units', 'meters'))
     1622  call nf90_err(nf90_put_att(ncid, levID, 'positive', 'up'))
     1623  call nf90_err(nf90_put_att(ncid, levID, 'standard_name', 'height'))
     1624
     1625  if (.not.allocated(coord)) allocate(coord(numxgrid))
     1626  do i = 1,numxgrid
     1627     coord(i) = outlon0 + (i-0.5)*dxout
     1628  enddo
     1629  call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgrid)))
     1630  deallocate(coord)
     1631  if (.not.allocated(coord)) allocate(coord(numygrid))
     1632  do i = 1,numygrid
     1633     coord(i) = outlat0 + (i-0.5)*dyout
     1634  enddo
     1635  call nf90_err(nf90_put_var(ncid, latID, coord(1:numygrid)))
     1636  deallocate(coord)
     1637  call nf90_err(nf90_put_var(ncid, levID, outheight(1:numzgrid)))
     1638
     1639  ! write time, one field per time - different to the others!
     1640  call nf90_err(nf90_put_var( ncid, timeID, itime, (/ 1 /)))
     1641 
     1642  dIDs = (/ londimid, latdimid, levdimid, timedimid, pointspecdimid, nageclassdimid /)
     1643
     1644  do ks=1,nspec
     1645    do kp=1,maxpointspec_act
     1646      do nage=1,nageclass
     1647 
     1648         write(anspec,'(i3.3)') ks
     1649
     1650 ! East Flux
     1651        call nf90_err(nf90_def_var(ncid,'flux_east_'//anspec, nf90_float, dIDs, &
     1652             fluxID))
     1653
     1654        do jy=0,numygrid-1
     1655          do ix=0,numxgrid-1
     1656            do kz=1, numzgrid
     1657               grid(ix,jy,kz)=flux(1,ix,jy,kz,ks,kp,nage)
     1658            end do
     1659          end do
     1660        end do
     1661
     1662        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1663            /areaeast(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
     1664           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1665
     1666 ! West Flux
     1667        call nf90_err(nf90_def_var(ncid,'flux_west_'//anspec, nf90_float, dIDs, &
     1668             fluxID))
     1669
     1670        do jy=0,numygrid-1
     1671          do ix=0,numxgrid-1
     1672            do kz=1, numzgrid
     1673               grid(ix,jy,kz)=flux(2,ix,jy,kz,ks,kp,nage)
     1674            end do
     1675          end do
     1676        end do
     1677
     1678        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1679            /areaeast(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
     1680           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1681
     1682 ! North Flux
     1683        call nf90_err(nf90_def_var(ncid,'flux_north_'//anspec, nf90_float, dIDs, &
     1684             fluxID))
     1685
     1686        do jy=0,numygrid-1
     1687          do ix=0,numxgrid-1
     1688            do kz=1, numzgrid
     1689               grid(ix,jy,kz)=flux(4,ix,jy,kz,ks,kp,nage)
     1690            end do
     1691          end do
     1692        end do
     1693
     1694        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1695            /areanorth(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
     1696           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1697
     1698 ! South Flux
     1699        call nf90_err(nf90_def_var(ncid,'flux_south_'//anspec, nf90_float, dIDs, &
     1700             fluxID))
     1701
     1702        do jy=0,numygrid-1
     1703          do ix=0,numxgrid-1
     1704            do kz=1, numzgrid
     1705               grid(ix,jy,kz)=flux(3,ix,jy,kz,ks,kp,nage)
     1706            end do
     1707          end do
     1708        end do
     1709
     1710        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1711            /areanorth(0:numxgrid-1,0:numygrid-1,1:numzgrid)/loutstep,&
     1712           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1713
     1714 ! Up Flux
     1715        call nf90_err(nf90_def_var(ncid,'flux_up_'//anspec, nf90_float, dIDs, &
     1716             fluxID))
     1717
     1718        do jy=0,numygrid-1
     1719          do ix=0,numxgrid-1
     1720            do kz=1, numzgrid
     1721               grid(ix,jy,kz)=flux(5,ix,jy,kz,ks,kp,nage)/area(ix,jy)
     1722            end do
     1723          end do
     1724        end do
     1725
     1726        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1727            /loutstep,&
     1728           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1729
     1730 ! Down Flux
     1731        call nf90_err(nf90_def_var(ncid,'flux_down_'//anspec, nf90_float, dIDs, &
     1732             fluxID))
     1733
     1734        do jy=0,numygrid-1
     1735          do ix=0,numxgrid-1
     1736            do kz=1, numzgrid
     1737               grid(ix,jy,kz)=flux(6,ix,jy,kz,ks,kp,nage)/area(ix,jy)
     1738            end do
     1739          end do
     1740        end do
     1741
     1742        call nf90_err(nf90_put_var(ncid,fluxid,1.e12*grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)&
     1743            /loutstep,&
     1744           (/ 1,1,1,1,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
     1745
     1746      end do
     1747    end do
     1748  end do
     1749
     1750  ! Close netCDF file
     1751  call nf90_err(nf90_close(ncid))
     1752
     1753  ! Reinitialization of grid
     1754  !*************************
     1755
     1756  do ks=1,nspec
     1757    do kp=1,maxpointspec_act
     1758      do jy=0,numygrid-1
     1759        do ix=0,numxgrid-1
     1760          do kz=1,numzgrid
     1761            do nage=1,nageclass
     1762              do i=1,6
     1763                flux(i,ix,jy,kz,ks,kp,nage)=0.
     1764              end do
     1765            end do
     1766          end do
     1767        end do
     1768      end do
     1769    end do
     1770  end do
     1771
     1772
     1773end subroutine fluxoutput_netcdf
     1774
    15221775end module netcdf_output_mod
    15231776

  • src/readspecies.f90

    r92fab65 r4138764  
    4242  implicit none
    4343
    44   integer :: i, pos_spec,j
     44  integer :: i, pos_spec,j,icheck_dow_hour
    4545  integer :: idow,ihour,id_spec
    4646  character(len=3) :: aspecnumb
    … …  
    229229    ohnconst(pos_spec)=pohnconst
    230230
     231
     232    icheck_dow_hour=0
    231233    do j=1,24     ! 24 hours, starting with 0-1 local time
    232234      area_hour(pos_spec,j)=parea_hour(j)
    233235      point_hour(pos_spec,j)=ppoint_hour(j)
     236      if (parea_hour(j).ne.1 .or. ppoint_hour(j).ne.1) icheck_dow_hour=1
    234237    end do
    235238    do j=1,7      ! 7 days of the week, starting with Monday
    236239      area_dow(pos_spec,j)=parea_dow(j)
    237240      point_dow(pos_spec,j)=ppoint_dow(j)
     241      if (parea_dow(j).ne.1 .or. ppoint_dow(j).ne.1) icheck_dow_hour=1
    238242    end do
    239243
    … …  
    35736120 continue
    358362
     363  if ((icheck_dow_hour.eq.1).and.(ldirect.lt.0)) then
     364    write(*,*) '#### FLEXPART MODEL WARNING                     ####'
     365    write(*,*) '#### The variation for an emission release      ####'
     366    write(*,*) '#### will have no effect in backward mode       ####' 
     367  endif
    359368
    36036922 close(unitspecies)

  • src/timemanager.f90

    ra803521 r4138764  
    8585#ifdef USE_NCF
    8686  use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
    87        &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
     87       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf,fluxoutput_netcdf
    8888#endif
    8989
    … …  
    433433        endif
    434434        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    435         if (iflux.eq.1) call fluxoutput(itime)
     435        if ((iflux.eq.1).and.(lnetcdfout.eq.0)) call fluxoutput(itime)
     436#ifdef USE_NCF
     437        if ((iflux.eq.1).and.(lnetcdfout.eq.1)) call fluxoutput_netcdf(itime)
     438#endif
    436439        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    437440 

  • src/timemanager_mpi.f90

    rf06b72a r4138764  
    556556        endif
    557557        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    558         if (iflux.eq.1) call fluxoutput(itime)
     558        if ((iflux.eq.1).and.(lnetcdfout.eq.0)) call fluxoutput(itime)
     559        if ((iflux.eq.1).and.(lnetcdfout.eq.1)) call fluxoutput_netcdf(itime)
    559560        if (mp_measure_time) call mpif_mtime('iotime',1)
    560561

  • src/calcmatrix.f90

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

  • src/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 raa939a9  
    147147! GFS
    148148   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
    149    integer :: nxshift=0 ! shift not fixed for the executable
     149  ! GFS 0.25
     150  ! integer,parameter :: nxmax=1441,nymax=721,nuvzmax=138,nwzmax=138,nzmax=138
     151  integer :: nxshift=0 ! shift not fixed for the executable
    150152
    151153

  • 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 r9ca6e38  
    44subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
    55
    6   !***********************************************************************
    7   !*                                                                     *
    8   !*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
    9   !*                                                                     *
    10   !***********************************************************************
    11   !*                                                                     *
    12   !*             AUTHOR:      G. WOTAWA                                  *
    13   !*             DATE:        1997-08-05                                 *
    14   !*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
    15   !*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
    16   !*                     qvh (on eta coordinates) in common block        *
    17   !*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
    18   !*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
    19   !*                     data with the ECMWF grib_api library            *
    20   !*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    21   !*                                 ECMWF grib_api                      *
    22   !                                                                      *
    23   !   Unified ECMWF and GFS builds                                       *
    24   !   Marian Harustak, 12.5.2017                                         *
    25   !     - Renamed routine from readwind to readwind_gfs                  *
    26   !*                                                                     *
    27   !***********************************************************************
    28   !*                                                                     *
    29   !* DESCRIPTION:                                                        *
    30   !*                                                                     *
    31   !* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
    32   !* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
    33   !*                                                                     *
    34   !* INPUT:                                                              *
    35   !* indj               indicates number of the wind field to be read in *
    36   !* n                  temporal index for meteorological fields (1 to 3)*
    37   !*                                                                     *
    38   !* IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
    39   !*                                                                     *
    40   !* wfname             File name of data to be read in                  *
    41   !* nx,ny,nuvz,nwz     expected field dimensions                        *
    42   !* nlev_ec            number of vertical levels ecmwf model            *
    43   !* uu,vv,ww           wind fields                                      *
    44   !* tt,qv              temperature and specific humidity                *
    45   !* ps                 surface pressure                                 *
    46   !*                                                                     *
    47   !***********************************************************************
     6!***********************************************************************
     7!*                                                                     *
     8!*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
     9!*                                                                     *
     10!***********************************************************************
     11!*                                                                     *
     12!*             AUTHOR:      G. WOTAWA                                  *
     13!*             DATE:        1997-08-05                                 *
     14!*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
     15!*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
     16!*                     qvh (on eta coordinates) in common block        *
     17!*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
     18!*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
     19!*                     data with the ECMWF grib_api library            *
     20!*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
     21!*                                 ECMWF grib_api                      *
     22!                                                                      *
     23!   Unified ECMWF and GFS builds                                       *
     24!   Marian Harustak, 12.5.2017                                         *
     25!     - Renamed routine from readwind to readwind_gfs                  *
     26!*                                                                     *
     27!***********************************************************************
     28!*                                                                     *
     29!* DESCRIPTION:                                                        *
     30!*                                                                     *
     31!* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
     32!* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
     33!*                                                                     *
     34!* INPUT:                                                              *
     35!* indj               indicates number of the wind field to be read in *
     36!* n                  temporal index for meteorological fields (1 to 3)*
     37!*                                                                     *
     38!* IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
     39!*                                                                     *
     40!* wfname             File name of data to be read in                  *
     41!* nx,ny,nuvz,nwz     expected field dimensions                        *
     42!* nlev_ec            number of vertical levels ecmwf model            *
     43!* uu,vv,ww           wind fields                                      *
     44!* tt,qv              temperature and specific humidity                *
     45!* ps                 surface pressure                                 *
     46!*                                                                     *
     47!***********************************************************************
    4848
    4949  use eccodes
    … …  
    5353  implicit none
    5454
    55   !HSO  new parameters for grib_api
     55!HSO  new parameters for grib_api
    5656  integer :: ifile
    5757  integer :: iret
    5858  integer :: igrib
    59   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
    60   !HSO end edits
     59  integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
     60!HSO end edits
    6161  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    6262  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    6464  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
    6565
    66   ! NCEP
     66! NCEP
    6767  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
    6868  real :: help, temp, ew
    … …  
    7272  real :: qvh2(0:nxmax-1,0:nymax-1)
    7373
    74   integer :: i179,i180,i181
    75 
    76   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    77   !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
     74  integer :: i180
     75
     76! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
     77!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
    7878
    7979  integer :: isec1(8),isec2(3)
     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
    … …  
    8788  logical :: hflswitch,strswitch
    8889
    89   !HSO  for grib api error messages
     90!HSO  for grib api error messages
    9091  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    9192  character(len=20) :: gribFunction = 'readwind_gfs'
    … …  
    100101
    101102
    102   ! OPENING OF DATA FILE (GRIB CODE)
    103 
    104   !HSO
     103! OPENING OF DATA FILE (GRIB CODE)
     104
     105!HSO
    105106  call grib_open_file(ifile,path(3)(1:length(3)) &
    106          //trim(wfname(indj)),'r',iret)
     107       //trim(wfname(indj)),'r',iret)
    107108  if (iret.ne.GRIB_SUCCESS) then
    108109    goto 888   ! ERROR DETECTED
    109110  endif
    110   !turn on support for multi fields messages
     111!turn on support for multi fields messages
    111112  call grib_multi_support_on
    112113
    … …  
    118119  numpclwch=0
    119120  ifield=0
    120 10   ifield=ifield+1
    121   !
    122   ! GET NEXT FIELDS
    123   !
     12110 ifield=ifield+1
     122!
     123! GET NEXT FIELDS
     124!
    124125  call grib_new_from_file(ifile,igrib,iret)
    125126  if (iret.eq.GRIB_END_OF_FILE)  then
    … …  
    129130  endif
    130131
    131   !first see if we read GRIB1 or GRIB2
     132!first see if we read GRIB1 or GRIB2
    132133  call grib_get_int(igrib,'editionNumber',gribVer,iret)
    133134!  call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    135136  if (gribVer.eq.1) then ! GRIB Edition 1
    136137
    137   !read the grib1 identifiers
    138   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    139 !  call grib_check(iret,gribFunction,gribErrorMsg)
    140   call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
    141 !  call grib_check(iret,gribFunction,gribErrorMsg)
    142   call grib_get_int(igrib,'level',isec1(8),iret)
    143 !  call grib_check(iret,gribFunction,gribErrorMsg)
     138!read the grib1 identifiers
     139
     140    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     141    call grib_check(iret,gribFunction,gribErrorMsg)
     142    call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
     143    call grib_check(iret,gribFunction,gribErrorMsg)
     144    call grib_get_int(igrib,'level',isec1(8),iret)
     145    call grib_check(iret,gribFunction,gribErrorMsg)
     146!JMA / SH: isec1(8) not evaluated any more below
     147!b/c with GRIB 2 this may be a real variable
     148    xsec18 = real(isec1(8))
    144149
    145150  else ! GRIB Edition 2
    146151
    147   !read the grib2 identifiers
    148   call grib_get_string(igrib,'shortName',shortname,iret)
    149 
    150   call grib_get_int(igrib,'discipline',discipl,iret)
    151 !  call grib_check(iret,gribFunction,gribErrorMsg)
    152   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    153 !  call grib_check(iret,gribFunction,gribErrorMsg)
    154   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    155 !  call grib_check(iret,gribFunction,gribErrorMsg)
    156   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    157 !  call grib_check(iret,gribFunction,gribErrorMsg)
    158   call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
    159        valSurf,iret)
    160 !  call grib_check(iret,gribFunction,gribErrorMsg)
    161  
    162 !  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
    163   !convert to grib1 identifiers
    164   isec1(6)=-1
    165   isec1(7)=-1
    166   isec1(8)=-1
    167   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
    168     isec1(6)=11          ! indicatorOfParameter
    169     isec1(7)=100         ! indicatorOfTypeOfLevel
    170     isec1(8)=valSurf/100 ! level, convert to hPa
    171   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
    172     isec1(6)=33          ! indicatorOfParameter
    173     isec1(7)=100         ! indicatorOfTypeOfLevel
    174     isec1(8)=valSurf/100 ! level, convert to hPa
    175   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
    176     isec1(6)=34          ! indicatorOfParameter
    177     isec1(7)=100         ! indicatorOfTypeOfLevel
    178     isec1(8)=valSurf/100 ! level, convert to hPa
    179   elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
    180     isec1(6)=39          ! indicatorOfParameter
    181     isec1(7)=100         ! indicatorOfTypeOfLevel
    182     isec1(8)=valSurf/100 ! level, convert to hPa
    183   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
    184     isec1(6)=52          ! indicatorOfParameter
    185     isec1(7)=100         ! indicatorOfTypeOfLevel
    186     isec1(8)=valSurf/100 ! level, convert to hPa
    187   elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
    188     isec1(6)=52          ! indicatorOfParameter
    189     isec1(7)=105         ! indicatorOfTypeOfLevel
    190     isec1(8)=2
    191   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
    192     isec1(6)=11          ! indicatorOfParameter
    193     isec1(7)=105         ! indicatorOfTypeOfLevel
    194     isec1(8)=2
    195   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
    196     isec1(6)=33          ! indicatorOfParameter
    197     isec1(7)=105         ! indicatorOfTypeOfLevel
    198     isec1(8)=10
    199   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
    200     isec1(6)=34          ! indicatorOfParameter
    201     isec1(7)=105         ! indicatorOfTypeOfLevel
    202     isec1(8)=10
    203   elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
    204     isec1(6)=153         ! indicatorOfParameter
    205     isec1(7)=100         ! indicatorOfTypeOfLevel
    206     isec1(8)=valSurf/100 ! level, convert to hPa
    207   elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
    208     isec1(6)=2           ! indicatorOfParameter
    209     isec1(7)=102         ! indicatorOfTypeOfLevel
    210     isec1(8)=0
    211   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
    212     isec1(6)=1           ! indicatorOfParameter
    213     isec1(7)=1           ! indicatorOfTypeOfLevel
    214     isec1(8)=0
    215   elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
    216     isec1(6)=66          ! indicatorOfParameter
    217     isec1(7)=1           ! indicatorOfTypeOfLevel
    218     isec1(8)=0
    219   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
    220     isec1(6)=11          ! indicatorOfParameter
    221     isec1(7)=107         ! indicatorOfTypeOfLevel
    222     isec1(8)=0.995       ! lowest sigma level
    223   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
    224     isec1(6)=33          ! indicatorOfParameter
    225     isec1(7)=107         ! indicatorOfTypeOfLevel
    226     isec1(8)=0.995       ! lowest sigma level
    227   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
    228     isec1(6)=34          ! indicatorOfParameter
    229     isec1(7)=107         ! indicatorOfTypeOfLevel
    230     isec1(8)=0.995       ! lowest sigma level
    231   elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
    232     isec1(6)=7           ! indicatorOfParameter
    233     isec1(7)=1           ! indicatorOfTypeOfLevel
    234     isec1(8)=0
    235   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
    236        .and.(discipl.eq.2)) then ! LSM
    237     isec1(6)=81          ! indicatorOfParameter
    238     isec1(7)=1           ! indicatorOfTypeOfLevel
    239     isec1(8)=0
    240   elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
    241     isec1(6)=221         ! indicatorOfParameter
    242     isec1(7)=1           ! indicatorOfTypeOfLevel
    243     isec1(8)=0
    244   elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
    245     isec1(6)=62          ! indicatorOfParameter
    246     isec1(7)=1           ! indicatorOfTypeOfLevel
    247     isec1(8)=0
    248   elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
    249     isec1(6)=63          ! indicatorOfParameter
    250     isec1(7)=1           ! indicatorOfTypeOfLevel
    251     isec1(8)=0
    252   endif
     152!read the grib2 identifiers
     153   
     154    call grib_get_string(igrib,'shortName',shortname,iret)
     155
     156    call grib_get_int(igrib,'discipline',discipl,iret)
     157    call grib_check(iret,gribFunction,gribErrorMsg)
     158    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     159    call grib_check(iret,gribFunction,gribErrorMsg)
     160    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     161    call grib_check(iret,gribFunction,gribErrorMsg)
     162    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     163    call grib_check(iret,gribFunction,gribErrorMsg)
     164!
     165    call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
     166         valSurf,iret)
     167    call grib_check(iret,gribFunction,gribErrorMsg)
     168
     169!    write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
     170!convert to grib1 identifiers
     171    isec1(6:8)=-1
     172    xsec18  =-1.0
     173!JMA / SH: isec1(8) not evaluated any more below
     174!b/c with GRIB 2 this may be a real variable
     175    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
     176      isec1(6)=11          ! indicatorOfParameter
     177      isec1(7)=100         ! indicatorOfTypeOfLevel
     178      xsec18=valSurf/100.0 ! level, convert to hPa
     179    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
     180      isec1(6)=33          ! indicatorOfParameter
     181      isec1(7)=100         ! indicatorOfTypeOfLevel
     182      xsec18=valSurf/100.0 ! level, convert to hPa
     183    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
     184      isec1(6)=34          ! indicatorOfParameter
     185      isec1(7)=100         ! indicatorOfTypeOfLevel
     186      xsec18=valSurf/100.0 ! level, convert to hPa
     187    elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
     188      isec1(6)=39          ! indicatorOfParameter
     189      isec1(7)=100         ! indicatorOfTypeOfLevel
     190      xsec18=valSurf/100.0 ! level, convert to hPa
     191    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
     192      isec1(6)=52          ! indicatorOfParameter
     193      isec1(7)=100         ! indicatorOfTypeOfLevel
     194      xsec18=valSurf/100.0 ! level, convert to hPa
     195    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
     196      isec1(6)=52          ! indicatorOfParameter
     197      isec1(7)=105         ! indicatorOfTypeOfLevel
     198      xsec18=real(2)
     199    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
     200      isec1(6)=11          ! indicatorOfParameter
     201      isec1(7)=105         ! indicatorOfTypeOfLevel
     202      xsec18=real(2)
     203    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
     204      isec1(6)=33          ! indicatorOfParameter
     205      isec1(7)=105         ! indicatorOfTypeOfLevel
     206      xsec18=real(10)
     207    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
     208      isec1(6)=34          ! indicatorOfParameter
     209      isec1(7)=105         ! indicatorOfTypeOfLevel
     210      xsec18=real(10)
     211    elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
     212      isec1(6)=153         ! indicatorOfParameter
     213      isec1(7)=100         ! indicatorOfTypeOfLevel
     214      xsec18=valSurf/100.0 ! level, convert to hPa
     215    elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
     216      isec1(6)=2           ! indicatorOfParameter
     217      isec1(7)=102         ! indicatorOfTypeOfLevel
     218      xsec18=real(0)
     219    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
     220      isec1(6)=1           ! indicatorOfParameter
     221      isec1(7)=1           ! indicatorOfTypeOfLevel
     222      xsec18=real(0)
     223    elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
     224      isec1(6)=66          ! indicatorOfParameter
     225      isec1(7)=1           ! indicatorOfTypeOfLevel
     226      xsec18=real(0)
     227    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
     228      isec1(6)=11          ! indicatorOfParameter
     229      isec1(7)=107         ! indicatorOfTypeOfLevel
     230      xsec18=0.995         ! lowest sigma level
     231    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
     232      isec1(6)=33          ! indicatorOfParameter
     233      isec1(7)=107         ! indicatorOfTypeOfLevel
     234      xsec18=0.995         ! lowest sigma level
     235    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
     236      isec1(6)=34          ! indicatorOfParameter
     237      isec1(7)=107         ! indicatorOfTypeOfLevel
     238      xsec18=0.995         ! lowest sigma level
     239    elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
     240      isec1(6)=7           ! indicatorOfParameter
     241      isec1(7)=1           ! indicatorOfTypeOfLevel
     242      xsec18=real(0)
     243    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
     244         .and.(discipl.eq.2)) then ! LSM
     245      isec1(6)=81          ! indicatorOfParameter
     246      isec1(7)=1           ! indicatorOfTypeOfLevel
     247      xsec18=real(0)
     248    elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
     249      isec1(6)=221         ! indicatorOfParameter
     250      isec1(7)=1           ! indicatorOfTypeOfLevel
     251      xsec18=real(0)
     252    elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
     253      isec1(6)=62          ! indicatorOfParameter
     254      isec1(7)=1           ! indicatorOfTypeOfLevel
     255      xsec18=real(0)
     256    elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
     257      isec1(6)=63          ! indicatorOfParameter
     258      isec1(7)=1           ! indicatorOfTypeOfLevel
     259      xsec18=real(0)
     260    endif
    253261
    254262  endif ! gribVer
    255263
    256264  if (isec1(6).ne.-1) then
    257   !  get the size and data of the values array
     265!  get the size and data of the values array
    258266    call grib_get_real4_array(igrib,'values',zsec4,iret)
    259 !    call grib_check(iret,gribFunction,gribErrorMsg)
     267    call grib_check(iret,gribFunction,gribErrorMsg)
    260268  endif
    261269
    262270  if(ifield.eq.1) then
    263271
    264   !get the required fields from section 2
    265   !store compatible to gribex input
    266   call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    267        isec2(2),iret)
    268 !  call grib_check(iret,gribFunction,gribErrorMsg)
    269   call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
    270        isec2(3),iret)
    271 !  call grib_check(iret,gribFunction,gribErrorMsg)
    272   call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    273        xauxin,iret)
    274 !  call grib_check(iret,gribFunction,gribErrorMsg)
    275   call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
    276        yauxin,iret)
    277 !  call grib_check(iret,gribFunction,gribErrorMsg)
    278   xaux=xauxin+real(nxshift)*dx
    279   yaux=yauxin
    280 
    281   ! CHECK GRID SPECIFICATIONS
     272!get the required fields from section 2
     273!store compatible to gribex input
     274    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
     275         isec2(2),iret)
     276    call grib_check(iret,gribFunction,gribErrorMsg)
     277    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
     278         isec2(3),iret)
     279    call grib_check(iret,gribFunction,gribErrorMsg)
     280    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
     281         xauxin,iret)
     282    call grib_check(iret,gribFunction,gribErrorMsg)
     283    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
     284         yauxin,iret)
     285    call grib_check(iret,gribFunction,gribErrorMsg)
     286    xaux=xauxin+real(nxshift)*dx
     287    yaux=yauxin
     288
     289! CHECK GRID SPECIFICATIONS
    282290
    283291    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    284292    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    285     if(xaux.eq.0.) xaux=-179.0     ! NCEP DATA
     293    if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
    286294    xaux0=xlon0
    287295    yaux0=ylat0
    … …  
    290298    if(xaux0.lt.0.) xaux0=xaux0+360.
    291299    if(yaux0.lt.0.) yaux0=yaux0+360.
    292     if(abs(xaux-xaux0).gt.eps) &
    293          stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    294     if(abs(yaux-yaux0).gt.eps) &
    295          stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
    296   endif
    297   !HSO end of edits
    298 
    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
     300    if(abs(xaux-xaux0).gt.eps) then
     301      write (*, *) xaux, xaux0
     302      stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
     303    endif
     304    if(abs(yaux-yaux0).gt.eps) then
     305      write (*, *) yaux, yaux0
     306      stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     307    end if
     308  endif
     309!HSO end of edits
     310
     311  i180=nint(180./dx)
    306312
    307313  if (isec1(6).ne.-1) then
    308314
    309   do j=0,nymin1
    310     do i=0,nxfield-1
    311       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
    323         endif
    324       endif
    325       if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
    326   ! U VELOCITY
    327          if((i.eq.0).and.(j.eq.0)) then
    328             do ii=1,nuvz
    329               if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
    330             end do
    331         endif
    332         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
    337         endif
    338       endif
    339       if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
    340   ! V VELOCITY
    341          if((i.eq.0).and.(j.eq.0)) then
    342             do ii=1,nuvz
    343               if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
    344             end do
    345         endif
    346         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
    351         endif
    352       endif
    353       if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
    354   ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
    355          if((i.eq.0).and.(j.eq.0)) then
    356             do ii=1,nuvz
    357               if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
    358             end do
    359         endif
    360         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
    365         endif
    366       endif
    367       if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
    368   ! SURFACE PRESSURE
    369         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
    374         endif
    375       endif
    376       if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
    377   ! W VELOCITY
    378          if((i.eq.0).and.(j.eq.0)) then
    379             do ii=1,nuvz
    380               if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
    381             end do
    382         endif
    383         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
    388         endif
    389       endif
    390       if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
    391   ! SNOW DEPTH
    392         help=zsec4(nxfield*(ny-j-1)+i+1)
    393         if(i.le.i180) then
    394           sd(i179+i,j,1,n)=help
    395         else
    396           sd(i-i181,j,1,n)=help
    397         endif
    398       endif
    399       if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
    400   ! MEAN SEA LEVEL PRESSURE
    401         help=zsec4(nxfield*(ny-j-1)+i+1)
    402         if(i.le.i180) then
    403           msl(i179+i,j,1,n)=help
    404         else
    405           msl(i-i181,j,1,n)=help
    406         endif
    407       endif
    408       if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
    409   ! TOTAL CLOUD COVER
    410         help=zsec4(nxfield*(ny-j-1)+i+1)
    411         if(i.le.i180) then
    412           tcc(i179+i,j,1,n)=help
    413         else
    414           tcc(i-i181,j,1,n)=help
    415         endif
    416       endif
    417       if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
    418            (isec1(8).eq.10)) then
    419   ! 10 M U VELOCITY
    420         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
    425         endif
    426       endif
    427       if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
    428            (isec1(8).eq.10)) then
    429   ! 10 M V VELOCITY
    430         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
    435         endif
    436       endif
    437       if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
    438            (isec1(8).eq.02)) then
    439   ! 2 M TEMPERATURE
    440         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
    445         endif
    446       endif
    447       if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
    448            (isec1(8).eq.02)) then
    449   ! 2 M DEW POINT TEMPERATURE
    450         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
    455         endif
    456       endif
    457       if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
    458   ! LARGE SCALE PREC.
    459         help=zsec4(nxfield*(ny-j-1)+i+1)
    460         if(i.le.i180) then
    461           lsprec(i179+i,j,1,n)=help
    462         else
    463           lsprec(i-i181,j,1,n)=help
    464         endif
    465       endif
    466       if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
    467   ! CONVECTIVE PREC.
    468         help=zsec4(nxfield*(ny-j-1)+i+1)
    469         if(i.le.i180) then
    470           convprec(i179+i,j,1,n)=help
    471         else
    472           convprec(i-i181,j,1,n)=help
    473         endif
    474       endif
    475       if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
    476   ! TOPOGRAPHY
    477         help=zsec4(nxfield*(ny-j-1)+i+1)
    478         if(i.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
    484         endif
    485       endif
    486       if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
    487   ! LAND SEA MASK
    488         help=zsec4(nxfield*(ny-j-1)+i+1)
    489         if(i.le.i180) then
    490           lsm(i179+i,j)=help
    491         else
    492           lsm(i-i181,j)=help
    493         endif
    494       endif
    495       if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
    496   ! MIXING HEIGHT
    497         help=zsec4(nxfield*(ny-j-1)+i+1)
    498         if(i.le.i180) then
    499           hmix(i179+i,j,1,n)=help
    500         else
    501           hmix(i-i181,j,1,n)=help
    502         endif
    503       endif
    504       if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
    505            (isec1(8).eq.02)) then
    506   ! 2 M RELATIVE HUMIDITY
    507         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
    512         endif
    513       endif
    514       if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
    515   ! TEMPERATURE LOWEST SIGMA LEVEL
    516         help=zsec4(nxfield*(ny-j-1)+i+1)
    517         if(i.le.i180) then
    518           tlev1(i179+i,j)=help
    519         else
    520           tlev1(i-i181,j)=help
    521         endif
    522       endif
    523       if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
    524   ! U VELOCITY LOWEST SIGMA LEVEL
    525         help=zsec4(nxfield*(ny-j-1)+i+1)
    526         if(i.le.i180) then
    527           ulev1(i179+i,j)=help
    528         else
    529           ulev1(i-i181,j)=help
    530         endif
    531       endif
    532       if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
    533   ! V VELOCITY LOWEST SIGMA LEVEL
    534         help=zsec4(nxfield*(ny-j-1)+i+1)
    535         if(i.le.i180) then
    536           vlev1(i179+i,j)=help
    537         else
    538           vlev1(i-i181,j)=help
    539         endif
    540       endif
     315! write (*, *) 'nxfield: ', nxfield, i180
     316
     317    do j=0,nymin1
     318      do i=0,nxfield-1
     319        if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
     320! TEMPERATURE
     321          if((i.eq.0).and.(j.eq.0)) then
     322            numpt=minloc(abs(xsec18*100.0-akz),dim=1)
     323          endif
     324          help=zsec4(nxfield*(ny-j-1)+i+1)
     325          if (help.eq.0) then
     326            write (*, *) 'i, j: ', i, j
     327            stop 'help == 0.0'
     328          endif
     329          if(i.lt.i180) then
     330            tth(i180+i,j,numpt,n)=help
     331          else
     332            tth(i-i180,j,numpt,n)=help
     333          endif
     334        endif
     335        if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
     336! U VELOCITY
     337          if((i.eq.0).and.(j.eq.0)) then
     338            numpu=minloc(abs(xsec18*100.0-akz),dim=1)
     339          endif
     340          help=zsec4(nxfield*(ny-j-1)+i+1)
     341          if(i.lt.i180) then
     342            uuh(i180+i,j,numpu)=help
     343          else
     344            uuh(i-i180,j,numpu)=help
     345          endif
     346        endif
     347        if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
     348! V VELOCITY
     349          if((i.eq.0).and.(j.eq.0)) then
     350            numpv=minloc(abs(xsec18*100.0-akz),dim=1)
     351          endif
     352          help=zsec4(nxfield*(ny-j-1)+i+1)
     353          if(i.lt.i180) then
     354            vvh(i180+i,j,numpv)=help
     355          else
     356            vvh(i-i180,j,numpv)=help
     357          endif
     358        endif
     359        if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
     360! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
     361          if((i.eq.0).and.(j.eq.0)) then
     362            numprh=minloc(abs(xsec18*100.0-akz),dim=1)
     363          endif
     364          help=zsec4(nxfield*(ny-j-1)+i+1)
     365          if(i.lt.i180) then
     366            qvh(i180+i,j,numprh,n)=help
     367          else
     368            qvh(i-i180,j,numprh,n)=help
     369          endif
     370        endif
     371        if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
     372! SURFACE PRESSURE
     373          help=zsec4(nxfield*(ny-j-1)+i+1)
     374          if(i.lt.i180) then
     375            ps(i180+i,j,1,n)=help
     376          else
     377            ps(i-i180,j,1,n)=help
     378          endif
     379        endif
     380        if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
     381! W VELOCITY
     382          if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
     383          help=zsec4(nxfield*(ny-j-1)+i+1)
     384          if(i.lt.i180) then
     385            wwh(i180+i,j,numpw)=help
     386          else
     387            wwh(i-i180,j,numpw)=help
     388          endif
     389        endif
     390        if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
     391! SNOW DEPTH
     392          help=zsec4(nxfield*(ny-j-1)+i+1)
     393          if(i.lt.i180) then
     394            sd(i180+i,j,1,n)=help
     395          else
     396            sd(i-i180,j,1,n)=help
     397          endif
     398        endif
     399        if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
     400! MEAN SEA LEVEL PRESSURE
     401          help=zsec4(nxfield*(ny-j-1)+i+1)
     402          if(i.lt.i180) then
     403            msl(i180+i,j,1,n)=help
     404          else
     405            msl(i-i180,j,1,n)=help
     406          endif
     407        endif
     408        if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
     409! TOTAL CLOUD COVER
     410          help=zsec4(nxfield*(ny-j-1)+i+1)
     411          if(i.lt.i180) then
     412            tcc(i180+i,j,1,n)=help
     413          else
     414            tcc(i-i180,j,1,n)=help
     415          endif
     416        endif
     417        if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
     418             (nint(xsec18).eq.10)) then
     419! 10 M U VELOCITY
     420          help=zsec4(nxfield*(ny-j-1)+i+1)
     421          if(i.lt.i180) then
     422            u10(i180+i,j,1,n)=help
     423          else
     424            u10(i-i180,j,1,n)=help
     425          endif
     426        endif
     427        if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
     428             (nint(xsec18).eq.10)) then
     429! 10 M V VELOCITY
     430          help=zsec4(nxfield*(ny-j-1)+i+1)
     431          if(i.lt.i180) then
     432            v10(i180+i,j,1,n)=help
     433          else
     434            v10(i-i180,j,1,n)=help
     435          endif
     436        endif
     437        if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
     438             (nint(xsec18).eq.2)) then
     439! 2 M TEMPERATURE
     440          help=zsec4(nxfield*(ny-j-1)+i+1)
     441          if(i.lt.i180) then
     442            tt2(i180+i,j,1,n)=help
     443          else
     444            tt2(i-i180,j,1,n)=help
     445          endif
     446        endif
     447        if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
     448             (nint(xsec18).eq.2)) then
     449! 2 M DEW POINT TEMPERATURE
     450          help=zsec4(nxfield*(ny-j-1)+i+1)
     451          if(i.lt.i180) then
     452            td2(i180+i,j,1,n)=help
     453          else
     454            td2(i-i180,j,1,n)=help
     455          endif
     456        endif
     457        if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
     458! LARGE SCALE PREC.
     459          help=zsec4(nxfield*(ny-j-1)+i+1)
     460          if(i.lt.i180) then
     461            lsprec(i180+i,j,1,n)=help
     462          else
     463            lsprec(i-i180,j,1,n)=help
     464          endif
     465        endif
     466        if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
     467! CONVECTIVE PREC.
     468          help=zsec4(nxfield*(ny-j-1)+i+1)
     469          if(i.lt.i180) then
     470            convprec(i180+i,j,1,n)=help
     471          else
     472            convprec(i-i180,j,1,n)=help
     473          endif
     474        endif
     475        if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
     476! TOPOGRAPHY
     477          help=zsec4(nxfield*(ny-j-1)+i+1)
     478          if(i.lt.i180) then
     479            oro(i180+i,j)=help
     480            excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     481          else
     482            oro(i-i180,j)=help
     483            excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
     484          endif
     485        endif
     486        if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
     487! LAND SEA MASK
     488          help=zsec4(nxfield*(ny-j-1)+i+1)
     489          if(i.lt.i180) then
     490            lsm(i180+i,j)=help
     491          else
     492            lsm(i-i180,j)=help
     493          endif
     494        endif
     495        if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
     496! MIXING HEIGHT
     497          help=zsec4(nxfield*(ny-j-1)+i+1)
     498          if(i.lt.i180) then
     499            hmix(i180+i,j,1,n)=help
     500          else
     501            hmix(i-i180,j,1,n)=help
     502          endif
     503        endif
     504        if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
     505             (nint(xsec18).eq.02)) then
     506! 2 M RELATIVE HUMIDITY
     507          help=zsec4(nxfield*(ny-j-1)+i+1)
     508          if(i.lt.i180) then
     509            qvh2(i180+i,j)=help
     510          else
     511            qvh2(i-i180,j)=help
     512          endif
     513        endif
     514        if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
     515! TEMPERATURE LOWEST SIGMA LEVEL
     516          help=zsec4(nxfield*(ny-j-1)+i+1)
     517          if(i.lt.i180) then
     518            tlev1(i180+i,j)=help
     519          else
     520            tlev1(i-i180,j)=help
     521          endif
     522        endif
     523        if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
     524! U VELOCITY LOWEST SIGMA LEVEL
     525          help=zsec4(nxfield*(ny-j-1)+i+1)
     526          if(i.lt.i180) then
     527            ulev1(i180+i,j)=help
     528          else
     529            ulev1(i-i180,j)=help
     530          endif
     531        endif
     532        if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
     533! V VELOCITY LOWEST SIGMA LEVEL
     534          help=zsec4(nxfield*(ny-j-1)+i+1)
     535          if(i.lt.i180) then
     536            vlev1(i180+i,j)=help
     537          else
     538            vlev1(i-i180,j)=help
     539          endif
     540        endif
    541541! 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
    559 
    560 
     542        if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg]
     543          if((i.eq.0).and.(j.eq.0)) then
     544            numpclwch=minloc(abs(xsec18*100.0-akz),dim=1)
     545          endif
     546          help=zsec4(nxfield*(ny-j-1)+i+1)
     547          if(i.lt.i180) then
     548            clwch(i180+i,j,numpclwch,n)=help
     549          else
     550            clwch(i-i180,j,numpclwch,n)=help
     551          endif
     552          readclouds=.true.
     553          sumclouds=.true.
     554        endif
     555
     556
     557      end do
    561558    end do
    562   end do
    563559
    564560  endif
    565561
    566562  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
    567   ! NCEP ISOBARIC LEVELS
     563! NCEP ISOBARIC LEVELS
    568564    iumax=iumax+1
    569565  endif
    … …  
    571567  call grib_release(igrib)
    572568  goto 10                      !! READ NEXT LEVEL OR PARAMETER
    573   !
    574   ! CLOSING OF INPUT DATA FILE
    575   !
    576 
    577   !HSO close grib file
    578 50   continue
     569!
     570! CLOSING OF INPUT DATA FILE
     571!
     572
     573!HSO close grib file
     57450 continue
    579575  call grib_close_file(ifile)
    580576
    581   ! SENS. HEAT FLUX
     577! SENS. HEAT FLUX
    582578  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
    583579  hflswitch=.false.    ! Heat flux not available
    584   ! SOLAR RADIATIVE FLUXES
     580! SOLAR RADIATIVE FLUXES
    585581  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
    586   ! EW SURFACE STRESS
     582! EW SURFACE STRESS
    587583  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
    588   ! NS SURFACE STRESS
     584! NS SURFACE STRESS
    589585  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
    590586  strswitch=.false.    ! stress not available
    591587
    592   ! CONVERT TP TO LSP (GRIB2 only)
     588! CONVERT TP TO LSP (GRIB2 only)
    593589  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
    613   endif
    614   !HSO end edits
    615 
    616 
    617   ! TRANSFORM RH TO SPECIFIC HUMIDITY
     590    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) )
     591  endif
     592!HSO end edits
     593
     594
     595! TRANSFORM RH TO SPECIFIC HUMIDITY
    618596
    619597  do j=0,ny-1
    … …  
    622600        help=qvh(i,j,k,n)
    623601        temp=tth(i,j,k,n)
     602        if (temp .eq. 0.0) then
     603          write (*, *) i, j, k, n
     604          temp = 273.0
     605        endif
    624606        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
    625607        elev=ew(temp)*help/100.0
    … …  
    629611  end do
    630612
    631   ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
    632   ! USING BOLTON'S (1980) FORMULA
    633   ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
     613! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
     614! USING BOLTON'S (1980) FORMULA
     615! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
    634616
    635617  do j=0,ny-1
    636618    do i=0,nxfield-1
    637         help=qvh2(i,j)
    638         temp=tt2(i,j,1,n)
    639         elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
    640         td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
    641         if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
     619      help=qvh2(i,j)
     620      temp=tt2(i,j,1,n)
     621      if (temp .eq. 0.0) then
     622        write (*, *) i, j, n
     623        temp = 273.0
     624      endif
     625      elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
     626      td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
     627      if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
    642628    end do
    643629  end do
    … …  
    653639
    654640
    655   ! For global fields, assign the leftmost data column also to the rightmost
    656   ! data column; if required, shift whole grid by nxshift grid points
    657   !*************************************************************************
     641! For global fields, assign the leftmost data column also to the rightmost
     642! data column; if required, shift whole grid by nxshift grid points
     643!*************************************************************************
    658644
    659645  if (xglobal) then
    … …  
    691677  do i=0,nxmin1
    692678    do j=0,nymin1
    693   ! Convert precip. from mm/s -> mm/hour
     679! Convert precip. from mm/s -> mm/hour
    694680      convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
    695681      lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
    … …  
    699685
    700686  if ((.not.hflswitch).or.(.not.strswitch)) then
    701   !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
    702   !    +  wfname(indj)
    703 
    704   ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
    705   !***************************************************************************
     687!  write(*,*) 'WARNING: No flux data contained in GRIB file ',
     688!    +  wfname(indj)
     689
     690! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
     691!***************************************************************************
    706692
    707693    do i=0,nxmin1
    … …  
    723709
    724710  return
    725 888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     711888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    726712  write(*,*) ' #### ',wfname(indj),'                    #### '
    727713  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
    728714  stop 'Execution terminated'
    729 999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
     715999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
    730716  write(*,*) ' #### ',wfname(indj),'                    #### '
    731717  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '

  • 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