Changeset 4553 for palm/trunk/UTIL/inifor/src/inifor_grid.f90

Timestamp:

Jun 3, 2020 4:34:15 PM (4 years ago)

Author:

eckhard

Message:

Fixed domain extend check, readablity and documentation improvements

File:

• palm/trunk/UTIL/inifor/src/inifor_grid.f90 (modified) (11 diffs)

palm/trunk/UTIL/inifor/src/inifor_grid.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Fixed domain extent check
+! Added and fixed code documentation
+! 
+! 
+! 4538 2020-05-18 13:45:35Z eckhard
 ! Modularize setup of PALM origin
 ! 
@@ -687 +692 @@
     averaging_width_ns = averaging_angle * EARTH_RADIUS
 
-    lonmin_tot = MIN(lam_centre - averaging_angle, lonmin_palm)
-    lonmax_tot = MAX(lam_centre + averaging_angle, lonmax_palm)
-    latmin_tot = MIN(phi_centre - averaging_angle, latmin_palm)
-    latmax_tot = MAX(phi_centre + averaging_angle, latmax_palm)
-
-    palm_domain_outside_cosmo = ANY(                                           &
-       (/ lonmin_tot,   -lonmax_tot,   latmin_tot,   -latmax_tot/) .LT.        &
-       (/ lonmin_cosmo, -lonmax_cosmo, latmin_cosmo, -latmax_cosmo/)           &
-    )
-
-    IF ( palm_domain_outside_cosmo )  THEN
-       message = 'PALM domain or geostrophic averaging domains extend ' //     &
-                 'outside COSMO domain.'
-       CALL inifor_abort( 'setup_parameters', message )
-    ENDIF
-
-
 !
 !-- Coriolis parameter
@@ -1077 +1065 @@
     latmax_palm = MAXVAL(palm_intermediate%clat)
 
+    lonmin_tot = MIN(lam_centre - averaging_angle, lonmin_palm)
+    lonmax_tot = MAX(lam_centre + averaging_angle, lonmax_palm)
+    latmin_tot = MIN(phi_centre - averaging_angle, latmin_palm)
+    latmax_tot = MAX(phi_centre + averaging_angle, latmax_palm)
+
+    palm_domain_outside_cosmo = ANY(                                           &
+       (/ lonmin_tot,   -lonmax_tot,   latmin_tot,   -latmax_tot/) .LT.        &
+       (/ lonmin_cosmo, -lonmax_cosmo, latmin_cosmo, -latmax_cosmo/)           &
+    )
+
+    IF ( palm_domain_outside_cosmo )  THEN
+       message = 'PALM domain or geostrophic averaging domains extend ' //     &
+                 'outside COSMO domain.'
+       CALL inifor_abort( 'setup_grids', message )
+    ENDIF
+
+
     CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid,   &
             x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
@@ -1268 +1273 @@
 
     IF (setup_volumetric)  THEN
-       ALLOCATE( intermediate_grid%h(0:intermediate_grid%nx,            &
-                                     0:intermediate_grid%ny,            &
-                                     0:intermediate_grid%nz) ) 
-       intermediate_grid%h(:,:,:) = - EARTH_RADIUS
+       ALLOCATE( intermediate_grid%intermediate_h(0:intermediate_grid%nx,      &
+                                                  0:intermediate_grid%ny,      &
+                                                  0:intermediate_grid%nz) ) 
+       intermediate_grid%intermediate_h(:,:,:) = - EARTH_RADIUS
 
 !
 !--    For w points, use hhl, for scalars use hfl
 !--    compute the full heights for the intermediate grids
-       CALL interpolate_2d(cosmo_h, intermediate_grid%h, intermediate_grid)
+       CALL interpolate_2d(cosmo_h, intermediate_grid%intermediate_h, intermediate_grid)
        CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
     ENDIF
@@ -1550 +1555 @@
 !> kind : kind of quantity to be averaged using this averaging grid.
 !>    Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
-!>    parent/COSMO columns for the region in get_cosmo_averaging_region() is
+!>    parent/COSMO columns for the region in get_latlon_averaging_region() is
 !>    independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
 !>    components, which are computed in the preprocessing step.
@@ -1595 +1600 @@
 !-- Find and store COSMO columns that fall into the coordinate range
 !-- given by avg_grid%clon, %clat
-    CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
+    CALL get_latlon_averaging_region(avg_grid, cosmo_grid)
 
     ALLOCATE (avg_grid%kkk(avg_grid%n_columns, avg_grid%nz, 2) )
@@ -1612 +1617 @@
           message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
                     "' is not supported. Use 'scalar', 'u', or 'v'."
-          CALL inifor_abort('get_cosmo_averaging_region', message)
+          CALL inifor_abort('get_latlon_averaging_region', message)
 
     END SELECT
 
 !
-!-- For level-besed averaging, compute average heights
+!-- For level-based averaging, compute average heights
     level_based_averaging = ( TRIM(cfg%averaging_mode) == 'level' )
     IF (level_based_averaging)  THEN
-       ALLOCATE(avg_grid%h(1,1,SIZE(avg_grid%cosmo_h, 3)) )
+       ALLOCATE(avg_grid%intermediate_h(1,1,SIZE(avg_grid%cosmo_h, 3)) )
  
-       CALL average_2d(avg_grid%cosmo_h, avg_grid%h(1,1,:),             &
+       CALL average_2d(avg_grid%cosmo_h, avg_grid%intermediate_h(1,1,:),       &
                        avg_grid%iii, avg_grid%jjj)
 
@@ -1636 +1641 @@
 
 
- SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_latlon_averaging_region() finds all mesocsale columns within the
+!> latitude-longitude reactagle given by the four values in avg_grid%lon(1:2)
+!> and %lat(1:2). The total number of all found columns is stored in
+!> avg_grid%n_columns, and their indices are stored in the sequential lists
+!> avg_grid%iii(:) and %jjj(:).
+!------------------------------------------------------------------------------!
+ SUBROUTINE get_latlon_averaging_region(avg_grid, cosmo_grid)
     TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
     TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
@@ -1663 +1677 @@
           message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
                     "' is not supported. Use 'scalar', 'u', or 'v'."
-          CALL inifor_abort('get_cosmo_averaging_region', message)
+          CALL inifor_abort('get_latlon_averaging_region', message)
 
     END SELECT
 
-!
-!-- FIXME: make dlon, dlat parameters of the grid_defintion type
     dlon = cosmo_lon(1) - cosmo_lon(0)
     dlat = cosmo_lat(1) - cosmo_lat(0)
@@ -1682 +1694 @@
               " and " //                                                    &
               TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax))
-    CALL report( 'get_cosmo_averaging_region', message )
+    CALL report( 'get_latlon_averaging_region', message )
 
     nx = imax - imin + 1
@@ -1700 +1712 @@
     ENDDO
 
- END SUBROUTINE get_cosmo_averaging_region
+ END SUBROUTINE get_latlon_averaging_region