Changeset 3081 for branches/2011/dev_NOC_2011_MERGE
- Timestamp:
- 2011-11-11T13:38:49+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
r3074 r3081 94 94 ENDIF 95 95 ! 96 SELECT CASE ( nhpg ) ! Hydr astatic pressure gradient computation96 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 97 97 CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate 98 98 CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation) … … 171 171 IF( ln_hpg_prj ) nhpg = 7 172 172 ! 173 ! ! Consi tency check173 ! ! Consistency check 174 174 ioptio = 0 175 175 IF( ln_hpg_zco ) ioptio = ioptio + 1 … … 1045 1045 !! 1046 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 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 variables for the correction term 1053 INTEGER :: jk1, jis, jid, jjs, jjd 1054 REAL(wp) :: zuijk, zvijk, pwes, pwed, pnss, pnsd, deps 1055 REAL(wp) :: rhdt1 1057 1056 REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 1058 1057 INTEGER :: bhitwe, bhitns … … 1073 1072 zcoef0 = - grav 1074 1073 znad = 0.0_wp 1075 IF( lk_vvl) znad = 1._wp1074 IF( lk_vvl ) znad = 1._wp 1076 1075 1077 1076 ! Save fsde3w and rhd … … 1079 1078 rhd_tmp(:,:,:) = rhd(:,:,:) 1080 1079 1081 ! Clean 3-D work arra ies1080 ! Clean 3-D work arrays 1082 1081 zhpi(:,:,:) = 0. 1083 1082 1084 !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here 1085 1083 ! how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here 1086 1084 ! Preparing vertical density profile for hybrid-sco coordinate 1087 1085 DO jj = 1, jpj 1088 1086 DO ji = 1, jpi 1089 1087 jk = mbathy(ji,jj) 1090 IF( jk <= 0) THEN; rhd(ji,jj,:) = 0._wp1091 1092 1093 DO jkk = jk+1, jpk1094 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 DO1088 IF( jk <= 0 ) THEN; rhd(ji,jj,:) = 0._wp 1089 ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 1090 ELSE IF(jk < jpkm1) THEN 1091 DO jkk = jk+1, jpk 1092 rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), & 1093 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1094 END DO 1097 1095 END IF 1098 1096 END DO … … 1104 1102 fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad 1105 1103 DO jk = 2, jpk 1106 fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)1104 fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk) 1107 1105 END DO 1108 1106 END DO … … 1118 1116 END DO 1119 1117 1120 ! !Construct the vertical density profile with the1121 ! !Constrained cubic spline interpolation1118 ! Construct the vertical density profile with the 1119 ! constrained cubic spline interpolation 1122 1120 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 1123 1121 … … 1126 1124 DO ji = 2, jpi 1127 1125 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)1126 bsp(ji,jj,1), csp(ji,jj,1), & 1127 dsp(ji,jj,1) ) * 0.5_wp * fsde3w(ji,jj,1) 1130 1128 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1131 1129 1132 !! assuming linear profile across the top half surface layer1130 ! assuming linear profile across the top half surface layer 1133 1131 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * rhdt1 1134 1132 END DO … … 1139 1137 DO jj = 2, jpj 1140 1138 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))1139 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1140 integ2(fsde3w(ji,jj,jk-1), fsde3w(ji,jj,jk),& 1141 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), & 1142 csp(ji,jj,jk-1), dsp(ji,jj,jk-1)) 1145 1143 END DO 1146 1144 END DO … … 1148 1146 1149 1147 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1150 1151 1148 DO jj = 2, jpjm1 1152 1149 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)1150 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 1151 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 1155 1152 END DO 1156 1153 END DO … … 1159 1156 DO jj = 2, jpjm1 1160 1157 DO ji = 2, jpim1 1161 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u w(ji,jj,jk)1162 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v w(ji,jj,jk)1158 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 1159 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 1163 1160 END DO 1164 1161 END DO 1165 1162 END DO 1166 1163 1167 ! Start pressure integration 1164 DO jk = 1, jpkm1 1165 DO jj = 2, jpjm1 1166 DO ji = 2, jpim1 1167 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 1168 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 1169 END DO 1170 END DO 1171 END DO 1168 1172 1169 1173 DO jk = 1, jpkm1 … … 1176 1180 1177 1181 !!!!! 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 IF( jk <= mbku(ji,jj) ) THEN 1183 IF( -fsde3w(ji+1,jj,mbku(ji,jj)) >= -fsde3w(ji,jj,mbku(ji,jj)) ) THEN 1184 jis = ji + 1; jid = ji 1185 ELSE 1186 jis = ji; jid = ji +1 1187 ENDIF 1188 1189 ! integrate the pressure on the shallow side 1190 jk1 = jk 1191 bhitwe = 0 1192 DO WHILE ( -fsde3w(jis,jj,jk1) > zuijk ) 1193 IF( jk1 == mbku(ji,jj) ) THEN 1194 bhitwe = 1 1195 EXIT 1196 ENDIF 1197 deps = min(fsde3w(jis,jj,jk1+1), -zuijk) 1198 pwes = pwes + & 1199 integ2(fsde3w(jis,jj,jk1), deps, & 1200 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1201 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1202 jk1 = jk1 + 1 1203 END DO 1204 1205 IF(bhitwe == 1) THEN 1206 zuijk = -fsde3w(jis,jj,jk1) 1207 ENDIF 1208 1209 ! integrate the pressure on the deep side 1210 jk1 = jk 1211 bhitwe = 0 1212 DO WHILE ( -fsde3w(jid,jj,jk1) < zuijk ) 1213 IF( jk1 == 1 ) THEN 1214 bhitwe = 1 1215 EXIT 1216 ENDIF 1217 deps = max(fsde3w(jid,jj,jk1-1), -zuijk) 1218 pwed = pwed + & 1219 integ2(deps, fsde3w(jid,jj,jk1), & 1220 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1221 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1222 jk1 = jk1 - 1 1223 END DO 1224 1225 IF( bhitwe == 1 ) THEN 1226 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), & 1228 bsp(jid,jj,1), csp(jid,jj,1), & 1229 dsp(jid,jj,1)) * deps 1230 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1231 pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps 1232 ENDIF 1233 1234 1235 dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1236 IF( lk_vvl ) THEN 1237 dpdx2 = zcoef0 / e1u(ji,jj) * & 1238 ( REAL(jis-jid, wp) * (pwes + pwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1239 ELSE 1240 dpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (pwes + pwed) 1241 ENDIF 1242 1243 ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * & 1244 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1182 1245 ENDIF 1183 1184 ! !integrate the pressure on the shallow side1185 jk1 = jk1186 bhitwe = 01187 IF( jk1 == mbathy(jis,jj) ) THEN1188 bhitwe = 11189 ELSE1190 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 + 11196 IF( jk1 == mbathy(jis,jj) ) THEN1197 bhitwe = 11198 EXIT1199 END IF1200 END DO1201 ENDIF1202 1203 IF( bhitwe /= 1 ) THEN1204 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 ELSE1209 zuijk = -fsde3w(jis,jj,jk1)1210 ENDIF1211 1212 ! !integrate the pressure on the deep side1213 jk1 = jk1214 bhitwe = 01215 IF( jk1 == 1 ) THEN1216 bhitwe = 11217 ELSE1218 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 - 11224 IF( jk1 == 1 ) THEN1225 bhitwe = 11226 EXIT1227 END IF1228 END DO1229 ENDIF1230 1231 IF( bhitwe /= 1 ) THEN1232 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 ELSE1237 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) ) * deps1241 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water1242 pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps1243 ENDIF1244 1245 IF( jid > jis ) THEN1246 pe = pwed; pw = pwes1247 ELSE1248 pe = pwes; pw = pwed1249 ENDIF1250 1251 dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))1252 IF( lk_vvl ) THEN1253 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj)))1254 ELSE1255 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe)1256 ENDIF1257 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 1246 1261 1247 !!!!! 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 1248 IF( jk <= mbkv(ji,jj) ) THEN 1249 IF( -fsde3w(ji,jj+1,mbkv(ji,jj)) >= -fsde3w(ji,jj,mbkv(ji,jj)) ) THEN 1250 jjs = jj + 1; jjd = jj 1251 ELSE 1252 jjs = jj ; jjd = jj + 1 1253 ENDIF 1254 1255 ! integrate the pressure on the shallow side 1256 jk1 = jk 1257 bhitns = 0 1258 DO WHILE ( -fsde3w(ji,jjs,jk1) > zvijk ) 1259 IF( jk1 == mbkv(ji,jj) ) THEN 1260 bhitns = 1 1261 EXIT 1262 ENDIF 1263 deps = min(fsde3w(ji,jjs,jk1+1), -zvijk) 1264 pnss = pnss + & 1265 integ2(fsde3w(ji,jjs,jk1), deps, & 1266 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1267 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1268 jk1 = jk1 + 1 1269 END DO 1270 1271 IF(bhitns == 1) THEN 1272 zvijk = -fsde3w(ji,jjs,jk1) 1273 ENDIF 1274 1275 ! integrate the pressure on the deep side 1276 jk1 = jk 1277 bhitns = 0 1278 DO WHILE ( -fsde3w(ji,jjd,jk1) < zvijk ) 1279 IF( jk1 == 1 ) THEN 1280 bhitns = 1 1281 EXIT 1282 ENDIF 1283 deps = max(fsde3w(ji,jjd,jk1-1), -zvijk) 1284 pnsd = pnsd + & 1285 integ2(deps, fsde3w(ji,jjd,jk1), & 1286 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1287 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1288 jk1 = jk1 - 1 1289 END DO 1290 1291 IF( bhitns == 1 ) THEN 1292 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), & 1294 bsp(ji,jjd,1), csp(ji,jjd,1), & 1295 dsp(ji,jjd,1) ) * deps 1296 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water 1297 pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps 1298 ENDIF 1299 1300 1301 dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1302 IF( lk_vvl ) THEN 1303 dpdy2 = zcoef0 / e2v(ji,jj) * & 1304 ( REAL(jjs-jjd, wp) * (pnss + pnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1305 ELSE 1306 dpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (pnss + pnsd ) 1307 ENDIF 1308 1309 va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*& 1310 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1267 1311 ENDIF 1268 1269 ! !integrate the pressure on the shallow side1270 jk1 = jk1271 bhitns = 01272 IF( jk1 == mbathy(ji,jjs) ) THEN1273 bhitns = 11274 ELSE1275 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 + 11281 IF( jk1 == mbathy(ji,jjs) ) THEN1282 bhitns = 11283 EXIT1284 END IF1285 END DO1286 ENDIF1287 1288 IF( bhitns /= 1 ) THEN1289 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 ELSE1294 zvijk = -fsde3w(ji,jjs,jk1)1295 ENDIF1296 1297 ! !integrate the pressure on the deep side1298 jk1 = jk1299 bhitns = 01300 IF( jk1 == 1 ) THEN1301 bhitns = 11302 ELSE1303 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 - 11309 IF( jk1 == 1 ) THEN1310 bhitns = 11311 EXIT1312 END IF1313 END DO1314 ENDIF1315 1316 IF( bhitns /= 1 ) THEN1317 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 ELSE1322 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) ) * deps1326 rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water1327 pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps1328 ENDIF1329 1330 IF( jjd > jjs ) THEN1331 pn = pnsd; ps = pnss1332 ELSE1333 pn = pnss; ps = pnsd1334 ENDIF1335 1336 dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))1337 if( lk_vvl ) then1338 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj)))1339 else1340 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn )1341 end if1342 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 1312 1346 1313 … … 1399 1366 Do jj = 1, jpj 1400 1367 1401 If (polynomial_type == 1) Then ! Constrained Cubic Spline1368 If (polynomial_type == 1) Then ! Constrained Cubic Spline 1402 1369 Do jk = 2, jpkm1-1 1403 1370 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) … … 1441 1408 End Do 1442 1409 1443 Else If (polynomial_type == 2) Then ! Linear1410 Else If (polynomial_type == 2) Then ! Linear 1444 1411 1445 1412 Do jk = 1, jpkm1-1
Note: See TracChangeset
for help on using the changeset viewer.
