Changeset 1960 for palm/trunk/SOURCE/check_parameters.f90
- Timestamp:
- Jul 12, 2016 4:34:24 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/check_parameters.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r1932 r1960 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Separate checks and calculations for humidity and passive scalar. Therefore, 22 ! a subroutine to calculate vertical gradients is introduced, in order to reduce 23 ! redundance. 24 ! Additional check - large-scale forcing is not implemented for passive scalars 25 ! so far. 22 26 ! 23 27 ! Former revisions: … … 1053 1057 ENDIF 1054 1058 1055 IF ( passive_scalar .AND. humidity ) THEN1056 message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &1057 'is not allowed simultaneously'1058 CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )1059 ENDIF1060 1061 1059 ! 1062 1060 !-- When plant canopy model is used, peform addtional checks … … 1094 1092 !-- Initial profiles for 1D and 3D model, respectively (u,v further below) 1095 1093 pt_init = pt_surface 1096 IF ( humidity ) THEN 1097 q_init = q_surface 1098 ENDIF 1099 IF ( ocean ) sa_init = sa_surface 1100 IF ( passive_scalar ) q_init = s_surface 1094 IF ( humidity ) q_init = q_surface 1095 IF ( ocean ) sa_init = sa_surface 1096 IF ( passive_scalar ) s_init = s_surface 1101 1097 1102 1098 ! … … 1282 1278 !-- Compute initial temperature profile using the given temperature gradients 1283 1279 IF ( .NOT. neutral ) THEN 1284 1285 i = 1 1286 gradient = 0.0_wp 1287 1288 IF ( .NOT. ocean ) THEN 1289 1290 pt_vertical_gradient_level_ind(1) = 0 1291 DO k = 1, nzt+1 1292 IF ( i < 11 ) THEN 1293 IF ( pt_vertical_gradient_level(i) < zu(k) .AND. & 1294 pt_vertical_gradient_level(i) >= 0.0_wp ) THEN 1295 gradient = pt_vertical_gradient(i) / 100.0_wp 1296 pt_vertical_gradient_level_ind(i) = k - 1 1297 i = i + 1 1298 ENDIF 1299 ENDIF 1300 IF ( gradient /= 0.0_wp ) THEN 1301 IF ( k /= 1 ) THEN 1302 pt_init(k) = pt_init(k-1) + dzu(k) * gradient 1303 ELSE 1304 pt_init(k) = pt_surface + dzu(k) * gradient 1305 ENDIF 1306 ELSE 1307 pt_init(k) = pt_init(k-1) 1308 ENDIF 1309 ENDDO 1310 1311 ELSE 1312 1313 pt_vertical_gradient_level_ind(1) = nzt+1 1314 DO k = nzt, 0, -1 1315 IF ( i < 11 ) THEN 1316 IF ( pt_vertical_gradient_level(i) > zu(k) .AND. & 1317 pt_vertical_gradient_level(i) <= 0.0_wp ) THEN 1318 gradient = pt_vertical_gradient(i) / 100.0_wp 1319 pt_vertical_gradient_level_ind(i) = k + 1 1320 i = i + 1 1321 ENDIF 1322 ENDIF 1323 IF ( gradient /= 0.0_wp ) THEN 1324 IF ( k /= nzt ) THEN 1325 pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient 1326 ELSE 1327 pt_init(k) = pt_surface - 0.5_wp * dzu(k+1) * gradient 1328 pt_init(k+1) = pt_surface + 0.5_wp * dzu(k+1) * gradient 1329 ENDIF 1330 ELSE 1331 pt_init(k) = pt_init(k+1) 1332 ENDIF 1333 ENDDO 1334 1335 ENDIF 1336 1337 ENDIF 1338 1339 ! 1340 !-- In case of no given temperature gradients, choose gradient of neutral 1341 !-- stratification 1342 IF ( pt_vertical_gradient_level(1) == -9999999.9_wp ) THEN 1343 pt_vertical_gradient_level(1) = 0.0_wp 1344 ENDIF 1345 1346 ! 1347 !-- Store temperature gradient at the top boundary for possible Neumann 1348 !-- boundary condition 1349 bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1) 1350 1351 ! 1352 !-- If required, compute initial humidity or scalar profile using the given 1353 !-- humidity/scalar gradient. In case of scalar transport, initially store 1354 !-- values of the scalar parameters on humidity parameters 1280 CALL init_vertical_profiles( pt_vertical_gradient_level_ind, & 1281 pt_vertical_gradient_level, & 1282 pt_vertical_gradient, pt_init, & 1283 pt_surface, bc_pt_t_val ) 1284 ENDIF 1285 ! 1286 !-- Compute initial humidity profile using the given humidity gradients 1287 IF ( humidity ) THEN 1288 CALL init_vertical_profiles( q_vertical_gradient_level_ind, & 1289 q_vertical_gradient_level, & 1290 q_vertical_gradient, q_init, & 1291 q_surface, bc_q_t_val ) 1292 ENDIF 1293 ! 1294 !-- Compute initial scalar profile using the given scalar gradients 1355 1295 IF ( passive_scalar ) THEN 1356 bc_q_b = bc_s_b 1357 bc_q_t = bc_s_t 1358 q_surface = s_surface 1359 q_surface_initial_change = s_surface_initial_change 1360 q_vertical_gradient = s_vertical_gradient 1361 q_vertical_gradient_level = s_vertical_gradient_level 1362 surface_waterflux = surface_scalarflux 1363 wall_humidityflux = wall_scalarflux 1364 ENDIF 1365 1366 IF ( humidity .OR. passive_scalar ) THEN 1367 1368 i = 1 1369 gradient = 0.0_wp 1370 1371 IF ( .NOT. ocean ) THEN 1372 1373 q_vertical_gradient_level_ind(1) = 0 1374 DO k = 1, nzt+1 1375 IF ( i < 11 ) THEN 1376 IF ( q_vertical_gradient_level(i) < zu(k) .AND. & 1377 q_vertical_gradient_level(i) >= 0.0_wp ) THEN 1378 gradient = q_vertical_gradient(i) / 100.0_wp 1379 q_vertical_gradient_level_ind(i) = k - 1 1380 i = i + 1 1381 ENDIF 1382 ENDIF 1383 IF ( gradient /= 0.0_wp ) THEN 1384 IF ( k /= 1 ) THEN 1385 q_init(k) = q_init(k-1) + dzu(k) * gradient 1386 ELSE 1387 q_init(k) = q_init(k-1) + dzu(k) * gradient 1388 ENDIF 1389 ELSE 1390 q_init(k) = q_init(k-1) 1391 ENDIF 1392 ! 1393 !-- Avoid negative humidities 1394 IF ( q_init(k) < 0.0_wp ) THEN 1395 q_init(k) = 0.0_wp 1396 ENDIF 1397 ENDDO 1398 1399 ELSE 1400 1401 q_vertical_gradient_level_ind(1) = nzt+1 1402 DO k = nzt, 0, -1 1403 IF ( i < 11 ) THEN 1404 IF ( q_vertical_gradient_level(i) > zu(k) .AND. & 1405 q_vertical_gradient_level(i) <= 0.0_wp ) THEN 1406 gradient = q_vertical_gradient(i) / 100.0_wp 1407 q_vertical_gradient_level_ind(i) = k + 1 1408 i = i + 1 1409 ENDIF 1410 ENDIF 1411 IF ( gradient /= 0.0_wp ) THEN 1412 IF ( k /= nzt ) THEN 1413 q_init(k) = q_init(k+1) - dzu(k+1) * gradient 1414 ELSE 1415 q_init(k) = q_surface - 0.5_wp * dzu(k+1) * gradient 1416 q_init(k+1) = q_surface + 0.5_wp * dzu(k+1) * gradient 1417 ENDIF 1418 ELSE 1419 q_init(k) = q_init(k+1) 1420 ENDIF 1421 ! 1422 !-- Avoid negative humidities 1423 IF ( q_init(k) < 0.0_wp ) THEN 1424 q_init(k) = 0.0_wp 1425 ENDIF 1426 ENDDO 1427 1428 ENDIF 1429 ! 1430 !-- In case of no given humidity gradients, choose zero gradient 1431 !-- conditions 1432 IF ( q_vertical_gradient_level(1) == -1.0_wp ) THEN 1433 q_vertical_gradient_level(1) = 0.0_wp 1434 ENDIF 1435 ! 1436 !-- Store humidity, rain water content and rain drop concentration 1437 !-- gradient at the top boundary for possile Neumann boundary condition 1438 bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1) 1439 ENDIF 1440 1296 CALL init_vertical_profiles( s_vertical_gradient_level_ind, & 1297 s_vertical_gradient_level, & 1298 s_vertical_gradient, s_init, & 1299 s_surface, bc_s_t_val ) 1300 ENDIF 1441 1301 ! 1442 1302 !-- If required, compute initial salinity profile using the given salinity 1443 1303 !-- gradients 1444 1304 IF ( ocean ) THEN 1445 1446 i = 1 1447 gradient = 0.0_wp 1448 1449 sa_vertical_gradient_level_ind(1) = nzt+1 1450 DO k = nzt, 0, -1 1451 IF ( i < 11 ) THEN 1452 IF ( sa_vertical_gradient_level(i) > zu(k) .AND. & 1453 sa_vertical_gradient_level(i) <= 0.0_wp ) THEN 1454 gradient = sa_vertical_gradient(i) / 100.0_wp 1455 sa_vertical_gradient_level_ind(i) = k + 1 1456 i = i + 1 1457 ENDIF 1458 ENDIF 1459 IF ( gradient /= 0.0_wp ) THEN 1460 IF ( k /= nzt ) THEN 1461 sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient 1462 ELSE 1463 sa_init(k) = sa_surface - 0.5_wp * dzu(k+1) * gradient 1464 sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient 1465 ENDIF 1466 ELSE 1467 sa_init(k) = sa_init(k+1) 1468 ENDIF 1469 ENDDO 1470 1305 CALL init_vertical_profiles( sa_vertical_gradient_level_ind, & 1306 sa_vertical_gradient_level, & 1307 sa_vertical_gradient, s_init, & 1308 sa_surface, -999.0_wp ) 1471 1309 ENDIF 1472 1310 … … 1860 1698 1861 1699 ! 1862 !-- In case of humidity or passive scalar, set boundary conditions for total 1863 !-- water content / scalar 1864 IF ( humidity .OR. passive_scalar ) THEN 1865 IF ( humidity ) THEN 1866 sq = 'q' 1867 ELSE 1868 sq = 's' 1869 ENDIF 1870 IF ( bc_q_b == 'dirichlet' ) THEN 1871 ibc_q_b = 0 1872 ELSEIF ( bc_q_b == 'neumann' ) THEN 1873 ibc_q_b = 1 1874 ELSE 1875 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 1876 '_b ="' // TRIM( bc_q_b ) // '"' 1877 CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 ) 1878 ENDIF 1879 IF ( bc_q_t == 'dirichlet' ) THEN 1880 ibc_q_t = 0 1881 ELSEIF ( bc_q_t == 'neumann' ) THEN 1882 ibc_q_t = 1 1883 ELSEIF ( bc_q_t == 'nested' ) THEN 1884 ibc_q_t = 3 1885 ELSE 1886 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 1887 '_t ="' // TRIM( bc_q_t ) // '"' 1888 CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 ) 1889 ENDIF 1700 !-- Set boundary conditions for total water content 1701 IF ( humidity ) THEN 1702 CALL check_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t, & 1703 'PA0071', 'PA0072', 'PA0073', 'PA0074', & 1704 surface_waterflux, constant_waterflux, & 1705 q_surface_initial_change ) 1890 1706 1891 1707 IF ( surface_waterflux == 9999999.9_wp ) THEN … … 1909 1725 ENDIF 1910 1726 1911 ! 1912 !-- A given surface humidity implies Dirichlet boundary condition for 1913 !-- humidity. In this case specification of a constant water flux is 1914 !-- forbidden. 1915 IF ( ibc_q_b == 0 .AND. constant_waterflux ) THEN 1916 message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // & 1917 '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // & 1918 'th prescribed surface flux' 1919 CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 ) 1920 ENDIF 1921 IF ( constant_waterflux .AND. q_surface_initial_change /= 0.0_wp ) THEN 1922 WRITE( message_string, * ) 'a prescribed surface flux is not allo', & 1923 'wed with ', sq, '_surface_initial_change (/=0) = ', & 1924 q_surface_initial_change 1925 CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 ) 1926 ENDIF 1927 1727 ENDIF 1728 1729 IF ( passive_scalar ) THEN 1730 CALL check_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t, & 1731 'PA0071', 'PA0072', 'PA0073', 'PA0074', & 1732 surface_scalarflux, constant_scalarflux, & 1733 s_surface_initial_change ) 1928 1734 ENDIF 1929 1735 ! … … 3850 3656 'the usage of large scale forcing from external file.' 3851 3657 CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 ) 3852 3658 ENDIF 3853 3659 3854 3660 IF ( large_scale_forcing .AND. ( .NOT. humidity ) ) THEN … … 3856 3662 'file LSF_DATA requires humidity = .T..' 3857 3663 CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 ) 3858 ENDIF 3664 ENDIF 3665 3666 IF ( large_scale_forcing .AND. passive_scalar ) THEN 3667 message_string = 'The usage of large scale forcing from external &'// & 3668 'file LSF_DATA is not implemented for psassive scalars' 3669 CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 ) 3670 ENDIF 3859 3671 3860 3672 IF ( large_scale_forcing .AND. topography /= 'flat' ) THEN … … 3864 3676 ENDIF 3865 3677 3866 IF ( large_scale_forcing .AND. ocean 3678 IF ( large_scale_forcing .AND. ocean ) THEN 3867 3679 message_string = 'The usage of large scale forcing from external &'// & 3868 3680 'file LSF_DATA is not implemented for ocean runs' … … 3943 3755 3944 3756 END SUBROUTINE check_dt_do 3757 3758 3759 3760 3761 !------------------------------------------------------------------------------! 3762 ! Description: 3763 ! ------------ 3764 !> Inititalizes the vertical profiles of scalar quantities. 3765 !------------------------------------------------------------------------------! 3766 3767 SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind, & 3768 vertical_gradient_level, & 3769 vertical_gradient, & 3770 pr_init, surface_value, bc_t_val ) 3771 3772 3773 IMPLICIT NONE 3774 3775 INTEGER(iwp) :: i !< counter 3776 INTEGER(iwp), DIMENSION(1:10) :: vertical_gradient_level_ind !< vertical grid indices for gradient levels 3777 3778 REAL(wp) :: bc_t_val !< model top gradient 3779 REAL(wp) :: gradient !< vertica gradient of the respective quantity 3780 REAL(wp) :: surface_value !< surface value of the respecitve quantity 3781 3782 3783 REAL(wp), DIMENSION(0:nz+1) :: pr_init !< initialisation profile 3784 REAL(wp), DIMENSION(1:10) :: vertical_gradient !< given vertical gradient 3785 REAL(wp), DIMENSION(1:10) :: vertical_gradient_level !< given vertical gradient level 3786 3787 i = 1 3788 gradient = 0.0_wp 3789 3790 IF ( .NOT. ocean ) THEN 3791 3792 vertical_gradient_level_ind(1) = 0 3793 DO k = 1, nzt+1 3794 IF ( i < 11 ) THEN 3795 IF ( vertical_gradient_level(i) < zu(k) .AND. & 3796 vertical_gradient_level(i) >= 0.0_wp ) THEN 3797 gradient = vertical_gradient(i) / 100.0_wp 3798 vertical_gradient_level_ind(i) = k - 1 3799 i = i + 1 3800 ENDIF 3801 ENDIF 3802 IF ( gradient /= 0.0_wp ) THEN 3803 IF ( k /= 1 ) THEN 3804 pr_init(k) = pr_init(k-1) + dzu(k) * gradient 3805 ELSE 3806 pr_init(k) = pr_init(k-1) + dzu(k) * gradient 3807 ENDIF 3808 ELSE 3809 pr_init(k) = pr_init(k-1) 3810 ENDIF 3811 ! 3812 !-- Avoid negative values 3813 IF ( pr_init(k) < 0.0_wp ) THEN 3814 pr_init(k) = 0.0_wp 3815 ENDIF 3816 ENDDO 3817 3818 ELSE 3819 3820 vertical_gradient_level_ind(1) = nzt+1 3821 DO k = nzt, 0, -1 3822 IF ( i < 11 ) THEN 3823 IF ( vertical_gradient_level(i) > zu(k) .AND. & 3824 vertical_gradient_level(i) <= 0.0_wp ) THEN 3825 gradient = vertical_gradient(i) / 100.0_wp 3826 vertical_gradient_level_ind(i) = k + 1 3827 i = i + 1 3828 ENDIF 3829 ENDIF 3830 IF ( gradient /= 0.0_wp ) THEN 3831 IF ( k /= nzt ) THEN 3832 pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient 3833 ELSE 3834 pr_init(k) = surface_value - 0.5_wp * dzu(k+1) * gradient 3835 pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient 3836 ENDIF 3837 ELSE 3838 pr_init(k) = pr_init(k+1) 3839 ENDIF 3840 ! 3841 !-- Avoid negative humidities 3842 IF ( pr_init(k) < 0.0_wp ) THEN 3843 pr_init(k) = 0.0_wp 3844 ENDIF 3845 ENDDO 3846 3847 ENDIF 3848 3849 ! 3850 !-- In case of no given gradients, choose zero gradient conditions 3851 IF ( vertical_gradient_level(1) == -1.0_wp ) THEN 3852 vertical_gradient_level(1) = 0.0_wp 3853 ENDIF 3854 ! 3855 !-- Store gradient at the top boundary for possible Neumann boundary condition 3856 bc_t_val = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1) 3857 3858 END SUBROUTINE init_vertical_profiles 3859 3860 3861 3862 !------------------------------------------------------------------------------! 3863 ! Description: 3864 ! ------------ 3865 !> Check the bottom and top boundary conditions for humidity and scalars. 3866 !------------------------------------------------------------------------------! 3867 3868 SUBROUTINE check_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, & 3869 err_nr_b, err_nr_t, err_nr_3, err_nr_4, & 3870 surface_flux, constant_flux, surface_initial_change ) 3871 3872 3873 IMPLICIT NONE 3874 3875 CHARACTER (LEN=1) :: sq !< 3876 CHARACTER (LEN=*) :: bc_b 3877 CHARACTER (LEN=*) :: bc_t 3878 CHARACTER (LEN=*) :: err_nr_b 3879 CHARACTER (LEN=*) :: err_nr_t 3880 CHARACTER (LEN=*) :: err_nr_3 3881 CHARACTER (LEN=*) :: err_nr_4 3882 3883 INTEGER(iwp) :: ibc_b 3884 INTEGER(iwp) :: ibc_t 3885 3886 LOGICAL :: constant_flux 3887 3888 REAL(wp) :: surface_flux 3889 REAL(wp) :: surface_initial_change 3890 3891 3892 IF ( bc_b == 'dirichlet' ) THEN 3893 ibc_b = 0 3894 ELSEIF ( bc_b == 'neumann' ) THEN 3895 ibc_b = 1 3896 ELSE 3897 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 3898 '_b ="' // TRIM( bc_b ) // '"' 3899 CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 ) 3900 ENDIF 3901 3902 IF ( bc_t == 'dirichlet' ) THEN 3903 ibc_t = 0 3904 ELSEIF ( bc_t == 'neumann' ) THEN 3905 ibc_t = 1 3906 ELSEIF ( bc_t == 'nested' ) THEN 3907 ibc_t = 3 3908 ELSE 3909 message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // & 3910 '_t ="' // TRIM( bc_t ) // '"' 3911 CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 ) 3912 ENDIF 3913 3914 ! 3915 !-- A given surface value implies Dirichlet boundary condition for 3916 !-- the respective quantity. In this case specification of a constant flux is 3917 !-- forbidden. 3918 IF ( ibc_b == 0 .AND. constant_flux ) THEN 3919 message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // & 3920 '= "' // TRIM( bc_b ) // '" is not allowed wi' // & 3921 'th prescribed surface flux' 3922 CALL message( 'check_parameters', err_nr_3, 1, 2, 0, 6, 0 ) 3923 ENDIF 3924 IF ( constant_waterflux .AND. surface_initial_change /= 0.0_wp ) THEN 3925 WRITE( message_string, * ) 'a prescribed surface flux is not allo', & 3926 'wed with ', sq, '_surface_initial_change (/=0) = ', & 3927 surface_initial_change 3928 CALL message( 'check_parameters', err_nr_4, 1, 2, 0, 6, 0 ) 3929 ENDIF 3930 3931 3932 END SUBROUTINE check_bc_scalars 3933 3934 3945 3935 3946 3936 END SUBROUTINE check_parameters
Note: See TracChangeset
for help on using the changeset viewer.
