Changes in / [6985a98:d1a8707] in flexpart.git
- Location:
- src
- Files:
-
- 1 deleted
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART.f90
rc8fc724 r6896b17 67 67 call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1)) 68 68 69 70 69 ! FLEXPART version string 71 70 flexversion_major = '10' ! Major version number, also used for species file names 72 flexversion='Version '//trim(flexversion_major)//'.1beta (201 6-11-02)'71 flexversion='Version '//trim(flexversion_major)//'.1beta (2017-01-01)' 73 72 verbosity=0 74 73 … … 385 384 end do 386 385 387 ! Inform whether output kernel is used or not388 !*********************************************389 if (lroot) then390 if (lnokernel) then391 write(*,*) "Concentrations are calculated without using kernel"392 else393 write(*,*) "Concentrations are calculated using kernel"394 end if395 end if396 397 398 386 ! Calculate particle trajectories 399 387 !******************************** … … 414 402 415 403 ! NIK 16.02.2005 416 do i=1,nspec 417 write(*,*) '**********************************************' 418 write(*,*) 'Scavenging statistics for species ', species(i), ':' 419 write(*,*) 'Total number of occurences of below-cloud scavenging', & 420 & tot_blc_count(i) 421 write(*,*) 'Total number of occurences of in-cloud scavenging', & 422 & tot_inc_count(i) 423 write(*,*) '**********************************************' 424 end do 425 404 write(*,*) '**********************************************' 405 write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count 406 write(*,*) 'Total number of occurences of in-cloud scavenging', tot_inc_count 407 write(*,*) '**********************************************' 408 426 409 write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE& 427 410 &XPART MODEL RUN!' -
src/FLEXPART_MPI.f90
rc8fc724 r79abee9 55 55 56 56 57 58 57 ! Initialize mpi 58 !********************* 59 59 call mpif_init 60 60 61 61 if (mp_measure_time) call mpif_mtime('flexpart',0) 62 62 63 ! Initialize arrays in com_mod 64 !***************************** 65 66 if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi) 67 63 ! Initialize arrays in com_mod 64 !***************************** 65 call com_mod_allocate_part(maxpart_mpi) 68 66 69 67 ! Generate a large number of random numbers … … 81 79 flexversion_major = '10' ! Major version number, also used for species file names 82 80 ! flexversion='Ver. 10 Beta MPI (2015-05-01)' 83 flexversion='Ver. '//trim(flexversion_major)//' .1beta MPI (2016-11-02)'81 flexversion='Ver. '//trim(flexversion_major)//' Beta MPI (2015-05-01)' 84 82 verbosity=0 85 83 … … 308 306 endif 309 307 310 if (.not.(lmpreader.and.lmp_use_reader)) then 311 do j=1, size(itra1) ! maxpart_mpi 312 itra1(j)=-999999999 313 end do 314 end if 308 do j=1, size(itra1) ! maxpart_mpi 309 itra1(j)=-999999999 310 end do 315 311 316 312 ! For continuation of previous run, read in particle positions … … 322 318 endif 323 319 ! readwind process skips this step 324 if ( .not.(lmpreader.and.lmp_use_reader)) call readpartpositions320 if (lmp_use_reader.and..not.lmpreader) call readpartpositions 325 321 else 326 322 if (verbosity.gt.0 .and. lroot) then … … 429 425 end do 430 426 431 ! Inform whether output kernel is used or not432 !*********************************************433 434 if (lroot) then435 if (lnokernel) then436 write(*,*) "Concentrations are calculated without using kernel"437 else438 write(*,*) "Concentrations are calculated using kernel"439 end if440 end if441 427 442 428 ! Calculate particle trajectories … … 462 448 ! NIK 16.02.2005 463 449 if (lroot) then 464 call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &450 call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, & 465 451 & mp_comm_used, mp_ierr) 466 call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &452 call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, & 467 453 & mp_comm_used, mp_ierr) 468 454 else 469 455 if (mp_partgroup_pid.ge.0) then ! Skip for readwind process 470 call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &456 call MPI_Reduce(tot_blc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, & 471 457 & mp_comm_used, mp_ierr) 472 call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &458 call MPI_Reduce(tot_inc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, & 473 459 & mp_comm_used, mp_ierr) 474 460 end if … … 476 462 477 463 if (lroot) then 478 do i=1,nspec 479 write(*,*) '**********************************************' 480 write(*,*) 'Scavenging statistics for species ', species(i), ':' 481 write(*,*) 'Total number of occurences of below-cloud scavenging', & 482 & tot_blc_count(i) 483 write(*,*) 'Total number of occurences of in-cloud scavenging', & 484 & tot_inc_count(i) 485 write(*,*) '**********************************************' 486 end do 464 write(*,*) '**********************************************' 465 write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count 466 write(*,*) 'Total number of occurences of in-cloud scavenging', tot_inc_count 467 write(*,*) '**********************************************' 487 468 488 469 write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE& -
src/boundcond_domainfill_mpi.f90
r16b61a5 r7999df47 21 21 22 22 subroutine boundcond_domainfill(itime,loutend) 23 ! i i 24 !***************************************************************************** 25 ! * 26 ! Particles are created by this subroutine continuously throughout the * 27 ! simulation at the boundaries of the domain-filling box. * 28 ! All particles carry the same amount of mass which alltogether comprises the* 29 ! mass of air within the box, which remains (more or less) constant. * 30 ! * 31 ! Author: A. Stohl * 32 ! * 33 ! 16 October 2002 * 34 ! * 35 !***************************************************************************** 36 ! * 37 ! Variables: * 38 ! * 39 ! nx_we(2) grid indices for western and eastern boundary of domain- * 40 ! filling trajectory calculations * 41 ! ny_sn(2) grid indices for southern and northern boundary of domain- * 42 ! filling trajectory calculations * 43 ! * 44 !***************************************************************************** 45 ! CHANGES 46 ! 08/2016 eso: MPI version: 47 ! 48 ! -Root process release particles and distributes to other processes. 49 ! Temporary arrays are used, also for the non-root (receiving) processes. 50 ! -The scheme can be improved by having all processes report numpart 51 ! (keeping track of how many particles have left the domain), so that 52 ! a proportional amount of new particles can be distributed (however 53 ! we have a separate function called from timemanager that will 54 ! redistribute particles among processes if there are imbalances) 55 !***************************************************************************** 23 ! i i 24 !***************************************************************************** 25 ! * 26 ! Particles are created by this subroutine continuously throughout the * 27 ! simulation at the boundaries of the domain-filling box. * 28 ! All particles carry the same amount of mass which alltogether comprises the* 29 ! mass of air within the box, which remains (more or less) constant. * 30 ! * 31 ! Author: A. Stohl * 32 ! * 33 ! 16 October 2002 * 34 ! * 35 !***************************************************************************** 36 ! * 37 ! Variables: * 38 ! * 39 ! nx_we(2) grid indices for western and eastern boundary of domain- * 40 ! filling trajectory calculations * 41 ! ny_sn(2) grid indices for southern and northern boundary of domain- * 42 ! filling trajectory calculations * 43 ! * 44 !***************************************************************************** 56 45 57 46 use point_mod … … 66 55 integer :: itime,in,indz,indzp,i,loutend 67 56 integer :: j,k,ix,jy,m,indzh,indexh,minpart,ipart,mmass 68 integer :: numactiveparticles, numpart_total, rel_counter 69 integer,allocatable,dimension(:) :: numrel_mpi !, numactiveparticles_mpi 57 integer :: numactiveparticles 70 58 71 59 real :: windl(2),rhol(2) … … 78 66 79 67 integer :: idummy = -11 80 integer :: mtag81 68 logical :: first_call=.true. 82 ! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if 83 ! needed. 84 real,parameter :: maxpartfract=0.1 85 integer :: tmp_size = int(maxpartfract*maxpart) 86 87 ! Use different seed for each process 69 70 ! Use different seed for each process 88 71 if (first_call) then 89 72 idummy=idummy+mp_seed … … 92 75 93 76 94 ! If domain-filling is global, no boundary conditions are needed95 !***************************************************************77 ! If domain-filling is global, no boundary conditions are needed 78 !*************************************************************** 96 79 97 80 if (gdomainfill) return … … 99 82 accmasst=0. 100 83 numactiveparticles=0 101 ! Keep track of active particles on each process 102 allocate(numrel_mpi(0:mp_partgroup_np-1)) 103 ! numactiveparticles_mpi(0:mp_partgroup_np-1) 104 105 ! New particles to be released on each process 106 numrel_mpi(:)=0 107 108 ! Terminate trajectories that have left the domain, if domain-filling 109 ! trajectory calculation domain is not global. Done for all processes 110 !******************************************************************** 84 85 ! Terminate trajectories that have left the domain, if domain-filling 86 ! trajectory calculation domain is not global 87 !******************************************************************** 111 88 112 89 do i=1,numpart … … 121 98 numactiveparticles+1 122 99 end do 123 ! numactiveparticles_mpi(mp_partid) = numactiveparticles 124 125 126 ! Collect number of active particles from all processes 127 ! call MPI_Allgather(numactiveparticles, 1, MPI_INTEGER, & 128 ! &numactiveparticles_mpi, 1, MPI_INTEGER, mp_comm_used, mp_ierr) 129 130 131 ! Total number of new releases 132 numpart_total = 0 133 134 135 ! This section only done by root process 136 !*************************************** 137 138 if (lroot) then 139 140 ! Use separate arrays for newly released particles 141 !************************************************* 142 143 allocate(itra1_tmp(tmp_size),npoint_tmp(tmp_size),nclass_tmp(tmp_size),& 144 & idt_tmp(tmp_size),itramem_tmp(tmp_size),itrasplit_tmp(tmp_size),& 145 & xtra1_tmp(tmp_size),ytra1_tmp(tmp_size),ztra1_tmp(tmp_size),& 146 & xmass1_tmp(tmp_size, maxspec)) 147 148 ! Initialize all particles as non-existent 149 itra1_tmp(:)=-999999999 150 151 ! Determine auxiliary variables for time interpolation 152 !***************************************************** 153 154 dt1=real(itime-memtime(1)) 155 dt2=real(memtime(2)-itime) 156 dtt=1./(dt1+dt2) 157 158 ! Initialize auxiliary variable used to search for vacant storage space 159 !********************************************************************** 160 161 minpart=1 162 163 !*************************************** 164 ! Western and eastern boundary condition 165 !*************************************** 166 167 ! Loop from south to north 168 !************************* 169 170 do jy=ny_sn(1),ny_sn(2) 171 172 ! Loop over western (index 1) and eastern (index 2) boundary 173 !*********************************************************** 174 175 do k=1,2 176 177 ! Loop over all release locations in a column 178 !******************************************** 179 180 do j=1,numcolumn_we(k,jy) 181 182 ! Determine, for each release location, the area of the corresponding boundary 183 !***************************************************************************** 184 185 if (j.eq.1) then 186 deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2. 187 else if (j.eq.numcolumn_we(k,jy)) then 188 ! deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+ 189 ! + zcolumn_we(k,jy,j))/2. 190 ! In order to avoid taking a very high column for very many particles, 191 ! use the deltaz from one particle below instead 192 deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2. 100 101 102 ! Determine auxiliary variables for time interpolation 103 !***************************************************** 104 105 dt1=real(itime-memtime(1)) 106 dt2=real(memtime(2)-itime) 107 dtt=1./(dt1+dt2) 108 109 ! Initialize auxiliary variable used to search for vacant storage space 110 !********************************************************************** 111 112 minpart=1 113 114 !*************************************** 115 ! Western and eastern boundary condition 116 !*************************************** 117 118 ! Loop from south to north 119 !************************* 120 121 do jy=ny_sn(1),ny_sn(2) 122 123 ! Loop over western (index 1) and eastern (index 2) boundary 124 !*********************************************************** 125 126 do k=1,2 127 128 ! Loop over all release locations in a column 129 !******************************************** 130 131 do j=1,numcolumn_we(k,jy) 132 133 ! Determine, for each release location, the area of the corresponding boundary 134 !***************************************************************************** 135 136 if (j.eq.1) then 137 deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2. 138 else if (j.eq.numcolumn_we(k,jy)) then 139 ! deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+ 140 ! + zcolumn_we(k,jy,j))/2. 141 ! In order to avoid taking a very high column for very many particles, 142 ! use the deltaz from one particle below instead 143 deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2. 144 else 145 deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2. 146 endif 147 if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then 148 boundarea=deltaz*111198.5/2.*dy 149 else 150 boundarea=deltaz*111198.5*dy 151 endif 152 153 154 ! Interpolate the wind velocity and density to the release location 155 !****************************************************************** 156 157 ! Determine the model level below the release position 158 !***************************************************** 159 160 do i=2,nz 161 if (height(i).gt.zcolumn_we(k,jy,j)) then 162 indz=i-1 163 indzp=i 164 goto 6 165 endif 166 end do 167 6 continue 168 169 ! Vertical distance to the level below and above current position 170 !**************************************************************** 171 172 dz1=zcolumn_we(k,jy,j)-height(indz) 173 dz2=height(indzp)-zcolumn_we(k,jy,j) 174 dz=1./(dz1+dz2) 175 176 ! Vertical and temporal interpolation 177 !************************************ 178 179 do m=1,2 180 indexh=memind(m) 181 do in=1,2 182 indzh=indz+in-1 183 windl(in)=uu(nx_we(k),jy,indzh,indexh) 184 rhol(in)=rho(nx_we(k),jy,indzh,indexh) 185 end do 186 187 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 188 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 189 end do 190 191 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 192 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 193 194 ! Calculate mass flux, divided by number of processes 195 !**************************************************** 196 197 fluxofmass=windx*rhox*boundarea*real(lsynctime)/mp_partgroup_np 198 199 200 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 201 ! if it is out of the domain, set accumulated mass flux to zero 202 !****************************************************************************** 203 204 if (k.eq.1) then 205 if (fluxofmass.ge.0.) then 206 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass 193 207 else 194 deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.208 acc_mass_we(k,jy,j)=0. 195 209 endif 196 if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then 197 boundarea=deltaz*111198.5/2.*dy 210 else 211 if (fluxofmass.le.0.) then 212 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass) 198 213 else 199 boundarea=deltaz*111198.5*dy214 acc_mass_we(k,jy,j)=0. 200 215 endif 201 202 203 ! Interpolate the wind velocity and density to the release location 204 !****************************************************************** 205 206 ! Determine the model level below the release position 207 !***************************************************** 208 209 do i=2,nz 210 if (height(i).gt.zcolumn_we(k,jy,j)) then 211 indz=i-1 212 indzp=i 213 goto 6 216 endif 217 accmasst=accmasst+acc_mass_we(k,jy,j) 218 219 ! If the accumulated mass exceeds half the mass that each particle shall carry, 220 ! one (or more) particle(s) is (are) released and the accumulated mass is 221 ! reduced by the mass of this (these) particle(s) 222 !****************************************************************************** 223 224 if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then 225 mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ & 226 xmassperparticle) 227 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- & 228 real(mmass)*xmassperparticle 229 else 230 mmass=0 231 endif 232 233 do m=1,mmass 234 do ipart=minpart,maxpart 235 236 ! If a vacant storage space is found, attribute everything to this array element 237 !***************************************************************************** 238 239 if (itra1(ipart).ne.itime) then 240 241 ! Assign particle positions 242 !************************** 243 244 xtra1(ipart)=real(nx_we(k)) 245 if (jy.eq.ny_sn(1)) then 246 ytra1(ipart)=real(jy)+0.5*ran1(idummy) 247 else if (jy.eq.ny_sn(2)) then 248 ytra1(ipart)=real(jy)-0.5*ran1(idummy) 249 else 250 ytra1(ipart)=real(jy)+(ran1(idummy)-.5) 251 endif 252 if (j.eq.1) then 253 ztra1(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- & 254 zcolumn_we(k,jy,1))/4. 255 else if (j.eq.numcolumn_we(k,jy)) then 256 ztra1(ipart)=(2.*zcolumn_we(k,jy,j)+ & 257 zcolumn_we(k,jy,j-1)+height(nz))/4. 258 else 259 ztra1(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* & 260 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1)) 261 endif 262 263 ! Interpolate PV to the particle position 264 !**************************************** 265 ixm=int(xtra1(ipart)) 266 jym=int(ytra1(ipart)) 267 ixp=ixm+1 268 jyp=jym+1 269 ddx=xtra1(ipart)-real(ixm) 270 ddy=ytra1(ipart)-real(jym) 271 rddx=1.-ddx 272 rddy=1.-ddy 273 p1=rddx*rddy 274 p2=ddx*rddy 275 p3=rddx*ddy 276 p4=ddx*ddy 277 do i=2,nz 278 if (height(i).gt.ztra1(ipart)) then 279 indzm=i-1 280 indzp=i 281 goto 26 282 endif 283 end do 284 26 continue 285 dz1=ztra1(ipart)-height(indzm) 286 dz2=height(indzp)-ztra1(ipart) 287 dz=1./(dz1+dz2) 288 do mm=1,2 289 indexh=memind(mm) 290 do in=1,2 291 indzh=indzm+in-1 292 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 293 +p2*pv(ixp,jym,indzh,indexh) & 294 +p3*pv(ixm,jyp,indzh,indexh) & 295 +p4*pv(ixp,jyp,indzh,indexh) 296 end do 297 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 298 end do 299 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 300 ylat=ylat0+ytra1(ipart)*dy 301 if (ylat.lt.0.) pvpart=-1.*pvpart 302 303 304 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 305 !***************************************************************************** 306 307 if (((ztra1(ipart).gt.3000.).and. & 308 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 309 nclass(ipart)=min(int(ran1(idummy)* & 310 real(nclassunc))+1,nclassunc) 311 numactiveparticles=numactiveparticles+1 312 numparticlecount=numparticlecount+1 313 npoint(ipart)=numparticlecount 314 idt(ipart)=mintime 315 itra1(ipart)=itime 316 itramem(ipart)=itra1(ipart) 317 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit 318 xmass1(ipart,1)=xmassperparticle 319 if (mdomainfill.eq.2) xmass1(ipart,1)= & 320 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9 321 else 322 goto 71 323 endif 324 325 326 ! Increase numpart, if necessary 327 !******************************* 328 329 numpart=max(numpart,ipart) 330 goto 73 ! Storage space has been found, stop searching 214 331 endif 215 332 end do 216 6 continue 217 218 ! Vertical distance to the level below and above current position 219 !**************************************************************** 220 221 dz1=zcolumn_we(k,jy,j)-height(indz) 222 dz2=height(indzp)-zcolumn_we(k,jy,j) 223 dz=1./(dz1+dz2) 224 225 ! Vertical and temporal interpolation 226 !************************************ 227 228 do m=1,2 229 indexh=memind(m) 230 do in=1,2 231 indzh=indz+in-1 232 windl(in)=uu(nx_we(k),jy,indzh,indexh) 233 rhol(in)=rho(nx_we(k),jy,indzh,indexh) 234 end do 235 236 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 237 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 238 end do 239 240 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 241 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 242 243 ! Calculate mass flux 244 !******************** 245 246 fluxofmass=windx*rhox*boundarea*real(lsynctime) 247 248 249 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 250 ! if it is out of the domain, set accumulated mass flux to zero 251 !****************************************************************************** 252 253 if (k.eq.1) then 254 if (fluxofmass.ge.0.) then 255 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass 256 else 257 acc_mass_we(k,jy,j)=0. 258 endif 259 else 260 if (fluxofmass.le.0.) then 261 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass) 262 else 263 acc_mass_we(k,jy,j)=0. 264 endif 265 endif 266 accmasst=accmasst+acc_mass_we(k,jy,j) 267 268 ! If the accumulated mass exceeds half the mass that each particle shall carry, 269 ! one (or more) particle(s) is (are) released and the accumulated mass is 270 ! reduced by the mass of this (these) particle(s) 271 !****************************************************************************** 272 273 if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then 274 mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ & 275 xmassperparticle) 276 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- & 277 real(mmass)*xmassperparticle 278 else 279 mmass=0 280 endif 281 282 do m=1,mmass 283 do ipart=minpart,maxpart 284 285 ! If a vacant storage space is found, attribute everything to this array element 286 ! TODO: for the MPI version this test can be removed, as all 287 ! elements in _tmp arrays are initialized to zero 288 !***************************************************************************** 289 290 if (itra1_tmp(ipart).ne.itime) then 291 292 ! Assign particle positions 293 !************************** 294 295 xtra1_tmp(ipart)=real(nx_we(k)) 296 if (jy.eq.ny_sn(1)) then 297 ytra1_tmp(ipart)=real(jy)+0.5*ran1(idummy) 298 else if (jy.eq.ny_sn(2)) then 299 ytra1_tmp(ipart)=real(jy)-0.5*ran1(idummy) 300 else 301 ytra1_tmp(ipart)=real(jy)+(ran1(idummy)-.5) 302 endif 303 if (j.eq.1) then 304 ztra1_tmp(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- & 305 zcolumn_we(k,jy,1))/4. 306 else if (j.eq.numcolumn_we(k,jy)) then 307 ztra1_tmp(ipart)=(2.*zcolumn_we(k,jy,j)+ & 308 zcolumn_we(k,jy,j-1)+height(nz))/4. 309 else 310 ztra1_tmp(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* & 311 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1)) 312 endif 313 314 ! Interpolate PV to the particle position 315 !**************************************** 316 ixm=int(xtra1_tmp(ipart)) 317 jym=int(ytra1_tmp(ipart)) 318 ixp=ixm+1 319 jyp=jym+1 320 ddx=xtra1_tmp(ipart)-real(ixm) 321 ddy=ytra1_tmp(ipart)-real(jym) 322 rddx=1.-ddx 323 rddy=1.-ddy 324 p1=rddx*rddy 325 p2=ddx*rddy 326 p3=rddx*ddy 327 p4=ddx*ddy 328 do i=2,nz 329 if (height(i).gt.ztra1_tmp(ipart)) then 330 indzm=i-1 331 indzp=i 332 goto 26 333 endif 334 end do 335 26 continue 336 dz1=ztra1_tmp(ipart)-height(indzm) 337 dz2=height(indzp)-ztra1_tmp(ipart) 338 dz=1./(dz1+dz2) 339 do mm=1,2 340 indexh=memind(mm) 341 do in=1,2 342 indzh=indzm+in-1 343 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 344 +p2*pv(ixp,jym,indzh,indexh) & 345 +p3*pv(ixm,jyp,indzh,indexh) & 346 +p4*pv(ixp,jyp,indzh,indexh) 347 end do 348 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 349 end do 350 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 351 ylat=ylat0+ytra1_tmp(ipart)*dy 352 if (ylat.lt.0.) pvpart=-1.*pvpart 353 354 355 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 356 !***************************************************************************** 357 358 if (((ztra1_tmp(ipart).gt.3000.).and. & 359 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 360 nclass_tmp(ipart)=min(int(ran1(idummy)* & 361 real(nclassunc))+1,nclassunc) 362 numactiveparticles=numactiveparticles+1 363 numparticlecount=numparticlecount+1 364 npoint_tmp(ipart)=numparticlecount 365 idt_tmp(ipart)=mintime 366 itra1_tmp(ipart)=itime 367 itramem_tmp(ipart)=itra1_tmp(ipart) 368 itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit 369 xmass1_tmp(ipart,1)=xmassperparticle 370 if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= & 371 xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9 372 else 373 goto 71 374 endif 375 376 377 ! Increase numpart, if necessary 378 !******************************* 379 380 numpart_total=max(numpart_total,ipart) 381 goto 73 ! Storage space has been found, stop searching 382 endif 383 end do 384 if (ipart.gt.tmp_size) & 385 stop 'boundcond_domainfill_mpi.f90: too many particles required' 386 73 minpart=ipart+1 387 71 continue 388 end do 389 390 333 if (ipart.gt.maxpart) & 334 stop 'boundcond_domainfill.f: too many particles required' 335 73 minpart=ipart+1 336 71 continue 391 337 end do 338 339 392 340 end do 393 341 end do 394 395 396 !***************************************** 397 ! Southern and northern boundary condition 398 !***************************************** 399 400 ! Loop from west to east 401 !*********************** 402 403 do ix=nx_we(1),nx_we(2) 404 405 ! Loop over southern (index 1) and northern (index 2) boundary 406 !************************************************************* 407 408 do k=1,2 409 ylat=ylat0+real(ny_sn(k))*dy 410 cosfact=cos(ylat*pi180) 411 412 ! Loop over all release locations in a column 413 !******************************************** 414 415 do j=1,numcolumn_sn(k,ix) 416 417 ! Determine, for each release location, the area of the corresponding boundary 418 !***************************************************************************** 419 420 if (j.eq.1) then 421 deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2. 422 else if (j.eq.numcolumn_sn(k,ix)) then 423 ! deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+ 424 ! + zcolumn_sn(k,ix,j))/2. 425 ! In order to avoid taking a very high column for very many particles, 426 ! use the deltaz from one particle below instead 427 deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2. 342 end do 343 344 345 !***************************************** 346 ! Southern and northern boundary condition 347 !***************************************** 348 349 ! Loop from west to east 350 !*********************** 351 352 do ix=nx_we(1),nx_we(2) 353 354 ! Loop over southern (index 1) and northern (index 2) boundary 355 !************************************************************* 356 357 do k=1,2 358 ylat=ylat0+real(ny_sn(k))*dy 359 cosfact=cos(ylat*pi180) 360 361 ! Loop over all release locations in a column 362 !******************************************** 363 364 do j=1,numcolumn_sn(k,ix) 365 366 ! Determine, for each release location, the area of the corresponding boundary 367 !***************************************************************************** 368 369 if (j.eq.1) then 370 deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2. 371 else if (j.eq.numcolumn_sn(k,ix)) then 372 ! deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+ 373 ! + zcolumn_sn(k,ix,j))/2. 374 ! In order to avoid taking a very high column for very many particles, 375 ! use the deltaz from one particle below instead 376 deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2. 377 else 378 deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2. 379 endif 380 if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then 381 boundarea=deltaz*111198.5/2.*cosfact*dx 382 else 383 boundarea=deltaz*111198.5*cosfact*dx 384 endif 385 386 387 ! Interpolate the wind velocity and density to the release location 388 !****************************************************************** 389 390 ! Determine the model level below the release position 391 !***************************************************** 392 393 do i=2,nz 394 if (height(i).gt.zcolumn_sn(k,ix,j)) then 395 indz=i-1 396 indzp=i 397 goto 16 398 endif 399 end do 400 16 continue 401 402 ! Vertical distance to the level below and above current position 403 !**************************************************************** 404 405 dz1=zcolumn_sn(k,ix,j)-height(indz) 406 dz2=height(indzp)-zcolumn_sn(k,ix,j) 407 dz=1./(dz1+dz2) 408 409 ! Vertical and temporal interpolation 410 !************************************ 411 412 do m=1,2 413 indexh=memind(m) 414 do in=1,2 415 indzh=indz+in-1 416 windl(in)=vv(ix,ny_sn(k),indzh,indexh) 417 rhol(in)=rho(ix,ny_sn(k),indzh,indexh) 418 end do 419 420 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 421 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 422 end do 423 424 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 425 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 426 427 ! Calculate mass flux, divided by number of processes 428 !**************************************************** 429 430 fluxofmass=windx*rhox*boundarea*real(lsynctime)/mp_partgroup_np 431 432 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 433 ! if it is out of the domain, set accumulated mass flux to zero 434 !****************************************************************************** 435 436 if (k.eq.1) then 437 if (fluxofmass.ge.0.) then 438 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass 428 439 else 429 deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.440 acc_mass_sn(k,ix,j)=0. 430 441 endif 431 if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then 432 boundarea=deltaz*111198.5/2.*cosfact*dx 442 else 443 if (fluxofmass.le.0.) then 444 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass) 433 445 else 434 boundarea=deltaz*111198.5*cosfact*dx446 acc_mass_sn(k,ix,j)=0. 435 447 endif 436 437 438 ! Interpolate the wind velocity and density to the release location 439 !****************************************************************** 440 441 ! Determine the model level below the release position 442 !***************************************************** 443 444 do i=2,nz 445 if (height(i).gt.zcolumn_sn(k,ix,j)) then 446 indz=i-1 447 indzp=i 448 goto 16 448 endif 449 accmasst=accmasst+acc_mass_sn(k,ix,j) 450 451 ! If the accumulated mass exceeds half the mass that each particle shall carry, 452 ! one (or more) particle(s) is (are) released and the accumulated mass is 453 ! reduced by the mass of this (these) particle(s) 454 !****************************************************************************** 455 456 if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then 457 mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ & 458 xmassperparticle) 459 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- & 460 real(mmass)*xmassperparticle 461 else 462 mmass=0 463 endif 464 465 do m=1,mmass 466 do ipart=minpart,maxpart 467 468 ! If a vacant storage space is found, attribute everything to this array element 469 !***************************************************************************** 470 471 if (itra1(ipart).ne.itime) then 472 473 ! Assign particle positions 474 !************************** 475 476 ytra1(ipart)=real(ny_sn(k)) 477 if (ix.eq.nx_we(1)) then 478 xtra1(ipart)=real(ix)+0.5*ran1(idummy) 479 else if (ix.eq.nx_we(2)) then 480 xtra1(ipart)=real(ix)-0.5*ran1(idummy) 481 else 482 xtra1(ipart)=real(ix)+(ran1(idummy)-.5) 483 endif 484 if (j.eq.1) then 485 ztra1(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- & 486 zcolumn_sn(k,ix,1))/4. 487 else if (j.eq.numcolumn_sn(k,ix)) then 488 ztra1(ipart)=(2.*zcolumn_sn(k,ix,j)+ & 489 zcolumn_sn(k,ix,j-1)+height(nz))/4. 490 else 491 ztra1(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* & 492 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1)) 493 endif 494 495 496 ! Interpolate PV to the particle position 497 !**************************************** 498 ixm=int(xtra1(ipart)) 499 jym=int(ytra1(ipart)) 500 ixp=ixm+1 501 jyp=jym+1 502 ddx=xtra1(ipart)-real(ixm) 503 ddy=ytra1(ipart)-real(jym) 504 rddx=1.-ddx 505 rddy=1.-ddy 506 p1=rddx*rddy 507 p2=ddx*rddy 508 p3=rddx*ddy 509 p4=ddx*ddy 510 do i=2,nz 511 if (height(i).gt.ztra1(ipart)) then 512 indzm=i-1 513 indzp=i 514 goto 126 515 endif 516 end do 517 126 continue 518 dz1=ztra1(ipart)-height(indzm) 519 dz2=height(indzp)-ztra1(ipart) 520 dz=1./(dz1+dz2) 521 do mm=1,2 522 indexh=memind(mm) 523 do in=1,2 524 indzh=indzm+in-1 525 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 526 +p2*pv(ixp,jym,indzh,indexh) & 527 +p3*pv(ixm,jyp,indzh,indexh) & 528 +p4*pv(ixp,jyp,indzh,indexh) 529 end do 530 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 531 end do 532 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 533 if (ylat.lt.0.) pvpart=-1.*pvpart 534 535 536 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 537 !***************************************************************************** 538 539 if (((ztra1(ipart).gt.3000.).and. & 540 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 541 nclass(ipart)=min(int(ran1(idummy)* & 542 real(nclassunc))+1,nclassunc) 543 numactiveparticles=numactiveparticles+1 544 numparticlecount=numparticlecount+1 545 npoint(ipart)=numparticlecount 546 idt(ipart)=mintime 547 itra1(ipart)=itime 548 itramem(ipart)=itra1(ipart) 549 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit 550 xmass1(ipart,1)=xmassperparticle 551 if (mdomainfill.eq.2) xmass1(ipart,1)= & 552 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9 553 else 554 goto 171 555 endif 556 557 558 ! Increase numpart, if necessary 559 !******************************* 560 numpart=max(numpart,ipart) 561 goto 173 ! Storage space has been found, stop searching 449 562 endif 450 563 end do 451 16 continue 452 453 ! Vertical distance to the level below and above current position 454 !**************************************************************** 455 456 dz1=zcolumn_sn(k,ix,j)-height(indz) 457 dz2=height(indzp)-zcolumn_sn(k,ix,j) 458 dz=1./(dz1+dz2) 459 460 ! Vertical and temporal interpolation 461 !************************************ 462 463 do m=1,2 464 indexh=memind(m) 465 do in=1,2 466 indzh=indz+in-1 467 windl(in)=vv(ix,ny_sn(k),indzh,indexh) 468 rhol(in)=rho(ix,ny_sn(k),indzh,indexh) 469 end do 470 471 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 472 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 473 end do 474 475 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 476 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 477 478 ! Calculate mass flux 479 !******************** 480 481 fluxofmass=windx*rhox*boundarea*real(lsynctime) 482 483 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 484 ! if it is out of the domain, set accumulated mass flux to zero 485 !****************************************************************************** 486 487 if (k.eq.1) then 488 if (fluxofmass.ge.0.) then 489 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass 490 else 491 acc_mass_sn(k,ix,j)=0. 492 endif 493 else 494 if (fluxofmass.le.0.) then 495 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass) 496 else 497 acc_mass_sn(k,ix,j)=0. 498 endif 499 endif 500 accmasst=accmasst+acc_mass_sn(k,ix,j) 501 502 ! If the accumulated mass exceeds half the mass that each particle shall carry, 503 ! one (or more) particle(s) is (are) released and the accumulated mass is 504 ! reduced by the mass of this (these) particle(s) 505 !****************************************************************************** 506 507 if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then 508 mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ & 509 xmassperparticle) 510 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- & 511 real(mmass)*xmassperparticle 512 else 513 mmass=0 514 endif 515 516 do m=1,mmass 517 do ipart=minpart,maxpart 518 519 ! If a vacant storage space is found, attribute everything to this array element 520 !***************************************************************************** 521 522 if (itra1_tmp(ipart).ne.itime) then 523 524 ! Assign particle positions 525 !************************** 526 527 ytra1_tmp(ipart)=real(ny_sn(k)) 528 if (ix.eq.nx_we(1)) then 529 xtra1_tmp(ipart)=real(ix)+0.5*ran1(idummy) 530 else if (ix.eq.nx_we(2)) then 531 xtra1_tmp(ipart)=real(ix)-0.5*ran1(idummy) 532 else 533 xtra1_tmp(ipart)=real(ix)+(ran1(idummy)-.5) 534 endif 535 if (j.eq.1) then 536 ztra1_tmp(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- & 537 zcolumn_sn(k,ix,1))/4. 538 else if (j.eq.numcolumn_sn(k,ix)) then 539 ztra1_tmp(ipart)=(2.*zcolumn_sn(k,ix,j)+ & 540 zcolumn_sn(k,ix,j-1)+height(nz))/4. 541 else 542 ztra1_tmp(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* & 543 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1)) 544 endif 545 546 547 ! Interpolate PV to the particle position 548 !**************************************** 549 ixm=int(xtra1_tmp(ipart)) 550 jym=int(ytra1_tmp(ipart)) 551 ixp=ixm+1 552 jyp=jym+1 553 ddx=xtra1_tmp(ipart)-real(ixm) 554 ddy=ytra1_tmp(ipart)-real(jym) 555 rddx=1.-ddx 556 rddy=1.-ddy 557 p1=rddx*rddy 558 p2=ddx*rddy 559 p3=rddx*ddy 560 p4=ddx*ddy 561 do i=2,nz 562 if (height(i).gt.ztra1_tmp(ipart)) then 563 indzm=i-1 564 indzp=i 565 goto 126 566 endif 567 end do 568 126 continue 569 dz1=ztra1_tmp(ipart)-height(indzm) 570 dz2=height(indzp)-ztra1_tmp(ipart) 571 dz=1./(dz1+dz2) 572 do mm=1,2 573 indexh=memind(mm) 574 do in=1,2 575 indzh=indzm+in-1 576 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 577 +p2*pv(ixp,jym,indzh,indexh) & 578 +p3*pv(ixm,jyp,indzh,indexh) & 579 +p4*pv(ixp,jyp,indzh,indexh) 580 end do 581 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 582 end do 583 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 584 if (ylat.lt.0.) pvpart=-1.*pvpart 585 586 587 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 588 !***************************************************************************** 589 590 if (((ztra1_tmp(ipart).gt.3000.).and. & 591 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 592 nclass_tmp(ipart)=min(int(ran1(idummy)* & 593 real(nclassunc))+1,nclassunc) 594 numactiveparticles=numactiveparticles+1 595 numparticlecount=numparticlecount+1 596 npoint_tmp(ipart)=numparticlecount 597 idt_tmp(ipart)=mintime 598 itra1_tmp(ipart)=itime 599 itramem_tmp(ipart)=itra1_tmp(ipart) 600 itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit 601 xmass1_tmp(ipart,1)=xmassperparticle 602 if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= & 603 xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9 604 else 605 goto 171 606 endif 607 608 609 ! Increase numpart, if necessary 610 !******************************* 611 numpart_total=max(numpart_total,ipart) 612 goto 173 ! Storage space has been found, stop searching 613 endif 614 end do 615 if (ipart.gt.tmp_size) & 616 stop 'boundcond_domainfill.f: too many particles required' 617 173 minpart=ipart+1 618 171 continue 619 end do 620 621 564 if (ipart.gt.maxpart) & 565 stop 'boundcond_domainfill.f: too many particles required' 566 173 minpart=ipart+1 567 171 continue 622 568 end do 569 570 623 571 end do 624 572 end do 625 626 627 ! xm=0.628 ! do i=1,numpart_total629 ! if (itra1_tmp(i).eq.itime) xm=xm+xmass1(i,1)630 ! end do631 632 !write(*,*) itime,numactiveparticles,numparticlecount,numpart,633 ! +xm,accmasst,xm+accmasst634 635 end if ! if lroot636 637 ! Distribute the number of particles to be released638 ! *************************************************639 call MPI_Bcast(numpart_total, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)640 641 do i=0, mp_partgroup_np-1642 numrel_mpi(i) = numpart_total/mp_partgroup_np643 if (i.lt.mod(numpart_total,mp_partgroup_np)) numrel_mpi(i) = numrel_mpi(i) + 1644 573 end do 645 574 646 ! Allocate temporary arrays for receiving processes 647 if (.not.lroot) then 648 allocate(itra1_tmp(numrel_mpi(mp_partid)),& 649 & npoint_tmp(numrel_mpi(mp_partid)),& 650 & nclass_tmp(numrel_mpi(mp_partid)),& 651 & idt_tmp(numrel_mpi(mp_partid)),& 652 & itramem_tmp(numrel_mpi(mp_partid)),& 653 & itrasplit_tmp(numrel_mpi(mp_partid)),& 654 & xtra1_tmp(numrel_mpi(mp_partid)),& 655 & ytra1_tmp(numrel_mpi(mp_partid)),& 656 & ztra1_tmp(numrel_mpi(mp_partid)),& 657 & xmass1_tmp(numrel_mpi(mp_partid),maxspec)) 658 659 ! Initialize all particles as non-existent 660 itra1_tmp(:)=-999999999 661 end if 662 663 ! Distribute particles 664 ! Keep track of released particles so far 665 rel_counter = 0 666 mtag = 1000 667 668 do i=0, mp_partgroup_np-1 669 670 ! For root process, nothing to do except update release count 671 if (i.eq.0) then 672 rel_counter = rel_counter + numrel_mpi(i) 673 cycle 674 end if 675 676 ! Send particles from root to non-root processes 677 if (lroot.and.numrel_mpi(i).gt.0) then 678 679 call MPI_SEND(nclass_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 680 &numrel_mpi(i),MPI_INTEGER,i,mtag+1*i,mp_comm_used,mp_ierr) 681 682 call MPI_SEND(npoint_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 683 &numrel_mpi(i),MPI_INTEGER,i,mtag+2*i,mp_comm_used,mp_ierr) 684 685 call MPI_SEND(itra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 686 &numrel_mpi(i),MPI_INTEGER,i,mtag+3*i,mp_comm_used,mp_ierr) 687 688 call MPI_SEND(idt_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 689 &numrel_mpi(i),MPI_INTEGER,i,mtag+4*i,mp_comm_used,mp_ierr) 690 691 call MPI_SEND(itramem_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 692 &numrel_mpi(i),MPI_INTEGER,i,mtag+5*i,mp_comm_used,mp_ierr) 693 694 call MPI_SEND(itrasplit_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 695 &numrel_mpi(i),MPI_INTEGER,i,mtag+6*i,mp_comm_used,mp_ierr) 696 697 call MPI_SEND(xtra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 698 &numrel_mpi(i),mp_dp,i,mtag+7*i,mp_comm_used,mp_ierr) 699 700 call MPI_SEND(ytra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 701 &numrel_mpi(i),mp_dp,i,mtag+8*i,mp_comm_used,mp_ierr) 702 703 call MPI_SEND(ztra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 704 &numrel_mpi(i),mp_sp,i,mtag+9*i,mp_comm_used,mp_ierr) 705 706 do j=1,nspec 707 call MPI_SEND(xmass1_tmp(rel_counter+1:rel_counter+numrel_mpi(i),j),& 708 &numrel_mpi(i),mp_sp,i,mtag+(9+j)*i,mp_comm_used,mp_ierr) 709 end do 710 711 ! Non-root processes issue receive requests 712 else if (i.eq.mp_partid.and.numrel_mpi(i).gt.0) then 713 call MPI_RECV(nclass_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 714 &MPI_INTEGER,id_root,mtag+1*i,mp_comm_used,mp_status,mp_ierr) 715 716 call MPI_RECV(npoint_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 717 &MPI_INTEGER,id_root,mtag+2*i,mp_comm_used,mp_status,mp_ierr) 718 719 call MPI_RECV(itra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 720 &MPI_INTEGER,id_root,mtag+3*i,mp_comm_used,mp_status,mp_ierr) 721 722 call MPI_RECV(idt_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 723 &MPI_INTEGER,id_root,mtag+4*i,mp_comm_used,mp_status,mp_ierr) 724 725 call MPI_RECV(itramem_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 726 &MPI_INTEGER,id_root,mtag+5*i,mp_comm_used,mp_status,mp_ierr) 727 728 call MPI_RECV(itrasplit_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 729 &MPI_INTEGER,id_root,mtag+6*i,mp_comm_used,mp_status,mp_ierr) 730 731 call MPI_RECV(xtra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 732 &mp_dp,id_root,mtag+7*i,mp_comm_used,mp_status,mp_ierr) 733 734 call MPI_RECV(ytra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 735 &mp_dp,id_root,mtag+8*i,mp_comm_used,mp_status,mp_ierr) 736 737 call MPI_RECV(ztra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 738 &mp_sp,id_root,mtag+9*i,mp_comm_used,mp_status,mp_ierr) 739 740 do j=1,nspec 741 call MPI_RECV(xmass1_tmp(1:numrel_mpi(i),j),numrel_mpi(i),& 742 &mp_sp,id_root,mtag+(9+j)*i,mp_comm_used,mp_status,mp_ierr) 743 744 end do 745 end if 746 rel_counter = rel_counter + numrel_mpi(i) 575 576 xm=0. 577 do i=1,numpart 578 if (itra1(i).eq.itime) xm=xm+xmass1(i,1) 747 579 end do 748 580 749 ! Find free storage space for the new particles. 750 ! This section is independent of the redistribution scheme used 751 ! ******************************************************************** 752 753 ! Keep track of released particles so far 754 minpart=1 755 756 ! The algorithm should be correct also for root process 757 do i=1, numrel_mpi(mp_partid) 758 do ipart=minpart, maxpart 759 if (itra1(ipart).ne.itime) then 760 itra1(ipart) = itra1_tmp(i) 761 npoint(ipart) = npoint_tmp(i) 762 nclass(ipart) = nclass_tmp(i) 763 idt(ipart) = idt_tmp(i) 764 itramem(ipart) = itramem_tmp(i) 765 itrasplit(ipart) = itrasplit_tmp(i) 766 xtra1(ipart) = xtra1_tmp(i) 767 ytra1(ipart) = ytra1_tmp(i) 768 ztra1(ipart) = ztra1_tmp(i) 769 xmass1(ipart,:) = xmass1_tmp(i,:) 770 ! Increase numpart, if necessary 771 numpart=max(numpart,ipart) 772 goto 200 ! Storage space has been found, stop searching 773 end if 774 end do 775 200 minpart=ipart+1 776 end do 777 778 ! If particles shall be dumped, then accumulated masses at the domain boundaries 779 ! must be dumped, too, to be used for later runs 780 !***************************************************************************** 781 581 !write(*,*) itime,numactiveparticles,numparticlecount,numpart, 582 ! +xm,accmasst,xm+accmasst 583 584 585 ! If particles shall be dumped, then accumulated masses at the domain boundaries 586 ! must be dumped, too, to be used for later runs 587 !***************************************************************************** 588 589 ! :TODO: eso parallelize 782 590 if ((ipout.gt.0).and.(itime.eq.loutend)) then 783 591 if (lroot) then 784 call mpif_mtime('iotime',0)785 592 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 786 593 form='unformatted') … … 788 595 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn 789 596 close(unitboundcond) 790 call mpif_mtime('iotime',1)791 597 end if 792 598 endif 793 599 794 ! Deallocate temporary arrays795 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,&796 & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp,numrel_mpi)797 ! numactiveparticles_mpi798 799 800 600 end subroutine boundcond_domainfill -
src/com_mod.f90
r5888298 r5888298 132 132 133 133 !ZHG SEP 2015 wheather or not to read clouds from GRIB 134 logical :: readclouds =.false.134 logical :: readclouds 135 135 !ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed) 136 logical :: sumclouds =.false.136 logical :: sumclouds 137 137 138 138 logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest … … 140 140 141 141 !NIK 16.02.2015 142 integer(selected_int_kind(16)), dimension(maxspec) :: tot_blc_count=0, & 143 &tot_inc_count=0 142 integer(selected_int_kind(16)) :: tot_blc_count=0, tot_inc_count=0 144 143 145 144 … … 577 576 real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn 578 577 !real outheight(maxzgrid),outheighthalf(maxzgrid) 579 580 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),& 581 & OHREA,ASSSPEC 578 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC 582 579 logical :: DRYBKDEP,WETBKDEP 583 580 … … 597 594 ! DRYDEPSPEC .true., if dry deposition is switched on for that species 598 595 ! WETDEP .true., if wet deposition is switched on 599 ! WETDEPSPEC .true., if wet deposition is switched on for that species600 596 ! OHREA .true., if OH reaction is switched on 601 597 ! ASSSPEC .true., if there are two species asscoiated … … 751 747 752 748 ! These variables are used to avoid having separate versions of 753 ! files in cases where differences with MPI version areminor (eso)749 ! files in cases where differences with MPI version is minor (eso) 754 750 !***************************************************************** 755 751 integer :: mpi_mode=0 ! .gt. 0 if running MPI version -
src/conccalc.f90
r9287c01 r9287c01 127 127 !***************************************************************************** 128 128 do ind=indz,indzp 129 rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind, memind(2)) &129 rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) & 130 130 +p2*rho(ixp,jy ,ind,2) & 131 131 +p3*rho(ix ,jyp,ind,2) & … … 181 181 !***************************************************************************** 182 182 183 if ( lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &183 if (((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & 184 184 (xl.gt.real(numxgrid-1)-0.5).or. & 185 (yl.gt.real(numygrid-1)-0.5)) ) then ! no kernel, direct attribution to grid cell185 (yl.gt.real(numygrid-1)-0.5)).or.(.not.usekernel)) then ! no kernel, direct attribution to grid cell 186 186 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 187 187 (jy.le.numygrid-1)) then -
src/conccalc_mpi.f90
r4c64400 r5f9d14a 194 194 !***************************************************************************** 195 195 196 if ( lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &196 if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & 197 197 (xl.gt.real(numxgrid-1)-0.5).or. & 198 198 (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell -
src/concoutput.f90
r16b61a5 r54cbd6c 635 635 !************************* 636 636 637 ! do ks=1,nspec638 ! do kp=1,maxpointspec_act639 ! do i=1,numreceptor640 ! creceptor(i,ks)=0.641 ! end do642 ! do jy=0,numygrid-1643 ! do ix=0,numxgrid-1644 ! do l=1,nclassunc645 ! do nage=1,nageclass646 ! do kz=1,numzgrid647 ! gridunc(ix,jy,kz,ks,kp,l,nage)=0.648 ! end do649 ! end do650 ! end do651 ! end do652 ! end do653 ! end do654 ! end do637 ! do ks=1,nspec 638 ! do kp=1,maxpointspec_act 639 ! do i=1,numreceptor 640 ! creceptor(i,ks)=0. 641 ! end do 642 ! do jy=0,numygrid-1 643 ! do ix=0,numxgrid-1 644 ! do l=1,nclassunc 645 ! do nage=1,nageclass 646 ! do kz=1,numzgrid 647 ! gridunc(ix,jy,kz,ks,kp,l,nage)=0. 648 ! end do 649 ! end do 650 ! end do 651 ! end do 652 ! end do 653 ! end do 654 ! end do 655 655 creceptor(:,:)=0. 656 656 gridunc(:,:,:,:,:,:,:)=0. -
src/concoutput_nest_mpi.f90
r16b61a5 r38b7917 105 105 ! Measure execution time 106 106 if (mp_measure_time) call mpif_mtime('iotime',0) 107 ! call cpu_time(mp_root_time_beg) 108 ! mp_root_wtime_beg = mpi_wtime() 109 ! end if 107 110 108 111 if (verbosity.eq.1) then -
src/concoutput_surf.f90
r16b61a5 r6a678e3 21 21 22 22 subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, & 23 drygridtotalunc)24 ! i i o o25 ! o26 !*****************************************************************************27 ! *28 ! Output of the concentration grid and the receptor concentrations. *29 ! *30 ! Author: A. Stohl *31 ! *32 ! 24 May 1995 *33 ! *34 ! 13 April 1999, Major update: if output size is smaller, dump output *35 ! in sparse matrix format; additional output of *36 ! uncertainty *37 ! *38 ! 05 April 2000, Major update: output of age classes; output for backward*39 ! runs is time spent in grid cell times total mass of *40 ! species. *41 ! *42 ! 17 February 2002, Appropriate dimensions for backward and forward runs *43 ! are now specified in file par_mod *44 ! *45 ! June 2006, write grid in sparse matrix with a single write command *46 ! in order to save disk space *47 ! *48 ! 2008 new sparse matrix format *49 ! *50 !*****************************************************************************51 ! *52 ! Variables: *53 ! outnum number of samples *54 ! ncells number of cells with non-zero concentrations *55 ! sparse .true. if in sparse matrix format, else .false. *56 ! tot_mu 1 for forward, initial mass mixing ration for backw. runs *57 ! *58 !*****************************************************************************23 drygridtotalunc) 24 ! i i o o 25 ! o 26 !***************************************************************************** 27 ! * 28 ! Output of the concentration grid and the receptor concentrations. * 29 ! * 30 ! Author: A. Stohl * 31 ! * 32 ! 24 May 1995 * 33 ! * 34 ! 13 April 1999, Major update: if output size is smaller, dump output * 35 ! in sparse matrix format; additional output of * 36 ! uncertainty * 37 ! * 38 ! 05 April 2000, Major update: output of age classes; output for backward* 39 ! runs is time spent in grid cell times total mass of * 40 ! species. * 41 ! * 42 ! 17 February 2002, Appropriate dimensions for backward and forward runs * 43 ! are now specified in file par_mod * 44 ! * 45 ! June 2006, write grid in sparse matrix with a single write command * 46 ! in order to save disk space * 47 ! * 48 ! 2008 new sparse matrix format * 49 ! * 50 !***************************************************************************** 51 ! * 52 ! Variables: * 53 ! outnum number of samples * 54 ! ncells number of cells with non-zero concentrations * 55 ! sparse .true. if in sparse matrix format, else .false. * 56 ! tot_mu 1 for forward, initial mass mixing ration for backw. runs * 57 ! * 58 !***************************************************************************** 59 59 60 60 use unc_mod … … 73 73 real :: outnum,densityoutrecept(maxreceptor),xl,yl 74 74 75 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),76 ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,77 ! + maxageclass)78 !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,79 ! + maxageclass)80 !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,81 ! + maxpointspec_act,maxageclass)82 !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,83 ! + maxpointspec_act,maxageclass),84 ! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,85 ! + maxpointspec_act,maxageclass),86 ! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,87 ! + maxpointspec_act,maxageclass)88 !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)89 !real sparse_dump_r(numxgrid*numygrid*numzgrid)90 !integer sparse_dump_i(numxgrid*numygrid*numzgrid)91 92 !real sparse_dump_u(numxgrid*numygrid*numzgrid)75 !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid), 76 ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act, 77 ! + maxageclass) 78 !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act, 79 ! + maxageclass) 80 !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec, 81 ! + maxpointspec_act,maxageclass) 82 !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, 83 ! + maxpointspec_act,maxageclass), 84 ! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec, 85 ! + maxpointspec_act,maxageclass), 86 ! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec, 87 ! + maxpointspec_act,maxageclass) 88 !real factor(0:numxgrid-1,0:numygrid-1,numzgrid) 89 !real sparse_dump_r(numxgrid*numygrid*numzgrid) 90 !integer sparse_dump_i(numxgrid*numygrid*numzgrid) 91 92 !real sparse_dump_u(numxgrid*numygrid*numzgrid) 93 93 real(dep_prec) :: auxgrid(nclassunc) 94 94 real(sp) :: gridtotal,gridsigmatotal,gridtotalunc … … 104 104 105 105 if (verbosity.eq.1) then 106 print*,'inside concoutput_surf '107 CALL SYSTEM_CLOCK(count_clock)108 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0109 endif 110 111 ! Determine current calendar date, needed for the file name112 !**********************************************************106 print*,'inside concoutput_surf ' 107 CALL SYSTEM_CLOCK(count_clock) 108 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 109 endif 110 111 ! Determine current calendar date, needed for the file name 112 !********************************************************** 113 113 114 114 jul=bdate+real(itime,kind=dp)/86400._dp … … 116 116 write(adate,'(i8.8)') jjjjmmdd 117 117 write(atime,'(i6.6)') ihmmss 118 !write(unitdates,'(a)') adate//atime119 120 open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')121 write(unitdates,'(a)') adate//atime122 close(unitdates)123 124 ! For forward simulations, output fields have dimension MAXSPEC,125 ! for backward simulations, output fields have dimension MAXPOINT.126 ! Thus, make loops either about nspec, or about numpoint127 !*****************************************************************118 !write(unitdates,'(a)') adate//atime 119 120 open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND') 121 write(unitdates,'(a)') adate//atime 122 close(unitdates) 123 124 ! For forward simulations, output fields have dimension MAXSPEC, 125 ! for backward simulations, output fields have dimension MAXPOINT. 126 ! Thus, make loops either about nspec, or about numpoint 127 !***************************************************************** 128 128 129 129 … … 144 144 145 145 if (verbosity.eq.1) then 146 print*,'concoutput_surf 2'147 CALL SYSTEM_CLOCK(count_clock)148 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0149 endif 150 151 !*******************************************************************152 ! Compute air density: sufficiently accurate to take it153 ! from coarse grid at some time154 ! Determine center altitude of output layer, and interpolate density155 ! data to that altitude156 !*******************************************************************146 print*,'concoutput_surf 2' 147 CALL SYSTEM_CLOCK(count_clock) 148 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 149 endif 150 151 !******************************************************************* 152 ! Compute air density: sufficiently accurate to take it 153 ! from coarse grid at some time 154 ! Determine center altitude of output layer, and interpolate density 155 ! data to that altitude 156 !******************************************************************* 157 157 158 158 do kz=1,numzgrid … … 166 166 (height(kzz).gt.halfheight)) goto 46 167 167 end do 168 46 kzz=max(min(kzz,nz),2)168 46 kzz=max(min(kzz,nz),2) 169 169 dz1=halfheight-height(kzz-1) 170 170 dz2=height(kzz)-halfheight … … 184 184 end do 185 185 186 do i=1,numreceptor187 xl=xreceptor(i)188 yl=yreceptor(i)189 iix=max(min(nint(xl),nxmin1),0)190 jjy=max(min(nint(yl),nymin1),0)191 densityoutrecept(i)=rho(iix,jjy,1,2)192 end do193 194 195 ! Output is different for forward and backward simulations196 do kz=1,numzgrid197 do jy=0,numygrid-1198 do ix=0,numxgrid-1199 if (ldirect.eq.1) then200 factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum201 else202 factor3d(ix,jy,kz)=real(abs(loutaver))/outnum203 endif204 end do205 end do206 end do207 208 !*********************************************************************209 ! Determine the standard deviation of the mean concentration or mixing210 ! ratio (uncertainty of the output) and the dry and wet deposition211 !*********************************************************************186 do i=1,numreceptor 187 xl=xreceptor(i) 188 yl=yreceptor(i) 189 iix=max(min(nint(xl),nxmin1),0) 190 jjy=max(min(nint(yl),nymin1),0) 191 densityoutrecept(i)=rho(iix,jjy,1,2) 192 end do 193 194 195 ! Output is different for forward and backward simulations 196 do kz=1,numzgrid 197 do jy=0,numygrid-1 198 do ix=0,numxgrid-1 199 if (ldirect.eq.1) then 200 factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum 201 else 202 factor3d(ix,jy,kz)=real(abs(loutaver))/outnum 203 endif 204 end do 205 end do 206 end do 207 208 !********************************************************************* 209 ! Determine the standard deviation of the mean concentration or mixing 210 ! ratio (uncertainty of the output) and the dry and wet deposition 211 !********************************************************************* 212 212 213 213 if (verbosity.eq.1) then 214 print*,'concoutput_surf 3 (sd)'215 CALL SYSTEM_CLOCK(count_clock)216 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0214 print*,'concoutput_surf 3 (sd)' 215 CALL SYSTEM_CLOCK(count_clock) 216 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 217 217 endif 218 218 gridtotal=0. … … 228 228 do ks=1,nspec 229 229 230 write(anspec,'(i3.3)') ks 231 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 232 if (ldirect.eq.1) then 233 open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// & 234 atime//'_'//anspec,form='unformatted') 235 else 236 open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// & 237 atime//'_'//anspec,form='unformatted') 238 endif 239 write(unitoutgrid) itime 230 write(anspec,'(i3.3)') ks 231 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 232 if (ldirect.eq.1) then 233 open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// & 234 atime//'_'//anspec,form='unformatted') 235 else 236 open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// & 237 atime//'_'//anspec,form='unformatted') 240 238 endif 241 242 if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio243 open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// & 244 atime//'_'//anspec,form='unformatted')245 246 write(unitoutgridppt) itime247 endif 248 249 do kp=1,maxpointspec_act250 do nage=1,nageclass 251 252 do jy=0,numygrid-1253 do ix=0,numxgrid-1 254 255 ! WET DEPOSITION 256 if ((WETDEP).and.(ldirect.gt.0)) then 257 do l=1,nclassunc258 auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)259 end do260 call mean(auxgrid,wetgrid(ix,jy), &261 wetgridsigma(ix,jy),nclassunc)262 ! Multiply by number of classes to get total concentration 263 wetgrid(ix,jy)=wetgrid(ix,jy) &264 *nclassunc265 wetgridtotal=wetgridtotal+wetgrid(ix,jy)266 ! Calculate standard deviation of the mean 267 wetgridsigma(ix,jy)= &268 wetgridsigma(ix,jy)* &269 sqrt(real(nclassunc))270 wetgridsigmatotal=wetgridsigmatotal+&271 wetgridsigma(ix,jy)272 endif273 274 ! DRY DEPOSITION 275 if ((DRYDEP).and.(ldirect.gt.0)) then 276 do l=1,nclassunc277 auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)278 end do279 call mean(auxgrid,drygrid(ix,jy), &280 drygridsigma(ix,jy),nclassunc)281 ! Multiply by number of classes to get total concentration 282 drygrid(ix,jy)=drygrid(ix,jy)* &283 nclassunc284 drygridtotal=drygridtotal+drygrid(ix,jy)285 ! Calculate standard deviation of the mean 286 drygridsigma(ix,jy)= &287 drygridsigma(ix,jy)* &288 sqrt(real(nclassunc))289 125 drygridsigmatotal=drygridsigmatotal+&290 drygridsigma(ix,jy)291 endif 292 293 ! CONCENTRATION OR MIXING RATIO 294 do kz=1,numzgrid 295 do l=1,nclassunc296 auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)297 end do298 call mean(auxgrid,grid(ix,jy,kz), &299 gridsigma(ix,jy,kz),nclassunc)300 ! Multiply by number of classes to get total concentration 301 grid(ix,jy,kz)= &302 grid(ix,jy,kz)*nclassunc303 gridtotal=gridtotal+grid(ix,jy,kz)304 ! Calculate standard deviation of the mean 305 gridsigma(ix,jy,kz)= &306 gridsigma(ix,jy,kz)* &307 sqrt(real(nclassunc))308 gridsigmatotal=gridsigmatotal+&309 gridsigma(ix,jy,kz)310 end do311 end do239 write(unitoutgrid) itime 240 endif 241 242 if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio 243 open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// & 244 atime//'_'//anspec,form='unformatted') 245 246 write(unitoutgridppt) itime 247 endif 248 249 do kp=1,maxpointspec_act 250 do nage=1,nageclass 251 252 do jy=0,numygrid-1 253 do ix=0,numxgrid-1 254 255 ! WET DEPOSITION 256 if ((WETDEP).and.(ldirect.gt.0)) then 257 do l=1,nclassunc 258 auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage) 259 end do 260 call mean(auxgrid,wetgrid(ix,jy), & 261 wetgridsigma(ix,jy),nclassunc) 262 ! Multiply by number of classes to get total concentration 263 wetgrid(ix,jy)=wetgrid(ix,jy) & 264 *nclassunc 265 wetgridtotal=wetgridtotal+wetgrid(ix,jy) 266 ! Calculate standard deviation of the mean 267 wetgridsigma(ix,jy)= & 268 wetgridsigma(ix,jy)* & 269 sqrt(real(nclassunc)) 270 wetgridsigmatotal=wetgridsigmatotal+ & 271 wetgridsigma(ix,jy) 272 endif 273 274 ! DRY DEPOSITION 275 if ((DRYDEP).and.(ldirect.gt.0)) then 276 do l=1,nclassunc 277 auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage) 278 end do 279 call mean(auxgrid,drygrid(ix,jy), & 280 drygridsigma(ix,jy),nclassunc) 281 ! Multiply by number of classes to get total concentration 282 drygrid(ix,jy)=drygrid(ix,jy)* & 283 nclassunc 284 drygridtotal=drygridtotal+drygrid(ix,jy) 285 ! Calculate standard deviation of the mean 286 drygridsigma(ix,jy)= & 287 drygridsigma(ix,jy)* & 288 sqrt(real(nclassunc)) 289 125 drygridsigmatotal=drygridsigmatotal+ & 290 drygridsigma(ix,jy) 291 endif 292 293 ! CONCENTRATION OR MIXING RATIO 294 do kz=1,numzgrid 295 do l=1,nclassunc 296 auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage) 297 end do 298 call mean(auxgrid,grid(ix,jy,kz), & 299 gridsigma(ix,jy,kz),nclassunc) 300 ! Multiply by number of classes to get total concentration 301 grid(ix,jy,kz)= & 302 grid(ix,jy,kz)*nclassunc 303 gridtotal=gridtotal+grid(ix,jy,kz) 304 ! Calculate standard deviation of the mean 305 gridsigma(ix,jy,kz)= & 306 gridsigma(ix,jy,kz)* & 307 sqrt(real(nclassunc)) 308 gridsigmatotal=gridsigmatotal+ & 309 gridsigma(ix,jy,kz) 312 310 end do 313 314 315 !******************************************************************* 316 ! Generate output: may be in concentration (ng/m3) or in mixing 317 ! ratio (ppt) or both 318 ! Output the position and the values alternated multiplied by 319 ! 1 or -1, first line is number of values, number of positions 320 ! For backward simulations, the unit is seconds, stored in grid_time 321 !******************************************************************* 322 323 if (verbosity.eq.1) then 324 print*,'concoutput_surf 4 (output)' 325 CALL SYSTEM_CLOCK(count_clock) 326 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 327 endif 328 329 ! Concentration output 330 !********************* 331 332 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 333 334 if (verbosity.eq.1) then 335 print*,'concoutput_surf (Wet deposition)' 336 CALL SYSTEM_CLOCK(count_clock) 337 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 338 endif 339 340 ! Wet deposition 341 sp_count_i=0 342 sp_count_r=0 343 sp_fact=-1. 344 sp_zer=.true. 345 if ((ldirect.eq.1).and.(WETDEP)) then 346 do jy=0,numygrid-1 347 do ix=0,numxgrid-1 348 ! concentraion greater zero 349 if (wetgrid(ix,jy).gt.smallnum) then 350 if (sp_zer.eqv..true.) then ! first non zero value 311 end do 312 end do 313 314 315 !******************************************************************* 316 ! Generate output: may be in concentration (ng/m3) or in mixing 317 ! ratio (ppt) or both 318 ! Output the position and the values alternated multiplied by 319 ! 1 or -1, first line is number of values, number of positions 320 ! For backward simulations, the unit is seconds, stored in grid_time 321 !******************************************************************* 322 323 if (verbosity.eq.1) then 324 print*,'concoutput_surf 4 (output)' 325 CALL SYSTEM_CLOCK(count_clock) 326 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 327 endif 328 329 ! Concentration output 330 !********************* 331 332 if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then 333 334 if (verbosity.eq.1) then 335 print*,'concoutput_surf (Wet deposition)' 336 CALL SYSTEM_CLOCK(count_clock) 337 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 338 endif 339 340 ! Wet deposition 341 sp_count_i=0 342 sp_count_r=0 343 sp_fact=-1. 344 sp_zer=.true. 345 if ((ldirect.eq.1).and.(WETDEP)) then 346 do jy=0,numygrid-1 347 do ix=0,numxgrid-1 348 ! concentraion greater zero 349 if (wetgrid(ix,jy).gt.smallnum) then 350 if (sp_zer.eqv..true.) then ! first non zero value 351 351 sp_count_i=sp_count_i+1 352 352 sparse_dump_i(sp_count_i)=ix+jy*numxgrid 353 353 sp_zer=.false. 354 354 sp_fact=sp_fact*(-1.) 355 356 357 358 359 360 361 355 endif 356 sp_count_r=sp_count_r+1 357 sparse_dump_r(sp_count_r)= & 358 sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy) 359 sparse_dump_u(sp_count_r)= & 360 1.e12*wetgridsigma(ix,jy)/area(ix,jy) 361 else ! concentration is zero 362 362 sp_zer=.true. 363 364 365 366 363 endif 364 end do 365 end do 366 else 367 367 sp_count_i=0 368 368 sp_count_r=0 369 370 371 372 373 369 endif 370 write(unitoutgrid) sp_count_i 371 write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i) 372 write(unitoutgrid) sp_count_r 373 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 374 374 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 375 375 376 377 378 379 380 381 ! Dry deposition382 383 384 385 386 387 388 389 390 376 if (verbosity.eq.1) then 377 print*,'concoutput_surf (Dry deposition)' 378 CALL SYSTEM_CLOCK(count_clock) 379 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 380 endif 381 ! Dry deposition 382 sp_count_i=0 383 sp_count_r=0 384 sp_fact=-1. 385 sp_zer=.true. 386 if ((ldirect.eq.1).and.(DRYDEP)) then 387 do jy=0,numygrid-1 388 do ix=0,numxgrid-1 389 if (drygrid(ix,jy).gt.smallnum) then 390 if (sp_zer.eqv..true.) then ! first non zero value 391 391 sp_count_i=sp_count_i+1 392 392 sparse_dump_i(sp_count_i)=ix+jy*numxgrid 393 393 sp_zer=.false. 394 394 sp_fact=sp_fact*(-1.) 395 396 397 398 399 395 endif 396 sp_count_r=sp_count_r+1 397 sparse_dump_r(sp_count_r)= & 398 sp_fact* & 399 1.e12*drygrid(ix,jy)/area(ix,jy) 400 400 sparse_dump_u(sp_count_r)= & 401 402 401 1.e12*drygridsigma(ix,jy)/area(ix,jy) 402 else ! concentration is zero 403 403 sp_zer=.true. 404 405 406 407 404 endif 405 end do 406 end do 407 else 408 408 sp_count_i=0 409 409 sp_count_r=0 410 411 412 413 414 410 endif 411 write(unitoutgrid) sp_count_i 412 write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i) 413 write(unitoutgrid) sp_count_r 414 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 415 415 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 416 416 417 418 419 420 421 422 423 ! Concentrations424 425 ! surf_only write only 1st layer426 427 428 429 430 417 if (verbosity.eq.1) then 418 print*,'concoutput_surf (Concentrations)' 419 CALL SYSTEM_CLOCK(count_clock) 420 WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0 421 endif 422 423 ! Concentrations 424 425 ! surf_only write only 1st layer 426 427 sp_count_i=0 428 sp_count_r=0 429 sp_fact=-1. 430 sp_zer=.true. 431 431 do kz=1,1 432 432 do jy=0,numygrid-1 … … 440 440 sp_fact=sp_fact*(-1.) 441 441 endif 442 sp_count_r=sp_count_r+1443 sparse_dump_r(sp_count_r)= &444 sp_fact* &445 grid(ix,jy,kz)* &446 factor3d(ix,jy,kz)/tot_mu(ks,kp)447 ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)448 ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp449 sparse_dump_u(sp_count_r)= &450 gridsigma(ix,jy,kz)* &451 factor3d(ix,jy,kz)/tot_mu(ks,kp)442 sp_count_r=sp_count_r+1 443 sparse_dump_r(sp_count_r)= & 444 sp_fact* & 445 grid(ix,jy,kz)* & 446 factor3d(ix,jy,kz)/tot_mu(ks,kp) 447 ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0) 448 ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp 449 sparse_dump_u(sp_count_r)= & 450 gridsigma(ix,jy,kz)* & 451 factor3d(ix,jy,kz)/tot_mu(ks,kp) 452 452 else ! concentration is zero 453 453 sp_zer=.true. 454 454 endif 455 456 457 458 459 460 461 455 end do 456 end do 457 end do 458 write(unitoutgrid) sp_count_i 459 write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i) 460 write(unitoutgrid) sp_count_r 461 write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r) 462 462 ! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r) 463 463 464 465 466 ! Mixing ratio output467 !********************468 469 470 471 ! Wet deposition472 473 474 475 476 477 478 464 endif ! concentration output 465 466 ! Mixing ratio output 467 !******************** 468 469 if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio 470 471 ! Wet deposition 472 sp_count_i=0 473 sp_count_r=0 474 sp_fact=-1. 475 sp_zer=.true. 476 if ((ldirect.eq.1).and.(WETDEP)) then 477 do jy=0,numygrid-1 478 do ix=0,numxgrid-1 479 479 if (wetgrid(ix,jy).gt.smallnum) then 480 480 if (sp_zer.eqv..true.) then ! first non zero value … … 484 484 sp_zer=.false. 485 485 sp_fact=sp_fact*(-1.) 486 487 488 489 490 491 492 493 486 endif 487 sp_count_r=sp_count_r+1 488 sparse_dump_r(sp_count_r)= & 489 sp_fact* & 490 1.e12*wetgrid(ix,jy)/area(ix,jy) 491 sparse_dump_u(sp_count_r)= & 492 1.e12*wetgridsigma(ix,jy)/area(ix,jy) 493 else ! concentration is zero 494 494 sp_zer=.true. 495 496 497 498 499 500 501 502 503 504 505 495 endif 496 end do 497 end do 498 else 499 sp_count_i=0 500 sp_count_r=0 501 endif 502 write(unitoutgridppt) sp_count_i 503 write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i) 504 write(unitoutgridppt) sp_count_r 505 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 506 506 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 507 507 508 508 509 ! Dry deposition510 511 512 513 514 515 516 509 ! Dry deposition 510 sp_count_i=0 511 sp_count_r=0 512 sp_fact=-1. 513 sp_zer=.true. 514 if ((ldirect.eq.1).and.(DRYDEP)) then 515 do jy=0,numygrid-1 516 do ix=0,numxgrid-1 517 517 if (drygrid(ix,jy).gt.smallnum) then 518 518 if (sp_zer.eqv..true.) then ! first non zero value … … 522 522 sp_zer=.false. 523 523 sp_fact=sp_fact*(-1) 524 525 526 527 528 529 530 531 524 endif 525 sp_count_r=sp_count_r+1 526 sparse_dump_r(sp_count_r)= & 527 sp_fact* & 528 1.e12*drygrid(ix,jy)/area(ix,jy) 529 sparse_dump_u(sp_count_r)= & 530 1.e12*drygridsigma(ix,jy)/area(ix,jy) 531 else ! concentration is zero 532 532 sp_zer=.true. 533 534 535 536 537 538 539 540 541 542 543 533 endif 534 end do 535 end do 536 else 537 sp_count_i=0 538 sp_count_r=0 539 endif 540 write(unitoutgridppt) sp_count_i 541 write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i) 542 write(unitoutgridppt) sp_count_r 543 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 544 544 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 545 545 546 546 547 ! Mixing ratios548 549 ! surf_only write only 1st layer550 551 552 553 554 547 ! Mixing ratios 548 549 ! surf_only write only 1st layer 550 551 sp_count_i=0 552 sp_count_r=0 553 sp_fact=-1. 554 sp_zer=.true. 555 555 do kz=1,1 556 556 do jy=0,numygrid-1 … … 563 563 sp_zer=.false. 564 564 sp_fact=sp_fact*(-1.) 565 566 567 568 569 570 571 572 573 574 575 576 565 endif 566 sp_count_r=sp_count_r+1 567 sparse_dump_r(sp_count_r)= & 568 sp_fact* & 569 1.e12*grid(ix,jy,kz) & 570 /volume(ix,jy,kz)/outnum* & 571 weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz) 572 sparse_dump_u(sp_count_r)= & 573 1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ & 574 outnum*weightair/weightmolar(ks)/ & 575 densityoutgrid(ix,jy,kz) 576 else ! concentration is zero 577 577 sp_zer=.true. 578 578 endif 579 579 end do 580 580 end do 581 581 end do 582 583 584 585 582 write(unitoutgridppt) sp_count_i 583 write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i) 584 write(unitoutgridppt) sp_count_r 585 write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r) 586 586 ! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r) 587 587 588 589 590 591 588 endif ! output for ppt 589 590 end do 591 end do 592 592 593 593 close(unitoutgridppt) … … 602 602 drygridtotal 603 603 604 ! Dump of receptor concentrations605 606 if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then607 write(unitoutreceptppt) itime608 do ks=1,nspec609 write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &610 weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)611 end do612 endif613 614 ! Dump of receptor concentrations615 616 if (numreceptor.gt.0) then617 write(unitoutrecept) itime618 do ks=1,nspec619 write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &620 i=1,numreceptor)621 end do622 endif623 624 625 626 ! Reinitialization of grid627 !*************************604 ! Dump of receptor concentrations 605 606 if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then 607 write(unitoutreceptppt) itime 608 do ks=1,nspec 609 write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* & 610 weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor) 611 end do 612 endif 613 614 ! Dump of receptor concentrations 615 616 if (numreceptor.gt.0) then 617 write(unitoutrecept) itime 618 do ks=1,nspec 619 write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, & 620 i=1,numreceptor) 621 end do 622 endif 623 624 625 626 ! Reinitialization of grid 627 !************************* 628 628 629 629 do ks=1,nspec 630 do kp=1,maxpointspec_act 631 do i=1,numreceptor 632 creceptor(i,ks)=0. 633 end do 634 do jy=0,numygrid-1 635 do ix=0,numxgrid-1 636 do l=1,nclassunc 637 do nage=1,nageclass 638 do kz=1,numzgrid 639 gridunc(ix,jy,kz,ks,kp,l,nage)=0. 640 end do 630 do kp=1,maxpointspec_act 631 do i=1,numreceptor 632 creceptor(i,ks)=0. 633 end do 634 do jy=0,numygrid-1 635 do ix=0,numxgrid-1 636 do l=1,nclassunc 637 do nage=1,nageclass 638 do kz=1,numzgrid 639 gridunc(ix,jy,kz,ks,kp,l,nage)=0. 641 640 end do 642 641 end do … … 645 644 end do 646 645 end do 646 end do 647 647 648 648 -
src/drydepokernel.f90
r1c0d5e6 r1c0d5e6 40 40 ! * 41 41 !***************************************************************************** 42 ! Changes:43 ! eso 10/2016: Added option to disregard kernel44 !45 !*****************************************************************************46 47 42 48 43 use unc_mod … … 52 47 implicit none 53 48 54 real(dep_prec), dimension(maxspec) :: deposit 55 real :: x,y,ddx,ddy,xl,yl,wx,wy,w 49 real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w 56 50 integer :: ix,jy,ixp,jyp,ks,nunc,nage,kp 57 51 … … 80 74 endif 81 75 82 ! If no kernel is used, direct attribution to grid cell83 !******************************************************84 85 if (lnokernel) then86 do ks=1,nspec87 if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then88 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &89 (jy.le.numygrid-1)) then90 drygridunc(ix,jy,ks,kp,nunc,nage)= &91 drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)92 end if93 end if94 end do95 else ! use kernel96 97 76 98 77 ! Determine mass fractions for four grid points 99 78 !********************************************** 100 do ks=1,nspec79 do ks=1,nspec 101 80 102 if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then81 if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then 103 82 104 83 if (.not.usekernel) then … … 108 87 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 109 88 (jy.le.numygrid-1)) then 110 111 89 w=wx*wy 90 drygridunc(ix,jy,ks,kp,nunc,nage)= & 112 91 drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 113 92 endif … … 137 116 endif ! deposit>0 138 117 139 end do 140 end if 118 end do 141 119 142 120 end subroutine drydepokernel -
src/drydepokernel_nest.f90
r4c64400 re200b7a 50 50 implicit none 51 51 52 real(dep_prec), dimension(maxspec) :: deposit 53 real :: x,y,ddx,ddy,xl,yl,wx,wy,w 52 real :: x,y,deposit(maxspec),ddx,ddy,xl,yl,wx,wy,w 54 53 integer :: ix,jy,ixp,jyp,ks,kp,nunc,nage 55 54 -
src/get_wetscav.f90
rc7e771d rc7e771d 72 72 integer(kind=1) :: clouds_v 73 73 integer :: ks, kp 74 integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count 75 74 integer :: inc_count, blc_count 76 75 ! integer :: n1,n2, icbot,ictop, indcloud !TEST 77 76 real :: S_i, act_temp, cl, cle ! in cloud scavenging … … 238 237 if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then 239 238 ! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then 240 blc_count (ks)=blc_count(ks)+1239 blc_count=blc_count+1 241 240 wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks) 242 241 … … 244 243 !********************************************************************************* 245 244 else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then 246 blc_count (ks)=blc_count(ks)+1245 blc_count=blc_count+1 247 246 248 247 !NIK 17.02.2015 … … 284 283 ! given in species file, or if gas and positive Henry's constant 285 284 if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then 286 inc_count (ks)=inc_count(ks)+1285 inc_count=inc_count+1 287 286 ! write(*,*) 'Incloud: ',inc_count 288 287 ! if negative coefficients (turned off) set to zero for use in equation -
src/gfs_mod.f90
r4c64400 r62e65c7 42 42 !********************************************* 43 43 44 ! integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64 45 integer,parameter :: nxmax=361,nymax=181,nuvzmax=64,nwzmax=64,nzmax=64 44 integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64 46 45 integer,parameter :: nxshift=0 ! for GFS or FNL 47 46 -
src/gridcheck_gfs.f90
r4c64400 r5f9d14a 421 421 !******************** 422 422 423 if (lroot) then 424 write(*,*) 425 write(*,*) 426 write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', & 427 nuvz,nwz 428 write(*,*) 429 write(*,'(a)') 'Mother domain:' 430 write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', & 431 xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx 432 write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range : ', & 433 ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy 434 write(*,*) 435 end if 423 write(*,*) 424 write(*,*) 425 write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', & 426 nuvz,nwz 427 write(*,*) 428 write(*,'(a)') 'Mother domain:' 429 write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', & 430 xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx 431 write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range : ', & 432 ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy 433 write(*,*) 434 436 435 437 436 ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL -
src/init_domainfill_mpi.f90
r16b61a5 r7999df47 47 47 ! * 48 48 !***************************************************************************** 49 ! MPI version:50 !51 ! -Root process allocates temporary arrays holding properties for52 ! all particles in the simulation53 ! -An index array is used to assign 1st particle to 1st process, 2nd particle54 ! to 2nd process and so on so that they are evenly distibuted geographically55 ! -Inititialization for domain-filling is done as in the serial code56 ! -Root process distributes particles evenly to other processes57 !58 !*****************************************************************************59 49 60 50 use point_mod … … 66 56 implicit none 67 57 68 ! ncolumn_mpi,numparttot_mpi ncolumn,numparttot per process 69 integer :: j,ix,jy,kz,ncolumn,numparttot,ncolumn_mpi,numparttot_mpi, arr_size 58 integer :: j,ix,jy,kz,ncolumn,numparttot 70 59 real :: gridarea(0:nymax-1),pp(nzmax),ylat,ylatp,ylatm,hzone 71 60 real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,fractus … … 77 66 78 67 integer :: idummy = -11 79 integer,allocatable,dimension(:) :: idx ! index array80 integer :: stride81 integer, parameter :: nullsize=082 68 logical :: first_call=.true. 83 69 84 ! Use different seed for each process ! TODO: not needed anymore 70 ! Use different seed for each process 85 71 if (first_call) then 86 72 idummy=idummy+mp_seed … … 116 102 117 103 118 ! This section only done by the root process119 !*******************************************120 if (lroot) then121 ! Arrays for particles to be released. Add a small number to npart(1) in case of122 ! round-off errors123 arr_size = npart(1) + mp_np124 allocate(itra1_tmp(arr_size),npoint_tmp(arr_size),nclass_tmp(arr_size),&125 & idt_tmp(arr_size),itramem_tmp(arr_size),itrasplit_tmp(arr_size),&126 & xtra1_tmp(arr_size),ytra1_tmp(arr_size),ztra1_tmp(arr_size),&127 & xmass1_tmp(arr_size, maxspec))128 129 ! Index array for particles. This is to avoid having particles130 ! near edges of domain all on one process.131 !****************************************************************************132 allocate(idx(npart(1)))133 stride = npart(1)/mp_partgroup_np134 135 jj=0136 do j=1, stride137 do i=0, mp_partgroup_np-1138 jj = jj+1139 if (jj.gt.npart(1)) exit140 idx(jj) = i*stride+j141 end do142 end do143 144 ! Add extra indices if npart(1) not evenly divisible by number of processes145 do i=1, mod(npart(1),mp_partgroup_np)146 jj = jj+1147 if (jj.gt.npart(1)) exit148 idx(jj) = jj149 end do150 151 ! Initialize all particles as non-existent152 itra1_tmp(:)=-999999999153 154 104 ! Calculate area of grid cell with formula M=2*pi*R*h*dx/360, 155 105 ! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90 156 106 !************************************************************ 157 107 158 do jy=ny_sn(1),ny_sn(2) ! loop about latitudes 159 ylat=ylat0+real(jy)*dy 160 ylatp=ylat+0.5*dy 161 ylatm=ylat-0.5*dy 162 if ((ylatm.lt.0).and.(ylatp.gt.0.)) then 163 hzone=1./dyconst 108 do jy=ny_sn(1),ny_sn(2) ! loop about latitudes 109 ylat=ylat0+real(jy)*dy 110 ylatp=ylat+0.5*dy 111 ylatm=ylat-0.5*dy 112 if ((ylatm.lt.0).and.(ylatp.gt.0.)) then 113 hzone=1./dyconst 114 else 115 cosfactp=cos(ylatp*pih)*r_earth 116 cosfactm=cos(ylatm*pih)*r_earth 117 if (cosfactp.lt.cosfactm) then 118 hzone=sqrt(r_earth**2-cosfactp**2)- & 119 sqrt(r_earth**2-cosfactm**2) 164 120 else 165 cosfactp=cos(ylatp*pih)*r_earth 166 cosfactm=cos(ylatm*pih)*r_earth 167 if (cosfactp.lt.cosfactm) then 168 hzone=sqrt(r_earth**2-cosfactp**2)- & 169 sqrt(r_earth**2-cosfactm**2) 170 else 171 hzone=sqrt(r_earth**2-cosfactm**2)- & 172 sqrt(r_earth**2-cosfactp**2) 173 endif 121 hzone=sqrt(r_earth**2-cosfactm**2)- & 122 sqrt(r_earth**2-cosfactp**2) 174 123 endif 175 gridarea(jy)=2.*pi*r_earth*hzone*dx/360. 176 end do 124 endif 125 gridarea(jy)=2.*pi*r_earth*hzone*dx/360. 126 end do 177 127 178 128 ! Do the same for the south pole 179 129 180 181 182 183 184 185 186 187 188 189 130 if (sglobal) then 131 ylat=ylat0 132 ylatp=ylat+0.5*dy 133 ylatm=ylat 134 cosfactm=0. 135 cosfactp=cos(ylatp*pih)*r_earth 136 hzone=sqrt(r_earth**2-cosfactm**2)- & 137 sqrt(r_earth**2-cosfactp**2) 138 gridarea(0)=2.*pi*r_earth*hzone*dx/360. 139 endif 190 140 191 141 ! Do the same for the north pole 192 142 193 194 195 196 197 198 199 200 201 202 143 if (nglobal) then 144 ylat=ylat0+real(nymin1)*dy 145 ylatp=ylat 146 ylatm=ylat-0.5*dy 147 cosfactp=0. 148 cosfactm=cos(ylatm*pih)*r_earth 149 hzone=sqrt(r_earth**2-cosfactp**2)- & 150 sqrt(r_earth**2-cosfactm**2) 151 gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360. 152 endif 203 153 204 154 … … 206 156 !********************************************************************* 207 157 208 209 210 211 212 213 colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)214 colmasstotal=colmasstotal+colmass(ix,jy)215 216 end do 158 colmasstotal=0. 159 do jy=ny_sn(1),ny_sn(2) ! loop about latitudes 160 do ix=nx_we(1),nx_we(2) ! loop about longitudes 161 pp(1)=rho(ix,jy,1,1)*r_air*tt(ix,jy,1,1) 162 pp(nz)=rho(ix,jy,nz,1)*r_air*tt(ix,jy,nz,1) 163 ! Each MPI process is assigned an equal share of particles 164 colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)/mp_partgroup_np 165 colmasstotal=colmasstotal+colmass(ix,jy) 166 217 167 end do 218 219 write(*,*) 'Atm. mass: ',colmasstotal 220 221 222 if (ipin.eq.0) numpart=0 168 end do 169 170 if (lroot) write(*,*) 'Atm. mass: ',colmasstotal 171 172 173 if (ipin.eq.0) numpart=0 223 174 224 175 ! Determine the particle positions 225 176 !********************************* 226 177 227 228 229 230 231 232 ncolumn=nint(0.999*real(npart(1))*colmass(ix,jy)/ &233 234 235 178 numparttot=0 179 numcolumn=0 180 do jy=ny_sn(1),ny_sn(2) ! loop about latitudes 181 ylat=ylat0+real(jy)*dy 182 do ix=nx_we(1),nx_we(2) ! loop about longitudes 183 ncolumn=nint(0.999*real(npart(1)/mp_partgroup_np)*colmass(ix,jy)/ & 184 colmasstotal) 185 if (ncolumn.eq.0) goto 30 186 if (ncolumn.gt.numcolumn) numcolumn=ncolumn 236 187 237 188 ! Calculate pressure at the altitudes of model surfaces, using the air density … … 239 190 !***************************************************************************** 240 191 241 242 243 244 245 246 247 248 249 250 192 do kz=1,nz 193 pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1) 194 end do 195 196 197 deltacol=(pp(1)-pp(nz))/real(ncolumn) 198 pnew=pp(1)+deltacol/2. 199 jj=0 200 do j=1,ncolumn 201 jj=jj+1 251 202 252 203 … … 257 208 258 209 259 260 261 262 263 264 265 266 267 268 269 210 if (ncolumn.gt.20) then 211 pnew=pnew-deltacol 212 else 213 pnew=pp(1)-ran1(idummy)*(pp(1)-pp(nz)) 214 endif 215 216 do kz=1,nz-1 217 if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then 218 dz1=pp(kz)-pnew 219 dz2=pnew-pp(kz+1) 220 dz=1./(dz1+dz2) 270 221 271 222 ! Assign particle position … … 273 224 ! Do the following steps only if particles are not read in from previous model run 274 225 !***************************************************************************** 275 276 xtra1_tmp(idx(numpart+jj))=real(ix)-0.5+ran1(idummy)277 if (ix.eq.0) xtra1_tmp(idx(numpart+jj))=ran1(idummy)278 if (ix.eq.nxmin1) xtra1_tmp(idx(numpart+jj))= &279 280 ytra1_tmp(idx(numpart+jj))=real(jy)-0.5+ran1(idummy)281 ztra1_tmp(idx(numpart+jj))=(height(kz)*dz2+height(kz+1)*dz1)*dz282 if (ztra1_tmp(idx(numpart+jj)).gt.height(nz)-0.5) &283 ztra1_tmp(idx(numpart+jj))=height(nz)-0.5226 if (ipin.eq.0) then 227 xtra1(numpart+jj)=real(ix)-0.5+ran1(idummy) 228 if (ix.eq.0) xtra1(numpart+jj)=ran1(idummy) 229 if (ix.eq.nxmin1) xtra1(numpart+jj)= & 230 real(nxmin1)-ran1(idummy) 231 ytra1(numpart+jj)=real(jy)-0.5+ran1(idummy) 232 ztra1(numpart+jj)=(height(kz)*dz2+height(kz+1)*dz1)*dz 233 if (ztra1(numpart+jj).gt.height(nz)-0.5) & 234 ztra1(numpart+jj)=height(nz)-0.5 284 235 285 236 286 237 ! Interpolate PV to the particle position 287 238 !**************************************** 288 ixm=int(xtra1_tmp(idx(numpart+jj)))289 jym=int(ytra1_tmp(idx(numpart+jj)))290 291 292 ddx=xtra1_tmp(idx(numpart+jj))-real(ixm)293 ddy=ytra1_tmp(idx(numpart+jj))-real(jym)294 295 296 297 298 299 300 301 if (height(i).gt.ztra1_tmp(idx(numpart+jj))) then302 303 304 305 306 307 6 308 dz1=ztra1_tmp(idx(numpart+jj))-height(indzm)309 dz2=height(indzp)-ztra1_tmp(idx(numpart+jj))310 311 312 313 314 315 316 317 318 319 239 ixm=int(xtra1(numpart+jj)) 240 jym=int(ytra1(numpart+jj)) 241 ixp=ixm+1 242 jyp=jym+1 243 ddx=xtra1(numpart+jj)-real(ixm) 244 ddy=ytra1(numpart+jj)-real(jym) 245 rddx=1.-ddx 246 rddy=1.-ddy 247 p1=rddx*rddy 248 p2=ddx*rddy 249 p3=rddx*ddy 250 p4=ddx*ddy 251 do i=2,nz 252 if (height(i).gt.ztra1(numpart+jj)) then 253 indzm=i-1 254 indzp=i 255 goto 6 256 endif 257 end do 258 6 continue 259 dz1=ztra1(numpart+jj)-height(indzm) 260 dz2=height(indzp)-ztra1(numpart+jj) 261 dz=1./(dz1+dz2) 262 do in=1,2 263 indzh=indzm+in-1 264 y1(in)=p1*pv(ixm,jym,indzh,1) & 265 +p2*pv(ixp,jym,indzh,1) & 266 +p3*pv(ixm,jyp,indzh,1) & 267 +p4*pv(ixp,jyp,indzh,1) 268 end do 269 pvpart=(dz2*y1(1)+dz1*y1(2))*dz 270 if (ylat.lt.0.) pvpart=-1.*pvpart 320 271 321 272 … … 323 274 !***************************************************************************** 324 275 325 if (((ztra1_tmp(idx(numpart+jj)).gt.3000.).and. &326 276 if (((ztra1(numpart+jj).gt.3000.).and. & 277 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 327 278 328 279 ! Assign certain properties to the particle 329 280 !****************************************** 330 nclass_tmp(idx(numpart+jj))=min(int(ran1(idummy)* & 331 real(nclassunc))+1,nclassunc) 332 numparticlecount=numparticlecount+1 333 npoint_tmp(idx(numpart+jj))=numparticlecount 334 idt_tmp(idx(numpart+jj))=mintime 335 itra1_tmp(idx(numpart+jj))=0 336 itramem_tmp(idx(numpart+jj))=0 337 itrasplit_tmp(idx(numpart+jj))=itra1_tmp(idx(numpart+jj))+ldirect* & 338 itsplit 339 xmass1_tmp(idx(numpart+jj),1)=colmass(ix,jy)/real(ncolumn) 340 if (mdomainfill.eq.2) xmass1_tmp(idx(numpart+jj),1)= & 341 xmass1_tmp(idx(numpart+jj),1)*pvpart*48./29.*ozonescale/10.**9 342 else 343 jj=jj-1 344 endif 281 nclass(numpart+jj)=min(int(ran1(idummy)* & 282 real(nclassunc))+1,nclassunc) 283 numparticlecount=numparticlecount+1 284 npoint(numpart+jj)=numparticlecount 285 idt(numpart+jj)=mintime 286 itra1(numpart+jj)=0 287 itramem(numpart+jj)=0 288 itrasplit(numpart+jj)=itra1(numpart+jj)+ldirect* & 289 itsplit 290 xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn) 291 if (mdomainfill.eq.2) xmass1(numpart+jj,1)= & 292 xmass1(numpart+jj,1)*pvpart*48./29.*ozonescale/10.**9 293 else 294 jj=jj-1 345 295 endif 346 296 endif 347 end do297 endif 348 298 end do 349 numparttot=numparttot+ncolumn350 if (ipin.eq.0) numpart=numpart+jj351 30 continue352 299 end do 300 numparttot=numparttot+ncolumn 301 if (ipin.eq.0) numpart=numpart+jj 302 30 continue 353 303 end do 304 end do 354 305 355 306 … … 357 308 !***************************************************** 358 309 359 ! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier 360 if (numpart.gt.maxpart) then 361 write(*,*) 'numpart too large: change source in init_atm_mass.f' 362 write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart 363 endif 364 365 366 xmassperparticle=colmasstotal/real(numparttot) 310 if (numpart.gt.maxpart) then 311 write(*,*) 'numpart too large: change source in init_atm_mass.f' 312 write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart 313 endif 314 315 316 xmassperparticle=colmasstotal/real(numparttot) 367 317 368 318 … … 370 320 !*********************************************** 371 321 372 ! do j=1,numpart 373 do j=1,npoint(1) 374 if ((xtra1_tmp(j).lt.0.).or.(xtra1_tmp(j).ge.real(nxmin1)).or. & 375 (ytra1_tmp(j).lt.0.).or.(ytra1_tmp(j).ge.real(nymin1))) then 376 itra1_tmp(j)=-999999999 377 endif 378 end do 322 do j=1,numpart 323 if ((xtra1(j).lt.0.).or.(xtra1(j).ge.real(nxmin1)).or. & 324 (ytra1(j).lt.0.).or.(ytra1(j).ge.real(nymin1))) then 325 itra1(j)=-999999999 326 endif 327 end do 379 328 380 329 … … 391 340 !**************************************************************************** 392 341 393 394 write(*,*) 'Total number of particles at model start: ',numpart395 write(*,*) 'Maximum number of particles per column: ',numcolumn396 write(*,*) 'If ',fractus,' <1, better use more particles'397 398 399 400 401 ncolumn=nint(0.999/fractus*real(npart(1))*colmass(ix,jy) &402 403 404 342 fractus=real(numcolumn)/real(nz) 343 write(*,*) 'Total number of particles at model start: ',numpart*mp_partgroup_np 344 write(*,*) 'Maximum number of particles per column: ',numcolumn*mp_partgroup_np 345 write(*,*) 'If ',fractus*mp_partgroup_np,' <1, better use more particles' 346 fractus=sqrt(max(fractus,1.))/2. 347 348 do jy=ny_sn(1),ny_sn(2) ! loop about latitudes 349 do ix=nx_we(1),nx_we(2) ! loop about longitudes 350 ncolumn=nint(0.999/fractus*real(npart(1)/mp_partgroup_np)*colmass(ix,jy) & 351 /colmasstotal) 352 if (ncolumn.gt.maxcolumn) stop 'maxcolumn too small' 353 if (ncolumn.eq.0) goto 80 405 354 406 355 … … 410 359 !************************************************************************ 411 360 412 413 414 415 361 if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn 362 if (ix.eq.nx_we(2)) numcolumn_we(2,jy)=ncolumn 363 if (jy.eq.ny_sn(1)) numcolumn_sn(1,ix)=ncolumn 364 if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn 416 365 417 366 ! Calculate pressure at the altitudes of model surfaces, using the air density … … 419 368 !***************************************************************************** 420 369 421 422 423 370 do kz=1,nz 371 pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1) 372 end do 424 373 425 374 ! Determine the reference starting altitudes 426 375 !******************************************* 427 376 428 429 430 431 432 433 434 435 436 437 438 377 deltacol=(pp(1)-pp(nz))/real(ncolumn) 378 pnew=pp(1)+deltacol/2. 379 do j=1,ncolumn 380 pnew=pnew-deltacol 381 do kz=1,nz-1 382 if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then 383 dz1=pp(kz)-pnew 384 dz2=pnew-pp(kz+1) 385 dz=1./(dz1+dz2) 386 zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz 387 if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5 439 388 440 389 ! Memorize vertical positions where particles are introduced … … 442 391 !*********************************************************** 443 392 444 445 446 447 393 if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition 394 if (ix.eq.nx_we(2)) zcolumn_we(2,jy,j)=zposition 395 if (jy.eq.ny_sn(1)) zcolumn_sn(1,ix,j)=zposition 396 if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition 448 397 449 398 ! Initialize mass that has accumulated at boundary to zero 450 399 !********************************************************* 451 400 452 acc_mass_we(1,jy,j)=0. 453 acc_mass_we(2,jy,j)=0. 454 acc_mass_sn(1,jy,j)=0. 455 acc_mass_sn(2,jy,j)=0. 456 endif 457 end do 401 acc_mass_we(1,jy,j)=0. 402 acc_mass_we(2,jy,j)=0. 403 acc_mass_sn(1,jy,j)=0. 404 acc_mass_sn(2,jy,j)=0. 405 endif 458 406 end do 459 80 continue460 407 end do 408 80 continue 461 409 end do 410 end do 462 411 463 412 ! If particles shall be read in to continue an existing run, … … 466 415 !*************************************************************************** 467 416 468 ! eso TODO: only needed for root process 469 if (ipin.eq.1) then 470 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 471 form='unformatted') 472 read(unitboundcond) numcolumn_we,numcolumn_sn, & 473 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn 474 close(unitboundcond) 475 endif 476 477 numpart = npart(1)/mp_partgroup_np 478 if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1 479 480 else ! Allocate dummy arrays for receiving processes 481 allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),& 482 & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),& 483 & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),& 484 & xmass1_tmp(nullsize, nullsize)) 485 486 end if ! end if(lroot) 487 488 489 ! Distribute particles to other processes (numpart is 'per-process', not total) 490 call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr) 491 ! eso TODO: xmassperparticle: not necessary to send 492 call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr) 493 call mpif_send_part_properties(numpart) 494 495 ! Deallocate the temporary arrays used for all particles 496 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,& 497 & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp) 417 ! :TODO: eso: parallelize 418 if (ipin.eq.1) then 419 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 420 form='unformatted') 421 read(unitboundcond) numcolumn_we,numcolumn_sn, & 422 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn 423 close(unitboundcond) 424 endif 425 426 498 427 499 428 -
src/makefile
rd1a8707 rd1a8707 423 423 par_mod.o point_mod.o unc_mod.o xmass_mod.o 424 424 timemanager_mpi.o: com_mod.o flux_mod.o mpi_mod.o oh_mod.o outg_mod.o \ 425 par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o425 par_mod.o point_mod.o unc_mod.o xmass_mod.o 426 426 unc_mod.o: par_mod.o 427 427 verttransform.o: cmapf_mod.o com_mod.o par_mod.o -
src/mean_mod.f90
r4c64400 r6a678e3 29 29 module procedure mean_sp 30 30 module procedure mean_dp 31 module procedure mean_mixed_dss 32 module procedure mean_mixed_dsd 33 31 module procedure mean_mixed_prec 34 32 end interface mean 35 33 … … 66 64 ! real(sp),parameter :: eps=1.0e-30 67 65 68 integer,intent(in) :: number69 66 real(sp), intent(in) :: x_sp(number) 70 67 real(sp), intent(out) ::xm,xs 68 integer,intent(in) :: number 71 69 real(sp) :: xl,xq,xaux 72 70 real(sp),parameter :: eps=1.0e-30 … … 118 116 implicit none 119 117 120 integer,intent(in) :: number121 118 real(dp), intent(in) :: x_dp(number) 122 119 real(dp), intent(out) ::xm,xs 120 integer,intent(in) :: number 123 121 real(dp) :: xl,xq,xaux 124 122 real(dp),parameter :: eps=1.0e-30 … … 144 142 end subroutine mean_dp 145 143 146 subroutine mean_mixed_ dss(x_dp,xm,xs,number)144 subroutine mean_mixed_prec(x_dp,xm,xs,number) 147 145 148 146 !***************************************************************************** … … 152 150 ! AUTHOR: Andreas Stohl, 25 January 1994 * 153 151 ! * 154 ! Mixed precision version ESO 2016 (dp in , sp out, sp out)*152 ! Mixed precision version ESO 2016 (dp input, sp output) * 155 153 !***************************************************************************** 156 154 ! * … … 170 168 implicit none 171 169 172 integer,intent(in) :: number173 170 real(dp), intent(in) :: x_dp(number) 174 171 real(sp), intent(out) ::xm,xs 172 integer,intent(in) :: number 175 173 real(sp) :: xl,xq,xaux 176 174 real(sp),parameter :: eps=1.0e-30 … … 194 192 endif 195 193 196 end subroutine mean_mixed_dss 197 198 subroutine mean_mixed_dsd(x_dp,xm,xs_dp,number) 199 200 !***************************************************************************** 201 ! * 202 ! This subroutine calculates mean and standard deviation of a given element.* 203 ! * 204 ! AUTHOR: Andreas Stohl, 25 January 1994 * 205 ! * 206 ! Mixed precision version ESO 2016 (dp in, sp out, dp out) * 207 !***************************************************************************** 208 ! * 209 ! Variables: * 210 ! x_dp(number) field of input data * 211 ! xm mean * 212 ! xs_dp standard deviation * 213 ! number number of elements of field x_dp * 214 ! * 215 ! Constants: * 216 ! eps tiny number * 217 ! * 218 !***************************************************************************** 219 220 use par_mod, only: sp,dp 221 222 implicit none 223 224 integer,intent(in) :: number 225 real(dp), intent(in) :: x_dp(number) 226 real(sp), intent(out) ::xm 227 real(dp), intent(out) ::xs_dp 228 real(dp) :: xl,xq,xaux 229 real(dp),parameter :: eps=1.0e-30_dp 230 integer :: i 231 232 xl=0._dp 233 xq=0._dp 234 do i=1,number 235 xl=xl+x_dp(i) 236 xq=xq+x_dp(i)*x_dp(i) 237 end do 238 239 xm=xl/real(number,kind=sp) 240 241 xaux=xq-xl*xl/real(number,kind=dp) 242 243 if (xaux.lt.eps) then 244 xs_dp=0._dp 245 else 246 xs_dp=sqrt(xaux/real(number-1,kind=dp)) 247 endif 248 249 end subroutine mean_mixed_dsd 250 194 end subroutine mean_mixed_prec 251 195 end module mean_mod -
src/mpi_mod.f90
r4c64400 r79abee9 90 90 ! MPI tags/requests for send/receive operation 91 91 integer :: tm1 92 integer, parameter :: nvar_async=26 92 integer, parameter :: nvar_async=26 !27 !29 :DBG: 93 93 !integer, dimension(:), allocatable :: tags 94 94 integer, dimension(:), allocatable :: reqs 95 96 ! Status array used for certain MPI operations (MPI_RECV)97 integer, dimension(MPI_STATUS_SIZE) :: mp_status98 95 99 96 … … 152 149 integer, private :: dat_lun 153 150 154 ! Temporary arrays for particles (allocated and deallocated as needed)155 integer, allocatable, dimension(:) :: nclass_tmp, npoint_tmp, itra1_tmp, idt_tmp, &156 & itramem_tmp, itrasplit_tmp157 real(kind=dp), allocatable, dimension(:) :: xtra1_tmp, ytra1_tmp158 real, allocatable, dimension(:) :: ztra1_tmp159 real, allocatable, dimension(:,:) :: xmass1_tmp160 161 ! mp_redist_fract Exchange particles between processes if relative difference162 ! is larger. A good value is between 0.0 and 0.5163 ! mp_maxpart_factor Allocate more memory per process than strictly needed164 ! (safety factor). Recommended value between 1.5 and 2.5165 ! mp_min_redist Do not redistribute particles if below this limit166 real, parameter :: mp_redist_fract=0.2, mp_maxpart_factor=1.5167 integer,parameter :: mp_min_redist=100000168 169 170 151 contains 171 152 … … 214 195 if (dep_prec==dp) then 215 196 mp_cp = MPI_REAL8 216 ! TODO: write info message for serial version as well197 ! TODO: write info message for serial version as well 217 198 if (lroot.and.verbosity>0) write(*,*) 'Using double precision for deposition fields' 218 199 else if (dep_prec==sp) then … … 251 232 write(*,FMT='(80("#"))') 252 233 end if 253 lmp_sync=.true. 234 lmp_sync=.true. ! :DBG: eso fix this... 254 235 end if 255 236 … … 330 311 331 312 ! Set maxpart per process 332 ! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution 333 maxpart_mpi=int(m p_maxpart_factor*real(maxpart)/real(mp_partgroup_np))313 if (mp_partid.lt.mod(maxpart,mp_partgroup_np)) addmaxpart=1 314 maxpart_mpi=int(maxpart/mp_partgroup_np)+addmaxpart 334 315 335 316 ! Do not allocate particle data arrays for readwind process … … 340 321 ! Set random seed for each non-root process 341 322 if (mp_pid.gt.0) then 323 ! if (mp_pid.ge.0) then 324 ! call system_clock(s) 342 325 s = 244 343 326 mp_seed = -abs(mod((s*181)*((mp_pid-83)*359), 104729)) 344 327 end if 328 if (mp_dev_mode) write(*,*) 'PID, mp_seed: ',mp_pid, mp_seed 345 329 if (mp_dbg_mode) then 330 ! :DBG: For debugging, set all seed to 0 and maxrand to e.g. 20 346 331 mp_seed=0 347 332 if (lroot) write(*,*) 'MPI: setting seed=0' … … 469 454 470 455 471 subroutine mpif_send_part_properties(num_part) 472 !*********************************************************************** 473 ! Distribute particle properties from root to all processes. 474 ! 475 ! Used by init_domainfill_mpi 476 ! 477 ! Variables: 478 ! 479 ! num_part input, number of particles per process (rounded up) 456 subroutine mpif_tm_send_vars 457 !*********************************************************************** 458 ! Distribute particle variables from pid0 to all processes. 459 ! Called from timemanager 460 ! *NOT IN USE* at the moment, but can be useful for debugging 480 461 ! 481 462 !*********************************************************************** … … 484 465 implicit none 485 466 486 integer,intent(in) :: num_part 487 integer :: i,jj, addone 488 489 ! Exit if too many particles 490 if (num_part.gt.maxpart_mpi) then 491 write(*,*) '#####################################################' 492 write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED ####' 493 write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE ####' 494 write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####' 495 write(*,*) '#### OR INCREASE MAXPART. ####' 496 write(*,*) '#####################################################' 497 ! call MPI_FINALIZE(mp_ierr) 498 stop 499 end if 500 467 integer :: i 501 468 502 469 ! Time for MPI communications … … 504 471 if (mp_measure_time) call mpif_mtime('commtime',0) 505 472 506 ! Distribute variables, send from pid 0 to other processes (including itself)507 !****************************************************************************508 509 call MPI_SCATTER(nclass_tmp,num_part,MPI_INTEGER,nclass,&510 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)511 if (mp_ierr /= 0) goto 600512 call MPI_SCATTER(npoint_tmp,num_part,MPI_INTEGER,npoint,&513 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)514 if (mp_ierr /= 0) goto 600515 call MPI_SCATTER(itra1_tmp,num_part,MPI_INTEGER,itra1,&516 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)517 if (mp_ierr /= 0) goto 600518 call MPI_SCATTER(idt_tmp,num_part,MPI_INTEGER,idt,&519 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)520 if (mp_ierr /= 0) goto 600521 call MPI_SCATTER(itramem_tmp,num_part,MPI_INTEGER,itramem,&522 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)523 if (mp_ierr /= 0) goto 600524 call MPI_SCATTER(itrasplit_tmp,num_part,MPI_INTEGER,itrasplit,&525 &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)526 if (mp_ierr /= 0) goto 600527 call MPI_SCATTER(xtra1_tmp,num_part,mp_dp,xtra1,&528 &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)529 if (mp_ierr /= 0) goto 600530 call MPI_SCATTER(ytra1_tmp,num_part,mp_dp,ytra1,&531 &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)532 if (mp_ierr /= 0) goto 600533 call MPI_SCATTER(ztra1_tmp,num_part,mp_sp,ztra1,&534 &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)535 if (mp_ierr /= 0) goto 600536 do i=1,nspec537 call MPI_SCATTER(xmass1_tmp(:,i),num_part,mp_sp,xmass1(:,i),&538 &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)539 if (mp_ierr /= 0) goto 600540 end do541 542 if (mp_measure_time) call mpif_mtime('commtime',1)543 544 goto 601545 546 600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr547 stop548 549 ! After the transfer of particles it it possible that the value of550 ! "numpart" is larger than the actual used, so we reduce it if there are551 ! invalid particles at the end of the arrays552 601 do i=num_part, 1, -1553 if (itra1(i).eq.-999999999) then554 numpart=numpart-1555 else556 exit557 end if558 end do559 560 561 !601 end subroutine mpif_send_part_properties562 end subroutine mpif_send_part_properties563 564 565 subroutine mpif_calculate_part_redist(itime)566 !***********************************************************************567 ! Calculate number of particles to redistribute between processes. This routine568 ! can be called at regular time intervals to keep a level number of569 ! particles on each process.570 !571 ! First, all processes report their local "numpart" to each other, which is572 ! stored in array "numpart_mpi(np)". This array is sorted from low to573 ! high values, along with a corresponding process ID array "idx_arr(np)".574 ! If the relative difference between the highest and lowest "numpart_mpi" value575 ! is above a threshold, particles are transferred from process idx_arr(np-1)576 ! to process idx_arr(0). Repeat for processes idx_arr(np-i) and idx_arr(i)577 ! for all valid i.578 ! Note: If np is an odd number, the process with the median579 ! number of particles is left unchanged580 !581 ! VARIABLES582 ! itime input, current time583 !***********************************************************************584 use com_mod585 586 implicit none587 588 integer, intent(in) :: itime589 real :: pmin,z590 integer :: i,jj,nn, num_part=1,m,imin, num_trans591 logical :: first_iter592 integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr593 real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this594 595 ! Exit if running with only 1 process596 !************************************************************************597 if (mp_np.eq.1) return598 599 ! All processes exchange information on number of particles600 !****************************************************************************601 allocate(numparticles_mpi(0:mp_partgroup_np-1), &602 &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))603 604 call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &605 & 1, MPI_INTEGER, mp_comm_used, mp_ierr)606 607 608 ! Sort from lowest to highest609 ! Initial guess: correct order610 sorted(:) = numparticles_mpi(:)611 do i=0, mp_partgroup_np-1612 idx_arr(i) = i613 end do614 615 ! For each successive element in index array, see if a lower value exists616 do i=0, mp_partgroup_np-2617 pmin=sorted(i)618 imin=idx_arr(i)619 m=i+1620 do jj=m, mp_partgroup_np-1621 if (pmin.le.sorted(jj)) cycle622 z=pmin623 pmin=sorted(jj)624 sorted(jj)=z625 626 nn=imin627 imin=idx_arr(jj)628 idx_arr(jj)=nn629 end do630 sorted(i)=pmin631 idx_arr(i)=imin632 end do633 634 ! For each pair of processes, transfer particles if the difference is above a635 ! limit, and numpart at sending process large enough636 637 m=mp_partgroup_np-1 ! index for last sorted process (most particles)638 do i=0,mp_partgroup_np/2-1639 num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))640 if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then641 if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&642 & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then643 ! DBG644 ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &645 ! &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi646 ! DBG647 call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)648 end if649 end if650 m=m-1651 end do652 653 deallocate(numparticles_mpi, idx_arr, sorted)654 655 end subroutine mpif_calculate_part_redist656 657 658 subroutine mpif_redist_part(itime, src_proc, dest_proc, num_trans)659 !***********************************************************************660 ! Transfer particle properties between two arbitrary processes.661 !662 ! VARIABLES663 !664 ! itime input, current time665 ! src_proc input, ID of source process666 ! dest_proc input, ID of destination process667 ! num_trans input, number of particles to transfer668 !669 !************************************************************************670 use com_mod !TODO: ,only: nclass etc671 672 implicit none673 674 integer, intent(in) :: itime, src_proc, dest_proc, num_trans675 integer :: ll, ul ! lower and upper indices in arrays676 integer :: arr_sz ! size of temporary arrays677 integer :: mtag ! MPI message tag678 integer :: i, j, minpart, ipart, maxnumpart679 680 ! Check for invalid input arguments681 !**********************************682 if (src_proc.eq.dest_proc) then683 write(*,*) '<mpi_mod::mpif_redist_part>: Error: &684 &src_proc.eq.dest_proc'685 stop686 end if687 688 ! Measure time for MPI communications689 !************************************690 if (mp_measure_time) call mpif_mtime('commtime',0)691 692 ! Send particles693 !***************694 if (mp_partid.eq.src_proc) then695 mtag = 2000696 697 ! Array indices for the transferred particles698 ll=numpart-num_trans+1699 ul=numpart700 701 if (mp_dev_mode) then702 write(*,FMT='(72("#"))')703 write(*,*) "Sending ", num_trans, "particles (from/to)", src_proc, dest_proc704 write(*,*) "numpart before: ", numpart705 end if706 707 call MPI_SEND(nclass(ll:ul),num_trans,&708 & MPI_INTEGER,dest_proc,mtag+1,mp_comm_used,mp_ierr)709 710 call MPI_SEND(npoint(ll:ul),num_trans,&711 & MPI_INTEGER,dest_proc,mtag+2,mp_comm_used,mp_ierr)712 713 call MPI_SEND(itra1(ll:ul),num_trans, &714 & MPI_INTEGER,dest_proc,mtag+3,mp_comm_used,mp_ierr)715 716 call MPI_SEND(idt(ll:ul),num_trans, &717 & MPI_INTEGER,dest_proc,mtag+4,mp_comm_used,mp_ierr)718 719 call MPI_SEND(itramem(ll:ul),num_trans, &720 & MPI_INTEGER,dest_proc,mtag+5,mp_comm_used,mp_ierr)721 722 call MPI_SEND(itrasplit(ll:ul),num_trans,&723 & MPI_INTEGER,dest_proc,mtag+6,mp_comm_used,mp_ierr)724 725 call MPI_SEND(xtra1(ll:ul),num_trans, &726 & mp_dp,dest_proc,mtag+7,mp_comm_used,mp_ierr)727 728 call MPI_SEND(ytra1(ll:ul),num_trans,&729 & mp_dp,dest_proc,mtag+8,mp_comm_used,mp_ierr)730 731 call MPI_SEND(ztra1(ll:ul),num_trans,&732 & mp_sp,dest_proc,mtag+9,mp_comm_used,mp_ierr)733 734 do j=1,nspec735 call MPI_SEND(xmass1(ll:ul,j),num_trans,&736 & mp_sp,dest_proc,mtag+(9+j),mp_comm_used,mp_ierr)737 end do738 739 ! Terminate transferred particles, update numpart740 itra1(ll:ul) = -999999999741 742 numpart = numpart-num_trans743 744 if (mp_dev_mode) then745 write(*,*) "numpart after: ", numpart746 write(*,FMT='(72("#"))')747 end if748 749 else if (mp_partid.eq.dest_proc) then750 751 ! Receive particles752 !******************753 mtag = 2000754 ! Array indices for the transferred particles755 ll=numpart+1756 ul=numpart+num_trans757 758 if (mp_dev_mode) then759 write(*,FMT='(72("#"))')760 write(*,*) "Receiving ", num_trans, "particles (from/to)", src_proc, dest_proc761 write(*,*) "numpart before: ", numpart762 end if763 764 ! We could receive the data directly at nclass(ll:ul) etc., but this leaves765 ! vacant spaces at lower indices. Using temporary arrays instead.766 arr_sz = ul-ll+1767 allocate(itra1_tmp(arr_sz),npoint_tmp(arr_sz),nclass_tmp(arr_sz),&768 & idt_tmp(arr_sz),itramem_tmp(arr_sz),itrasplit_tmp(arr_sz),&769 & xtra1_tmp(arr_sz),ytra1_tmp(arr_sz),ztra1_tmp(arr_sz),&770 & xmass1_tmp(arr_sz, maxspec))771 772 call MPI_RECV(nclass_tmp,num_trans,&773 & MPI_INTEGER,src_proc,mtag+1,mp_comm_used,mp_status,mp_ierr)774 775 call MPI_RECV(npoint_tmp,num_trans,&776 & MPI_INTEGER,src_proc,mtag+2,mp_comm_used,mp_status,mp_ierr)777 778 call MPI_RECV(itra1_tmp,num_trans, &779 & MPI_INTEGER,src_proc,mtag+3,mp_comm_used,mp_status,mp_ierr)780 781 call MPI_RECV(idt_tmp,num_trans, &782 & MPI_INTEGER,src_proc,mtag+4,mp_comm_used,mp_status,mp_ierr)783 784 call MPI_RECV(itramem_tmp,num_trans, &785 & MPI_INTEGER,src_proc,mtag+5,mp_comm_used,mp_status,mp_ierr)786 787 call MPI_RECV(itrasplit_tmp,num_trans,&788 & MPI_INTEGER,src_proc,mtag+6,mp_comm_used,mp_status,mp_ierr)789 790 call MPI_RECV(xtra1_tmp,num_trans, &791 & mp_dp,src_proc,mtag+7,mp_comm_used,mp_status,mp_ierr)792 793 call MPI_RECV(ytra1_tmp,num_trans,&794 & mp_dp,src_proc,mtag+8,mp_comm_used,mp_status,mp_ierr)795 796 call MPI_RECV(ztra1_tmp,num_trans,&797 & mp_sp,src_proc,mtag+9,mp_comm_used,mp_status,mp_ierr)798 799 do j=1,nspec800 call MPI_RECV(xmass1_tmp(:,j),num_trans,&801 & mp_sp,src_proc,mtag+(9+j),mp_comm_used,mp_status,mp_ierr)802 end do803 804 ! This is the maximum value numpart can possibly have after the transfer805 maxnumpart=numpart+num_trans806 807 ! Search for vacant space and copy from temporary storage808 !********************************************************809 minpart=1810 do i=1, num_trans811 ! Take into acount that we may have transferred invalid particles812 if (itra1_tmp(i).ne.itime) cycle813 do ipart=minpart,maxnumpart814 if (itra1(ipart).ne.itime) then815 itra1(ipart) = itra1_tmp(i)816 npoint(ipart) = npoint_tmp(i)817 nclass(ipart) = nclass_tmp(i)818 idt(ipart) = idt_tmp(i)819 itramem(ipart) = itramem_tmp(i)820 itrasplit(ipart) = itrasplit_tmp(i)821 xtra1(ipart) = xtra1_tmp(i)822 ytra1(ipart) = ytra1_tmp(i)823 ztra1(ipart) = ztra1_tmp(i)824 xmass1(ipart,:) = xmass1_tmp(i,:)825 goto 200 ! Storage space has been found, stop searching826 end if827 ! :TODO: add check for if too many particles requiried828 end do829 200 minpart=ipart+1830 end do831 ! Increase numpart, if necessary832 numpart=max(numpart,ipart)833 834 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&835 &itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)836 837 if (mp_dev_mode) then838 write(*,*) "numpart after: ", numpart839 write(*,FMT='(72("#"))')840 end if841 842 else843 ! This routine should only be called by the two participating processes844 write(*,*) "ERROR: wrong process has entered mpi_mod::mpif_redist_part"845 stop846 return847 end if848 849 ! Measure time for MPI communications850 !************************************851 if (mp_measure_time) call mpif_mtime('commtime',1)852 853 end subroutine mpif_redist_part854 855 856 subroutine mpif_tm_send_vars857 !***********************************************************************858 ! Distribute particle variables from pid0 to all processes.859 ! Called from timemanager860 ! *NOT IN USE* at the moment, but can be useful for debugging861 !862 !***********************************************************************863 use com_mod864 865 implicit none866 867 integer :: i868 869 ! Time for MPI communications870 !****************************871 if (mp_measure_time) call mpif_mtime('commtime',0)872 873 473 ! Distribute variables, send from pid 0 to other processes 874 474 !********************************************************************** … … 877 477 if (lroot) then 878 478 call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,& 879 & 479 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 880 480 if (mp_ierr /= 0) goto 600 881 481 call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,& 882 & 482 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 883 483 if (mp_ierr /= 0) goto 600 884 484 call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,& 885 & 485 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 886 486 if (mp_ierr /= 0) goto 600 887 487 call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,& 888 & 488 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 889 489 if (mp_ierr /= 0) goto 600 890 490 call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,& 891 & 491 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 892 492 if (mp_ierr /= 0) goto 600 893 493 894 494 ! int2 895 495 call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,MPI_IN_PLACE,& 896 & 496 &numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr) 897 497 if (mp_ierr /= 0) goto 600 898 498 899 499 ! real 900 500 call MPI_SCATTER(uap,numpart_mpi,mp_sp,MPI_IN_PLACE,& 901 & 501 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 902 502 if (mp_ierr /= 0) goto 600 903 503 call MPI_SCATTER(ucp,numpart_mpi,mp_sp,MPI_IN_PLACE,& 904 & 504 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 905 505 if (mp_ierr /= 0) goto 600 906 506 call MPI_SCATTER(uzp,numpart_mpi,mp_sp,MPI_IN_PLACE,& 907 & 507 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 908 508 if (mp_ierr /= 0) goto 600 909 509 910 510 call MPI_SCATTER(us,numpart_mpi,mp_sp,MPI_IN_PLACE,& 911 & 511 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 912 512 if (mp_ierr /= 0) goto 600 913 513 call MPI_SCATTER(vs,numpart_mpi,mp_sp,MPI_IN_PLACE,& 914 & 514 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 915 515 if (mp_ierr /= 0) goto 600 916 516 call MPI_SCATTER(ws,numpart_mpi,mp_sp,MPI_IN_PLACE,& 917 & 517 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 918 518 if (mp_ierr /= 0) goto 600 919 519 920 520 call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,MPI_IN_PLACE,& 921 & 521 &numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr) 922 522 if (mp_ierr /= 0) goto 600 923 523 call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,MPI_IN_PLACE,& 924 & 524 &numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr) 925 525 if (mp_ierr /= 0) goto 600 926 526 call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,MPI_IN_PLACE,& 927 & 527 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 928 528 if (mp_ierr /= 0) goto 600 929 529 930 530 do i=1,nspec 931 531 call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,MPI_IN_PLACE,& 932 & 532 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 933 533 if (mp_ierr /= 0) goto 600 934 534 end do … … 937 537 ! integers 938 538 call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,npoint,& 939 & 539 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 940 540 if (mp_ierr /= 0) goto 600 941 541 call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,idt,& 942 & 542 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 943 543 if (mp_ierr /= 0) goto 600 944 544 call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,itra1,& 945 & 545 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 946 546 if (mp_ierr /= 0) goto 600 947 547 call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,nclass,& 948 & 548 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 949 549 if (mp_ierr /= 0) goto 600 950 550 call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,itramem,& 951 & 551 &numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr) 952 552 if (mp_ierr /= 0) goto 600 953 553 954 554 ! int2 955 555 call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,cbt,& 956 & 556 &numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr) 957 557 if (mp_ierr /= 0) goto 600 958 558 959 559 ! reals 960 560 call MPI_SCATTER(uap,numpart_mpi,mp_sp,uap,& 961 & 561 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 962 562 if (mp_ierr /= 0) goto 600 963 563 call MPI_SCATTER(ucp,numpart_mpi,mp_sp,ucp,& 964 & 564 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 965 565 if (mp_ierr /= 0) goto 600 966 566 call MPI_SCATTER(uzp,numpart_mpi,mp_sp,uzp,& 967 & 567 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 968 568 if (mp_ierr /= 0) goto 600 969 569 970 570 call MPI_SCATTER(us,numpart_mpi,mp_sp,us,& 971 & 571 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 972 572 if (mp_ierr /= 0) goto 600 973 573 call MPI_SCATTER(vs,numpart_mpi,mp_sp,vs,& 974 & 574 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 975 575 if (mp_ierr /= 0) goto 600 976 576 call MPI_SCATTER(ws,numpart_mpi,mp_sp,ws,& 977 & 577 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 978 578 if (mp_ierr /= 0) goto 600 979 579 980 580 call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,xtra1,& 981 & 581 &numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr) 982 582 if (mp_ierr /= 0) goto 600 983 583 call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,ytra1,& 984 & 584 &numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr) 985 585 if (mp_ierr /= 0) goto 600 986 586 call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,ztra1,& 987 & 587 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 988 588 if (mp_ierr /= 0) goto 600 989 589 990 590 do i=1,nspec 991 591 call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,xmass1(:,i),& 992 & 592 &numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr) 993 593 if (mp_ierr /= 0) goto 600 994 594 end do … … 1545 1145 1546 1146 ! cloud water/ice: 1547 1548 ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)1549 ! if (mp_ierr /= 0) goto 6001550 1551 1552 1147 if (readclouds_nest(i)) then 1148 ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) 1149 ! if (mp_ierr /= 0) goto 600 1150 call MPI_Bcast(ctwcn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr) 1151 if (mp_ierr /= 0) goto 600 1152 end if 1553 1153 1554 1154 ! 2D fields … … 1798 1398 integer :: d2s1 = nxmax*nymax 1799 1399 integer :: d2s2 = nxmax*nymax*maxspec 1800 !integer :: d1_size1 = maxwf1400 !integer :: d1_size1 = maxwf 1801 1401 1802 1402 ! integer :: d3s1,d3s2,d2s1,d2s2 … … 2045 1645 if (dest.eq.id_read) cycle 2046 1646 do k=1, numbnests 2047 i=dest*nvar_async 2048 call MPI_Isend(uun(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2049 if (mp_ierr /= 0) goto 600 1647 i=dest*nvar_async 1648 call MPI_Isend(uun(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1649 if (mp_ierr /= 0) goto 600 1650 i=i+1 1651 call MPI_Isend(vvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1652 if (mp_ierr /= 0) goto 600 1653 i=i+1 1654 call MPI_Isend(wwn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1655 if (mp_ierr /= 0) goto 600 1656 i=i+1 1657 call MPI_Isend(ttn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1658 if (mp_ierr /= 0) goto 600 1659 i=i+1 1660 call MPI_Isend(rhon(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1661 if (mp_ierr /= 0) goto 600 1662 i=i+1 1663 call MPI_Isend(drhodzn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1664 if (mp_ierr /= 0) goto 600 1665 i=i+1 1666 call MPI_Isend(tthn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1667 if (mp_ierr /= 0) goto 600 1668 i=i+1 1669 call MPI_Isend(qvhn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1670 if (mp_ierr /= 0) goto 600 1671 i=i+1 1672 call MPI_Isend(qvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1673 if (mp_ierr /= 0) goto 600 1674 i=i+1 1675 call MPI_Isend(pvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1676 if (mp_ierr /= 0) goto 600 1677 i=i+1 1678 call MPI_Isend(cloudsn(:,:,:,mind,k),d3s1,MPI_INTEGER1,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1679 i=i+1 1680 if (mp_ierr /= 0) goto 600 1681 call MPI_Isend(cloudshn(:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1682 if (mp_ierr /= 0) goto 600 1683 i=i+1 1684 call MPI_Isend(vdepn(:,:,:,mind,k),d2s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1685 if (mp_ierr /= 0) goto 600 1686 i=i+1 1687 call MPI_Isend(psn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1688 if (mp_ierr /= 0) goto 600 1689 i=i+1 1690 call MPI_Isend(sdn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1691 if (mp_ierr /= 0) goto 600 1692 i=i+1 1693 ! 15 1694 call MPI_Isend(tccn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1695 if (mp_ierr /= 0) goto 600 1696 i=i+1 1697 call MPI_Isend(tt2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1698 if (mp_ierr /= 0) goto 600 1699 i=i+1 1700 call MPI_Isend(td2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1701 if (mp_ierr /= 0) goto 600 1702 i=i+1 1703 call MPI_Isend(lsprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1704 if (mp_ierr /= 0) goto 600 1705 i=i+1 1706 call MPI_Isend(convprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1707 if (mp_ierr /= 0) goto 600 1708 i=i+1 1709 call MPI_Isend(ustarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1710 if (mp_ierr /= 0) goto 600 1711 i=i+1 1712 call MPI_Isend(wstarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1713 if (mp_ierr /= 0) goto 600 1714 i=i+1 1715 call MPI_Isend(hmixn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1716 if (mp_ierr /= 0) goto 600 1717 i=i+1 1718 call MPI_Isend(tropopausen(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1719 if (mp_ierr /= 0) goto 600 1720 i=i+1 1721 call MPI_Isend(olin(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1722 if (mp_ierr /= 0) goto 600 1723 ! 25 1724 1725 ! Send cloud water if it exists. Increment counter always (as on receiving end) 1726 if (readclouds) then 2050 1727 i=i+1 2051 call MPI_Isend(vvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2052 if (mp_ierr /= 0) goto 600 2053 i=i+1 2054 call MPI_Isend(wwn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 1728 call MPI_Isend(ctwcn(:,:,mind,k),d2s1,mp_sp,dest,tm1,& 1729 &MPI_COMM_WORLD,reqs(i),mp_ierr) 2055 1730 if (mp_ierr /= 0) goto 600 2056 i=i+1 2057 call MPI_Isend(ttn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2058 if (mp_ierr /= 0) goto 600 2059 i=i+1 2060 call MPI_Isend(rhon(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2061 if (mp_ierr /= 0) goto 600 2062 i=i+1 2063 call MPI_Isend(drhodzn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2064 if (mp_ierr /= 0) goto 600 2065 i=i+1 2066 call MPI_Isend(tthn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2067 if (mp_ierr /= 0) goto 600 2068 i=i+1 2069 call MPI_Isend(qvhn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2070 if (mp_ierr /= 0) goto 600 2071 i=i+1 2072 call MPI_Isend(qvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2073 if (mp_ierr /= 0) goto 600 2074 i=i+1 2075 call MPI_Isend(pvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2076 if (mp_ierr /= 0) goto 600 2077 i=i+1 2078 call MPI_Isend(cloudsn(:,:,:,mind,k),d3s1,MPI_INTEGER1,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2079 i=i+1 2080 if (mp_ierr /= 0) goto 600 2081 call MPI_Isend(cloudshn(:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2082 if (mp_ierr /= 0) goto 600 2083 i=i+1 2084 call MPI_Isend(vdepn(:,:,:,mind,k),d2s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2085 if (mp_ierr /= 0) goto 600 2086 i=i+1 2087 call MPI_Isend(psn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2088 if (mp_ierr /= 0) goto 600 2089 i=i+1 2090 call MPI_Isend(sdn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2091 if (mp_ierr /= 0) goto 600 2092 i=i+1 2093 ! 15 2094 call MPI_Isend(tccn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2095 if (mp_ierr /= 0) goto 600 2096 i=i+1 2097 call MPI_Isend(tt2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2098 if (mp_ierr /= 0) goto 600 2099 i=i+1 2100 call MPI_Isend(td2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2101 if (mp_ierr /= 0) goto 600 2102 i=i+1 2103 call MPI_Isend(lsprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2104 if (mp_ierr /= 0) goto 600 2105 i=i+1 2106 call MPI_Isend(convprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2107 if (mp_ierr /= 0) goto 600 2108 i=i+1 2109 call MPI_Isend(ustarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2110 if (mp_ierr /= 0) goto 600 2111 i=i+1 2112 call MPI_Isend(wstarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2113 if (mp_ierr /= 0) goto 600 2114 i=i+1 2115 call MPI_Isend(hmixn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2116 if (mp_ierr /= 0) goto 600 2117 i=i+1 2118 call MPI_Isend(tropopausen(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2119 if (mp_ierr /= 0) goto 600 2120 i=i+1 2121 call MPI_Isend(olin(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr) 2122 if (mp_ierr /= 0) goto 600 2123 ! 25 2124 2125 ! Send cloud water if it exists. Increment counter always (as on receiving end) 2126 if (readclouds) then 2127 i=i+1 2128 call MPI_Isend(ctwcn(:,:,mind,k),d2s1,mp_sp,dest,tm1,& 2129 &MPI_COMM_WORLD,reqs(i),mp_ierr) 2130 if (mp_ierr /= 0) goto 600 2131 end if 2132 end do 1731 end if 2133 1732 end do 1733 end do 2134 1734 2135 1735 if (mp_measure_time) call mpif_mtime('commtime',1) … … 2210 1810 do k=1, numbnests 2211 1811 ! Get MPI tags/requests for communications 2212 2213 2214 2215 2216 2217 2218 2219 2220 2221 2222 2223 2224 2225 2226 2227 2228 2229 2230 2231 2232 2233 2234 2235 2236 2237 2238 2239 2240 2241 2242 2243 2244 2245 2246 2247 2248 2249 2250 2251 2252 2253 2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268 2269 2270 2271 2272 2273 2274 2275 2276 2277 2278 2279 2280 2281 2282 2283 2284 2285 2286 2287 2288 2289 2290 2291 2292 2293 2294 2295 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 2309 2310 1812 j=mp_pid*nvar_async 1813 call MPI_Irecv(uun(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1814 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1815 if (mp_ierr /= 0) goto 600 1816 j=j+1 1817 call MPI_Irecv(vvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1818 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1819 if (mp_ierr /= 0) goto 600 1820 j=j+1 1821 call MPI_Irecv(wwn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1822 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1823 if (mp_ierr /= 0) goto 600 1824 j=j+1 1825 call MPI_Irecv(ttn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1826 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1827 if (mp_ierr /= 0) goto 600 1828 j=j+1 1829 call MPI_Irecv(rhon(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1830 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1831 if (mp_ierr /= 0) goto 600 1832 j=j+1 1833 call MPI_Irecv(drhodzn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1834 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1835 if (mp_ierr /= 0) goto 600 1836 j=j+1 1837 call MPI_Irecv(tthn(:,:,:,mind,k),d3s2,mp_sp,id_read,MPI_ANY_TAG,& 1838 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1839 if (mp_ierr /= 0) goto 600 1840 j=j+1 1841 call MPI_Irecv(qvhn(:,:,:,mind,k),d3s2,mp_sp,id_read,MPI_ANY_TAG,& 1842 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1843 if (mp_ierr /= 0) goto 600 1844 j=j+1 1845 call MPI_Irecv(qvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1846 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1847 if (mp_ierr /= 0) goto 600 1848 j=j+1 1849 call MPI_Irecv(pvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,& 1850 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1851 if (mp_ierr /= 0) goto 600 1852 j=j+1 1853 call MPI_Irecv(cloudsn(:,:,:,mind,k),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,& 1854 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1855 if (mp_ierr /= 0) goto 600 1856 j=j+1 1857 call MPI_Irecv(cloudshn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1858 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1859 if (mp_ierr /= 0) goto 600 1860 j=j+1 1861 call MPI_Irecv(vdepn(:,:,:,mind,k),d2s2,mp_sp,id_read,MPI_ANY_TAG,& 1862 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1863 if (mp_ierr /= 0) goto 600 1864 j=j+1 1865 call MPI_Irecv(psn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1866 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1867 if (mp_ierr /= 0) goto 600 1868 j=j+1 1869 call MPI_Irecv(sdn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1870 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1871 if (mp_ierr /= 0) goto 600 1872 j=j+1 1873 call MPI_Irecv(tccn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1874 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1875 if (mp_ierr /= 0) goto 600 1876 j=j+1 1877 call MPI_Irecv(tt2n(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1878 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1879 if (mp_ierr /= 0) goto 600 1880 j=j+1 1881 call MPI_Irecv(td2n(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1882 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1883 if (mp_ierr /= 0) goto 600 1884 j=j+1 1885 call MPI_Irecv(lsprecn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1886 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1887 if (mp_ierr /= 0) goto 600 1888 j=j+1 1889 call MPI_Irecv(convprecn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1890 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1891 if (mp_ierr /= 0) goto 600 1892 call MPI_Irecv(ustarn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1893 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1894 if (mp_ierr /= 0) goto 600 1895 j=j+1 1896 call MPI_Irecv(wstarn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1897 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1898 if (mp_ierr /= 0) goto 600 1899 j=j+1 1900 call MPI_Irecv(hmixn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1901 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1902 if (mp_ierr /= 0) goto 600 1903 j=j+1 1904 call MPI_Irecv(tropopausen(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1905 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1906 if (mp_ierr /= 0) goto 600 1907 j=j+1 1908 call MPI_Irecv(olin(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,& 1909 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1910 if (mp_ierr /= 0) goto 600 2311 1911 2312 1912 ! Post request for clwc. These data are possibly not sent, request must then be cancelled 2313 1913 ! For now assume that data at all steps either have or do not have water 2314 2315 2316 2317 2318 2319 2320 1914 if (readclouds) then 1915 j=j+1 1916 call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,& 1917 &MPI_COMM_WORLD,reqs(j),mp_ierr) 1918 if (mp_ierr /= 0) goto 600 1919 end if 1920 end do 2321 1921 2322 1922 if (mp_measure_time) call mpif_mtime('commtime',1) … … 2340 1940 ! 2341 1941 ! TODO 2342 ! NB: take into account nested wind fields by using a separate 2343 ! subroutine mpif_gf_request_nest (and separate request arrays for the 2344 ! nested variables) 1942 ! take into account nested wind fields 2345 1943 ! 2346 1944 !******************************************************************************* … … 2751 2349 write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR CONCCALC:',& 2752 2350 & mp_conccalc_time_total 2753 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR VERTTRANSFORM:',& 2754 ! & mp_vt_wtime_total 2755 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',& 2756 ! & mp_vt_time_total 2757 ! NB: the 'flush' function is possibly a gfortran-specific extension, 2758 ! comment out if it gives problems 2759 ! call flush() 2351 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR VERTTRANSFORM:',& 2352 ! & mp_vt_wtime_total 2353 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',& 2354 ! & mp_vt_time_total 2355 ! NB: the 'flush' function is possibly a gfortran-specific extension 2356 call flush() 2760 2357 end if 2761 2358 end do … … 2791 2388 2792 2389 2793 end subroutine mpif_finalize2794 2795 2796 subroutine get_lun(my_lun)2390 end subroutine mpif_finalize 2391 2392 2393 subroutine get_lun(my_lun) 2797 2394 !*********************************************************************** 2798 2395 ! get_lun: … … 2800 2397 !*********************************************************************** 2801 2398 2802 implicit none2803 2804 integer, intent(inout) :: my_lun2805 integer, save :: free_lun=1002806 logical :: exists, iopen2807 2808 !*********************************************************************** 2809 2810 loop1: do2811 inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen)2812 if (exists .and. .not.iopen) exit loop12813 free_lun = free_lun+12814 end do loop12815 my_lun = free_lun2816 2817 end subroutine get_lun2818 2819 2820 subroutine write_data_dbg(array_in, array_name, tstep, ident)2399 implicit none 2400 2401 integer, intent(inout) :: my_lun 2402 integer, save :: free_lun=100 2403 logical :: exists, iopen 2404 2405 !*********************************************************************** 2406 2407 loop1: do 2408 inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen) 2409 if (exists .and. .not.iopen) exit loop1 2410 free_lun = free_lun+1 2411 end do loop1 2412 my_lun = free_lun 2413 2414 end subroutine get_lun 2415 2416 2417 subroutine write_data_dbg(array_in, array_name, tstep, ident) 2821 2418 !*********************************************************************** 2822 2419 ! Write one-dimensional arrays to file (for debugging purposes) 2823 2420 !*********************************************************************** 2824 implicit none2825 2826 real, intent(in), dimension(:) :: array_in2827 integer, intent(in) :: tstep2828 integer :: lios2829 character(LEN=*), intent(in) :: ident, array_name2830 2831 character(LEN=8) :: c_ts2832 character(LEN=40) :: fn_1, fn_22833 2834 !*********************************************************************** 2835 2836 write(c_ts, FMT='(I8.8,BZ)') tstep2837 fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)2838 write(c_ts, FMT='(I2.2,BZ)') mp_np2839 fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat'2840 2841 call get_lun(dat_lun)2842 open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', &2843 FORM='UNFORMATTED', STATUS='REPLACE')2844 write(UNIT=dat_lun, IOSTAT=lios) array_in2845 close(UNIT=dat_lun)2846 2847 end subroutine write_data_dbg2421 implicit none 2422 2423 real, intent(in), dimension(:) :: array_in 2424 integer, intent(in) :: tstep 2425 integer :: lios 2426 character(LEN=*), intent(in) :: ident, array_name 2427 2428 character(LEN=8) :: c_ts 2429 character(LEN=40) :: fn_1, fn_2 2430 2431 !*********************************************************************** 2432 2433 write(c_ts, FMT='(I8.8,BZ)') tstep 2434 fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident) 2435 write(c_ts, FMT='(I2.2,BZ)') mp_np 2436 fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat' 2437 2438 call get_lun(dat_lun) 2439 open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', & 2440 FORM='UNFORMATTED', STATUS='REPLACE') 2441 write(UNIT=dat_lun, IOSTAT=lios) array_in 2442 close(UNIT=dat_lun) 2443 2444 end subroutine write_data_dbg 2848 2445 2849 2446 … … 2902 2499 clwc(:,:,:,li:ui)=0.0 2903 2500 ciwc(:,:,:,li:ui)=0.0 2904 2501 2905 2502 ! 2D fields 2906 2503 … … 2919 2516 tropopause(:,:,:,li:ui)=10000. 2920 2517 oli(:,:,:,li:ui)=0.01 2921 2518 2922 2519 end subroutine set_fields_synthetic 2923 2520 -
src/netcdf_output_mod.f90
r4c64400 rf28aa0a 19 19 !********************************************************************** 20 20 21 22 21 !***************************************************************************** 23 22 ! * … … 25 24 ! residence time and wet and dry deposition output. * 26 25 ! * 27 ! - writeheader_netcdf generates file including all information previously *26 ! - writeheader_netcdf generates file including all information previously * 28 27 ! stored in separate header files * 29 ! - concoutput_netcdf write concentration output and wet and dry deposition *28 ! - concoutput_netcdf write concentration output and wet and dry deposition * 30 29 ! * 31 30 ! Author: D. Brunner * … … 894 893 gridsigmatotal=0. 895 894 gridtotalunc=0. 896 wetgridtotal=0. _dep_prec897 wetgridsigmatotal=0. _dep_prec898 wetgridtotalunc=0. _dep_prec899 drygridtotal=0. _dep_prec900 drygridsigmatotal=0. _dep_prec901 drygridtotalunc=0. _dep_prec895 wetgridtotal=0. 896 wetgridsigmatotal=0. 897 wetgridtotalunc=0. 898 drygridtotal=0. 899 drygridsigmatotal=0. 900 drygridtotalunc=0. 902 901 903 902 do ks=1,nspec … … 923 922 wetgridsigma(ix,jy),nclassunc) 924 923 ! Multiply by number of classes to get total concentration 925 wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc,kind= sp)924 wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc,kind=dep_prec) 926 925 wetgridtotal=wetgridtotal+wetgrid(ix,jy) 927 926 ! Calculate standard deviation of the mean … … 947 946 drygridsigma(ix,jy),nclassunc) 948 947 ! Multiply by number of classes to get total concentration 949 drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc ,kind=sp)948 drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc) 950 949 drygridtotal=drygridtotal+drygrid(ix,jy) 951 950 ! Calculate standard deviation of the mean 952 951 drygridsigma(ix,jy)= & 953 952 drygridsigma(ix,jy)* & 954 sqrt(real(nclassunc , kind=dep_prec))953 sqrt(real(nclassunc)) 955 954 drygridsigmatotal=drygridsigmatotal+ & 956 955 drygridsigma(ix,jy) … … 1055 1054 if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ & 1056 1055 wetgridtotal 1057 if (drygridtotal.gt.0.) drygridtotalunc= real(drygridsigmatotal/ &1058 drygridtotal , kind=dep_prec)1056 if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ & 1057 drygridtotal 1059 1058 1060 1059 ! Dump of receptor concentrations … … 1299 1298 wetgridsigma(ix,jy)= & 1300 1299 wetgridsigma(ix,jy)* & 1301 sqrt(real(nclassunc ,kind=dep_prec))1300 sqrt(real(nclassunc)) 1302 1301 endif 1303 1302 … … 1320 1319 drygridsigma(ix,jy)= & 1321 1320 drygridsigma(ix,jy)* & 1322 sqrt(real(nclassunc ,kind=dep_prec))1321 sqrt(real(nclassunc)) 1323 1322 endif 1324 1323 -
src/outg_mod.f90
r4c64400 r6a678e3 22 22 module outg_mod 23 23 24 use par_mod, only: dep_prec , sp24 use par_mod, only: dep_prec 25 25 26 26 implicit none … … 39 39 real,allocatable, dimension (:,:,:) :: factor3d 40 40 real,allocatable, dimension (:,:,:) :: grid 41 real( sp),allocatable, dimension (:,:) :: wetgrid42 real( sp),allocatable, dimension (:,:) :: drygrid41 real(dep_prec),allocatable, dimension (:,:) :: wetgrid 42 real(dep_prec),allocatable, dimension (:,:) :: drygrid 43 43 real,allocatable, dimension (:,:,:) :: gridsigma 44 44 real(dep_prec),allocatable, dimension (:,:) :: drygridsigma -
src/par_mod.f90
rd1a8707 rd1a8707 58 58 59 59 integer,parameter :: dep_prec=sp 60 61 !****************************************************************62 ! Set to T to disable use of kernel for concentrations/deposition63 !****************************************************************64 65 logical, parameter :: lnokernel=.false.66 60 67 61 !*********************************************************** … … 192 186 !************************************************** 193 187 194 integer,parameter :: maxpart=10000000 195 integer,parameter :: maxspec=4 196 188 integer,parameter :: maxpart=20000000 189 integer,parameter :: maxspec=2 197 190 real,parameter :: minmass=0.0001 198 191 -
src/readavailable.f90
r4c64400 rdb712a8 236 236 do i=2,numbwf 237 237 idiff=abs(wftime(i)-wftime(i-1)) 238 if (idiff.gt.idiffmax .and.lroot) then238 if (idiff.gt.idiffmax) then 239 239 write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO' 240 240 write(*,*) 'WIND FIELDS IS TOO BIG FOR TRANSPORT CALCULATION.& 241 241 &' 242 242 write(*,*) 'THEREFORE, TRAJECTORIES HAVE TO BE SKIPPED.' 243 else if (idiff.gt.idiffnorm .and.lroot) then243 else if (idiff.gt.idiffnorm) then 244 244 write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO' 245 245 write(*,*) 'WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION' -
src/readcommand.f90
r8ee24a5 r8ee24a5 170 170 if (index(line,'LDIRECT') .eq. 0) then 171 171 old = .false. 172 if (lroot) write(*,*) 'COMMAND in old short format, & 173 &please update to namelist format' 172 write(*,*) 'COMMAND in old short format, please update to namelist format' 174 173 else 175 174 old = .true. 176 if (lroot) write(*,*) 'COMMAND in old long format, & 177 &please update to namelist format' 175 write(*,*) 'COMMAND in old long format, please update to namelist format' 178 176 endif 179 177 rewind(unitcommand) -
src/readreleases.f90
rc8fc724 r8ee24a5 218 218 rewind(unitreleases) 219 219 220 if (nspec.gt.maxspec) goto 994221 222 220 ! allocate arrays of matching size for number of species (namelist output) 223 221 deallocate(mass) … … 276 274 do i=1,maxspec 277 275 DRYDEPSPEC(i)=.false. 278 WETDEPSPEC(i)=.false.279 276 end do 280 277 … … 370 367 &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.))) then 371 368 WETDEP=.true. 372 WETDEPSPEC(i)=.true.373 369 if (lroot) then 374 370 write (*,*) ' Below-cloud scavenging: ON' … … 382 378 if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.)) then 383 379 WETDEP=.true. 384 WETDEPSPEC(i)=.true.385 380 if (lroot) then 386 381 write (*,*) ' In-cloud scavenging: ON' … … 402 397 endif 403 398 404 end do ! end loop over species399 end do 405 400 406 401 if (WETDEP.or.DRYDEP) DEP=.true. … … 593 588 endif 594 589 595 596 590 ! Determine the release rate (particles per second) and total number 597 591 ! of particles released during the simulation … … 653 647 endif 654 648 655 if (lroot) then656 write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))657 end if658 659 649 return 660 650 -
src/releaseparticles_mpi.f90
r16b61a5 r7e52e2e 21 21 22 22 subroutine releaseparticles(itime) 23 ! o24 !*****************************************************************************25 ! *26 ! This subroutine releases particles from the release locations. *27 ! *28 ! It searches for a "vacant" storage space and assigns all particle *29 ! information to that space. A space is vacant either when no particle *30 ! is yet assigned to it, or when it's particle is expired and, thus, *31 ! the storage space is made available to a new particle. *32 ! *33 ! Author: A. Stohl *34 ! *35 ! 29 June 2002 *36 ! *37 ! CHANGES *38 ! 12/2014 eso: MPI version *39 ! *40 !*****************************************************************************41 ! *42 ! Variables: *43 ! itime [s] current time *44 ! ireleasestart, ireleaseend start and end times of all releases *45 ! npart(maxpoint) number of particles to be released in total *46 ! numrel number of particles to be released during this time *47 ! step *48 ! addoneextra particle assigned to MPI process if *49 ! mod(npart(i),mp_partgroup_np) .ne. 0) *50 !*****************************************************************************23 ! o 24 !***************************************************************************** 25 ! * 26 ! This subroutine releases particles from the release locations. * 27 ! * 28 ! It searches for a "vacant" storage space and assigns all particle * 29 ! information to that space. A space is vacant either when no particle * 30 ! is yet assigned to it, or when it's particle is expired and, thus, * 31 ! the storage space is made available to a new particle. * 32 ! * 33 ! Author: A. Stohl * 34 ! * 35 ! 29 June 2002 * 36 ! * 37 ! CHANGES * 38 ! 12/2014 eso: MPI version * 39 ! * 40 !***************************************************************************** 41 ! * 42 ! Variables: * 43 ! itime [s] current time * 44 ! ireleasestart, ireleaseend start and end times of all releases * 45 ! npart(maxpoint) number of particles to be released in total * 46 ! numrel number of particles to be released during this time * 47 ! step * 48 ! addpart extra particle assigned to MPI process if * 49 ! mod(npart(i),mp_partgroup_np) .ne. 0) * 50 !***************************************************************************** 51 51 52 52 use point_mod … … 59 59 implicit none 60 60 61 !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)61 !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint) 62 62 real :: xaux,yaux,zaux,rfraction 63 63 real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4 … … 73 73 74 74 integer :: idummy = -7 75 !save idummy,xmasssave76 !data idummy/-7/,xmasssave/maxpoint*0./75 !save idummy,xmasssave 76 !data idummy/-7/,xmasssave/maxpoint*0./ 77 77 78 78 logical :: first_call=.true. 79 79 80 ! Use different seed for each process.81 !****************************************************************************80 ! Use different seed for each process. 81 !**************************************************************************** 82 82 if (first_call) then 83 83 idummy=idummy+mp_seed … … 87 87 mind2=memind(2) 88 88 89 ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time)90 !*****************************************************************************89 ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time) 90 !***************************************************************************** 91 91 92 92 julmonday=juldate(19000101,0) ! this is a Monday … … 97 97 98 98 99 ! For every release point, check whether we are in the release time interval100 !***************************************************************************99 ! For every release point, check whether we are in the release time interval 100 !*************************************************************************** 101 101 102 102 minpart=1 … … 105 105 (itime.le.ireleaseend(i))) then 106 106 107 ! Determine the local day and time108 !*********************************107 ! Determine the local day and time 108 !********************************* 109 109 110 110 xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx ! longitude needed to determine local time … … 124 124 endif 125 125 126 ! Calculate a species- and time-dependent correction factor, distinguishing between127 ! area (those with release starting at surface) and point (release starting above surface) sources128 ! Also, calculate an average time correction factor (species independent)129 !*****************************************************************************126 ! Calculate a species- and time-dependent correction factor, distinguishing between 127 ! area (those with release starting at surface) and point (release starting above surface) sources 128 ! Also, calculate an average time correction factor (species independent) 129 !***************************************************************************** 130 130 average_timecorrect=0. 131 131 do k=1,nspec … … 139 139 average_timecorrect=average_timecorrect/real(nspec) 140 140 141 ! Determine number of particles to be released this time; at start and at end of release,142 ! only half the particles are released143 !*****************************************************************************141 ! Determine number of particles to be released this time; at start and at end of release, 142 ! only half the particles are released 143 !***************************************************************************** 144 144 145 145 if (ireleasestart(i).ne.ireleaseend(i)) then … … 149 149 (itime.eq.ireleaseend(i))) rfraction=rfraction/2. 150 150 151 ! Take the species-average time correction factor in order to scale the152 ! number of particles released this time153 ! Also scale by number of MPI processes154 !**********************************************************************151 ! Take the species-average time correction factor in order to scale the 152 ! number of particles released this time 153 ! Also scale by number of MPI processes 154 !********************************************************************** 155 155 156 156 rfraction=rfraction*average_timecorrect … … 158 158 rfraction=rfraction+xmasssave(i) ! number to be released at this time 159 159 160 ! number to be released for one process160 ! number to be released for one process 161 161 if (mp_partid.lt.mod(int(rfraction),mp_partgroup_np)) then 162 162 addone=1 … … 180 180 numrel=npart(i)/mp_partgroup_np+addone 181 181 endif 182 182 183 183 xaux=xpoint2(i)-xpoint1(i) 184 184 yaux=ypoint2(i)-ypoint1(i) … … 187 187 do ipart=minpart,maxpart_mpi ! search for free storage space 188 188 189 ! If a free storage space is found, attribute everything to this array element190 !*****************************************************************************189 ! If a free storage space is found, attribute everything to this array element 190 !***************************************************************************** 191 191 192 192 if (itra1(ipart).ne.itime) then 193 193 194 ! Particle coordinates are determined by using a random position within the release volume195 !*****************************************************************************196 197 ! Determine horizontal particle position198 !***************************************194 ! Particle coordinates are determined by using a random position within the release volume 195 !***************************************************************************** 196 197 ! Determine horizontal particle position 198 !*************************************** 199 199 200 200 xtra1(ipart)=xpoint1(i)+ran1(idummy)*xaux … … 207 207 ytra1(ipart)=ypoint1(i)+ran1(idummy)*yaux 208 208 209 ! Assign mass to particle: Total mass divided by total number of particles.210 ! Time variation has partly been taken into account already by a species-average211 ! correction factor, by which the number of particles released this time has been212 ! scaled. Adjust the mass per particle by the species-dependent time correction factor213 ! divided by the species-average one214 !*****************************************************************************209 ! Assign mass to particle: Total mass divided by total number of particles. 210 ! Time variation has partly been taken into account already by a species-average 211 ! correction factor, by which the number of particles released this time has been 212 ! scaled. Adjust the mass per particle by the species-dependent time correction factor 213 ! divided by the species-average one 214 !***************************************************************************** 215 215 do k=1,nspec 216 xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &217 *timecorrect(k)/average_timecorrect218 ! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k219 ! Assign certain properties to particle220 !**************************************216 xmass1(ipart,k)=xmass(i,k)/real(npart(i)) & 217 *timecorrect(k)/average_timecorrect 218 ! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k 219 ! Assign certain properties to particle 220 !************************************** 221 221 end do 222 222 nclass(ipart)=min(int(ran1(idummy)*real(nclassunc))+1, & … … 234 234 235 235 236 ! Determine vertical particle position237 !*************************************236 ! Determine vertical particle position 237 !************************************* 238 238 239 239 ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux 240 240 241 ! Interpolation of topography and density242 !****************************************243 244 ! Determine the nest we are in245 !*****************************241 ! Interpolation of topography and density 242 !**************************************** 243 244 ! Determine the nest we are in 245 !***************************** 246 246 247 247 ngrid=0 … … 257 257 43 continue 258 258 259 ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation260 !*****************************************************************************259 ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation 260 !***************************************************************************** 261 261 262 262 if (ngrid.gt.0) then … … 294 294 endif 295 295 296 ! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters297 !*****************************************************************************296 ! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters 297 !***************************************************************************** 298 298 if (kindz(i).eq.3) then 299 299 presspart=ztra1(ipart) … … 337 337 endif 338 338 339 ! If release positions are given in meters above sea level, subtract the340 ! topography from the starting height341 !***********************************************************************339 ! If release positions are given in meters above sea level, subtract the 340 ! topography from the starting height 341 !*********************************************************************** 342 342 343 343 if (kindz(i).eq.2) ztra1(ipart)=ztra1(ipart)-topo … … 348 348 349 349 350 ! For special simulations, multiply particle concentration air density;351 ! Simply take the 2nd field in memory to do this (accurate enough)352 !***********************************************************************353 !AF IND_SOURCE switches between different units for concentrations at the source354 !Af NOTE that in backward simulations the release of particles takes place at the355 !Af receptor and the sampling at the source.356 !Af 1="mass"357 !Af 2="mass mixing ratio"358 !Af IND_RECEPTOR switches between different units for concentrations at the receptor359 !Af 1="mass"360 !Af 2="mass mixing ratio"361 362 !Af switches for the releasefile:363 !Af IND_REL = 1 : xmass * rho364 !Af IND_REL = 0 : xmass * 1365 366 !Af ind_rel is defined in readcommand.f350 ! For special simulations, multiply particle concentration air density; 351 ! Simply take the 2nd field in memory to do this (accurate enough) 352 !*********************************************************************** 353 !AF IND_SOURCE switches between different units for concentrations at the source 354 !Af NOTE that in backward simulations the release of particles takes place at the 355 !Af receptor and the sampling at the source. 356 !Af 1="mass" 357 !Af 2="mass mixing ratio" 358 !Af IND_RECEPTOR switches between different units for concentrations at the receptor 359 !Af 1="mass" 360 !Af 2="mass mixing ratio" 361 362 !Af switches for the releasefile: 363 !Af IND_REL = 1 : xmass * rho 364 !Af IND_REL = 0 : xmass * 1 365 366 !Af ind_rel is defined in readcommand.f 367 367 368 368 if (ind_rel .eq. 1) then 369 369 370 ! Interpolate the air density371 !****************************370 ! Interpolate the air density 371 !**************************** 372 372 373 373 do ii=2,nz … … 403 403 404 404 405 ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density406 !********************************************************************405 ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density 406 !******************************************************************** 407 407 408 408 do k=1,nspec … … 416 416 endif 417 417 end do 418 if (ipart.gt.maxpart _mpi) goto 996418 if (ipart.gt.maxpart) goto 996 419 419 420 420 34 minpart=ipart+1 421 421 end do 422 endif422 endif 423 423 end do 424 424 … … 426 426 return 427 427 428 996 continue428 996 continue 429 429 write(*,*) '#####################################################' 430 430 write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####' -
src/timemanager_mpi.f90
r16b61a5 r3b80e98 104 104 105 105 logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete 106 integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 106 integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 !,mind 107 107 ! integer :: ksp 108 108 integer :: ip … … 155 155 156 156 do itime=0,ideltas,lsynctime 157 158 157 159 158 ! Computation of wet deposition, OH reaction and mass transfer … … 167 166 !******************************************************************** 168 167 169 if (mp_d bg_mode) write(*,*) 'myid, itime: ',mp_pid,itime168 if (mp_dev_mode) write(*,*) 'myid, itime: ',mp_pid,itime 170 169 171 170 if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then … … 331 330 call releaseparticles(itime) 332 331 endif 333 334 335 ! Check if particles should be redistributed among processes336 !***********************************************************337 call mpif_calculate_part_redist(itime)338 332 339 333 … … 554 548 ! Decide whether to write an estimate of the number of particles released, 555 549 ! or exact number (require MPI reduce operation) 556 if (mp_d bg_mode) then550 if (mp_dev_mode) then 557 551 numpart_tot_mpi = numpart 558 552 else … … 561 555 562 556 if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader).and.& 563 &.not.mp_d bg_mode) then557 &.not.mp_dev_mode) then 564 558 call MPI_Reduce(numpart, numpart_tot_mpi, 1, MPI_INTEGER, MPI_SUM, id_root, & 565 559 & mp_comm_used, mp_ierr) … … 567 561 568 562 !CGZ-lifetime: output species lifetime 569 if (lroot.or.mp_d bg_mode) then563 if (lroot.or.mp_dev_mode) then 570 564 ! write(*,*) 'Overview species lifetime in days', & 571 565 ! real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0)) … … 577 571 ! end if 578 572 end if 579 580 ! Write number of particles for all processes581 if (mp_dev_mode) write(*,*) "PID, itime, numpart", mp_pid,itime,numpart582 583 573 584 574 45 format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES: Uncertainty: ',3f7.3) … … 870 860 endif 871 861 deallocate(gridunc) 872 deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass) 873 if (allocated(checklifetime)) deallocate(checklifetime) 862 deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass, checklifetime) 874 863 deallocate(ireleasestart,ireleaseend,npart,kindz) 875 864 deallocate(xmasssave) -
src/unc_mod.f90
r4c64400 r6a678e3 26 26 implicit none 27 27 28 real,allocatable ,dimension (:,:,:,:,:,:,:) :: gridunc28 real,allocatable ,dimension (:,:,:,:,:,:,:) :: gridunc 29 29 real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn 30 30 real(dep_prec),allocatable, dimension (:,:,:,:,:,:) :: drygridunc -
src/wetdepo.f90
rc8fc724 r0539b8f 66 66 integer :: itage,nage,il,interp_time, n 67 67 integer :: ks, kp 68 integer (selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count68 integer :: blc_count, inc_count 69 69 real :: grfraction(3),wetscav 70 70 real :: wetdeposit(maxspec),restmass … … 83 83 !************************ 84 84 85 blc_count (:)=086 inc_count (:)=085 blc_count=0 86 inc_count=0 87 87 88 88 do jpart=1,numpart … … 104 104 105 105 do ks=1,nspec ! loop over species 106 107 if (WETDEPSPEC(ks).eqv..false.) cycle108 106 109 107 !************************************************** … … 166 164 167 165 ! count the total number of below-cloud and in-cloud occurences: 168 tot_blc_count (1:nspec)=tot_blc_count(1:nspec)+blc_count(1:nspec)169 tot_inc_count (1:nspec)=tot_inc_count(1:nspec)+inc_count(1:nspec)166 tot_blc_count=tot_blc_count+blc_count 167 tot_inc_count=tot_inc_count+inc_count 170 168 171 169 end subroutine wetdepo -
src/wetdepokernel.f90
r1c0d5e6 r1c0d5e6 40 40 ! * 41 41 !***************************************************************************** 42 ! Changes:43 ! eso 10/2016: Added option to disregard kernel44 !45 !*****************************************************************************46 42 47 43 use unc_mod … … 77 73 endif 78 74 79 ! If no kernel is used, direct attribution to grid cell80 !******************************************************81 75 82 if (lnokernel) then83 do ks=1,nspec84 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &85 (jy.le.numygrid-1)) then86 wetgridunc(ix,jy,ks,kp,nunc,nage)= &87 wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)88 end if89 end do90 else ! use kernel91 92 76 ! Determine mass fractions for four grid points 93 77 !********************************************** … … 128 112 endif 129 113 end do 130 end if131 114 132 115 end subroutine wetdepokernel
Note: See TracChangeset
for help on using the changeset viewer.