Changeset 3099
- Timestamp:
- 2011-11-14T17:35:33+01:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3081 r3099 1020 1020 !! 1021 1021 !! ** Method : s-coordinate case. 1022 !! Reformulate the horizontal hydrostatical pressure gradient1023 !! term using Pressure Jacobian. A new correction term1024 !! is developed to eliminate the sigma-coordinate error.1022 !! A Pressure-Jacobian horizontal pressure gradient method 1023 !! based on the constrained cubic-spline interpolation for 1024 !! all vertical coordinate systems 1025 1025 !! 1026 1026 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend … … 1030 1030 1031 1031 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1032 USE oce , ONLY: fsde3w_tmp=> ta ! (ta,sa) used as 3D workspace1033 USE oce , ONLY: rhd_tmp=> sa1032 USE oce , ONLY: zdeptht => ta ! (ta,sa) used as 3D workspace 1033 USE oce , ONLY: zrhh => sa 1034 1034 USE wrk_nemo, ONLY: zhpi => wrk_3d_3 1035 1035 USE wrk_nemo, ONLY: zu => wrk_3d_4 … … 1052 1052 !! The local variables for the correction term 1053 1053 INTEGER :: jk1, jis, jid, jjs, jjd 1054 REAL(wp) :: zuijk, zvijk, pwes, pwed, pnss, pnsd,deps1055 REAL(wp) :: rhdt11056 REAL(wp) :: dpdx1, dpdx2, dpdy1,dpdy21057 INTEGER :: bhitwe,bhitns1054 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 1055 REAL(wp) :: zrhdt1 1056 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 1057 INTEGER :: zbhitwe, zbhitns 1058 1058 !!---------------------------------------------------------------------- 1059 1059 … … 1074 1074 IF( lk_vvl ) znad = 1._wp 1075 1075 1076 ! Save fsde3w and rhd1077 fsde3w_tmp(:,:,:) = fsde3w(:,:,:)1078 rhd_tmp(:,:,:) = rhd(:,:,:)1079 1080 1076 ! Clean 3-D work arrays 1081 1077 zhpi(:,:,:) = 0. 1082 1078 1083 ! how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here1084 1079 ! Preparing vertical density profile for hybrid-sco coordinate 1085 1080 DO jj = 1, jpj 1086 1081 DO ji = 1, jpi 1087 1082 jk = mbathy(ji,jj) 1088 IF( jk <= 0 ) THEN; rhd(ji,jj,:) = 0._wp1089 ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)1083 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 1084 ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 1090 1085 ELSE IF(jk < jpkm1) THEN 1091 1086 DO jkk = jk+1, jpk 1092 rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), &1087 zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), & 1093 1088 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1094 1089 END DO 1095 END 1090 ENDIF 1096 1091 END DO 1097 1092 END DO … … 1099 1094 DO jj = 1, jpj 1100 1095 DO ji = 1, jpi 1101 fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)1102 fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad1096 zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 1097 zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 1103 1098 DO jk = 2, jpk 1104 fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)1099 zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 1105 1100 END DO 1106 1101 END DO … … 1110 1105 DO jj = 1, jpj 1111 1106 DO ji = 1, jpi 1112 fsp(ji,jj,jk) = rhd(ji,jj,jk)1113 xsp(ji,jj,jk) = fsde3w(ji,jj,jk)1107 fsp(ji,jj,jk) = zrhh(ji,jj,jk) 1108 xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 1114 1109 END DO 1115 1110 END DO … … 1123 1118 DO jj = 2, jpj 1124 1119 DO ji = 2, jpi 1125 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), &1120 zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), & 1126 1121 bsp(ji,jj,1), csp(ji,jj,1), & 1127 dsp(ji,jj,1) ) * 0.5_wp * fsde3w(ji,jj,1)1128 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water1122 dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 1123 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1129 1124 1130 1125 ! assuming linear profile across the top half surface layer 1131 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * rhdt11126 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1 1132 1127 END DO 1133 1128 END DO … … 1138 1133 DO ji = 2, jpi 1139 1134 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1140 integ2( fsde3w(ji,jj,jk-1), fsde3w(ji,jj,jk),&1135 integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),& 1141 1136 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), & 1142 1137 csp(ji,jj,jk-1), dsp(ji,jj,jk-1)) … … 1174 1169 DO jj = 2, jpjm1 1175 1170 DO ji = 2, jpim1 1176 pwes = 0._wp;pwed = 0._wp1177 pnss = 0._wp;pnsd = 0._wp1171 zpwes = 0._wp; zpwed = 0._wp 1172 zpnss = 0._wp; zpnsd = 0._wp 1178 1173 zuijk = zu(ji,jj,jk) 1179 1174 zvijk = zv(ji,jj,jk) … … 1181 1176 !!!!! for u equation 1182 1177 IF( jk <= mbku(ji,jj) ) THEN 1183 IF( - fsde3w(ji+1,jj,mbku(ji,jj)) >= -fsde3w(ji,jj,mbku(ji,jj)) ) THEN1178 IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN 1184 1179 jis = ji + 1; jid = ji 1185 1180 ELSE … … 1189 1184 ! integrate the pressure on the shallow side 1190 1185 jk1 = jk 1191 bhitwe = 01192 DO WHILE ( - fsde3w(jis,jj,jk1) > zuijk )1186 zbhitwe = 0 1187 DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 1193 1188 IF( jk1 == mbku(ji,jj) ) THEN 1194 bhitwe = 11189 zbhitwe = 1 1195 1190 EXIT 1196 1191 ENDIF 1197 deps = min(fsde3w(jis,jj,jk1+1), -zuijk)1198 pwes =pwes + &1199 integ2( fsde3w(jis,jj,jk1),deps, &1192 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk) 1193 zpwes = zpwes + & 1194 integ2(zdeptht(jis,jj,jk1), zdeps, & 1200 1195 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1201 1196 csp(jis,jj,jk1), dsp(jis,jj,jk1)) … … 1203 1198 END DO 1204 1199 1205 IF( bhitwe == 1) THEN1206 zuijk = - fsde3w(jis,jj,jk1)1200 IF(zbhitwe == 1) THEN 1201 zuijk = -zdeptht(jis,jj,jk1) 1207 1202 ENDIF 1208 1203 1209 1204 ! integrate the pressure on the deep side 1210 1205 jk1 = jk 1211 bhitwe = 01212 DO WHILE ( - fsde3w(jid,jj,jk1) < zuijk )1206 zbhitwe = 0 1207 DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 1213 1208 IF( jk1 == 1 ) THEN 1214 bhitwe = 11209 zbhitwe = 1 1215 1210 EXIT 1216 1211 ENDIF 1217 deps = max(fsde3w(jid,jj,jk1-1), -zuijk)1218 pwed =pwed + &1219 integ2( deps, fsde3w(jid,jj,jk1), &1212 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk) 1213 zpwed = zpwed + & 1214 integ2(zdeps, zdeptht(jid,jj,jk1), & 1220 1215 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1221 1216 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) … … 1223 1218 END DO 1224 1219 1225 IF( bhitwe == 1 ) THEN1226 deps = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad)1227 rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), &1220 IF( zbhitwe == 1 ) THEN 1221 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 1222 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), & 1228 1223 bsp(jid,jj,1), csp(jid,jj,1), & 1229 dsp(jid,jj,1)) * deps1230 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water1231 pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) *deps1224 dsp(jid,jj,1)) * zdeps 1225 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1226 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1232 1227 ENDIF 1233 1228 1234 1235 dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1229 ! update the momentum trends in u direction 1230 1231 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1236 1232 IF( lk_vvl ) THEN 1237 dpdx2 = zcoef0 / e1u(ji,jj) * &1238 ( REAL(jis-jid, wp) * ( pwes +pwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )1233 zdpdx2 = zcoef0 / e1u(ji,jj) * & 1234 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1239 1235 ELSE 1240 dpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (pwes +pwed)1236 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1241 1237 ENDIF 1242 1238 1243 ua(ji,jj,jk) = ua(ji,jj,jk) + ( dpdx1 +dpdx2) * &1239 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1244 1240 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1245 1241 ENDIF … … 1247 1243 !!!!! for v equation 1248 1244 IF( jk <= mbkv(ji,jj) ) THEN 1249 IF( - fsde3w(ji,jj+1,mbkv(ji,jj)) >= -fsde3w(ji,jj,mbkv(ji,jj)) ) THEN1245 IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN 1250 1246 jjs = jj + 1; jjd = jj 1251 1247 ELSE … … 1255 1251 ! integrate the pressure on the shallow side 1256 1252 jk1 = jk 1257 bhitns = 01258 DO WHILE ( - fsde3w(ji,jjs,jk1) > zvijk )1253 zbhitns = 0 1254 DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 1259 1255 IF( jk1 == mbkv(ji,jj) ) THEN 1260 bhitns = 11256 zbhitns = 1 1261 1257 EXIT 1262 1258 ENDIF 1263 deps = min(fsde3w(ji,jjs,jk1+1), -zvijk)1264 pnss =pnss + &1265 integ2( fsde3w(ji,jjs,jk1),deps, &1259 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk) 1260 zpnss = zpnss + & 1261 integ2(zdeptht(ji,jjs,jk1), zdeps, & 1266 1262 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1267 1263 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) … … 1269 1265 END DO 1270 1266 1271 IF( bhitns == 1) THEN1272 zvijk = - fsde3w(ji,jjs,jk1)1267 IF(zbhitns == 1) THEN 1268 zvijk = -zdeptht(ji,jjs,jk1) 1273 1269 ENDIF 1274 1270 1275 1271 ! integrate the pressure on the deep side 1276 1272 jk1 = jk 1277 bhitns = 01278 DO WHILE ( - fsde3w(ji,jjd,jk1) < zvijk )1273 zbhitns = 0 1274 DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 1279 1275 IF( jk1 == 1 ) THEN 1280 bhitns = 11276 zbhitns = 1 1281 1277 EXIT 1282 1278 ENDIF 1283 deps = max(fsde3w(ji,jjd,jk1-1), -zvijk)1284 pnsd =pnsd + &1285 integ2( deps, fsde3w(ji,jjd,jk1), &1279 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk) 1280 zpnsd = zpnsd + & 1281 integ2(zdeps, zdeptht(ji,jjd,jk1), & 1286 1282 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1287 1283 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) … … 1289 1285 END DO 1290 1286 1291 IF( bhitns == 1 ) THEN1292 deps = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad)1293 rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), &1287 IF( zbhitns == 1 ) THEN 1288 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 1289 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), & 1294 1290 bsp(ji,jjd,1), csp(ji,jjd,1), & 1295 dsp(ji,jjd,1) ) * deps1296 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water1297 pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) *deps1291 dsp(ji,jjd,1) ) * zdeps 1292 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1293 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1298 1294 ENDIF 1299 1295 1300 1301 dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1296 ! update the momentum trends in v direction 1297 1298 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1302 1299 IF( lk_vvl ) THEN 1303 dpdy2 = zcoef0 / e2v(ji,jj) * &1304 ( REAL(jjs-jjd, wp) * ( pnss +pnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )1300 zdpdy2 = zcoef0 / e2v(ji,jj) * & 1301 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1305 1302 ELSE 1306 dpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (pnss +pnsd )1303 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1307 1304 ENDIF 1308 1305 1309 va(ji,jj,jk) = va(ji,jj,jk) + ( dpdy1 +dpdy2)*&1306 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1310 1307 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1311 1308 ENDIF … … 1316 1313 END DO 1317 1314 1318 ! Restore fsde3w and rhd1319 fsde3w(:,:,:) = fsde3w_tmp(:,:,:)1320 rhd(:,:,:) = rhd_tmp(:,:,:)1321 1322 1315 ! 1323 1316 IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) ) & 1324 1317 CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 1325 1318 ! 1319 XXXXXXX 1326 1320 END SUBROUTINE hpg_prj 1327 1321 … … 1362 1356 dsp = 0.0_wp 1363 1357 1358 DO ji = 1, jpi 1359 DO jj = 1, jpj 1360 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1361 DO jk = 2, jpkm1-1 1362 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1363 zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1364 zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1365 zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1366 1367 zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1368 1369 IF(zdf1 * zdf2 <= 0._wp) THEN 1370 zdf(jk) = 0._wp 1371 ELSE 1372 zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1373 ENDIF 1374 END DO 1375 1376 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1377 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1378 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1379 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1380 & 0.5_wp * zdf(jpkm1 - 1) 1381 1382 DO jk = 1, jpkm1 - 1 1383 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1384 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1385 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1386 zddf1 = -2._wp * ztmp1 + ztmp2 1387 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1388 zddf2 = 2._wp * ztmp1 - ztmp2 1389 1390 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1391 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1392 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1393 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1394 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1395 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1396 & xsp(ji,jj,jk)**2 ) 1397 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1398 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1399 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1400 END DO 1401 1402 ELSE IF (polynomial_type == 2) THEN ! Linear 1403 1404 DO jk = 1, jpkm1-1 1405 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1406 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1407 1408 dsp(ji,jj,jk) = 0._wp 1409 csp(ji,jj,jk) = 0._wp 1410 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1411 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1412 END DO 1413 1414 ELSE 1415 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1416 ENDIF 1417 1418 END DO 1419 END DO 1364 1420 1365 Do ji = 1, jpi 1366 Do jj = 1, jpj 1367 1368 If (polynomial_type == 1) Then ! Constrained Cubic Spline 1369 Do jk = 2, jpkm1-1 1370 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1371 zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1372 zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1373 zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1374 1375 zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1376 1377 If(zdf1 * zdf2 <= 0._wp) Then 1378 zdf(jk) = 0._wp 1379 Else 1380 zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1381 End If 1382 End Do 1383 1384 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - & 1385 & 0.5_wp * zdf(2) 1386 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1387 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1388 & 0.5_wp * zdf(jpkm1 - 1) 1389 1390 Do jk = 1, jpkm1 - 1 1391 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1392 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1393 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1394 zddf1 = -2._wp * ztmp1 + ztmp2 1395 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1396 zddf2 = 2._wp * ztmp1 - ztmp2 1397 1398 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1399 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1400 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1401 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1402 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1403 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1404 & xsp(ji,jj,jk)**2 ) 1405 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1406 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1407 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1408 End Do 1409 1410 Else If (polynomial_type == 2) Then ! Linear 1411 1412 Do jk = 1, jpkm1-1 1413 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1414 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1415 1416 dsp(ji,jj,jk) = 0._wp 1417 csp(ji,jj,jk) = 0._wp 1418 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1419 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1420 End Do 1421 1422 Else 1423 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1424 End If 1425 1426 End Do 1427 End Do 1428 1429 End Subroutine cspline 1421 END SUBROUTINE cspline 1430 1422 1431 1423 … … 1451 1443 IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 1452 1444 f = 0.5_wp * (fl + fr) 1453 1445 ELSE 1454 1446 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1455 END 1447 ENDIF 1456 1448 1457 1449 END FUNCTION interp1 … … 1480 1472 !! *** ROUTINE interp1 *** 1481 1473 !! 1482 !! ** Purpose : deriavtive of a cubic spline function 1474 !! ** Purpose : Calculate the first order of deriavtive of 1475 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1483 1476 !! 1484 !! ** Method : 1477 !! ** Method : f=dy/dx=b+2*c*x+3*d*x^2 1485 1478 !! 1486 1479 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
