Changeset 4366 for palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
- Timestamp:
- Jan 9, 2020 8:12:43 AM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r4360 r4366 26 26 ! ----------------- 27 27 ! $Id$ 28 ! vector version for calculation of Obukhov length via Newton iteration added 29 ! 30 ! 4360 2020-01-07 11:25:50Z suehring 28 31 ! Calculation of diagnostic-only 2-m potential temperature moved to 29 32 ! diagnostic_output_quantities. … … 112 115 constant_waterflux, coupling_mode, & 113 116 debug_output_timestep, & 114 humidity, 117 humidity, loop_optimization, & 115 118 ibc_e_b, ibc_pt_b, indoor_model, & 116 119 land_surface, large_scale_forcing, lsf_surf, message_string, & … … 859 862 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 860 863 864 LOGICAL, DIMENSION(surf%ns) :: convergence_reached !< convergence switch for vectorization 865 861 866 REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0 862 867 f_d_ol, & !< Derivative of f … … 866 871 ol_u !< Upper bound of L for Newton iteration 867 872 873 REAL(wp), DIMENSION(surf%ns) :: ol_old_vec !< temporary array required for vectorization 874 REAL(wp), DIMENSION(surf%ns) :: z_mo_vec !< temporary array required for vectorization 875 868 876 ! 869 877 !-- Evaluate bulk Richardson number (calculation depends on … … 926 934 ENDIF 927 935 928 ! 929 !-- Calculate the Obukhov length using Newton iteration 930 !$OMP PARALLEL DO PRIVATE(i, j, z_mo) & 931 !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) 932 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) & 933 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) & 934 !$ACC PRESENT(surf) 935 DO m = 1, surf%ns 936 937 i = surf%i(m) 938 j = surf%j(m) 939 940 z_mo = surf%z_mo(m) 941 942 ! 943 !-- Store current value in case the Newton iteration fails 944 ol_old = surf%ol(m) 945 946 ! 947 !-- Ensure that the bulk Richardson number and the Obukhov 948 !-- length have the same sign 949 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. & 950 ABS( surf%ol(m) ) == ol_max ) THEN 951 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 952 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 953 ENDIF 936 IF ( loop_optimization == 'cache' ) THEN 937 ! 938 !-- Calculate the Obukhov length using Newton iteration 939 !$OMP PARALLEL DO PRIVATE(i, j, z_mo) & 940 !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) 941 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) & 942 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) & 943 !$ACC PRESENT(surf) 944 DO m = 1, surf%ns 945 946 i = surf%i(m) 947 j = surf%j(m) 948 949 z_mo = surf%z_mo(m) 950 951 ! 952 !-- Store current value in case the Newton iteration fails 953 ol_old = surf%ol(m) 954 955 ! 956 !-- Ensure that the bulk Richardson number and the Obukhov 957 !-- length have the same sign 958 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN 959 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 960 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 961 ENDIF 962 ! 963 !-- Iteration to find Obukhov length 964 iter = 0 965 DO 966 iter = iter + 1 967 ! 968 !-- In case of divergence, use the value of the previous time step 969 IF ( iter > 1000 ) THEN 970 surf%ol(m) = ol_old 971 EXIT 972 ENDIF 973 974 ol_m = surf%ol(m) 975 ol_l = ol_m - 0.001_wp * ol_m 976 ol_u = ol_m + 0.001_wp * ol_m 977 978 979 IF ( ibc_pt_b /= 1 ) THEN 980 ! 981 !-- Calculate f = Ri - [...]/[...]^2 = 0 982 f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) ) & 983 - psi_h( z_mo / ol_m ) & 984 + psi_h( surf%z0h(m) / ol_m ) & 985 ) / & 986 ( LOG( z_mo / surf%z0(m) ) & 987 - psi_m( z_mo / ol_m ) & 988 + psi_m( surf%z0(m) / ol_m ) & 989 )**2 990 991 ! 992 !-- Calculate df/dL 993 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) ) & 994 - psi_h( z_mo / ol_u ) & 995 + psi_h( surf%z0h(m) / ol_u ) & 996 ) / & 997 ( LOG( z_mo / surf%z0(m) ) & 998 - psi_m( z_mo / ol_u ) & 999 + psi_m( surf%z0(m) / ol_u ) & 1000 )**2 & 1001 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) & 1002 - psi_h( z_mo / ol_l ) & 1003 + psi_h( surf%z0h(m) / ol_l ) & 1004 ) / & 1005 ( LOG( z_mo / surf%z0(m) ) & 1006 - psi_m( z_mo / ol_l ) & 1007 + psi_m( surf%z0(m) / ol_l ) & 1008 )**2 & 1009 ) / ( ol_u - ol_l ) 1010 ELSE 1011 ! 1012 !-- Calculate f = Ri - 1 /[...]^3 = 0 1013 f = surf%rib(m) - ( z_mo / ol_m ) / ( LOG( z_mo / surf%z0(m) ) & 1014 - psi_m( z_mo / ol_m ) & 1015 + psi_m( surf%z0(m) / ol_m ) & 1016 )**3 1017 1018 ! 1019 !-- Calculate df/dL 1020 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) & 1021 - psi_m( z_mo / ol_u ) & 1022 + psi_m( surf%z0(m) / ol_u ) & 1023 )**3 & 1024 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1025 - psi_m( z_mo / ol_l ) & 1026 + psi_m( surf%z0(m) / ol_l ) & 1027 )**3 & 1028 ) / ( ol_u - ol_l ) 1029 ENDIF 1030 ! 1031 !-- Calculate new L 1032 surf%ol(m) = ol_m - f / f_d_ol 1033 1034 ! 1035 !-- Ensure that the bulk Richardson number and the Obukhov 1036 !-- length have the same sign and ensure convergence. 1037 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1038 ! 1039 !-- If unrealistic value occurs, set L to the maximum 1040 !-- value that is allowed 1041 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1042 surf%ol(m) = ol_max 1043 EXIT 1044 ENDIF 1045 ! 1046 !-- Assure that Obukhov length does not become zero. If the limit is 1047 !-- reached, exit the loop. 1048 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1049 surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) ) 1050 EXIT 1051 ENDIF 1052 ! 1053 !-- Check for convergence 1054 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT 1055 1056 ENDDO 1057 ENDDO 1058 1059 ! 1060 !-- Vector Version 1061 ELSE 1062 ! 1063 !-- Calculate the Obukhov length using Newton iteration 1064 !-- First set arrays required for vectorization 1065 DO m = 1, surf%ns 1066 1067 z_mo_vec(m) = surf%z_mo(m) 1068 1069 ! 1070 !-- Store current value in case the Newton iteration fails 1071 ol_old_vec(m) = surf%ol(m) 1072 1073 ! 1074 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign 1075 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN 1076 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 1077 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 1078 ENDIF 1079 1080 ENDDO 1081 954 1082 ! 955 1083 !-- Iteration to find Obukhov length 1084 convergence_reached(:) = .FALSE. 956 1085 iter = 0 957 1086 DO 1087 958 1088 iter = iter + 1 959 1089 ! 960 !-- In case of divergence, use the value of the previous time step1090 !-- In case of divergence, use the value(s) of the previous time step 961 1091 IF ( iter > 1000 ) THEN 962 surf%ol(m) = ol_old 1092 DO m = 1, surf%ns 1093 IF ( .NOT. convergence_reached(m) ) surf%ol(1:surf%ns) = ol_old 1094 ENDDO 963 1095 EXIT 964 1096 ENDIF 965 1097 966 ol_m = surf%ol(m) 967 ol_l = ol_m - 0.001_wp * ol_m 968 ol_u = ol_m + 0.001_wp * ol_m 969 970 971 IF ( ibc_pt_b /= 1 ) THEN 972 ! 973 !-- Calculate f = Ri - [...]/[...]^2 = 0 974 f = surf%rib(m) - ( z_mo / ol_m ) * ( & 975 LOG( z_mo / surf%z0h(m) ) & 976 - psi_h( z_mo / ol_m ) & 977 + psi_h( surf%z0h(m) / & 978 ol_m ) & 979 ) & 980 / ( LOG( z_mo / surf%z0(m) ) & 981 - psi_m( z_mo / ol_m ) & 982 + psi_m( surf%z0(m) / ol_m ) & 983 )**2 984 985 ! 986 !-- Calculate df/dL 987 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / & 988 surf%z0h(m) ) & 989 - psi_h( z_mo / ol_u ) & 990 + psi_h( surf%z0h(m) / ol_u ) & 991 ) & 992 / ( LOG( z_mo / surf%z0(m) ) & 993 - psi_m( z_mo / ol_u ) & 994 + psi_m( surf%z0(m) / ol_u ) & 995 )**2 & 996 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) & 997 - psi_h( z_mo / ol_l ) & 998 + psi_h( surf%z0h(m) / ol_l ) & 999 ) & 1000 / ( LOG( z_mo / surf%z0(m) ) & 1001 - psi_m( z_mo / ol_l ) & 1002 + psi_m( surf%z0(m) / ol_l ) & 1003 )**2 & 1004 ) / ( ol_u - ol_l ) 1005 ELSE 1006 ! 1007 !-- Calculate f = Ri - 1 /[...]^3 = 0 1008 f = surf%rib(m) - ( z_mo / ol_m ) / & 1009 ( LOG( z_mo / surf%z0(m) ) & 1010 - psi_m( z_mo / ol_m ) & 1011 + psi_m( surf%z0(m) / ol_m ) & 1012 )**3 1013 1014 ! 1015 !-- Calculate df/dL 1016 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) & 1017 - psi_m( z_mo / ol_u ) & 1018 + psi_m( surf%z0(m) / ol_u ) & 1019 )**3 & 1020 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1021 - psi_m( z_mo / ol_l ) & 1022 + psi_m( surf%z0(m) / ol_l ) & 1023 )**3 & 1024 ) / ( ol_u - ol_l ) 1025 ENDIF 1026 ! 1027 !-- Calculate new L 1028 surf%ol(m) = ol_m - f / f_d_ol 1029 1030 ! 1031 !-- Ensure that the bulk Richardson number and the Obukhov 1032 !-- length have the same sign and ensure convergence. 1033 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1034 ! 1035 !-- If unrealistic value occurs, set L to the maximum 1036 !-- value that is allowed 1037 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1038 surf%ol(m) = ol_max 1039 EXIT 1040 ENDIF 1041 ! 1042 !-- Assure that Obukhov length does not become zero. If the limit is 1043 !-- reached, exit the loop. 1044 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1045 surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) ) 1046 EXIT 1047 ENDIF 1048 ! 1049 !-- Check for convergence 1050 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN 1051 EXIT 1052 ELSE 1053 CYCLE 1054 ENDIF 1055 1056 ENDDO 1057 ENDDO 1098 1099 DO m = 1, surf%ns 1100 1101 IF ( convergence_reached(m) ) CYCLE 1102 1103 ol_m = surf%ol(m) 1104 ol_l = ol_m - 0.001_wp * ol_m 1105 ol_u = ol_m + 0.001_wp * ol_m 1106 1107 1108 IF ( ibc_pt_b /= 1 ) THEN 1109 ! 1110 !-- Calculate f = Ri - [...]/[...]^2 = 0 1111 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1112 - psi_h( z_mo_vec(m) / ol_m ) & 1113 + psi_h( surf%z0h(m) / ol_m ) & 1114 ) / & 1115 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1116 - psi_m( z_mo_vec(m) / ol_m ) & 1117 + psi_m( surf%z0(m) / ol_m ) & 1118 )**2 1119 1120 ! 1121 !-- Calculate df/dL 1122 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1123 - psi_h( z_mo_vec(m) / ol_u ) & 1124 + psi_h( surf%z0h(m) / ol_u ) & 1125 ) / & 1126 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1127 - psi_m( z_mo_vec(m) / ol_u ) & 1128 + psi_m( surf%z0(m) / ol_u ) & 1129 )**2 & 1130 + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1131 - psi_h( z_mo_vec(m) / ol_l ) & 1132 + psi_h( surf%z0h(m) / ol_l ) & 1133 ) / & 1134 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1135 - psi_m( z_mo_vec(m) / ol_l ) & 1136 + psi_m( surf%z0(m) / ol_l ) & 1137 )**2 & 1138 ) / ( ol_u - ol_l ) 1139 ELSE 1140 ! 1141 !-- Calculate f = Ri - 1 /[...]^3 = 0 1142 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1143 - psi_m( z_mo_vec(m) / ol_m ) & 1144 + psi_m( surf%z0(m) / ol_m ) & 1145 )**3 1146 1147 ! 1148 !-- Calculate df/dL 1149 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1150 - psi_m( z_mo_vec(m) / ol_u ) & 1151 + psi_m( surf%z0(m) / ol_u ) & 1152 )**3 & 1153 + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1154 - psi_m( z_mo_vec(m) / ol_l ) & 1155 + psi_m( surf%z0(m) / ol_l ) & 1156 )**3 & 1157 ) / ( ol_u - ol_l ) 1158 ENDIF 1159 ! 1160 !-- Calculate new L 1161 surf%ol(m) = ol_m - f / f_d_ol 1162 1163 ! 1164 !-- Ensure that the bulk Richardson number and the Obukhov 1165 !-- length have the same sign and ensure convergence. 1166 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1167 1168 ! 1169 !-- Check for convergence 1170 !-- This check does not modify surf%ol, therefore this is done first 1171 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN 1172 convergence_reached(m) = .TRUE. 1173 ENDIF 1174 ! 1175 !-- If unrealistic value occurs, set L to the maximum allowed value 1176 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1177 surf%ol(m) = ol_max 1178 convergence_reached(m) = .TRUE. 1179 ENDIF 1180 1181 ENDDO 1182 ! 1183 !-- Assure that Obukhov length does not become zero 1184 DO m = 1, surf%ns 1185 IF ( convergence_reached(m) ) CYCLE 1186 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1187 surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) ) 1188 convergence_reached(m) = .TRUE. 1189 ENDIF 1190 ENDDO 1191 1192 IF ( ALL( convergence_reached ) ) EXIT 1193 1194 ENDDO ! end of iteration loop 1195 1196 ENDIF ! end of vector branch 1058 1197 1059 1198 END SUBROUTINE calc_ol
Note: See TracChangeset
for help on using the changeset viewer.
