Changeset 1709 for palm/trunk
- Timestamp:
- Nov 4, 2015 2:47:01 PM (9 years ago)
- Location:
- palm/trunk
- Files:
-
- 8 edited
-
DOC/palm2doxygen.config (modified) (1 diff)
-
SOURCE/flow_statistics.f90 (modified) (3 diffs)
-
SOURCE/init_1d_model.f90 (modified) (4 diffs)
-
SOURCE/land_surface_model.f90 (modified) (25 diffs)
-
SOURCE/radiation_model.f90 (modified) (12 diffs)
-
SOURCE/read_3d_binary.f90 (modified) (4 diffs)
-
SOURCE/surface_layer_fluxes.f90 (modified) (9 diffs)
-
SOURCE/write_3d_binary.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/DOC/palm2doxygen.config
r1684 r1709 812 812 # exclude all test directories for example use the pattern */test/* 813 813 814 EXCLUDE_PATTERNS = */message.f90 */cpulog.f90 814 EXCLUDE_PATTERNS = */message.f90 */cpulog.f90 */tridia_solver.f90 815 815 816 816 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names -
palm/trunk/SOURCE/flow_statistics.f90
r1701 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Updated output of Obukhov length 22 22 ! 23 23 ! Former revisions: … … 1554 1554 ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0 1555 1555 ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0 1556 1557 ! 1558 !-- Calculate obukhov length 1559 IF ( ts_value(5,sr) /= 0.0_wp ) THEN 1560 ! IF ( most_method == 'circular' ) THEN 1561 ! ts_value(22,sr) = ts_value(4,sr)**2 / & 1562 ! ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) 1563 ! ELSE 1564 ts_value(22,sr) = hom(nzb,1,114,sr) 1565 ! ENDIF 1556 1557 IF ( .NOT. neutral ) THEN 1558 ts_value(22,sr) = hom(nzb,1,114,sr) ! L 1566 1559 ELSE 1567 ts_value(22,sr) = 1 0000.0_wp1560 ts_value(22,sr) = 1.0E10_wp 1568 1561 ENDIF 1569 1562 … … 3503 3496 ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0 3504 3497 3505 IF ( ts_value(5,sr) /= 0.0_wp) THEN3506 ts_value(22,sr) = hom(nzb,1,114,sr) 3498 IF ( .NOT. neutral ) THEN 3499 ts_value(22,sr) = hom(nzb,1,114,sr) ! L 3507 3500 ELSE 3508 ts_value(22,sr) = 1 0000.0_wp3501 ts_value(22,sr) = 1.0E10_wp 3509 3502 ENDIF 3510 3503 -
palm/trunk/SOURCE/init_1d_model.f90
r1698 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Set initial time step to 10 s to avoid instability of the 1d model for small 22 ! grid spacings 22 23 ! 23 24 ! Former revisions: … … 81 82 !> 82 83 !> @todo harmonize code with new surface_layer_fluxes module 84 !> @bug 1D model crashes when using small grid spacings in the order of 1 m 83 85 !------------------------------------------------------------------------------! 84 86 SUBROUTINE init_1d_model … … 262 264 !-- Determine the time step at the start of a 1D-simulation and 263 265 !-- determine and printout quantities used for run control 264 CALL timestep_1d266 dt_1d = 10.0_wp 265 267 CALL run_control_1d 266 268 … … 864 866 USE pegrid 865 867 866 USE control_parameters, 868 USE control_parameters, & 867 869 ONLY: message_string 868 870 -
palm/trunk/SOURCE/land_surface_model.f90
r1698 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Renamed pt_1 and qv_1 to pt1 and qv1. 22 ! Bugfix: set initial values for t_surface_p in case of restart runs 23 ! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky' 24 ! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky' 25 ! Added todo action 22 26 ! 23 27 ! Former revisions: … … 87 91 !> 88 92 !> @todo Consider partial absorption of the net shortwave radiation by the 89 !> s urfacelayer.93 !> skin layer. 90 94 !> @todo Allow for water surfaces 91 95 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, … … 93 97 !> @todo Implement surface runoff model (required when performing long-term LES 94 98 !> with considerable precipitation. 99 !> @todo Fix crashes with radiation_scheme == 'constant' 95 100 !> 96 !> @note No time step criterion is required as long as the soil layers do not97 !> 101 !> @note No time step criterion is required as long as the soil layers do not 102 !> become too thin. 98 103 !------------------------------------------------------------------------------! 99 104 MODULE land_surface_model_mod … … 124 129 USE radiation_model_mod, & 125 130 ONLY: force_radiation_call, radiation_scheme, rad_net, rad_sw_in, & 126 rad_lw_out, sigma_sb131 rad_lw_out, rad_lw_out_change_0, sigma_sb 127 132 128 133 #if defined ( __rrtmg ) 129 134 USE radiation_model_mod, & 130 ONLY: rrtm_ lwuflx_dt, rrtm_idrv135 ONLY: rrtm_idrv 131 136 #endif 132 137 … … 275 280 rad_net_l, & !< local copy of rad_net (net radiation at surface) 276 281 r_a, & !< aerodynamic resistance 277 r_a_av, & !< aver gae of r_a282 r_a_av, & !< average of r_a 278 283 r_canopy, & !< canopy resistance 279 284 r_soil, & !< soil resistance … … 635 640 INTEGER(iwp) :: k !< running index 636 641 637 REAL(wp) :: pt _1 !< potential temperature at first grid level642 REAL(wp) :: pt1 !< potential temperature at first grid level 638 643 639 644 … … 788 793 !-- Calculate surface temperature 789 794 t_surface = pt_surface * exn 790 t_surface_p = t_surface791 795 792 796 ! … … 798 802 799 803 IF ( cloud_physics ) THEN 800 pt _1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)804 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 801 805 ELSE 802 pt _1 = pt(k+1,j,i)806 pt1 = pt(k+1,j,i) 803 807 ENDIF 804 808 805 809 ! 806 810 !-- Assure that r_a cannot be zero at model start 807 IF ( pt _1 == pt(k,j,i) ) pt_1 = pt_1 + 1.0E-10_wp811 IF ( pt1 == pt(k,j,i) ) pt1 = pt1 + 1.0E-10_wp 808 812 809 813 us(j,i) = 0.1_wp 810 ts(j,i) = (pt _1 - pt(k,j,i)) / r_a(j,i)814 ts(j,i) = (pt1 - pt(k,j,i)) / r_a(j,i) 811 815 shf(j,i) = - us(j,i) * ts(j,i) 812 816 ENDDO … … 899 903 m_soil_p = m_soil 900 904 m_liq_eb_p = m_liq_eb 905 t_surface_p = t_surface 906 907 901 908 902 909 !-- Store initial profiles of t_soil and m_soil (assuming they are … … 926 933 dots_unit(dots_soil+6:dots_soil+7) = "s/m" 927 934 928 RETURN929 935 930 936 END SUBROUTINE init_lsm … … 964 970 lambda_surface, & !< Current value of lambda_surface 965 971 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 966 pt _1,& !< potential temperature at first grid level967 qv _1!< specific humidity at first grid level972 pt1, & !< potential temperature at first grid level 973 qv1 !< specific humidity at first grid level 968 974 969 975 ! … … 989 995 !-- equivalent to the ECMWF formulation using drag coefficients 990 996 IF ( cloud_physics ) THEN 991 pt _1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)992 qv _1 = q(k+1,j,i) - ql(k+1,j,i)997 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 998 qv1 = q(k+1,j,i) - ql(k+1,j,i) 993 999 ELSE 994 pt_1 = pt(k+1,j,i) 995 qv_1 = q(k+1,j,i) 996 ENDIF 997 998 r_a(j,i) = (pt_1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1000 pt1 = pt(k+1,j,i) 1001 qv1 = q(k+1,j,i) 1002 ENDIF 1003 1004 r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 1005 1006 ! 1007 !-- Make sure that the resistance does not drop to zero 1008 IF ( ABS(r_a(j,i)) < 1.0E-10_wp ) r_a(j,i) = 1.0E-10_wp 999 1009 1000 1010 ! … … 1005 1015 !-- night) 1006 1016 IF ( radiation_scheme /= 'constant' ) THEN 1007 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / &1008 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) &1017 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) / & 1018 (0.81_wp * (0.004_wp * rad_sw_in(k,j,i) & 1009 1019 + 1.0_wp)) ) 1010 1020 ELSE 1011 1021 f1 = 1.0_wp 1012 1022 ENDIF 1023 1013 1024 1014 1025 ! … … 1040 1051 ! 1041 1052 !-- Calculate vapour pressure 1042 e = qv _1 * surface_pressure / 0.622_wp1053 e = qv1 * surface_pressure / 0.622_wp 1043 1054 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1044 1055 ELSE … … 1082 1093 !-- In case of dewfall, set evapotranspiration to zero 1083 1094 !-- All super-saturated water is then removed from the air 1084 IF ( humidity .AND. q_s <= qv _1 ) THEN1095 IF ( humidity .AND. q_s <= qv1 ) THEN 1085 1096 r_canopy(j,i) = 0.0_wp 1086 1097 r_soil(j,i) = 0.0_wp … … 1108 1119 !-- reservoir are not allowed to take up water from the super-saturated 1109 1120 !-- air. 1110 IF ( humidity ) THEN 1111 IF ( q_s <= qv_1 ) THEN 1112 IF ( .NOT. dewfall ) THEN 1113 f_qsws_veg = 0.0_wp 1114 f_qsws_soil = 0.0_wp 1115 f_qsws_liq = 0.0_wp 1116 ENDIF 1117 ENDIF 1121 IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall ) THEN 1122 f_qsws_veg = 0.0_wp 1123 f_qsws_soil = 0.0_wp 1124 f_qsws_liq = 0.0_wp 1118 1125 ENDIF 1119 1126 … … 1133 1140 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1134 1141 1135 1136 1142 IF ( humidity ) THEN 1137 1143 #if defined ( __rrtmg ) 1138 1144 ! 1139 1145 !-- Numerator of the prognostic equation 1140 coef_1 = rad_net_l(j,i) + r rtm_lwuflx_dt(0,nzb)&1146 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1141 1147 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1142 + f_shf * pt _1 + f_qsws * ( qv_1 - q_s&1148 + f_shf * pt1 + f_qsws * ( qv1 - q_s & 1143 1149 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1144 1150 * t_soil(nzb_soil,j,i) … … 1146 1152 ! 1147 1153 !-- Denominator of the prognostic equation 1148 coef_2 = r rtm_lwuflx_dt(0,nzb) + f_qsws * dq_s_dt&1154 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt & 1149 1155 + lambda_surface + f_shf / exn 1150 1156 #else … … 1154 1160 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1155 1161 * t_surface(j,i) ** 4 & 1156 + f_shf * pt _1 + f_qsws * ( qv_1&1162 + f_shf * pt1 + f_qsws * ( qv1 & 1157 1163 - q_s + dq_s_dt * t_surface(j,i) ) & 1158 1164 + lambda_surface * t_soil(nzb_soil,j,i) 1159 1165 1160 1166 ! 1161 !-- 1162 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws&1163 1167 !-- Denominator of the prognostic equation 1168 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1169 * dq_s_dt + lambda_surface + f_shf / exn 1164 1170 1165 1171 #endif … … 1169 1175 ! 1170 1176 !-- Numerator of the prognostic equation 1171 coef_1 = rad_net_l(j,i) + r rtm_lwuflx_dt(0,nzb)&1177 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1172 1178 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1173 + f_shf * pt _1 + lambda_surface&1179 + f_shf * pt1 + lambda_surface & 1174 1180 * t_soil(nzb_soil,j,i) 1175 1181 1176 1182 ! 1177 1183 !-- Denominator of the prognostic equation 1178 coef_2 = r rtm_lwuflx_dt(0,nzb) + lambda_surface + f_shf / exn1184 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn 1179 1185 #else 1180 1186 … … 1182 1188 !-- Numerator of the prognostic equation 1183 1189 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1184 * t_surface(j,i) ** 4 + f_shf * pt _1&1190 * t_surface(j,i) ** 4 + f_shf * pt1 & 1185 1191 + lambda_surface * t_soil(nzb_soil,j,i) 1186 1192 … … 1191 1197 #endif 1192 1198 ENDIF 1193 1194 1199 1195 1200 tend = 0.0_wp … … 1238 1243 !-- Calculate fluxes 1239 1244 #if defined ( __rrtmg ) 1240 rad_net_l(j,i) = rad_net_l(j,i) + r rtm_lwuflx_dt(0,nzb)&1245 rad_net_l(j,i) = rad_net_l(j,i) + rad_lw_out_change_0(j,i) & 1241 1246 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1242 - r rtm_lwuflx_dt(0,nzb) * t_surface_p(j,i)1247 - rad_lw_out_change_0(j,i) * t_surface_p(j,i) 1243 1248 1244 1249 IF ( rrtm_idrv == 1 ) THEN 1245 1250 rad_net(j,i) = rad_net_l(j,i) 1246 1251 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & 1247 + r rtm_lwuflx_dt(0,nzb)&1252 + rad_lw_out_change_0(j,i) & 1248 1253 * ( t_surface_p(j,i) - t_surface(j,i) ) 1249 1254 ENDIF … … 1252 1257 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1253 1258 * t_surface(j,i)**3 * t_surface_p(j,i) 1254 rad_net(j,i) = rad_net_l(j,i)1255 1259 #endif 1260 1256 1261 1257 1262 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1258 1263 - t_soil(nzb_soil,j,i)) 1259 1264 1260 shf_eb(j,i) = - f_shf * ( pt _1 - pt(k,j,i) )1265 shf_eb(j,i) = - f_shf * ( pt1 - pt(k,j,i) ) 1261 1266 1262 1267 shf(j,i) = shf_eb(j,i) / rho_cp 1263 1268 1264 1269 IF ( humidity ) THEN 1265 qsws_eb(j,i) = - f_qsws * ( qv _1 - q_s + dq_s_dt&1270 qsws_eb(j,i) = - f_qsws * ( qv1 - q_s + dq_s_dt & 1266 1271 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1267 1272 1268 1273 qsws(j,i) = qsws_eb(j,i) / rho_lv 1269 1274 1270 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv _1 - q_s&1275 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv1 - q_s & 1271 1276 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1272 1277 * t_surface_p(j,i) ) 1273 1278 1274 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv _1 - q_s&1279 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s & 1275 1280 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1276 1281 * t_surface_p(j,i) ) 1277 1282 1278 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv _1 - q_s&1283 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv1 - q_s & 1279 1284 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1280 1285 * t_surface_p(j,i) ) … … 1285 1290 r_s(j,i) = 1.0E10_wp 1286 1291 ELSE 1287 r_s(j,i) = - rho_lv * ( qv _1 - q_s + dq_s_dt &1292 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1288 1293 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1289 1294 / qsws_eb(j,i) - r_a(j,i) -
palm/trunk/SOURCE/radiation_model.f90
r1702 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting 22 ! corrections 22 23 ! 23 24 ! Former revisions: … … 197 198 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 198 199 alpha, & !< surface broadband albedo (used for clear-sky scheme) 200 rad_lw_out_change_0, & !< change in LW out due to change in surface temperature 199 201 rad_net, & !< net radiation at the surface 200 202 rad_net_av !< average of rad_net … … 363 365 radiation_clearsky, radiation_rrtmg, radiation_scheme, & 364 366 radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out, & 365 rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, & 366 rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av, & 367 rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, & 368 skip_time_do_radiation, sw_radiation, time_radiation, time_utc_init 367 rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr, rad_lw_cs_hr_av, & 368 rad_lw_hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, & 369 rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, & 370 rad_sw_hr_av, sigma_sb, skip_time_do_radiation, sw_radiation, & 371 time_radiation, time_utc_init 369 372 370 373 371 374 #if defined ( __rrtmg ) 372 PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv, & 373 rrtm_lwuflx_dt 375 PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv 374 376 #endif 375 377 … … 390 392 ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) ) 391 393 rad_net = 0.0_wp 394 ENDIF 395 396 ! 397 !-- Allocate array for storing the surface net radiation 398 IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) ) THEN 399 ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) ) 400 rad_lw_out_change_0 = 0.0_wp 392 401 ENDIF 393 402 … … 654 663 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 655 664 IF ( radiation_scheme == 'clear-sky' ) THEN 656 665 CALL radiation_clearsky 657 666 ELSEIF ( radiation_scheme == 'rrtmg' ) THEN 658 667 CALL radiation_rrtmg … … 679 688 INTEGER(iwp) :: i, j, k !< loop indices 680 689 REAL(wp) :: exn, & !< Exner functions at surface 681 exn _1,& !< Exner functions at first grid level682 pt _1!< potential temperature at first grid level690 exn1, & !< Exner functions at first grid level 691 pt1 !< potential temperature at first grid level 683 692 684 693 ! … … 696 705 !-- Calculate radiation fluxes and net radiation (rad_net) for each grid 697 706 !-- point 698 DO i = nxl , nxr699 DO j = nys , nyn707 DO i = nxlg, nxrg 708 DO j = nysg, nyng 700 709 k = nzb_s_inner(j,i) 701 710 702 exn _1 = (hyp(k+1) / 100000.0_wp )**0.286_wp711 exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp 703 712 704 713 rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0) … … 707 716 708 717 IF ( cloud_physics ) THEN 709 pt _1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)710 rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt _1 * exn_1)**4718 pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i) 719 rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4 711 720 ELSE 712 rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn _1)**4721 rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4 713 722 ENDIF 714 723 … … 718 727 ENDDO 719 728 ENDDO 720 721 CALL exchange_horiz_2d( rad_lw_in, nbgp )722 CALL exchange_horiz_2d( rad_lw_out, nbgp )723 CALL exchange_horiz_2d( rad_sw_in, nbgp )724 CALL exchange_horiz_2d( rad_sw_out, nbgp )725 CALL exchange_horiz_2d( rad_net, nbgp )726 727 RETURN728 729 729 730 END SUBROUTINE radiation_clearsky … … 901 902 rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day 902 903 ENDDO 904 905 ! 906 !-- Save change in LW heating rate 907 rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb) 903 908 904 909 ENDIF … … 953 958 954 959 CALL exchange_horiz_2d( rad_net, nbgp ) 960 CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp ) 955 961 #endif 956 962 … … 1347 1353 ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) ) 1348 1354 1349 rrtm_lwuflx c_dt = 0.0_wp1355 rrtm_lwuflx_dt = 0.0_wp 1350 1356 rrtm_lwuflxc_dt = 0.0_wp 1351 1357 -
palm/trunk/SOURCE/read_3d_binary.f90
r1692 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added rad_lw_out_change_0, increased binary_version 22 22 ! 23 23 ! Former revisions: … … 124 124 ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out, & 125 125 rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, & 126 rad_lw_out_av, rad_ sw_in, rad_sw_in_av, rad_sw_out,&127 rad_sw_out _av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,&128 rad_sw_hr _av126 rad_lw_out_av, rad_lw_out_change_0, rad_sw_in, rad_sw_in_av, & 127 rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, & 128 rad_sw_hr, rad_sw_hr_av 129 129 130 130 USE random_function_mod, & … … 348 348 !-- First compare the version numbers 349 349 READ ( 13 ) version_on_file 350 binary_version = '4. 1'350 binary_version = '4.2' 351 351 IF ( TRIM( version_on_file ) /= TRIM( binary_version ) ) THEN 352 352 WRITE( message_string, * ) 'version mismatch concerning data ', & … … 877 877 rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 878 878 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 879 880 CASE ( 'rad_lw_out_change_0' ) 881 IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) ) THEN 882 ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) ) 883 ENDIF 884 IF ( k == 1 ) READ ( 13 ) tmp_2d 885 rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)& 886 = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 879 887 880 888 CASE ( 'rad_lw_cs_hr' ) -
palm/trunk/SOURCE/surface_layer_fluxes.f90
r1706 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix: division by zero could occur when calculating rib at low wind speeds 22 ! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for 23 ! neutral = .T. 22 24 ! 23 25 ! Former revisions: … … 96 98 !> that this method cannot be used in case of roughness heterogeneity 97 99 !> 98 !> @todo limiting of uv_total might be required in some cases99 100 !> @todo (re)move large_scale_forcing actions 100 101 !> @todo check/optimize OpenMP and OpenACC directives … … 208 209 CALL calc_scaling_parameters 209 210 211 CALL calc_uv_total 212 210 213 IF ( .NOT. neutral ) THEN 211 214 CALL calc_ol … … 220 223 !-- number to calculate the Obukhov length 221 224 ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN 225 226 CALL calc_uv_total 222 227 223 228 IF ( .NOT. neutral ) THEN … … 273 278 !-- Allocate field for storing the horizontal velocity 274 279 ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) ) 280 281 282 ! 283 !-- In case of runs with neutral statification, set Obukhov length to a 284 !-- large value 285 IF ( neutral ) ol = 1.0E10_wp 275 286 276 287 IF ( most_method == 'lookup' ) THEN … … 384 395 ! Description: 385 396 ! ------------ 397 !> Compute the absolute value of the horizontal velocity (relative to the 398 !> surface). This is required by all methods 399 !------------------------------------------------------------------------------! 400 SUBROUTINE calc_uv_total 401 402 IMPLICIT NONE 403 404 405 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 406 !$acc kernels loop 407 DO i = nxl, nxr 408 DO j = nys, nyn 409 410 k = nzb_s_inner(j,i) 411 uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) & 412 - u(k,j,i) - u(k,j,i+1) ) )**2 + & 413 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) & 414 - v(k,j,i) - v(k,j+1,i) ) )**2 ) 415 416 ! 417 !-- For too small values of the local wind, MOST does not work. A 418 !-- threshold value is thus set if required 419 ! uv_total(j,i) = MAX(0.01_wp,uv_total(j,i)) 420 421 ENDDO 422 ENDDO 423 424 ! 425 !-- Values of uv_total need to be exchanged at the ghost boundaries 426 !$acc update host( uv_total ) 427 CALL exchange_horiz_2d( uv_total ) 428 !$acc update device( uv_total ) 429 430 END SUBROUTINE calc_uv_total 431 432 433 !------------------------------------------------------------------------------! 434 ! Description: 435 ! ------------ 386 436 !> Calculate the Obukhov length (L) and Richardson flux number (z/L) 387 437 !------------------------------------------------------------------------------! … … 394 444 395 445 REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number 396 446 397 447 REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0 398 448 f_d_ol, & !< Derivative of f 399 449 ol_l, & !< Lower bound of L for Newton iteration 400 450 ol_m, & !< Previous value of L for Newton iteration 401 ol_old, & !< Previous time step value of L 451 ol_old, & !< Previous time step value of L 402 452 ol_u !< Upper bound of L for Newton iteration 403 453 404 !405 !-- Compute the absolute value of the horizontal velocity (relative to the406 !-- surface). This is required by all methods407 !$OMP PARALLEL DO PRIVATE( k, z_mo )408 !$acc kernels loop409 DO i = nxl, nxr410 DO j = nys, nyn411 412 k = nzb_s_inner(j,i)413 414 uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) &415 - u(k,j,i) - u(k,j,i+1) ) )**2 + &416 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) &417 - v(k,j,i) - v(k,j+1,i) ) )**2 )418 419 !420 !-- For too small values of the local wind, MOST does not work. A421 !-- threshold value is thus set if required422 ! uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))423 424 ENDDO425 ENDDO426 427 !428 !-- Values of uv_total need to be exchanged at the ghost boundaries429 !$acc update host( uv_total )430 CALL exchange_horiz_2d( uv_total )431 !$acc update device( uv_total )432 454 433 455 IF ( TRIM( most_method ) /= 'circular' ) THEN … … 447 469 IF ( humidity ) THEN 448 470 rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) ) & 449 / ( uv_total(j,i)**2 * vpt(k+1,j,i) )471 / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp ) 450 472 ELSE 451 473 rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) ) & 452 / ( uv_total(j,i)**2 * pt(k+1,j,i) )474 / ( uv_total(j,i)**2 * pt(k+1,j,i) + 1.0E-20_wp ) 453 475 ENDIF 454 476 ELSE … … 462 484 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp & 463 485 * pt(k+1,j,i) * qsws(j,i) ) & 464 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2 ) 486 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2& 487 + 1.0E-20_wp) 465 488 ELSE 466 489 rib(j,i) = - g * z_mo * shf(j,i) & 467 / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 ) 490 / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 & 491 + 1.0E-20_wp ) 468 492 ENDIF 469 493 ENDIF -
palm/trunk/SOURCE/write_3d_binary.f90
r1692 r1709 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added rad_lw_out_change_0, increased binary_version 22 22 ! 23 23 ! Former revisions: … … 114 114 USE radiation_model_mod, & 115 115 ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, & 116 rad_lw_out, rad_lw_out_av, rad_lw_ cs_hr, rad_lw_cs_hr_av,&117 rad_lw_ hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,&118 rad_sw_ out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,&119 rad_sw_ hr_av116 rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr, & 117 rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, & 118 rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, & 119 rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av 120 120 121 121 USE random_function_mod, & … … 139 139 ! 140 140 !-- Write arrays. 141 binary_version = '4. 1'141 binary_version = '4.2' 142 142 143 143 WRITE ( 14 ) binary_version … … 297 297 WRITE ( 14 ) 'rad_lw_out_av '; WRITE ( 14 ) rad_lw_out_av 298 298 ENDIF 299 IF ( ALLOCATED( rad_lw_out_change_0 ) ) THEN 300 WRITE ( 14 ) 'rad_lw_out_change_0 ' 301 WRITE ( 14 ) rad_lw_out_change_0 302 ENDIF 299 303 IF ( ALLOCATED( rad_lw_cs_hr ) ) THEN 300 304 WRITE ( 14 ) 'rad_lw_cs_hr '; WRITE ( 14 ) rad_lw_cs_hr
Note: See TracChangeset
for help on using the changeset viewer.
