Changeset bb579a9 in flexpart.git


Ignore:
Timestamp:
Sep 13, 2017, 11:55:16 AM (7 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
08a38b5
Parents:
aa8c34a (diff), dd6a4aa (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 'ctbto' into dev

Files:
8 added
7 deleted
28 edited
1 moved

Legend:

Unmodified
Added
Removed
  • src/FLEXPART.f90

    rc8fc724 rd8eed02  
    3333  !                                                                            *
    3434  !*****************************************************************************
     35  ! Changes:                                                                   *
     36  !   Unified ECMWF and GFS builds                                             *
     37  !   Marian Harustak, 12.5.2017                                               *
     38  !     - Added detection of metdata format using gributils routines           *
     39  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
     40  !       detected metdata format                                              *
     41  !     - Passed metdata format down to timemanager                            *
     42  !*****************************************************************************
    3543  !                                                                            *
    3644  ! Variables:                                                                 *
     
    4654  use netcdf_output_mod, only: writeheader_netcdf
    4755  use random_mod, only: gasdev1
     56  use class_gribfile
    4857
    4958  implicit none
     
    5261  integer :: idummy = -320
    5362  character(len=256) :: inline_options  !pathfile, flexversion, arg2
     63  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
     64  integer :: detectformat
     65
    5466
    5567
     
    7082  ! FLEXPART version string
    7183  flexversion_major = '10' ! Major version number, also used for species file names
    72   flexversion='Version '//trim(flexversion_major)//'.1beta (2016-11-02)'
     84  flexversion='Version '//trim(flexversion_major)//'.2beta (2017-08-01)'
    7385  verbosity=0
    7486
     
    172184  call readavailable
    173185
     186  ! Detect metdata format
     187  !**********************
     188
     189  metdata_format = detectformat()
     190
     191  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     192    print *,'ECMWF metdata detected'
     193  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
     194    print *,'NCEP metdata detected'
     195  else
     196    print *,'Unknown metdata format'
     197    return
     198  endif
     199
     200
     201
    174202  ! If nested wind fields are used, allocate arrays
    175203  !************************************************
     
    188216  endif
    189217
    190   call gridcheck
     218  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     219    call gridcheck_ecmwf
     220  else
     221    call gridcheck_gfs
     222  end if
    191223
    192224  if (verbosity.gt.1) then   
     
    411443  endif
    412444
    413   call timemanager
     445  call timemanager(metdata_format)
    414446
    415447! NIK 16.02.2005
  • src/FLEXPART_MPI.f90

    r5184a7c rd8eed02  
    3333  !                                                                            *
    3434  !*****************************************************************************
     35  ! Changes:                                                                   *
     36  !   Unified ECMWF and GFS builds                                             *
     37  !   Marian Harustak, 12.5.2017                                               *
     38  !     - Added detection of metdata format using gributils routines           *
     39  !     - Distinguished calls to ecmwf/gfs gridcheck versions based on         *
     40  !       detected metdata format                                              *
     41  !     - Passed metdata format down to timemanager                            *
     42  !*****************************************************************************
    3543  !                                                                            *
    3644  ! Variables:                                                                 *
     
    4755  use netcdf_output_mod, only: writeheader_netcdf
    4856  use random_mod, only: gasdev1
     57  use class_gribfile
    4958
    5059  implicit none
     
    5362  integer :: idummy = -320
    5463  character(len=256) :: inline_options  !pathfile, flexversion, arg2
     64  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
     65  integer :: detectformat
     66
    5567
    5668
     
    8092  ! FLEXPART version string
    8193  flexversion_major = '10' ! Major version number, also used for species file names
    82   flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
     94  flexversion='Ver. '//trim(flexversion_major)//'.2beta MPI (2017-08-01)'
    8395  verbosity=0
    8496
     
    197209  call readavailable
    198210
     211  ! Detect metdata format
     212  !**********************
     213
     214  metdata_format = detectformat()
     215
     216  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     217    print *,'ECMWF metdata detected'
     218  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
     219    print *,'NCEP metdata detected'
     220  else
     221    print *,'Unknown metdata format'
     222    return
     223  endif
     224
     225
     226
    199227  ! If nested wind fields are used, allocate arrays
    200228  !************************************************
     
    213241  endif
    214242
    215   call gridcheck
     243  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     244    call gridcheck_ecmwf
     245  else
     246    call gridcheck_gfs
     247  end if
     248
    216249
    217250  if (verbosity.gt.1 .and. lroot) then   
     
    456489
    457490
    458   call timemanager
     491  call timemanager(metdata_format)
    459492
    460493
  • src/calcmatrix.f90

    re200b7a rd8eed02  
    2020!**********************************************************************
    2121
    22 subroutine calcmatrix(lconv,delt,cbmf)
     22subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
    2323  !                        o    i    o
    2424  !*****************************************************************************
     
    3030  !  Petra Seibert, Bernd C. Krueger, 2000-2001                                *
    3131  !                                                                            *
     32  !*****************************************************************************
     33  ! Changes:                                                                   *
    3234  !  changed by C. Forster, November 2003 - February 2004                      *
    3335  !  array fmassfrac(nconvlevmax,nconvlevmax) represents                       *
    3436  !  the convective redistribution matrix for the particles                    *
    3537  !                                                                            *
     38  !   Unified ECMWF and GFS builds                                             *
     39  !   Marian Harustak, 12.5.2017                                               *
     40  !     - Merged calcmatrix and calcmatrix_gfs into one routine using if-then  *
     41  !       for meteo-type dependent code                                        *
     42  !*****************************************************************************
     43  !                                                                            *
    3644  !  lconv        indicates whether there is convection in this cell, or not   *
    3745  !  delt         time step for convection [s]                                 *
    3846  !  cbmf         cloud base mass flux                                         *
     47  !  metdata_format format of metdata (ecmwf/gfs)                              *
    3948  !                                                                            *
    4049  !*****************************************************************************
     
    4352  use com_mod
    4453  use conv_mod
     54  use class_gribfile
    4555
    4656  implicit none
    4757
    4858  real :: rlevmass,summe
     59  integer :: metdata_format
    4960
    5061  integer :: iflag, k, kk, kuvz
     
    7788  do kuvz = 2,nuvz
    7889    k = kuvz-1
     90    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    7991    pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
    8092    phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
     93    else
     94      phconv(kuvz) =  0.5*(pconv(kuvz)+pconv(k))
     95    endif
    8196    dpr(k) = phconv(k) - phconv(kuvz)
    8297    qsconv(k) = f_qvsat( pconv(k), tconv(k) )
     
    85100    do kk=1,nconvlev
    86101      fmassfrac(k,kk)=0.
    87     enddo
    88   enddo
     102    end do
     103  end do
    89104
    90105
  • src/calcpar.f90

    re200b7a r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine calcpar(n,uuh,vvh,pvh)
     22subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
    2323  !                   i  i   i   o
    2424  !*****************************************************************************
     
    3737  !     new variables in call to richardson                                    *
    3838  !                                                                            *
    39   !*****************************************************************************
    40   !  Changes, Bernd C. Krueger, Feb. 2001:
    41   !   Variables tth and qvh (on eta coordinates) in common block
     39  !   CHANGE 17/11/2005 Caroline Forster NCEP GFS version                      *
     40  !                                                                            *
     41  !   Changes, Bernd C. Krueger, Feb. 2001:                                    *
     42  !    Variables tth and qvh (on eta coordinates) in common block              *
     43  !                                                                            *
     44  !   Unified ECMWF and GFS builds                                             *
     45  !   Marian Harustak, 12.5.2017                                               *
     46  !     - Merged calcpar and calcpar_gfs into one routine using if-then        *
     47  !       for meteo-type dependent code                                        *
     48  !*****************************************************************************
     49
    4250  !*****************************************************************************
    4351  !                                                                            *
    4452  ! Variables:                                                                 *
    4553  ! n                  temporal index for meteorological fields (1 to 3)       *
     54  ! uuh                                                                        *
     55  ! vvh                                                                        *
     56  ! pvh                                                                        *
     57  ! metdata_format     format of metdata (ecmwf/gfs)                           *
    4658  !                                                                            *
    4759  ! Constants:                                                                 *
     
    5668  use par_mod
    5769  use com_mod
     70  use class_gribfile
    5871
    5972  implicit none
    6073
    61   integer :: n,ix,jy,i,kz,lz,kzmin
     74  integer :: metdata_format
     75  integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start
    6276  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
    6377  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
    64   real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
     78  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
    6579  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    6680  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
     
    111125  !***********************************************
    112126
     127      if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
     128        ! NCEP version: find first level above ground
     129        llev = 0
     130        do i=1,nuvz
     131          if (ps(ix,jy,1,n).lt.akz(i)) llev=i
     132        end do
     133        llev = llev+1
     134        if (llev.gt.nuvz) llev = nuvz-1
     135        ! NCEP version
     136
     137        ! calculate inverse Obukhov length scale with tth(llev)
    113138      ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
    114            tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
     139            tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
     140      else
     141        llev=0
     142        ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
     143            tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
     144      end if
     145
    115146      if (ol.ne.0.) then
    116147        oli(ix,jy,1,n)=1./ol
     
    130161      end do
    131162
     163      if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
     164        ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
    132165      call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
    133166           ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
    134            td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
     167             td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
     168      else
     169        call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
     170             ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
     171             td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,metdata_format)
     172      end if
    135173
    136174      if(lsubgrid.eq.1) then
     
    173211  !******************************************************
    174212
    175   ! 1) Calculate altitudes of ECMWF model levels
    176   !*********************************************
     213  ! 1) Calculate altitudes of model levels
     214  !***************************************
    177215
    178216      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
     
    180218      pold=ps(ix,jy,1,n)
    181219      zold=0.
    182       do kz=2,nuvz
     220      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     221        loop_start=2
     222      else
     223        loop_start=llev
     224      end if
     225      do kz=loop_start,nuvz
    183226        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)  ! pressure on model layers
    184227        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
     
    199242  !************************************************************************
    200243
    201       do kz=1,nuvz
     244      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     245        loop_start=1
     246      else
     247        loop_start=llev
     248      end if
     249
     250      do kz=loop_start,nuvz
    202251        if (zlev(kz).ge.altmin) then
    203252          kzmin=kz
  • src/calcpar_nests.f90

    rdb712a8 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
     22subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
    2323  !                         i  i    i    o
    2424  !*****************************************************************************
     
    3939  !                                                                            *
    4040  !*****************************************************************************
    41   !  Changes, Bernd C. Krueger, Feb. 2001:
    42   !   Variables tth and qvh (on eta coordinates) in common block
     41  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
     42  !   Variables tth and qvh (on eta coordinates) in common block               *
     43  !                                                                            *
     44  !   Unified ECMWF and GFS builds                                             *
     45  !   Marian Harustak, 12.5.2017                                               *
     46  !     - Added passing of metdata_format as it was needed by called routines  *
    4347  !*****************************************************************************
    4448  !                                                                            *
    4549  ! Variables:                                                                 *
    4650  ! n                  temporal index for meteorological fields (1 to 3)       *
     51  ! metdata_format     format of metdata (ecmwf/gfs)                           *
    4752  !                                                                            *
    4853  ! Constants:                                                                 *
     
    6065  implicit none
    6166
     67  integer :: metdata_format
    6268  integer :: n,ix,jy,i,l,kz,lz,kzmin
    63   real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
     69  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
    6470  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
    6571  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
     
    110116      ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
    111117           td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
    112            sshfn(ix,jy,1,n,l),akm,bkm)
     118           sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,metdata_format)
    113119      if (ol.ne.0.) then
    114120        olin(ix,jy,1,n,l)=1./ol
     
    131137           qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
    132138           tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
    133            wstarn(ix,jy,1,n,l),hmixplus)
     139           wstarn(ix,jy,1,n,l),hmixplus,metdata_format)
    134140
    135141      if(lsubgrid.eq.1) then
  • src/convmix.f90

    r8a65cb0 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine convmix(itime)
     22subroutine convmix(itime,metdata_format)
    2323  !                     i
    2424  !**************************************************************
     
    3131  !CHANGES by A. Stohl:
    3232  !  various run-time optimizations - February 2005
     33  ! CHANGES by C. Forster, November 2005, NCEP GFS version
     34  !      in the ECMWF version convection is calculated on the
     35  !      original eta-levels
     36  !      in the GFS version convection is calculated on the
     37  !      FLEXPART levels
     38  !
     39  !   Unified ECMWF and GFS builds                                             
     40  !   Marian Harustak, 12.5.2017                                             
     41  !     - Merged convmix and convmix_gfs into one routine using if-then           
     42  !       for meteo-type dependent code                                       
    3343  !**************************************************************
    3444
     
    3747  use com_mod
    3848  use conv_mod
     49  use class_gribfile
    3950
    4051  implicit none
     
    4455  integer :: jy, kpart, ktop, ngrid,kz
    4556  integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
     57  integer :: metdata_format
    4658
    4759  ! itime [s]                 current time
     
    104116
    105117    ngrid=0
     118    if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    106119    do j=numbnests,1,-1
    107120      if ( x.gt.xln(j)+eps .and. x.lt.xrn(j)-eps .and. &
     
    111124      endif
    112125    end do
     126    else
     127      do j=numbnests,1,-1
     128        if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. &
     129             y.gt.yln(j) .and. y.lt.yrn(j) ) then
     130          ngrid=j
     131          goto 23
     132        endif
     133      end do
     134    endif
    113135 23   continue
    114136
     
    167189      td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt
    168190!!$      do kz=1,nconvlev+1      !old
     191      if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    169192        do kz=1,nuvz-1           !bugfix
    170193        tconv(kz)=(tth(ix,jy,kz+1,mind1)*dt2+ &
     
    173196             qvh(ix,jy,kz+1,mind2)*dt1)*dtt
    174197      end do
     198      else
     199        do kz=1,nuvz-1           !bugfix
     200          pconv(kz)=(pplev(ix,jy,kz,mind1)*dt2+ &
     201              pplev(ix,jy,kz,mind2)*dt1)*dtt
     202          tconv(kz)=(tt(ix,jy,kz,mind1)*dt2+ &
     203              tt(ix,jy,kz,mind2)*dt1)*dtt
     204          qconv(kz)=(qv(ix,jy,kz,mind1)*dt2+ &
     205              qv(ix,jy,kz,mind2)*dt1)*dtt
     206        end do
     207      end if
    175208
    176209  ! Calculate translocation matrix
    177       call calcmatrix(lconv,delt,cbaseflux(ix,jy))
     210      call calcmatrix(lconv,delt,cbaseflux(ix,jy),metdata_format)
    178211      igrold = igr
    179212      ktop = 0
     
    252285  ! calculate translocation matrix
    253286  !*******************************
    254         call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest))
     287        call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest),metdata_format)
    255288        igrold = igr
    256289        ktop = 0
  • src/getfields.f90

    rd1a8707 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine getfields(itime,nstop)
     22subroutine getfields(itime,nstop,metdata_format)
    2323!                       i     o
    2424!*****************************************************************************
     
    3838!                                                                            *
    3939!*****************************************************************************
    40 !  Changes, Bernd C. Krueger, Feb. 2001:
    41 !   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
    42 !   Function of nstop extended.
     40!  Changes, Bernd C. Krueger, Feb. 2001:                                     *
     41!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.        *
     42!   Function of nstop extended.                                              *
     43!                                                                            *
     44!   Unified ECMWF and GFS builds                                             *
     45!   Marian Harustak, 12.5.2017                                               *
     46!     - Added passing of metdata_format as it was needed by called routines  *
    4347!*****************************************************************************
    4448!                                                                            *
     
    5862! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
    5963! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
     64! metdata_format     format of metdata (ecmwf/gfs)                           *
    6065!                                                                            *
    6166! Constants:                                                                 *
     
    6772  use par_mod
    6873  use com_mod
     74  use class_gribfile
    6975
    7076  implicit none
    7177
    7278  integer :: indj,itime,nstop,memaux
     79  integer :: metdata_format
    7380
    7481  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     
    125132    do indj=indmin,numbwf-1
    126133      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
    127         call readwind(indj+1,memind(2),uuh,vvh,wwh)
     134        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     135          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
     136        else
     137          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
     138        end if
    128139        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
    129         call calcpar(memind(2),uuh,vvh,pvh)
    130         call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
    131         call verttransform(memind(2),uuh,vvh,wwh,pvh)
     140        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
     141        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
     142        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     143          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
     144        else
     145          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
     146        end if
    132147        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
    133148        memtime(2)=wftime(indj+1)
     
    152167           (ldirect*wftime(indj+1).gt.ldirect*itime)) then
    153168        memind(1)=1
    154         call readwind(indj,memind(1),uuh,vvh,wwh)
     169        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     170          call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
     171        else
     172          call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
     173        end if
    155174        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
    156         call calcpar(memind(1),uuh,vvh,pvh)
    157         call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
    158         call verttransform(memind(1),uuh,vvh,wwh,pvh)
     175        call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
     176        call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
     177        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     178          call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
     179        else
     180          call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
     181        end if
    159182        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
    160183        memtime(1)=wftime(indj)
    161184        memind(2)=2
    162         call readwind(indj+1,memind(2),uuh,vvh,wwh)
     185        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     186          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
     187        else
     188          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
     189        end if
    163190        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
    164         call calcpar(memind(2),uuh,vvh,pvh)
    165         call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
    166         call verttransform(memind(2),uuh,vvh,wwh,pvh)
     191        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
     192        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
     193        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     194          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
     195        else
     196          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
     197        end if
    167198        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
    168199        memtime(2)=wftime(indj+1)
  • src/getfields_mpi.f90

    r78e62dc r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine getfields(itime,nstop,memstat)
     22subroutine getfields(itime,nstop,memstat,metdata_format)
    2323!                       i     o       o
    2424!*****************************************************************************
     
    5858!    memstat=0:      no new fields to be read
    5959!         
     60!   Unified ECMWF and GFS builds                                             
     61!   Marian Harustak, 12.5.2017                                               
     62!     - Added passing of metdata_format as it was needed by called routines 
     63!         
    6064!*****************************************************************************
    6165!                                                                            *
     
    7781! tt(0:nxmax,0:nymax,nuvzmax,numwfmem)   temperature [K]                     *
    7882! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
     83! metdata_format     format of metdata (ecmwf/gfs)                           *
    7984!                                                                            *
    8085! Constants:                                                                 *
     
    8792  use com_mod
    8893  use mpi_mod, only: lmpreader,lmp_use_reader,lmp_sync
     94  use class_gribfile
    8995
    9096  implicit none
    9197
     98  integer :: metdata_format
    9299  integer :: indj,itime,nstop,memaux,mindread
    93100! mindread: index of where to read 3rd field
     
    204211      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
    205212        if (lmpreader.or..not. lmp_use_reader) then
    206           call readwind(indj+1,mindread,uuh,vvh,wwh)
     213          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     214            call readwind_ecmwf(indj+1,mindread,uuh,vvh,wwh)
     215          else
     216            call readwind_gfs(indj+1,mindread,uuh,vvh,wwh)
     217          end if
    207218          call readwind_nests(indj+1,mindread,uuhn,vvhn,wwhn)
    208           call calcpar(mindread,uuh,vvh,pvh)
    209           call calcpar_nests(mindread,uuhn,vvhn,pvhn)
    210           call verttransform(mindread,uuh,vvh,wwh,pvh)
     219          call calcpar(mindread,uuh,vvh,pvh,metdata_format)
     220          call calcpar_nests(mindread,uuhn,vvhn,pvhn,metdata_format)
     221          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     222            call verttransform_ecmwf(mindread,uuh,vvh,wwh,pvh)
     223          else
     224            call verttransform_gfs(mindread,uuh,vvh,wwh,pvh)
     225          end if
    211226          call verttransform_nests(mindread,uuhn,vvhn,wwhn,pvhn)
    212227        end if
     
    230245        memind(1)=1
    231246        if (lmpreader.or..not.lmp_use_reader) then
    232           call readwind(indj,memind(1),uuh,vvh,wwh)
     247          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     248            call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
     249          else
     250            call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
     251          end if
    233252          call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
    234           call calcpar(memind(1),uuh,vvh,pvh)
    235           call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
    236           call verttransform(memind(1),uuh,vvh,wwh,pvh)
     253          call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
     254          call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
     255          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     256            call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
     257          else
     258            call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
     259          end if
    237260          call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
    238261        end if
     
    240263        memind(2)=2
    241264        if (lmpreader.or..not.lmp_use_reader) then
    242           call readwind(indj+1,memind(2),uuh,vvh,wwh)
     265          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     266            call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
     267          else
     268            call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
     269          end if
    243270          call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
    244           call calcpar(memind(2),uuh,vvh,pvh)
    245           call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
    246           call verttransform(memind(2),uuh,vvh,wwh,pvh)
     271          call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
     272          call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
     273          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     274            call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
     275          else
     276            call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
     277          end if
    247278          call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
    248279        end if
  • src/gridcheck_ecmwf.f90

    rdb712a8 r61e07ba  
    2020!**********************************************************************
    2121
    22 subroutine gridcheck
     22subroutine gridcheck_ecmwf
    2323
    2424  !**********************************************************************
     
    3737  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    3838  !                                 ECMWF grib_api                      *
     39  !                                                                     *
     40  !   Unified ECMWF and GFS builds                                      *
     41  !   Marian Harustak, 12.5.2017                                        *
     42  !     - Renamed from gridcheck to gridcheck_ecmwf                     *
    3943  !                                                                     *
    4044  !**********************************************************************
     
    117121  ! OPENING OF DATA FILE (GRIB CODE)
    118122  !
    119 5   call grib_open_file(ifile,path(3)(1:length(3)) &
    120          //trim(wfname(ifn)),'r',iret)
     1235 call grib_open_file(ifile,path(3)(1:length(3)) &
     124       //trim(wfname(ifn)),'r',iret)
    121125  if (iret.ne.GRIB_SUCCESS) then
    122126    goto 999   ! ERROR DETECTED
     
    127131  gotGrid=0
    128132  ifield=0
    129 10   ifield=ifield+1
     13310 ifield=ifield+1
    130134
    131135  !
     
    145149  if (gribVer.eq.1) then ! GRIB Edition 1
    146150
    147   !print*,'GRiB Edition 1'
    148   !read the grib2 identifiers
    149   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    150   call grib_check(iret,gribFunction,gribErrorMsg)
    151   call grib_get_int(igrib,'level',isec1(8),iret)
    152   call grib_check(iret,gribFunction,gribErrorMsg)
    153 
    154   !change code for etadot to code for omega
    155   if (isec1(6).eq.77) then
    156     isec1(6)=135
    157   endif
    158 
    159   !print*,isec1(6),isec1(8)
     151    !print*,'GRiB Edition 1'
     152    !read the grib2 identifiers
     153    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
     154    call grib_check(iret,gribFunction,gribErrorMsg)
     155    call grib_get_int(igrib,'level',isec1(8),iret)
     156    call grib_check(iret,gribFunction,gribErrorMsg)
     157
     158    !change code for etadot to code for omega
     159    if (isec1(6).eq.77) then
     160      isec1(6)=135
     161    endif
     162
     163    !print*,isec1(6),isec1(8)
    160164
    161165  else
    162166
    163   !print*,'GRiB Edition 2'
    164   !read the grib2 identifiers
    165   call grib_get_int(igrib,'discipline',discipl,iret)
    166   call grib_check(iret,gribFunction,gribErrorMsg)
    167   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    168   call grib_check(iret,gribFunction,gribErrorMsg)
    169   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    170   call grib_check(iret,gribFunction,gribErrorMsg)
    171   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    172   call grib_check(iret,gribFunction,gribErrorMsg)
    173   call grib_get_int(igrib,'level',valSurf,iret)
    174   call grib_check(iret,gribFunction,gribErrorMsg)
    175   call grib_get_int(igrib,'paramId',parId,iret)
    176   call grib_check(iret,gribFunction,gribErrorMsg)
    177 
    178   !print*,discipl,parCat,parNum,typSurf,valSurf
    179 
    180   !convert to grib1 identifiers
    181   isec1(6)=-1
    182   isec1(7)=-1
    183   isec1(8)=-1
    184   isec1(8)=valSurf     ! level
    185   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    186     isec1(6)=130         ! indicatorOfParameter
    187   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
    188     isec1(6)=131         ! indicatorOfParameter
    189   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
    190     isec1(6)=132         ! indicatorOfParameter
    191   elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    192     isec1(6)=133         ! indicatorOfParameter
    193 !ZHG FOR CLOUDS FROM GRIB
    194   elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
    195     isec1(6)=246         ! indicatorOfParameter
    196   elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
    197     isec1(6)=247         ! indicatorOfParameter
    198 !ZHG end
    199 ! ESO qc(=clwc+ciwc)
    200   elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
    201     isec1(6)=201031      ! indicatorOfParameter
    202   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    203     isec1(6)=134         ! indicatorOfParameter
    204   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    205     isec1(6)=135         ! indicatorOfParameter
    206   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    207     isec1(6)=135         ! indicatorOfParameter
    208   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    209     isec1(6)=151         ! indicatorOfParameter
    210   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
    211     isec1(6)=165         ! indicatorOfParameter
    212   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
    213     isec1(6)=166         ! indicatorOfParameter
    214   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
    215     isec1(6)=167         ! indicatorOfParameter
    216   elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
    217     isec1(6)=168         ! indicatorOfParameter
    218   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    219     isec1(6)=141         ! indicatorOfParameter
    220   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    221     isec1(6)=164         ! indicatorOfParameter
    222   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    223     isec1(6)=142         ! indicatorOfParameter
    224   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    225     isec1(6)=143         ! indicatorOfParameter
    226   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    227     isec1(6)=146         ! indicatorOfParameter
    228   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    229     isec1(6)=176         ! indicatorOfParameter
    230   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    231     isec1(6)=180         ! indicatorOfParameter
    232   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    233     isec1(6)=181         ! indicatorOfParameter
    234   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    235     isec1(6)=129         ! indicatorOfParameter
    236   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    237     isec1(6)=160         ! indicatorOfParameter
    238   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    239        (typSurf.eq.1)) then ! LSM
    240     isec1(6)=172         ! indicatorOfParameter
    241   else
    242     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
    243          parCat,parNum,typSurf
    244   endif
    245   if(parId .ne. isec1(6) .and. parId .ne. 77) then
    246     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
    247 !    stop
    248   endif
    249 
    250   endif
     167    !print*,'GRiB Edition 2'
     168    !read the grib2 identifiers
     169    call grib_get_int(igrib,'discipline',discipl,iret)
     170    call grib_check(iret,gribFunction,gribErrorMsg)
     171    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     172    call grib_check(iret,gribFunction,gribErrorMsg)
     173    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     174    call grib_check(iret,gribFunction,gribErrorMsg)
     175    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     176    call grib_check(iret,gribFunction,gribErrorMsg)
     177    call grib_get_int(igrib,'level',valSurf,iret)
     178    call grib_check(iret,gribFunction,gribErrorMsg)
     179    call grib_get_int(igrib,'paramId',parId,iret)
     180    call grib_check(iret,gribFunction,gribErrorMsg)
     181
     182    !print*,discipl,parCat,parNum,typSurf,valSurf
     183
     184    !convert to grib1 identifiers
     185    isec1(6)=-1
     186    isec1(7)=-1
     187    isec1(8)=-1
     188    isec1(8)=valSurf     ! level
     189    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     190      isec1(6)=130         ! indicatorOfParameter
     191    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
     192      isec1(6)=131         ! indicatorOfParameter
     193    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
     194      isec1(6)=132         ! indicatorOfParameter
     195    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     196      isec1(6)=133         ! indicatorOfParameter
     197      !ZHG FOR CLOUDS FROM GRIB
     198    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     199      isec1(6)=246         ! indicatorOfParameter
     200    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     201      isec1(6)=247         ! indicatorOfParameter
     202      !ZHG end
     203      ! ESO qc(=clwc+ciwc)
     204    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
     205      isec1(6)=201031      ! indicatorOfParameter
     206    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     207      isec1(6)=134         ! indicatorOfParameter
     208    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     209      isec1(6)=135         ! indicatorOfParameter
     210    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
     211      isec1(6)=135         ! indicatorOfParameter
     212    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     213      isec1(6)=151         ! indicatorOfParameter
     214    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
     215      isec1(6)=165         ! indicatorOfParameter
     216    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
     217      isec1(6)=166         ! indicatorOfParameter
     218    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
     219      isec1(6)=167         ! indicatorOfParameter
     220    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
     221      isec1(6)=168         ! indicatorOfParameter
     222    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     223      isec1(6)=141         ! indicatorOfParameter
     224    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     225      isec1(6)=164         ! indicatorOfParameter
     226    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     227      isec1(6)=142         ! indicatorOfParameter
     228    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     229      isec1(6)=143         ! indicatorOfParameter
     230    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     231      isec1(6)=146         ! indicatorOfParameter
     232    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     233      isec1(6)=176         ! indicatorOfParameter
     234    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     235      isec1(6)=180         ! indicatorOfParameter
     236    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     237      isec1(6)=181         ! indicatorOfParameter
     238    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     239      isec1(6)=129         ! indicatorOfParameter
     240    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     241      isec1(6)=160         ! indicatorOfParameter
     242    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
     243         (typSurf.eq.1)) then ! LSM
     244      isec1(6)=172         ! indicatorOfParameter
     245    else
     246      print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     247           parCat,parNum,typSurf
     248    endif
     249    if(parId .ne. isec1(6) .and. parId .ne. 77) then
     250      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     251      !    stop
     252    endif
     253
     254  endif
     255
     256  CALL grib_get_int(igrib,'numberOfPointsAlongAParallel', &
     257       isec2(2),iret)
     258  ! nx=isec2(2)
     259  ! WRITE(*,*) nx,nxmax
     260  IF (isec2(2).GT.nxmax) THEN
     261    WRITE(*,*) 'FLEXPART error: Too many grid points in x direction.'
     262    WRITE(*,*) 'Reduce resolution of wind fields.'
     263    WRITE(*,*) 'Or change parameter settings in file ecmwf_mod.'
     264    WRITE(*,*) nx,nxmax
     265!    STOP
     266  ENDIF
    251267
    252268  !get the size and data of the values array
     
    258274  if (ifield.eq.1) then
    259275
    260   !HSO  get the required fields from section 2 in a gribex compatible manner
     276    !HSO  get the required fields from section 2 in a gribex compatible manner
    261277    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    262278         isec2(2),iret)
     
    272288    call grib_check(iret,gribFunction,gribErrorMsg)
    273289
    274   !  get the size and data of the vertical coordinate array
     290    !  get the size and data of the vertical coordinate array
    275291    call grib_get_real4_array(igrib,'pv',zsec2,iret)
    276292    call grib_check(iret,gribFunction,gribErrorMsg)
     
    279295    ny=isec2(3)
    280296    nlev_ec=isec2(12)/2-1
    281    endif
     297  endif
    282298
    283299  !HSO  get the second part of the grid dimensions only from GRiB1 messages
     
    308324    dyconst=180./(dy*r_earth*pi)
    309325    gotGrid=1
    310   ! Check whether fields are global
    311   ! If they contain the poles, specify polar stereographic map
    312   ! projections using the stlmbr- and stcm2p-calls
    313   !***********************************************************
     326    ! Check whether fields are global
     327    ! If they contain the poles, specify polar stereographic map
     328    ! projections using the stlmbr- and stcm2p-calls
     329    !***********************************************************
    314330
    315331    xauxa=abs(xaux2+dx-360.-xaux1)
     
    332348    if (xglobal.and.xauxa.lt.0.001) then
    333349      sglobal=.true.               ! field contains south pole
    334   ! Enhance the map scale by factor 3 (*2=6) compared to north-south
    335   ! map scale
     350      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
     351      ! map scale
    336352      sizesouth=6.*(switchsouth+90.)/dy
    337353      call stlmbr(southpolemap,-90.,0.)
     
    346362    if (xglobal.and.xauxa.lt.0.001) then
    347363      nglobal=.true.               ! field contains north pole
    348   ! Enhance the map scale by factor 3 (*2=6) compared to north-south
    349   ! map scale
     364      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
     365      ! map scale
    350366      sizenorth=6.*(90.-switchnorth)/dy
    351367      call stlmbr(northpolemap,90.,0.)
     
    363379
    364380  if (nx.gt.nxmax) then
    365    write(*,*) 'FLEXPART error: Too many grid points in x direction.'
     381    write(*,*) 'FLEXPART error: Too many grid points in x direction.'
    366382    write(*,*) 'Reduce resolution of wind fields.'
    367383    write(*,*) 'Or change parameter settings in file par_mod.'
     
    371387
    372388  if (ny.gt.nymax) then
    373    write(*,*) 'FLEXPART error: Too many grid points in y direction.'
     389    write(*,*) 'FLEXPART error: Too many grid points in y direction.'
    374390    write(*,*) 'Reduce resolution of wind fields.'
    375391    write(*,*) 'Or change parameter settings in file par_mod.'
     
    410426  !
    411427
    412 30   call grib_close_file(ifile)
     42830 call grib_close_file(ifile)
    413429
    414430  !error message if no fields found with correct first longitude in it
     
    469485
    470486  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
    471                         ! by trajectory model
     487  ! by trajectory model
    472488  !do 8940 i=1,244
    473489  !   write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip
     
    482498    k=nlev_ec+1+numskip+i
    483499    akm(nwz-i+1)=zsec2(j)
    484   !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
     500    !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
    485501    bkm(nwz-i+1)=zsec2(k)
    486502  end do
     
    498514  bkz(1)=1.0
    499515  do i=1,nuvz
    500      akz(i+1)=0.5*(akm(i+1)+akm(i))
    501      bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
     516    akz(i+1)=0.5*(akm(i+1)+akm(i))
     517    bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
    502518  end do
    503519  nuvz=nuvz+1
     
    539555    if (pint.lt.5000.) goto 96
    540556  end do
    541 96   nconvlev=i
     55796 nconvlev=i
    542558  if (nconvlev.gt.nconvlevmax-1) then
    543559    nconvlev=nconvlevmax-1
     
    548564  return
    549565
    550 999   write(*,*)
     566999 write(*,*)
    551567  write(*,*) ' ###########################################'// &
    552568       '###### '
     
    568584  endif
    569585
    570 end subroutine gridcheck
    571 
     586end subroutine gridcheck_ecmwf
     587
  • src/gridcheck_gfs.f90

    r4c64400 rd8eed02  
    2020!**********************************************************************
    2121
    22 subroutine gridcheck
     22subroutine gridcheck_gfs
    2323
    2424  !**********************************************************************
     
    3838  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    3939  !                                 ECMWF grib_api                      *
     40  !                                                                     *
     41  !   Unified ECMWF and GFS builds                                      *
     42  !   Marian Harustak, 12.5.2017                                        *
     43  !     - Renamed routine from gridcheck to gridcheck_gfs               *
    4044  !                                                                     *
    4145  !**********************************************************************
     
    227231       yaux2in,iret)
    228232  call grib_check(iret,gribFunction,gribErrorMsg)
     233
     234  ! Fix for flexpart.eu ticket #48
     235  if (xaux2in.lt.0) xaux2in = 359.0
    229236
    230237  xaux1=xaux1in
     
    537544  endif
    538545
    539 end subroutine gridcheck
     546end subroutine gridcheck_gfs
  • src/makefile

    r6985a98 rd8eed02  
    1212#       'make -j ecmwf gcc=4.9',
    1313#    also set environment variable LD_LIBRARY_PATH to point to compiler libraries
     14#
     15#    Makefile was modified to produce unified executable for both ECMWF and GFS meteo data formats
     16#    gributils were included to detect format of meteo data
     17#
     18#    Cpp directives USE_MPIINPLACE were added to three source files. The effect of these directives
     19#    are to enable the MPI_IN_PLACE option only if compiled with a -DUSE_MPIINPLACE directive.
     20#    Otherwise, a safer option (which requires the allocation of another array) is used by default.
     21#    In makefile added the -x f95-cpp-input flag for compiling of cpp directives.
    1422#
    1523#  USAGE
    16 #    Compile serial FLEXPART (ECMWF)
    17 #      make [-j] ecmwf
    18 #
    19 #    Compile parallel FLEXPART (ECMWF)
    20 #      make [-j] ecmwf-mpi
     24#    Compile serial FLEXPART
     25#      make [-j] serial
     26#
     27#    Compile parallel FLEXPART
     28#      make [-j] mpi
    2129#     
    22 #    Compile for debugging parallel FLEXPART (ECMWF)
    23 #      make [-j] ecmwf-mpi-dbg
    24 #
    25 #    Compile serial FLEXPART (GFS)
    26 #      make [-j] gfs
    27 #
    28 #    Compile parallel FLEXPART (GFS)
    29 #      make [-j] gfs-mpi
     30#    Compile for debugging parallel FLEXPART
     31#      make [-j] mpi-dbg
    3032#
    3133################################################################################
    3234
    3335## PROGRAMS
    34 FLEXPART-ECMWF-MPI      = FP_ecmwf_MPI
    35 FLEXPART-ECMWF-MPI-DBG  = DBG_FP_ecmwf_MPI
    36 FLEXPART-ECMWF          = FP_ecmwf_gfortran
    37 FLEXPART-GFS            = FP_gfs_gfortran
    38 FLEXPART-GFS-MPI        = FP_gfs_MPI
     36# Unified executable names
     37# The same executable is used for both ECMWF and GFS metdata
     38
     39# Parallel processing executable
     40FLEXPART-MPI = FLEXPART_MPI
     41
     42# Parallel processing executable with debugging info
     43FLEXPART-MPI-DBG = DBG_FLEXPART_MPI
     44
     45# Serial processing executable
     46FLEXPART-SERIAL = FLEXPART
     47
    3948
    4049ifeq ($(gcc), 4.9)
     
    6170
    6271
     72# path to gributils used to detect meteodata format
     73VPATH = gributils/
     74
     75
    6376## OPTIMIZATION LEVEL
    6477O_LEV = 2 # [0,1,2,3,g,s,fast]
     
    6881LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
    6982
    70 FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
    71 
    72 DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
    73 
    74 LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
     83FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER)  #-Warray-bounds -fcheck=all # -march=native
     84
     85DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace   -Wall  -fdump-core $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero  -Warray-bounds -fcheck=all
     86
     87LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
    7588LDDEBUG  = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
    7689
     
    8295xmass_mod.o             flux_mod.o \
    8396point_mod.o             outg_mod.o \
    84 mean_mod.o              random_mod.o
     97mean_mod.o              random_mod.o \
     98class_gribfile_mod.o
    8599
    86100MPI_MODOBJS = \
     
    99113        redist.o                \
    100114        concoutput_surf.o       concoutput_surf_nest.o  \
    101         getfields.o
     115        getfields.o \
     116        readwind_ecmwf.o
    102117
    103118## For MPI version
     
    112127        redist_mpi.o            \
    113128        concoutput_surf_mpi.o   concoutput_surf_nest_mpi.o      \
    114         getfields_mpi.o
    115 
    116 ### WINDFIELDS
    117 ## For ECMWF (serial) version:
    118 OBJECTS_ECMWF = \
    119         calcpar.o          readwind.o \
    120         richardson.o       verttransform.o \
    121         obukhov.o          gridcheck.o  \
    122         convmix.o          calcmatrix.o \
    123         ecmwf_mod.o
    124 
    125 
    126 ## For ECMWF MPI version:
    127 OBJECTS_ECMWF_MPI = \
    128         gridcheck.o        readwind_mpi.o \
    129         calcpar.o          \
    130         richardson.o       verttransform.o \
    131         obukhov.o          \
    132         convmix.o          calcmatrix.o \
    133         ecmwf_mod.o
    134 
    135 ## For GFS (serial) version:
    136 OBJECTS_GFS = \
    137         calcpar_gfs.o          readwind_gfs.o \
    138         richardson_gfs.o       verttransform_gfs.o \
    139         obukhov_gfs.o          gridcheck_gfs.o  \
    140         convmix_gfs.o          calcmatrix_gfs.o \
    141         gfs_mod.o
     129        getfields_mpi.o \
     130        readwind_ecmwf_mpi.o
    142131
    143132OBJECTS = \
     
    154143ew.o                    readreleases.o  \
    155144readdepo.o              get_vdep_prob.o   \
    156 get_wetscav.o   \
     145get_wetscav.o           readwind_gfs.o \
    157146psim.o                  outgrid_init.o  \
    158 outgrid_init_nest.o     \
     147outgrid_init_nest.o     calcmatrix.o \
    159148photo_O1D.o             readlanduse.o \
    160149interpol_wind.o         readoutgrid.o \
    161150interpol_all.o          readpaths.o \
    162 getrb.o                 \
    163 getrc.o                 \
     151getrb.o                 obukhov.o \
     152getrc.o                 convmix.o \
    164153getvdep.o               readspecies.o \
    165 interpol_misslev.o      \
    166 scalev.o \
    167 pbl_profile.o           readOHfield.o\
    168 juldate.o               \
     154interpol_misslev.o      richardson.o \
     155scalev.o                verttransform_ecmwf.o \
     156pbl_profile.o           readOHfield.o \
     157juldate.o               verttransform_gfs.o \
    169158interpol_vdep.o         interpol_rain.o \
    170159hanna.o                 wetdepokernel.o \
    171                   wetdepo.o \
     160calcpar.o               wetdepo.o \
    172161hanna_short.o           windalign.o \
    173 hanna1.o                        \
    174                         gridcheck_nests.o \
     162hanna1.o                gridcheck_ecmwf.o \
     163gridcheck_gfs.o         gridcheck_nests.o \
    175164readwind_nests.o        calcpar_nests.o \
    176165verttransform_nests.o   interpol_all_nests.o \
    177166interpol_wind_nests.o   interpol_misslev_nests.o \
    178167interpol_vdep_nests.o   interpol_rain_nests.o \
    179 readageclasses.o        \
     168readageclasses.o        detectformat.o \
    180169calcfluxes.o            fluxoutput.o \
    181170qvsat.o                 skplin.o \
     
    201190%.o: %.mod
    202191
    203 ecmwf: $(FLEXPART-ECMWF)
    204 ecmwf: FC := $(F90)
    205 
    206 ecmwf-mpi: $(FLEXPART-ECMWF-MPI)
    207 ecmwf-mpi: FC := $(MPIF90)
    208 
    209 ecmwf-mpi-dbg: $(FLEXPART-ECMWF-MPI-DBG)
    210 ecmwf-mpi-dbg: FFLAGS := $(DBGFLAGS)
    211 ecmwf-mpi-dbg: LDFLAGS:= $(LDDEBUG)
    212 ecmwf-mpi-dbg: FC := $(MPIF90)
    213 
    214 gfs: $(FLEXPART-GFS)
    215 gfs: FC := $(F90)
    216 
    217 gfs-mpi: $(FLEXPART-GFS-MPI)
    218 gfs-mpi: FC := $(MPIF90)
    219 
    220 #all: $(FLEXPART-ECMWF)
    221 #all: $(FLEXPART-ECMWF-MPI)
    222 
    223 ## This allows for switching between ECMWF/GFS without editing source code.
    224 wind_mod = ecmwf_mod.o # default wind
    225 # ifeq ($(MAKECMDGOALS),ecmwf)
    226 # wind_mod = ecmwf_mod.o
    227 # endif
    228 # ifeq ($(MAKECMDGOALS),ecmwf-mpi)
    229 # wind_mod = ecmwf_mod.o
    230 # endif
    231 # ifeq ($(MAKECMDGOALS),ecmwf-mpi-dbg)
    232 # wind_mod = ecmwf_mod.o
    233 # endif
    234 
    235 ifeq ($(MAKECMDGOALS),gfs)
    236 wind_mod = gfs_mod.o
    237 endif
    238 ifeq ($(MAKECMDGOALS),gfs-mpi)
    239 wind_mod = gfs_mod.o
    240 endif
    241 
    242 $(FLEXPART-ECMWF): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF)
    243         +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF) $(LDFLAGS)
    244 
    245 $(FLEXPART-ECMWF-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_ECMWF_MPI)
     192# serial executable
     193serial: $(FLEXPART-SERIAL)
     194serial: FC := $(F90)
     195
     196# parallel processing executable
     197mpi: $(FLEXPART-MPI)
     198mpi: FC := $(MPIF90)
     199
     200# parallel processing with debugging info
     201mpi-dbg: $(FLEXPART-MPI-DBG)
     202mpi-dbg: FFLAGS := $(DBGFLAGS)
     203mpi-dbg: LDFLAGS:= $(LDDEBUG)
     204mpi-dbg: FC := $(MPIF90)
     205
     206$(FLEXPART-SERIAL): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL)
     207        +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(LDFLAGS)
     208
     209$(FLEXPART-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI)
    246210        +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
    247         $(OBJECTS_ECMWF_MPI) $(LDFLAGS)
    248 #       +$(FC) -o $@ *.o $(LDFLAGS)
    249 
    250 $(FLEXPART-ECMWF-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
    251         $(OBJECTS_ECMWF_MPI)
     211        $(LDFLAGS)
     212
     213$(FLEXPART-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI)
    252214        +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
    253         $(OBJECTS_ECMWF_MPI) $(LDFLAGS)
    254 
    255 $(FLEXPART-GFS): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS)
    256         +$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS) $(LDFLAGS)
    257 
    258 $(FLEXPART-GFS-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_GFS)
    259         +$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
    260         $(OBJECTS_GFS) $(LDFLAGS)
     215        $(LDFLAGS)
    261216
    262217%.o: %.f90
     
    267222
    268223cleanall:
    269         \rm -f *.o *.mod $(FLEXPART-ECMWF-MPI) $(FLEXPART-ECMWF-MPI-DBG) $(FLEXPART-ECMWF) \
    270         $(FLEXPART-GFS-MPI) $(FLEXPART-GFS)
     224        \rm -f *.o *.mod $(FLEXPART-MPI) $(FLEXPART-MPI-DBG) $(FLEXPART-SERIAL)
     225
    271226
    272227.SUFFIXES = $(SUFFIXES) .f90
     
    282237        random_mod.o
    283238calcfluxes.o: com_mod.o flux_mod.o outg_mod.o par_mod.o
    284 calcmatrix.o: com_mod.o conv_mod.o par_mod.o
    285 calcmatrix_gfs.o: com_mod.o conv_mod.o par_mod.o
    286 calcpar.o: com_mod.o par_mod.o
    287 calcpar_gfs.o: com_mod.o par_mod.o
     239calcmatrix.o: com_mod.o conv_mod.o par_mod.o class_gribfile_mod.o
     240calcpar.o: com_mod.o par_mod.o class_gribfile_mod.o
    288241calcpar_nests.o: com_mod.o par_mod.o
    289242calcpv.o: com_mod.o par_mod.o
     
    311264conv_mod.o: par_mod.o
    312265convect43c.o: conv_mod.o par_mod.o
    313 convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o
    314 convmix_gfs.o: com_mod.o conv_mod.o par_mod.o
     266convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o class_gribfile_mod.o
    315267coordtrafo.o: com_mod.o par_mod.o point_mod.o
     268detectformat.o: com_mod.o par_mod.o class_gribfile_mod.o
    316269distance.o: par_mod.o
    317270distance2.o: par_mod.o
     
    319272drydepokernel_nest.o: com_mod.o par_mod.o unc_mod.o
    320273erf.o: par_mod.o
    321 FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o
     274FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o class_gribfile_mod.o
    322275FLEXPART_MPI.o: com_mod.o conv_mod.o mpi_mod.o par_mod.o point_mod.o \
    323         random_mod.o netcdf_output_mod.o
     276        random_mod.o netcdf_output_mod.o class_gribfile_mod.o
    324277fluxoutput.o: com_mod.o flux_mod.o outg_mod.o par_mod.o
    325278get_settling.o: com_mod.o par_mod.o
    326 getfields.o: com_mod.o par_mod.o
    327 getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o
     279getfields.o: com_mod.o par_mod.o class_gribfile_mod.o
     280getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o class_gribfile_mod.o
    328281gethourlyOH.o: com_mod.o oh_mod.o par_mod.o
    329282getrb.o: par_mod.o
     
    331284getvdep.o: com_mod.o par_mod.o
    332285getvdep_nests.o: com_mod.o par_mod.o
    333 gridcheck.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
     286gridcheck_ecmwf.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
    334287gridcheck_emos.o: com_mod.o conv_mod.o par_mod.o
    335288gridcheck_fnl.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
     
    366319mpi_mod.o: com_mod.o par_mod.o unc_mod.o
    367320netcdf_output_mod.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
    368 obukhov.o: par_mod.o
    369 obukhov_gfs.o: par_mod.o
     321obukhov.o: par_mod.o class_gribfile_mod.o
    370322ohreaction.o: com_mod.o oh_mod.o par_mod.o
    371323openouttraj.o: com_mod.o par_mod.o point_mod.o
     
    374326outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o
    375327outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
    376 par_mod.o : $(wind_mod)
    377328part0.o: par_mod.o
    378329partdep.o: par_mod.o
     
    402353readreleases.o: com_mod.o par_mod.o point_mod.o xmass_mod.o
    403354readspecies.o: com_mod.o par_mod.o
    404 readwind.o: com_mod.o par_mod.o
     355readwind_ecmwf.o: com_mod.o par_mod.o
    405356readwind_emos.o: com_mod.o par_mod.o
    406357readwind_gfs.o: com_mod.o par_mod.o
    407358readwind_gfs_emos.o: com_mod.o par_mod.o
    408 readwind_mpi.o: com_mod.o mpi_mod.o par_mod.o
     359readwind_ecmwf_mpi.o: com_mod.o mpi_mod.o par_mod.o
    409360readwind_nests.o: com_mod.o par_mod.o
    410361readwind_nests_emos.o: com_mod.o par_mod.o
     
    415366releaseparticles_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o \
    416367        random_mod.o xmass_mod.o
    417 richardson.o: par_mod.o
    418 richardson_gfs.o: par_mod.o
     368richardson.o: par_mod.o class_gribfile_mod.o
    419369scalev.o: par_mod.o
    420370shift_field.o: par_mod.o
     
    425375        par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o
    426376unc_mod.o: par_mod.o
    427 verttransform.o: cmapf_mod.o com_mod.o par_mod.o
     377verttransform_ecmwf.o: cmapf_mod.o com_mod.o par_mod.o
    428378verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
    429379verttransform_nests.o: com_mod.o par_mod.o
  • src/mpi_mod.f90

    raa8c34a rbb579a9  
    5959!                                                                            *
    6060!*****************************************************************************
     61!                                                                            *
     62! Modification by DJM, 2017-05-09 - added #ifdef USE_MPIINPLACE cpp          *
     63! directive to mpif_tm_reduce_grid() to insure that MPI_IN_PLACE is          *
     64! used in the MPI_Reduce() only if specifically compiled with that           *
     65! directive.                                                                 *
     66!                                                                            *
     67!*****************************************************************************
    6168
    6269  use mpi
     
    24252432
    24262433! 2) Using in-place reduction
     2434
     2435!!!--------------------------------------------------------------------
     2436!!! DJM - 2017-05-09 change - MPI_IN_PLACE option for MPI_Reduce() causes
     2437!!! severe numerical problems in some cases.  I'm guessing there are memory
     2438!!! overrun issues in this complex code, but have so far failed to identify
     2439!!! a specific problem.  And/or, when one searches the Internet for this
     2440!!! issue, there is "some" hint that the implementation may be buggy. 
     2441!!!
     2442!!! At this point, with the following CPP USE_MPIINPLACE directive, the
     2443!!! default behaviour will be to not use the MPI_IN_PLACE option.
     2444!!! Users will have to compile with -DUSE_MPIINPLACE if they want that option.
     2445!!! Introduction of the CPP directives also requires that the code be compiled
     2446!!! with the "-x f95-cpp-input" option.
     2447!!!
     2448!!! Modification of this section requires the addition of array gridunc0, which
     2449!!! requires an allocation in outgrid_init.f90 and initial declaration in
     2450!!! unc_mod.f90.
     2451!!!--------------------------------------------------------------------
     2452
     2453#ifdef USE_MPIINPLACE
     2454
    24272455    if (lroot) then
    24282456      call MPI_Reduce(MPI_IN_PLACE, gridunc, grid_size3d, mp_sp, MPI_SUM, id_root, &
     
    24332461           & mp_comm_used, mp_ierr)
    24342462    end if
     2463
     2464#else
     2465
     2466      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
     2467           & mp_comm_used, mp_ierr)
     2468      if (lroot) gridunc = gridunc0
     2469
     2470#endif
    24352471
    24362472    if ((WETDEP).and.(ldirect.gt.0)) then
  • src/obukhov.f90

    re200b7a r6ecb30a  
    2020!**********************************************************************
    2121
    22 real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm)
     22real function obukhov(ps,tsurf,tdsurf,tlev,ustar,hf,akm,bkm,plev,metdata_format)
    2323
    2424  !********************************************************************
     
    2727  !                       Date:   1994-06-27                          *
    2828  !                                                                   *
    29   !  Update: A. Stohl, 2000-09-25, avoid division by zero by          *
    30   !  setting ustar to minimum value                                   *
     29  !     This program calculates Obukhov scale height from surface     *
     30  !     meteorological data and sensible heat flux.                   *
    3131  !                                                                   *
    3232  !********************************************************************
    3333  !                                                                   *
    34   !     This program calculates Obukhov scale height from surface     *
    35   !     meteorological data and sensible heat flux.                   *
     34  !  Update: A. Stohl, 2000-09-25, avoid division by zero by          *
     35  !  setting ustar to minimum value                                   *
     36  !  CHANGE: 17/11/2005 Caroline Forster NCEP GFS version             *
     37  !                                                                   *
     38  !   Unified ECMWF and GFS builds                                    *
     39  !   Marian Harustak, 12.5.2017                                      *
     40  !     - Merged obukhov and obukhov_gfs into one routine using       *
     41  !       if-then for meteo-type dependent code                       *
    3642  !                                                                   *
    3743  !********************************************************************
     
    4753  !     akm     ECMWF vertical discretization parameter               *
    4854  !     bkm     ECMWF vertical discretization parameter               *
     55  !     plev                                                          *
     56  !     metdata_format format of metdata (ecmwf/gfs)                  *
    4957  !                                                                   *
    5058  !********************************************************************
    5159
    5260  use par_mod
     61  use class_gribfile
    5362
    5463  implicit none
    5564
     65  integer :: metdata_format
    5666  real :: akm(nwzmax),bkm(nwzmax)
    5767  real :: ps,tsurf,tdsurf,tlev,ustar,hf,e,ew,tv,rhoa,plev
     
    6272  tv=tsurf*(1.+0.378*e/ps)               ! virtual temperature
    6373  rhoa=ps/(r_air*tv)                      ! air density
     74  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    6475  ak1=(akm(1)+akm(2))/2.
    6576  bk1=(bkm(1)+bkm(2))/2.
    6677  plev=ak1+bk1*ps                        ! Pressure level 1
     78  end if
    6779  theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature
    6880  if (ustar.le.0.) ustar=1.e-8
  • src/outgrid_init.f90

    r6b22af9 r6ecb30a  
    1919! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
    2020!**********************************************************************
     21!                                                                     *
     22! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to     *
     23! enable allocation of a gridunc0 array if required by MPI code in    *
     24! mpi_mod.f90                                                         *
     25!                                                                     *
     26!**********************************************************************
     27
    2128
    2229subroutine outgrid_init
     
    215222! Extra field for totals at MPI root process
    216223  if (lroot.and.mpi_mode.gt.0) then
    217     ! allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
    218     !      maxpointspec_act,nclassunc,maxageclass),stat=stat)
    219     ! if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0'
     224
     225#ifdef USE_MPIINPLACE
     226#else
     227    ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(),
     228    ! then an aux array is needed for parallel grid reduction
     229    allocate(gridunc0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
     230         maxpointspec_act,nclassunc,maxageclass),stat=stat)
     231    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc0'
     232#endif
    220233    if (ldirect.gt.0) then
    221234      allocate(wetgridunc0(0:numxgrid-1,0:numygrid-1,maxspec, &
  • src/par_mod.f90

    r6985a98 rd8eed02  
    3737
    3838module par_mod
    39 
    40 !************************************************************************
    41 ! wind_mod: is gfs_mod.f90 for target gfs, ecmwf_mod.f90 for target ecmwf
    42 !************************************************************************
    43   use wind_mod
    4439
    4540  implicit none
     
    146141  !*********************************************
    147142 
    148   ! Moved to ecmwf_mod.f90 (for ECMWF) / gfs_mod.f90 (for GFS)
     143!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new
     144  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
     145
     146!  integer,parameter :: nxmax=181,nymax=91,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
     147
     148!  INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data
     149
     150!  integer,parameter :: nxshift=359  ! for ECMWF
     151  integer,parameter :: nxshift=0     ! for GFS
     152
     153  !*********************************************
     154  ! Maximum dimensions of the nested input grids
     155  !*********************************************
     156
     157  integer,parameter :: maxnests=1,nxmaxn=451,nymaxn=226
     158
     159  ! nxmax,nymax        maximum dimension of wind fields in x and y
     160  !                    direction, respectively
     161  ! nuvzmax,nwzmax     maximum dimension of (u,v) and (w) wind fields in z
     162  !                    direction (for fields on eta levels)
     163  ! nzmax              maximum dimension of wind fields in z direction
     164  !                    for the transformed Cartesian coordinates
     165  ! nxshift            for global grids (in x), the grid can be shifted by
     166  !                    nxshift grid points, in order to accomodate nested
     167  !                    grids, and output grids overlapping the domain "boundary"
     168  !                    nxshift must not be negative; "normal" setting would be 0
     169
    149170 
    150171  integer,parameter :: nconvlevmax = nuvzmax-1
     
    269290  integer,parameter ::  icmv=-9999
    270291
    271 ! Parameters for testing
    272 !*******************************************
    273 !  integer :: verbosity=0
    274292
    275293end module par_mod
  • src/readwind_gfs.f90

    r8a65cb0 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine readwind(indj,n,uuh,vvh,wwh)
     22subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
    2323
    2424  !***********************************************************************
     
    3838  !*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    3939  !*                                 ECMWF grib_api                      *
     40  !                                                                      *
     41  !   Unified ECMWF and GFS builds                                       *
     42  !   Marian Harustak, 12.5.2017                                         *
     43  !     - Renamed routine from readwind to readwind_gfs                  *
    4044  !*                                                                     *
    4145  !***********************************************************************
     
    716720  stop 'Execution terminated'
    717721
    718 end subroutine readwind
     722end subroutine readwind_gfs
  • src/richardson.f90

    rc2162ce r6ecb30a  
    2121
    2222subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
    23        akz,bkz,hf,tt2,td2,h,wst,hmixplus)
     23       akz,bkz,hf,tt2,td2,h,wst,hmixplus,metdata_format)
    2424  !                        i    i    i     i    i    i    i
    2525  ! i   i  i   i   i  o  o     o
     
    4141  !     Meteor. 81, 245-269.                                                  *
    4242  !                                                                           *
     43  !****************************************************************************
     44  !                                                                           *
    4345  !     Update: 1999-02-01 by G. Wotawa                                       *
    4446  !                                                                           *
     
    4648  !     instead of first model level.                                         *
    4749  !     New input variables tt2, td2 introduced.                              *
     50  !                                                                           *
     51  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  *
     52  !                                                                           *
     53  !     Unified ECMWF and GFS builds                                          *
     54  !     Marian Harustak, 12.5.2017                                            *
     55  !       - Merged richardson and richardson_gfs into one routine using       *
     56  !         if-then for meteo-type dependent code                             *
    4857  !                                                                           *
    4958  !****************************************************************************
     
    5564  ! tv                         virtual temperature                            *
    5665  ! wst                        convective velocity scale                      *
     66  ! metdata_format             format of metdata (ecmwf/gfs)                  *
    5767  !                                                                           *
    5868  ! Constants:                                                                *
     
    6272
    6373  use par_mod
     74  use class_gribfile
    6475
    6576  implicit none
    6677
    67   integer :: i,k,nuvz,iter
     78  integer :: metdata_format
     79  integer :: i,k,nuvz,iter,llev,loop_start
    6880  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
    6981  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
     
    7789  iter=0
    7890
     91  if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
     92    ! NCEP version: find first model level above ground
     93    !**************************************************
     94
     95     llev = 0
     96     do i=1,nuvz
     97       if (psurf.lt.akz(i)) llev=i
     98     end do
     99     llev = llev+1
     100    ! sec llev should not be 1!
     101     if (llev.eq.1) llev = 2
     102     if (llev.gt.nuvz) llev = nuvz-1
     103    ! NCEP version
     104  end if
     105
     106
    79107  ! Compute virtual temperature and virtual potential temperature at
    80108  ! reference level (2 m)
     
    95123  ! Integrate z up to one level above zt
    96124  !*************************************
    97 
    98   do k=2,nuvz
     125  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
     126    loop_start=2
     127  else
     128    loop_start=llev
     129  end if
     130  do k=loop_start,nuvz
    99131    pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
    100132    tv=ttlev(k)*(1.+0.608*qvlev(k))
  • src/timemanager.f90

    ra76d954 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine timemanager
     22subroutine timemanager(metdata_format)
    2323
    2424  !*****************************************************************************
     
    5050  !   For compatibility with MPI version,                                      *
    5151  !   variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
     52  !  Unified ECMWF and GFS builds                                              *
     53  !   Marian Harustak, 12.5.2017                                               *
     54  !   - Added passing of metdata_format as it was needed by called routines    *
    5255  !*****************************************************************************
    5356  !                                                                            *
     
    8386  ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) =                           *
    8487  !                    spatial positions of trajectories                       *
     88  ! metdata_format     format of metdata (ecmwf/gfs)                           *
    8589  !                                                                            *
    8690  ! Constants:                                                                 *
     
    102106  implicit none
    103107
     108  integer :: metdata_format
    104109  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
    105110! integer :: ksp
     
    194199           write (*,*) 'timemanager> call convmix -- backward'
    195200        endif         
    196       call convmix(itime)
     201      call convmix(itime,metdata_format)
    197202        if (verbosity.gt.1) then
    198203          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     
    207212           write (*,*) 'timemanager> call getfields'
    208213    endif
    209     call getfields(itime,nstop1)
     214    call getfields(itime,nstop1,metdata_format)
    210215        if (verbosity.gt.1) then
    211216          CALL SYSTEM_CLOCK(count_clock)
     
    266271       write (*,*) 'timemanager> call convmix -- forward'
    267272     endif   
    268      call convmix(itime)
     273     call convmix(itime,metdata_format)
    269274   endif
    270275
  • src/timemanager_mpi.f90

    rb5127f9 rd8eed02  
    2020!**********************************************************************
    2121
    22 subroutine timemanager
     22subroutine timemanager(metdata_format)
    2323
    2424!*****************************************************************************
     
    5050!   MPI version                                                              *
    5151!   Variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
     52!  Unified ECMWF and GFS builds                                              *
     53!   Marian Harustak, 12.5.2017                                               *
     54!   - Added passing of metdata_format as it was needed by called routines    *
    5255!*****************************************************************************
    5356!                                                                            *
     
    8386! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) =                           *
    8487!                    spatial positions of trajectories                       *
     88! metdata_format     format of metdata (ecmwf/gfs)                           *
    8589!                                                                            *
    8690! Constants:                                                                 *
     
    103107  implicit none
    104108
     109  integer :: metdata_format
    105110  logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
    106111  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
     
    205210      endif
    206211! readwind process skips this step
    207       if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime)
     212      if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime,metdata_format)
    208213    endif
    209214
     
    218223    if (mp_measure_time) call mpif_mtime('getfields',0)
    219224
    220     call getfields(itime,nstop1,memstat)
     225    call getfields(itime,nstop1,memstat,metdata_format)
    221226
    222227    if (mp_measure_time) call mpif_mtime('getfields',1)
     
    346351        write (*,*) 'timemanager> call convmix -- forward'
    347352      endif
    348       call convmix(itime)
     353      call convmix(itime,metdata_format)
    349354    endif
    350355
  • src/unc_mod.f90

    r4c64400 r6ecb30a  
    1919! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
    2020!**********************************************************************
     21!                                                                     *
     22! DJM - 2017-05-09 - added #ifdef USE_MPIINPLACE cpp directive to     *
     23! enable declaration of a gridunc0 array if required by MPI code in   *
     24! mpi_mod.f90                                                         *
     25!                                                                     *
     26!**********************************************************************
    2127
    2228module unc_mod
     
    2733
    2834  real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc
     35#ifdef USE_MPIINPLACE
     36#else
     37  ! If MPI_IN_PLACE option is not used in mpi_mod.f90::mpif_tm_reduce_grid(),
     38  ! then an aux array is needed for parallel grid reduction
     39  real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc0
     40#endif
    2941  real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
    3042  real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
  • src/verttransform_gfs.f90

    r4fbe7a5 r6ecb30a  
    2020!**********************************************************************
    2121
    22 subroutine verttransform(n,uuh,vvh,wwh,pvh)
     22subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
    2323  !                      i  i   i   i   i
    2424  !*****************************************************************************
     
    4747  !  Changes, Bernd C. Krueger, Feb. 2001:
    4848  !   Variables tth and qvh (on eta coordinates) from common block
     49  !
     50  !   Unified ECMWF and GFS builds                                     
     51  !   Marian Harustak, 12.5.2017                                       
     52  !     - Renamed routine from verttransform to verttransform_gfs
     53  !
    4954  !*****************************************************************************
    5055  !                                                                            *
     
    531536
    532537
    533 end subroutine verttransform
     538end subroutine verttransform_gfs
  • src/writeheader.f90

    rf13406c r6ecb30a  
    3535  !*****************************************************************************
    3636  !                                                                            *
     37  !  Modified to remove TRIM around the output of flexversion so that          *
     38  !  it will be a constant length (defined in com_mod.f90) in output header    *
     39  !                                                                            *
     40  !     Don Morton, Boreal Scientific Computing                                *
     41  !     07 May 2017                                                            *
     42  !                                                                            *
     43  !*****************************************************************************
     44  !                                                                            *
    3745  ! Variables:                                                                 *
    3846  !                                                                            *
     
    6775
    6876  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime, trim(flexversion)
     77    write(unitheader) ibdate,ibtime, flexversion
    7078  else
    71     write(unitheader) iedate,ietime, trim(flexversion)
     79    write(unitheader) iedate,ietime, flexversion
    7280  endif
    7381
  • src/writeheader_nest.f90

    r4fbe7a5 r6ecb30a  
    3535  !*****************************************************************************
    3636  !                                                                            *
     37  !  Modified to remove TRIM around the output of flexversion so that          *
     38  !  it will be a constant length (defined in com_mod.f90) in output header    *
     39  !                                                                            *
     40  !     Don Morton, Boreal Scientific Computing                                *
     41  !     07 May 2017                                                            *
     42  !                                                                            *
     43  !*****************************************************************************
     44  !                                                                            *
    3745  ! Variables:                                                                 *
    3846  !                                                                            *
     
    6775
    6876  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime, trim(flexversion)
     77    write(unitheader) ibdate,ibtime, flexversion
    7078  else
    71     write(unitheader) iedate,ietime, trim(flexversion)
     79    write(unitheader) iedate,ietime, flexversion
    7280  endif
    7381
  • src/writeheader_nest_surf.f90

    rf13406c r6ecb30a  
    3535  !*****************************************************************************
    3636  !                                                                            *
     37  !  Modified to remove TRIM around the output of flexversion so that          *
     38  !  it will be a constant length (defined in com_mod.f90) in output header    *
     39  !                                                                            *
     40  !     Don Morton, Boreal Scientific Computing                                *
     41  !     07 May 2017                                                            *
     42  !                                                                            *
     43  !*****************************************************************************
     44  !                                                                            *
    3745  ! Variables:                                                                 *
    3846  !                                                                            *
     
    6775
    6876  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime,trim(flexversion)
     77    write(unitheader) ibdate,ibtime,flexversion
    7078  else
    71     write(unitheader) iedate,ietime,trim(flexversion)
     79    write(unitheader) iedate,ietime,flexversion
    7280  endif
    7381
  • src/writeheader_surf.f90

    rf13406c r6ecb30a  
    3535  !*****************************************************************************
    3636  !                                                                            *
     37  !  Modified to remove TRIM around the output of flexversion so that          *
     38  !  it will be a constant length (defined in com_mod.f90) in output header    *
     39  !                                                                            *
     40  !     Don Morton, Boreal Scientific Computing                                *
     41  !     07 May 2017                                                            *
     42  !                                                                            *
     43  !*****************************************************************************
     44  !                                                                            *
    3745  ! Variables:                                                                 *
    3846  !                                                                            *
     
    6775
    6876  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime, trim(flexversion)
     77    write(unitheader) ibdate,ibtime, flexversion
    7078  else
    71     write(unitheader) iedate,ietime, trim(flexversion)
     79    write(unitheader) iedate,ietime, flexversion
    7280  endif
    7381
  • tests/NILU/basic_tests

    rca350ba rdd6a4aa  
    11#WORKSPACE=/xnilu_wrk/flexbuild
    22#WORKSPACE=/home/ignacio/repos/
    3 FP_exec=$WORKSPACE/src/FP_ecmwf_gfortran
     3FP_exec=$WORKSPACE/src/FLEXPART
    44path_flextest=$WORKSPACE/flextest/
    55declare -a test_names=('1' 'HelloWorld' 'Fwd1' 'Fwd2' 'Bwd1' 'Volc' '2')
  • tests/NILU/run_tests

    rca350ba rdd6a4aa  
    99#FP_exec=/home/ignacio/repos/flexpart/src/FP_ecmwf_gfortran
    1010#FP_exec=/xnilu_wrk/flexbuild/tests/NILU/FP_ecmwf_gfortran
    11 FP_exec=$WORKSPACE/src/FP_ecmwf_gfortran
     11#FP_exec=$WORKSPACE/src/FP_ecmwf_gfortran
     12# CTBTO unified executable
     13FP_exec=$WORKSPACE/src/FLEXPART
    1214#path_flextest=/home/ignacio/repos/flextest/
    1315#path_flextest=/xnilu_wrk/flexbuild/flextest/
  • src/coordtrafo.f90

    re200b7a r1d072c0  
    4646
    4747  integer :: i,j,k
     48  real :: yrspc ! small real number relative to x
    4849
    4950  if (numpoint.eq.0) goto 30
     
    6465  ! CHECK IF RELEASE POINTS ARE WITHIN DOMAIN
    6566  !******************************************
    66 
     67 
     68  yrspc = spacing(real(nymin1,kind=sp))
     69 
    6770  do i=1,numpoint
    6871    if (sglobal.and.(ypoint1(i).lt.1.e-6)) ypoint1(i)=1.e-6
    69     if (nglobal.and.(ypoint2(i).gt.real(nymin1)-1.e-5)) &
    70          ypoint2(i)=real(nymin1)-1.e-5
    71   if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1)-1.e-6) &
    72        .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1)-1.e-6) &
     72    if (nglobal.and.(ypoint2(i).gt.real(nymin1,kind=dp)-1.e-5)) &
     73         ypoint2(i)=real(nymin1,kind=dp)-10*yrspc
     74    if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1,kind=dp)-1.e-6) &
     75       .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1,kind=dp)-yrspc) &
    7376       .or.((.not.xglobal).and.((xpoint1(i).lt.1.e-6).or. &
    74        (xpoint1(i).ge.real(nxmin1)-1.e-6).or.(xpoint2(i).lt.1.e-6).or. &
    75        (xpoint2(i).ge.real(nxmin1)-1.e-6)))) then
     77       (xpoint1(i).ge.real(nxmin1,kind=dp)-1.e-6).or.(xpoint2(i).lt.1.e-6).or. &
     78       (xpoint2(i).ge.real(nxmin1,kind=dp)-1.e-6)))) then
    7679      write(*,*) ' NOTICE: RELEASE POINT OUT OF DOMAIN DETECTED.'
    7780      write(*,*) ' IT IS REMOVED NOW ... '
    78       if (i.ge.1000) then
     81      if (i.le.1000) then
    7982         write(*,*) ' COMMENT: ',compoint(i)
    8083      else
  • src/readspecies.f90

    rc8fc724 raa8c34a  
    6767  character(len=16) :: pspecies
    6868  real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
    69   real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
     69  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pkao
    7070  real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
    71   integer :: readerror
     71  integer :: readerror, pspec_ass
    7272
    7373! declare namelist
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG