- Timestamp:
- Jun 18, 2019, 1:59:41 PM (5 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug
- Children:
- 79d7d08
- Parents:
- 6741557 (diff), 941db73 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Location:
- src
- Files:
-
- 6 added
- 27 edited
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART.f90
r2753a5c r2eefa58 68 68 integer :: detectformat 69 69 70 71 72 ! Initialize arrays in com_mod73 !*****************************74 call com_mod_allocate_part(maxpart)75 70 76 77 71 ! Generate a large number of random numbers 78 72 !****************************************** … … 143 137 write(*,*) 'call readpaths' 144 138 endif 145 call readpaths (pathfile)139 call readpaths 146 140 147 141 if (verbosity.gt.1) then !show clock info … … 171 165 endif 172 166 endif 167 168 ! Initialize arrays in com_mod 169 !***************************** 170 call com_mod_allocate_part(maxpart) 171 173 172 174 173 ! Read the age classes to be used … … 453 452 endif 454 453 454 if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo' 455 455 call timemanager(metdata_format) 456 456 457 457 458 if (verbosity.gt.0) then … … 468 469 endif 469 470 end do 470 write (*,*) 'timemanager> call wetdepo'471 471 endif 472 472 -
src/FLEXPART_MPI.f90
r20963b1 r2eefa58 77 77 if (mp_measure_time) call mpif_mtime('flexpart',0) 78 78 79 ! Initialize arrays in com_mod 80 !***************************** 81 82 if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi) 83 84 79 85 80 ! Generate a large number of random numbers 86 81 !****************************************** … … 151 146 write(*,*) 'call readpaths' 152 147 endif 153 call readpaths (pathfile)148 call readpaths 154 149 155 150 if (verbosity.gt.1) then !show clock info … … 179 174 endif 180 175 endif 176 177 ! Initialize arrays in com_mod 178 !***************************** 179 180 if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi) 181 181 182 182 … … 413 413 end if ! (mpif_pid == 0) 414 414 415 if (mp_measure_time) call mpif_mtime('iotime', 0)415 if (mp_measure_time) call mpif_mtime('iotime',1) 416 416 417 417 if (verbosity.gt.0 .and. lroot) then -
src/com_mod.f90
re9e0f06 r2eefa58 18 18 19 19 implicit none 20 21 20 22 21 23 !**************************************************************** … … 69 71 70 72 real :: ctl,fine 71 integer :: ifine,iout,ipout,ipin,iflux,mdomainfill 73 integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac 72 74 integer :: mquasilag,nested_output,ind_source,ind_receptor 73 75 integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only … … 82 84 ! iout output options: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both 83 85 ! ipout particle dump options: 0 no, 1 every output interval, 2 only at end 86 ! ipoutfac increase particle dump interval by factor (default 1) 84 87 ! ipin read in particle positions from dumped file from a previous run 85 88 ! fine real(ifine) … … 121 124 ! lnetcdfout 1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input) 122 125 126 integer :: linversionout 127 ! linversionout 1 for one grid_time output file for each release containing all timesteps 128 123 129 integer :: nageclass,lage(maxageclass) 124 130 … … 128 134 129 135 logical :: gdomainfill 130 131 136 ! gdomainfill .T., if domain-filling is global, .F. if not 132 137 … … 174 179 real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass) 175 180 real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass) 176 real :: rm(maxspec),dryvel(maxspec) 181 real :: rm(maxspec),dryvel(maxspec),kao(maxspec) 177 182 real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec) 178 183 … … 359 364 real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice [kg/kg] 360 365 real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !combined [m3/m3] 361 366 ! RLT add pressure and dry air density 367 real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem) 368 real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem) 362 369 real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem) 363 370 real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem) … … 381 388 ! uupol,vvpol [m/s] wind components in polar stereographic projection 382 389 ! tt [K] temperature data 390 ! prs air pressure 383 391 ! qv specific humidity data 384 392 ! pv (pvu) potential vorticity … … 651 659 real :: receptorarea(maxreceptor) 652 660 real :: creceptor(maxreceptor,maxspec) 661 real, allocatable, dimension(:,:) :: creceptor0 653 662 character(len=16) :: receptorname(maxreceptor) 654 663 integer :: numreceptor … … 673 682 real, allocatable, dimension(:,:) :: xmass1 674 683 real, allocatable, dimension(:,:) :: xscav_frac1 684 685 ! Variables used for writing out interval averages for partoutput 686 !**************************************************************** 687 688 integer, allocatable, dimension(:) :: npart_av 689 real, allocatable, dimension(:) :: part_av_cartx,part_av_carty,part_av_cartz,part_av_z,part_av_topo 690 real, allocatable, dimension(:) :: part_av_pv,part_av_qv,part_av_tt,part_av_rho,part_av_tro,part_av_hmix 691 real, allocatable, dimension(:) :: part_av_uu,part_av_vv,part_av_energy 675 692 676 693 ! eso: Moved from timemanager … … 780 797 & idt(nmpart),itramem(nmpart),itrasplit(nmpart),& 781 798 & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),& 782 & xmass1(nmpart, maxspec),& 783 & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime 799 & xmass1(nmpart, maxspec)) ! ,& 800 ! & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime 801 802 if (ipout.eq.3) then 803 allocate(npart_av(nmpart),part_av_cartx(nmpart),part_av_carty(nmpart),& 804 & part_av_cartz(nmpart),part_av_z(nmpart),part_av_topo(nmpart)) 805 allocate(part_av_pv(nmpart),part_av_qv(nmpart),part_av_tt(nmpart),& 806 & part_av_rho(nmpart),part_av_tro(nmpart),part_av_hmix(nmpart)) 807 allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart)) 808 end if 784 809 785 810 786 811 allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),& 787 812 & vs(nmpart),ws(nmpart),cbt(nmpart)) 788 813 789 814 end subroutine com_mod_allocate_part 790 815 -
src/concoutput.f90
r20963b1 r2eefa58 72 72 real :: sp_fact 73 73 real :: outnum,densityoutrecept(maxreceptor),xl,yl 74 ! RLT 75 real :: densitydryrecept(maxreceptor) 76 real :: factor_dryrecept(maxreceptor) 74 77 75 78 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), … … 106 109 character(LEN=8),save :: file_stat='REPLACE' 107 110 logical :: ldates_file 111 logical :: lexist 108 112 integer :: ierr 109 113 character(LEN=100) :: dates_char … … 203 207 densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ & 204 208 rho(iix,jjy,kzz-1,mind)*dz2)/dz 209 ! RLT 210 densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ & 211 rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz 205 212 end do 206 213 end do … … 214 221 !densityoutrecept(i)=rho(iix,jjy,1,2) 215 222 densityoutrecept(i)=rho(iix,jjy,1,mind) 223 ! RLT 224 densitydryrecept(i)=rho_dry(iix,jjy,1,mind) 216 225 end do 217 226 227 ! RLT 228 ! conversion factor for output relative to dry air 229 factor_drygrid=densityoutgrid/densitydrygrid 230 factor_dryrecept=densityoutrecept/densitydryrecept 218 231 219 232 ! Output is different for forward and backward simulations … … 353 366 ! Concentration output 354 367 !********************* 355 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5) .or.(iout.eq.6)) then368 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 356 369 357 370 ! Wet deposition … … 614 627 end do 615 628 629 ! RLT Aug 2017 630 ! Write out conversion factor for dry air 631 inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist) 632 if (lexist) then 633 ! open and append 634 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',& 635 status='old',action='write',access='append') 636 else 637 ! create new 638 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',& 639 status='new',action='write') 640 endif 641 sp_count_i=0 642 sp_count_r=0 643 sp_fact=-1. 644 sp_zer=.true. 645 do kz=1,numzgrid 646 do jy=0,numygrid-1 647 do ix=0,numxgrid-1 648 if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then 649 if (sp_zer.eqv..true.) then ! first value not equal to one 650 sp_count_i=sp_count_i+1 651 sparse_dump_i(sp_count_i)= & 652 ix+jy*numxgrid+kz*numxgrid*numygrid 653 sp_zer=.false. 654 sp_fact=sp_fact*(-1.) 655 endif 656 sp_count_r=sp_count_r+1 657 sparse_dump_r(sp_count_r)= & 658 sp_fact*factor_drygrid(ix,jy,kz) 659 else ! factor is one 660 sp_zer=.true. 661 endif 662 end do 663 end do 664 end do 665 write(unitoutfactor) sp_count_i 666 write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i) 667 write(unitoutfactor) sp_count_r 668 write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r) 669 close(unitoutfactor) 670 671 616 672 if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal 617 673 if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ & … … 640 696 endif 641 697 642 698 ! RLT Aug 2017 699 ! Write out conversion factor for dry air 700 if (numreceptor.gt.0) then 701 inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist) 702 if (lexist) then 703 ! open and append 704 open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',& 705 status='old',action='write',access='append') 706 else 707 ! create new 708 open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',& 709 status='new',action='write') 710 endif 711 write(unitoutfactor) itime 712 write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor) 713 close(unitoutfactor) 714 endif 643 715 644 716 ! Reinitialization of grid -
src/concoutput_mpi.f90
r20963b1 r6741557 364 364 ! Concentration output 365 365 !********************* 366 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5) .or.(iout.eq.6)) then366 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 367 367 368 368 ! Wet deposition -
src/concoutput_nest.f90
r7d02c2f r2eefa58 70 70 real :: sp_fact 71 71 real :: outnum,densityoutrecept(maxreceptor),xl,yl 72 ! RLT 73 real :: densitydryrecept(maxreceptor) 74 real :: factor_dryrecept(maxreceptor) 72 75 73 76 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), … … 96 99 character :: adate*8,atime*6 97 100 character(len=3) :: anspec 101 logical :: lexist 98 102 integer :: mind 99 103 ! mind eso:added to ensure identical results between 2&3-fields versions … … 164 168 densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ & 165 169 rho(iix,jjy,kzz-1,mind)*dz2)/dz 170 ! RLT 171 densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ & 172 rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz 166 173 end do 167 174 end do … … 175 182 !densityoutrecept(i)=rho(iix,jjy,1,2) 176 183 densityoutrecept(i)=rho(iix,jjy,1,mind) 184 ! RLT 185 densitydryrecept(i)=rho_dry(iix,jjy,1,mind) 177 186 end do 178 187 188 ! RLT 189 ! conversion factor for output relative to dry air 190 factor_drygrid=densityoutgrid/densitydrygrid 191 factor_dryrecept=densityoutrecept/densitydryrecept 179 192 180 193 ! Output is different for forward and backward simulations … … 551 564 end do 552 565 566 ! RLT Aug 2017 567 ! Write out conversion factor for dry air 568 inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist) 569 if (lexist) then 570 ! open and append 571 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',& 572 status='old',action='write',access='append') 573 else 574 ! create new 575 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',& 576 status='new',action='write') 577 endif 578 sp_count_i=0 579 sp_count_r=0 580 sp_fact=-1. 581 sp_zer=.true. 582 do kz=1,numzgrid 583 do jy=0,numygridn-1 584 do ix=0,numxgridn-1 585 if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then 586 if (sp_zer.eqv..true.) then ! first value not equal to one 587 sp_count_i=sp_count_i+1 588 sparse_dump_i(sp_count_i)= & 589 ix+jy*numxgridn+kz*numxgridn*numygridn 590 sp_zer=.false. 591 sp_fact=sp_fact*(-1.) 592 endif 593 sp_count_r=sp_count_r+1 594 sparse_dump_r(sp_count_r)= & 595 sp_fact*factor_drygrid(ix,jy,kz) 596 else ! factor is one 597 sp_zer=.true. 598 endif 599 end do 600 end do 601 end do 602 write(unitoutfactor) sp_count_i 603 write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i) 604 write(unitoutfactor) sp_count_r 605 write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r) 606 close(unitoutfactor) 553 607 554 608 -
src/concoutput_surf.f90
r16b61a5 r2eefa58 72 72 real :: sp_fact 73 73 real :: outnum,densityoutrecept(maxreceptor),xl,yl 74 ! RLT 75 real :: densitydryrecept(maxreceptor) 76 real :: factor_dryrecept(maxreceptor) 74 77 75 78 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), … … 101 104 character :: adate*8,atime*6 102 105 character(len=3) :: anspec 106 logical :: lexist 103 107 104 108 … … 180 184 densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ & 181 185 rho(iix,jjy,kzz-1,2)*dz2)/dz 186 ! RLT 187 densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ & 188 rho_dry(iix,jjy,kzz-1,2)*dz2)/dz 182 189 end do 183 190 end do … … 190 197 jjy=max(min(nint(yl),nymin1),0) 191 198 densityoutrecept(i)=rho(iix,jjy,1,2) 199 ! RLT 200 densitydryrecept(i)=rho_dry(iix,jjy,1,2) 192 201 end do 193 202 203 ! RLT 204 ! conversion factor for output relative to dry air 205 factor_drygrid=densityoutgrid/densitydrygrid 206 factor_dryrecept=densityoutrecept/densitydryrecept 194 207 195 208 ! Output is different for forward and backward simulations … … 596 609 end do 597 610 611 ! RLT Aug 2017 612 ! Write out conversion factor for dry air 613 inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist) 614 if (lexist) then 615 ! open and append 616 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',& 617 status='old',action='write',access='append') 618 else 619 ! create new 620 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',& 621 status='new',action='write') 622 endif 623 sp_count_i=0 624 sp_count_r=0 625 sp_fact=-1. 626 sp_zer=.true. 627 do kz=1,1 628 do jy=0,numygrid-1 629 do ix=0,numxgrid-1 630 if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then 631 if (sp_zer.eqv..true.) then ! first value not equal to one 632 sp_count_i=sp_count_i+1 633 sparse_dump_i(sp_count_i)= & 634 ix+jy*numxgrid+kz*numxgrid*numygrid 635 sp_zer=.false. 636 sp_fact=sp_fact*(-1.) 637 endif 638 sp_count_r=sp_count_r+1 639 sparse_dump_r(sp_count_r)= & 640 sp_fact*factor_drygrid(ix,jy,kz) 641 else ! factor is one 642 sp_zer=.true. 643 endif 644 end do 645 end do 646 end do 647 write(unitoutfactor) sp_count_i 648 write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i) 649 write(unitoutfactor) sp_count_r 650 write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r) 651 close(unitoutfactor) 652 653 598 654 if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal 599 655 if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ & … … 622 678 endif 623 679 624 680 ! RLT Aug 2017 681 ! Write out conversion factor for dry air 682 if (numreceptor.gt.0) then 683 inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist) 684 if (lexist) then 685 ! open and append 686 open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',& 687 status='old',action='write',access='append') 688 else 689 ! create new 690 open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',& 691 status='new',action='write') 692 endif 693 write(unitoutfactor) itime 694 write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor) 695 close(unitoutfactor) 696 endif 625 697 626 698 ! Reinitialization of grid -
src/concoutput_surf_nest.f90
r6a678e3 r2eefa58 70 70 real :: sp_fact 71 71 real :: outnum,densityoutrecept(maxreceptor),xl,yl 72 ! RLT 73 real :: densitydryrecept(maxreceptor) 74 real :: factor_dryrecept(maxreceptor) 72 75 73 76 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), … … 96 99 character :: adate*8,atime*6 97 100 character(len=3) :: anspec 98 101 logical :: lexist 99 102 100 103 ! Determine current calendar date, needed for the file name … … 159 162 densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ & 160 163 rho(iix,jjy,kzz-1,2)*dz2)/dz 164 ! RLT 165 densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ & 166 rho_dry(iix,jjy,kzz-1,2)*dz2)/dz 161 167 end do 162 168 end do … … 169 175 jjy=max(min(nint(yl),nymin1),0) 170 176 densityoutrecept(i)=rho(iix,jjy,1,2) 177 ! RLT 178 densitydryrecept(i)=rho_dry(iix,jjy,1,2) 171 179 end do 172 180 181 ! RLT 182 ! conversion factor for output relative to dry air 183 factor_drygrid=densityoutgrid/densitydrygrid 184 factor_dryrecept=densityoutrecept/densitydryrecept 173 185 174 186 ! Output is different for forward and backward simulations … … 317 329 write(unitoutgrid) sp_count_r 318 330 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 319 write(unitoutgrid) sp_count_r320 write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)331 ! write(unitoutgrid) sp_count_r 332 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 321 333 322 334 ! Dry deposition … … 354 366 write(unitoutgrid) sp_count_r 355 367 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 356 write(unitoutgrid) sp_count_r357 write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)368 ! write(unitoutgrid) sp_count_r 369 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 358 370 359 371 … … 399 411 write(unitoutgrid) sp_count_r 400 412 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 401 write(unitoutgrid) sp_count_r402 write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)413 ! write(unitoutgrid) sp_count_r 414 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 403 415 else 404 416 … … 440 452 write(unitoutgrid) sp_count_r 441 453 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 442 write(unitoutgrid) sp_count_r443 write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)454 ! write(unitoutgrid) sp_count_r 455 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 444 456 endif ! surf_only 445 457 … … 487 499 write(unitoutgridppt) sp_count_r 488 500 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 489 write(unitoutgridppt) sp_count_r490 write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)501 ! write(unitoutgridppt) sp_count_r 502 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 491 503 492 504 … … 526 538 write(unitoutgridppt) sp_count_r 527 539 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 528 write(unitoutgridppt) sp_count_r529 write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)540 ! write(unitoutgridppt) sp_count_r 541 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 530 542 531 543 … … 570 582 write(unitoutgridppt) sp_count_r 571 583 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 572 write(unitoutgridppt) sp_count_r573 write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)584 ! write(unitoutgridppt) sp_count_r 585 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 574 586 else 575 587 … … 611 623 write(unitoutgridppt) sp_count_r 612 624 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 613 write(unitoutgridppt) sp_count_r614 write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)625 ! write(unitoutgridppt) sp_count_r 626 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 615 627 endif ! surf_only 616 628 … … 624 636 625 637 end do 638 639 ! RLT Aug 2017 640 ! Write out conversion factor for dry air 641 inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist) 642 if (lexist) then 643 ! open and append 644 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',& 645 status='old',action='write',access='append') 646 else 647 ! create new 648 open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',& 649 status='new',action='write') 650 endif 651 sp_count_i=0 652 sp_count_r=0 653 sp_fact=-1. 654 sp_zer=.true. 655 do kz=1,1 656 do jy=0,numygridn-1 657 do ix=0,numxgridn-1 658 if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then 659 if (sp_zer.eqv..true.) then ! first value not equal to one 660 sp_count_i=sp_count_i+1 661 sparse_dump_i(sp_count_i)= & 662 ix+jy*numxgridn+kz*numxgridn*numygridn 663 sp_zer=.false. 664 sp_fact=sp_fact*(-1.) 665 endif 666 sp_count_r=sp_count_r+1 667 sparse_dump_r(sp_count_r)= & 668 sp_fact*factor_drygrid(ix,jy,kz) 669 else ! factor is one 670 sp_zer=.true. 671 endif 672 end do 673 end do 674 end do 675 write(unitoutfactor) sp_count_i 676 write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i) 677 write(unitoutfactor) sp_count_r 678 write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r) 679 close(unitoutfactor) 626 680 627 681 -
src/getfields.f90
r6ecb30a r2eefa58 87 87 real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) 88 88 real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) 89 ! RLT added partial pressure water vapor 90 real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem) 91 integer :: kz, ix 92 character(len=100) :: rowfmt 89 93 90 94 integer :: indmin = 1 … … 153 157 40 indmin=indj 154 158 155 if (WETBKDEP) then156 157 endif159 if (WETBKDEP) then 160 call writeprecip(itime,memind(1)) 161 endif 158 162 159 163 else … … 204 208 60 indmin=indj 205 209 206 if (WETBKDEP) then 207 call writeprecip(itime,memind(1)) 208 endif 209 210 endif 210 if (WETBKDEP) then 211 call writeprecip(itime,memind(1)) 212 endif 213 214 end if 215 216 ! RLT calculate dry air density 217 pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv) 218 rho_dry=(prs-pwater)/(r_air*tt) 219 220 ! test density 221 ! write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))' 222 ! if(itime.eq.0) then 223 ! open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write') 224 ! do kz=1,nzmax 225 ! do ix=1,nxmax 226 ! write(500,fmt=rowfmt) rho_dry(ix,:,kz,1) 227 ! end do 228 ! end do 229 ! close(500) 230 ! open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write') 231 ! do kz=1,nzmax 232 ! do ix=1,nxmax 233 ! write(500,fmt=rowfmt) rho(ix,:,kz,1) 234 ! end do 235 ! end do 236 ! close(500) 237 ! endif 211 238 212 239 lwindinterv=abs(memtime(2)-memtime(1)) -
src/init_domainfill.f90
rb5127f9 r0a94e13 86 86 endif 87 87 endif 88 89 ! Exit here if resuming a run from particle dump 90 !*********************************************** 91 if (gdomainfill.and.ipin.ne.0) return 88 92 89 93 ! Do not release particles twice (i.e., not at both in the leftmost and rightmost … … 414 418 !*************************************************************************** 415 419 416 if ( ipin.eq.1) then420 if ((ipin.eq.1).and.(.not.gdomainfill)) then 417 421 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 418 422 form='unformatted') -
src/init_domainfill_mpi.f90
rb5127f9 r328fdf9 110 110 endif 111 111 112 ! Exit here if resuming a run from particle dump 113 !*********************************************** 114 if (gdomainfill.and.ipin.ne.0) return 115 112 116 ! Do not release particles twice (i.e., not at both in the leftmost and rightmost 113 117 ! grid cell) for a global domain … … 213 217 colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy) 214 218 colmasstotal=colmasstotal+colmass(ix,jy) 215 216 219 end do 217 220 end do … … 466 469 467 470 ! eso TODO: only needed for root process 468 if ( ipin.eq.1) then471 if ((ipin.eq.1).and.(.not.gdomainfill)) then 469 472 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 470 473 form='unformatted') … … 474 477 endif 475 478 476 numpart = numpart/mp_partgroup_np 477 if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1 478 479 else ! Allocate dummy arrays for receiving processes 480 allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),& 481 & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),& 482 & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),& 483 & xmass1_tmp(nullsize, nullsize)) 479 if (ipin.eq.0) then 480 numpart = numpart/mp_partgroup_np 481 if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1 482 end if 483 484 else ! Allocate dummy arrays for receiving processes 485 if (ipin.eq.0) then 486 allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),& 487 & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),& 488 & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),& 489 & xmass1_tmp(nullsize, nullsize)) 490 end if 484 491 485 end if ! end if(lroot) 492 end if ! end if(lroot) 493 486 494 487 495 488 496 ! Distribute particles to other processes (numpart is 'per-process', not total) 489 call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr) 490 ! eso TODO: xmassperparticle: not necessary to send 491 call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr)492 call mpif_send_part_properties(numpart)497 ! Only if not restarting from previous run 498 if (ipin.eq.0) then 499 call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr) 500 call mpif_send_part_properties(npart(1)/mp_partgroup_np) 493 501 494 502 ! Deallocate the temporary arrays used for all particles 495 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&503 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,& 496 504 & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp) 505 end if 497 506 498 507 -
src/makefile
r7123c70 r2eefa58 118 118 OBJECTS_SERIAL = \ 119 119 releaseparticles.o partoutput.o \ 120 partoutput_average.o \ 120 121 conccalc.o \ 121 122 init_domainfill.o concoutput.o \ … … 127 128 redist.o \ 128 129 concoutput_surf.o concoutput_surf_nest.o \ 130 concoutput_inversion_nest.o \ 131 concoutput_inversion.o \ 129 132 getfields.o \ 130 133 readwind_ecmwf.o … … 132 135 ## For MPI version 133 136 OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \ 134 conccalc_mpi.o \137 partoutput_average_mpi.o conccalc_mpi.o \ 135 138 init_domainfill_mpi.o concoutput_mpi.o \ 136 139 timemanager_mpi.o FLEXPART_MPI.o \ … … 149 152 advance.o initialize.o \ 150 153 writeheader.o writeheader_txt.o \ 151 writeprecip.o \154 partpos_average.o writeprecip.o \ 152 155 writeheader_surf.o assignland.o\ 153 156 part0.o gethourlyOH.o\ … … 199 202 drydepokernel_nest.o zenithangle.o \ 200 203 ohreaction.o getvdep_nests.o \ 201 initial_cond_calc.o initial_cond_output.o \204 initial_cond_calc.o initial_cond_output.o initial_cond_output_inversion.o \ 202 205 dynamic_viscosity.o get_settling.o \ 203 206 initialize_cbl_vel.o re_initialize_particle.o \ … … 271 274 conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o 272 275 concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o 276 concoutput_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o 273 277 concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \ 274 278 unc_mod.o mean_mod.o 275 279 concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o 280 concoutput_inversion_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o 276 281 concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \ 277 282 unc_mod.o mean_mod.o … … 320 325 initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o 321 326 initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o 327 initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o 322 328 initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o 323 329 initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o … … 348 354 part0.o: par_mod.o 349 355 partdep.o: par_mod.o 356 partpos_average.o: com_mod.o par_mod.o 350 357 partoutput.o: com_mod.o par_mod.o 358 partoutput_average.o: com_mod.o par_mod.o 359 partoutput_average_mpi.o: com_mod.o par_mod.o mpi_mod.o 351 360 partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o 352 361 partoutput_short.o: com_mod.o par_mod.o -
src/mpi_mod.f90
r0ecc1fe r0c8c7f2 88 88 ! Variables for MPI processes in the 'particle' group 89 89 integer, allocatable, dimension(:) :: mp_partgroup_rank 90 integer, allocatable, dimension(:) :: npart_per_process 90 91 integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np 91 92 … … 125 126 ! mp_time_barrier Measure MPI barrier time 126 127 ! mp_exact_numpart Use an extra MPI communication to give the exact number of particles 127 ! to standard output (this does *not* otherwise affect the simulation)128 ! to standard output (this does not otherwise affect the simulation) 128 129 logical, parameter :: mp_dbg_mode = .false. 129 130 logical, parameter :: mp_dev_mode = .false. … … 190 191 ! mp_np number of running processes, decided at run-time 191 192 !*********************************************************************** 192 use par_mod, only: maxpart, numwfmem, dep_prec 193 use com_mod, only: mpi_mode, verbosity 193 use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec 194 use com_mod, only: mpi_mode, verbosity, creceptor0 194 195 195 196 implicit none … … 337 338 338 339 ! Set maxpart per process 339 ! eso08/2016: Increase maxpart per process, in case of unbalanced distribution340 ! ESO 08/2016: Increase maxpart per process, in case of unbalanced distribution 340 341 maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np)) 341 342 if (mp_np == 1) maxpart_mpi = maxpart … … 365 366 end if 366 367 368 ! Allocate array for number of particles per process 369 allocate(npart_per_process(0:mp_partgroup_np-1)) 370 371 ! Write whether MPI_IN_PLACE is used or not 372 #ifdef USE_MPIINPLACE 373 if (lroot) write(*,*) 'Using MPI_IN_PLACE operations' 374 #else 375 if (lroot) allocate(creceptor0(maxreceptor,maxspec)) 376 if (lroot) write(*,*) 'Not using MPI_IN_PLACE operations' 377 #endif 367 378 goto 101 368 379 … … 559 570 ! invalid particles at the end of the arrays 560 571 561 601 do i=num _part, 1, -1572 601 do i=numpart, 1, -1 562 573 if (itra1(i).eq.-999999999) then 563 574 numpart=numpart-1 … … 598 609 integer :: i,jj,nn, num_part=1,m,imin, num_trans 599 610 logical :: first_iter 600 integer,allocatable,dimension(:) :: numparticles_mpi,idx_arr611 integer,allocatable,dimension(:) :: idx_arr 601 612 real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this 602 613 … … 607 618 ! All processes exchange information on number of particles 608 619 !**************************************************************************** 609 allocate(numparticles_mpi(0:mp_partgroup_np-1), & 610 &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1)) 611 612 call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, & 620 allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1)) 621 622 call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, & 613 623 & 1, MPI_INTEGER, mp_comm_used, mp_ierr) 614 624 … … 616 626 ! Sort from lowest to highest 617 627 ! Initial guess: correct order 618 sorted(:) = n umparticles_mpi(:)628 sorted(:) = npart_per_process(:) 619 629 do i=0, mp_partgroup_np-1 620 630 idx_arr(i) = i 621 631 end do 632 633 ! Do not rebalance particles for ipout=3 634 if (ipout.eq.3) return 622 635 623 636 ! For each successive element in index array, see if a lower value exists … … 645 658 m=mp_partgroup_np-1 ! index for last sorted process (most particles) 646 659 do i=0,mp_partgroup_np/2-1 647 num_trans = n umparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))660 num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i)) 648 661 if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then 649 if ( n umparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&650 & real(num_trans)/real(n umparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then662 if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.& 663 & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then 651 664 ! DBG 652 ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, n umparticles_mpi', &653 ! &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, n umparticles_mpi665 ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', & 666 ! &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process 654 667 ! DBG 655 668 call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2) … … 659 672 end do 660 673 661 deallocate( numparticles_mpi,idx_arr, sorted)674 deallocate(idx_arr, sorted) 662 675 663 676 end subroutine mpif_calculate_part_redist … … 1961 1974 if (readclouds) then 1962 1975 j=j+1 1963 call MPI_Irecv(ctwc(:,:,mind),d2s1 ,mp_sp,id_read,MPI_ANY_TAG,&1976 call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,& 1964 1977 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1965 1978 if (mp_ierr /= 0) goto 600 … … 2326 2339 if (readclouds) then 2327 2340 j=j+1 2328 call MPI_Irecv(ctwcn(:,:,mind,k),d2s1 ,mp_sp,id_read,MPI_ANY_TAG,&2341 call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,& 2329 2342 &MPI_COMM_WORLD,reqs(j),mp_ierr) 2330 2343 if (mp_ierr /= 0) goto 600 … … 2462 2475 end if 2463 2476 2477 ! Receptor concentrations 2478 if (lroot) then 2479 call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, & 2480 & mp_comm_used,mp_ierr) 2481 if (mp_ierr /= 0) goto 600 2482 else 2483 call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, & 2484 & mp_comm_used,mp_ierr) 2485 end if 2486 2464 2487 #else 2465 2488 2466 2489 call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, & 2467 2490 & mp_comm_used, mp_ierr) 2491 if (mp_ierr /= 0) goto 600 2468 2492 if (lroot) gridunc = gridunc0 2493 2494 call MPI_Reduce(creceptor, creceptor0,rcpt_size,mp_sp,MPI_SUM,id_root, & 2495 & mp_comm_used,mp_ierr) 2496 if (mp_ierr /= 0) goto 600 2497 if (lroot) creceptor = creceptor0 2469 2498 2470 2499 #endif … … 2482 2511 end if 2483 2512 2484 ! Receptor concentrations2485 if (lroot) then2486 call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &2487 & mp_comm_used,mp_ierr)2488 if (mp_ierr /= 0) goto 6002489 else2490 call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &2491 & mp_comm_used,mp_ierr)2492 end if2493 2513 2494 2514 if (mp_measure_time) call mpif_mtime('commtime',1) … … 2700 2720 end if 2701 2721 2702 2703 2704 2705 2706 2707 2708 2709 2710 2711 2712 2713 2714 2722 case ('readwind') 2723 if (imode.eq.0) then 2724 call cpu_time(mp_readwind_time_beg) 2725 mp_readwind_wtime_beg = mpi_wtime() 2726 else 2727 call cpu_time(mp_readwind_time_end) 2728 mp_readwind_wtime_end = mpi_wtime() 2729 2730 mp_readwind_time_total = mp_readwind_time_total + & 2731 &(mp_readwind_time_end - mp_readwind_time_beg) 2732 mp_readwind_wtime_total = mp_readwind_wtime_total + & 2733 &(mp_readwind_wtime_end - mp_readwind_wtime_beg) 2734 end if 2715 2735 2716 2736 case ('commtime') … … 2788 2808 write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',& 2789 2809 & mp_getfields_time_total 2790 write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&2791 & mp_readwind_wtime_total2792 write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&2793 & mp_readwind_time_total2810 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',& 2811 ! & mp_readwind_wtime_total 2812 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',& 2813 ! & mp_readwind_time_total 2794 2814 write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',& 2795 2815 & mp_io_wtime_total -
src/netcdf_output_mod.f90
r4ad96c5 r0a94e13 93 93 character(len=255), parameter :: institution = 'NILU' 94 94 95 integer :: tpointer 95 integer :: tpointer=0 96 96 character(len=255) :: ncfname, ncfnamen 97 97 -
src/outg_mod.f90
r4c64400 r2eefa58 37 37 real,allocatable, dimension (:,:,:) :: areanorth 38 38 real,allocatable, dimension (:,:,:) :: densityoutgrid 39 real,allocatable, dimension (:,:,:) :: densitydrygrid ! added RLT 40 real,allocatable, dimension (:,:,:) :: factor_drygrid ! added RLT 39 41 real,allocatable, dimension (:,:,:) :: factor3d 40 42 real,allocatable, dimension (:,:,:) :: grid -
src/outgrid_init.f90
rd2a5a83 r2eefa58 268 268 if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc' 269 269 allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, & 270 0:max(numygrid,numygridn)-1,numzgrid),stat=stat) 271 if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc' 272 ! RLT 273 allocate(densitydrygrid(0:max(numxgrid,numxgridn)-1, & 274 0:max(numygrid,numygridn)-1,numzgrid),stat=stat) 275 if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc' 276 allocate(factor_drygrid(0:max(numxgrid,numxgridn)-1, & 270 277 0:max(numygrid,numygridn)-1,numzgrid),stat=stat) 271 278 if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc' -
src/par_mod.f90
rdf96ea65 r2eefa58 30 30 ! Update 15 August 2013 IP * 31 31 ! * 32 ! ESO 2016: *33 ! GFS specific parameters moved to gfs_mod.f90 *34 ! ECMWF specific parameters moved to ecmwf_mod.f90 *35 32 ! * 36 33 !******************************************************************************* … … 80 77 real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81 81 78 real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4 79 ! additional constants RLT Aug-2017 80 real,parameter :: rgas=8.31447 81 real,parameter :: r_water=461.495 82 82 83 83 ! pi number "pi" … … 89 89 ! kappa exponent of formula for potential temperature 90 90 ! vonkarman von Karman constant 91 ! rgas universal gas constant [J/mol/K] 92 ! r_water specific gas constant for water vapor [J/kg/K] 91 93 92 94 real,parameter :: karman=0.40, href=15., convke=2.0 … … 280 282 281 283 integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1 282 integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93 284 integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105 283 285 integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96 284 286 integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1 … … 290 292 integer,parameter :: unitboundcond=89 291 293 integer,parameter :: unittmp=101 294 ! RLT 295 integer,parameter :: unitoutfactor=102 292 296 293 297 !****************************************************** -
src/partoutput.f90
rd2a5a83 r0a94e13 71 71 !************************************** 72 72 73 if (ipout.eq.1 ) then73 if (ipout.eq.1.or.ipout.eq.3) then 74 74 open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// & 75 75 atime,form='unformatted') -
src/partoutput_mpi.f90
rd2a5a83 r0a94e13 78 78 !************************************** 79 79 80 if (ipout.eq.1 ) then80 if (ipout.eq.1.or.ipout.eq.3) then 81 81 open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// & 82 82 atime,form='unformatted',status=file_stat,position='append') -
src/readcommand.f90
r20963b1 r2eefa58 50 50 ! ipin 1 continue simulation with dumped particle data, 0 no * 51 51 ! ipout 0 no particle dump, 1 every output time, 3 only at end* 52 ! ipoutfac increase particle dump interval by factor (default 1) * 52 53 ! itsplit [s] time constant for particle splitting * 53 54 ! loutaver [s] concentration output is an average over loutaver * … … 97 98 iout, & 98 99 ipout, & 100 ipoutfac, & 99 101 lsubgrid, & 100 102 lconvection, & … … 112 114 surf_only, & 113 115 cblflag, & 116 linversionout, & 114 117 ohfields_path 115 118 … … 129 132 iout=3 130 133 ipout=0 134 ipoutfac=1 131 135 lsubgrid=1 132 136 lconvection=1 … … 144 148 surf_only=0 145 149 cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine 150 linversionout=0 146 151 ohfields_path="../../flexin/" 147 152 … … 411 416 !********************************************************************** 412 417 413 if ((iout.lt.1).or.(iout.gt. 6)) then418 if ((iout.lt.1).or.(iout.gt.5)) then 414 419 write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' 415 420 write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR #### ' … … 471 476 write(*,*) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR ####' 472 477 write(*,*) '#### EACH RELEASE IS FORBIDDEN ! ####' 478 stop 479 endif 480 481 ! Inversion output format only for backward runs 482 !***************************************************************************** 483 484 if ((linversionout.eq.1).and.(ldirect.eq.1)) then 485 write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' 486 write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR ####' 487 write(*,*) '#### BACKWARD RUNS ####' 473 488 stop 474 489 endif … … 507 522 !**************************************************************** 508 523 509 if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2) ) then524 if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then 510 525 write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' 511 write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!#### '526 write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3! #### ' 512 527 stop 513 528 endif -
src/readpaths.f90
r62e65c7 r2eefa58 20 20 !********************************************************************** 21 21 22 subroutine readpaths !(pathfile)22 subroutine readpaths 23 23 24 24 !***************************************************************************** -
src/readwind_gfs.f90
rdb91eb7 r5d74ed9 83 83 84 84 ! NCEP 85 integer :: numpt,numpu,numpv,numpw,numprh 85 integer :: numpt,numpu,numpv,numpw,numprh,numpclwch 86 86 real :: help, temp, ew 87 87 real :: elev … … 134 134 numpw=0 135 135 numprh=0 136 numpclwch=0 136 137 ifield=0 137 138 10 ifield=ifield+1 … … 557 558 endif 558 559 ! SEC & IP 12/2018 read GFS clouds 559 if(isec1(6).eq.153) then !! CLWCR Cloud liquid water content [kg/kg] 560 clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1) 560 if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg] 561 if((i.eq.0).and.(j.eq.0)) then 562 do ii=1,nuvz 563 if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii 564 end do 565 endif 566 help=zsec4(nxfield*(ny-j-1)+i+1) 567 if(i.le.i180) then 568 clwch(i179+i,j,numpclwch,n)=help 569 else 570 clwch(i-i181,j,numpclwch,n)=help 571 endif 561 572 readclouds=.true. 562 573 sumclouds=.true. 574 ! readclouds=.false. 575 ! sumclouds=.false. 563 576 endif 564 577 -
src/releaseparticles.f90
r75a4ded r7873bf7 114 114 average_timecorrect=0. 115 115 do k=1,nspec 116 if (zpoint1(i).gt.0.5) then ! point source 116 if(abs(xpoint2(i)-xpoint1(i)).lt.1.E-4.and.abs(ypoint2(i)-ypoint1(i)).lt.1.E-4) then 117 ! if (zpoint1(i).gt.0.5) then ! point source 117 118 timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek) 118 119 else ! area source -
src/timemanager.f90
rc7d1052 r2eefa58 408 408 #endif 409 409 else 410 if (linversionout.eq.1) then 411 call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc) 412 if (verbosity.eq.1) then 413 print*,'called concoutput_inversion' 414 call system_clock(count_clock) 415 write(*,*) 'system clock',count_clock - count_clock0 416 endif 417 else 410 418 call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc) 419 endif 411 420 if (verbosity.eq.1) then 412 421 print*,'called concoutput_surf ' … … 422 431 call concoutput_nest(itime,outnum) 423 432 else 433 if(linversionout.eq.1) then 434 call concoutput_inversion_nest(itime,outnum) 435 else 424 436 call concoutput_surf_nest(itime,outnum) 437 endif 425 438 endif 426 439 else … … 451 464 45 format(i13,' Seconds simulated: ',i13, ' Particles: Uncertainty: ',3f7.3) 452 465 46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles') 453 if (ipout.ge.1) call partoutput(itime) ! dump particle positions 466 if (ipout.ge.1) then 467 if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions 468 if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions 469 endif 454 470 loutnext=loutnext+loutstep 455 471 loutstart=loutnext-loutaver/2 … … 609 625 ! write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j) 610 626 627 ! Calculate average position for particle dump output 628 !**************************************************** 629 630 if (ipout.eq.3) call partpos_average(itime,j) 631 632 611 633 ! Calculate the gross fluxes across layer interfaces 612 634 !*************************************************** … … 730 752 if (ipout.eq.2) call partoutput(itime) ! dump particle positions 731 753 732 if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field 754 if (linit_cond.ge.1) then 755 if(linversionout.eq.1) then 756 call initial_cond_output_inversion(itime) ! dump initial cond. field 757 else 758 call initial_cond_output(itime) ! dump initial cond. fielf 759 endif 760 endif 733 761 734 762 !close(104) -
src/timemanager_mpi.f90
r20963b1 r0c8c7f2 113 113 integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 114 114 ! integer :: ksp 115 integer :: ip 115 integer :: ip,irec 116 116 integer :: loutnext,loutstart,loutend 117 117 integer :: ix,jy,ldeltat,itage,nage,idummy … … 129 129 ! Measure time spent in timemanager 130 130 if (mp_measure_time) call mpif_mtime('timemanager',0) 131 131 132 132 133 ! First output for time 0 … … 532 533 end if 533 534 534 else ! :TODO: check for zeroing in the netcdf module535 else 535 536 call concoutput_surf_nest(itime,outnum) 536 537 end if … … 593 594 46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles') 594 595 if (ipout.ge.1) then 596 if (mp_measure_time) call mpif_mtime('iotime',0) 597 irec=0 595 598 do ip=0, mp_partgroup_np-1 596 if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions 599 if (ip.eq.mp_partid) then 600 if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions 601 if (ipout.eq.3) call partoutput_average(itime,irec) ! dump particle positions 602 endif 603 if (ipout.eq.3) irec=irec+npart_per_process(ip) 597 604 call mpif_mpi_barrier 598 605 end do 606 if (mp_measure_time) call mpif_mtime('iotime',1) 599 607 end if 600 608 … … 757 765 if (mp_measure_time) call mpif_mtime('advance',1) 758 766 767 ! Calculate average position for particle dump output 768 !**************************************************** 769 770 if (ipout.eq.3) call partpos_average(itime,j) 771 759 772 760 773 ! Calculate the gross fluxes across layer interfaces … … 895 908 do ip=0, mp_partgroup_np-1 896 909 if (ip.eq.mp_partid) then 897 !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid910 if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid 898 911 call partoutput(itime) ! dump particle positions 899 912 end if -
src/verttransform_ecmwf.f90
r79e0349 r2eefa58 73 73 use com_mod 74 74 use cmapf_mod, only: cc2gll 75 ! use mpi_mod76 75 77 76 implicit none … … 82 81 real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev 83 82 real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv 83 ! RLT added pressure 84 real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: prsh 84 85 real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv 85 86 real,dimension(0:nymax-1) :: cosf … … 219 220 !************************* 220 221 221 222 222 do jy=0,nymin1 223 223 do ix=0,nxmin1 … … 230 230 wzlev(:,:,1)=0. 231 231 rhoh(:,:,1)=pold/(r_air*tvold) 232 ! RLT add pressure 233 prsh(:,:,1)=ps(:,:,1,n) 232 234 233 235 … … 237 239 do kz=2,nuvz 238 240 pint=akz(kz)+bkz(kz)*ps(:,:,1,n) 241 ! RLT add pressure 242 prsh(:,:,kz)=pint 239 243 tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n)) 240 244 rhoh(:,:,kz)=pint(:,:)/(r_air*tv) … … 288 292 pv(:,:,1,n)=pvh(:,:,1) 289 293 rho(:,:,1,n)=rhoh(:,:,1) 294 ! RLT add pressure 295 prs(:,:,1,n)=prsh(:,:,1) 290 296 291 297 uu(:,:,nz,n)=uuh(:,:,nuvz) … … 301 307 pv(:,:,nz,n)=pvh(:,:,nuvz) 302 308 rho(:,:,nz,n)=rhoh(:,:,nuvz) 303 309 ! RLT 310 prs(:,:,nz,n)=prsh(:,:,nuvz) 304 311 305 312 kmin=2 … … 321 328 pv(ix,jy,iz,n)=pv(ix,jy,nz,n) 322 329 rho(ix,jy,iz,n)=rho(ix,jy,nz,n) 330 ! RLT 331 prs(ix,jy,iz,n)=prs(ix,jy,nz,n) 323 332 else 324 333 innuvz: do kz=idx(ix,jy),nuvz … … 354 363 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz 355 364 rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz 365 ! RLT add pressure 366 prs(ix,jy,iz,n)=(prsh(ix,jy,kz-1)*dz2+prsh(ix,jy,kz)*dz1)/dz 356 367 endif 357 368 enddo … … 654 665 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 655 666 656 do kz=nz, 1,-1 !go Bottom up!667 do kz=nz,2,-1 !go Bottom up! 657 668 if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud 658 669 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) -
src/verttransform_gfs.f90
rdb91eb7 r437c545 548 548 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 549 549 550 do kz=nz, 1,-1 !go Bottom up!550 do kz=nz,2,-1 !go Bottom up! 551 551 if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud 552 552 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 553 553 clouds(ix,jy,kz,n)=1 ! is a cloud 554 554 if (lsp.ge.convp) then 555 clouds(ix,jy,kz,n)=3 ! lsp in-cloud555 clouds(ix,jy,kz,n)=3 ! lsp in-cloud 556 556 else 557 557 clouds(ix,jy,kz,n)=2 ! convp in-cloud
Note: See TracChangeset
for help on using the changeset viewer.