Changeset 2753a5c in flexpart.git for src

Timestamp:

Dec 20, 2018, 2:55:22 PM (5 years ago)

Author:

Ignacio Pisso <ip@…>

Branches:

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

Children:

7123c70

Parents:

df96ea65 (diff), 4ad96c5 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

resolved merging conflicts with GFS branch

Location:

src

Files:

• FLEXPART.f90 (modified) (1 diff)
• com_mod.f90 (modified) (1 diff)
• get_wetscav.f90 (modified) (7 diffs)
• netcdf_output_mod.f90 (modified) (5 diffs)
• readspecies.f90 (modified) (8 diffs)
• timemanager.f90 (modified) (1 diff)
• makefile.gfs (deleted)
• par_mod.f90 (modified) (2 diffs)
• pathnames (modified) (1 diff)

src/FLEXPART.f90

@@ -455 +455 @@
   call timemanager(metdata_format)
 
+  if (verbosity.gt.0) then
 ! NIK 16.02.2005 
-  do i=1,nspec
-    write(*,*) '**********************************************'
-    write(*,*) 'Scavenging statistics for species ', species(i), ':'
-    write(*,*) 'Total number of occurences of below-cloud scavenging', &
-         & tot_blc_count(i)
-    write(*,*) 'Total number of occurences of in-cloud    scavenging', &
-         & tot_inc_count(i)
-    write(*,*) '**********************************************'
-  end do
+    do i=1,nspec
+      if (tot_inc_count(i).gt.0) then
+         write(*,*) '**********************************************'
+         write(*,*) 'Scavenging statistics for species ', species(i), ':'
+         write(*,*) 'Total number of occurences of below-cloud scavenging', &
+           & tot_blc_count(i)
+         write(*,*) 'Total number of occurences of in-cloud    scavenging', &
+           & tot_inc_count(i)
+         write(*,*) '**********************************************'
+      endif
+    end do
+    write (*,*) 'timemanager> call wetdepo'
+  endif
   
   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&

src/com_mod.f90

@@ -174 +174 @@
   real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
   real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
-  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
+  real :: rm(maxspec),dryvel(maxspec)
   real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
-  ! se  it is possible to associate a species with a second one to make transfer from gas to aerosol
-  integer :: spec_ass(maxspec)
 
   real :: area_hour(maxspec,24),point_hour(maxspec,24)

src/get_wetscav.f90

@@ -151 +151 @@
       if (height(il).gt.ztra1(jpart)) then
         hz=il-1
-!        goto 26
         exit
       endif
     end do
-!26  continue
-
 
     if (ngrid.eq.0) then
@@ -203 +200 @@
 
 
-!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 
+!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp - 2 and 3 not used removed by SE
 ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
 ! for now they are treated the same
     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
-    grfraction(2)=max(0.05,cc*(lfr(i)))
-    grfraction(3)=max(0.05,cc*(cfr(j)))
-
 
 ! 2) Computation of precipitation rate in sub-grid cell
 !******************************************************
     prec(1)=(lsp+convp)/grfraction(1)
-    prec(2)=(lsp)/grfraction(2)
-    prec(3)=(convp)/grfraction(3)
-
 
 ! 3) Computation of scavenging coefficients for all species
 !    Computation of wet deposition
 !**********************************************************
-
 
       if (ngrid.gt.0) then
@@ -237 +227 @@
 !******************************************************************
         if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
-          !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then 
           blc_count(ks)=blc_count(ks)+1
           wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
@@ -271 +260 @@
           endif
           
-!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
-
         endif ! gas or particle
 !      endif ! positive below-cloud scavenging parameters given in Species file
@@ -285 +272 @@
         if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then 
           inc_count(ks)=inc_count(ks)+1
-!          write(*,*) 'Incloud: ',inc_count
 ! if negative coefficients (turned off) set to zero for use in equation
           if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
@@ -300 +286 @@
 ! sec test
 !           cl=1E6*1E-7*prec(1)**0.3 !Sec GFS new
-            cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new
+            cl=1E6*2E-7*prec(1)**0.36 !Sec ECMWF new, is also suitable for GFS
 !           cl=2E-7*prec(1)**0.36 !Andreas
 !           cl=1.6E-6*prec(1)**0.36 !Henrik
@@ -322 +308 @@
           if (dquer(ks).gt.0.) then
             S_i= frac_act/cl
-!           write(*,*) 'Si: ',S_i
-
 ! GAS
 !****
           else
-
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-!REPLACE to switch old/ new scheme 
-          ! S_i=frac_act/cle
             S_i=1/cle
           endif ! gas or particle
 
 ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
-!OLD 
-          if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
+!SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
-          else
-!SEC wetscav fix
-             wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
-!            wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
-          endif
         endif ! positive in-cloud scavenging parameters given in Species file
       endif !incloud

src/netcdf_output_mod.f90

@@ -67 +67 @@
                        ccn_aero,in_aero, & ! wetc_in,wetd_in, &
                        reldiff,henry,f0,density,dquer,dsigma,dryvel,&
-!                       weightmolar,ohreact,spec_ass,kao,vsetaver,&
-                       weightmolar,ohcconst,ohdconst,spec_ass,kao,vsetaver,&
+                       weightmolar,ohcconst,ohdconst,vsetaver,&
                        ! for concoutput_netcdf and concoutput_nest_netcdf
                        nxmin1,nymin1,nz,oro,oron,rho,rhon,&
@@ -114 +113 @@
   logical, parameter :: write_area = .false.
 
+  ! coordinate transformation from internal to world coord
+  real :: xp1,yp1,xp2,yp2
 contains
 
@@ -512 +513 @@
         call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
-        call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
-        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
 
         if (lnest) then
@@ -535 +534 @@
         call nf90_err(nf90_put_att(ncid, sID, 'ohcconst', ohcconst(i)))
         call nf90_err(nf90_put_att(ncid, sID, 'ohdconst', ohdconst(i)))
-        call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i))) 
         call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
-        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
 
         if (lnest) then
@@ -673 +670 @@
        call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
        call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
-       call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
+       xp1=xpoint1(i)*dx+xlon0
+       yp1=ypoint1(i)*dy+ylat0
+       xp2=xpoint2(i)*dx+xlon0
+       yp2=ypoint2(i)*dy+ylat0
+       call nf90_err(nf90_put_var(ncid, rellng1ID, xp1, (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellng2ID, xp2, (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellat1ID, yp1, (/i/)))
+       call nf90_err(nf90_put_var(ncid, rellat2ID, yp2, (/i/)))
        call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
        call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))

src/readspecies.f90

@@ -67 +67 @@
   character(len=16) :: pspecies
   real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
-  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pkao
+  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
   real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
-  integer :: readerror, pspec_ass
+  real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
+  integer :: readerror
 
 ! declare namelist
@@ -76 +77 @@
        pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
        preldiff, phenry, pf0, pdensity, pdquer, &
-       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
+       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
+       parea_dow, parea_hour, ppoint_dow, ppoint_hour
 
   pspecies="" ! read failure indicator value
@@ -96 +98 @@
   pohdconst=-9.9E-09
   pohnconst=2.0
-  pspec_ass=-9
-  pkao=-99.99
   pweightmolar=-999.9
 
+  do j=1,24           ! initialize everything to no variation
+    area_hour(pos_spec,j)=1.
+    point_hour(pos_spec,j)=1.
+  end do
+  do j=1,7
+    area_dow(pos_spec,j)=1.
+    point_dow(pos_spec,j)=1.
+  end do
+
+  if (readerror.ne.0) then ! text format input
 ! Open the SPECIES file and read species names and properties
 !************************************************************
@@ -164 +174 @@
     read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
 !  write(*,*) ohnconst(pos_spec)
-    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
-!  write(*,*) spec_ass(pos_spec)
-    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
-!       write(*,*) kao(pos_spec)
+
+! Read in daily and day-of-week variation of emissions, if available
+!*******************************************************************
+
+    read(unitspecies,*,end=22)
+    do j=1,24     ! 24 hours, starting with 0-1 local time
+      read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
+    end do
+    read(unitspecies,*)
+    do j=1,7      ! 7 days of the week, starting with Monday
+      read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
+    end do
 
     pspecies=species(pos_spec)
@@ -188 +206 @@
     pohdconst=ohdconst(pos_spec)
     pohnconst=ohnconst(pos_spec)
-    pspec_ass=spec_ass(pos_spec)
-    pkao=kao(pos_spec)
-
-  else
+
+
+    do j=1,24     ! 24 hours, starting with 0-1 local time
+      parea_hour(j)=area_hour(pos_spec,j)
+      ppoint_hour(j)=point_hour(pos_spec,j)
+    end do
+    do j=1,7      ! 7 days of the week, starting with Monday
+      parea_dow(j)=area_dow(pos_spec,j)
+      ppoint_dow(j)=point_dow(pos_spec,j)
+    end do
+
+  else ! namelist available
 
     species(pos_spec)=pspecies
@@ -212 +238 @@
     ohdconst(pos_spec)=pohdconst
     ohnconst(pos_spec)=pohnconst
-    spec_ass(pos_spec)=pspec_ass
-    kao(pos_spec)=pkao
+
+    do j=1,24     ! 24 hours, starting with 0-1 local time
+      area_hour(pos_spec,j)=parea_hour(j)
+      point_hour(pos_spec,j)=ppoint_hour(j)
+    end do
+    do j=1,7      ! 7 days of the week, starting with Monday
+      area_dow(pos_spec,j)=parea_dow(j)
+      point_dow(pos_spec,j)=ppoint_dow(j)
+    end do
 
   endif
@@ -303 +336 @@
   end if
 
-  if (spec_ass(pos_spec).gt.0) then
-    spec_found=.FALSE.
-    do j=1,pos_spec-1
-      if (spec_ass(pos_spec).eq.specnum(j)) then
-        spec_ass(pos_spec)=j
-        spec_found=.TRUE.
-        ASSSPEC=.TRUE.
-      endif
-    end do
-    if (spec_found.eqv..false.) then
-      goto 997
-    endif
-  endif
-
-  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
   if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
 
@@ -330 +348 @@
 
 
-! Read in daily and day-of-week variation of emissions, if available
-!*******************************************************************
-! HSO: This is not yet implemented as namelist parameters
-
-  do j=1,24           ! initialize everything to no variation
-    area_hour(i,j)=1.
-    point_hour(i,j)=1.
-  end do
-  do j=1,7
-    area_dow(i,j)=1.
-    point_dow(i,j)=1.
-  end do
-
-  if (readerror.ne.0) then ! text format input
-
-    read(unitspecies,*,end=22)
-    do j=1,24     ! 24 hours, starting with 0-1 local time
-      read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
-    end do
-    read(unitspecies,*)
-    do j=1,7      ! 7 days of the week, starting with Monday
-      read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
-    end do
-
   endif
 

src/timemanager.f90

@@ -449 +449 @@
 
         !write(*,46) float(itime)/3600,itime,numpart
-45      format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES:    Uncertainty: ',3f7.3)
+45      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
 46      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions

src/par_mod.f90

@@ -148 +148 @@
   ! ECMWF
 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=359 ! 1.0 degree 92 level
+!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0 ! 1.0 degree 138 level
 !   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 1.0 degree 138 level
 ! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359  ! 0.5 degree 138 level
@@ -153 +154 @@
 
 ! GFS
-   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0
+   integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
+   integer :: nxshift=0
 
 

src/pathnames

@@ -1 +1 @@
 ../options/
-../output/
+./output/
 /
 /xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global