Changes in / [d9f0585:d404d98] in flexpart.git
- Location:
- src
- Files:
-
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART.f90
rc8fc724 rdb712a8 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 (2016-11-02)'71 flexversion='Version '//trim(flexversion_major)//'.0beta (2015-05-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 r0f7835d 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_mpi 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_mpi) & 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_mpi 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_mpi) & 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
rc8fc724 r62e65c7 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 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),& 580 & OHREA,ASSSPEC 578 logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC 581 579 582 580 ! numxgrid,numygrid number of grid points in x,y-direction … … 595 593 ! DRYDEPSPEC .true., if dry deposition is switched on for that species 596 594 ! WETDEP .true., if wet deposition is switched on 597 ! WETDEPSPEC .true., if wet deposition is switched on for that species598 595 ! OHREA .true., if OH reaction is switched on 599 596 ! ASSSPEC .true., if there are two species asscoiated … … 746 743 747 744 ! These variables are used to avoid having separate versions of 748 ! files in cases where differences with MPI version areminor (eso)745 ! files in cases where differences with MPI version is minor (eso) 749 746 !***************************************************************** 750 747 integer :: mpi_mode=0 ! .gt. 0 if running MPI version -
src/conccalc.f90
r4c64400 rfdc0f03 126 126 !***************************************************************************** 127 127 do ind=indz,indzp 128 rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind, memind(2)) &128 rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) & 129 129 +p2*rho(ixp,jy ,ind,2) & 130 130 +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 185 (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell -
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 rdb712a8 626 626 !************************* 627 627 628 ! do ks=1,nspec629 ! do kp=1,maxpointspec_act630 ! do i=1,numreceptor631 ! creceptor(i,ks)=0.632 ! end do633 ! do jy=0,numygrid-1634 ! do ix=0,numxgrid-1635 ! do l=1,nclassunc636 ! do nage=1,nageclass637 ! do kz=1,numzgrid638 ! gridunc(ix,jy,kz,ks,kp,l,nage)=0.639 ! end do640 ! end do641 ! end do642 ! end do643 ! end do644 ! end do645 ! end do628 ! do ks=1,nspec 629 ! do kp=1,maxpointspec_act 630 ! do i=1,numreceptor 631 ! creceptor(i,ks)=0. 632 ! end do 633 ! do jy=0,numygrid-1 634 ! do ix=0,numxgrid-1 635 ! do l=1,nclassunc 636 ! do nage=1,nageclass 637 ! do kz=1,numzgrid 638 ! gridunc(ix,jy,kz,ks,kp,l,nage)=0. 639 ! end do 640 ! end do 641 ! end do 642 ! end do 643 ! end do 644 ! end do 645 ! end do 646 646 creceptor(:,:)=0. 647 647 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
r4c64400 re200b7a 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 ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & 105 84 (jy.le.numygrid-1)) then 106 107 drygridunc(ix,jy,ks,kp,nunc,nage)= &108 drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w109 continue110 85 w=wx*wy 86 drygridunc(ix,jy,ks,kp,nunc,nage)= & 87 drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w 88 continue 89 endif 111 90 112 91 if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. & … … 133 112 endif 134 113 135 end do 136 end if 114 end do 137 115 138 116 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/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 r0f7835d 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 354 355 356 ! Check whether numpart is really smaller than maxpart 357 ! *****************************************************358 359 ! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier 360 if (numpart.gt.maxpart) then361 362 write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart363 364 365 366 304 end do 305 306 307 ! Check whether numpart is really smaller than maxpart per process 308 !***************************************************************** 309 310 if (numpart.gt.maxpart_mpi) then 311 write(*,*) 'numpart too large: change source in init_atm_mass.f' 312 write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart_mpi 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
r4c64400 r341f4b7 419 419 par_mod.o point_mod.o unc_mod.o xmass_mod.o 420 420 timemanager_mpi.o: com_mod.o flux_mod.o mpi_mod.o oh_mod.o outg_mod.o \ 421 par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o421 par_mod.o point_mod.o unc_mod.o xmass_mod.o 422 422 unc_mod.o: par_mod.o 423 423 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 r861805a 122 122 logical, parameter :: mp_dev_mode = .false. 123 123 logical, parameter :: mp_dbg_out = .false. 124 logical, parameter :: mp_time_barrier=. true.124 logical, parameter :: mp_time_barrier=.false. 125 125 logical, parameter :: mp_measure_time=.false. 126 126 logical, parameter :: mp_exact_numpart=.true. … … 251 251 write(*,FMT='(80("#"))') 252 252 end if 253 lmp_sync=.true. 253 lmp_sync=.true. ! :DBG: eso fix this... 254 254 end if 255 255 … … 261 261 !********************************************************************** 262 262 263 ! id_read = min(mp_np-1, 1)264 263 id_read = mp_np-1 265 264 … … 331 330 ! Set maxpart per process 332 331 ! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution 333 maxpart_mpi=int(mp_maxpart_factor* real(maxpart)/real(mp_partgroup_np))332 maxpart_mpi=int(mp_maxpart_factor*maxpart/mp_partgroup_np) 334 333 335 334 ! Do not allocate particle data arrays for readwind process … … 487 486 integer :: i,jj, addone 488 487 489 ! Exit if too many particles490 if (num_part.gt.maxpart_mpi) then491 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 stop499 end if500 501 502 488 ! Time for MPI communications 503 489 !**************************** … … 541 527 542 528 if (mp_measure_time) call mpif_mtime('commtime',1) 529 write(*,*) "PID ", mp_partid, "ending MPI_Scatter operation" 543 530 544 531 goto 601 … … 548 535 549 536 ! After the transfer of particles it it possible that the value of 550 ! "numpart" is larger than the actual used, so we reduce it if there are537 ! "numpart" is larger than the actual, so we reduce it if there are 551 538 ! invalid particles at the end of the arrays 552 539 601 do i=num_part, 1, -1 … … 641 628 if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.& 642 629 & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then 643 ! 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 630 call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2) 648 631 end if … … 678 661 integer :: i, j, minpart, ipart, maxnumpart 679 662 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 663 ! Measure time for MPI communications 689 664 !************************************ … … 699 674 ul=numpart 700 675 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 if676 ! if (mp_dev_mode) then 677 ! write(*,FMT='(72("#"))') 678 ! write(*,*) "Sending ", num_trans, "particles (from/to)", src_proc, dest_proc 679 ! write(*,*) "numpart before: ", numpart 680 ! end if 706 681 707 682 call MPI_SEND(nclass(ll:ul),num_trans,& … … 742 717 numpart = numpart-num_trans 743 718 744 if (mp_dev_mode) then745 write(*,*) "numpart after: ", numpart746 write(*,FMT='(72("#"))')747 end if719 ! if (mp_dev_mode) then 720 ! write(*,*) "numpart after: ", numpart 721 ! write(*,FMT='(72("#"))') 722 ! end if 748 723 749 724 else if (mp_partid.eq.dest_proc) then … … 756 731 ul=numpart+num_trans 757 732 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 if733 ! if (mp_dev_mode) then 734 ! write(*,FMT='(72("#"))') 735 ! write(*,*) "Receiving ", num_trans, "particles (from/to)", src_proc, dest_proc 736 ! write(*,*) "numpart before: ", numpart 737 ! end if 763 738 764 739 ! We could receive the data directly at nclass(ll:ul) etc., but this leaves … … 810 785 do i=1, num_trans 811 786 ! Take into acount that we may have transferred invalid particles 812 if (itra1_tmp( i).ne.itime) cycle787 if (itra1_tmp(minpart).ne.itime) goto 200 813 788 do ipart=minpart,maxnumpart 814 789 if (itra1(ipart).ne.itime) then 815 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,:) 790 itra1(ipart) = itra1_tmp(minpart) 791 npoint(ipart) = npoint_tmp(minpart) 792 nclass(ipart) = nclass_tmp(minpart) 793 idt(ipart) = idt_tmp(minpart) 794 itramem(ipart) = itramem_tmp(minpart) 795 itrasplit(ipart) = itrasplit_tmp(minpart) 796 xtra1(ipart) = xtra1_tmp(minpart) 797 ytra1(ipart) = ytra1_tmp(minpart) 798 ztra1(ipart) = ztra1_tmp(minpart) 799 xmass1(ipart,:) = xmass1_tmp(minpart,:) 800 ! Increase numpart, if necessary 801 numpart=max(numpart,ipart) 825 802 goto 200 ! Storage space has been found, stop searching 826 803 end if 827 ! :TODO: add check for if too many particles requiried828 804 end do 829 200 minpart= ipart+1805 200 minpart=minpart+1 830 806 end do 831 ! Increase numpart, if necessary 832 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) then 838 write(*,*) "numpart after: ", numpart 839 write(*,FMT='(72("#"))') 840 end if 807 808 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,& 809 & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp) 810 811 ! if (mp_dev_mode) then 812 ! write(*,*) "numpart after: ", numpart 813 ! write(*,FMT='(72("#"))') 814 ! end if 841 815 842 816 else … … 2340 2314 ! 2341 2315 ! 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) 2316 ! take into account nested wind fields 2345 2317 ! 2346 2318 !******************************************************************************* … … 2755 2727 ! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',& 2756 2728 ! & 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() 2729 ! NB: the 'flush' function is possibly a gfortran-specific extension 2730 call flush() 2760 2731 end if 2761 2732 end do -
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
r4c64400 r6b3dad4 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 !*********************************************************** -
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
r4c64400 rc04b739 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 r05cf28d 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. … … 579 574 endif 580 575 581 582 576 ! Determine the release rate (particles per second) and total number 583 577 ! of particles released during the simulation … … 639 633 endif 640 634 641 if (lroot) then642 write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))643 end if644 645 635 return 646 636 -
src/readspecies.f90
rc8fc724 r62e65c7 197 197 weta_gas(pos_spec)=pweta_gas 198 198 wetb_gas(pos_spec)=pwetb_gas 199 crain_aero (pos_spec)=pcrain_aero200 csnow_aero (pos_spec)=pcsnow_aero199 crain_aero=pcrain_aero 200 csnow_aero=pcsnow_aero 201 201 ccn_aero(pos_spec)=pccn_aero 202 202 in_aero(pos_spec)=pin_aero -
src/releaseparticles_mpi.f90
r16b61a5 r861805a 416 416 endif 417 417 end do 418 ! ESO TODO: Here we could use dynamic allocation and increase array sizes 418 419 if (ipart.gt.maxpart_mpi) goto 996 419 420 -
src/timemanager_mpi.f90
r16b61a5 r861805a 578 578 end if 579 579 580 ! Write number ofparticles for all processes580 ! Write particles for all processes 581 581 if (mp_dev_mode) write(*,*) "PID, itime, numpart", mp_pid,itime,numpart 582 582 … … 870 870 endif 871 871 deallocate(gridunc) 872 deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass) 873 if (allocated(checklifetime)) deallocate(checklifetime) 872 deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass, checklifetime) 874 873 deallocate(ireleasestart,ireleaseend,npart,kindz) 875 874 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 r05cf28d 90 90 real :: frac_act, liq_frac, dquer_m 91 91 92 integer (selected_int_kind(16)), dimension(nspec):: blc_count, inc_count92 integer :: blc_count, inc_count 93 93 real :: Si_dummy, wetscav_dummy 94 94 logical :: readclouds_this_nest … … 107 107 !************************ 108 108 109 blc_count (:)=0110 inc_count (:)=0109 blc_count=0 110 inc_count=0 111 111 112 112 do jpart=1,numpart … … 256 256 do ks=1,nspec ! loop over species 257 257 wetdeposit(ks)=0. 258 wetscav=0. 259 260 ! Cycle loop if wet deposition for the species is off 261 if (WETDEPSPEC(ks).eqv..false.) cycle 258 wetscav=0. 262 259 263 260 if (ngrid.gt.0) then … … 277 274 if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then 278 275 ! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then 279 blc_count (ks)=blc_count(ks)+1276 blc_count=blc_count+1 280 277 wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks) 281 278 … … 283 280 !********************************************************************************* 284 281 else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then 285 blc_count (ks)=blc_count(ks)+1282 blc_count=blc_count+1 286 283 287 284 !NIK 17.02.2015 … … 323 320 ! given in species file, or if gas and positive Henry's constant 324 321 if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then 325 inc_count (ks)=inc_count(ks)+1322 inc_count=inc_count+1 326 323 ! if negative coefficients (turned off) set to zero for use in equation 327 324 if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0. … … 435 432 436 433 ! count the total number of below-cloud and in-cloud occurences: 437 tot_blc_count (1:nspec)=tot_blc_count(1:nspec)+blc_count(1:nspec)438 tot_inc_count (1:nspec)=tot_inc_count(1:nspec)+inc_count(1:nspec)434 tot_blc_count=tot_blc_count+blc_count 435 tot_inc_count=tot_inc_count+inc_count 439 436 440 437 end subroutine wetdepo -
src/wetdepokernel.f90
r4c64400 re200b7a 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 !********************************************** … … 123 107 endif 124 108 end do 125 end if126 109 127 110 end subroutine wetdepokernel
Note: See TracChangeset
for help on using the changeset viewer.