Changeset 660bcb7 in flexpart.git for src/timemanager.f90

Timestamp:

Nov 24, 2014, 1:36:35 PM (9 years ago)

Author:

Anne Fouilloux <annefou@…>

Branches:

NetCDF

Children:

3073eaf

Parents:

4bf4a69

Message:

NETCDF outputs from svn trunk (hasod: ADD: netcdf module file)
I did not take changes in advance.f90; it will be added later.
I also changed OPENs where status was set to new and access is set to APPEND (had problems on abel.uio.no with intel compilers 2013.sp1.3)
It should contains technical changes only and results should be identical.

File:

• src/timemanager.f90 (modified) (16 diffs)

src/timemanager.f90

@@ -43 +43 @@
   !     call convection BEFORE new fields are read in BWD mode                 *
   !  Changes Caroline Forster, Feb 2005                                        *
-  !new interface between flexpart and convection scheme                        *
-  !Emanuel's latest subroutine convect43c.f is used                            *
+  !   new interface between flexpart and convection scheme                     *
+  !   Emanuel's latest subroutine convect43c.f is used                         *
+  !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
+  !   added netcdf output code                                                 *
   !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
-  ! DEP                .true. if either wet or dry deposition is switched on   *
+  ! dep                .true. if either wet or dry deposition is switched on   *
   ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
-  ! DRYDEP             .true. if dry deposition is switched on                 *
+  ! drydep             .true. if dry deposition is switched on                 *
   ! ideltas [s]        modelling period                                        *
   ! itime [s]          actual temporal position of calculation                 *
@@ -70 +72 @@
   ! prob               probability of absorption at ground due to dry          *
   !                    deposition                                              *
-  ! WETDEP             .true. if wet deposition is switched on                 *
+  ! wetdep             .true. if wet deposition is switched on                 *
   ! weight             weight for each concentration sample (1/2 or 1)         *
   ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
@@ -92 +94 @@
   use par_mod
   use com_mod
+  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
 
   implicit none
@@ -129 +132 @@
 
 
+  itime=0
   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
   write(*,46) float(itime)/3600,itime,numpart
@@ -141 +145 @@
   do itime=0,ideltas,lsynctime
 
-
   ! Computation of wet deposition, OH reaction and mass transfer
   ! between two species every lsynctime seconds
@@ -152 +155 @@
   !********************************************************************
 
-    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
         if (verbosity.gt.0) then
            write (*,*) 'timemanager> call wetdepo'
@@ -159 +162 @@
     endif
 
-    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
+    if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
          call ohreaction(itime,lsynctime,loutnext)
 
-    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
        stop 'associated species not yet implemented!'
   !     call transferspec(itime,lsynctime,loutnext)
@@ -240 +243 @@
   !***********************************************************************
 
-    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
+    if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
       do ks=1,nspec
       do kp=1,maxpointspec_act
@@ -350 +353 @@
         if ((iout.le.3.).or.(iout.eq.5)) then
           if (surf_only.ne.1) then 
-          call concoutput(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
+            if (lnetcdfout.eq.1) then 
+              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            else 
+              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            endif
           else  
-  if (verbosity.eq.1) then
-     print*,'call concoutput_surf '
-     CALL SYSTEM_CLOCK(count_clock)
-     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
-          call concoutput_surf(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
-  if (verbosity.eq.1) then
-     print*,'called concoutput_surf '
-     CALL SYSTEM_CLOCK(count_clock)
-     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
+            if (verbosity.eq.1) then
+             print*,'call concoutput_surf '
+             call system_clock(count_clock)
+             write(*,*) 'system clock',count_clock - count_clock0   
+            endif
+            if (lnetcdfout.eq.1) then 
+              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            else
+              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+              if (verbosity.eq.1) then
+                print*,'called concoutput_surf '
+                call system_clock(count_clock)
+                write(*,*) 'system clock',count_clock - count_clock0   
+              endif
+            endif
           endif
 
-          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
-          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
+          if (nested_output .eq. 1) then
+            if (lnetcdfout.eq.0) then
+              if (surf_only.ne.1) then
+                call concoutput_nest(itime,outnum)
+              else 
+                call concoutput_surf_nest(itime,outnum)
+              endif
+            else
+              if (surf_only.ne.1) then
+                call concoutput_nest_netcdf(itime,outnum)
+              else 
+                call concoutput_surf_nest_netcdf(itime,outnum)
+              endif
+            endif
+          endif
           outnum=0.
         endif
@@ -474 +496 @@
   !****************************
 
-        xold=xtra1(j)
-        yold=ytra1(j)
+        xold=real(xtra1(j))
+        yold=real(ytra1(j))
         zold=ztra1(j)
 
@@ -515 +537 @@
             endif
 
-            if (DRYDEPSPEC(ks)) then        ! dry deposition
+            if (drydepspec(ks)) then        ! dry deposition
               drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
               xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
@@ -526 +548 @@
             endif
 
-
             if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
@@ -538 +559 @@
           if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
             itra1(j)=-999999999
+            if (verbosity.gt.0) then
+              print*,'terminated particle ',j,' for small mass'
+            endif
           endif
 
   !        Sabine Eckhardt, June 2008
   !        don't create depofield for backward runs
-          if (DRYDEP.AND.(ldirect.eq.1)) then
-            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
-                 real(ytra1(j)),nage,kp)
-            if (nested_output.eq.1) call drydepokernel_nest( &
-                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
-                 nage,kp)
+          if (drydep.AND.(ldirect.eq.1)) then
+            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+            if (nested_output.eq.1) then
+              call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+            endif
           endif
 
@@ -554 +577 @@
 
           if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
-            if (linit_cond.ge.1) &
-                 call initial_cond_calc(itime+lsynctime,j)
+            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
             itra1(j)=-999999999
+            if (verbosity.gt.0) then
+              print*,'terminated particle ',j,' for age'
+            endif
           endif
         endif
@@ -584 +609 @@
 
   if (iflux.eq.1) then
-      deallocate(flux)
+    deallocate(flux)
   endif
-  if (OHREA.eqv..TRUE.) then
-      deallocate(OH_field,OH_field_height)
+  if (ohrea.eqv..TRUE.) then
+    deallocate(OH_field,OH_field_height)
   endif
   if (ldirect.gt.0) then
-  deallocate(drygridunc,wetgridunc)
+    deallocate(drygridunc,wetgridunc)
   endif
   deallocate(gridunc)
@@ -597 +622 @@
   deallocate(xmasssave)
   if (nested_output.eq.1) then
-     deallocate(orooutn, arean, volumen)
-     if (ldirect.gt.0) then
-     deallocate(griduncn,drygriduncn,wetgriduncn)
-     endif
+    deallocate(orooutn, arean, volumen)
+    if (ldirect.gt.0) then
+      deallocate(griduncn,drygriduncn,wetgriduncn)
+    endif
   endif
   deallocate(outheight,outheighthalf)
-  deallocate(oroout, area, volume)
+  deallocate(oroout,area,volume)
 
 end subroutine timemanager