Changeset 38b7917 in flexpart.git for src/init_domainfill.f90

Timestamp:

Mar 8, 2016, 9:54:38 AM (8 years ago)

Author:

Espen Sollum ATMOS <eso@…>

Branches:

master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie

Children:

9b53903

Parents:

db712a8

Message:

Parallelization of domain fill option (save/restart not implemented yet)

File:

• src/init_domainfill.f90 (modified) (22 diffs)

src/init_domainfill.f90

@@ -21 +21 @@
 
 subroutine init_domainfill
-  !
-  !*****************************************************************************
-  !                                                                            *
-  ! Initializes particles equally distributed over the first release location  *
-  ! specified in file RELEASES. This box is assumed to be the domain for doing *
-  ! domain-filling trajectory calculations.                                    *
-  ! All particles carry the same amount of mass which alltogether comprises the*
-  ! mass of air within the box.                                                *
-  !                                                                            *
-  !     Author: A. Stohl                                                       *
-  !                                                                            *
-  !     15 October 2002                                                        *
-  !                                                                            *
-  !*****************************************************************************
-  !                                                                            *
-  ! Variables:                                                                 *
-  !                                                                            *
-  ! numparticlecount    consecutively counts the number of particles released  *
-  ! nx_we(2)       grid indices for western and eastern boundary of domain-    *
-  !                filling trajectory calculations                             *
-  ! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
-  !                filling trajectory calculations                             *
-  !                                                                            *
-  !*****************************************************************************
+!
+!*****************************************************************************
+!                                                                            *
+! Initializes particles equally distributed over the first release location  *
+! specified in file RELEASES. This box is assumed to be the domain for doing *
+! domain-filling trajectory calculations.                                    *
+! All particles carry the same amount of mass which alltogether comprises the*
+! mass of air within the box.                                                *
+!                                                                            *
+!     Author: A. Stohl                                                       *
+!                                                                            *
+!     15 October 2002                                                        *
+!                                                                            *
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+!                                                                            *
+! numparticlecount    consecutively counts the number of particles released  *
+! nx_we(2)       grid indices for western and eastern boundary of domain-    *
+!                filling trajectory calculations                             *
+! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
+!                filling trajectory calculations                             *
+!                                                                            *
+!*****************************************************************************
 
   use point_mod
@@ -65 +65 @@
 
 
-  ! Determine the release region (only full grid cells), over which particles
-  ! shall be initialized
-  ! Use 2 fields for west/east and south/north boundary
-  !**************************************************************************
+! Determine the release region (only full grid cells), over which particles
+! shall be initialized
+! Use 2 fields for west/east and south/north boundary
+!**************************************************************************
 
   nx_we(1)=max(int(xpoint1(1)),0)
@@ -75 +75 @@
   ny_sn(2)=min((int(ypoint2(1))+1),nymin1)
 
-  ! For global simulations (both global wind data and global domain-filling),
-  ! set a switch, such that no boundary conditions are used
-  !**************************************************************************
+! For global simulations (both global wind data and global domain-filling),
+! set a switch, such that no boundary conditions are used
+!**************************************************************************
   if (xglobal.and.sglobal.and.nglobal) then
     if ((nx_we(1).eq.0).and.(nx_we(2).eq.nxmin1).and. &
@@ -87 +87 @@
   endif
 
-  ! Do not release particles twice (i.e., not at both in the leftmost and rightmost
-  ! grid cell) for a global domain
-  !*****************************************************************************
+! Do not release particles twice (i.e., not at both in the leftmost and rightmost
+! grid cell) for a global domain
+!*****************************************************************************
   if (xglobal) nx_we(2)=min(nx_we(2),nx-2)
 
 
-  ! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
-  ! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
-  !************************************************************
+! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
+! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
+!************************************************************
 
   do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
@@ -117 +117 @@
   end do
 
-  ! Do the same for the south pole
+! Do the same for the south pole
 
   if (sglobal) then
@@ -130 +130 @@
   endif
 
-  ! Do the same for the north pole
+! Do the same for the north pole
 
   if (nglobal) then
@@ -144 +144 @@
 
 
-  ! Calculate total mass of each grid column and of the whole atmosphere
-  !*********************************************************************
+! Calculate total mass of each grid column and of the whole atmosphere
+!*********************************************************************
 
   colmasstotal=0.
@@ -157 +157 @@
   end do
 
-           write(*,*) 'Atm. mass: ',colmasstotal
+  write(*,*) 'Atm. mass: ',colmasstotal
 
 
   if (ipin.eq.0) numpart=0
 
-  ! Determine the particle positions
-  !*********************************
+! Determine the particle positions
+!*********************************
 
   numparttot=0
@@ -175 +175 @@
       if (ncolumn.gt.numcolumn) numcolumn=ncolumn
 
-  ! Calculate pressure at the altitudes of model surfaces, using the air density
-  ! information, which is stored as a 3-d field
-  !*****************************************************************************
+! Calculate pressure at the altitudes of model surfaces, using the air density
+! information, which is stored as a 3-d field
+!*****************************************************************************
 
       do kz=1,nz
@@ -191 +191 @@
 
 
-  ! For columns with many particles (i.e. around the equator), distribute
-  ! the particles equally, for columns with few particles (i.e. around the
-  ! poles), distribute the particles randomly
-  !***********************************************************************
+! For columns with many particles (i.e. around the equator), distribute
+! the particles equally, for columns with few particles (i.e. around the
+! poles), distribute the particles randomly
+!***********************************************************************
 
 
@@ -209 +209 @@
             dz=1./(dz1+dz2)
 
-  ! Assign particle position
-  !*************************
-  ! Do the following steps only if particles are not read in from previous model run
-  !*****************************************************************************
+! Assign particle position
+!*************************
+! Do the following steps only if particles are not read in from previous model run
+!*****************************************************************************
             if (ipin.eq.0) then
               xtra1(numpart+jj)=real(ix)-0.5+ran1(idummy)
@@ -224 +224 @@
 
 
-  ! Interpolate PV to the particle position
-  !****************************************
+! Interpolate PV to the particle position
+!****************************************
               ixm=int(xtra1(numpart+jj))
               jym=int(ytra1(numpart+jj))
@@ -260 +260 @@
 
 
-  ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
-  !*****************************************************************************
+! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
+!*****************************************************************************
 
               if (((ztra1(numpart+jj).gt.3000.).and. &
                    (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
 
-  ! Assign certain properties to the particle
-  !******************************************
+! Assign certain properties to the particle
+!******************************************
                 nclass(numpart+jj)=min(int(ran1(idummy)* &
                      real(nclassunc))+1,nclassunc)
@@ -293 +293 @@
   end do
 
-
-  ! Check whether numpart is really smaller than maxpart
-  !*****************************************************
-
+  write(*,*) 'init_domainfill> ncolumn: ', ncolumn 
+  write(*,*) 'init_domainfill> numcolumn: ', numcolumn 
+  write(*,*) 'init_domainfill> ny_sn(1),ny_sn(2): ', ny_sn(1),ny_sn(2)
+  write(*,*) 'init_domainfill> nx_we(1),nx_we(2): ', nx_we(1),nx_we(2)
+
+
+! Check whether numpart is really smaller than maxpart
+!*****************************************************
+
+! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier
   if (numpart.gt.maxpart) then
     write(*,*) 'numpart too large: change source in init_atm_mass.f'
@@ -306 +312 @@
 
 
-  ! Make sure that all particles are within domain
-  !***********************************************
+! Make sure that all particles are within domain
+!***********************************************
 
   do j=1,numpart
@@ -319 +325 @@
 
 
-  ! For boundary conditions, we need fewer particle release heights per column,
-  ! because otherwise it takes too long until enough mass has accumulated to
-  ! release a particle at the boundary (would take dx/u seconds), leading to
-  ! relatively large position errors of the order of one grid distance.
-  ! It's better to release fewer particles per column, but to do so more often.
-  ! Thus, use on the order of nz starting heights per column.
-  ! We thus repeat the above to determine fewer starting heights, that are
-  ! used furtheron in subroutine boundcond_domainfill.f.
-  !****************************************************************************
+! For boundary conditions, we need fewer particle release heights per column,
+! because otherwise it takes too long until enough mass has accumulated to
+! release a particle at the boundary (would take dx/u seconds), leading to
+! relatively large position errors of the order of one grid distance.
+! It's better to release fewer particles per column, but to do so more often.
+! Thus, use on the order of nz starting heights per column.
+! We thus repeat the above to determine fewer starting heights, that are
+! used furtheron in subroutine boundcond_domainfill.f.
+!****************************************************************************
 
   fractus=real(numcolumn)/real(nz)
@@ -343 +349 @@
 
 
-  ! Memorize how many particles per column shall be used for all boundaries
-  ! This is further used in subroutine boundcond_domainfill.f
-  ! Use 2 fields for west/east and south/north boundary
-  !************************************************************************
+! Memorize how many particles per column shall be used for all boundaries
+! This is further used in subroutine boundcond_domainfill.f
+! Use 2 fields for west/east and south/north boundary
+!************************************************************************
 
       if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn
@@ -353 +359 @@
       if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn
 
-  ! Calculate pressure at the altitudes of model surfaces, using the air density
-  ! information, which is stored as a 3-d field
-  !*****************************************************************************
+! Calculate pressure at the altitudes of model surfaces, using the air density
+! information, which is stored as a 3-d field
+!*****************************************************************************
 
       do kz=1,nz
@@ -361 +367 @@
       end do
 
-  ! Determine the reference starting altitudes
-  !*******************************************
+! Determine the reference starting altitudes
+!*******************************************
 
       deltacol=(pp(1)-pp(nz))/real(ncolumn)
@@ -374 +380 @@
             dz=1./(dz1+dz2)
             zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
-           if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
-
-  ! Memorize vertical positions where particles are introduced
-  ! This is further used in subroutine boundcond_domainfill.f
-  !***********************************************************
+            if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
+
+! Memorize vertical positions where particles are introduced
+! This is further used in subroutine boundcond_domainfill.f
+!***********************************************************
 
             if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition
@@ -385 +391 @@
             if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition
 
-  ! Initialize mass that has accumulated at boundary to zero
-  !*********************************************************
+! Initialize mass that has accumulated at boundary to zero
+!*********************************************************
 
             acc_mass_we(1,jy,j)=0.
@@ -399 +405 @@
   end do
 
-  ! If particles shall be read in to continue an existing run,
-  ! then the accumulated masses at the domain boundaries must be read in, too.
-  ! This overrides any previous calculations.
-  !***************************************************************************
+! If particles shall be read in to continue an existing run,
+! then the accumulated masses at the domain boundaries must be read in, too.
+! This overrides any previous calculations.
+!***************************************************************************
 
   if (ipin.eq.1) then