Changeset f75967d in flexpart.git


Ignore:
Timestamp:
Jan 5, 2016, 12:20:01 PM (3 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
dev, release-10, univie
Children:
d8107c2
Parents:
d6a0977
Message:

Fixed a bug with number of occurances of wet scavenging written to stdout

Location:
src
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • src/FLEXPART_MPI.f90

    rd6a0977 rf75967d  
    436436    stop
    437437  endif
    438    
    439438
    440439
    441440  call timemanager
     441
    442442
    443443! NIK 16.02.2005
    444444  if (lroot) then
    445     call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
     445    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
    446446         & mp_comm_used, mp_ierr)
    447     call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
     447    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
    448448         & mp_comm_used, mp_ierr)
    449449  else
    450450    if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
    451       call MPI_Reduce(tot_blc_count, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
     451      call MPI_Reduce(tot_blc_count, tot_blc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
    452452           & mp_comm_used, mp_ierr)
    453       call MPI_Reduce(tot_inc_count, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
     453      call MPI_Reduce(tot_inc_count, tot_inc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
    454454           & mp_comm_used, mp_ierr)
    455455    end if
  • src/com_mod.f90

    rd6a0977 rf75967d  
    132132!ZHG SEP 2015 wheather or not to read clouds from GRIB
    133133  logical :: readclouds
     134!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
     135  logical :: sumclouds
     136 
    134137
    135138!NIK 16.02.2015
    136   integer :: tot_blc_count, tot_inc_count
     139  integer :: tot_blc_count=0, tot_inc_count=0
    137140
    138141
     
    746749  !*****************************************************************
    747750  integer :: mpi_mode=0 ! .gt. 0 if running MPI version
    748   logical :: lroot=.true. ! true if serial version, or if MPI and root process
     751  logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
    749752
    750753  contains
  • src/timemanager.f90

    rd6a0977 rf75967d  
    138138!ZHG 2015
    139139!CGZ-lifetime: set lifetime to 0
    140   checklifetime(:,:)=0
    141   species_lifetime(:,:)=0
    142   print*, 'Initialized lifetime'
     140  ! checklifetime(:,:)=0
     141  ! species_lifetime(:,:)=0
     142  ! print*, 'Initialized lifetime'
    143143!CGZ-lifetime: set lifetime to 0
    144144 
     
    422422        !CGZ-lifetime: output species lifetime
    423423!ZHG
    424         write(*,*) 'Overview species lifetime in days', &
    425              real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
    426         write(*,*) 'all info:',species_lifetime
     424        ! write(*,*) 'Overview species lifetime in days', &
     425        !      real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
     426        ! write(*,*) 'all info:',species_lifetime
    427427!ZHG
    428428        !CGZ-lifetime: output species lifetime
     
    594594!ZHG 2015
    595595                  !CGZ-lifetime: Check mass fraction left/save lifetime
    596                    if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
     596                   ! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
    597597                       !Mass below 1% of initial >register lifetime
    598                        checklifetime(j,ks)=abs(itra1(j)-itramem(j))
    599                        species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
    600                        species_lifetime(ks,2)= species_lifetime(ks,2)+1
    601                    endif
     598                       ! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
     599                       ! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
     600                       ! species_lifetime(ks,2)= species_lifetime(ks,2)+1
     601                   ! endif
    602602                   !CGZ-lifetime: Check mass fraction left/save lifetime
    603603!ZHG 2015
  • src/timemanager_mpi.f90

    rd6a0977 rf75967d  
    151151
    152152!CGZ-lifetime: set lifetime to 0
    153   checklifetime(:,:)=0
    154   species_lifetime(:,:)=0
    155   print*, 'Initialized lifetime'
     153  ! checklifetime(:,:)=0
     154  ! species_lifetime(:,:)=0
     155  ! print*, 'Initialized lifetime'
    156156!CGZ-lifetime: set lifetime to 0
    157157
     
    564564       
    565565        !CGZ-lifetime: output species lifetime
    566         ! if (lroot) then
     566        if (lroot) then
    567567        !   write(*,*) 'Overview species lifetime in days', &
    568568        !        real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
     
    573573        !     write (*,*) 'timemanager> starting simulation'
    574574        !   end if
    575         ! end if
     575        end if
    576576
    57757745      format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES:    Uncertainty: ',3f7.3)
     
    757757                   
    758758                   !CGZ-lifetime: Check mass fraction left/save lifetime
    759                    if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
     759                   ! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
    760760                       !Mass below 1% of initial >register lifetime
    761                        checklifetime(j,ks)=abs(itra1(j)-itramem(j))
    762 
    763                        species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
    764                        species_lifetime(ks,2)= species_lifetime(ks,2)+1
    765                    endif
     761                   !     checklifetime(j,ks)=abs(itra1(j)-itramem(j))
     762
     763                   !     species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
     764                   !     species_lifetime(ks,2)= species_lifetime(ks,2)+1
     765                   ! endif
    766766                   !CGZ-lifetime: Check mass fraction left/save lifetime
    767767                   
  • src/verttransform.f90

    rd6a0977 rf75967d  
    568568
    569569!*********************************************************************************** 
    570 if (readclouds) then !HG METHOD
     570  if (readclouds) then !HG METHOD
    571571! The method is loops all grids vertically and constructs the 3D matrix for clouds
    572572! Cloud top and cloud bottom gid cells are assigned as well as the total column
     
    577577! to include future cloud processing by non-precipitating-clouds.
    578578!***********************************************************************************
    579   if (readclouds) write(*,*) 'using cloud water from ECMWF'
    580   clw(:,:,:,n)=0
    581   icloud_stats(:,:,:,n)=0
    582   clouds(:,:,:,n)=0
    583   do jy=0,nymin1
    584    do ix=0,nxmin1
    585     lsp=lsprec(ix,jy,1,n)
    586     convp=convprec(ix,jy,1,n)
    587     prec=lsp+convp
    588     tot_cloud_h=0
    589     ! Find clouds in the vertical
    590     do kz=1, nz-1 !go from top to bottom
    591      if (clwc(ix,jy,kz,n).gt.0) then     
    592       ! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
    593       clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
    594       tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
    595       icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
    596       icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
     579    write(*,*) 'using cloud water from ECMWF'
     580    clw(:,:,:,n)=0
     581    icloud_stats(:,:,:,n)=0
     582    clouds(:,:,:,n)=0
     583    do jy=0,nymin1
     584      do ix=0,nxmin1
     585        lsp=lsprec(ix,jy,1,n)
     586        convp=convprec(ix,jy,1,n)
     587        prec=lsp+convp
     588        tot_cloud_h=0
     589! Find clouds in the vertical
     590        do kz=1, nz-1 !go from top to bottom
     591          if (clwc(ix,jy,kz,n).gt.0) then     
     592! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
     593            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
     594            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
     595            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
     596            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
    597597!ZHG 2015 extra for testing
    598       clh(ix,jy,kz,n)=height(kz+1)-height(kz)
    599       icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
    600       icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
     598            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
     599            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
     600            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
    601601!ZHG
    602     endif
    603     end do
    604 
    605    ! If Precipitation. Define removal type in the vertical
    606   if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    607 
    608     do kz=nz,1,-1 !go Bottom up!
    609      if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
    610         cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
    611         clouds(ix,jy,kz,n)=1                               ! is a cloud
    612         if (lsp.ge.convp) then
    613            clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
    614         else
    615           clouds(ix,jy,kz,n)=2                             ! convp in-cloud
    616         endif                                              ! convective or large scale
    617      elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
    618         if (lsp.ge.convp) then
    619           clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
    620         else
    621           clouds(ix,jy,kz,n)=4                             ! convp dominated washout
    622         endif                                              ! convective or large scale
    623      endif
    624 
    625      if (height(kz).ge. 19000) then                        ! set a max height for removal
    626         clouds(ix,jy,kz,n)=0
    627      endif !clw>0
    628     end do !nz
    629    endif ! precipitation
    630   end do
    631  end do
     602          endif
     603        end do
     604
     605! If Precipitation. Define removal type in the vertical
     606        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     607
     608          do kz=nz,1,-1 !go Bottom up!
     609            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
     610              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
     611              clouds(ix,jy,kz,n)=1                               ! is a cloud
     612              if (lsp.ge.convp) then
     613                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
     614              else
     615                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
     616              endif                                              ! convective or large scale
     617            elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
     618              if (lsp.ge.convp) then
     619                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
     620              else
     621                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
     622              endif                                              ! convective or large scale
     623            endif
     624
     625            if (height(kz).ge. 19000) then                        ! set a max height for removal
     626              clouds(ix,jy,kz,n)=0
     627            endif !clw>0
     628          end do !nz
     629        endif ! precipitation
     630      end do
     631    end do
    632632!**************************************************************************
    633  else       ! use old definitions
     633  else       ! use old definitions
    634634!**************************************************************************
    635635!   create a cloud and rainout/washout field, clouds occur where rh>80%
    636636!   total cloudheight is stored at level 0
    637  if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
    638   do jy=0,nymin1
    639     do ix=0,nxmin1
     637! if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
     638    write(*,*) 'using cloud water from Parameterization'
     639    do jy=0,nymin1
     640      do ix=0,nxmin1
    640641! OLD METHOD
    641       rain_cloud_above(ix,jy)=0
    642       lsp=lsprec(ix,jy,1,n)
    643       convp=convprec(ix,jy,1,n)
    644       cloudsh(ix,jy,n)=0
    645       do kz_inv=1,nz-1
    646         kz=nz-kz_inv+1
    647         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    648         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    649         clouds(ix,jy,kz,n)=0
    650         if (rh.gt.0.8) then ! in cloud
    651           if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    652             rain_cloud_above(ix,jy)=1
    653             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    654                  height(kz)-height(kz-1)
    655             if (lsp.ge.convp) then
    656               clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    657             else
    658               clouds(ix,jy,kz,n)=2 ! convp dominated rainout
     642        rain_cloud_above(ix,jy)=0
     643        lsp=lsprec(ix,jy,1,n)
     644        convp=convprec(ix,jy,1,n)
     645        cloudsh(ix,jy,n)=0
     646        do kz_inv=1,nz-1
     647          kz=nz-kz_inv+1
     648          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     649          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     650          clouds(ix,jy,kz,n)=0
     651          if (rh.gt.0.8) then ! in cloud
     652            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     653              rain_cloud_above(ix,jy)=1
     654              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
     655                   height(kz)-height(kz-1)
     656              if (lsp.ge.convp) then
     657                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
     658              else
     659                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
     660              endif
     661            else ! no precipitation
     662              clouds(ix,jy,kz,n)=1 ! cloud
    659663            endif
    660           else ! no precipitation
    661             clouds(ix,jy,kz,n)=1 ! cloud
    662           endif
    663         else ! no cloud
    664           if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
    665             if (lsp.ge.convp) then
    666               clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    667             else
    668               clouds(ix,jy,kz,n)=4 ! convp dominated washout
     664          else ! no cloud
     665            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
     666              if (lsp.ge.convp) then
     667                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     668              else
     669                clouds(ix,jy,kz,n)=4 ! convp dominated washout
     670              endif
    669671            endif
    670672          endif
    671         endif
    672       end do
     673        end do
    673674!END OLD METHOD
    674675      end do
    675      end do
    676 endif !readclouds
     676    end do
     677  endif !readclouds
    677678
    678679     !********* TEST ***************
     
    838839    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
    839840    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
    840   end subroutine verttransform
    841 
     841end subroutine verttransform
     842
  • src/wetdepo.f90

    rd6a0977 rf75967d  
    7373  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
    7474  integer :: ks, kp
    75   integer :: n1,n2, icbot,ictop, indcloud !TEST
     75!  integer :: n1,n2, icbot,ictop, indcloud !TEST
    7676  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    7777  real :: clouds_h ! cloud height for the specific grid point
     
    105105  ! Loop over all particles
    106106  !************************
    107 
    108107
    109108  blc_count=0
     
    185184    n=memind(2)
    186185    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    187     n=memind(1)
     186         n=memind(1)
    188187
    189188    if (ngrid.eq.0) then
     
    199198  ! scavenging is done
    200199
    201      if (clouds_v.le.1) goto 20
     200    if (clouds_v.le.1) goto 20
    202201 
    203202  ! 1) Parameterization of the the area fraction of the grid cell where the
     
    311310! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
    312311        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then
     312          inc_count=inc_count+1
    313313! if negative coefficients (turned off) set to zero for use in equation
    314314          if (weta_in(ks).lt.0.) weta_in(ks)=0.
     
    355355           !OLD
    356356          if (readclouds) then
    357            wetscav=S_i*(prec(1)/3.6E6)
     357            wetscav=S_i*(prec(1)/3.6E6)
    358358          else
    359           wetscav=S_i*(prec(1)/3.6E6)/clouds_h
     359            wetscav=S_i*(prec(1)/3.6E6)/clouds_h
    360360          endif
    361361
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG