Index: palm/trunk/UTIL/inifor/src/inifor.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3183)
+++ palm/trunk/UTIL/inifor/src/inifor.f90	(revision 3395)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Added main loop support for computation of geostrophic winds and surface
+!     pressure
+! Cleaned up terminal output, show some meVssages only with --debug option
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching
 ! Renamend initial-condition mode variable 'mode' to 'ic_mode'
@@ -56,10 +61,12 @@
                fini_file_lists, preprocess, origin_lon, origin_lat,            &
                output_file, io_group_list, output_var_table,                   &
-               cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,             &
-               average_imin, average_imax, average_jmin, average_jmax
-
+               cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3,                 &
+               averaging_width_ns, averaging_width_ew, phi_n, lambda_n,        &
+               lam_centre, phi_centre
     USE io
     USE transform,                                                             &
-        ONLY:  average_profile, average_2d, interpolate_2d, interpolate_3d
+        ONLY:  average_profile, interpolate_2d, interpolate_3d,                &
+               geostrophic_winds, extrapolate_density, extrapolate_pressure,   &
+               get_surface_pressure
     USE types
     
@@ -70,7 +77,22 @@
     INTEGER                                 ::  iter
     REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) ::  output_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_arr
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east
+    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west
+    REAL(dp), POINTER, DIMENSION(:)         ::  internal_arr
+    REAL(dp), POINTER, DIMENSION(:)         ::  ug_vg_arr
     TYPE(nc_var), POINTER                   ::  output_var
     TYPE(io_group), POINTER                 ::  group
     TYPE(container), ALLOCATABLE            ::  input_buffer(:)
+    LOGICAL, SAVE                           ::  ug_vg_have_been_computed = .FALSE.
+    LOGICAL, SAVE                           ::  debugging_output = .TRUE.
     
 !> \mainpage About INIFOR
@@ -199,12 +221,152 @@
  CALL run_control('time', 'alloc')
                       
-
                       CALL average_profile(                                    &
                          input_buffer(output_var % input_id) % array(:,:,:),   &
-                         output_arr(:,:,:), average_imin, average_imax,        &
-                         average_jmin, average_jmax,                           &
-                         output_var % intermediate_grid,                       &
-                         output_var % grid)
- CALL run_control('time', 'comp')
+                         output_arr(0,0,:),                                    &
+                         output_var % averaging_grid)
+
+                      IF ( TRIM(output_var % name) ==                          &
+                           'surface_forcing_surface_pressure' )  THEN
+
+                         IF ( cfg % p0_is_set )  THEN
+                            output_arr(0,0,1) = p0
+                         ELSE
+                            CALL get_surface_pressure(                         &
+                               output_arr(0,0,:), rho_centre,                  &
+                               output_var % averaging_grid)
+                         END IF
+
+                      END IF
+ CALL run_control('time', 'comp')
+
+                   CASE ( 'internal profile' )
+
+                      message = "Averaging of internal profile for variable '" //&
+                         TRIM(output_var % name) // "' is not supported."
+
+                      SELECT CASE (TRIM(output_var % name))
+
+                      CASE('internal_density_centre')
+                         ALLOCATE( rho_centre( 1:output_var % grid % nz ) )
+                         internal_arr => rho_centre
+
+                      CASE('internal_density_north')
+                         ALLOCATE( rho_north( 1:output_var % grid % nz ) )
+                         internal_arr => rho_north
+
+                      CASE('internal_density_south')
+                         ALLOCATE( rho_south( 1:output_var % grid % nz ) )
+                         internal_arr => rho_south
+
+                      CASE('internal_density_east')
+                         ALLOCATE( rho_east( 1:output_var % grid % nz) )
+                         internal_arr => rho_east
+
+                      CASE('internal_density_west')
+                         ALLOCATE( rho_west( 1:output_var % grid % nz ) )
+                         internal_arr => rho_west
+
+                      CASE('internal_pressure_north')
+                         ALLOCATE( p_north( 1:output_var % grid % nz ) )
+                         internal_arr => p_north
+
+                      CASE('internal_pressure_south')
+                         ALLOCATE( p_south( 1:output_var % grid % nz ) )
+                         internal_arr => p_south
+
+                      CASE('internal_pressure_east')
+                         ALLOCATE( p_east( 1:output_var % grid % nz) )
+                         internal_arr => p_east
+
+                      CASE('internal_pressure_west')
+                         ALLOCATE( p_west( 1:output_var % grid % nz ) )
+                         internal_arr => p_west
+
+                      CASE DEFAULT
+                         CALL abort('main loop', message)
+
+                      END SELECT
+ CALL run_control('time', 'alloc')
+
+
+                      CALL average_profile(                                 &
+                         input_buffer(output_var % input_id) % array(:,:,:),&
+                         internal_arr(:),                                   &
+                         output_var % averaging_grid)
+
+                      SELECT CASE (TRIM(output_var % name))
+
+                      CASE('internal_density_centre',                          &
+                           'internal_density_north',                           &
+                           'internal_density_south',                           &
+                           'internal_density_east',                            &
+                           'internal_density_west')
+                         CALL extrapolate_density(internal_arr,                &
+                                                  output_var % averaging_grid)
+
+                      CASE('internal_pressure_north')
+                         CALL extrapolate_pressure(internal_arr, rho_north,    &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_south')
+                         CALL extrapolate_pressure(internal_arr, rho_south,    &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_east')
+                         CALL extrapolate_pressure(internal_arr, rho_east,     &
+                                                   output_var % averaging_grid)
+
+                      CASE('internal_pressure_west')
+                         CALL extrapolate_pressure(internal_arr, rho_west,     &
+                                                   output_var % averaging_grid)
+
+                      CASE DEFAULT
+                         CALL abort('main loop', message)
+
+                      END SELECT
+
+                      IF (.TRUE.)  THEN
+                         ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
+                         output_arr(1,1,:) = internal_arr(:)
+                      END IF
+ CALL run_control('time', 'comp')
+
+
+                   ! This case gets called twice, the first time for ug, the
+                   ! second time for vg. We compute ug and vg at the first call
+                   ! and keep vg (and ug for that matter) around for the second
+                   ! call.
+                   CASE ( 'geostrophic winds' )
+
+                      IF (.NOT. ug_vg_have_been_computed )  THEN
+                         ALLOCATE( ug_arr(output_var % grid % nz) )
+                         ALLOCATE( vg_arr(output_var % grid % nz) )
+
+                         IF ( cfg % ug_is_set )  THEN
+                            ug_arr = cfg % ug
+                            vg_arr = cfg % vg
+                         ELSE
+                            CALL geostrophic_winds( p_north, p_south, p_east,     &
+                                                    p_west, rho_centre, f3,      &
+                                                    averaging_width_ew,           &
+                                                    averaging_width_ns,           &
+                                                    phi_n, lambda_n, phi_centre, lam_centre, &
+                                                    ug_arr, vg_arr )
+                         END IF
+
+                         ug_vg_have_been_computed = .TRUE.
+
+                      END IF
+
+                      ! Prepare output of geostrophic winds
+                      SELECT CASE(TRIM(output_var % name))
+                      CASE ('ls_forcing_ug')
+                         ug_vg_arr => ug_arr
+                      CASE ('ls_forcing_vg')
+                         ug_vg_arr => vg_arr
+                      END SELECT
+
+                      ALLOCATE( output_arr(1,1,output_var % grid % nz) )
+                      output_arr(1,1,:) = ug_vg_arr(:)
 
                    CASE ( 'average scalar' )
@@ -222,9 +384,9 @@
                       SELECT CASE (TRIM(output_var % name))
 
-                      CASE('ls_forcing_ug')
-                          output_arr(1, 1, :) = ug
-
-                      CASE('ls_forcing_vg')
-                          output_arr(1, 1, :) = vg
+                      !CASE('ls_forcing_ug')
+                      !    output_arr(1, 1, :) = ug
+
+                      !CASE('ls_forcing_vg')
+                      !    output_arr(1, 1, :) = vg
 
                       CASE('nudging_tau')
@@ -256,10 +418,12 @@
 !- Section 2.3: Write current time step of current variable
 !------------------------------------------------------------------------------
-                   message = "Writing variable '" // TRIM(output_var%name) // "'."
-                   CALL report('main loop', message)
-                   CALL update_output(output_var, output_arr, iter, output_file)
+                   IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
+                      message = "Writing variable '" // TRIM(output_var%name) // "'."
+                      CALL report('main loop', message)
+                      CALL update_output(output_var, output_arr, iter, output_file)
  CALL run_control('time', 'write')
-
-                   DEALLOCATE(output_arr)
+                   END IF
+
+                   IF (ALLOCATED(output_arr))  DEALLOCATE(output_arr)
  CALL run_control('time', 'alloc')
 
@@ -267,4 +431,18 @@
 
              END DO ! ouput variables
+
+             IF ( group % kind == 'thermodynamics' )  THEN
+                DEALLOCATE( rho_centre )
+                IF ( .NOT. cfg % ug_is_set )  THEN
+                   DEALLOCATE( rho_north )
+                   DEALLOCATE( rho_south )
+                   DEALLOCATE( rho_east )
+                   DEALLOCATE( rho_west )
+                   DEALLOCATE( p_north )
+                   DEALLOCATE( p_south )
+                   DEALLOCATE( p_east )
+                   DEALLOCATE( p_west )
+                END IF
+             END IF
 
              IF ( group % kind == 'running average' .OR. &
@@ -273,5 +451,5 @@
                 ! accumulated COSMO-DE quantities (precipitation).
              ELSE
-                CALL report('main loop', 'Deallocating input buffer')
+                CALL report('main loop', 'Deallocating input buffer', cfg % debug)
                 DEALLOCATE(input_buffer)
              END IF
@@ -281,5 +459,5 @@
 
           IF (ALLOCATED(input_buffer))  THEN
-             CALL report('main loop', 'Deallocating input buffer')
+             CALL report('main loop', 'Deallocating input buffer', cfg % debug)
              DEALLOCATE(input_buffer)
           END IF
@@ -294,5 +472,5 @@
           END IF
 
-          CALL report('main loop', message)
+          CALL report('main loop', message, cfg % debug)
 
        END IF ! IO group % to_be_processed
