Changeset db712a8 in flexpart.git


Ignore:
Timestamp:
Mar 3, 2016, 12:34:56 PM (8 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
38b7917
Parents:
b0434e1
Message:

Completed handling of nested wind fields with cloud water (for wet deposition).

Location:
src
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • src/FLEXPART.f90

    rb0434e1 rdb712a8  
    315315  endif
    316316
    317   !! testing !!
    318   ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')
    319   ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')
    320 
    321 
    322317  ! Write basic information on the simulation to a file "header"
    323318  ! and open files that are to be kept open throughout the simulation

  • src/FLEXPART_MPI.f90

    rb0434e1 rdb712a8  
    348348  endif
    349349
    350   !! testing !!
    351   ! open(999,file=trim(path(1))//'OH_FIELDS/jscalar_50N.txt',action='write',status='new')
    352   ! open(998,file=trim(path(1))//'OH_FIELDS/jscalar_50S.txt',action='write',status='new')
    353 
    354 
    355350  ! Write basic information on the simulation to a file "header"
    356351  ! and open files that are to be kept open throughout the simulation

  • src/calcpar_nests.f90

    re200b7a rdb712a8  
    103103      ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
    104104           td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
     105      if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8
    105106
    106107  ! 2) Calculation of inverse Obukhov length scale
    … …  
    135136        subsceff=min(excessoron(ix,jy,l),hmixplus)
    136137      else
    137         subsceff=0
     138        subsceff=0.0
    138139      endif
    139140  !
    … …  
    149150
    150151      if (DRYDEP) then
    151         z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
    152         z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
     152        ! z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
     153        ! z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
     154        z0(7)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
    153155
    154156  ! Calculate relative humidity at surface
    … …  
    228230  end do
    229231
    230 
    231     call calcpv_nests(l,n,uuhn,vvhn,pvhn)
     232  ! Calculation of potential vorticity on 3-d grid
     233  !***********************************************
     234
     235  call calcpv_nests(l,n,uuhn,vvhn,pvhn)
    232236
    233237  end do

  • src/calcpv.f90

    r8a65cb0 rdb712a8  
    7070  enddo
    7171
    72 ! :DBG: halts with SIGFPE if compiling with -ffpe-trap
    73 !  where(ppml.le.0) ppml=10000.0
    7472!  ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa
    75   ppmk=(100000./ppml)**kappa
     73  ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa
    7674
    7775  do jy=0,nymin1

  • src/calcpv_nests.f90

    r4fbe7a5 rdb712a8  
    6868    enddo
    6969  enddo
    70   ppmk=(100000./ppml)**kappa
     70!  ppmk=(100000./ppml)**kappa
     71  ppmk(0:nxn(l)-1,0:nyn(l)-1,1:nuvz)=(100000./ppml(0:nxn(l)-1,0:nyn(l)-1,1:nuvz))**kappa
    7172
    7273  do jy=0,nyn(l)-1

  • src/com_mod.f90

    rb0434e1 rdb712a8  
    134134  logical :: readclouds
    135135!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
    136   logical :: sumclouds
     136  logical :: sumclouds
     137
     138  logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
    137139 
    138140
    … …  
    370372!ZHG Sep 2015 
    371373   real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
    372    real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
    373374   real :: clw4(0:nxmax-1,0:nymax-1,numwfmem) ! eso: =icloud_stats(:,:,4,:)
    374375
    … …  
    487488
    488489  real,allocatable,dimension(:,:,:,:,:) :: uun, vvn, wwn, ttn, qvn, pvn,&
    489        & rhon, drhodzn, tthn, qvhn
    490   integer,allocatable,dimension(:,:,:,:) :: cloudsnh
     490       & rhon, drhodzn, tthn, qvhn, clwcn, ciwcn, clwn, clwchn, ciwchn
     491  real,allocatable,dimension(:,:,:,:) :: clw4n
     492  integer,allocatable,dimension(:,:,:,:) :: cloudshn
    491493  integer(kind=1),allocatable,dimension(:,:,:,:,:) :: cloudsn
    492494
    … …  
    788790    allocate(qvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
    789791    allocate(pvn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
    790     allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,0:nzmax,numwfmem,numbnests))
    791     allocate(cloudsnh(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests))
     792    allocate(clwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     793    allocate(ciwcn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     794    allocate(clwn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     795
     796    allocate(cloudsn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
     797    allocate(cloudshn(0:nxmaxn-1,0:nymaxn-1,numwfmem,numbnests))
    792798    allocate(rhon(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
    793799    allocate(drhodzn(0:nxmaxn-1,0:nymaxn-1,nzmax,numwfmem,numbnests))
    794800    allocate(tthn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
    795801    allocate(qvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
    796 
     802    allocate(clwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
     803    allocate(ciwchn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,numwfmem,numbnests))
     804    allocate(clw4n(0:nxmax-1,0:nymax-1,numwfmem,numbnests))
     805
     806!    clw4n(:,:,:,:)=0.
     807    clwcn(:,:,:,:,:)=0.
     808    ciwcn(:,:,:,:,:)=0.
     809    clwchn(:,:,:,:,:)=0.
     810    ciwchn(:,:,:,:,:)=0.
     811    ! clwn(:,:,:,:,:)=0.
     812
     813    ! cloudsn(:,:,:,:,:)=0
     814    ! cloudshn(:,:,:,:)=0
    797815   
    798816  end subroutine com_mod_allocate_nests

  • src/concoutput.f90

    r6a678e3 rdb712a8  
    626626!*************************
    627627
    628   do ks=1,nspec
    629     do kp=1,maxpointspec_act
    630       do i=1,numreceptor
    631         creceptor(i,ks)=0.
    632       end do
    633       do jy=0,numygrid-1
    634         do ix=0,numxgrid-1
    635           do l=1,nclassunc
    636             do nage=1,nageclass
    637               do kz=1,numzgrid
    638                 gridunc(ix,jy,kz,ks,kp,l,nage)=0.
    639               end do
    640             end do
    641           end do
    642         end do
    643       end do
    644     end do
    645   end do
     628  ! do ks=1,nspec
     629  !   do kp=1,maxpointspec_act
     630  !     do i=1,numreceptor
     631  !       creceptor(i,ks)=0.
     632  !     end do
     633  !     do jy=0,numygrid-1
     634  !       do ix=0,numxgrid-1
     635  !         do l=1,nclassunc
     636  !           do nage=1,nageclass
     637  !             do kz=1,numzgrid
     638  !               gridunc(ix,jy,kz,ks,kp,l,nage)=0.
     639  !             end do
     640  !           end do
     641  !         end do
     642  !       end do
     643  !     end do
     644  !   end do
     645  ! end do
     646  creceptor(:,:)=0.
     647  gridunc(:,:,:,:,:,:,:)=0.
    646648
    647649

  • src/concoutput_nest.f90

    r6a678e3 rdb712a8  
    545545  !*************************
    546546
    547   do ks=1,nspec
    548   do kp=1,maxpointspec_act
    549     do i=1,numreceptor
    550       creceptor(i,ks)=0.
    551     end do
    552     do jy=0,numygridn-1
    553       do ix=0,numxgridn-1
    554         do l=1,nclassunc
    555           do nage=1,nageclass
    556             do kz=1,numzgrid
    557               griduncn(ix,jy,kz,ks,kp,l,nage)=0.
    558             end do
    559           end do
    560         end do
    561       end do
    562     end do
    563   end do
    564   end do
     547  ! do ks=1,nspec
     548  ! do kp=1,maxpointspec_act
     549  !   do i=1,numreceptor
     550  !     creceptor(i,ks)=0.
     551  !   end do
     552  !   do jy=0,numygridn-1
     553  !     do ix=0,numxgridn-1
     554  !       do l=1,nclassunc
     555  !         do nage=1,nageclass
     556  !           do kz=1,numzgrid
     557  !             griduncn(ix,jy,kz,ks,kp,l,nage)=0.
     558  !           end do
     559  !         end do
     560  !       end do
     561  !     end do
     562  !   end do
     563  ! end do
     564  ! end do
     565  creceptor(:,:)=0.
     566  griduncn(:,:,:,:,:,:,:)=0.
    565567
    566568

  • src/ecmwf_mod.f90

    rb0434e1 rdb712a8  
    4141  !*********************************************
    4242
    43   integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new
     43!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new
     44  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
    4445  integer,parameter :: nxshift=359
     46!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !test BUG
     47!  integer,parameter :: nxshift=0
    4548
    4649!
    … …  
    4952  !*********************************************
    5053
    51   integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181 !ECMWF
     54  integer,parameter :: maxnests=1,nxmaxn=361,nymaxn=181
     55!  integer,parameter :: maxnests=1,nxmaxn=86,nymaxn=31
    5256
    5357

  • src/gridcheck.f90

    rb0434e1 rdb712a8  
    194194  elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
    195195    isec1(6)=246         ! indicatorOfParameter
    196     ! readclouds=.true.
    197     ! sumclouds=.false.
    198196  elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
    199197    isec1(6)=247         ! indicatorOfParameter
    … …  
    202200  elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
    203201    isec1(6)=201031      ! indicatorOfParameter
    204     ! readclouds=.true.
    205     ! sumclouds=.true.
    206202  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    207203    isec1(6)=134         ! indicatorOfParameter
    … …  
    366362  endif ! gotGrid
    367363
     364  if (nx.gt.nxmax) then
     365   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
     366    write(*,*) 'Reduce resolution of wind fields.'
     367    write(*,*) 'Or change parameter settings in file par_mod.'
     368    write(*,*) nx,nxmax
     369    stop
     370  endif
     371
     372  if (ny.gt.nymax) then
     373   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
     374    write(*,*) 'Reduce resolution of wind fields.'
     375    write(*,*) 'Or change parameter settings in file par_mod.'
     376    write(*,*) ny,nymax
     377    stop
     378  endif
     379
    368380  k=isec1(8)
    369381  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
    … …  
    410422  nwz =iwmax
    411423  if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
    412 
    413   if (nx.gt.nxmax) then
    414    write(*,*) 'FLEXPART error: Too many grid points in x direction.'
    415     write(*,*) 'Reduce resolution of wind fields.'
    416     write(*,*) 'Or change parameter settings in file par_mod.'
    417     write(*,*) nx,nxmax
    418     stop
    419   endif
    420 
    421   if (ny.gt.nymax) then
    422    write(*,*) 'FLEXPART error: Too many grid points in y direction.'
    423     write(*,*) 'Reduce resolution of wind fields.'
    424     write(*,*) 'Or change parameter settings in file par_mod.'
    425     write(*,*) ny,nymax
    426     stop
    427   endif
    428424
    429425  if (nuvz+1.gt.nuvzmax) then

  • src/gridcheck_nests.f90

    rb0434e1 rdb712a8  
    172172  elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
    173173    isec1(6)=246         ! indicatorOfParameter
    174     ! readclouds=.true.
    175     ! sumclouds=.false.
    176174  elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
    177175    isec1(6)=247         ! indicatorOfParameter
    … …  
    180178  elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
    181179    isec1(6)=201031      ! indicatorOfParameter
    182     ! readclouds=.true.
    183     ! sumclouds=.true.
    184180  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    185181    isec1(6)=134         ! indicatorOfParameter
    186182  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    187183    isec1(6)=135         ! indicatorOfParameter
    188   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90
    189     isec1(6)=135         ! indicatorOfParameter    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " "
     184  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridcheck.f90
     185    isec1(6)=135         ! indicatorOfParameter    !
    190186  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    191187    isec1(6)=151         ! indicatorOfParameter
    … …  
    258254  endif ! ifield
    259255
     256  if (nxn(l).gt.nxmaxn) then
     257  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
     258  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
     259  write(*,*) 'for nesting level ',l
     260  write(*,*) 'Or change parameter settings in file par_mod.'
     261  write(*,*) nxn(l),nxmaxn
     262  stop
     263  endif
     264
     265  if (nyn(l).gt.nymaxn) then
     266  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
     267  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
     268  write(*,*) 'for nesting level ',l
     269  write(*,*) 'Or change parameter settings in file par_mod.'
     270  write(*,*) nyn(l),nymaxn
     271  stop
     272  endif
     273
    260274  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    261275  if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!!
    … …  
    333347  if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
    334348
    335   if (nxn(l).gt.nxmaxn) then
    336   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
    337   write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
    338   write(*,*) 'for nesting level ',l
    339   write(*,*) 'Or change parameter settings in file par_mod.'
    340   write(*,*) nxn(l),nxmaxn
    341   stop
    342   endif
    343 
    344   if (nyn(l).gt.nymaxn) then
    345   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
    346   write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
    347   write(*,*) 'for nesting level ',l
    348   write(*,*) 'Or change parameter settings in file par_mod.'
    349   write(*,*) nyn(l),nymaxn
    350   stop
    351   endif
    352 
    353349  if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
    354350  write(*,*) 'FLEXPART error: Nested wind fields have too many'// &

  • src/initialize.f90

    r8a65cb0 rdb712a8  
    9999  jyp=jy+1
    100100
    101   h=max(hmix(ix ,jy ,1,memind(1)), &
     101  h=max(hmix(ix ,jy,1,memind(1)), &
    102102       hmix(ixp,jy ,1,memind(1)), &
    103103       hmix(ix ,jyp,1,memind(1)), &
    … …  
    157157    wp=rannumb(nrand+2)
    158158    if (.not.turbswitch) then     ! modified by mc
     159      wp=wp*sigw
     160    else if (cblflag.eq.1) then   ! modified by mc
     161      if(-h/ol.gt.5) then
     162!if (ol.lt.0.) then
     163!if (ol.gt.0.) then !by mc : only for test correct is lt.0
     164        call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol)
     165      else
    159166        wp=wp*sigw
    160     else if (cblflag.eq.1) then   ! modified by mc
    161         if(-h/ol.gt.5) then
    162         !if (ol.lt.0.) then
    163         !if (ol.gt.0.) then !by mc : only for test correct is lt.0
    164             call initialize_cbl_vel(idummy,zt,ust,wst,h,sigw,wp,ol)
    165         else
    166             wp=wp*sigw
    167         end if           
     167      end if
    168168    end if
    169169

  • src/interpol_rain_nests.f90

    r8a65cb0 rdb712a8  
    8888  !***********************************************************
    8989
    90   if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
    91        xt=real(nxn(ngrid)-1)-0.0001
    92   if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
    93        yt=real(nyn(ngrid)-1)-0.0001
     90  ! if (xt.ge.(real(nxn(ngrid)-1)-0.0001)) &
     91  !      xt=real(nxn(ngrid)-1)-0.0001
     92  ! if (yt.ge.(real(nyn(ngrid)-1)-0.0001)) &
     93  !      yt=real(nyn(ngrid)-1)-0.0001
     94
     95! ESO make it consistent with interpol_rain
     96  if (xt.ge.(real(nxn(ngrid)-1))) xt=real(nxn(ngrid)-1)-0.00001
     97  if (yt.ge.(real(nyn(ngrid)-1))) yt=real(nyn(ngrid)-1)-0.00001
    9498
    9599
    … …  
    105109  ix=int(xt)
    106110  jy=int(yt)
     111
    107112  ixp=ix+1
    108113  jyp=jy+1

  • src/makefile

    rb0434e1 rdb712a8  
    388388readdepo.o: com_mod.o par_mod.o
    389389readlanduse.o: com_mod.o par_mod.o
    390 readlanduse_int1.o: com_mod.o par_mod.o
     390#readlanduse_int1.o: com_mod.o par_mod.o
    391391readOHfield.o: com_mod.o oh_mod.o par_mod.o
    392392readoutgrid.o: com_mod.o outg_mod.o par_mod.o

  • src/mpi_mod.f90

    r6a678e3 rdb712a8  
    860860!   step
    861861!
    862 ! TODO
    863 !   GFS version
    864862!
    865863!*******************************************************************************
    … …  
    922920
    923921! The non-reader processes need to know if cloud water were read.
    924 ! TODO: only at first step or always?
    925922    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
    926923    if (mp_ierr /= 0) goto 600
    … …  
    10171014    if (mp_ierr /= 0) goto 600
    10181015
    1019     call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    1020     if (mp_ierr /= 0) goto 600
    1021 
    10221016    if (mp_measure_time) call mpif_mtime('commtime',1)
    10231017
    … …  
    10451039!   step
    10461040!
    1047 ! TODO
    1048 !   Transfer cloud water/ice if and when available for nested
    10491041!
    10501042!***********************************************************************
    … …  
    10601052
    10611053! Common array sizes used for communications
    1062     integer :: d3_size1 = nxmaxn*nymaxn*nzmax*maxnests
    1063     integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax*maxnests
    1064     integer :: d2_size1 = nxmaxn*nymaxn*maxnests
    1065     integer :: d2_size2 = nxmaxn*nymaxn*maxspec*maxnests
    1066     integer :: d2_size3 = nxmaxn*nymaxn*maxnests
     1054    integer :: d3_size1 = nxmaxn*nymaxn*nzmax
     1055    integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax
     1056    integer :: d2_size1 = nxmaxn*nymaxn
     1057    integer :: d2_size2 = nxmaxn*nymaxn*maxspec
     1058    integer :: d2_size3 = nxmaxn*nymaxn
    10671059
    10681060    integer :: d3s1,d3s2,d2s1,d2s2
    … …  
    11061098!**********************************************************************
    11071099
     1100! The non-reader processes need to know if cloud water were read.
     1101    call MPI_Bcast(readclouds_nest,maxnests,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
     1102    if (mp_ierr /= 0) goto 600
     1103
    11081104! Static fields/variables sent only at startup
    11091105    if (first_call) then
    … …  
    11201116
    11211117! MPI prefers contiguous arrays for sending (else a buffer is created),
    1122 ! hence the loop
    1123 
     1118! hence the loop over nests
     1119!**********************************************************************
    11241120    do i=1, numbnests
    11251121! 3D fields
    … …  
    11451141      if (mp_ierr /= 0) goto 600
    11461142      call MPI_Bcast(cloudsn(:,:,:,li:ui,i),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
    1147       if (mp_ierr /= 0) goto 600
     1143      if (mp_ierr /= 0) goto 600
     1144
     1145! cloud water/ice:
     1146    if (readclouds_nest(i)) then
     1147      ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
     1148      ! if (mp_ierr /= 0) goto 600
     1149      call MPI_Bcast(clw4n(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
     1150      if (mp_ierr /= 0) goto 600
     1151    end if
    11481152
    11491153! 2D fields
    1150       call MPI_Bcast(cloudsnh(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
     1154      call MPI_Bcast(cloudshn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    11511155      if (mp_ierr /= 0) goto 600
    11521156      call MPI_Bcast(vdepn(:,:,:,li:ui,i),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    … …  
    12061210!   mind    -- index where to place new fields
    12071211!
    1208 ! TODO
    12091212!
    12101213!*******************************************************************************
    … …  
    13921395!   memstat -- input, used to resolve windfield index being received
    13931396!
    1394 ! TODO
    13951397!
    13961398!*******************************************************************************
    … …  
    20192021! In the implementation with 3 fields, the processes may have posted
    20202022! MPI_Irecv requests that should be cancelled here
    2021 !! TODO:
    20222023! if (.not.lmp_sync) then
    20232024!   r=mp_pid*nvar_async
    … …  
    21052106!   
    21062107!
    2107 ! TODO
    21082108!
    21092109!*******************************************************************************
    … …  
    21322132    uu(:,:,:,li:ui) = 10.0
    21332133    vv(:,:,:,li:ui) = 0.0
    2134     uupol(:,:,:,li:ui) = 10.0 ! TODO check if ok
     2134    uupol(:,:,:,li:ui) = 10.0
    21352135    vvpol(:,:,:,li:ui)=0.0
    21362136    ww(:,:,:,li:ui)=0.
    … …  
    21642164    tropopause(:,:,:,li:ui)=10000.
    21652165    oli(:,:,:,li:ui)=0.01
    2166     z0=1.0
    21672166 
    21682167  end subroutine set_fields_synthetic

  • src/netcdf_output_mod.f90

    r6a678e3 rdb712a8  
    7979  implicit none
    8080
     81  private
     82
     83  public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
     84       &concoutput_nest_netcdf,concoutput_surf_netcdf
     85
    8186!  include 'netcdf.inc'
    8287
    … …  
    97102  real,parameter :: eps=nxmax/3.e5
    98103
    99   private:: writemetadata, output_units, nf90_err
     104!  private:: writemetadata, output_units, nf90_err
    100105
    101106  ! switch output of release point information on/off
    … …  
    628633
    629634  ! volume
    630   call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:)))
     635  if (lnest) then
     636    call nf90_err(nf90_put_var(ncid, volID, volumen(:,:,:)))
     637  else
     638    call nf90_err(nf90_put_var(ncid, volID, volume(:,:,:)))
     639  end if
    631640
    632641  ! area
    633   call nf90_err(nf90_put_var(ncid, areaID, area(:,:)))
     642  if (lnest) then
     643    call nf90_err(nf90_put_var(ncid, areaID, arean(:,:)))
     644  else
     645    call nf90_err(nf90_put_var(ncid, areaID, area(:,:)))
     646  end if
    634647
    635648  if (write_releases.eqv..true.) then

  • src/par_mod.f90

    rb0434e1 rdb712a8  
    139139  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new
    140140  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
    141   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    142141  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    143   !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
    144142
    145143!  integer,parameter :: nxshift=359 ! for ECMWF
    … …  
    216214  !**************************************************
    217215
    218   integer,parameter :: maxpart=40000000
     216  integer,parameter :: maxpart=400000
    219217  integer,parameter :: maxspec=6
    220   integer,parameter :: minmass=0.0 !0.0001
     218  real,parameter :: minmass=0.0 !0.0001
    221219
    222220  ! maxpart                 Maximum number of particles

  • src/readavailable.f90

    r5f9d14a rdb712a8  
    283283  stop
    284284
    285 999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
     285999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### '
    286286  write(*,'(a)') '     '//path(4)(1:length(4))
    287287  write(*,*) ' #### CANNOT BE OPENED           #### '

  • src/readwind_nests.f90

    rb0434e1 rdb712a8  
    347347! -add check for if one of clwc/ciwc missing (error),
    348348!    also if all 3 cw fields present, use qc and disregard the others
    349 ! -use same flags readclouds/sumclouds as in mother grid? this assumes
    350 !    that both the nested and mother grids contain CW in same format
    351349        if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
    352           clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    353           ! readclouds=.true.
    354           ! sumclouds=.false.
     350          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
     351          readclouds_nest(l)=.true.
     352          sumclouds_nest(l)=.false.
    355353        endif
    356354        if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
    357           ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     355          ciwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    358356        endif
    359357!ZHG end
    360358!ESO read qc (=clwc+ciwc)
    361359        if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
    362           clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    363           ! readclouds=.true.
    364           ! sumclouds=.true.
     360          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
     361          readclouds_nest(l)=.true.
     362          sumclouds_nest(l)=.true.
    365363        endif
    366364

  • src/timemanager_mpi.f90

    r6a678e3 rdb712a8  
    237237    if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then
    238238      call mpif_gf_send_vars(memstat)
    239       call mpif_gf_send_vars_nest(memstat)
     239      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
    240240! Version 2  (lmp_sync=.false., see below) is also used whenever 2 new fields are
    241241! read (as at first time step), in which case async send/recv is impossible.
    242242    else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
    243243      call mpif_gf_send_vars(memstat)
    244       call mpif_gf_send_vars_nest(memstat)
     244      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
    245245    end if
    246246
    … …  
    260260! Issued at start of each new field period.
    261261      if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
    262 ! TODO: z0(7) changes with time, so should be dimension (numclass,2) to
    263 ! allow transfer of the future value in the background
    264         call MPI_Bcast(z0,numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    265262        call mpif_gf_request
    266263      end if

  • src/verttransform.f90

    rb0434e1 rdb712a8  
    6767  use com_mod
    6868  use cmapf_mod, only: cc2gll
    69 ! eso TODO:
    70 ! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below)
    71 ! as this routine should not be dependent on MPI
    7269!  use mpi_mod
    73 ! :TODO
    7470
    7571  implicit none
    7672
    77   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    78   integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1)
    79   real :: f_qvsat,pressure,cosf(0:nymax-1)
    80   real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec
    81   real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax)
    82   real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax)
    83   real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1)
    84   real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi
     73  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
     74  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh
     75
     76  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
     77  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
     78  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
     79  real,dimension(0:nymax-1) :: cosf
     80
     81  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
     82
     83  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
     84  real :: f_qvsat,pressure,rh,lsp,convp,prec
     85  real :: ew,dz1,dz2,dz
    8586  real :: xlon,ylat,xlonr,dzdx,dzdy
    8687  real :: dzdx1,dzdx2,dzdy1,dzdy2
    8788  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    88   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    89   real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    90   real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
    91   real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    92   real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)
    9389  real,parameter :: const=r_air/ga
    94 !  integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2
    95   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
     90  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
    9691
    9792  logical :: init = .true.
    9893
    9994  !ZHG SEP 2014 tests 
    100   integer :: cloud_ver,cloud_min, cloud_max
    101   integer ::teller(5), convpteller=0, lspteller=0
    102   real :: cloud_col_wat, cloud_water
     95  ! integer :: cloud_ver,cloud_min, cloud_max
     96  ! integer ::teller(5), convpteller=0, lspteller=0
     97  ! real :: cloud_col_wat, cloud_water
    10398  !ZHG 2015 temporary variables for testing
    104   real :: rcw(0:nxmax-1,0:nymax-1)
    105   real :: rpc(0:nxmax-1,0:nymax-1)
     99  ! real :: rcw(0:nxmax-1,0:nymax-1)
     100  ! real :: rpc(0:nxmax-1,0:nymax-1)
    106101  ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
    107102  ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
    108   CHARACTER(LEN=3)  :: aspec
    109   integer :: virr=0
     103  ! CHARACTER(LEN=3)  :: aspec
     104  ! integer :: virr=0
    110105  real :: tot_cloud_h
     106  real :: dbg_height(nzmax)
    111107!ZHG
    112108
    … …  
    143139
    144140
    145     tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
     141    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ &
    146142         ps(ixm,jym,1,n))
    147143    pold(ixm,jym)=ps(ixm,jym,1,n)
    … …  
    182178    init=.false.
    183179
     180    dbg_height = height
     181
    184182  endif
    185183
    … …  
    191189  do jy=0,nymin1
    192190    do ix=0,nxmin1
    193       tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
     191      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
    194192           ps(ix,jy,1,n))
    195193    enddo
    … …  
    250248  qv(:,:,1,n)=qvh(:,:,1,n)
    251249!hg adding the cloud water
    252   clwc(:,:,1,n)=clwch(:,:,1,n)
    253   if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)   
     250  if (readclouds) then
     251    clwc(:,:,1,n)=clwch(:,:,1,n)
     252    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
     253  end if
    254254!hg
    255255  pv(:,:,1,n)=pvh(:,:,1)
    256256  rho(:,:,1,n)=rhoh(:,:,1)
     257
    257258  uu(:,:,nz,n)=uuh(:,:,nuvz)
    258259  vv(:,:,nz,n)=vvh(:,:,nuvz)
    259260  tt(:,:,nz,n)=tth(:,:,nuvz,n)
    260261  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
    261 
    262262!hg adding the cloud water
    263   clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
    264   if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
     263  if (readclouds) then
     264    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
     265    if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
     266  end if
    265267!hg
    266268  pv(:,:,nz,n)=pvh(:,:,nuvz)
    … …  
    279281          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
    280282!hg adding the cloud water
    281           clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
    282           if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     283          if (readclouds) then
     284            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
     285            if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     286          end if
    283287!hg
    284288          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
    … …  
    309313               +qvh(ix,jy,kz,n)*dz1)/dz
    310314!hg adding the cloud water
    311           clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
    312           if (.not.sumclouds) &
    313                &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
     315          if (readclouds) then
     316            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
     317            if (.not.sumclouds) &
     318                 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
     319          end if
    314320!hg
    315321          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    … …  
    601607            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
    602608!ZHG 2015 extra for testing
    603             clh(ix,jy,kz,n)=height(kz+1)-height(kz)
     609!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
    604610            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
    605611            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]

  • src/verttransform_nests.f90

    rb0434e1 rdb712a8  
    2121
    2222subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
    23   !                            i   i    i    i   i
    24   !*****************************************************************************
    25   !                                                                            *
    26   !     This subroutine transforms temperature, dew point temperature and      *
    27   !     wind components from eta to meter coordinates.                         *
    28   !     The vertical wind component is transformed from Pa/s to m/s using      *
    29   !     the conversion factor pinmconv.                                        *
    30   !     In addition, this routine calculates vertical density gradients        *
    31   !     needed for the parameterization of the turbulent velocities.           *
    32   !     It is similar to verttransform, but makes the transformations for      *
    33   !     the nested grids.                                                      *
    34   !                                                                            *
    35   !     Author: A. Stohl, G. Wotawa                                            *
    36   !                                                                            *
    37   !     12 August 1996                                                         *
    38   !     Update: 16 January 1998                                                *
    39   !                                                                            *
    40   !     Major update: 17 February 1999                                         *
    41   !     by G. Wotawa                                                           *
    42   !                                                                            *
    43   !     - Vertical levels for u, v and w are put together                      *
    44   !     - Slope correction for vertical velocity: Modification of calculation  *
    45   !       procedure                                                            *
    46   !                                                                            *
    47   !*****************************************************************************
    48   !  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
    49   !   Variables tthn and qvhn (on eta coordinates) from common block
    50   !*****************************************************************************
    51   ! Sabine Eckhardt, March 2007
    52   ! add the variable cloud for use with scavenging - descr. in com_mod
    53   !*****************************************************************************
    54   !                                                                            *
    55   ! Variables:                                                                 *
    56   ! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
    57   ! uun                             wind components in x-direction [m/s]       *
    58   ! vvn                             wind components in y-direction [m/s]       *
    59   ! wwn                             wind components in z-direction [deltaeta/s]*
    60   ! ttn                             temperature [K]                            *
    61   ! pvn                             potential vorticity (pvu)                  *
    62   ! psn                             surface pressure [Pa]                      *
    63   !                                                                            *
    64   !*****************************************************************************
     23!                            i   i    i    i   i
     24!*****************************************************************************
     25!                                                                            *
     26!     This subroutine transforms temperature, dew point temperature and      *
     27!     wind components from eta to meter coordinates.                         *
     28!     The vertical wind component is transformed from Pa/s to m/s using      *
     29!     the conversion factor pinmconv.                                        *
     30!     In addition, this routine calculates vertical density gradients        *
     31!     needed for the parameterization of the turbulent velocities.           *
     32!     It is similar to verttransform, but makes the transformations for      *
     33!     the nested grids.                                                      *
     34!                                                                            *
     35!     Author: A. Stohl, G. Wotawa                                            *
     36!                                                                            *
     37!     12 August 1996                                                         *
     38!     Update: 16 January 1998                                                *
     39!                                                                            *
     40!     Major update: 17 February 1999                                         *
     41!     by G. Wotawa                                                           *
     42!                                                                            *
     43!     - Vertical levels for u, v and w are put together                      *
     44!     - Slope correction for vertical velocity: Modification of calculation  *
     45!       procedure                                                            *
     46!                                                                            *
     47!*****************************************************************************
     48!  Changes, Bernd C. Krueger, Feb. 2001:       (marked "C-cv")
     49!   Variables tthn and qvhn (on eta coordinates) from common block
     50!*****************************************************************************
     51! Sabine Eckhardt, March 2007
     52! add the variable cloud for use with scavenging - descr. in com_mod
     53!*****************************************************************************
     54! ESO, 2016
     55! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than
     56!  the actual field dimensions
     57!*****************************************************************************
     58!                                                                            *
     59! Variables:                                                                 *
     60! nxn,nyn,nuvz,nwz                field dimensions in x,y and z direction    *
     61! uun                             wind components in x-direction [m/s]       *
     62! vvn                             wind components in y-direction [m/s]       *
     63! wwn                             wind components in z-direction [deltaeta/s]*
     64! ttn                             temperature [K]                            *
     65! pvn                             potential vorticity (pvu)                  *
     66! psn                             surface pressure [Pa]                      *
     67!                                                                            *
     68!*****************************************************************************
    6569
    6670  use par_mod
    … …  
    6973  implicit none
    7074
    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
     75  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: uuhn,vvhn,pvhn
     76  real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
     77
     78  real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
     79  real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
     80  real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
     81  real,dimension(0:nymaxn-1) :: cosf
     82
     83  integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
     84
     85  integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
     86  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
     87
     88  real :: ew,dz1,dz2,dz
    7789  real :: dzdx,dzdy
    7890  real :: dzdx1,dzdx2,dzdy1,dzdy2
    79   real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    80   real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    81   real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    82   real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
    8391  real,parameter :: const=r_air/ga
    84 
    85 
    86   ! Loop over all nests
    87   !********************
     92  real :: tot_cloud_h
     93  integer :: nxm1, nym1
     94
     95!  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
     96
     97! Loop over all nests
     98!********************
    8899
    89100  do l=1,numbnests
    90 
    91   ! Loop over the whole grid
    92   !*************************
    93 
    94   do jy=0,nyn(l)-1
    95     do ix=0,nxn(l)-1
    96 
    97       tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    98       psn(ix,jy,1,n,l))
    99       pold=psn(ix,jy,1,n,l)
    100       uvwzlev(ix,jy,1)=0.
    101       wzlev(1)=0.
    102       rhoh(1)=pold/(r_air*tvold)
    103 
    104 
    105   ! Compute heights of eta levels
    106   !******************************
    107 
    108       do kz=2,nuvz
    109         pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)
    110         tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
    111         rhoh(kz)=pint/(r_air*tv)
    112 
    113         if (abs(tv-tvold).gt.0.2) then
    114           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    115           (tv-tvold)/log(tv/tvold)
    116         else
    117           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    118         endif
    119 
    120         tvold=tv
    121         pold=pint
     101    nxm1=nxn(l)-1
     102    nym1=nyn(l)-1
     103
     104! Loop over the whole grid
     105!*************************
     106
     107    do jy=0,nyn(l)-1
     108      do ix=0,nxn(l)-1
     109        tvold(ix,jy)=tt2n(ix,jy,1,n,l)*(1.+0.378*ew*(td2n(ix,jy,1,n,l))/ &
     110             psn(ix,jy,1,n,l))
    122111      end do
    123 
    124 
    125       do kz=2,nwz-1
    126         wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
    127       end do
    128       wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    129 
    130   ! pinmconv=(h2-h1)/(p2-p1)
    131 
    132       pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    133       ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    134       (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    135       do kz=2,nz-1
    136         pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    137         ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
    138         (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
    139       end do
    140       pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    141       ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    142       (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    143 
    144 
    145   ! Levels, where u,v,t and q are given
    146   !************************************
    147 
    148       uun(ix,jy,1,n,l)=uuhn(ix,jy,1,l)
    149       vvn(ix,jy,1,n,l)=vvhn(ix,jy,1,l)
    150       ttn(ix,jy,1,n,l)=tthn(ix,jy,1,n,l)
    151       qvn(ix,jy,1,n,l)=qvhn(ix,jy,1,n,l)
    152       pvn(ix,jy,1,n,l)=pvhn(ix,jy,1,l)
    153       rhon(ix,jy,1,n,l)=rhoh(1)
    154       uun(ix,jy,nz,n,l)=uuhn(ix,jy,nuvz,l)
    155       vvn(ix,jy,nz,n,l)=vvhn(ix,jy,nuvz,l)
    156       ttn(ix,jy,nz,n,l)=tthn(ix,jy,nuvz,n,l)
    157       qvn(ix,jy,nz,n,l)=qvhn(ix,jy,nuvz,n,l)
    158       pvn(ix,jy,nz,n,l)=pvhn(ix,jy,nuvz,l)
    159       rhon(ix,jy,nz,n,l)=rhoh(nuvz)
    160       kmin=2
    161       do iz=2,nz-1
    162         do kz=kmin,nuvz
    163           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     112    end do
     113
     114    pold(0:nxm1,0:nym1)=psn(0:nxm1,0:nym1,1,n,l)
     115    uvzlev(0:nxm1,0:nym1,1)=0.
     116    wzlev(0:nxm1,0:nym1,1)=0.
     117    rhohn(0:nxm1,0:nym1,1)=pold(0:nxm1,0:nym1)/(r_air*tvold(0:nxm1,0:nym1))
     118
     119! Compute heights of eta levels
     120!******************************
     121
     122    do kz=2,nuvz
     123      pint(0:nxm1,0:nym1)=akz(kz)+bkz(kz)*psn(0:nxm1,0:nym1,1,n,l)
     124      tv(0:nxm1,0:nym1)=tthn(0:nxm1,0:nym1,kz,n,l)*(1.+0.608*qvhn(0:nxm1,0:nym1,kz,n,l))
     125      rhohn(0:nxm1,0:nym1,kz)=pint(0:nxm1,0:nym1)/(r_air*tv(0:nxm1,0:nym1))
     126
     127      where (abs(tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1)).gt.0.2)
     128        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
     129             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))* &
     130             (tv(0:nxm1,0:nym1)-tvold(0:nxm1,0:nym1))/&
     131             &log(tv(0:nxm1,0:nym1)/tvold(0:nxm1,0:nym1))
     132      elsewhere
     133        uvzlev(0:nxm1,0:nym1,kz)=uvzlev(0:nxm1,0:nym1,kz-1)+const*&
     134             &log(pold(0:nxm1,0:nym1)/pint(0:nxm1,0:nym1))*tv(0:nxm1,0:nym1)
     135      endwhere
     136
     137      tvold(0:nxm1,0:nym1)=tv(0:nxm1,0:nym1)
     138      pold(0:nxm1,0:nym1)=pint(0:nxm1,0:nym1)
     139
     140    end do
     141
     142    do kz=2,nwz-1
     143      wzlev(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)+uvzlev(0:nxm1,0:nym1,kz))/2.
     144    end do
     145    wzlev(0:nxm1,0:nym1,nwz)=wzlev(0:nxm1,0:nym1,nwz-1)+ &
     146         uvzlev(0:nxm1,0:nym1,nuvz)-uvzlev(0:nxm1,0:nym1,nuvz-1)
     147
     148
     149    pinmconv(0:nxm1,0:nym1,1)=(uvzlev(0:nxm1,0:nym1,2))/ &
     150         ((aknew(2)+bknew(2)*psn(0:nxm1,0:nym1,1,n,l))- &
     151         (aknew(1)+bknew(1)*psn(0:nxm1,0:nym1,1,n,l)))
     152    do kz=2,nz-1
     153      pinmconv(0:nxm1,0:nym1,kz)=(uvzlev(0:nxm1,0:nym1,kz+1)-uvzlev(0:nxm1,0:nym1,kz-1))/ &
     154           ((aknew(kz+1)+bknew(kz+1)*psn(0:nxm1,0:nym1,1,n,l))- &
     155           (aknew(kz-1)+bknew(kz-1)*psn(0:nxm1,0:nym1,1,n,l)))
     156    end do
     157    pinmconv(0:nxm1,0:nym1,nz)=(uvzlev(0:nxm1,0:nym1,nz)-uvzlev(0:nxm1,0:nym1,nz-1))/ &
     158         ((aknew(nz)+bknew(nz)*psn(0:nxm1,0:nym1,1,n,l))- &
     159         (aknew(nz-1)+bknew(nz-1)*psn(0:nxm1,0:nym1,1,n,l)))
     160
     161! Levels, where u,v,t and q are given
     162!************************************
     163
     164    uun(0:nxm1,0:nym1,1,n,l)=uuhn(0:nxm1,0:nym1,1,l)
     165    vvn(0:nxm1,0:nym1,1,n,l)=vvhn(0:nxm1,0:nym1,1,l)
     166    ttn(0:nxm1,0:nym1,1,n,l)=tthn(0:nxm1,0:nym1,1,n,l)
     167    qvn(0:nxm1,0:nym1,1,n,l)=qvhn(0:nxm1,0:nym1,1,n,l)
     168    if (readclouds_nest(l)) then
     169      clwcn(0:nxm1,0:nym1,1,n,l)=clwchn(0:nxm1,0:nym1,1,n,l)
     170      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,1,n,l)=ciwchn(0:nxm1,0:nym1,1,n,l)
     171    end if
     172    pvn(0:nxm1,0:nym1,1,n,l)=pvhn(0:nxm1,0:nym1,1,l)
     173    rhon(0:nxm1,0:nym1,1,n,l)=rhohn(0:nxm1,0:nym1,1)
     174
     175    uun(0:nxm1,0:nym1,nz,n,l)=uuhn(0:nxm1,0:nym1,nuvz,l)
     176    vvn(0:nxm1,0:nym1,nz,n,l)=vvhn(0:nxm1,0:nym1,nuvz,l)
     177    ttn(0:nxm1,0:nym1,nz,n,l)=tthn(0:nxm1,0:nym1,nuvz,n,l)
     178    qvn(0:nxm1,0:nym1,nz,n,l)=qvhn(0:nxm1,0:nym1,nuvz,n,l)
     179    if (readclouds_nest(l)) then
     180      clwcn(0:nxm1,0:nym1,nz,n,l)=clwchn(0:nxm1,0:nym1,nuvz,n,l)
     181      if (.not.sumclouds_nest(l)) ciwcn(0:nxm1,0:nym1,nz,n,l)=ciwchn(0:nxm1,0:nym1,nuvz,n,l)
     182    end if
     183    pvn(0:nxm1,0:nym1,nz,n,l)=pvhn(0:nxm1,0:nym1,nuvz,l)
     184    rhon(0:nxm1,0:nym1,nz,n,l)=rhohn(0:nxm1,0:nym1,nuvz)
     185
     186
     187    kmin=2
     188    idx=kmin
     189    do iz=2,nz-1
     190      do jy=0,nyn(l)-1
     191        do ix=0,nxn(l)-1
     192          if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
    164193            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    165194            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    166195            ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
    167196            qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
     197!hg adding the cloud water
     198            if (readclouds_nest(l)) then
     199              clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
     200              if (.not.sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
     201            end if
     202!hg
    168203            pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
    169204            rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
    170             goto 30
     205          else
     206            innuvz: do kz=idx(ix,jy),nuvz
     207              if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
     208                   (height(iz).le.uvzlev(ix,jy,kz))) then
     209                idx(ix,jy)=kz
     210                exit innuvz
     211              endif
     212            enddo innuvz
    171213          endif
    172           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    173           (height(iz).le.uvwzlev(ix,jy,kz))) then
    174            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    175            dz2=uvwzlev(ix,jy,kz)-height(iz)
    176            dz=dz1+dz2
    177            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    178            uuhn(ix,jy,kz,l)*dz1)/dz
    179            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    180            vvhn(ix,jy,kz,l)*dz1)/dz
    181            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    182            tthn(ix,jy,kz,n,l)*dz1)/dz
    183            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    184            qvhn(ix,jy,kz,n,l)*dz1)/dz
    185            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    186            pvhn(ix,jy,kz,l)*dz1)/dz
    187            rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    188            kmin=kz
    189            goto 30
     214        enddo
     215      enddo
     216      do jy=0,nyn(l)-1
     217        do ix=0,nxn(l)-1
     218          if(height(iz).le.uvzlev(ix,jy,nuvz)) then
     219            kz=idx(ix,jy)
     220            dz1=height(iz)-uvzlev(ix,jy,kz-1)
     221            dz2=uvzlev(ix,jy,kz)-height(iz)
     222            dz=dz1+dz2
     223            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+uuhn(ix,jy,kz,l)*dz1)/dz
     224            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+vvhn(ix,jy,kz,l)*dz1)/dz
     225            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
     226                 +tthn(ix,jy,kz,n,l)*dz1)/dz
     227            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
     228                 +qvhn(ix,jy,kz,n,l)*dz1)/dz
     229!hg adding the cloud water
     230            if (readclouds_nest(l)) then
     231              clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2+clwchn(ix,jy,kz,n,l)*dz1)/dz
     232              if (.not.sumclouds_nest(l)) &
     233                   &ciwcn(ix,jy,iz,n,l)=(ciwchn(ix,jy,kz-1,n,l)*dz2+ciwchn(ix,jy,kz,n,l)*dz1)/dz
     234            end if
     235!hg
     236            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+pvhn(ix,jy,kz,l)*dz1)/dz
     237            rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)*dz2+rhohn(ix,jy,kz)*dz1)/dz
    190238          endif
     239        enddo
     240      enddo
     241    enddo
     242
     243
     244
     245!         do kz=kmin,nuvz
     246!           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     247!             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
     248!             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
     249!             ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
     250!             qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
     251!             pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
     252!             rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
     253!             goto 30
     254!           endif
     255!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     256!           (height(iz).le.uvwzlev(ix,jy,kz))) then
     257!            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     258!            dz2=uvwzlev(ix,jy,kz)-height(iz)
     259!            dz=dz1+dz2
     260!            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
     261!            uuhn(ix,jy,kz,l)*dz1)/dz
     262!            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
     263!            vvhn(ix,jy,kz,l)*dz1)/dz
     264!            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
     265!            tthn(ix,jy,kz,n,l)*dz1)/dz
     266!            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
     267!            qvhn(ix,jy,kz,n,l)*dz1)/dz
     268!            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
     269!            pvhn(ix,jy,kz,l)*dz1)/dz
     270!            rhon(ix,jy,iz,n,l)=(rhohn(kz-1)*dz2+rhohn(kz)*dz1)/dz
     271!            kmin=kz
     272!            goto 30
     273!           endif
     274!         end do
     275! 30      continue
     276!       end do
     277
     278
     279! Levels, where w is given
     280!*************************
     281
     282    wwn(0:nxm1,0:nym1,1,n,l)=wwhn(0:nxm1,0:nym1,1,l)*pinmconv(0:nxm1,0:nym1,1)
     283    wwn(0:nxm1,0:nym1,nz,n,l)=wwhn(0:nxm1,0:nym1,nwz,l)*pinmconv(0:nxm1,0:nym1,nz)
     284    kmin=2
     285    idx=kmin
     286    do iz=2,nz
     287      do jy=0,nym1
     288        do ix=0,nxm1
     289          inn:         do kz=idx(ix,jy),nwz
     290            if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
     291                 height(iz).le.wzlev(ix,jy,kz)) then
     292              idx(ix,jy)=kz
     293              exit inn
     294            endif
     295          enddo inn
     296        enddo
     297      enddo
     298      do jy=0,nym1
     299        do ix=0,nxm1
     300          kz=idx(ix,jy)
     301          dz1=height(iz)-wzlev(ix,jy,kz-1)
     302          dz2=wzlev(ix,jy,kz)-height(iz)
     303          dz=dz1+dz2
     304          wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(ix,jy,kz-1)*dz2 &
     305               +wwhn(ix,jy,kz,l)*pinmconv(ix,jy,kz)*dz1)/dz
     306        enddo
     307      enddo
     308    enddo
     309
     310!       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
     311!       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
     312!       kmin=2
     313!       do iz=2,nz
     314!         do kz=kmin,nwz
     315!           if ((height(iz).gt.wzlev(kz-1)).and. &
     316!           (height(iz).le.wzlev(kz))) then
     317!             dz1=height(iz)-wzlev(kz-1)
     318!             dz2=wzlev(kz)-height(iz)
     319!             dz=dz1+dz2
     320!             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
     321!             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
     322!             kmin=kz
     323!             goto 40
     324!           endif
     325!         end do
     326! 40      continue
     327!       end do
     328
     329! Compute density gradients at intermediate levels
     330!*************************************************
     331
     332    drhodzn(0:nxm1,0:nym1,1,n,l)=(rhon(0:nxm1,0:nym1,2,n,l)-rhon(0:nxm1,0:nym1,1,n,l))/ &
     333         (height(2)-height(1))
     334    do kz=2,nz-1
     335      drhodzn(0:nxm1,0:nym1,kz,n,l)=(rhon(0:nxm1,0:nym1,kz+1,n,l)-rhon(0:nxm1,0:nym1,kz-1,n,l))/ &
     336           (height(kz+1)-height(kz-1))
     337    end do
     338    drhodzn(0:nxm1,0:nym1,nz,n,l)=drhodzn(0:nxm1,0:nym1,nz-1,n,l)
     339
     340
     341! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
     342! (height(2)-height(1))
     343! do kz=2,nz-1
     344!   drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
     345!   rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     346! end do
     347! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
     348
     349
     350
     351!   end do
     352! end do
     353
     354
     355!****************************************************************
     356! Compute slope of eta levels in windward direction and resulting
     357! vertical wind correction
     358!****************************************************************
     359
     360    do jy=1,nyn(l)-2
     361!    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     362      cosf(jy)=1./cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     363    end do
     364
     365    kmin=2
     366    idx=kmin
     367    do iz=2,nz-1
     368      do jy=1,nyn(l)-2
     369        do ix=1,nxn(l)-2
     370
     371          inneta: do kz=idx(ix,jy),nz
     372            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
     373                 (height(iz).le.uvzlev(ix,jy,kz))) then
     374              idx(ix,jy)=kz
     375              exit inneta
     376            endif
     377          enddo inneta
     378        enddo
     379      enddo
     380
     381      do jy=1,nyn(l)-2
     382        do ix=1,nxn(l)-2
     383          kz=idx(ix,jy)
     384          dz1=height(iz)-uvzlev(ix,jy,kz-1)
     385          dz2=uvzlev(ix,jy,kz)-height(iz)
     386          dz=dz1+dz2
     387          ix1=ix-1
     388          jy1=jy-1
     389          ixp=ix+1
     390          jyp=jy+1
     391
     392          dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
     393          dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
     394          dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
     395
     396          dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
     397          dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
     398          dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
     399
     400          wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*uun(ix,jy,iz,n,l)*dxconst*xresoln(l)*cosf(jy)+&
     401               &dzdy*vvn(ix,jy,iz,n,l)*dyconst*yresoln(l))
     402
    191403        end do
    192 30      continue
     404
    193405      end do
    194 
    195 
    196   ! Levels, where w is given
    197   !*************************
    198 
    199       wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
    200       wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
    201       kmin=2
    202       do iz=2,nz
    203         do kz=kmin,nwz
    204           if ((height(iz).gt.wzlev(kz-1)).and. &
    205           (height(iz).le.wzlev(kz))) then
    206             dz1=height(iz)-wzlev(kz-1)
    207             dz2=wzlev(kz)-height(iz)
    208             dz=dz1+dz2
    209             wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
    210             +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
    211             kmin=kz
    212             goto 40
    213           endif
     406    end do
     407
     408
     409!   do jy=1,nyn(l)-2
     410!     do ix=1,nxn(l)-2
     411!       kmin=2
     412!       do iz=2,nz-1
     413
     414!         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf(jy)
     415!         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
     416
     417!         do kz=kmin,nz
     418!           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     419!                (height(iz).le.uvwzlev(ix,jy,kz))) then
     420!             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     421!             dz2=uvwzlev(ix,jy,kz)-height(iz)
     422!             dz=dz1+dz2
     423!             kl=kz-1
     424!             klp=kz
     425!             kmin=kz
     426!             goto 47
     427!           endif
     428!         end do
     429
     430! 47      ix1=ix-1
     431!         jy1=jy-1
     432!         ixp=ix+1
     433!         jyp=jy+1
     434
     435!         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
     436!         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
     437!         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
     438
     439!         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
     440!         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
     441!         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
     442
     443!         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
     444
     445!       end do
     446
     447!     end do
     448!   end do
     449
     450
     451!write (*,*) 'initializing nested cloudsn, n:',n
     452!   create a cloud and rainout/washout field, cloudsn occur where rh>80%
     453
     454
     455!*********************************************************************************** 
     456    if (readclouds_nest(l)) then !HG METHOD
     457! The method is loops all grids vertically and constructs the 3D matrix for clouds
     458! Cloud top and cloud bottom gid cells are assigned as well as the total column
     459! cloud water. For precipitating grids, the type and whether it is in or below
     460! cloud scavenging are assigned with numbers 2-5 (following the old metod).
     461! Distinction is done for lsp and convp though they are treated the same in regards
     462! to scavenging. Also clouds that are not precipitating are defined which may be
     463! to include future cloud processing by non-precipitating-clouds.
     464!***********************************************************************************
     465      write(*,*) 'Nested ECMWF fields: using cloud water'
     466      clwn(0:nxm1,0:nym1,:,n,l)=0.0
     467!    icloud_stats(0:nxm1,0:nym1,:,n)=0.0
     468      clw4n(0:nxm1,0:nym1,n,l)=0.0
     469      cloudsn(0:nxm1,0:nym1,:,n,l)=0
     470! If water/ice are read separately into clwc and ciwc, store sum in clwcn
     471      if (.not.sumclouds_nest(l)) then
     472        clwcn(0:nxm1,0:nym1,:,n,l) = clwcn(0:nxm1,0:nym1,:,n,l) + ciwcn(0:nxm1,0:nym1,:,n,l)
     473      end if
     474      do jy=0,nym1
     475        do ix=0,nxm1
     476          lsp=lsprecn(ix,jy,1,n,l)
     477          convp=convprecn(ix,jy,1,n,l)
     478          prec=lsp+convp
     479          tot_cloud_h=0
     480! Find clouds in the vertical
     481          do kz=1, nz-1 !go from top to bottom
     482            if (clwcn(ix,jy,kz,n,l).gt.0) then     
     483! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
     484              clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)*rhon(ix,jy,kz,n,l))*(height(kz+1)-height(kz))
     485              tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
     486              clw4n(ix,jy,n,l) = clw4n(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
     487!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
     488!           icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
     489              cloudh_min=min(height(kz+1),height(kz))
     490!ZHG 2015 extra for testing
     491!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
     492!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
     493!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
     494!ZHG
     495            endif
     496          end do
     497
     498! If Precipitation. Define removal type in the vertical
     499          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     500
     501            do kz=nz,1,-1 !go Bottom up!
     502              if (clwn(ix,jy,kz,n,l).gt.0.0) then ! is in cloud
     503                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
     504                cloudsn(ix,jy,kz,n,l)=1                               ! is a cloud
     505                if (lsp.ge.convp) then
     506                  cloudsn(ix,jy,kz,n,l)=3                            ! lsp in-cloud
     507                else
     508                  cloudsn(ix,jy,kz,n,l)=2                             ! convp in-cloud
     509                endif                                              ! convective or large scale
     510              else if((clwn(ix,jy,kz,n,l).le.0.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
     511                if (lsp.ge.convp) then
     512                  cloudsn(ix,jy,kz,n,l)=5                             ! lsp dominated washout
     513                else
     514                  cloudsn(ix,jy,kz,n,l)=4                             ! convp dominated washout
     515                endif                                              ! convective or large scale
     516              endif
     517
     518              if (height(kz).ge. 19000) then                        ! set a max height for removal
     519                cloudsn(ix,jy,kz,n,l)=0
     520              endif !clw>0
     521            end do !nz
     522          endif ! precipitation
    214523        end do
    215 40      continue
    216524      end do
    217 
    218   ! Compute density gradients at intermediate levels
    219   !*************************************************
    220 
    221       drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    222       (height(2)-height(1))
    223       do kz=2,nz-1
    224         drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    225         rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     525!********************************************************************
     526    else ! old method:
     527!********************************************************************
     528      write(*,*) 'Nested fields: using cloud water from Parameterization'
     529      do jy=0,nyn(l)-1
     530        do ix=0,nxn(l)-1
     531          rain_cloud_above(ix,jy)=0
     532          lsp=lsprecn(ix,jy,1,n,l)
     533          convp=convprecn(ix,jy,1,n,l)
     534          cloudshn(ix,jy,n,l)=0
     535          do kz_inv=1,nz-1
     536            kz=nz-kz_inv+1
     537            pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
     538            rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
     539            cloudsn(ix,jy,kz,n,l)=0
     540            if (rh.gt.0.8) then ! in cloud
     541              if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
     542                rain_cloud_above(ix,jy)=1
     543                cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
     544                     height(kz)-height(kz-1)
     545                if (lsp.ge.convp) then
     546                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
     547                else
     548                  cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
     549                endif
     550              else ! no precipitation
     551                cloudsn(ix,jy,kz,n,l)=1 ! cloud
     552              endif
     553            else ! no cloud
     554              if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
     555                if (lsp.ge.convp) then
     556                  cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
     557                else
     558                  cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
     559                endif
     560              endif
     561            endif
     562          end do
     563        end do
    226564      end do
    227       drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
    228 
    229     end do
    230   end do
    231 
    232 
    233   !****************************************************************
    234   ! Compute slope of eta levels in windward direction and resulting
    235   ! vertical wind correction
    236   !****************************************************************
    237 
    238   do jy=1,nyn(l)-2
    239     cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    240     do ix=1,nxn(l)-2
    241 
    242       kmin=2
    243       do iz=2,nz-1
    244 
    245         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
    246         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    247 
    248         do kz=kmin,nz
    249           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    250                (height(iz).le.uvwzlev(ix,jy,kz))) then
    251             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    252             dz2=uvwzlev(ix,jy,kz)-height(iz)
    253             dz=dz1+dz2
    254             kl=kz-1
    255             klp=kz
    256             kmin=kz
    257             goto 47
    258           endif
    259         end do
    260 
    261 47      ix1=ix-1
    262         jy1=jy-1
    263         ixp=ix+1
    264         jyp=jy+1
    265 
    266         dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
    267         dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
    268         dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
    269 
    270         dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
    271         dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
    272         dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
    273 
    274         wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdx*ui+dzdy*vi)
    275 
    276       end do
    277 
    278     end do
    279   end do
    280 
    281 
    282   !write (*,*) 'initializing nested cloudsn, n:',n
    283   !   create a cloud and rainout/washout field, cloudsn occur where rh>80%
    284   write(*,*) 'Nested fields: using cloud water from Parameterization'
    285   do jy=0,nyn(l)-1
    286     do ix=0,nxn(l)-1
    287       rain_cloud_above=0
    288       lsp=lsprecn(ix,jy,1,n,l)
    289       convp=convprecn(ix,jy,1,n,l)
    290       cloudsnh(ix,jy,n,l)=0
    291       do kz_inv=1,nz-1
    292          kz=nz-kz_inv+1
    293          pressure=rhon(ix,jy,kz,n,l)*r_air*ttn(ix,jy,kz,n,l)
    294          rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
    295          cloudsn(ix,jy,kz,n,l)=0
    296          if (rh.gt.0.8) then ! in cloud
    297             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then
    298                rain_cloud_above=1
    299                cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    300                height(kz)-height(kz-1)
    301                if (lsp.ge.convp) then
    302                   cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
    303                else
    304                   cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
    305                endif
    306             else ! no precipitation
    307                   cloudsn(ix,jy,kz,n,l)=1 ! cloud
    308             endif
    309          else ! no cloud
    310             if (rain_cloud_above.eq.1) then ! scavenging
    311                if (lsp.ge.convp) then
    312                   cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
    313                else
    314                   cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
    315                endif
    316             endif
    317          endif
    318       end do
    319     end do
    320   end do
    321 
    322   end do
     565    end if
     566
     567  end do ! end loop over nests
    323568
    324569end subroutine verttransform_nests

  • src/wetdepo.f90

    r41d8574 rdb712a8  
    2121
    2222subroutine wetdepo(itime,ltsample,loutnext)
    23   !                  i      i        i
    24   !*****************************************************************************
    25   !                                                                            *
    26   ! Calculation of wet deposition using the concept of scavenging coefficients.*
    27   ! For lack of detailed information, washout and rainout are jointly treated. *
    28   ! It is assumed that precipitation does not occur uniformly within the whole *
    29   ! grid cell, but that only a fraction of the grid cell experiences rainfall. *
    30   ! This fraction is parameterized from total cloud cover and rates of large   *
    31   ! scale and convective precipitation.                                        *
    32   !                                                                            *
    33   !    Author: A. Stohl                                                        *
    34   !                                                                            *
    35   !    1 December 1996                                                         *
    36   !                                                                            *
    37   ! Correction by Petra Seibert, Sept 2002:                                    *
    38   ! use centred precipitation data for integration                             *
    39   ! Code may not be correct for decay of deposition!                           *
    40   !                                                                            *
    41   !*****************************************************************************
    42   !                                                                            *
    43   ! Variables:                                                                 *
    44   ! cc [0-1]           total cloud cover                                       *
    45   ! convp [mm/h]       convective precipitation rate                           *
    46   ! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
    47   ! ix,jy              indices of output grid cell for each particle           *
    48   ! itime [s]          actual simulation time [s]                              *
    49   ! jpart              particle index                                          *
    50   ! ldeltat [s]        interval since radioactive decay was computed           *
    51   ! lfr, cfr           area fraction covered by precipitation for large scale  *
    52   !                    and convective precipitation (dependent on prec. rate)  *
    53   ! loutnext [s]       time for which gridded deposition is next output        *
    54   ! loutstep [s]       interval at which gridded deposition is output          *
    55   ! lsp [mm/h]         large scale precipitation rate                          *
    56   ! ltsample [s]       interval over which mass is deposited                   *
    57   ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
    58   ! wetdeposit         mass that is wet deposited                              *
    59   ! wetgrid            accumulated deposited mass on output grid               *
    60   ! wetscav            scavenging coefficient                                  *
    61   !                                                                            *
    62   ! Constants:                                                                 *
    63   !                                                                            *
    64   !*****************************************************************************
     23!                  i      i        i
     24!*****************************************************************************
     25!                                                                            *
     26! Calculation of wet deposition using the concept of scavenging coefficients.*
     27! For lack of detailed information, washout and rainout are jointly treated. *
     28! It is assumed that precipitation does not occur uniformly within the whole *
     29! grid cell, but that only a fraction of the grid cell experiences rainfall. *
     30! This fraction is parameterized from total cloud cover and rates of large   *
     31! scale and convective precipitation.                                        *
     32!                                                                            *
     33!    Author: A. Stohl                                                        *
     34!                                                                            *
     35!    1 December 1996                                                         *
     36!                                                                            *
     37! Correction by Petra Seibert, Sept 2002:                                    *
     38! use centred precipitation data for integration                             *
     39! Code may not be correct for decay of deposition!                           *
     40!                                                                            *
     41!*****************************************************************************
     42!                                                                            *
     43! Variables:                                                                 *
     44! cc [0-1]           total cloud cover                                       *
     45! convp [mm/h]       convective precipitation rate                           *
     46! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
     47! ix,jy              indices of output grid cell for each particle           *
     48! itime [s]          actual simulation time [s]                              *
     49! jpart              particle index                                          *
     50! ldeltat [s]        interval since radioactive decay was computed           *
     51! lfr, cfr           area fraction covered by precipitation for large scale  *
     52!                    and convective precipitation (dependent on prec. rate)  *
     53! loutnext [s]       time for which gridded deposition is next output        *
     54! loutstep [s]       interval at which gridded deposition is output          *
     55! lsp [mm/h]         large scale precipitation rate                          *
     56! ltsample [s]       interval over which mass is deposited                   *
     57! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
     58! wetdeposit         mass that is wet deposited                              *
     59! wetgrid            accumulated deposited mass on output grid               *
     60! wetscav            scavenging coefficient                                  *
     61!                                                                            *
     62! Constants:                                                                 *
     63!                                                                            *
     64!*****************************************************************************
    6565
    6666  use point_mod
    … …  
    7171
    7272  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
     73  integer :: ngrid,itage,nage,hz,il,interp_time, n
     74  integer(kind=1) :: clouds_v
    7475  integer :: ks, kp
    7576!  integer :: n1,n2, icbot,ictop, indcloud !TEST
    … …  
    7980  real :: wetdeposit(maxspec),restmass
    8081  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    81   !save lfr,cfr
     82!save lfr,cfr
    8283
    8384  real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
    … …  
    9192  integer :: blc_count, inc_count
    9293  real    :: Si_dummy, wetscav_dummy
    93 
    94 
    95 
    96   ! Compute interval since radioactive decay of deposited mass was computed
    97   !************************************************************************
     94  logical :: readclouds_this_nest
     95
     96
     97! Compute interval since radioactive decay of deposited mass was computed
     98!************************************************************************
    9899
    99100  if (itime.le.loutnext) then
    … …  
    103104  endif
    104105
    105   ! Loop over all particles
    106   !************************
     106! Loop over all particles
     107!************************
    107108
    108109  blc_count=0
    … …  
    118119    endif
    119120
    120   ! Determine age class of the particle
     121! Determine age class of the particle
    121122    itage=abs(itra1(jpart)-itramem(jpart))
    122123    do nage=1,nageclass
    … …  
    126127
    127128
    128   ! Determine which nesting level to be used
    129   !*****************************************
     129! Determine which nesting level to be used
     130!*****************************************
    130131
    131132    ngrid=0
    132133    do j=numbnests,1,-1
    133134      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    134       (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
     135           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    135136        ngrid=j
    136137        goto 23
    … …  
    140141
    141142
    142   ! Determine nested grid coordinates
    143   !**********************************
     143! Determine nested grid coordinates
     144!**********************************
    144145
    145146    if (ngrid.gt.0) then
    … …  
    148149      ix=int(xtn)
    149150      jy=int(ytn)
     151      if (readclouds_nest(ngrid)) then
     152        readclouds_this_nest=.true.
     153      else
     154        readclouds_this_nest=.false.
     155      end if
    150156    else
    151157      ix=int(xtra1(jpart))
    … …  
    154160
    155161
    156   ! Interpolate large scale precipitation, convective precipitation and
    157   ! total cloud cover
    158   ! Note that interpolated time refers to itime-0.5*ltsample [PS]
    159   !********************************************************************
     162! Interpolate large scale precipitation, convective precipitation and
     163! total cloud cover
     164! Note that interpolated time refers to itime-0.5*ltsample [PS]
     165!********************************************************************
    160166    interp_time=nint(itime-0.5*ltsample)
    161    
     167
    162168    if (ngrid.eq.0) then
    163169      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    164       1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    165       memtime(1),memtime(2),interp_time,lsp,convp,cc)
     170           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
     171           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    166172    else
    167173      call interpol_rain_nests(lsprecn,convprecn,tccn, &
    168       nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    169       memtime(1),memtime(2),interp_time,lsp,convp,cc)
     174           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
     175           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    170176    endif
    171177
    … …  
    173179    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    174180
    175   ! get the level were the actual particle is in
     181! get the level were the actual particle is in
    176182    do il=2,nz
    177183      if (height(il).gt.ztra1(jpart)) then
    178184        hz=il-1
    179         goto 26
     185!        goto 26
     186        exit
    180187      endif
    181188    end do
    182 26  continue
     189!26  continue
    183190
    184191    n=memind(2)
    … …  
    190197      clouds_h=cloudsh(ix,jy,n)
    191198    else
    192       ! new removal not implemented for nests yet
    193199      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    194       clouds_h=cloudsnh(ix,jy,n,ngrid)
    195     endif
    196 
    197   ! if there is no precipitation or the particle is above the clouds no
    198   ! scavenging is done
     200      clouds_h=cloudshn(ix,jy,n,ngrid)
     201    endif
     202
     203! if there is no precipitation or the particle is above the clouds no
     204! scavenging is done
    199205
    200206    if (clouds_v.le.1) goto 20
    201  
    202   ! 1) Parameterization of the the area fraction of the grid cell where the
    203   !    precipitation occurs: the absolute limit is the total cloud cover, but
    204   !    for low precipitation rates, an even smaller fraction of the grid cell
    205   !    is used. Large scale precipitation occurs over larger areas than
    206   !    convective precipitation.
    207   !**************************************************************************
     207
     208! 1) Parameterization of the the area fraction of the grid cell where the
     209!    precipitation occurs: the absolute limit is the total cloud cover, but
     210!    for low precipitation rates, an even smaller fraction of the grid cell
     211!    is used. Large scale precipitation occurs over larger areas than
     212!    convective precipitation.
     213!**************************************************************************
    208214
    209215    if (lsp.gt.20.) then
    … …  
    232238
    233239
    234   !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
    235   ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
    236   ! for now they are treated the same
     240!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
     241! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
     242! for now they are treated the same
    237243    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
    238244    grfraction(2)=max(0.05,cc*(lfr(i)))
    … …  
    240246
    241247
    242   ! 2) Computation of precipitation rate in sub-grid cell
    243   !******************************************************
     248! 2) Computation of precipitation rate in sub-grid cell
     249!******************************************************
    244250    prec(1)=(lsp+convp)/grfraction(1)
    245251    prec(2)=(lsp)/grfraction(2)
    … …  
    247253
    248254
    249   ! 3) Computation of scavenging coefficients for all species
    250   !    Computation of wet deposition
    251   !**********************************************************
     255! 3) Computation of scavenging coefficients for all species
     256!    Computation of wet deposition
     257!**********************************************************
    252258
    253259    do ks=1,nspec      ! loop over species
    … …  
    255261      wetscav=0.   
    256262
    257       !ZHG test if it nested?
     263!ZHG test if it nested?
    258264      if (ngrid.gt.0) then
    259265        act_temp=ttn(ix,jy,hz,n,ngrid)
    … …  
    261267        act_temp=tt(ix,jy,hz,n)
    262268      endif
    263      
    264 
    265   !***********************
    266   ! BELOW CLOUD SCAVENGING
    267   !*********************** 
     269
     270
     271!***********************
     272! BELOW CLOUD SCAVENGING
     273!*********************** 
    268274      if (clouds_v.ge.4) then !below cloud
    269275
    … …  
    274280          if (dquer(ks) .le. 0.) then  !is gas
    275281            wetscav=weta(ks)*prec(1)**wetb(ks)
    276            
     282
    277283!AEROSOL
    278284          else !is particle
    … …  
    281287! for particles larger than 10 um use the largest size defined in the parameterizations (10um)
    282288            dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
    283             if (act_temp .ge. 273 .and. weta(ks).gt.0.)  then !Rain
    284               ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
    285               ! the below-cloud scavenging (rain efficienty)
    286               ! parameter A (=weta) from SPECIES file
     289            if (act_temp .ge. 273. .and. weta(ks).gt.0.)  then !Rain
     290! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
     291! the below-cloud scavenging (rain efficienty)
     292! parameter A (=weta) from SPECIES file
    287293              wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
    288294                   (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
    289295
    290             elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.)  then ! Snow
    291               ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
    292               ! the below-cloud scavenging (Snow efficiency)
    293               ! parameter B (=wetb) from SPECIES file
     296            elseif (act_temp .lt. 273. .and. wetb(ks).gt.0.)  then ! Snow
     297! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
     298! the below-cloud scavenging (Snow efficiency)
     299! parameter B (=wetb) from SPECIES file
    294300              wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
    295301                   (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
    … …  
    304310
    305311
    306   !********************
    307   ! IN CLOUD SCAVENGING
    308       !******************************************************
     312!********************
     313! IN CLOUD SCAVENGING
     314!********************
    309315      if (clouds_v.lt.4) then ! In-cloud
    310316! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
    … …  
    315321          if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
    316322
    317           !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
    318           if (readclouds) then                  !icloud_stats(ix,jy,4,n) has units kg/m2
    319 !            cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
    320 ! ESO: stop using icloud_stats
    321             cl =clw4(ix,jy,n)*(grfraction(1)/cc)
     323!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
     324! nested fields
     325          if (ngrid.gt.0.and.readclouds_this_nest) then
     326            cl=clw4n(ix,jy,n,ngrid)*(grfraction(1)/cc)
     327          else if (ngrid.eq.0.and.readclouds) then
     328            cl=clw4(ix,jy,n)*(grfraction(1)/cc)
    322329          else                                  !parameterize cloudwater m2/m3
    323             !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
     330!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
    324331            cl=1.6E-6*prec(1)**0.36
    325332          endif
    326333
    327             !ZHG: Calculate the partition between liquid and water phase water.
    328             if (act_temp .le. 253) then
    329               liq_frac=0
    330             else if (act_temp .ge. 273) then
    331               liq_frac=1
    332             else
    333               liq_frac =((act_temp-273)/(273-253))**2
    334             endif
    335            ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
    336             frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
    337  
    338   !ZHG Use the activated fraction and the liqid water to calculate the washout
    339 
    340   ! AEROSOL
    341           !**************************************************
     334!ZHG: Calculate the partition between liquid and water phase water.
     335          if (act_temp .le. 253.) then
     336            liq_frac=0
     337          else if (act_temp .ge. 273.) then
     338            liq_frac=1
     339          else
     340            liq_frac =((act_temp-273.)/(273.-253.))**2.
     341          end if
     342! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
     343          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
     344
     345!ZHG Use the activated fraction and the liqid water to calculate the washout
     346
     347! AEROSOL
     348!**************************************************
    342349          if (dquer(ks).gt. 0.) then ! is particle
    343350
    344351            S_i= frac_act/cl
    345352
    346           !*********************
    347   ! GAS
     353!*********************
     354! GAS
    348355          else ! is gas
    349                
     356
    350357            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    351             !REPLACE to switch old/ new scheme
    352             ! S_i=frac_act/cle
     358!REPLACE to switch old/ new scheme
     359! S_i=frac_act/cle
    353360            S_i=1/cle
    354361          endif ! gas or particle
    355362
    356   ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
    357            !OLD
    358           if (readclouds) then
     363! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
     364!OLD
     365          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
    359366            wetscav=S_i*(prec(1)/3.6E6)
    360367          else
    … …  
    362369          endif
    363370
    364 !ZHG 2015 TEST
    365 !          Si_dummy=frac_act/2E-7*prec(1)**0.36
    366 !           wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h
    367 !           if (clouds_v.lt.4) then
    368 !           talltest=talltest+1
    369 !if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')
    370 !if(talltest .lt. 100001)  write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
    371 !if(talltest .lt. 100001)  write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl
    372 !if(talltest .eq. 100001) CLOSE(199)
    373 !if(talltest .eq. 100001) STOP
    374 !
    375 !write(*,*)  'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36',  2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h
    376 !write(*,*)  'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp
    377 !write(*,*)  wetscav, wetscav_dummy
    378 !write(*,*) cc, grfraction(1), cc/grfraction(1)
    379 
    380 !write(*,*)  'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)
    381 
    382 
    383 !write(*,*) '**************************************************'
    384 !write(*,*)  'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)
    385 !write(*,*)  'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy
    386 !write(*,*)  'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1)
    387 
    388 
    389 !           write(*,*) 'PRECIPITATION ,cl  ECMWF , cl PARAMETIZED, clouds_v, lat' &
    390 !                      ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
    391 
    392 !endif
    393371
    394372        endif ! positive in-cloud scavenging parameters given in Species file
    395373      endif !incloud
    396374!END ZHG TEST
    397      
    398   !**************************************************
    399   ! CALCULATE DEPOSITION
    400   !**************************************************
    401   !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    402   !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     375
     376!**************************************************
     377! CALCULATE DEPOSITION
     378!**************************************************
     379!     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
     380!     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
    403381
    404382      if (wetscav.gt.0.) then
    … …  
    421399      if (restmass .gt. smallnum) then
    422400        xmass1(jpart,ks)=restmass
    423   !   depostatistic
    424   !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    425   !   depostatistic
     401!   depostatistic
     402!   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
     403!   depostatistic
    426404      else
    427405        xmass1(jpart,ks)=0.
    428406      endif
    429   !   Correct deposited mass to the last time step when radioactive decay of
    430   !   gridded deposited mass was calculated
     407!   Correct deposited mass to the last time step when radioactive decay of
     408!   gridded deposited mass was calculated
    431409      if (decay(ks).gt.0.) then
    432410        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
    … …  
    436414    end do !all species
    437415
    438   ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
    439   ! Add the wet deposition to accumulated amount on output grid and nested output grid
    440   !*****************************************************************************
     416! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
     417! Add the wet deposition to accumulated amount on output grid and nested output grid
     418!*****************************************************************************
    441419
    442420    if (ldirect.eq.1) then
    … …  
    450428  end do ! all particles
    451429
    452   ! count the total number of below-cloud and in-cloud occurences:
     430! count the total number of below-cloud and in-cloud occurences:
    453431  tot_blc_count=tot_blc_count+blc_count
    454432  tot_inc_count=tot_inc_count+inc_count

  • src/wetdepokernel_nest.f90

    re200b7a rdb712a8  
    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,nage,kp)
     23  !                           i    i       i i i    i
    2524  !*****************************************************************************
    2625  !                                                                            *
    … …  
    5554  integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage
    5655
    57 
     56  real :: dbg_dx, dbg_dy, dbg_xoutshiftn, dbg_youtshiftn, dbg_dxoutn, dbg_dyoutn,dbg_t
    5857
    5958  xl=(x*dx+xoutshiftn)/dxoutn
    6059  yl=(y*dy+youtshiftn)/dyoutn
    61   ix=int(xl)
    62   jy=int(yl)
     60
     61  ! old:
     62  ! ix=int(xl)
     63  ! jy=int(yl)
     64
     65  ! ESO: for xl,yl in range <-.5,-1> we get ix,jy=0 and thus negative
     66  ! wx,wy as the int function rounds upwards for negative numbers.
     67  ! Either use the floor function, or (perhaps more correctly?) use "(xl.gt.-0.5)"
     68  ! in place of "(ix.ge.0)" and similar for the upper boundary.
     69
     70  ! new:
     71  ix=floor(xl)
     72  jy=floor(yl)
     73
    6374  ddx=xl-real(ix)                   ! distance to left cell border
    6475  ddy=yl-real(jy)                   ! distance to lower cell border
     76
    6577
    6678  if (ddx.gt.0.5) then
    … …  
    8092  endif
    8193
    82 
    83   ! Determine mass fractions for four grid points
    84   !**********************************************
     94! Determine mass fractions for four grid points
     95!**********************************************
    8596
    8697  do ks=1,nspec
    87 
    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
     98    if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
     99         (jy.le.numygridn-1)) then
     100      w=wx*wy
    91101      wetgriduncn(ix,jy,ks,kp,nunc,nage)= &
    92102           wetgriduncn(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
    93   endif
     103    endif
    94104
    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)
     105    if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgridn-1).and. &
     106         (jyp.le.numygridn-1)) then
     107      w=(1.-wx)*(1.-wy)
    98108      wetgriduncn(ixp,jyp,ks,kp,nunc,nage)= &
    99109           wetgriduncn(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    100   endif
     110    endif
    101111
    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
     112    if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgridn-1).and. &
     113         (jy.le.numygridn-1)) then
     114      w=(1.-wx)*wy
    105115      wetgriduncn(ixp,jy,ks,kp,nunc,nage)= &
    106116           wetgriduncn(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
    107   endif
     117    endif
    108118
    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)
     119    if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgridn-1).and. &
     120         (jyp.le.numygridn-1)) then
     121      w=wx*(1.-wy)
    112122      wetgriduncn(ix,jyp,ks,kp,nunc,nage)= &
    113123           wetgriduncn(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
    114   endif
     124    endif
    115125
    116126  end do
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG