import rv4_0 source code from ftp site

AOD_PM_ml.f90

@@ -53,7 +53,7 @@ module AOD_PM_ml
  !// subroutines
   public ::   AOD_calc 
 
-  real, dimension(MAXLIMAX,MAXLJMAX,KMAX_MID) :: kext 
+  real :: kext 
 
   contains
 
@@ -83,7 +83,6 @@ subroutine AOD_calc (i,j,debug)
 
 
   AOD(i,j)     = 0.0
-  kext(i,j,:)  = 0.0
 
   do k =  KCHEMTOP, KMAX_MID   !_______________ vertical layer loop
 
@@ -97,26 +96,27 @@ subroutine AOD_calc (i,j,debug)
 !..=>  xn_2d(ispec,k) * species(ispec)%molwt * 1.e6 / AVOG  [g/m3]
 !.. ===========================================================================
 
-    do n = 1, size(AOD_GROUP)
+   kext  = 0.0
+   do n = 1, size(AOD_GROUP)
       itot = AOD_GROUP(n)
 
-      kext(i,j,k) = kext(i,j,k) +   &
+      kext= kext +   &
                     xn_2d(itot,k) * species(itot)%molwt * species(itot)%ExtC        
     enddo
 
-     kext(i,j,k) = kext(i,j,k) * 1.0e6 / AVOG 
+     kext = kext * 1.0e6 / AVOG 
 
 !     if(debug .and. (k == 18 .or. k == KMAX_MID) )  &
-!            write(6,'(a17,i4,es15.3)') '> Ext. coeff', k, kext(i,j,k)
+!            write(6,'(a17,i4,es15.3)') '> Ext. coeff', k, kext
 
 !.. Aerosol extinction optical depth : integral over all vertical layers
 !.. [1/m} * [m]
 
-      AOD(i,j) = AOD(i,j) + kext(i,j,k) * (z_bnd(i,j,k)-z_bnd(i,j,k+1))
+      AOD(i,j) = AOD(i,j) + kext * (z_bnd(i,j,k)-z_bnd(i,j,k+1))
 
 !      if(debug .and. (k == 18 .or. k == KMAX_MID) )  & 
 !      write(6,'(a25,i4,2es15.4,2f8.1)') '>> Kext AOD for layer', k,  &
-!                kext(i,j,k), AOD(i,j), z_bnd(i,j,k), z_bnd(i,j,k+1)
+!                kext, AOD(i,j), z_bnd(i,j,k), z_bnd(i,j,k+1)
 
   enddo                        !_______________ vertical layer loop
 

AOTnPOD_ml.f90

@@ -145,9 +145,13 @@ subroutine Calc_AOTx(iO3cl,iLC, aot, debug_flag, debug_txt )
              return
     end if
 
-   !If night, or outside growing season, we simply exit with aot=0
-    if ( vego3_outputs(iO3cl)%defn == "EU" .and.  (current_date%hour < 9 .or. &
-         current_date%hour > 21 )) then ! 8-20 CET, assuming summertime
+   ! If night, or outside growing season, we simply exit with aot=0
+   ! EU AOT for 8:00 -- 20:00 CET, is really 8:00 -- 19:59 CET 
+   ! Or:        7:00 -- 18:59 UTC
+   ! (nb hour is integer value)
+
+    if ( vego3_outputs(iO3cl)%defn == "EU" .and.  &
+         (current_date%hour  <  7 .or. current_date%hour  > 18 )) then
           return
 
     else if ( Grid%izen >= AOT_HORIZON ) then  !UN or MM use daylight

Advection_ml.f90

@@ -70,7 +70,7 @@ Module Advection_ml
   use CheckStop_ml,      only : CheckStop
   use Convection_ml,     only : convection_pstar
   use GridValues_ml,     only : GRIDWIDTH_M,xm2,xmd,xm2ji,xmdji, &
-                                carea,xm_i, Pole_included,dA,dB
+                                carea,xm_i, Pole_Singular,dA,dB
   use Io_ml,             only : datewrite
   use ModelConstants_ml, only : KMAX_BND,KMAX_MID,NMET, nstep, nmax, &
                   dt_advec, dt_advec_inv,  PT,KCHEMTOP, NPROCX,NPROCY,NPROC, &
@@ -89,7 +89,7 @@ Module Advection_ml
 
   INCLUDE 'mpif.h'
   INTEGER STATUS(MPI_STATUS_SIZE)
-  real :: MPIbuff(KMAX_MID*max(gimax,gjmax))
+  real,allocatable :: MPIbuff(:)
   integer, private, parameter :: NADVS      =  3
 
   real, private, save, dimension(KMAX_BND)  ::  dhs1, dhs1i, dhs2i
@@ -98,16 +98,17 @@ Module Advection_ml
   real, private, save, dimension(9,2:KMAX_MID,0:1)  ::  alfnew
   real, private, save, dimension(3)  ::  alfbegnew,alfendnew
 
-  real, private,save, dimension(MAXLJMAX,KMAX_MID,NMET) :: uw,ue
+  real, private,save,allocatable, dimension(:,:,:) :: uw,ue
 
-  real, private,save, dimension(MAXLIMAX,KMAX_MID,NMET) :: vs,vn
+  real, private,save,allocatable, dimension(:,:,:) :: vs,vn
 
   integer, public, parameter :: ADVEC_TYPE = 1 ! Divides by advected p*
 ! integer, public, parameter :: ADVEC_TYPE = 2 ! Divides by "meteorologically" 
                                                ! advected p* 
 
   public :: assign_dtadvec
   public :: assign_nmax
+  public :: alloc_adv_arrays
   public :: vgrid
   public :: advecdiff
   public :: advecdiff_poles
@@ -163,6 +164,8 @@ subroutine assign_dtadvec(GRIDWIDTH_M)
 
    if(me==0)write(*,fmt="(a,F8.1,a)")' advection time step (dt_advec) set to: ',dt_advec,' seconds'
 
+   call alloc_adv_arrays!should be moved elsewhere
+
   end subroutine assign_dtadvec
 
   !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -176,7 +179,9 @@ subroutine assign_nmax(metstep)
 
     call CheckStop(mod(3600*metstep,nint(dt_advec)).ne.0, "3600*metstep/dt_advec must be an integer")
 
-    nmax = (3600*metstep)/dt_advec
+   ! Use nint for safety anyway:
+
+    nmax = nint(  (3600*metstep)/dt_advec )
 
     if (me .eq. 0) then
 !      write(6,*)
@@ -684,8 +689,13 @@ subroutine advecdiff_poles
     call Code_timer(tim_before)
 
     if(firstcall)then
-      if(NPROCY>2.and.me==0.and.Pole_included==1)write(*,*)&
+      if(NPROCY>2.and.me==0.and.Pole_Singular>1)then
+         write(*,*)&
           'COMMENT: Advection routine will work faster if NDY = 2 (or 1)'
+      elseif(NPROCY>1.and.me==0.and.Pole_Singular==1)then
+         write(*,*)&
+          'COMMENT: Advection routine will work faster if NDY = 1'
+      endif
     endif
 
     if(KCHEMTOP==2)then
@@ -1394,16 +1404,6 @@ subroutine advvk(xn_adv,ps3d,sdot,dt_s)
         fc3 = fc1*fc2
         n1k = 0
         if(fc1.lt.0)n1k=1
-!pw bug corrected 29/8-2002 (emep1.2beta):
-!       zzfl1 = alfnew(1,k,n1k)*fc1           &
-!             + alfnew(2,k,n1k)*fc2           &
-!             + alfnew(3,k,n1k)*fc3
-!       zzfl2 = alfnew(4,k,n1k)*fc1           &
-!             + alfnew(5,k,n1k)*fc2           &
-!             + alfnew(6,k,n1k)*fc3
-!       zzfl3 = alfnew(7,k,n1k)*fc1           &
-!             + alfnew(8,k,n1k)*fc2           &
-!             + alfnew(9,k,n1k)*fc3
         zzfl1 = alfnew(1,k+1,n1k)*fc1         &
               + alfnew(2,k+1,n1k)*fc2         &
               + alfnew(3,k+1,n1k)*fc3
@@ -3545,4 +3545,16 @@ end subroutine adv_int
 !  moved to Convection_ml.f90
 
 
+
+  subroutine alloc_adv_arrays
+
+    !allocate the arrays once
+    allocate(MPIbuff(KMAX_MID*max(gimax,gjmax)))
+    allocate(uw(MAXLJMAX,KMAX_MID,NMET),ue(MAXLJMAX,KMAX_MID,NMET))
+    allocate(vs(MAXLIMAX,KMAX_MID,NMET),vn(MAXLIMAX,KMAX_MID,NMET))
+
+
+  end subroutine alloc_adv_arrays
+
+
 end module Advection_ml

Aero_Vds_ml.f90

@@ -3,7 +3,7 @@ module Aero_Vds_ml
 !==============================================================================
   use PhysicalConstants_ml, only : FREEPATH, VISCO, BOLTZMANN, PI, GRAV, ROWATER
   use My_Aerosols_ml,       only : NSIZE
-  use ModelConstants_ml,    only : DEBUG_VDS
+  use ModelConstants_ml,    only : DEBUG_VDS, MasterProc
 
   ! DESCRIPTION
   ! Calculates laminar sub-layer resistance (rb) and gravitational settling 
@@ -44,6 +44,13 @@ module Aero_Vds_ml
   public  :: RuijgrokWetSO4
   public  :: Wesely1985
 
+  real, public, parameter, dimension(NSIZE) ::   &
+                 ! diam   = (/ 0.33e-6, 4.0e-6, 8.5e-6 /),  &
+                 !Mc: diam   = (/ 0.33e-6, 1.7e-6, 8.5e-6 /),  &
+                  diam   = (/ 0.33e-6, 3.0e-6, 4.0e-6, 4.5e-6 ,22e-6 /),  &
+                 ! sigma  = (/ 1.8, 2.0, 2.2 /),                    &
+                  sigma  = (/ 1.8, 2.0, 2.0, 2.2 ,2.0/),                    &
+                  PMdens = (/ 1600.0, 2200.0, 2200.0, 2600.0, 800.0/) ! kg/m3
 contains
 
    !------------------------------------------------------------------------
@@ -62,13 +69,13 @@ function SettlingVelocity(tsK,roa) result(Vs)
    ! and dp=1.7 for coarse
    ! Extra 'giant' size is used for sea salt only
 
-     real, parameter, dimension(NSIZE) ::   &
-                ! diam   = (/ 0.33e-6, 4.0e-6, 8.5e-6 /),  &
-                !Mc: diam   = (/ 0.33e-6, 1.7e-6, 8.5e-6 /),  &
-                 diam   = (/ 0.33e-6, 2.5e-6, 8.5e-6 /),  &
-                ! sigma  = (/ 1.8, 2.0, 2.2 /),                    &
-                 sigma  = (/ 1.8, 1.8, 2.2 /),                    &
-                 PMdens = (/ 1600.0, 2200.0, 2200.0 /) ! kg/m3
+!     real, parameter, dimension(NSIZE) ::   &
+!                ! diam   = (/ 0.33e-6, 4.0e-6, 8.5e-6 /),  &
+!                !Mc: diam   = (/ 0.33e-6, 1.7e-6, 8.5e-6 /),  &
+!                 diam   = (/ 0.33e-6, 2.5e-6, 4.0e-6, 4.5e-6 ,22e-6 /),  &
+!                ! sigma  = (/ 1.8, 2.0, 2.2 /),                    &
+!                 sigma  = (/ 1.8, 1.8, 2.0, 2.2 ,2.0/),                    &
+!                 PMdens = (/ 1600.0, 1600.0, 2200.0, 2600.0, 800.0/) ! kg/m3
      real, parameter :: one2three = 1.0/3.0
      integer :: imod 
      real    :: lnsig2, dg, & 
@@ -95,7 +102,11 @@ function SettlingVelocity(tsK,roa) result(Vs)
        !... Settling velocity for poly-disperse 
         vs(imod) = vs_help*(exp(8.0*lnsig2)+1.246*knut*exp(3.5*lnsig2)) ! A31, k=3
 
-        if (DEBUG_VDS) write(6,'(a19,i3,f8.3)') "** Settling Vd **",imod,vs(imod)*100.0
+        if (DEBUG_VDS.and.MasterProc ) &
+            write(6,'(a,i3,es12.3,f10.3,5es12.3,3f9.2,f9.3)') &
+             "** Settling Vd ** ", imod, roa, tsK, &
+             dg,knut,Di_help,vs_help,Di, lnsig2, &
+             1.0e6*diam(imod), PMdens(imod), sigma(imod), vs(imod)*100.0
 
      end do !imod
    end function SettlingVelocity

AirEmis_ml.f90

@@ -46,7 +46,7 @@ module AirEmis_ml
    implicit none
    private
 
-   real, public, dimension(KCHEMTOP:KMAX_MID,MAXLIMAX,MAXLJMAX), save :: &
+   real, public, dimension(:,:,:), save,allocatable :: &
                   airn                 & ! aircraft  NOx emissions
                  ,airlig                 ! lightning NOx emissions
 
@@ -90,7 +90,7 @@ subroutine lightning()
                                                      ! molecules/cm3/s 
 
 
-      character*20 fname
+      character(len=20) :: fname
 
       data ygrdum / 85.76058712, 80.26877907, 74.74454037, &
                 69.21297617, 63.67863556, 58.14295405, &
@@ -113,6 +113,11 @@ subroutine lightning()
       secmonth = 1.
       flux(:,:,:) = 0.
 
+
+      if(.not.allocated(airlig))then
+         allocate(airlig(KCHEMTOP:KMAX_MID,MAXLIMAX,MAXLJMAX))
+      endif
+
 ! --- Read Emission data received from DLR 
 
      if(me == 0)then