Changeset 37


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

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

Location:
branches/petra/src
Files:
1 added
17 edited

Legend:

Unmodified
Added
Removed

  • branches/petra/src/FLEXPART.f90

    r36 r37  
    6666
    6767  ! FLEXPART version string
    68   flexversion='Version 9.2_beta.r34 (2015-02-12)'
     68  flexversion='Version 9.2_beta.r3? (2015-02-25)'
    6969  verbosity=0
    7070

  • branches/petra/src/com_mod.f90

    r36 r37  
    3232!        Modifications: 15 August 2013 IP,
    3333!        2/2015 PS, add incremental deposition switch
     34!        2/2015 PS, change variables for cloud diagnostics
    3435!                                                                              *
    3536!*******************************************************************************
    … …  
    365366  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    366367  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
    367   integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    368   integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
    369 
    370 
    371368  ! uu,vv,ww [m/2]       wind components in x,y and z direction
    372369  ! uupol,vvpol [m/s]    wind components in polar stereographic projection
    … …  
    377374  ! drhodz [kg/m2]       vertical air density gradient
    378375  ! tth,qvh              tth,qvh on original eta levels
    379   ! clouds:   no cloud, no precipitation   0
    380   !      cloud, no precipitation      1
    381   !      rainout  conv/lsp dominated  2/3
    382   !      washout  conv/lsp dominated  4/5
    383376
    384377  ! pplev for the GFS version
    … …  
    406399  real :: oli(0:nxmax-1,0:nymax-1,1,2)
    407400  real :: diffk(0:nxmax-1,0:nymax-1,1,2)
     401  integer, dimension (0:nxmax-1,0:nymax-1,2) :: icloudbot, icloudthck
    408402
    409403  ! ps                   surface pressure
    … …  
    426420  ! oli [m]              inverse Obukhov length (1/L)
    427421  ! diffk [m2/s]         diffusion coefficient at reference height
     422  ! icloudbot (m)        cloud bottom height
     423  ! icloudthck (m)       cloud thickness   
    428424
    429425
    … …  
    484480  real :: qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    485481  real :: pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    486   integer(kind=1) :: cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,2,maxnests)
    487   integer :: cloudsnh(0:nxmaxn-1,0:nymaxn-1,2,maxnests)
    488482  real :: rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    489483  real :: drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    … …  
    514508  real :: diffkn(0:nxmaxn-1,0:nymaxn-1,1,2,maxnests)
    515509  real :: vdepn(0:nxmaxn-1,0:nymaxn-1,maxspec,2,maxnests)
     510  integer, dimension (0:nxmax-1,0:nymax-1,2,maxnests) :: icloudbotn, icloudthckn
    516511
    517512

  • branches/petra/src/getfields.f90

    r4 r37  
    4141  !   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
    4242  !   Function of nstop extended.
     43  !
     44  !  3/2015, PS: put in verbosity output
     45  !
    4346  !*****************************************************************************
    4447  !                                                                            *
    … …  
    8386  integer :: indmin = 1
    8487
     88  call verboutput(verbosity,1,'  getfields')
    8589
    8690  ! Check, if wind fields are available for the current time step
    … …  
    125129    do indj=indmin,numbwf-1
    126130       if (ldirect*wftime(indj+1).gt.ldirect*itime) then
     131           call verboutput(verbosity,1,'    call readwind')
    127132          call readwind(indj+1,memind(2),uuh,vvh,wwh)
     133           call verboutput(verbosity,1,'    call readwind_nests')
    128134          call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
     135           call verboutput(verbosity,1,'    call calcpar')
    129136          call calcpar(memind(2),uuh,vvh,pvh)
     137           call verboutput(verbosity,1,'    call calcpar_nests')
    130138          call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
     139           call verboutput(verbosity,1,'    call verttransform')
    131140          call verttransform(memind(2),uuh,vvh,wwh,pvh)
     141           call verboutput(verbosity,1,'    call verttransform_nests')
    132142          call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
    133143          memtime(2)=wftime(indj+1)
    … …  
    148158             (ldirect*wftime(indj+1).gt.ldirect*itime)) then
    149159           memind(1)=1
     160            call verboutput(verbosity,1,'    call readwind')
    150161           call readwind(indj,memind(1),uuh,vvh,wwh)
     162            call verboutput(verbosity,1,'    call readwind_nests')
    151163           call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
     164            call verboutput(verbosity,1,'    call calcpar')
    152165           call calcpar(memind(1),uuh,vvh,pvh)
     166            call verboutput(verbosity,1,'    call calcpar_nests')
    153167           call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
     168            call verboutput(verbosity,1,'    call verttransform')
    154169           call verttransform(memind(1),uuh,vvh,wwh,pvh)
     170            call verboutput(verbosity,1,'    call verttransform_nests')
    155171           call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
    156172           memtime(1)=wftime(indj)
    157173           memind(2)=2
     174            call verboutput(verbosity,1,'    call readwind')
    158175           call readwind(indj+1,memind(2),uuh,vvh,wwh)
     176            call verboutput(verbosity,1,'    call readwind_nests')
    159177           call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
     178            call verboutput(verbosity,1,'    call calcpar')
    160179           call calcpar(memind(2),uuh,vvh,pvh)
     180            call verboutput(verbosity,1,'    call calcpar_nests')
    161181           call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
     182            call verboutput(verbosity,1,'    call verttransform')
    162183           call verttransform(memind(2),uuh,vvh,wwh,pvh)
     184            call verboutput(verbosity,1,'    call verttransform_nests')
    163185           call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
    164186           memtime(2)=wftime(indj+1)
    … …  
    175197  if (lwindinterv.gt.idiffmax) nstop=3
    176198
     199  call verboutput(verbosity,1,'    leaving getfields')
     200
    177201end subroutine getfields

  • branches/petra/src/gridcheck.f90

    r27 r37  
    3232  !             LAST UPDATE: 1997-10-10                                 *
    3333  !                                                                     *
     34  !**********************************************************************
     35  !                                                                     *
    3436  !             Update:      1999-02-08, global fields allowed, A. Stohl*
    3537  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
    … …  
    3739  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
    3840  !                                 ECMWF grib_api                      *
    39   !                                                                     *
     41  !
     42  !   3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
     43  ! 
     44  !
    4045  !**********************************************************************
    4146  !                                                                     *
    … …  
    144149  if (gribVer.eq.1) then ! GRIB Edition 1
    145150
    146   !print*,'GRiB Edition 1'
    147   !read the grib2 identifiers
    148   call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    149   call grib_check(iret,gribFunction,gribErrorMsg)
    150   call grib_get_int(igrib,'level',isec1(8),iret)
    151   call grib_check(iret,gribFunction,gribErrorMsg)
    152 
    153   !change code for etadot to code for omega
    154   if (isec1(6).eq.77) then
    155     isec1(6)=135
    156   endif
    157 
    158   !print*,isec1(6),isec1(8)
    159 
    160   else
    161 
    162   !print*,'GRiB Edition 2'
    163   !read the grib2 identifiers
    164   call grib_get_int(igrib,'discipline',discipl,iret)
    165   call grib_check(iret,gribFunction,gribErrorMsg)
    166   call grib_get_int(igrib,'parameterCategory',parCat,iret)
    167   call grib_check(iret,gribFunction,gribErrorMsg)
    168   call grib_get_int(igrib,'parameterNumber',parNum,iret)
    169   call grib_check(iret,gribFunction,gribErrorMsg)
    170   call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    171   call grib_check(iret,gribFunction,gribErrorMsg)
    172   call grib_get_int(igrib,'level',valSurf,iret)
    173   call grib_check(iret,gribFunction,gribErrorMsg)
    174   call grib_get_int(igrib,'paramId',parId,iret)
    175   call grib_check(iret,gribFunction,gribErrorMsg)
    176 
    177   !print*,discipl,parCat,parNum,typSurf,valSurf
    178 
    179   !convert to grib1 identifiers
    180   isec1(6)=-1
    181   isec1(7)=-1
    182   isec1(8)=-1
    183   isec1(8)=valSurf     ! level
    184   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    185     isec1(6)=130         ! indicatorOfParameter
    186   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
    187     isec1(6)=131         ! indicatorOfParameter
    188   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
    189     isec1(6)=132         ! indicatorOfParameter
    190   elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    191     isec1(6)=133         ! indicatorOfParameter
    192   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    193     isec1(6)=134         ! indicatorOfParameter
    194   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    195     isec1(6)=135         ! indicatorOfParameter
    196   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    197     isec1(6)=135         ! indicatorOfParameter
    198   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    199     isec1(6)=151         ! indicatorOfParameter
    200   elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
    201     isec1(6)=165         ! indicatorOfParameter
    202   elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
    203     isec1(6)=166         ! indicatorOfParameter
    204   elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
    205     isec1(6)=167         ! indicatorOfParameter
    206   elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
    207     isec1(6)=168         ! indicatorOfParameter
    208   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    209     isec1(6)=141         ! indicatorOfParameter
    210   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    211     isec1(6)=164         ! indicatorOfParameter
    212   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    213     isec1(6)=142         ! indicatorOfParameter
    214   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    215     isec1(6)=143         ! indicatorOfParameter
    216   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    217     isec1(6)=146         ! indicatorOfParameter
    218   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    219     isec1(6)=176         ! indicatorOfParameter
    220   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    221     isec1(6)=180         ! indicatorOfParameter
    222   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    223     isec1(6)=181         ! indicatorOfParameter
    224   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    225     isec1(6)=129         ! indicatorOfParameter
    226   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    227     isec1(6)=160         ! indicatorOfParameter
    228   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    229        (typSurf.eq.1)) then ! LSM
    230     isec1(6)=172         ! indicatorOfParameter
    231   else
    232     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
    233          parCat,parNum,typSurf
    234   endif
    235   if(parId .ne. isec1(6) .and. parId .ne. 77) then
    236     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
    237 !    stop
    238   endif
    239 
    240   endif
     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)
     164
     165  else ! GRIB 2
     166
     167    !print*,'GRiB Edition 2'
     168    ! read the grib2 identifiers
     169
     170    call grib_get_int(igrib,'discipline',discipl,iret)
     171    call grib_check(iret,gribFunction,gribErrorMsg)
     172    call grib_get_int(igrib,'parameterCategory',parCat,iret)
     173    call grib_check(iret,gribFunction,gribErrorMsg)
     174    call grib_get_int(igrib,'parameterNumber',parNum,iret)
     175    call grib_check(iret,gribFunction,gribErrorMsg)
     176    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
     177    call grib_check(iret,gribFunction,gribErrorMsg)
     178    call grib_get_int(igrib,'level',valSurf,iret)
     179    call grib_check(iret,gribFunction,gribErrorMsg)
     180    call grib_get_int(igrib,'paramId',parId,iret)
     181    call grib_check(iret,gribFunction,gribErrorMsg)
     182
     183    !print*,discipl,parCat,parNum,typSurf,valSurf
     184
     185    ! convert to grib1 identifiers
     186
     187    isec1(6)=-1
     188    isec1(7)=-1
     189    isec1(8)=-1
     190    isec1(8)=valSurf     ! level
     191    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     192      isec1(6)=130         ! indicatorOfParameter
     193    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
     194      isec1(6)=131         ! indicatorOfParameter
     195    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
     196      isec1(6)=132         ! indicatorOfParameter
     197    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     198      isec1(6)=133         ! indicatorOfParameter
     199    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     200      isec1(6)=134         ! indicatorOfParameter
     201    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta
     202      isec1(6)=135         ! indicatorOfParameter
     203    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
     204      isec1(6)=135         ! indicatorOfParameter
     205    elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
     206      isec1(6)=135         ! indicatorOfParameter
     207    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     208      isec1(6)=151         ! indicatorOfParameter
     209    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
     210      isec1(6)=165         ! indicatorOfParameter
     211    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
     212      isec1(6)=166         ! indicatorOfParameter
     213    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
     214      isec1(6)=167         ! indicatorOfParameter
     215    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
     216      isec1(6)=168         ! indicatorOfParameter
     217    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     218      isec1(6)=141         ! indicatorOfParameter
     219    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     220      isec1(6)=164         ! indicatorOfParameter
     221    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     222      isec1(6)=142         ! indicatorOfParameter
     223    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     224      isec1(6)=143         ! indicatorOfParameter
     225    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     226      isec1(6)=146         ! indicatorOfParameter
     227    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     228      isec1(6)=176         ! indicatorOfParameter
     229    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     230      isec1(6)=180         ! indicatorOfParameter
     231    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     232      isec1(6)=181         ! indicatorOfParameter
     233    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     234      isec1(6)=129         ! indicatorOfParameter
     235    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     236      isec1(6)=160         ! indicatorOfParameter
     237    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
     238         (typSurf.eq.1)) then ! LSM
     239      isec1(6)=172         ! indicatorOfParameter
     240    else
     241      print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     242           parCat,parNum,typSurf
     243    endif
     244    if(parId .ne. isec1(6) .and. parId .ne. 77) then
     245      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     246  !    stop
     247    endif
     248
     249  endif ! GRIB 1 / GRIB 2 cases
    241250
    242251  !get the size and data of the values array

  • branches/petra/src/interpol_rain.f90

    r24 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    2020!**********************************************************************
    2121
    22 subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
    23        ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
     22subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx,ny, &
     23       memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
     24       icbot,ictop)
    2425  !                          i   i   i    i    i     i   i
    2526  !i    i    i  i    i     i      i      i     o     o     o
    … …  
    4041  !                                                                           *
    4142  !     30 August 1996                                                        *
     43  !****************************************************************************
     44  ! Petra Seibert, 2011/2012-2015:
     45  !  Add interpolation of cloud bottom and cloud thickness
     46  !  for fix to SE's new wet scavenging scheme
    4247  !                                                                           *
    4348  !****************************************************************************
    … …  
    6671  !                                                                           *
    6772  !****************************************************************************
     73 
     74  use com_mod, only: icloudbot, icloudthck
     75  use par_mod, only: icmv
    6876
    6977  implicit none
    7078
    7179  integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
    72   integer :: itime,itime1,itime2,level,indexh
     80  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
     81  integer :: icbot,icthck,ictop,ipsum
    7382  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    7483  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    7584  real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
    76   real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
     85  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2)
    7786  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
    78 
    7987
    8088
    … …  
    112120  !***********************
    113121
     122  memind_loop: &
    114123  do m=1,2
     124
    115125    indexh=memind(m)
    116126
    … …  
    127137         + p3*yy3(ix ,jyp,level,indexh) &
    128138         + p4*yy3(ixp,jyp,level,indexh)
    129   end do
     139
     140
     141! PS clouds:
     142    ip1=1
     143    ip2=1
     144    ip3=1
     145    ip4=1
     146    if (icloudbot(ix ,jy ,indexh) .eq. icmv) ip1=0
     147    if (icloudbot(ixp,jy ,indexh) .eq. icmv) ip2=0
     148    if (icloudbot(ix ,jyp,indexh) .eq. icmv) ip3=0
     149    if (icloudbot(ixp,jyp,indexh) .eq. icmv) ip4=0
     150    ipsum = ip1+ip2+ip3+ip4
     151    if (ipsum .eq. 0) then
     152      cbot(m)=icmv
     153    else
     154      cbot(m)=(ip1*p1*icloudbot(ix ,jy ,indexh) &
     155             + ip2*p2*icloudbot(ixp,jy ,indexh) &
     156             + ip3*p3*icloudbot(ix ,jyp,indexh) &
     157             + ip4*p4*icloudbot(ixp,jyp,indexh)) / ipsum
     158    endif
     159
     160    ip1=1
     161    ip2=1
     162    ip3=1
     163    ip4=1
     164    if (icloudthck(ix ,jy ,indexh) .eq. icmv) ip1=0
     165    if (icloudthck(ixp,jy ,indexh) .eq. icmv) ip2=0
     166    if (icloudthck(ix ,jyp,indexh) .eq. icmv) ip3=0
     167    if (icloudthck(ixp,jyp,indexh) .eq. icmv) ip4=0
     168    ipsum = ip1+ip2+ip3+ip4
     169    if (ipsum .eq. 0) then
     170      cthck(m)=icmv
     171    else
     172      cthck(m)=(ip1*p1*icloudthck(ix ,jy ,indexh) &
     173              + ip2*p2*icloudthck(ixp,jy ,indexh) &
     174              + ip3*p3*icloudthck(ix ,jyp,indexh) &
     175              + ip4*p4*icloudthck(ixp,jyp,indexh)) / ipsum
     176    endif
     177! PS end clouds
     178
     179  enddo memind_loop
    130180
    131181
    … …  
    142192  yint3=(y3(1)*dt2+y3(2)*dt1)/dt
    143193
     194! PS clouds:
     195  icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt )
     196  if (nint(cbot(1)) .eq. icmv) icbot=cbot(2)
     197  if (nint(cbot(2)) .eq. icmv) icbot=cbot(1)
     198
     199  icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt )
     200  if (nint(cthck(1)) .eq. icmv) icthck=cthck(2)
     201  if (nint(cthck(2)) .eq. icmv) icthck=cthck(1)
     202
     203  if (icbot .ne. icmv .and. icthck .ne. icmv) then
     204    ictop = icbot + icthck ! convert cloud thickness to cloud top
     205  else
     206    icbot=icmv
     207    ictop=icmv
     208  endif
     209! PS end clouds
     210
    144211
    145212end subroutine interpol_rain

  • branches/petra/src/interpol_rain_nests.f90

    r4 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    2121
    2222subroutine interpol_rain_nests(yy1,yy2,yy3,nxmaxn,nymaxn,nzmax, &
    23        maxnests,ngrid,nxn,nyn,memind,xt,yt,level,itime1,itime2,itime, &
    24        yint1,yint2,yint3)
     23       maxnests,ngrid,nxn,nyn, &
     24       memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
     25       icbot,ictop)
    2526  !                                i   i   i    i      i      i
    2627  !   i       i    i   i    i    i  i    i     i      i      i
    … …  
    4445  !                                                                           *
    4546  !     15 March 2000                                                         *
     47  !                                                                           *
     48  !****************************************************************************
     49  ! Petra Seibert, 2011/2012-2015:
     50  !  Add interpolation of cloud bottom and cloud thickness
     51  !  for fix to SE's new wet scavenging scheme
    4652  !                                                                           *
    4753  !****************************************************************************
    … …  
    7177  !****************************************************************************
    7278
     79  use com_mod, only: icloudbotn, icloudthckn
     80  use par_mod, only: icmv
     81
    7382  implicit none
    7483
    7584  integer :: maxnests,ngrid
    76   integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2)
    77   integer :: m,ix,jy,ixp,jyp,itime,itime1,itime2,level,indexh
     85  integer :: nxn(maxnests),nyn(maxnests),nxmaxn,nymaxn,nzmax,memind(2),m
     86  integer :: ix,jy,ixp,jyp
     87  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
     88  integer :: icbot,ictop,icthck,ipsum
    7889  real :: yy1(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    7990  real :: yy2(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    8091  real :: yy3(0:nxmaxn-1,0:nymaxn-1,nzmax,2,maxnests)
    81   real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
     92  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),cbot(2),cthck(2)
    8293  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
    83 
    84 
    8594
    8695  ! If point at border of grid -> small displacement into grid
    8796  !***********************************************************
    8897
    89   if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
    90        xt=real(nxn(ngrid)-1)-0.0001
    91   if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
    92        yt=real(nyn(ngrid)-1)-0.0001
     98  if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) xt=real(nxn(ngrid)-1)-0.0001
     99  if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) yt=real(nyn(ngrid)-1)-0.0001
    93100
    94101
    … …  
    119126  !***********************
    120127
     128  memind_loop: &
    121129  do m=1,2
     130
    122131    indexh=memind(m)
    123132
    … …  
    134143         + p3*yy3(ix ,jyp,level,indexh,ngrid) &
    135144         + p4*yy3(ixp,jyp,level,indexh,ngrid)
    136   end do
     145
     146! PS clouds:
     147    ip1=1
     148    ip2=1
     149    ip3=1
     150    ip4=1
     151    if (icloudbotn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
     152    if (icloudbotn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
     153    if (icloudbotn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
     154    if (icloudbotn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
     155    ipsum = ip1+ip2+ip3+ip4
     156    if (ipsum .eq. 0) then
     157      cbot(m)=icmv
     158    else
     159      cbot(m)=(ip1*p1*icloudbotn(ix ,jy ,indexh,ngrid) &
     160             + ip2*p2*icloudbotn(ixp,jy ,indexh,ngrid) &
     161             + ip3*p3*icloudbotn(ix ,jyp,indexh,ngrid) &
     162             + ip4*p4*icloudbotn(ixp,jyp,indexh,ngrid)) / ipsum
     163    endif
     164
     165    ip1=1
     166    ip2=1
     167    ip3=1
     168    ip4=1
     169    if (icloudthckn(ix ,jy ,indexh,ngrid) .eq. icmv) ip1=0
     170    if (icloudthckn(ixp,jy ,indexh,ngrid) .eq. icmv) ip2=0
     171    if (icloudthckn(ix ,jyp,indexh,ngrid) .eq. icmv) ip3=0
     172    if (icloudthckn(ixp,jyp,indexh,ngrid) .eq. icmv) ip4=0
     173    ipsum = ip1+ip2+ip3+ip4
     174    if (ipsum .eq. 0) then
     175      cthck(m)=icmv
     176    else
     177      cthck(m)=(ip1*p1*icloudthckn(ix ,jy ,indexh,ngrid) &
     178              + ip2*p2*icloudthckn(ixp,jy ,indexh,ngrid) &
     179              + ip3*p3*icloudthckn(ix ,jyp,indexh,ngrid) &
     180              + ip4*p4*icloudthckn(ixp,jyp,indexh,ngrid)) / ipsum
     181    endif
     182! PS end clouds
     183
     184  enddo memind_loop
    137185
    138186
    … …  
    150198
    151199
     200! PS clouds:
     201  icbot = nint( (cbot(1)*dt2 + cbot(2)*dt1)/dt )
     202  if (nint(cbot(1)) .eq. icmv) icbot=cbot(2)
     203  if (nint(cbot(2)) .eq. icmv) icbot=cbot(1)
     204
     205  icthck = nint( (cthck(1)*dt2 + cthck(2)*dt1)/dt )
     206  if (nint(cthck(1)) .eq. icmv) icthck=cthck(2)
     207  if (nint(cthck(2)) .eq. icmv) icthck=cthck(1)
     208
     209  if (icbot .ne. icmv .and. icthck .ne. icmv) then
     210    ictop = icbot + icthck ! convert cloud thickness to cloud top
     211  else
     212    icbot=icmv
     213    ictop=icmv
     214  endif
     215! PS end clouds
     216
     217
    152218end subroutine interpol_rain_nests

  • branches/petra/src/par_mod.f90

    r36 r37  
    103103
    104104  integer,parameter :: idiffnorm=10800, idiffmax=2*idiffnorm, minstep=1
     105  integer,parameter :: itagekernmin=10800
    105106
    106107  ! idiffnorm [s]           normal time interval between two wind fields
    107108  ! idiffmax [s]            maximum time interval between two wind fields
    108109  ! minstep [s]             minimum time step to be used within FLEXPART
    109 
     110  ! itagekernmin [s]        minimum particle age [s] for using output kernel
    110111
    111112  !*****************************************************************
    … …  
    265266
    266267!******************************************************
    267 ! integer code for missing values, used in wet scavenging (PS, 2012)
     268! for use in wet scavenging (PS, 2012, 2015)
    268269!******************************************************
    269270
    270     !  integer icmv
    271     !  parameter(icmv=-9999)
    272       integer,parameter ::  icmv=-9999
     271  integer,parameter :: icmv=-9999 ! integer code for missing values
     272  ! some cloud heights used when trying to construct reasonable cloud:
     273  integer,parameter :: icloudtopconvmin = 6000
     274  integer,parameter :: icloudtopmin = 3000
     275  integer,parameter :: icloudbot1 = 500, icloudtop1 =  8000
     276  integer,parameter :: icloudbot2 = 0,   icloudtop2 = 10000
     277 
     278  real, parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
     279  real, parameter :: rhmininit = 0.90 ! standard condition for presence of clouds
     280! PS       note that original by Sabine Eckhart was 80%
     281! PS       however, for T<-20 C we consider saturation over ice
     282! PS       so I think 90% should be enough         
     283  real, parameter :: rhminx = 0.30 ! extreme condition for presence of clouds
    273284
    274285! Parameters for testing

  • branches/petra/src/readcommand.f90

    r36 r37  
    235235    read(unitcommand,*,iostat=icmdstat) linit_cond
    236236    if (icmdstat .gt. 0) &
    237       print*, 'readcommand: linit_cond not read properly',icmdstat,linit_cond
    238     if (old) call skplin(3,unitcommand)
    239     read(unitcommand,*,iostat=icmdstat) surf_only
     237      print*, 'readcommand: linit_cond not present or no read properly,&
     238       & - check values: ',icmdstat,linit_cond
     239    if (old) call skplin(3,unitcommand)
     240    read(unitcommand,*,iostat=icmdstat) icmdstat,surf_only
    240241    if (icmdstat .gt. 0) &
    241       print*, 'readcommand: linit_cond not read properly',icmdstat,surf_only
    242     if (old) call skplin(3,unitcommand)
    243     read(unitcommand,*,iostat=icmdstat) ldep_incr
     242      print*, 'readcommand: surf_only not present or not read properly&
     243       & - check values: ',icmdstat,surf_only
     244    if (old) call skplin(3,unitcommand)
     245    read(unitcommand,*,iostat=icmdstat) icmdstat,ldep_incr
    244246    if (icmdstat .gt. 0) &
    245       print*, 'readcommand: linit_cond not read properly',icmdstat, ldep_incr
     247      print*, 'readcommand: ldep_incr not present or not read properly&
     248       & - check values: ',icmdstat, ldep_incr
    246249    close(unitcommand)
    247250
    … …  
    416419
    417420
    418   ! For backward runs one releasefield for all releases makes no sense,
    419   ! For quasilag and domainfill ioutputforechrelease is forbidden
    420   !*****************************************************************************
    421 
    422   if ((ldirect .lt. 0) .and. (ioutputforeachrelease .eq. 0)) then
    423       write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND:     ####'
    424       write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
    425       write(*,*) '#### MUST BE SET TO ONE!                     ####'
    426       stop
    427   endif
    428 
    429 
    430   ! For backward runs one releasefield for all releases makes no sense,
     421  ! For domainfill runs one releasefield for all releases makes no sense,
    431422  ! and is "forbidden"
    432423  !*****************************************************************************

  • branches/petra/src/readspecies.f90

    r36 r37  
    2020!**********************************************************************
    2121
    22 subroutine readspecies(id_spec,pos_spec)
     22subroutine readspecies(id_spec,id_pos)
    2323
    2424  !*****************************************************************************
    2525  !                                                                            *
    2626  !     This routine reads names and physical constants of chemical species/   *
    27   !     radionuclides given in the parameter pos_spec                          *
     27  !     radionuclides given in the parameter id_pos                          *
    2828  !                                                                            *
    2929  !   Author: A. Stohl                                                         *
    … …  
    3232  !
    3333  !   Changes:                                                                 *
     34  !                                                                            *
    3435  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
    3536  !                                                                            *
    3637  !   HSO, 13 August 2013: added optional namelist input
    3738  !   PS, 2/2015: access= -> position=
     39  !   PS, 4/2015: quick fix for bug associated with end=22
    3840  !                                                                            *
    3941  !*****************************************************************************
    … …  
    6062  implicit none
    6163
    62   integer :: i, pos_spec,j
     64  integer :: i, id_pos,j
    6365  integer :: idow,ihour,id_spec
    6466  character(len=3) :: aspecnumb
    … …  
    6971  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
    7072  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
    71   integer :: readerror
     73  integer :: ios
    7274
    7375  ! declare namelist
    … …  
    100102  ! Open the SPECIES file and read species names and properties
    101103  !************************************************************
    102   specnum(pos_spec)=id_spec
    103   write(aspecnumb,'(i3.3)') specnum(pos_spec)
     104  specnum(id_pos)=id_spec
     105  write(aspecnumb,'(i3.3)') specnum(id_pos)
    104106  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
    105   !write(*,*) 'reading SPECIES',specnum(pos_spec)
     107  !write(*,*) 'reading SPECIES',specnum(id_pos)
    106108
    107109  ASSSPEC=.FALSE.
    108110
    109111  ! try namelist input
    110   read(unitspecies,species_params,iostat=readerror)
     112  read(unitspecies,species_params,iostat=ios)
    111113  close(unitspecies)
    112114   
    113   if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
    114 
    115     readerror=1
    116 
    117     open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
     115  if ((pweightmolar.eq.-789.0).or.(ios.ne.0)) then ! no namelist found
     116
     117    open(unitspecies,file=path(1) (1:length(1))//'SPECIES/SPECIES_'//aspecnumb, status='old', err=998)
    118118
    119119    do i=1,6
    120120      read(unitspecies,*)
    121121    end do
    122 
    123     read(unitspecies,'(a10)',end=22) species(pos_spec)
    124   !  write(*,*) species(pos_spec)
    125     read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
    126   !  write(*,*) decay(pos_spec)
    127     read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
    128   !  write(*,*) weta(pos_spec)
    129     read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    130   !  write(*,*) wetb(pos_spec)
     122    read(unitspecies,'(a10)',end=20) species(id_pos)
     123  !  write(*,*) species(id_pos)
     124    read(unitspecies,'(f18.1)',end=20) decay(id_pos)
     125  !  write(*,*) decay(id_pos)
     126    read(unitspecies,'(e18.1)',end=20) weta(id_pos)
     127  !  write(*,*) weta(id_pos)
     128    read(unitspecies,'(f18.2)',end=20) wetb(id_pos)
     129  !  write(*,*) wetb(id_pos)
    131130
    132131  !*** NIK 31.01.2013: including in-cloud scavening parameters
    133    read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
    134   !  write(*,*) weta_in(pos_spec)
    135    read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
    136   !  write(*,*) wetb_in(pos_spec)
    137    read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
    138   !  write(*,*) wetc_in(pos_spec)
    139    read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
    140   !  write(*,*) wetd_in(pos_spec)
    141 
    142     read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    143   !  write(*,*) reldiff(pos_spec)
    144     read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
    145   !  write(*,*) henry(pos_spec)
    146     read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
    147   !  write(*,*) f0(pos_spec)
    148     read(unitspecies,'(e18.1)',end=22) density(pos_spec)
    149   !  write(*,*) density(pos_spec)
    150     read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
    151   !  write(*,*) dquer(pos_spec)
    152     read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
    153   !  write(*,*) dsigma(pos_spec)
    154     read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
    155   !  write(*,*) dryvel(pos_spec)
    156     read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
    157   !  write(*,*) weightmolar(pos_spec)
    158     read(unitspecies,'(e18.1)',end=22) ohreact(pos_spec)
    159   !  write(*,*) ohreact(pos_spec)
    160     read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
    161   !  write(*,*) spec_ass(pos_spec)
    162     read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
    163   !       write(*,*) kao(pos_spec)
    164 
    165     pspecies=species(pos_spec)
    166     pdecay=decay(pos_spec)
    167     pweta=weta(pos_spec)
    168     pwetb=wetb(pos_spec)
    169     pweta_in=weta_in(pos_spec)
    170     pwetb_in=wetb_in(pos_spec)
    171     pwetc_in=wetc_in(pos_spec)
    172     pwetd_in=wetd_in(pos_spec)
    173     preldiff=reldiff(pos_spec)
    174     phenry=henry(pos_spec)
    175     pf0=f0(pos_spec)
    176     pdensity=density(pos_spec)
    177     pdquer=dquer(pos_spec)
    178     pdsigma=dsigma(pos_spec)
    179     pdryvel=dryvel(pos_spec)
    180     pweightmolar=weightmolar(pos_spec)
    181     pohreact=ohreact(pos_spec)
    182     pspec_ass=spec_ass(pos_spec)
    183     pkao=kao(pos_spec)
     132   read(unitspecies,'(e18.1)',end=20) weta_in(id_pos)
     133  !  write(*,*) weta_in(id_pos)
     134   read(unitspecies,'(f18.2)',end=20) wetb_in(id_pos)
     135  !  write(*,*) wetb_in(id_pos)
     136   read(unitspecies,'(f18.2)',end=20) wetc_in(id_pos)
     137  !  write(*,*) wetc_in(id_pos)
     138   read(unitspecies,'(f18.2)',end=20) wetd_in(id_pos)
     139  !  write(*,*) wetd_in(id_pos)
     140
     141    read(unitspecies,'(f18.1)',end=20) reldiff(id_pos)
     142  !  write(*,*) reldiff(id_pos)
     143    read(unitspecies,'(e18.1)',end=20) henry(id_pos)
     144  !  write(*,*) henry(id_pos)
     145    read(unitspecies,'(f18.1)',end=20) f0(id_pos)
     146  !  write(*,*) f0(id_pos)
     147    read(unitspecies,'(e18.1)',end=20) density(id_pos)
     148  !  write(*,*) density(id_pos)
     149    read(unitspecies,'(e18.1)',end=20) dquer(id_pos)
     150  !  write(*,*) dquer(id_pos)
     151    read(unitspecies,'(e18.1)',end=20) dsigma(id_pos)
     152  !  write(*,*) dsigma(id_pos)
     153    read(unitspecies,'(f18.2)',end=20) dryvel(id_pos)
     154  !  write(*,*) dryvel(id_pos)
     155    read(unitspecies,'(f18.2)',end=20) weightmolar(id_pos)
     156  !  write(*,*) weightmolar(id_pos)
     157    read(unitspecies,'(e18.1)',end=20) ohreact(id_pos)
     158  !  write(*,*) ohreact(id_pos)
     159    read(unitspecies,'(i18)',end=20) spec_ass(id_pos)
     160  !  write(*,*) spec_ass(id_pos)
     161    read(unitspecies,'(f18.2)',end=20) kao(id_pos)
     162  !       write(*,*) kao(id_pos)
     163
     164    pspecies=species(id_pos)
     165    pdecay=decay(id_pos)
     166    pweta=weta(id_pos)
     167    pwetb=wetb(id_pos)
     168    pweta_in=weta_in(id_pos)
     169    pwetb_in=wetb_in(id_pos)
     170    pwetc_in=wetc_in(id_pos)
     171    pwetd_in=wetd_in(id_pos)
     172    preldiff=reldiff(id_pos)
     173    phenry=henry(id_pos)
     174    pf0=f0(id_pos)
     175    pdensity=density(id_pos)
     176    pdquer=dquer(id_pos)
     177    pdsigma=dsigma(id_pos)
     178    pdryvel=dryvel(id_pos)
     179    pweightmolar=weightmolar(id_pos)
     180    pohreact=ohreact(id_pos)
     181    pspec_ass=spec_ass(id_pos)
     182    pkao=kao(id_pos)
    184183
    185184  else
    186185
    187     species(pos_spec)=pspecies
    188     decay(pos_spec)=pdecay
    189     weta(pos_spec)=pweta
    190     wetb(pos_spec)=pwetb
    191     weta_in(pos_spec)=pweta_in
    192     wetb_in(pos_spec)=pwetb_in
    193     wetc_in(pos_spec)=pwetc_in
    194     wetd_in(pos_spec)=pwetd_in
    195     reldiff(pos_spec)=preldiff
    196     henry(pos_spec)=phenry
    197     f0(pos_spec)=pf0
    198     density(pos_spec)=pdensity
    199     dquer(pos_spec)=pdquer
    200     dsigma(pos_spec)=pdsigma
    201     dryvel(pos_spec)=pdryvel
    202     weightmolar(pos_spec)=pweightmolar
    203     ohreact(pos_spec)=pohreact
    204     spec_ass(pos_spec)=pspec_ass
    205     kao(pos_spec)=pkao
    206 
    207   endif
    208 
    209   i=pos_spec
    210 
    211   if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
    212    if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
    213   endif
    214 
    215   if (spec_ass(pos_spec).gt.0) then
     186    species(id_pos)=pspecies
     187    decay(id_pos)=pdecay
     188    weta(id_pos)=pweta
     189    wetb(id_pos)=pwetb
     190    weta_in(id_pos)=pweta_in
     191    wetb_in(id_pos)=pwetb_in
     192    wetc_in(id_pos)=pwetc_in
     193    wetd_in(id_pos)=pwetd_in
     194    reldiff(id_pos)=preldiff
     195    henry(id_pos)=phenry
     196    f0(id_pos)=pf0
     197    density(id_pos)=pdensity
     198    dquer(id_pos)=pdquer
     199    dsigma(id_pos)=pdsigma
     200    dryvel(id_pos)=pdryvel
     201    weightmolar(id_pos)=pweightmolar
     202    ohreact(id_pos)=pohreact
     203    spec_ass(id_pos)=pspec_ass
     204    kao(id_pos)=pkao
     205
     206  endif
     207
     208  if ((weta(id_pos).gt.0).and.(henry(id_pos).le.0)) then
     209   if (dquer(id_pos).le.0) goto 996 ! no particle, no henry set
     210  endif
     211
     212  if (spec_ass(id_pos).gt.0) then
    216213    spec_found=.FALSE.
    217     do j=1,pos_spec-1
    218       if (spec_ass(pos_spec).eq.specnum(j)) then
    219         spec_ass(pos_spec)=j
     214    do j=1,id_pos-1
     215      if (spec_ass(id_pos).eq.specnum(j)) then
     216        spec_ass(id_pos)=j
    220217        spec_found=.TRUE.
    221218        ASSSPEC=.TRUE.
    … …  
    227224  endif
    228225
    229   if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
    230   if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
    231 
    232   if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
     226  if (dsigma(id_pos).eq.1.) dsigma(id_pos)=1.0001   ! avoid floating exception
     227  if (dsigma(id_pos).eq.0.) dsigma(id_pos)=1.0001   ! avoid floating exception
     228
     229  if ((reldiff(id_pos).gt.0.).and.(density(id_pos).gt.0.)) then
    233230    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
    234231    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
    … …  
    246243
    247244  do j=1,24           ! initialize everything to no variation
    248     area_hour(i,j)=1.
    249     point_hour(i,j)=1.
     245    area_hour(id_pos,j)=1.
     246    point_hour(id_pos,j)=1.
    250247  end do
    251248  do j=1,7
    252     area_dow(i,j)=1.
    253     point_dow(i,j)=1.
     249    area_dow(id_pos,j)=1.
     250    point_dow(id_pos,j)=1.
    254251  end do
    255252
    256   if (readerror.ne.0) then ! text format input
    257 
    258     read(unitspecies,*,end=22)
     253  if (ios.ne.0) then ! text format input
     254
     255    read(unitspecies,*,iostat=ios)
     256    if (ios .ne. 0) goto 22 ! ifort does not like err=
    259257    do j=1,24     ! 24 hours, starting with 0-1 local time
    260       read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
     258      read(unitspecies,*) ihour,area_hour(id_pos,j),point_hour(id_pos,j)
    261259    end do
    262260    read(unitspecies,*)
    263261    do j=1,7      ! 7 days of the week, starting with Monday
    264       read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
     262      read(unitspecies,*) idow,area_dow(id_pos,j),point_dow(id_pos,j)
    265263    end do
    266264

  • branches/petra/src/readwind.f90

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

  • branches/petra/src/readwind_nests.f90

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

  • branches/petra/src/timemanager.f90

    r30 r37  
    137137  !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    138138  write(*,46) float(itime)/3600,itime,numpart
    139   if (verbosity.gt.0) then
    140     write (*,*) 'timemanager> starting simulation'
    141     if (verbosity.gt.1) then
    142       CALL SYSTEM_CLOCK(count_clock)
    143       WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
    144     endif     
     139  call verboutput(verbosity,1,'timemanager> starting simulation')
     140  if (verbosity.gt.1) then
     141    CALL SYSTEM_CLOCK(count_clock)
     142    WRITE(*,*) 'timemanager> SYSTEM CLOCK', (count_clock-count_clock0) / real(count_rate)
    145143  endif
    146144

  • branches/petra/src/verttransform.f90

    r24 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    5050  ! added the variable cloud for use with scavenging - descr. in com_mod
    5151  !*****************************************************************************
     52  ! Petra Seibert, Feb 2015: Add quick fix from 2013
     53  !*****************************************************************************
    5254  !                                                                            *
    5355  ! Variables:                                                                 *
    5456  ! nx,ny,nz                        field dimensions in x,y and z direction    *
    55   ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
    5657  ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
    5758  ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
    … …  
    6970  implicit none
    7071
    71   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    72   integer :: rain_cloud_above,kz_inv
     72  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
     73  integer :: icloudtop
    7374  real :: f_qvsat,pressure
    74   real :: rh,lsp,convp
     75  real :: rh,lsp,convp,prec,rhmin
    7576  real :: rhoh(nuvzmax),pinmconv(nzmax)
    7677  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    … …  
    8586  real,parameter :: const=r_air/ga
    8687
    87   logical :: init = .true.
     88  logical :: init = .true., lconvectprec, lsearch
    8889
    8990
    … …  
    477478
    478479
    479   !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
    480   !   create a cloud and rainout/washout field, clouds occur where rh>80%
    481   !   total cloudheight is stored at level 0
     480  ! write (*,*) 'diagnosing clouds, n:',n,nymin1,nxmin1,nz
     481  jy_loop: &
    482482  do jy=0,nymin1
    483483    do ix=0,nxmin1
    484       rain_cloud_above=0
     484   
    485485      lsp=lsprec(ix,jy,1,n)
    486486      convp=convprec(ix,jy,1,n)
    487       cloudsh(ix,jy,n)=0
    488       do kz_inv=1,nz-1
    489          kz=nz-kz_inv+1
    490          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    491          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    492          clouds(ix,jy,kz,n)=0
    493          if (rh.gt.0.8) then ! in cloud
    494             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    495                rain_cloud_above=1
    496                cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    497                height(kz)-height(kz-1)
    498                if (lsp.ge.convp) then
    499                   clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    500                else
    501                   clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    502                endif
    503             else ! no precipitation
    504                   clouds(ix,jy,kz,n)=1 ! cloud
     487      prec=lsp+convp
     488      if (lsp .gt. convp) then !  prectype='lsp'
     489        lconvectprec = .false.
     490      else ! prectype='cp '
     491        lconvectprec = .true.
     492      endif
     493      icloudbot(ix,jy,n)=icmv
     494      icloudtop=icmv ! this is just a local variable
     495      rhmin=rhmininit ! just initialise in a way that cond is true
     496      lsearch=.true.
     497
     498      cloudsearch_loop: &
     499      do while (rhmin .ge. rhminx .and. lsearch)
     500        ! give up for < rhminx
     501
     502        kz_loop: &
     503        do kz_inv=1,nz-1
     504          kz=nz-kz_inv+1
     505          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     506          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     507          if (rh .gt. rhmin) then
     508            if (icloudbot(ix,jy,n) .eq. icmv) then
     509              icloudtop=nint(height(kz)) ! use int to save memory
    505510            endif
    506          else ! no cloud
    507             if (rain_cloud_above.eq.1) then ! scavenging
    508                if (lsp.ge.convp) then
    509                   clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    510                else
    511                   clouds(ix,jy,kz,n)=4 ! convp dominated washout
    512                endif
    513             endif
    514          endif
    515       end do
    516     end do
    517   end do
     511            icloudbot(ix,jy,n)=nint(height(kz))
     512          endif
     513        end do kz_loop
     514
     515        ! PS try to get a cloud thicker than 50 m
     516        ! PS if there is at least precmin mm/h
     517        if (prec .gt. precmin .and. &
     518          ( icloudbot(ix,jy,n) .eq. icmv .or. &
     519            icloudtop-icloudbot(ix,jy,n) .lt. 50)) then
     520          rhmin = rhmin - 0.05
     521        else
     522          lsearch = .false.
     523        endif
     524       
     525      enddo cloudsearch_loop
     526
     527      ! PS implement a rough fix for badly represented convection
     528      ! PS is based on looking at a limited set of comparison data
     529      if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. &
     530         icloudbot(ix,jy,n) .lt. icloudtopmin .and. prec .gt. precmin) then
     531        if (convp .lt. 0.1) then
     532          icloudbot(ix,jy,n) = icloudbot1
     533          icloudtop =          icloudtop1
     534        else
     535          icloudbot(ix,jy,n) = icloudbot2
     536          icloudtop =          icloudtop2
     537        endif
     538      endif
     539
     540      if (icloudtop .ne. icmv) then
     541        icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
     542      else
     543        icloudthck(ix,jy,n) = icmv ! no cloud found whatsoever
     544      endif
     545
     546      ! PS get rid of too thin clouds     
     547      if (icloudthck(ix,jy,n) .lt. 50) then
     548        icloudbot(ix,jy,n)=icmv
     549        icloudthck(ix,jy,n)=icmv
     550      endif
     551
     552    enddo
     553  enddo jy_loop
    518554
    519555end subroutine verttransform

  • branches/petra/src/verttransform_nests.f90

    r24 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    5252  ! add the variable cloud for use with scavenging - descr. in com_mod
    5353  !*****************************************************************************
     54  ! Petra Seibert, Feb 2015: Add quick fix from 2013
     55  !*****************************************************************************
    5456  !                                                                            *
    5557  ! Variables:                                                                 *
    … …  
    6971  implicit none
    7072
    71   integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp
    72   integer :: rain_cloud_above,kz_inv
    73   real :: f_qvsat,pressure,rh,lsp,convp
    74   real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
    75   real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
    76   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
     73  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
     74  integer :: icloudtop
     75  real :: f_qvsat,pressure
     76  real :: rh,lsp,convp,prec,rhmin
     77  real :: rhoh(nuvzmax),pinmconv(nzmax)
     78  real :: wzlev(nwzmax),uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
     79  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    7780  real :: dzdx,dzdy
    78   real :: dzdx1,dzdx2,dzdy1,dzdy2
     81  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    7982  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    8083  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    … …  
    8386  real,parameter :: const=r_air/ga
    8487
     88  logical :: lconvectprec, lsearch
    8589
    8690  ! Loop over all nests
    8791  !********************
    8892
     93  numbnest_loop: &
    8994  do l=1,numbnests
    9095
    … …  
    143148
    144149
    145   ! Levels, where u,v,t and q are given
    146   !************************************
     150  ! Levels where u,v,t and q are given
     151  !***********************************
    147152
    148153      uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
    … …  
    194199
    195200
    196   ! Levels, where w is given
    197   !*************************
     201  ! Levels where w is given
     202  !************************
    198203
    199204      wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
    … …  
    275280
    276281      end do
    277 
    278282    end do
    279283  end do
    280284
    281 
    282   !write (*,*) 'initializing nested cloudsn, n:',n
    283   !   create a cloud and rainout/washout field, cloudsn occur where rh>80%
     285  ! write (*,*) 'diagnosing nested clouds, n:',n
     286  jy_loop: &
    284287  do jy=0,nyn(l)-1
    285288    do ix=0,nxn(l)-1
    286       rain_cloud_above=0
     289
    287290      lsp=lsprecn(ix,jy,1,n,l)
    288291      convp=convprecn(ix,jy,1,n,l)
    289       cloudsnh(ix,jy,n,l)=0
    290       do kz_inv=1,nz-1
    291          kz=nz-kz_inv+1
    292          pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
    293          rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
    294          cloudsn(ix,jy,kz,n,l)=0
    295          if (rh.gt.0.8) then ! in cloud
    296             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
    297                rain_cloud_above=1
    298                cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    299                height(kz)-height(kz-1)
    300                if (lsp.ge.convp) then
    301                   cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
    302                else
    303                   cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
    304                endif
    305             else ! no precipitation
    306                   cloudsn(ix,jy,kz,n,l)=1 ! cloud
     292      prec=lsp+convp
     293      if (lsp .gt. convp) then !  prectype='lsp'
     294        lconvectprec = .false.
     295      else ! prectype='cp '
     296        lconvectprec = .true.
     297      endif
     298      icloudbotn(ix,jy,n,l)=icmv
     299      icloudtop=icmv ! this is just a local variable
     300      rhmin=rhmininit 
     301      lsearch=.true.
     302
     303      cloudsearch_loop: &
     304      do while (rhmin .ge. rhminx .and. lsearch) 
     305        ! give up for < rhminx
     306
     307        kz_loop: &
     308        do kz_inv=1,nz-1
     309           kz=nz-kz_inv+1
     310           pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
     311           rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
     312           if (rh .gt. rhmin) then
     313            if (icloudbotn(ix,jy,n,l) .eq. icmv) then
     314              icloudtop=nint(height(kz)) ! use int to save memory
    307315            endif
    308          else ! no cloud
    309             if (rain_cloud_above.eq.1) then ! scavenging
    310                if (lsp.ge.convp) then
    311                   cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
    312                else
    313                   cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
    314                endif
    315             endif
    316          endif
    317       end do
    318     end do
    319   end do
    320 
    321   end do
     316            icloudbotn(ix,jy,n,l)=nint(height(kz))
     317          endif
     318        end do kz_loop
     319
     320        ! PS try to get a cloud thicker than 50 m
     321        ! PS if there is at least precmin mm/h
     322        if (prec .gt. precmin .and. &
     323          ( icloudbotn(ix,jy,n,l) .eq. icmv .or. &
     324            icloudtop-icloudbotn(ix,jy,n,l) .lt. 50)) then
     325          rhmin = rhmin - 0.05
     326        else
     327          lsearch = .false.
     328        endif
     329       
     330      enddo cloudsearch_loop
     331     
     332      ! PS implement a rough fix for badly represented convection
     333      ! PS is based on looking at a limited set of comparison data
     334      if (lconvectprec .and. icloudtop .lt. icloudtopconvmin .or. &
     335         icloudbotn(ix,jy,n,l) .lt. icloudtopmin .and. prec .gt. precmin) then
     336        if (convp .lt. 0.1) then
     337          icloudbotn(ix,jy,n,l) = icloudbot1
     338          icloudtop =             icloudtop1
     339        else
     340          icloudbotn(ix,jy,n,l) = icloudbot2
     341          icloudtop =             icloudtop2
     342        endif
     343      endif
     344      if (icloudtop .ne. icmv) then
     345        icloudthckn(ix,jy,n,l) = icloudtop-icloudbotn(ix,jy,n,l)
     346      else
     347        icloudthckn(ix,jy,n,l) = icmv
     348      endif
     349      ! PS get rid of too thin clouds     
     350      if (icloudthck(ix,jy,n) .lt. 50) then
     351        icloudbotn(ix,jy,n,l)=icmv
     352        icloudthckn(ix,jy,n,l)=icmv
     353      endif
     354
     355    enddo
     356  enddo jy_loop
     357 
     358  enddo numbnest_loop
    322359
    323360end subroutine verttransform_nests

  • branches/petra/src/wetdepo.f90

    r24 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    3939  ! Code may not be correct for decay of deposition!                           *
    4040  !                                                                            *
     41  ! PS, 2/2015: implement wet depo quick fix from 2011/12 in this version
     42  !  it implements interpolated and improved clouds
     43  !  Also, certain deficiencies for thin clouds fixed also
     44  !  Pass itage to wetdepokernel
     45  !                                                                            *
    4146  !*****************************************************************************
    4247  !                                                                            *
    … …  
    7176
    7277  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
    74   integer :: ks, kp
     78  integer :: ngrid,itage,nage,kz,il,interp_time,n
     79  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
    7580  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    76   real :: clouds_h ! cloud height for the specific grid point
    77   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
     81  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav,wetscavold,precsub,f
    7882  real :: wetdeposit(maxspec),restmass
    7983  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    97101  !************************
    98102
     103  particle_loop: &
    99104  do jpart=1,numpart
    100105
    101106    if (itra1(jpart).eq.-999999999) goto 20
    102     if(ldirect.eq.1)then
     107    if (ldirect.eq.1) then
    103108      if (itra1(jpart).gt.itime) goto 20
    104109    else
    105110      if (itra1(jpart).lt.itime) goto 20
    106111    endif
     112   
    107113  ! Determine age class of the particle
    108114    itage=abs(itra1(jpart)-itramem(jpart))
    … …  
    123129        goto 23
    124130      endif
    125     end do
     131    enddo
    12613223  continue
    127133
    … …  
    131137
    132138    if (ngrid.gt.0) then
    133       xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
    134       ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
     139      xtn=real(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
     140      ytn=real(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
    135141      ix=int(xtn)
    136142      jy=int(ytn)
    … …  
    148154
    149155    if (ngrid.eq.0) then
    150       call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    151       1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    152       memtime(1),memtime(2),interp_time,lsp,convp,cc)
     156      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax,1,nx,ny,memind,&
     157      real(xtra1(jpart)),real(ytra1(jpart)), &
     158      1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop)
    153159    else
    154160      call interpol_rain_nests(lsprecn,convprecn,tccn, &
    155       nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    156       memtime(1),memtime(2),interp_time,lsp,convp,cc)
    157     endif
    158 
    159     if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
     161      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn, &
     162      1,memtime(1),memtime(2),interp_time,lsp,convp,cc,icbot,ictop)
     163    endif
     164
     165! PS 2012/2015: subtract a small value, eg 0.01 mm/h,
     166! to remove spurious precip; replaces previous code
     167        prec = lsp+convp
     168        precsub = 0.01
     169        if (prec .lt. precsub) then
     170          goto 20
     171        else
     172          f = (prec-precsub)/prec
     173          lsp = f*lsp
     174          convp = f*convp
     175        endif
     176
    160177  ! get the level were the actual particle is in
    161178    do il=2,nz
    162179      if (height(il).gt.ztra1(jpart)) then
    163         hz=il-1
     180        kz=il-1
    164181        goto 26
    165182      endif
    166     end do
     183    enddo
    16718426  continue
    168185
    169186    n=memind(2)
    170187    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    171     n=memind(1)
     188      n=memind(1)
    172189
    173190  ! if there is no precipitation or the particle is above the clouds no
    174191  ! scavenging is done
    175 
    176     if (ngrid.eq.0) then
    177       clouds_v=clouds(ix,jy,hz,n)
    178       clouds_h=cloudsh(ix,jy,n)
    179     else
    180       clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    181       clouds_h=cloudsnh(ix,jy,n,ngrid)
    182     endif
    183   !write(*,*) 'there is
    184   !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
    185     if (clouds_v.le.1) goto 20
    186   !write (*,*) 'there is scavenging'
     192  ! PS: part of 2011/2012/2015 fix, replaces previous code
     193 
     194    if (ztra1(jpart) .le. float(ictop)) then
     195      if (ztra1(jpart) .gt. float(icbot)) then
     196        indcloud = 2 ! in-cloud
     197      else
     198        indcloud = 1 ! below-cloud
     199      endif
     200    elseif (ictop .eq. icmv) then
     201      indcloud = 0 ! no cloud found, use old scheme
     202    else
     203      goto 20 ! above cloud
     204    endif
     205
    187206
    188207  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    228247  !**********************************************************
    229248
     249    species_loop: &
    230250    do ks=1,nspec      ! loop over species
     251   
    231252      wetdeposit(ks)=0.
    232253      wetscav=0.   
    233254
    234255  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
    235 
    236 
    237   !****i*******************
    238   ! BELOW CLOUD SCAVENGING
    239   !*********************** 
    240       if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime
    241  !       for aerosols and not highliy soluble substances weta=5E-6
    242         wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient
    243       endif !end below-cloud
    244 
    245 
    246   !********************
    247   ! IN CLOUD SCAVENGING
    248   !********************
    249       if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime
    250 
    251         if (ngrid.gt.0) then
    252           act_temp=ttn(ix,jy,hz,n,ngrid)
    253         else
    254           act_temp=tt(ix,jy,hz,n)
    255         endif
     256      if (weta(ks).gt.0.) then
     257        ! positive below-cloud coefficient from SPECIES file
     258       
     259        if (indcloud .eq. 1) then ! below-cloud scavenging
     260
     261         ! for aerosols and not highly soluble substances weta=5E-6
     262           wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient
     263
     264        elseif (indcloud .eq. 2) then ! in-cloud scavenging
     265
     266          if (ngrid.gt.0) then
     267            act_temp=ttn(ix,jy,kz,n,ngrid)
     268          else
     269            act_temp=tt(ix,jy,kz,n)
     270          endif
    256271
    257272! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    … …  
    259274! wetb_in=0.36 (default)
    260275! wetc_in=0.9 (default)
    261 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
     276! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0 - no scaling)
    262277        cl=weta_in(ks)*prec**wetb_in(ks)
    263278        if (dquer(ks).gt.0) then ! is particle
    264279          S_i=wetc_in(ks)/cl
    265280        else ! is gas
    266           cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    267           S_i=1/cle
     281          cle=(1.-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
     282          S_i=1./cle
    268283        endif
    269         wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
    270   !     write(*,*) 'in. wetscav:'
    271   !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
    272        
    273       endif ! end in-cloud
    274 
    275 
     284        wetscav=S_i*prec/3.6E6/wetd_in(ks)/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
     285! PS - prevent extrem values of wetscav:
     286        wetscavold = 2.e-5*prec**0.8     
     287        wetscav = min( 0.1*wetscav, wetscavold ) ! here we sacle the above wetscav. This is now twice, also in wetd_in.
     288
     289      else ! PS: no cloud diagnosed, old scheme,
     290
     291! PS    using with fixed a,b for simplicity, one may wish to change!!
     292        wetscav = 2.e-5*prec**0.8 ! was before: 1.e-4*prec**0.62
     293
     294      endif ! end regime cases
    276295     
    277   !**************************************************
    278   ! CALCULATE DEPOSITION
    279   !**************************************************
    280   !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    281   !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
    282 
    283       if (wetscav.gt.0.) then
    284         wetdeposit(ks)=xmass1(jpart,ks)* &
    285         (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
    286       else ! if no scavenging
    287         wetdeposit(ks)=0.
    288       endif
    289 
    290 
     296  ! calculate deposition
     297      wetdeposit(ks)=xmass1(jpart,ks)* &
     298        (1.-exp(-wetscav*abs(ltsample)))*grfraction
    291299      restmass = xmass1(jpart,ks)-wetdeposit(ks)
    292300      if (ioutputforeachrelease.eq.1) then
    … …  
    297305      if (restmass .gt. smallnum) then
    298306        xmass1(jpart,ks)=restmass
    299   !   depostatistic
    300   !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    301   !   depostatistic
     307    !   depostatistic:
     308   !!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    302309      else
    303310        xmass1(jpart,ks)=0.
    … …  
    305312  !   Correct deposited mass to the last time step when radioactive decay of
    306313  !   gridded deposited mass was calculated
    307       if (decay(ks).gt.0.) then
     314      if (decay(ks).gt.0.) &
    308315        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
    309       endif
    310 
    311 
    312 
    313     end do !all species
    314 
    315   ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
    316   ! Add the wet deposition to accumulated amount on output grid and nested output grid
    317   !*****************************************************************************
     316
     317    else ! weta(k) <= 0
     318
     319      wetdeposit(ks)=0.
     320   
     321    endif
     322
     323  enddo species_loop
     324
     325 ! Sabine Eckhardt, June 2008: write deposition only in forward runs
     326 ! add the wet deposition from this step to accumulated amount
     327 ! on output grid and nested output grid
    318328
    319329    if (ldirect.eq.1) then
    320       call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
    321       real(ytra1(jpart)),nage,kp)
     330      call wetdepokernel(nclass(jpart),wetdeposit, &
     331        real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
    322332      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
    323         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
    324     endif
    325 
    326 20  continue
    327   end do ! all particles
     333        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),itage,nage,kp)
     334    endif
     335
     33620  continue ! jump here for particles not to be treated
     337  enddo particle_loop
    328338
    329339end subroutine wetdepo

  • branches/petra/src/wetdepokernel.f90

    r4 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    2020!**********************************************************************
    2121
    22 subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
    23   !                          i      i    i i  i
     22subroutine wetdepokernel(nunc,deposit,x,y,itage,nage,kp)
     23  !                          i      i    i  i    i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    2626  !     Attribution of the deposition from an individual particle to the       *
    27   !     deposition fields using a uniform kernel with bandwidths dxout and dyout.*
     27  !     deposition fields using a uniform kernel with bandwidths dxout         *
     28  !     and dyout.                                                             *
    2829  !                                                                            *
    2930  !     Author: A. Stohl                                                       *
    3031  !                                                                            *
    3132  !     26 December 1996                                                       *
     33  !
     34  !  PS, 2/2015: do not use kernel for itage < 3 h
     35  !   same as for concentration in conccalc.f90
    3236  !                                                                            *
    3337  !*****************************************************************************
    … …  
    4852
    4953  real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
    50   integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp
     54  integer :: ix,jy,ixp,jyp,nunc,nage,ks,kp,itage
    5155
    5256  xl=(x*dx+xoutshift)/dxout
    … …  
    7478
    7579
    76   ! Determine mass fractions for four grid points
    77   !**********************************************
     80  if (itage .lt. itagekernmin) then
     81    ! no kernel, direct attribution to grid cell
     82   
     83    do ks=1,nspec
     84      if (ix.ge.0 .and. jy.ge.0 .and. &
     85          ix.le.numxgrid-1 .and. jy.le.numygrid-1) &
     86        wetgridunc(ix,jy,ks,kp,nunc,nage)= &
     87          wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
     88    enddo
     89       
     90  else
     91    ! Determine mass fractions for four grid points & distribute
    7892
    79   do ks=1,nspec
     93    do ks=1,nspec
    8094
    81   if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
    82        (jy.le.numygrid-1)) then
    83     w=wx*wy
    84       wetgridunc(ix,jy,ks,kp,nunc,nage)= &
    85            wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
     95      if (ix.ge.0 .and. jy.ge.0 .and. &
     96          ix.le.numxgrid-1 .and. jy.le.numygrid-1) then
     97        w=wx*wy
     98        wetgridunc(ix,jy,ks,kp,nunc,nage)= &
     99          wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
     100      endif
     101
     102      if ((ixp.ge.0).and.(jyp.ge.0).and.&
     103         (ixp.le.numxgrid-1).and. (jyp.le.numygrid-1)) then
     104        w=(1.-wx)*(1.-wy)
     105        wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
     106          wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
     107      endif
     108
     109      if ((ixp.ge.0).and.(jy.ge.0).and. &
     110          (ixp.le.numxgrid-1).and.(jy.le.numygrid-1)) then
     111        w=(1.-wx)*wy
     112        wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
     113           wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
     114      endif
     115
     116      if ((ix.ge.0).and.(jyp.ge.0).and. &
     117          (ix.le.numxgrid-1).and.(jyp.le.numygrid-1)) then
     118        w=wx*(1.-wy)
     119        wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
     120          wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
     121      endif
     122
     123    end do
     124
    86125  endif
    87126
    88   if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
    89        (jyp.le.numygrid-1)) then
    90     w=(1.-wx)*(1.-wy)
    91       wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
    92            wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    93   endif
    94 
    95   if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
    96        (jy.le.numygrid-1)) then
    97     w=(1.-wx)*wy
    98       wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
    99            wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    100   endif
    101 
    102   if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
    103        (jyp.le.numygrid-1)) then
    104     w=wx*(1.-wy)
    105       wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
    106            wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    107   endif
    108   end do
    109 
    110127end subroutine wetdepokernel

  • branches/petra/src/wetdepokernel_nest.f90

    r4 r37  
    11!**********************************************************************
    2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
     2! Copyright 1998-2015                                                 *
    33! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
    44! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
    … …  
    2020!**********************************************************************
    2121
    22 subroutine wetdepokernel_nest &
    23        (nunc,deposit,x,y,nage,kp)
    24   !        i      i    i i  i    i
     22subroutine wetdepokernel_nest (nunc,deposit,x,y,itage,nage,kp)
     23  !                             i      i    i i   i    i   i
    2524  !*****************************************************************************
    2625  !                                                                            *
    … …  
    3534  !      2 September 2004: Adaptation from wetdepokernel.                      *
    3635  !                                                                            *
     36  !
     37  !  PS, 2/2015: do not use kernel for itage < 3 h
     38  !   same as for concentration in conccalc.f90
    3739  !                                                                            *
    3840  !*****************************************************************************
    … …  
    5355
    5456  real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w
    55   integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage
     57  integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage,itage
    5658
    5759
    … …  
    8082  endif
    8183
     84  if (itage .lt. itagekernmin) then 
     85   ! no kernel, direct attribution to grid cell
     86   
     87    do ks=1,nspec
     88      if (ix.ge.0 .and. jy.ge.0 .and. &
     89          ix.le.numxgridn-1 .and. jy.le.numygridn-1) &
     90        wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
     91          wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)
     92    enddo
     93       
     94  else
     95    ! Determine mass fractions for four grid points & distribute
    8296
    83   ! Determine mass fractions for four grid points
    84   !**********************************************
     97    do ks=1,nspec
    8598
    86   do ks=1,nspec
     99      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
     100           (jy.le.numygridn-1)) then
     101        w=wx*wy
     102        wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
     103          wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
     104      endif
    87105
    88   if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
    89        (jy.le.numygridn-1)) then
    90     w=wx*wy
    91       wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
    92            wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
     106      if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &
     107           (jyp.le.numygridn-1)) then
     108        w=(1.-wx)*(1.-wy)
     109        wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &
     110          wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
     111      endif
     112
     113      if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &
     114           (jy.le.numygridn-1)) then
     115        w=(1.-wx)*wy
     116        wetgriduncn(ixp,jy,ks,kp,nunc,nage)= &
     117          wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
     118      endif
     119
     120      if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &
     121           (jyp.le.numygridn-1)) then
     122        w=wx*(1.-wy)
     123        wetgriduncn(ix,jyp,ks,kp,nunc,nage)= &
     124          wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
     125      endif
     126
     127    end do
     128
    93129  endif
    94 
    95   if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &
    96        (jyp.le.numygridn-1)) then
    97     w=(1.-wx)*(1.-wy)
    98       wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &
    99            wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    100   endif
    101 
    102   if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &
    103        (jy.le.numygridn-1)) then
    104     w=(1.-wx)*wy
    105       wetgriduncn(ixp,jy,ks,kp,nunc,nage)= &
    106            wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    107   endif
    108 
    109   if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &
    110        (jyp.le.numygridn-1)) then
    111     w=wx*(1.-wy)
    112       wetgriduncn(ix,jyp,ks,kp,nunc,nage)= &
    113            wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    114   endif
    115 
    116   end do
    117130end subroutine wetdepokernel_nest
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG