- Timestamp:
- 2011-11-10T13:43:36+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
r2715 r3074 27 27 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 28 !! hpg_rot : s-coordinate (ROTated axes scheme) 29 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 29 30 !!---------------------------------------------------------------------- 30 31 USE oce ! ocean dynamics and tracers … … 52 53 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 53 54 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme) 55 LOGICAL , PUBLIC :: ln_hpg_prj = .FALSE. !: s-coordinate (Pressure Jacobian scheme) 54 56 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient 55 57 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 56 58 57 INTEGER :: nhpg = 0 ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)59 INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 58 60 59 61 !! * Substitutions … … 100 102 CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 101 103 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 104 CASE ( 7 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 105 END SELECT 103 106 ! … … 129 132 !! 130 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, & 131 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma , ln_dynhpg_imp 134 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj, & 135 & rn_gamma , ln_dynhpg_imp 132 136 !!---------------------------------------------------------------------- 133 137 ! … … 147 151 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 152 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 153 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 149 154 WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma 150 155 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 151 156 ENDIF 152 157 ! 153 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) & 154 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 158 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) ) & 159 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 160 & the standard jacobian formulation hpg_sco or & 161 & the pressure jacobian formulation hpg_prj') 155 162 ! 156 163 ! ! Set nhpg from ln_hpg_... flags … … 162 169 IF( ln_hpg_djc ) nhpg = 5 163 170 IF( ln_hpg_rot ) nhpg = 6 171 IF( ln_hpg_prj ) nhpg = 7 164 172 ! 165 173 ! ! Consitency check … … 172 180 IF( ln_hpg_djc ) ioptio = ioptio + 1 173 181 IF( ln_hpg_rot ) ioptio = ioptio + 1 182 IF( ln_hpg_prj ) ioptio = ioptio + 1 174 183 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 175 184 ! … … 1005 1014 END SUBROUTINE hpg_rot 1006 1015 1016 1017 SUBROUTINE hpg_prj( kt ) 1018 !!--------------------------------------------------------------------- 1019 !! *** ROUTINE hpg_prj *** 1020 !! 1021 !! ** Method : s-coordinate case. 1022 !! Reformulate the horizontal hydrostatical pressure gradient 1023 !! term using Pressure Jacobian. A new correction term 1024 !! is developed to eliminate the sigma-coordinate error. 1025 !! 1026 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 1027 !! - Save the trend (l_trddyn=T) 1028 !! 1029 !!---------------------------------------------------------------------- 1030 1031 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1032 USE oce , ONLY: fsde3w_tmp => ta ! (ta,sa) used as 3D workspace 1033 USE oce , ONLY: rhd_tmp => sa 1034 USE wrk_nemo, ONLY: zhpi => wrk_3d_3 1035 USE wrk_nemo, ONLY: zu => wrk_3d_4 1036 USE wrk_nemo, ONLY: zv => wrk_3d_5 1037 USE wrk_nemo, ONLY: fsp => wrk_3d_6 1038 USE wrk_nemo, ONLY: xsp => wrk_3d_7 1039 USE wrk_nemo, ONLY: asp => wrk_3d_8 1040 USE wrk_nemo, ONLY: bsp => wrk_3d_9 1041 USE wrk_nemo, ONLY: csp => wrk_3d_10 1042 USE wrk_nemo, ONLY: dsp => wrk_3d_11 1043 !! 1044 !!---------------------------------------------------------------------- 1045 !! 1046 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 1047 INTEGER, INTENT(in) :: kt ! ocean time-step index 1048 !! 1049 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 1050 REAL(wp) :: zcoef0, znad ! temporary scalars 1051 !! 1052 !! The local varialbes for the correction term 1053 INTEGER :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd 1054 REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs 1055 REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps 1056 REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2 1057 REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 1058 INTEGER :: bhitwe, bhitns 1059 !!---------------------------------------------------------------------- 1060 1061 IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 1062 CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable') ; RETURN 1063 ENDIF 1064 1065 IF( kt == nit000 ) THEN 1066 IF(lwp) WRITE(numout,*) 1067 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 1068 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 1069 ENDIF 1070 1071 !!---------------------------------------------------------------------- 1072 ! Local constant initialization 1073 zcoef0 = - grav 1074 znad = 0.0_wp 1075 IF(lk_vvl) znad = 1._wp 1076 1077 ! Save fsde3w and rhd 1078 fsde3w_tmp(:,:,:) = fsde3w(:,:,:) 1079 rhd_tmp(:,:,:) = rhd(:,:,:) 1080 1081 ! Clean 3-D work arraies 1082 zhpi(:,:,:) = 0. 1083 1084 !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here 1085 1086 ! Preparing vertical density profile for hybrid-sco coordinate 1087 DO jj = 1, jpj 1088 DO ji = 1, jpi 1089 jk = mbathy(ji,jj) 1090 IF(jk <= 0) THEN; rhd(ji,jj,:) = 0._wp 1091 ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 1092 ELSE IF(jk < jpkm1) THEN 1093 DO jkk = jk+1, jpk 1094 rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),& 1095 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1096 END DO 1097 END IF 1098 END DO 1099 END DO 1100 1101 DO jj = 1, jpj 1102 DO ji = 1, jpi 1103 fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 1104 fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad 1105 DO jk = 2, jpk 1106 fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk) 1107 END DO 1108 END DO 1109 END DO 1110 1111 DO jk = 1, jpkm1 1112 DO jj = 1, jpj 1113 DO ji = 1, jpi 1114 fsp(ji,jj,jk) = rhd(ji,jj,jk) 1115 xsp(ji,jj,jk) = fsde3w(ji,jj,jk) 1116 END DO 1117 END DO 1118 END DO 1119 1120 ! ! Construct the vertical density profile with the 1121 ! !Constrained cubic spline interpolation 1122 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 1123 1124 ! Calculate the hydrostatic pressure at T(ji,jj,1) 1125 DO jj = 2, jpj 1126 DO ji = 2, jpi 1127 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1128 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1129 * 0.5_wp * fsde3w(ji,jj,1) 1130 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1131 1132 ! ! assuming linear profile across the top half surface layer 1133 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * rhdt1 1134 END DO 1135 END DO 1136 1137 ! Calculate the pressure at T(ji,jj,2:jpkm1) 1138 DO jk = 2, jpkm1 1139 DO jj = 2, jpj 1140 DO ji = 2, jpi 1141 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1142 integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),& 1143 asp(ji,jj,jk-1),bsp(ji,jj,jk-1),& 1144 csp(ji,jj,jk-1),dsp(ji,jj,jk-1)) 1145 END DO 1146 END DO 1147 END DO 1148 1149 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1150 1151 DO jj = 2, jpjm1 1152 DO ji = 2, jpim1 1153 zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad) 1154 zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad) 1155 END DO 1156 END DO 1157 1158 DO jk = 2, jpkm1 1159 DO jj = 2, jpjm1 1160 DO ji = 2, jpim1 1161 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk) 1162 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk) 1163 END DO 1164 END DO 1165 END DO 1166 1167 ! Start pressure integration 1168 1169 DO jk = 1, jpkm1 1170 DO jj = 2, jpjm1 1171 DO ji = 2, jpim1 1172 pwes = 0._wp; pwed = 0._wp 1173 pnss = 0._wp; pnsd = 0._wp 1174 zuijk = zu(ji,jj,jk) 1175 zvijk = zv(ji,jj,jk) 1176 1177 !!!!! for u equation 1178 IF( -fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN 1179 jis = ji + 1; jid = ji 1180 ELSE 1181 jis = ji; jid = ji +1 1182 ENDIF 1183 1184 ! !integrate the pressure on the shallow side 1185 jk1 = jk 1186 bhitwe = 0 1187 IF( jk1 == mbathy(jis,jj) ) THEN 1188 bhitwe = 1 1189 ELSE 1190 DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk ) 1191 pwes = pwes + & 1192 integ2(fsde3w(jis,jj,jk1), fsde3w(jis,jj,jk1+1),& 1193 asp(jis,jj,jk1) , bsp(jis,jj,jk1), & 1194 csp(jis,jj,jk1) , dsp(jis,jj,jk1)) 1195 jk1 = jk1 + 1 1196 IF( jk1 == mbathy(jis,jj) ) THEN 1197 bhitwe = 1 1198 EXIT 1199 END IF 1200 END DO 1201 ENDIF 1202 1203 IF( bhitwe /= 1 ) THEN 1204 pwes = pwes + & 1205 integ2(fsde3w(jis,jj,jk1), -zuijk, & 1206 asp(jis,jj,jk1), bsp(jis,jj,jk1),& 1207 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1208 ELSE 1209 zuijk = -fsde3w(jis,jj,jk1) 1210 ENDIF 1211 1212 ! !integrate the pressure on the deep side 1213 jk1 = jk 1214 bhitwe = 0 1215 IF( jk1 == 1 ) THEN 1216 bhitwe = 1 1217 ELSE 1218 DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk ) 1219 pwed = pwed + & 1220 integ2(fsde3w(jid,jj,jk1-1), fsde3w(jid,jj,jk1),& 1221 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1222 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1)) 1223 jk1 = jk1 - 1 1224 IF( jk1 == 1 ) THEN 1225 bhitwe = 1 1226 EXIT 1227 END IF 1228 END DO 1229 ENDIF 1230 1231 IF( bhitwe /= 1 ) THEN 1232 pwed = pwed + & 1233 integ2(-zuijk, fsde3w(jid,jj,jk1),& 1234 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1235 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1)) 1236 ELSE 1237 deps = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad) 1238 rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), & 1239 bsp(jid,jj,1), csp(jid,jj,1), & 1240 dsp(jid,jj,1) ) * deps 1241 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1242 pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps 1243 ENDIF 1244 1245 IF( jid > jis ) THEN 1246 pe = pwed; pw = pwes 1247 ELSE 1248 pe = pwes; pw = pwed 1249 ENDIF 1250 1251 dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1252 IF( lk_vvl ) THEN 1253 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj))) 1254 ELSE 1255 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe) 1256 ENDIF 1257 1258 ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * & 1259 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1260 1261 !!!!! for v equation 1262 1263 IF( -fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN 1264 jjs = jj + 1; jjd = jj 1265 ELSE 1266 jjs = jj ; jjd = jj + 1 1267 ENDIF 1268 1269 ! !integrate the pressure on the shallow side 1270 jk1 = jk 1271 bhitns = 0 1272 IF( jk1 == mbathy(ji,jjs) ) THEN 1273 bhitns = 1 1274 ELSE 1275 DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk ) 1276 pnss = pnss + & 1277 integ2(fsde3w(ji,jjs,jk1), fsde3w(ji,jjs,jk1+1),& 1278 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1279 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1280 jk1 = jk1 + 1 1281 IF( jk1 == mbathy(ji,jjs) ) THEN 1282 bhitns = 1 1283 EXIT 1284 END IF 1285 END DO 1286 ENDIF 1287 1288 IF( bhitns /= 1 ) THEN 1289 pnss = pnss + & 1290 integ2(fsde3w(ji,jjs,jk1), -zvijk, & 1291 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1292 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1293 ELSE 1294 zvijk = -fsde3w(ji,jjs,jk1) 1295 ENDIF 1296 1297 ! !integrate the pressure on the deep side 1298 jk1 = jk 1299 bhitns = 0 1300 IF( jk1 == 1 ) THEN 1301 bhitns = 1 1302 ELSE 1303 DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk ) 1304 pnsd = pnsd + & 1305 integ2(fsde3w(ji,jjd,jk1-1), fsde3w(ji,jjd,jk1), & 1306 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1307 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1308 jk1 = jk1 - 1 1309 IF( jk1 == 1 ) THEN 1310 bhitns = 1 1311 EXIT 1312 END IF 1313 END DO 1314 ENDIF 1315 1316 IF( bhitns /= 1 ) THEN 1317 pnsd = pnsd + & 1318 integ2(-zvijk, fsde3w(ji,jjd,jk1), & 1319 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1320 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1321 ELSE 1322 deps = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad) 1323 rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), & 1324 bsp(ji,jjd,1), csp(ji,jjd,1), & 1325 dsp(ji,jjd,1) ) * deps 1326 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1327 pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps 1328 ENDIF 1329 1330 IF( jjd > jjs ) THEN 1331 pn = pnsd; ps = pnss 1332 ELSE 1333 pn = pnss; ps = pnsd 1334 ENDIF 1335 1336 dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1337 if( lk_vvl ) then 1338 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj))) 1339 else 1340 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn ) 1341 end if 1342 1343 va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)* & 1344 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1345 1346 1347 END DO 1348 END DO 1349 END DO 1350 1351 ! Restore fsde3w and rhd 1352 fsde3w(:,:,:) = fsde3w_tmp(:,:,:) 1353 rhd(:,:,:) = rhd_tmp(:,:,:) 1354 1355 ! 1356 IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) ) & 1357 CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 1358 ! 1359 END SUBROUTINE hpg_prj 1360 1361 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 1362 !!---------------------------------------------------------------------- 1363 !! *** ROUTINE cspline *** 1364 !! 1365 !! ** Purpose : constrained cubic spline interpolation 1366 !! 1367 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1368 !! Reference: K.W. Brodlie, A review of mehtods for curve and function 1369 !! drawing, 1980 1370 !! 1371 !!---------------------------------------------------------------------- 1372 IMPLICIT NONE 1373 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1374 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1375 ! the interpoated function 1376 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1377 ! 2: Linear 1378 1379 ! Local Variables 1380 INTEGER :: ji, jj, jk ! dummy loop indices 1381 INTEGER :: jpi, jpj, jpkm1 1382 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 1383 REAL(wp) :: zdxtmp1, zdxtmp2, zalpha 1384 REAL(wp) :: zdf(size(fsp,3)) 1385 !!---------------------------------------------------------------------- 1386 1387 jpi = size(fsp,1) 1388 jpj = size(fsp,2) 1389 jpkm1 = size(fsp,3) - 1 1390 1391 ! Clean output arrays 1392 asp = 0.0_wp 1393 bsp = 0.0_wp 1394 csp = 0.0_wp 1395 dsp = 0.0_wp 1396 1397 1398 Do ji = 1, jpi 1399 Do jj = 1, jpj 1400 1401 If (polynomial_type == 1) Then !Constrained Cubic Spline 1402 Do jk = 2, jpkm1-1 1403 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1404 zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1405 zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1406 zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1407 1408 zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1409 1410 If(zdf1 * zdf2 <= 0._wp) Then 1411 zdf(jk) = 0._wp 1412 Else 1413 zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1414 End If 1415 End Do 1416 1417 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - & 1418 & 0.5_wp * zdf(2) 1419 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1420 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1421 & 0.5_wp * zdf(jpkm1 - 1) 1422 1423 Do jk = 1, jpkm1 - 1 1424 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1425 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1426 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1427 zddf1 = -2._wp * ztmp1 + ztmp2 1428 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1429 zddf2 = 2._wp * ztmp1 - ztmp2 1430 1431 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1432 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1433 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1434 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1435 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1436 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1437 & xsp(ji,jj,jk)**2 ) 1438 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1439 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1440 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1441 End Do 1442 1443 Else If (polynomial_type == 2) Then !Linear 1444 1445 Do jk = 1, jpkm1-1 1446 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1447 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1448 1449 dsp(ji,jj,jk) = 0._wp 1450 csp(ji,jj,jk) = 0._wp 1451 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1452 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1453 End Do 1454 1455 Else 1456 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1457 End If 1458 1459 End Do 1460 End Do 1461 1462 End Subroutine cspline 1463 1464 1465 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1466 !!---------------------------------------------------------------------- 1467 !! *** ROUTINE interp1 *** 1468 !! 1469 !! ** Purpose : 1-d linear interpolation 1470 !! 1471 !! ** Method : 1472 !! interpolation is straight forward 1473 !! extrapolation is also permitted (no value limit) 1474 !! 1475 !! H.Liu, Jan 2009, POL 1476 !!---------------------------------------------------------------------- 1477 IMPLICIT NONE 1478 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1479 REAL(wp) :: f ! result of the interpolation (extrapolation) 1480 REAL(wp) :: zdeltx 1481 !!---------------------------------------------------------------------- 1482 1483 zdeltx = xr - xl 1484 IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 1485 f = 0.5_wp * (fl + fr) 1486 ELSE 1487 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1488 END IF 1489 1490 END FUNCTION interp1 1491 1492 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1493 !!---------------------------------------------------------------------- 1494 !! *** ROUTINE interp1 *** 1495 !! 1496 !! ** Purpose : 1-d constrained cubic spline interpolation 1497 !! 1498 !! ** Method : cubic spline interpolation 1499 !! 1500 !!---------------------------------------------------------------------- 1501 IMPLICIT NONE 1502 REAL(wp), INTENT(in) :: x, a, b, c, d 1503 REAL(wp) :: f ! value from the interpolation 1504 !!---------------------------------------------------------------------- 1505 1506 f = a + x* ( b + x * ( c + d * x ) ) 1507 1508 END FUNCTION interp2 1509 1510 1511 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1512 !!---------------------------------------------------------------------- 1513 !! *** ROUTINE interp1 *** 1514 !! 1515 !! ** Purpose : deriavtive of a cubic spline function 1516 !! 1517 !! ** Method : 1518 !! 1519 !!---------------------------------------------------------------------- 1520 IMPLICIT NONE 1521 REAL(wp), INTENT(in) :: x, a, b, c, d 1522 REAL(wp) :: f ! value from the interpolation 1523 !!---------------------------------------------------------------------- 1524 1525 f = b + x * ( 2._wp * c + 3._wp * d * x) 1526 1527 END FUNCTION interp3 1528 1529 1530 FUNCTION integ2(xl, xr, a, b, c, d) RESULT(f) 1531 !!---------------------------------------------------------------------- 1532 !! *** ROUTINE interp1 *** 1533 !! 1534 !! ** Purpose : 1-d constrained cubic spline integration 1535 !! 1536 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1537 !! 1538 !!---------------------------------------------------------------------- 1539 IMPLICIT NONE 1540 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1541 REAL(wp) :: za1, za2, za3 1542 REAL(wp) :: f ! integration result 1543 !!---------------------------------------------------------------------- 1544 1545 za1 = 0.5_wp * b 1546 za2 = c / 3.0_wp 1547 za3 = 0.25_wp * d 1548 1549 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1550 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1551 1552 END FUNCTION integ2 1553 1554 1007 1555 !!====================================================================== 1008 1556 END MODULE dynhpg
Note: See TracChangeset
for help on using the changeset viewer.
