Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
May 7, 2020 3:58:16 PM (4 years ago)
Author:
eckhard
Message:

fixed constant-density pressure extrapolation, respect integer working precision

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/UTIL/inifor/src/inifor_transform.f90

    r4481 r4523  
    2626! -----------------
    2727! $Id$
     28! bugfix: pressure extrapolation
     29! respect integer working precision (iwp) specified in inifor_defs.f90
     30! moved fill_water_cells() routine here
     31! remove unused routine, appropriately renamed constand_density_pressure()
     32!
     33!
     34! 4481 2020-03-31 18:55:54Z maronga
    2835! Use PALM's working precision
    2936! Improved coding style
    … …  
    111118    USE inifor_control
    112119    USE inifor_defs,                                                           &
    113         ONLY: BETA, G, P_SL, PI, RD, T_SL, TO_DEGREES, TO_RADIANS, wp
     120        ONLY: BETA, G, P_SL, PI, RD, T_SL, TO_DEGREES, TO_RADIANS, WATER_ID,   &
     121              iwp, wp
    114122    USE inifor_types
    115123    USE inifor_util,                                                           &
    … …  
    126134    REAL(wp), INTENT(OUT)             ::  out_arr(:)
    127135
    128     INTEGER :: k, l, nz
     136    INTEGER(iwp) :: k, l, nz
    129137
    130138    nz = UBOUND(out_arr, 1)
    … …  
    179187    REAL(wp), INTENT(OUT)             ::  out_arr(0:,0:,:)
    180188
    181     INTEGER :: i, j, k, l, nz
     189    INTEGER(iwp) :: i, j, k, l, nz
    182190
    183191    nz = UBOUND(out_arr, 3)
    … …  
    241249    TYPE(nc_var), INTENT(IN), OPTIONAL ::  ncvar
    242250
    243     INTEGER ::  i, j, k, l
     251    INTEGER(iwp) ::  i, j, k, l
    244252
    245253!
    … …  
    247255    IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) )  THEN
    248256        message = "Output array for '" // TRIM(ncvar%name) // "' has ' more levels (" // &
    249            TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
    250            TRIM(str(UBOUND(invar, 3))) // ")."
     257           TRIM(str(UBOUND(outvar, 3, kind=iwp))) // ") than input variable ("//&
     258           TRIM(str(UBOUND(invar, 3, kind=iwp))) // ")."
    251259        CALL inifor_abort('interpolate_2d', message)
    252260    ENDIF
    … …  
    277285!------------------------------------------------------------------------------!
    278286 SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
    279     REAL(wp), INTENT(IN)              ::  in_arr(0:,0:,0:)
    280     REAL(wp), INTENT(OUT)             ::  out_arr(0:)
    281     INTEGER, INTENT(IN), DIMENSION(:) ::  ii, jj
    282 
    283     INTEGER  ::  i, j, k, l
    284     REAL(wp) ::  ni
     287    REAL(wp), INTENT(IN)                   ::  in_arr(0:,0:,0:)
     288    REAL(wp), INTENT(OUT)                  ::  out_arr(0:)
     289    INTEGER(iwp), INTENT(IN), DIMENSION(:) ::  ii, jj
     290
     291    INTEGER(iwp) ::  i, j, k, l
     292    REAL(wp)     ::  ni
    285293
    286294    IF (SIZE(ii) /= SIZE(jj))  THEN
    287295       message = "Length of 'ii' and 'jj' index lists do not match." //     &
    288           NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //   &
    289           NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
     296          NEW_LINE(' ') // "ii has " // str(SIZE(ii, kind=iwp)) // " elements, " //   &
     297          NEW_LINE(' ') // "jj has " // str(SIZE(jj, kind=iwp)) // "."
    290298       CALL inifor_abort('average_2d', message)
    291299    ENDIF
    … …  
    330338    REAL(wp), DIMENSION(:,:,:), INTENT(OUT) ::  palm_array
    331339    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  intermediate_array
    332     INTEGER ::  nx, ny, nlev
     340    INTEGER(iwp) ::  nx, ny, nlev
    333341
    334342    nx = palm_intermediate%nx
    … …  
    365373    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
    366374
    367     INTEGER ::  i_source, j_source, k_profile, k_source, l, m
     375    INTEGER(iwp) ::  i_source, j_source, k_profile, k_source, l, m
    368376
    369377    REAL ::  ni_columns
    … …  
    422430    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
    423431
    424     INTEGER ::  i_source, j_source, l, nz, nlev
     432    INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
    425433
    426434    REAL(wp) ::  ni_columns
    … …  
    468476    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
    469477
    470     INTEGER ::  i_source, j_source, l, nz, nlev
     478    INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
    471479
    472480    REAL(wp)                            ::  ni_columns
    … …  
    511519
    512520
    513 
    514 
    515 !------------------------------------------------------------------------------!
    516 ! Description:
    517 ! ------------
    518 !> Extrapolates density linearly from the level 'k_min' downwards.
    519 !------------------------------------------------------------------------------!
    520  SUBROUTINE extrapolate_density(rho, avg_grid)
    521     REAL(wp), DIMENSION(:), INTENT(INOUT) ::  rho
    522     TYPE(grid_definition), INTENT(IN)     ::  avg_grid
    523 
    524     REAL(wp) ::  drhodz, dz, zk, rhok
    525     INTEGER  ::  k_min
    526 
    527     k_min  = avg_grid%k_min
    528     zk     = avg_grid%z(k_min)
    529     rhok   = rho(k_min)
    530     dz     = avg_grid%z(k_min + 1) - avg_grid%z(k_min)
    531     drhodz = (rho(k_min + 1) - rho(k_min)) / dz
    532 
    533     rho(1:k_min-1) = rhok + drhodz * (avg_grid%z(1:k_min-1) - zk)
    534 
    535  END SUBROUTINE extrapolate_density
    536 
    537 
    538 !------------------------------------------------------------------------------!
    539 ! Description:
    540 ! ------------
    541 !> Driver for extrapolating pressure from PALM level k_min downwards
    542 !------------------------------------------------------------------------------!
    543  SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
     521!------------------------------------------------------------------------------!
     522! Description:
     523! ------------
     524!> Computes average soil values in COSMO-DE water cells from neighbouring
     525!> non-water cells. This is done as a preprocessing step for the COSMO-DE
     526!> soil input arrays, which contain unphysical values for water cells.
     527!>
     528!> This routine computes the average of up to all nine neighbouring cells
     529!> or keeps the original value, if not at least one non-water neightbour
     530!> is available.
     531!>
     532!> By repeatedly applying this step, soil data can be extrapolated into
     533!> 'water' regions occupying multiple cells, one cell per iteration.
     534!>
     535!> Input parameters:
     536!> -----------------
     537!> soiltyp : 2d map of COSMO-DE soil types
     538!> nz : number of layers in the COSMO-DE soil
     539!> niter : number iterations
     540!>
     541!> Output parameters:
     542!> ------------------
     543!> array : the soil array (i.e. water content or temperature)
     544!------------------------------------------------------------------------------!
     545 SUBROUTINE fill_water_cells(soiltyp, array, nz, niter)
     546    INTEGER(iwp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp
     547    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT)  :: array
     548    INTEGER(iwp), INTENT(IN)                   :: nz, niter
     549
     550    REAL(wp), DIMENSION(nz)                    :: column
     551    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  :: old_soiltyp, new_soiltyp
     552    INTEGER(iwp)                               :: l, i, j, nx, ny, n_cells, ii, jj, iter
     553    INTEGER(iwp), DIMENSION(8)                 :: di, dj
     554
     555    nx = SIZE(array, 1)
     556    ny = SIZE(array, 2)
     557    di = (/ -1, -1, -1, 0,  0,  1, 1, 1 /)
     558    dj = (/ -1,  0,  1, -1, 1, -1, 0, 1 /)
     559
     560    ALLOCATE(old_soiltyp(SIZE(soiltyp,1), &
     561                         SIZE(soiltyp,2) ))
     562
     563    ALLOCATE(new_soiltyp(SIZE(soiltyp,1), &
     564                         SIZE(soiltyp,2) ))
     565
     566    old_soiltyp(:,:) = soiltyp(:,:,1)
     567    new_soiltyp(:,:) = soiltyp(:,:,1)
     568
     569    DO  iter = 1, niter
     570
     571       DO  j = 1, ny
     572       DO  i = 1, nx
     573       
     574          IF (old_soiltyp(i,j) == WATER_ID)  THEN
     575
     576             n_cells = 0
     577             column(:) = 0.0_wp
     578             DO  l = 1, SIZE(di)
     579
     580                ii = MIN(nx, MAX(1_iwp, i + di(l)))
     581                jj = MIN(ny, MAX(1_iwp, j + dj(l)))
     582
     583                IF (old_soiltyp(ii,jj) .NE. WATER_ID)  THEN
     584                   n_cells = n_cells + 1
     585                   column(:) = column(:) + array(ii,jj,:)
     586                ENDIF
     587
     588             ENDDO
     589
     590!
     591!--          Overwrite if at least one non-water neighbour cell is available
     592             IF (n_cells > 0)  THEN
     593                array(i,j,:) = column(:) / n_cells
     594                new_soiltyp(i,j) = 0
     595             ENDIF
     596
     597          ENDIF
     598
     599       ENDDO
     600       ENDDO
     601
     602       old_soiltyp(:,:) = new_soiltyp(:,:)
     603
     604    ENDDO
     605
     606    DEALLOCATE(old_soiltyp, new_soiltyp)
     607
     608 END SUBROUTINE fill_water_cells
     609
     610
     611!------------------------------------------------------------------------------!
     612! Description:
     613! ------------
     614!> Takes the averaged pressure profile <p> and sets the lowest entry to the
     615!> extrapolated pressure at the surface.
     616!------------------------------------------------------------------------------!
     617 SUBROUTINE get_surface_pressure(p, rho, avg_grid)
    544618    REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho
    545619    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p
    546620    TYPE(grid_definition), INTENT(IN)     ::  avg_grid
    547621
    548     REAL(wp) ::  drhodz, dz, zk, rhok
    549     INTEGER  ::  k, k_min
     622    REAL(wp)     ::  drhodz, dz, zk, rhok
     623    INTEGER(iwp) :: k_min
    550624
    551625    k_min = avg_grid%k_min
    … …  
    553627    rhok  = rho(k_min)
    554628    dz    = avg_grid%z(k_min + 1) - avg_grid%z(k_min)
    555     drhodz = 0.5_wp * (rho(k_min + 1) - rho(k_min)) / dz
    556 
    557     DO  k = 1, k_min-1
    558        p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz,         &
    559                                         avg_grid%z(k), G)
    560     ENDDO
    561 
    562  END SUBROUTINE extrapolate_pressure
    563 
    564 
    565 !------------------------------------------------------------------------------!
    566 ! Description:
    567 ! ------------
    568 !> Takes the averaged pressure profile <p> and sets the lowest entry to the
    569 !> extrapolated pressure at the surface.
    570 !------------------------------------------------------------------------------!
    571  SUBROUTINE get_surface_pressure(p, rho, avg_grid)
    572     REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho
    573     REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p
    574     TYPE(grid_definition), INTENT(IN)     ::  avg_grid
    575 
    576     REAL(wp) ::  drhodz, dz, zk, rhok
    577     INTEGER  ::  k_min
    578 
    579     k_min = avg_grid%k_min
    580     zk    = avg_grid%z(k_min)
    581     rhok  = rho(k_min)
    582     dz    = avg_grid%z(k_min + 1) - avg_grid%z(k_min)
    583     drhodz = 0.5_wp * (rho(k_min + 1) - rho(k_min)) / dz
    584 
    585     p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz,            &
    586                                      0.0_wp, G)
     629    drhodz = ( rho(k_min + 1) - rho(k_min) ) / dz
     630
     631    p(1) = linear_density_pressure( p(k_min), zk, rhok, drhodz,                &
     632                                    z = 0.0_wp, g=G )
    587633
    588634 END SUBROUTINE get_surface_pressure
    589635
    590636
    591  FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
     637 FUNCTION linear_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
    592638
    593639    REAL(wp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
    594640    REAL(wp) ::  p
    595641
    596     p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) )
    597 
    598  END FUNCTION constant_density_pressure
     642    p = pk + ( zk - z ) * g * ( rhok - 0.5_wp * drhodz * (zk - z) )
     643
     644 END FUNCTION linear_density_pressure
    599645
    600646!-----------------------------------------------------------------------------!
    … …  
    769815    REAL(wp), INTENT(OUT) ::  phi(0:,0:), lam(0:,0:)
    770816
    771     INTEGER ::  i, j
     817    INTEGER(iwp) ::  i, j
    772818   
    773819    IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
    … …  
    816862    REAL(wp), INTENT(IN)                  :: angle !< rotation angle [deg]
    817863
    818     INTEGER  :: i
    819     REAL(wp) :: sine, cosine, v_rot(2), rotation(2,2)
     864    INTEGER(iwp) :: i
     865    REAL(wp)     :: sine, cosine, v_rot(2), rotation(2,2)
    820866
    821867    sine = SIN(angle * TO_RADIANS)
    … …  
    912958                                       palm_ii, palm_jj)
    913959
    914     REAL(wp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
    915     REAL(wp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
    916     REAL(wp)                                   ::  cosmo_dxi, cosmo_dyi
    917     INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_ii, palm_jj
    918 
    919     REAL(wp) ::  lonpos, latpos, lon0, lat0
    920     INTEGER  ::  i, j
     960    REAL(wp), DIMENSION(0:), INTENT(IN)            ::  cosmo_lat, cosmo_lon
     961    REAL(wp), DIMENSION(0:,0:), INTENT(IN)         ::  palm_clat, palm_clon
     962    REAL(wp)                                       ::  cosmo_dxi, cosmo_dyi
     963    INTEGER(iwp), DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_ii, palm_jj
     964
     965    REAL(wp)     ::  lonpos, latpos, lon0, lat0
     966    INTEGER(iwp) ::  i, j
    921967
    922968    lon0 = cosmo_lon(0)
    … …  
    9681014!> column of the given palm grid.
    9691015!------------------------------------------------------------------------------!
    970  SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
    971                                                             palm_intermediate )
     1016 SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,            &
     1017                                                         palm_intermediate )
    9721018    TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
    9731019    TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
    9741020
    975     INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
    976     LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
    977                  point_is_in_current_cell
    978     REAL(wp) ::  current_height, column_base, column_top, h_top, h_bottom, &
    979                  weight
     1021    INTEGER(iwp) ::  i, j, k, nx, ny, nz, nlev, k_intermediate
     1022    LOGICAL      ::  point_is_below_grid, point_is_above_grid,                 &
     1023                     point_is_in_current_cell
     1024    REAL(wp)     ::  current_height, column_base, column_top, h_top, h_bottom, &
     1025                     weight
    9801026
    9811027    nx   = palm_grid%nx
    … …  
    10911137    LOGICAL                                      ::  level_based_averaging
    10921138
    1093     INTEGER           ::  i, j, k_palm, k_intermediate, l, nlev
     1139    INTEGER(iwp)      ::  i, j, k_palm, k_intermediate, l, nlev
    10941140    LOGICAL           ::  point_is_below_grid, point_is_above_grid,         &
    10951141                          point_is_in_current_cell
    … …  
    12591305    palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
    12601306       
    1261     REAL(wp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
    1262     REAL(wp)                                   ::  cosmo_dxi, cosmo_dyi
    1263     REAL(wp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
    1264     INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)  ::  palm_ii, palm_jj
     1307    REAL(wp), DIMENSION(0:), INTENT(IN)           ::  cosmo_lat, cosmo_lon
     1308    REAL(wp)                                      ::  cosmo_dxi, cosmo_dyi
     1309    REAL(wp), DIMENSION(0:,0:), INTENT(IN)        ::  palm_clat, palm_clon
     1310    INTEGER(iwp), DIMENSION(0:,0:,1:), INTENT(IN) ::  palm_ii, palm_jj
    12651311
    12661312    REAL(wp), DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_w_horiz
    12671313
    1268     REAL(wp) ::  wlambda, wphi
    1269     INTEGER  ::  i, j
     1314    REAL(wp)     ::  wlambda, wphi
     1315    INTEGER(iwp) ::  i, j
    12701316
    12711317    cosmo_dxi = 1.0_wp / (cosmo_lon(1) - cosmo_lon(0))
    … …  
    13191365    REAL(wp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
    13201366    REAL(wp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  u_centre, v_centre
    1321     INTEGER ::  nx, ny
     1367    INTEGER(iwp) ::  nx, ny
    13221368
    13231369    nx = UBOUND(u_face, 1)
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |