Emit stuff at one grid point

DYNAMICO_emission/src/emission.F90

@@ -11,6 +11,8 @@ MODULE emission
   
   REAL, PARAMETER :: oneday = 86400. ! hard-coded                                                                                                                                                         
   INTEGER, PARAMETER :: log_unit = 15
+
+  INTEGER :: emission_ij ! index of cell where volcanic emissions occur
   
   !  TYPE(t_field),POINTER :: f_write2d(:), f_write_llm(:), f_write_llmp1(:)
   
@@ -27,13 +29,52 @@ CONTAINS
     USE getin_mod, ONLY : getin
     USE physics_interface_mod, ONLY : inout => physics_inout
     REAL(rstd) :: unjours
-
-    WRITE(*,*) 'init_emission called'
     
     unjours = 86400.
     CALL getin('unjours', unjours)
-    
+
+    CALL init_emission(inout%lon, inout%lat)
   END SUBROUTINE init_physics
+
+  SUBROUTINE init_emission(lon,lat)
+    USE icosa,              ONLY : rstd
+    USE getin_mod,          ONLY : getin
+    USE spherical_geom_mod, ONLY : dist_lonlat ! distance on the sphere
+    REAL(rstd), INTENT(IN) :: lon(:), lat(:)
+    REAL(rstd) :: emission_lon, emission_lat, dist, dist_min, dist_max
+    INTEGER :: ij, ij_min
+    
+    WRITE(*,*) 'init_emission called'
+
+    emission_lon = 0.
+    CALL getin('emission_lon', emission_lon)
+    emission_lat = 0.
+    CALL getin('emission_lat', emission_lat)
+    dist_max = 1e-4
+    CALL getin('emission_dist_max', dist_max)
+
+    dist_min = 4. ! distances on the unit sphere are always less than pi
+    ij_min = 1
+    DO ij=1, SIZE(lon)
+       CALL dist_lonlat(emission_lon, emission_lat, lon(ij), lat(ij), dist)
+       IF(dist<dist_min) THEN
+          dist_min = dist
+          ij_min = ij
+       END IF
+    END DO
+
+    PRINT *, 'emission lon, lat    :', emission_lon, emission_lat
+    PRINT *, 'closest  lon, lat    :', lon(ij_min), lat(ij_min)
+    PRINT *, 'dist to closest cell :', ij, dist_min
+
+    IF(dist_min < dist_max) THEN
+       emission_ij = ij_min
+    ELSE
+       ! if the closest point is too far from the target, this means that
+       ! this points is on another MPI process
+       emission_ij = -1
+    END IF
+  END SUBROUTINE init_emission
   
   SUBROUTINE physics
     USE mpipara,  ONLY : is_mpi_master
@@ -66,8 +107,25 @@ CONTAINS
     !           &        inout%p, play, pphi,                   &
     !           &        inout%ulon,  inout%ulat,  inout%temp,  &
     !           &        inout%dulon, inout%dulat, inout%dtemp, dps)
+
+    IF( emission_ij>0 ) THEN
+       inout%dq(:,:,:)=0.
+       CALL emit(llm, time, inout%p(emission_ij, :), inout%geopot(emission_ij,:), &
+            &    inout%dq(emission_ij,:,1) )
+    END IF
     
   END SUBROUTINE physics
+
+  SUBROUTINE emit(llm, time, p, geopot, dq)
+    USE icosa,    ONLY : rstd
+    INTEGER, INTENT(IN)     :: llm           ! number of model layers
+    REAL(rstd), INTENT(IN)  :: time          ! current time in seconds since start of simulation
+    REAL(rstd), INTENT(IN)  :: p(:), &   ! pressure at interfaces
+         &                     geopot(:) ! geopotential at interfaces
+    REAL(rstd), INTENT(OUT) :: dq(:)
+    ! uniform emission on the whole column
+    dq(:) = 1.
+  END SUBROUTINE emit
   
   SUBROUTINE compute_play(ngrid, llm, plev, play)
     INTEGER, INTENT(IN) :: ngrid, llm