Index: /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain.F	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain.F	(revision 22)
@@ -921,4 +921,18 @@
       grid%ntsloc        = 0
       grid%ntsloc_domain = 0
+
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+      ALLOCATE( grid%latlidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%lonlidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%namelidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%desclidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%ilidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%jlidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%id_lidarloc( grid%max_lidar_locs ) )
+      ALLOCATE( grid%lidar_filename( grid%max_lidar_locs ) )
+      grid%nlidarloc        = 0
+      grid%nlidarloc_domain = 0
+#endif
 
 #ifdef DM_PARALLEL
@@ -1592,4 +1606,31 @@
         DEALLOCATE( grid%ts_filename )
       ENDIF 
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+      IF ( ASSOCIATED( grid%ilidarloc ) ) THEN
+        DEALLOCATE( grid%ilidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%jlidarloc ) ) THEN
+        DEALLOCATE( grid%jlidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%id_lidarloc ) ) THEN
+        DEALLOCATE( grid%id_lidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%latlidarloc ) ) THEN
+        DEALLOCATE( grid%latlidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%lonlidarloc ) ) THEN
+        DEALLOCATE( grid%lonlidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%namelidarloc ) ) THEN
+        DEALLOCATE( grid%namelidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%desclidarloc ) ) THEN
+        DEALLOCATE( grid%desclidarloc )
+      ENDIF 
+      IF ( ASSOCIATED( grid%lidar_filename ) ) THEN
+        DEALLOCATE( grid%lidar_filename )
+      ENDIF 
+#endif
       DEALLOCATE( grid )
       NULLIFY( grid )
Index: /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain_type.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain_type.F	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/frame/module_domain_type.F	(revision 22)
@@ -222,4 +222,18 @@
       LOGICAL :: have_displayed_alloc_stats   ! used in module_alloc_space to display alloc stats; only do it once.
 
+
+#ifdef WRFMEAS
+! L. Fita, LMD. May 2014
+      ! lidar location information
+      INTEGER :: nlidarloc, nlidarloc_domain
+      INTEGER :: next_lidar_time
+      INTEGER, POINTER, DIMENSION(:) :: ilidarloc, jlidarloc, id_lidarloc
+      REAL, POINTER, DIMENSION(:) :: latlidarloc, lonlidarloc
+      CHARACTER (LEN=5), POINTER, DIMENSION(:) :: namelidarloc
+      CHARACTER (LEN=25), POINTER, DIMENSION(:) :: desclidarloc
+      CHARACTER (LEN=256), POINTER, DIMENSION(:) :: lidar_filename
+      LOGICAL :: have_calculated_lidarlocs
+#endif
+
    END TYPE domain
 END MODULE module_domain_type
Index: /branches/LMDZ_WRFmeas/WRFV3/share/Makefile
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/Makefile	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/Makefile	(revision 22)
@@ -36,4 +36,5 @@
         output_wrf.o                    \
         wrf_timeseries.o                \
+        wrf_lidar.o                     \
         wrf_ext_write_field.o           \
         wrf_ext_read_field.o            \
@@ -42,4 +43,5 @@
         wrf_bdyin.o             	\
         wrf_tsin.o                      \
+        wrf_lidarin.o                   \
         landread.o			\
 	setfeenv.o
@@ -82,5 +84,5 @@
 		../frame/module_wrf_error.o
 
-start_domain.o: start_domain_em.int wrf_timeseries.o ../frame/module_domain.o ../frame/module_configure.o
+start_domain.o: start_domain_em.int wrf_timeseries.o wrf_lidar.o ../frame/module_domain.o ../frame/module_configure.o
 
 module_bc.o: ../frame/module_configure.o ../frame/module_state_description.o \
@@ -109,4 +111,11 @@
 
 wrf_timeseries.o: wrf_tsin.o \
+		module_model_constants.o \
+		module_llxy.o \
+		../frame/module_domain.o \
+		../frame/module_configure.o \
+		../frame/module_dm.o
+
+wrf_lidar.o: wrf_lidarin.o \
 		module_model_constants.o \
 		module_llxy.o \
@@ -293,4 +302,5 @@
                                 ../frame/module_configure.o ../frame/module_io.o module_io_wrf.o module_bc_time_utilities.o
 wrf_tsin.o	                : ../frame/module_domain.o
+wrf_lidarin.o	                : ../frame/module_domain.o
 
 
Index: /branches/LMDZ_WRFmeas/WRFV3/share/README.WRFmeas_lidar
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/README.WRFmeas_lidar	(revision 22)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/README.WRFmeas_lidar	(revision 22)
@@ -0,0 +1,317 @@
+WRF Model Version 3.3 (April 6, 2011)
+http://wrf-model.org/users/users.php
+
+------------------------
+WRF PUBLIC DOMAIN NOTICE
+------------------------
+
+WRF was developed at the National Center for Atmospheric Research
+(NCAR) which is operated by the University Corporation for
+Atmospheric Research (UCAR). NCAR and UCAR make no proprietary
+claims, either statutory or otherwise, to this version and
+release of WRF and consider WRF to be in the public domain for
+use by any person or entity for any purpose without any fee or
+charge. UCAR requests that any WRF user include this notice on
+any partial or full copies of WRF. WRF is provided on an "AS
+IS" basis and any warranties, either express or implied,
+including but not limited to implied warranties of
+non-infringement, originality, merchantability and fitness for a
+particular purpose, are disclaimed. In no event shall
+UCAR be liable for any damages, whatsoever, whether direct,
+indirect, consequential or special, that arise out of or in
+connection with the access, use or performance of WRF, including
+infringement actions.
+
+
+======================================
+This is the main directory for the WRF Version 3 source code release.
+======================================
+
+V3.3 Release Notes:
+-------------------
+
+Version 3.3 is released on April 6, 2011 (rev 4896).
+
+- For more information on WRF V3.3 release, visit WRF User's home pages
+  http://www.mmm.ucar.edu/wrf/users/, and
+  http://www.dtcenter.org/wrf-nmm/users/, and read the online User's Guide.
+
+
+======================================
+
+V3.2.1 Release Notes (8/18/10):
+------------------
+
+This is a bug fix release. The detailed updates for WRF-ARW can be found at
+http://www.mmm.ucar.edu/wrf/users/wrfv3.2/updates-3.2.1.html, and
+for WRF-NMM at
+http://www.dtcenter.org/wrf-nmm/users/model/wrfv3/updates.php
+
+See online User's Guides for latest information.
+
+======================================
+
+V3.2 Release Notes:
+-------------------
+
+Version 3.2 is released on March 31, 2010.
+
+- For more information on WRF V3.2 release, visit WRF User's home pages
+  http://www.mmm.ucar.edu/wrf/users/, and 
+  http://www.dtcenter.org/wrf-nmm/users/, and read the online User's Guide.
+
+
+======================================
+
+V3.1.1 Release Notes (7/31/09):
+------------------
+
+This is a bug fix release. The detailed updates for WRF-ARW can be found at
+http://www.mmm.ucar.edu/wrf/users/wrfv3.1/updates-3.1.1.html, and
+for WRF-NMM at
+http://www.dtcenter.org/wrf-nmm/users/model/wrfv3/updates.php
+
+See online User's Guides for latest information.
+
+======================================
+
+V3.1 Release Notes:
+-------------------
+
+Version 3.1 is released on April 9, 2009.
+
+- For more information on WRF V3.1 release, visit WRF Users home page
+  http://www.mmm.ucar.edu/wrf/users/, and read the online User's Guide.
+- WRF V3 executable will work with V3.0 wrfinput/wrfbdy (but requires 
+  inserting a new namelist variable use_baseparam_fr_nml in &dynamics). As
+  always, rerunning the new programs is recommended.
+- WRF V3.1 has changed a number of namelists. For example, ucmcall has
+  become sf_urban_physics, pd_moist(scalar/tke/chem) has become 
+  moist_adv_opt (scalar/tke/chem_adv_opt). Helpful messages will be printed
+  to aid you.
+
+
+V3.0.1.1 Release Notes:
+-----------------------
+
+Version 3.0.1.1 has only limited bug fixes compared to version 3.0.1.
+The detailed updates for WRF-ARW can be found at
+http://www.mmm.ucar.edu/wrf/users/wrfv3/updates-3.0.1.1.html.
+
+
+V3.0.1 Release Notes:
+---------------------
+
+This is a bug fix release. The detailed updates for WRF-ARW can be found at
+http://www.mmm.ucar.edu/wrf/users/wrfv3/updates-3.0.1.html.
+
+Online User's Guides have also been updated.
+
+======================================
+
+V3.0 Release Notes:
+------------------
+
+- For directions on compiling WRF, see below or Users Web page.
+- For more information on WRF V3 release, visit WRF Users home page
+  http://www.mmm.ucar.edu/wrf/users/
+- WRF V3 works with WPS, and SI is no longer supported. 
+  Please see User' Guide for WPS.
+- WRF V3 executable does not work with wrfinput/wrfbdy produced by previous versions.
+- WRF V3 has removed one obsolete namelist variable, dyn_opt. Edit with caution.
+
+======================================
+
+WRF update history:
+
+- V2.0.1: May 21, 2004
+- V2.0.2: June 3, 2004
+- V2.0.3: Nov 12, 2004 
+  V2.0.3.1: Dec 3, 2004
+- V2.1: August 4, 2005
+- V2.1.1: Nov 8, 2005
+- V2.1.2: Jan 27, 2006
+- V2.2: Dec 21, 2006
+- V2.2.1: Nov 1, 2007
+- V3.0: April 4, 2008
+- V3.0.1: August 5, 2008
+- V3.0.1.1: August 22, 2008
+- V3.1: April 9, 2009
+- V3.1.1: July 31, 2009
+- V3.2: March 31, 2010
+- V3.2.1: August 18, 2010
+- V3.3: April 6, 2011
+
+======================================
+
+How to compile and run?
+-----------------------
+
+- In WRFV3 directory, type 'configure' - this will create a configure.wrf
+  file that has appropriate compile options for the supported computers.
+
+  Note: WRF requires netCDF library. If your netCDF library is installed in
+        some odd directory, set environment variable NETCDF before you type
+        'configure'. For example, 
+
+        setenv NETCDF /usr/local/netcdf-pgi
+
+- Type 'compile case_name' where you can find the case_names by simply typing
+  compile. The following are available:
+
+      compile em_b_wave
+      compile em_grav2d_x
+      compile em_heldsuarez
+      compile em_hill2d_x
+      compile em_les
+      compile em_quarter_ss
+      compile em_real
+      compile em_seabreeze2d_x
+      compile em_squall2d_x
+      compile em_squall2d_y
+      compile em_scm_xy
+      compile em_tropical_cyclone
+
+- If sucessful, this will create either real.exe or ideal.exe and wrf.exe
+  in directory main/, and the appropriate executables will be linked into
+  the test directoires under test/case_name, or run/.
+
+- cd to the appropriate test or run directory to run ideal/real/wrf.
+
+- If it is one of the idealized cases (b_wave, hill2d_x, grav2d_x, quarter_ss,
+  squall2d_x, squall2d_y, em_les or em_heldsuarez), cd the the appropriate directory, type
+
+  ./ideal.exe 
+
+  to produce wrfinput_d01 file for wrf model. Then type
+
+  ./wrf.exe 
+
+  to run. 
+
+- If it is real-data case (real), place files from WPS (met_em.*)
+  in the appropriate directory, type
+
+  ./real.exe
+
+  to produce wrfbdy_d01 and wrfinput_d01. Then type
+
+  ./wrf.exe
+
+  to run.
+
+- If you use mpich, type
+
+  mpirun -np number-of-processors wrf.exe
+
+- For information on how to make nested runs, visit
+  http://www.mmm.ucar.edu/wrf/users/
+
+
+======================================
+
+What is in WRF V3?
+
+- Advanced Research WRF (ARW) solver: Eulerian mass, hydrostatic and non-hydrostatic
+  * Arakawa C-grid staggering
+  * Runge-Kutta 2nd and 3rd order timestep options
+  * scalar-conserving flux form for prognostic variables
+  * 2nd to 6th order advection options (horizontal and vertical)
+  * time-split small step for acoustic modes
+  * small step horizontally explicit, vertically implicit
+  * divergence damping option and vertical time off-centering
+  * external-mode filtering option for mass model
+  * hydrostatic option via namelist option
+  * positive-definite and monotonic advection for moisture, scalar, tke and chemical tracers
+  * global modeling capability on latitude-longitude grid
+  * digital filter initialization
+
+- Two-way nesting:
+  * multiple domains and multiple nest levels
+  * supports integer nest grid ratio
+  * feedback option for both odd and even nest grid ratios 
+  * smoothing options
+
+- One-way nesting 
+
+- Moving nest 
+  * Specified move
+  * Automatic move using a mid-level vortex-following algorithm
+
+- Physics options:
+  * microphysics (Kessler/ WRF Single Moment 3, 5 and 6 classes / Lin et al./ Mibrandt 2-moment /
+    Eta Ferrier / Thompson / Goddard / 2-moment Morrison / WRF Double Moment 5 and 6 classes /
+    SBU-Lin 5-classes )
+  * cumulus parameterization (Kain-Fritsch with shallow convection / 
+    Betts-Miller-Janjic / Grell-Devenyi ensemble / New Grell (with shallow convection option) /
+    Tiedtke (with shallow conv and momentum transport) / NSAS (with shallow conv and momentum transport) /
+    SAS (with shallow conv for ARW) / Zhang-McFarlane (with momentum transport) )
+  * UW shallow convection
+  * planetary boundary layer (Yosei University / Mellor-Yamada-Janjic / ACM2 / QNSE / MYNN / BouLac /
+    UW / TEMF )
+  * slab soil model (5-layer thermal diffusion / Noah land-surface model (4 levels) /
+    RUC LSM (6 levels) / Pleim-Xu (2 levels)
+  * Urban canopy model, BEP multi-layer and BEM  (works with Noah LSM, BEP and 
+    BEM also requires MYJ or BouLac)
+  * longwave radiation (RRTM / CAM / RRTMG / new Goddard )
+  * shortwave radiation (simple MM5 scheme / old Goddard / CAM / RRTMG / new Goddard), 
+    terrain slope and shading effect
+  * sub-grid turbulence (constant K diffusion/ 2-D Smagorinsky/ predicted TKE /
+    2-D, 6th order diffusion / Nonlinear Backscatter Anisotropic (NBA) sub-grid turbulence stress for LES )
+  * Rayleigh damping for w at the upper boundary layer
+  * gravity wave drag
+  * land-use categories determine surface properties; support for 24 category USGS and 
+    20 category MODIS (Noah LSM only in V3.1)
+  * Options for modifying SST, sea ice, vegetation fraction, albedo, and deep soil temp
+    for long simulations
+  * fractional sea ice option for polar regions; modified snow/ice physics
+  * single-column ocean mixed layer model
+  * drag and enthalpy flux formulation for hurricane applications
+
+- Nudging:
+  * three-dimensional and surface analysis nudging 
+  * observation nudging 
+  * spectral nudging
+
+- Software
+
+  * Hierarchical software architecture that insulates scientific code
+    (Model Layer) from computer architecture (Driver Layer)
+  * Multi-level parallelism supporting shared-memory (OpenMP), distributed-memory (MPI), 
+    and hybrid share/distributed modes of execution
+  * Active data registry: defines and manages model state fields, I/O,
+    nesting, configuration, and numerous other aspects of WRF through a single file, 
+    called the Registry
+  * Two-way nesting:
+      Easy to extend: forcing and feedback of new fields specified by
+        editing a single table in the Registry
+      Efficient: 5-8% overhead on 64 processes of IBM
+      Moving nests.
+  * Enhanced I/O options:
+      NetCDF and Parallel HDF5 formats
+      Nine auxiliary input and history output streams separately controllable through the namelist
+      Special input streams for analysis and observation nudging
+      Output file names and time-stamps specifiable through namelist
+      Special output stream for 3DVAR
+  * Efficient execution on a range of computing platforms:
+      IBM SP systems, (e.g. NCAR "bluevista","blueice","bluefire" Power5-based system)
+      IBM Blue Gene
+      SGI Origin and Altix
+      Linux/Intel
+         IA64 MPP (HP Superdome, SGI Altix, NCSA Teragrid systems)
+         IA64 SMP
+         x86_64 (e.g. TACC's "Ranger", NOAA/GSD "wJet" )
+	 PGI, Intel, Pathscale, gfortran, g95 compilers supported
+      Sun Solaris (single threaded and SMP)
+      Cray X1, X1e (vector), XT3/4 (Opteron)
+      Mac Intel/ppc, PGI/ifort/g95
+      NEC SX/8
+      HP-UX
+      Fujitsu VPP 5000
+  * RSL_LITE: communication layer, scalable to very large domains, supports nesting.
+  * I/O: NetCDF, parallel NetCDF (Argonne), HDF5, GRIB, raw binary, Quilting (asynchronous I/O), MCEL (coupling)
+  * ESMF Time Management, including exact arithmetic for fractional
+    time steps (no drift).
+  * ESMF integration - WRF can be run as an ESMF component.
+  * Improved documentation, both on-line (web based browsing tools) and in-line
+
Index: /branches/LMDZ_WRFmeas/WRFV3/share/input_wrf.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/input_wrf.F	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/input_wrf.F	(revision 22)
@@ -1287,4 +1287,8 @@
 #if (DA_CORE != 1)
     CALL wrf_tsin( grid , ierr )
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+    CALL wrf_lidarin( grid , ierr )
+#endif
 #endif
 
Index: /branches/LMDZ_WRFmeas/WRFV3/share/mediation_integrate.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/mediation_integrate.F	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/mediation_integrate.F	(revision 22)
@@ -171,4 +171,8 @@
    ! Compute time series variables
    CALL calc_ts(grid)
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+    CALL calc_lidar(grid)
+#endif
 
    RETURN
@@ -349,4 +353,13 @@
         INTEGER, INTENT(INOUT) :: ierr
      END SUBROUTINE wrf_tsin
+
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+    SUBROUTINE wrf_lidarin ( grid , ierr )
+        USE module_domain
+        TYPE ( domain ), INTENT(INOUT) :: grid
+        INTEGER, INTENT(INOUT) :: ierr
+     END SUBROUTINE wrf_lidarin
+#endif
 
    END INTERFACE
@@ -467,4 +480,8 @@
        !   from the tslist file
        CALL wrf_tsin( nest , ierr )
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+       CALL wrf_lidarin( nest , ierr )
+#endif
 #endif
      END IF
@@ -1159,4 +1176,8 @@
    ! Write out time series
    CALL write_ts( grid )
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+   CALL write_lidar( grid )
+#endif
 
    RETURN
Index: /branches/LMDZ_WRFmeas/WRFV3/share/start_domain.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/start_domain.F	(revision 21)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/start_domain.F	(revision 22)
@@ -27,4 +27,11 @@
          TYPE (domain) :: grid
       END SUBROUTINE calc_ts_locations
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+      SUBROUTINE calc_lidar_locations( grid )
+         USE module_domain
+         TYPE (domain) :: grid
+      END SUBROUTINE calc_lidar_locations
+#endif
 #endif
    END INTERFACE
@@ -56,4 +63,8 @@
 #if (DA_CORE != 1)
    CALL calc_ts_locations( grid )
+#if WRFMEAS
+! L. Fita, LMD. May 2014
+   CALL calc_lidar_locations( grid )
+#endif
 #endif
 
Index: /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidar.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidar.F	(revision 22)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidar.F	(revision 22)
@@ -0,0 +1,530 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! L. Fita, LMD. May 2014. Module to compute lidar retrievals (values on the vertical 
+!   column) at each time-step. Based on wrf_ts.F
+! This routine prints out the current value of variables at all specified
+!   time series locations that are within the current patch.
+!
+! Michael G. Duda -- 25 August 2005
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+SUBROUTINE calc_lidar_locations( grid )
+
+   USE module_domain, ONLY : domain, get_ijk_from_grid, domain_clock_get
+   USE module_configure, ONLY : model_config_rec, grid_config_rec_type, model_to_grid_config_rec
+   USE module_dm, ONLY : wrf_dm_min_real
+   USE module_llxy
+   USE module_state_description
+
+   IMPLICIT NONE
+
+   ! Arguments
+   TYPE (domain), INTENT(INOUT) :: grid
+
+   ! Externals
+   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+! L. Fita, LMD. May 2014
+! Already defined as function by wrf_timeseries.F
+   INTEGER, EXTERNAL :: get_unused_unit
+
+   ! Local variables
+   INTEGER :: nlidarloc_temp
+   INTEGER :: i, k, iunit
+   REAL :: lidar_rx, lidar_ry, lidar_xlat, lidar_xlong, lidar_hgt
+   REAL :: known_lat, known_lon
+   CHARACTER (LEN=132) :: message
+   TYPE (PROJ_INFO) :: lidar_proj
+   TYPE (grid_config_rec_type) :: config_flags
+
+   INTEGER :: ids, ide, jds, jde, kds, kde,        &
+              ims, ime, jms, jme, kms, kme,        &
+              ips, ipe, jps, jpe, kps, kpe,        &
+              imsx, imex, jmsx, jmex, kmsx, kmex,  &
+              ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
+              imsy, imey, jmsy, jmey, kmsy, kmey,  &
+              ipsy, ipey, jpsy, jpey, kpsy, kpey
+   CHARACTER (LEN=50)                                    :: SimStartTime
+
+   IF ( grid%nlidarloc .LE. 0 ) RETURN
+
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+   IF ( grid%dfi_stage == DFI_FST ) THEN
+#endif
+      CALL get_ijk_from_grid ( grid ,                               &
+                               ids, ide, jds, jde, kds, kde,        &
+                               ims, ime, jms, jme, kms, kme,        &
+                               ips, ipe, jps, jpe, kps, kpe,        &
+                               imsx, imex, jmsx, jmex, kmsx, kmex,  &
+                               ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
+                               imsy, imey, jmsy, jmey, kmsy, kmey,  &
+                               ipsy, ipey, jpsy, jpey, kpsy, kpey )
+   
+      CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
+   
+      ! Set up map transformation structure
+      CALL map_init(lidar_proj)
+   
+      IF (ips <= 1 .AND. 1 <= ipe .AND. &
+          jps <= 1 .AND. 1 <= jpe) THEN
+         known_lat = grid%xlat(1,1)
+         known_lon = grid%xlong(1,1)
+      ELSE
+         known_lat = 9999.
+         known_lon = 9999.
+      END IF
+      known_lat = wrf_dm_min_real(known_lat)
+      known_lon = wrf_dm_min_real(known_lon)
+   
+      ! Mercator
+      IF (config_flags%map_proj == PROJ_MERC) THEN
+         CALL map_set(PROJ_MERC, lidar_proj,               &
+                      truelat1 = config_flags%truelat1, &
+                      lat1     = known_lat,             &
+                      lon1     = known_lon,             &
+                      knowni   = 1.,                    &
+                      knownj   = 1.,                    &
+                      dx       = config_flags%dx)
+   
+      ! Lambert conformal
+      ELSE IF (config_flags%map_proj == PROJ_LC) THEN
+         CALL map_set(PROJ_LC, lidar_proj,                  &
+                      truelat1 = config_flags%truelat1,  &
+                      truelat2 = config_flags%truelat2,  &
+                      stdlon   = config_flags%stand_lon, &
+                      lat1     = known_lat,              &
+                      lon1     = known_lon,              &
+                      knowni   = 1.,                     &
+                      knownj   = 1.,                     &
+                      dx       = config_flags%dx)
+   
+      ! Polar stereographic
+      ELSE IF (config_flags%map_proj == PROJ_PS) THEN
+         CALL map_set(PROJ_PS, lidar_proj,                  &
+                      truelat1 = config_flags%truelat1,  &
+                      stdlon   = config_flags%stand_lon, &
+                      lat1     = known_lat,              &
+                      lon1     = known_lon,              &
+                      knowni   = 1.,                     &
+                      knownj   = 1.,                     &
+                      dx       = config_flags%dx)
+   
+#if (EM_CORE == 1)
+      ! Cassini (global ARW)
+      ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
+         CALL map_set(PROJ_CASSINI, lidar_proj,                            &
+                      latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
+                      loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), & 
+                      lat1     = known_lat,                             &
+                      lon1     = known_lon,                             &
+                      lat0     = config_flags%pole_lat,                 &
+                      lon0     = config_flags%pole_lon,                 &
+                      knowni   = 1.,                                    &
+                      knownj   = 1.,                                    &
+                      stdlon   = config_flags%stand_lon)
+#endif
+
+      ! Rotated latitude-longitude
+      ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
+         CALL map_set(PROJ_ROTLL, lidar_proj,                      &
+! I have no idea how this should work for NMM nested domains
+                      ixdim    = grid%e_we-1,                   &
+                      jydim    = grid%e_sn-1,                   &
+                      phi      = real(grid%e_sn-2)*grid%dy/2.0, &
+                      lambda   = real(grid%e_we-2)*grid%dx,     &
+                      lat1     = config_flags%cen_lat,          &
+                      lon1     = config_flags%cen_lon,          &
+                      latinc   = grid%dy,                       &
+                      loninc   = grid%dx,                       &
+                      stagger  = HH)
+   
+      END IF
+   
+      ! Determine lidar locations for domain
+      IF (.NOT. grid%have_calculated_lidarlocs) THEN
+         grid%have_calculated_lidarlocs = .TRUE.
+         WRITE(message, '(A43,I3)') 'Computing lidar locations for domain ', grid%id
+         CALL wrf_message(message)
+   
+         nlidarloc_temp = 0
+         DO k=1,grid%nlidarloc
+   
+            IF (config_flags%map_proj == 0) THEN  ! For idealized cases, no map transformation needed
+               lidar_rx = grid%latlidarloc(k)           ! NB: (x,y) = (lat,lon) rather than (x,y) = (lon,lat)
+               lidar_ry = grid%lonlidarloc(k)
+            ELSE
+               CALL latlon_to_ij(lidar_proj, grid%latlidarloc(k), grid%lonlidarloc(k), lidar_rx, lidar_ry)
+            END IF
+            
+
+            nlidarloc_temp = nlidarloc_temp + 1
+            grid%ilidarloc(nlidarloc_temp) = NINT(lidar_rx)
+            grid%jlidarloc(nlidarloc_temp) = NINT(lidar_ry)
+            grid%id_lidarloc(nlidarloc_temp) = k
+   
+            ! Is point outside of domain (or on the edge of domain)?
+            IF (grid%ilidarloc(nlidarloc_temp) < ids .OR. grid%ilidarloc(nlidarloc_temp) > ide .OR. &
+                grid%jlidarloc(nlidarloc_temp) < jds .OR. grid%jlidarloc(nlidarloc_temp) > jde) THEN
+               nlidarloc_temp = nlidarloc_temp - 1
+   
+            END IF
+   
+         END DO
+   
+         grid%next_lidar_time = 1
+   
+         grid%nlidarloc_domain = nlidarloc_temp
+   
+         DO k=1,grid%nlidarloc_domain
+   
+            ! If location is outside of patch, we need to get lat/lon of TS grid cell from another patch
+            IF (grid%ilidarloc(k) < ips .OR. grid%ilidarloc(k) > ipe .OR. &
+                grid%jlidarloc(k) < jps .OR. grid%jlidarloc(k) > jpe) THEN
+               lidar_xlat  = 1.E30
+               lidar_xlong = 1.E30
+               lidar_hgt   = 1.E30
+            ELSE
+               lidar_xlat  = grid%xlat(grid%ilidarloc(k),grid%jlidarloc(k))
+               lidar_xlong = grid%xlong(grid%ilidarloc(k),grid%jlidarloc(k))
+#if (EM_CORE == 1)
+               lidar_hgt   = grid%ht(grid%ilidarloc(k),grid%jlidarloc(k))
+#endif
+            END IF
+#if DM_PARALLEL
+            lidar_xlat  = wrf_dm_min_real(lidar_xlat)
+            lidar_xlong = wrf_dm_min_real(lidar_xlong)
+            lidar_hgt   = wrf_dm_min_real(lidar_hgt)
+#endif
+   
+            IF ( wrf_dm_on_monitor() ) THEN
+
+               iunit = get_unused_unit()
+               IF ( iunit <= 0 ) THEN
+                  CALL wrf_error_fatal('Error in calc_lidar_locations: could not find a free Fortran unit.')
+               END IF
+               CALL domain_clock_get(grid, simulationStartTimeStr = SimStartTime)
+
+               WRITE(grid%lidar_filename(k),'(A)') TRIM(grid%namelidarloc(grid%id_lidarloc(k)))//'.LIDAR.d00'
+               i = LEN_TRIM(grid%lidar_filename(k))
+               WRITE(grid%lidar_filename(k)(i-1:i),'(I2.2)') grid%id
+               OPEN(UNIT=iunit, FILE=TRIM(grid%lidar_filename(k)), FORM='FORMATTED', STATUS='REPLACE')
+#if (EM_CORE == 1)
+               WRITE(UNIT=iunit, &
+                     FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A32,A19)') &
+                     grid%desclidarloc(grid%id_lidarloc(k))//' ', grid%id, grid%id_lidarloc(k), &
+                     ' '//grid%namelidarloc(grid%id_lidarloc(k)), &
+                     ' (', grid%latlidarloc(grid%id_lidarloc(k)), ',', grid%lonlidarloc(grid%id_lidarloc(k)), ') (', &
+                     grid%ilidarloc(k), ',', grid%jlidarloc(k), ') (', &
+                     lidar_xlat, ',', lidar_xlong, ') ', &
+                     lidar_hgt,' meters. simulation start time: ',TRIM(SimStartTime)
+#else
+               WRITE(UNIT=iunit, &
+                     FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A32,A19)') &
+                     grid%desclidarloc(grid%id_lidarloc(k))//' ', grid%id, grid%id_lidarloc(k), &
+                     ' '//grid%namelidarloc(grid%id_lidarloc(k)), &
+                     ' (', grid%latlidarloc(grid%id_lidarloc(k)), ',', grid%lonlidarloc(grid%id_lidarloc(k)), ') (', &
+                     grid%ilidarloc(k), ',', grid%jlidarloc(k), ') (', &
+                     lidar_xlat, ',', lidar_xlong, ') ', &
+                     lidar_hgt,' meters. simulation start time: ',TRIM(SimStartTime)
+#endif
+               CLOSE(UNIT=iunit)
+            END IF
+         END DO
+   
+      END IF
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+   END IF
+#endif
+
+END SUBROUTINE calc_lidar_locations
+
+
+SUBROUTINE calc_lidar( grid )
+
+   USE module_domain
+   USE module_model_constants
+
+   IMPLICIT NONE
+
+   ! Arguments
+   TYPE (domain), INTENT(INOUT) :: grid
+
+   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+
+   ! Local variables
+   INTEGER :: i, k, mm, n, ix, iy, rc
+   REAL :: xtime_minutes
+   REAL, ALLOCATABLE, DIMENSION(:)                       :: earth_u, earth_v,        &
+     output_t, output_qv, output_qc, output_qr, output_qs, output_qh, output_qi,     &
+     output_qg
+   REAL, ALLOCATABLE, DIMENSION(:) :: p8w
+
+   ! Parameter lidar_model_level:  
+       ! TRUE to output T, Q, and wind at lowest model level
+       ! FALSE to output T and Q at 2-m and wind at 10-m diagnostic levels:
+   LOGICAL, PARAMETER :: lidar_model_level = .FALSE.  
+
+   IF ( grid%nlidarloc_domain .LE. 0 ) RETURN
+
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+   IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
+#endif
+
+   n = grid%next_lidar_time
+
+   ALLOCATE(p8w(grid%sm32:grid%em32))
+   ALLOCATE(earth_u(grid%sm32:grid%em32))
+   ALLOCATE(earth_v(grid%sm32:grid%em32))
+   ALLOCATE(output_t(grid%sm32:grid%em32))
+   ALLOCATE(output_qv(grid%sm32:grid%em32))
+   ALLOCATE(output_qc(grid%sm32:grid%em32))
+   ALLOCATE(output_qr(grid%sm32:grid%em32))
+   ALLOCATE(output_qs(grid%sm32:grid%em32))
+   ALLOCATE(output_qh(grid%sm32:grid%em32))
+   ALLOCATE(output_qi(grid%sm32:grid%em32))
+   ALLOCATE(output_qg(grid%sm32:grid%em32))
+
+   DO i=1,grid%nlidarloc_domain
+
+      ix = grid%ilidarloc(i)
+      iy = grid%jlidarloc(i)
+  
+      IF (grid%sp31 <= ix .AND. ix <= grid%ep31 .AND. &
+          grid%sp33 <= iy .AND. iy <= grid%ep33) THEN
+          DO k=grid%sm32, grid%em32-1
+            !
+            ! Output from the lowest model computational level:
+            !
+#if (EM_CORE == 1)
+            earth_u(k) = grid%u_2(ix,k,iy)*grid%cosa(ix,iy)-grid%v_2(ix,k,iy)*grid%sina(ix,iy)
+            earth_v(k) = grid%v_2(ix,k,iy)*grid%cosa(ix,iy)+grid%u_2(ix,k,iy)*grid%sina(ix,iy)
+            output_t(k) = grid%t_2(ix,k,iy) + 300.
+#else
+            earth_u(k) = grid%u_2(ix,k,iy)*grid%cosa(ix,iy)-grid%v_2(ix,k,iy)*grid%sina(ix,iy)
+            earth_v(k) = grid%v_2(ix,k,iy)*grid%cosa(ix,iy)+grid%u_2(ix,k,iy)*grid%sina(ix,iy)
+            output_t(k) = grid%t(ix,k,iy) + 300.
+#endif
+            output_qv(k) = grid%moist(ix,k,iy,P_QV)
+            output_qc(k) = grid%moist(ix,k,iy,P_QC)
+            output_qr(k) = grid%moist(ix,k,iy,P_QR)
+            output_qs(k) = grid%moist(ix,k,iy,P_QS)
+            output_qh(k) = grid%moist(ix,k,iy,P_QH)
+            output_qi(k) = grid%moist(ix,k,iy,P_QI)
+            output_qg(k) = grid%moist(ix,k,iy,P_QG)
+          END DO
+
+      
+         CALL domain_clock_get( grid, minutesSinceSimulationStart=xtime_minutes )
+         grid%lidar_hour(n,i) = xtime_minutes / 60.
+         grid%lidar_z(n,:,i) = ( grid%ph_2(ix,:,iy) + grid%phb(ix,:,iy) ) / g
+         grid%lidar_p(n,:,i) = grid%p(ix,:,iy) + grid%pb(ix,:,iy) * 0.01
+         grid%lidar_u(n,:,i) = earth_u
+         grid%lidar_v(n,:,i) = earth_v
+         grid%lidar_w(n,:,i) = grid%w_2(ix,:,iy)
+         grid%lidar_t(n,:,i) = output_t
+         grid%lidar_qv(n,:,i) = output_qv
+         grid%lidar_qc(n,:,i) = output_qc
+         grid%lidar_qr(n,:,i) = output_qr
+         grid%lidar_qs(n,:,i) = output_qs
+         grid%lidar_qh(n,:,i) = output_qh
+         grid%lidar_qi(n,:,i) = output_qi
+         grid%lidar_qg(n,:,i) = output_qg
+         grid%lidar_dens(n,:,i) = 1./grid%alt(ix,:,iy)
+         grid%lidar_cldfra(n,:,i) = grid%cldfra(ix,:,iy)
+         grid%lidar_drydens(n,i) = grid%mu_2(ix,iy) + grid%mub(ix,iy) * 0.01
+         grid%lidar_psfc(n,i) = grid%psfc(ix,iy)
+#if (EM_CORE == 1)
+         grid%lidar_rainc(n,i)  = grid%rainc(ix,iy)
+         grid%lidar_rainnc(n,i) = grid%rainnc(ix,iy)
+#endif   
+      ELSE
+         grid%lidar_hour(n,i) = 1.E30
+         grid%lidar_u(n,:,i) = 1.E30
+         grid%lidar_v(n,:,i) = 1.E30
+         grid%lidar_w(n,:,i) = 1.E30
+         grid%lidar_t(n,:,i) = 1.E30
+         grid%lidar_qv(n,:,i) = 1.E30
+         grid%lidar_qc(n,:,i) = 1.E30
+         grid%lidar_qr(n,:,i) = 1.E30
+         grid%lidar_qs(n,:,i) = 1.E30
+         grid%lidar_qh(n,:,i) = 1.E30
+         grid%lidar_qi(n,:,i) = 1.E30
+         grid%lidar_qg(n,:,i) = 1.E30
+         grid%lidar_dens(n,:,i) = 1.E30
+         grid%lidar_cldfra(n,:,i) = 1.E30
+         grid%lidar_psfc(n,i) = 1.E30
+#if (EM_CORE == 1)
+         grid%lidar_rainc(n,i)  = 1.E30
+         grid%lidar_rainnc(n,i) = 1.E30
+#endif
+   
+      END IF
+   END DO
+
+   DEALLOCATE(p8w)
+ 
+   grid%next_lidar_time = grid%next_lidar_time + 1
+
+   IF ( grid%next_lidar_time > grid%lidar_buf_size ) CALL write_lidar(grid)
+
+END SUBROUTINE calc_lidar
+
+
+SUBROUTINE write_lidar( grid )
+
+   USE module_domain, ONLY : domain
+   USE module_dm, ONLY : wrf_dm_min_reals
+   USE module_state_description
+
+   IMPLICIT NONE
+
+   ! Arguments
+   TYPE (domain), INTENT(INOUT) :: grid
+
+   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+! L. Fita, LMD. May 2014
+! Already defined as function by wrf_timeseries.F
+   INTEGER, EXTERNAL :: get_unused_unit
+
+   ! Local variables
+   INTEGER :: i, n, k, ix, iy, iunit
+   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: lidar_buf
+
+   IF ( grid%nlidarloc_domain .LE. 0 ) RETURN
+
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+   IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
+#endif
+
+#ifdef DM_PARALLEL
+   ALLOCATE(lidar_buf(grid%lidar_buf_size,grid%sm32:grid%em32,grid%max_lidar_locs))
+
+   lidar_buf(:,1,:) = grid%lidar_hour(:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,1,:),grid%lidar_hour(:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_z(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_z(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_p(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_p(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_u(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_u(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_v(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_v(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_w(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_w(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_t(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_t(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qv(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qv(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qc(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qc(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qr(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qr(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qs(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qs(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qh(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qh(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qi(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qi(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_qg(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_qg(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_dens(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_dens(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,:,:) = grid%lidar_cldfra(:,:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,:,:),grid%lidar_cldfra(:,:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,1,:) = grid%lidar_drydens(:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,1,:),grid%lidar_drydens(:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,1,:) = grid%lidar_psfc(:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,1,:),grid%lidar_psfc(:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+#if (EM_CORE == 1)
+   lidar_buf(:,1,:) = grid%lidar_rainc(:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,1,:),grid%lidar_rainc(:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+
+   lidar_buf(:,1,:) = grid%lidar_rainnc(:,:)
+   CALL wrf_dm_min_reals(lidar_buf(:,1,:),grid%lidar_rainnc(:,:),grid%lidar_buf_size*grid%max_lidar_locs)
+#endif
+
+   DEALLOCATE(lidar_buf)
+#endif
+
+   IF ( wrf_dm_on_monitor() ) THEN
+
+      iunit = get_unused_unit()
+      IF ( iunit <= 0 ) THEN
+         CALL wrf_error_fatal('Error in write_lidar: could not find a free Fortran unit.')
+      END IF
+
+      DO i=1,grid%nlidarloc_domain
+
+         ix = grid%ilidarloc(i)
+         iy = grid%jlidarloc(i)
+
+         OPEN(UNIT=iunit, FILE=TRIM(grid%lidar_filename(i)), STATUS='unknown', POSITION='append', FORM='formatted')
+
+         DO n=1,grid%next_lidar_time - 1
+#if (EM_CORE == 1)
+            WRITE(UNIT=iunit,FMT='(a8,1x,i2,f13.6,i5,i5,i5,1x,4(e13.5,1x))') &
+                            'new_time',grid%id, grid%lidar_hour(n,i), &
+                            grid%id_lidarloc(i), ix, iy, &
+                            grid%lidar_psfc(n,i), &
+                            grid%lidar_rainc(n,i), &
+                            grid%lidar_rainnc(n,i), &
+                            grid%lidar_drydens(n,i)
+#else
+            WRITE(UNIT=iunit,FMT='(a8,1x,i2,f13.6,i5,i5,i5,1x,2(e13.5,1x))') &
+                            'new_time',grid%id, grid%lidar_hour(n,i), &
+                            grid%id_lidarloc(i), ix, iy, &
+                            grid%lidar_psfc(n,i), &
+                            grid%lidar_drydens(n,i)
+#endif
+            WRITE(UNIT=iunit, FMT='(5x,a3,1x,15(a14,1x),a10)') 'k', 'z [m]', 'p [hPa]',   &
+              'u [ms-1]', 'v [ms-1]', 'w [ms-1]', 't_pot [k]', 'qv [kgkg-1]',             &
+              'qc [kgkg-1]', 'qr [kgkg-1]', 'qs [kgkg-1]', 'qh [kgkg-1]', 'qi [kgkg-1]',  &
+              'qg [kgkg-1]', 'dens [kg m-3]', 'cldfra [1]', '__________'
+            DO k=grid%sm32,grid%em32-1
+              WRITE(UNIT=iunit,FMT='(5x,i3,1x,15(e13.5,1x))') k, &
+                              grid%lidar_z(n,k,i), &
+                              grid%lidar_p(n,k,i), &
+                              grid%lidar_u(n,k,i), &
+                              grid%lidar_v(n,k,i), &
+                              grid%lidar_w(n,k,i), &
+                              grid%lidar_t(n,k,i), &
+                              grid%lidar_qv(n,k,i), &
+                              grid%lidar_qc(n,k,i), &
+                              grid%lidar_qr(n,k,i), &
+                              grid%lidar_qs(n,k,i), &
+                              grid%lidar_qh(n,k,i), &
+                              grid%lidar_qi(n,k,i), &
+                              grid%lidar_qg(n,k,i), &
+                              grid%lidar_dens(n,k,i), &
+                              grid%lidar_cldfra(n,k,i)
+           END DO
+         END DO
+
+         CLOSE(UNIT=iunit)
+
+      END DO
+
+   END IF
+
+   grid%next_lidar_time = 1
+
+END SUBROUTINE write_lidar
+
+! L. Fita, LMD. May 2014
+! Already defined as function by wrf_timeseries.F
+! No SUBROUTINE calc_p8w(grid, ix, iy, p8w, k_start, k_end)
Index: /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidarin.F
===================================================================
--- /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidarin.F	(revision 22)
+++ /branches/LMDZ_WRFmeas/WRFV3/share/wrf_lidarin.F	(revision 22)
@@ -0,0 +1,117 @@
+SUBROUTINE wrf_lidarin ( grid , ierr )
+
+    USE module_domain
+    USE module_utility
+
+    IMPLICIT NONE
+
+#include <wrf_io_flags.h>
+#include <wrf_status_codes.h>
+
+    TYPE(domain), INTENT(INOUT) :: grid
+    INTEGER, INTENT(INOUT) :: ierr
+
+    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+    INTEGER, EXTERNAL :: get_unused_unit
+
+    REAL, ALLOCATABLE, DIMENSION(:) :: latlidarlocs, lonlidarlocs
+    INTEGER :: istatus, iunit
+    LOGICAL :: exists
+    CHARACTER (LEN=256) :: errmess
+
+    ierr = 0
+
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+    IF ( grid%dfi_opt == DFI_NODFI .OR. (grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage == DFI_SETUP) ) THEN
+#endif
+
+       grid%nlidarloc = 0
+       grid%have_calculated_lidarlocs = .FALSE.
+   
+       IF ( grid%max_lidar_locs <= 0 ) RETURN
+   
+       IF ( wrf_dm_on_monitor() ) THEN
+   
+          INQUIRE(FILE='lidarlist', EXIST=exists)
+
+          IF (exists) THEN
+
+             iunit = get_unused_unit()
+             IF ( iunit <= 0 ) THEN
+                CALL wrf_error_fatal('Error in wrf_lidarin: could not find a free Fortran unit.')
+             END IF
+
+             ! Input time series locations
+             OPEN(UNIT=iunit, FILE='lidarlist', FORM='formatted', STATUS='old', IOSTAT=istatus)
+   
+             IF (istatus == 0) THEN
+   
+                ! Ignore first three lines, which constitute a header
+                READ(UNIT=iunit, FMT='(1X)')
+                READ(UNIT=iunit, FMT='(1X)')
+                READ(UNIT=iunit, FMT='(1X)')
+   
+                ! Read in time series locations
+                istatus = 0
+                DO WHILE (istatus == 0)
+                   READ(UNIT=iunit, FMT='(A25,1X,A5,1X,F7.3,1X,F8.3)', IOSTAT=istatus)            &
+                        grid%desclidarloc(grid%nlidarloc+1), grid%namelidarloc(grid%nlidarloc+1), &
+                        grid%latlidarloc(grid%nlidarloc+1), grid%lonlidarloc(grid%nlidarloc+1)
+                   IF (istatus == 0) grid%nlidarloc = grid%nlidarloc + 1
+                   IF (istatus > 0) THEN
+                      WRITE(errmess, FMT='(I4)') grid%nlidarloc + 3   ! Three extra for the header of the file
+                      CALL wrf_message('Error in lidarlist, line '//trim(errmess))
+                      EXIT    ! (technically unecessary, as we will exit the loop anyway)
+                   END IF
+                   IF ( grid%nlidarloc >= grid%max_lidar_locs ) THEN
+                      IF ( istatus == 0 ) THEN                 ! Assume there were more lines in the file
+                         WRITE(errmess, FMT='(A,I4,A)') 'Ignoring all time series locations beyond #', &
+                                                    grid%nlidarloc,'. Increase max_lidar_locs in namelist.input'
+                         CALL wrf_message(trim(errmess))
+                      END IF
+                      EXIT
+                   END IF
+                END DO
+      
+                CLOSE(iunit)
+   
+             END IF
+
+          END IF  ! lidarlist file exists
+   
+       END IF
+
+#ifdef DM_PARALLEL
+       CALL wrf_dm_bcast_integer(grid%nlidarloc, 1)
+       CALL wrf_dm_bcast_real(grid%latlidarloc, grid%max_lidar_locs)
+       CALL wrf_dm_bcast_real(grid%lonlidarloc, grid%max_lidar_locs)
+#endif
+#if ((EM_CORE == 1) && (DA_CORE != 1))
+    END IF
+#endif
+
+END SUBROUTINE wrf_lidarin
+
+
+! L. Fita, LMD. May 2014
+! Already defined as function by wrf_timeseries.F
+!
+!!INTEGER FUNCTION get_unused_unit()
+!!
+!!    IMPLICIT NONE
+!!
+!!    INTEGER, PARAMETER :: min_unit_number = 30
+!!    INTEGER, PARAMETER :: max_unit_number = 99 
+!!
+!!    LOGICAL :: opened
+!!
+!!    DO get_unused_unit = min_unit_number, max_unit_number
+!!       INQUIRE(UNIT=get_unused_unit, OPENED=opened)
+!!       IF ( .NOT. opened ) RETURN
+!!    END DO
+!!
+!!    get_unused_unit = -1
+!!
+!!    RETURN
+!!    
+!!END FUNCTION get_unused_unit
