- Timestamp:
- 2011-09-28T10:25:27+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r2715 r2873 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_imp134 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj, rn_gamma , ln_dynhpg_imp 132 135 !!---------------------------------------------------------------------- 133 136 ! … … 147 150 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 151 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 152 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 149 153 WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma 150 154 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 151 155 ENDIF 152 156 ! 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') 157 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) ) & 158 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require:& 159 & the standard jacobian formulation hpg_sco or & 160 & the pressure jacobian formulation hpg_prj') 155 161 ! 156 162 ! ! Set nhpg from ln_hpg_... flags … … 162 168 IF( ln_hpg_djc ) nhpg = 5 163 169 IF( ln_hpg_rot ) nhpg = 6 170 IF( ln_hpg_prj ) nhpg = 7 164 171 ! 165 172 ! ! Consitency check … … 172 179 IF( ln_hpg_djc ) ioptio = ioptio + 1 173 180 IF( ln_hpg_rot ) ioptio = ioptio + 1 181 IF( ln_hpg_prj ) ioptio = ioptio + 1 174 182 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 175 183 ! … … 1005 1013 END SUBROUTINE hpg_rot 1006 1014 1015 1016 SUBROUTINE hpg_prj( kt ) 1017 !!--------------------------------------------------------------------- 1018 !! *** ROUTINE hpg_prj *** 1019 !! 1020 !! ** Method : s-coordinate case. 1021 !! Reformulate the horizontal hydrostatical pressure gradient 1022 !! term using Pressure Jacobian. A new correction term 1023 !! is developed to eliminate the sigma-coordinate error. 1024 !! 1025 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 1026 !! - Save the trend (l_trddyn=T) 1027 !! 1028 !!---------------------------------------------------------------------- 1029 1030 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1031 USE oce , ONLY: fsde3w_tmp => ta ! (ta,sa) used as 3D workspace 1032 USE oce , ONLY: rhd_tmp => sa 1033 USE wrk_nemo, ONLY: zhpi => wrk_3d_3 1034 USE wrk_nemo, ONLY: zu => wrk_3d_4 1035 USE wrk_nemo, ONLY: zv => wrk_3d_5 1036 USE wrk_nemo, ONLY: fsp => wrk_3d_6 1037 USE wrk_nemo, ONLY: xsp => wrk_3d_7 1038 USE wrk_nemo, ONLY: asp => wrk_3d_8 1039 USE wrk_nemo, ONLY: bsp => wrk_3d_9 1040 USE wrk_nemo, ONLY: csp => wrk_3d_10 1041 USE wrk_nemo, ONLY: dsp => wrk_3d_11 1042 !! 1043 !!---------------------------------------------------------------------- 1044 !! 1045 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 1046 INTEGER, INTENT(in) :: kt ! ocean time-step index 1047 !! 1048 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 1049 REAL(wp) :: zcoef0, znad ! temporary scalars 1050 !! 1051 !! The local varialbes for the correction term 1052 INTEGER :: jke, jkw, jkn, jks, jk1 1053 REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs 1054 REAL(wp) :: zuijk, zvijk, pe, pw, pn, ps, dze, dzw, dzn, dzs 1055 REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2 1056 REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 1057 INTEGER :: bhitwe, bhitns 1058 !!---------------------------------------------------------------------- 1059 1060 IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 1061 CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable') ; RETURN 1062 ENDIF 1063 1064 IF( kt == nit000 ) THEN 1065 IF(lwp) WRITE(numout,*) 1066 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 1067 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 1068 ENDIF 1069 1070 !!---------------------------------------------------------------------- 1071 ! Local constant initialization 1072 zcoef0 = - grav 1073 znad = 0.0_wp 1074 IF(lk_vvl) znad = 1._wp 1075 1076 ! Save fsde3w and rhd 1077 fsde3w_tmp(:,:,:) = fsde3w(:,:,:) 1078 rhd_tmp(:,:,:) = rhd(:,:,:) 1079 1080 ! Clean 3-D work arraies 1081 zhpi(:,:,:) = 0. 1082 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 ! Constrained cubic spline interpolation 1121 1122 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 1123 1124 1125 1126 ! Calculate the pressure at T(ji,jj,1) 1127 DO jj = 2, jpj 1128 DO ji = 2, jpi 1129 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1130 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1131 * 0.5_wp * fsde3w(ji,jj,1) 1132 rhdt1 = max(rhdt1, 0.0_wp) 1133 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * rhdt1 1134 END DO 1135 END DO 1136 1137 1138 1139 ! print*, 'max&min fse3w=',maxval(fse3w), minval(fse3w) 1140 ! print*, 'max&min rhd===',maxval(rhd), minval(rhd) 1141 1142 ! Calculate the pressure at T(ji,jj,2:jpkm1) 1143 DO jk = 2, jpkm1 1144 DO jj = 2, jpj 1145 DO ji = 2, jpi 1146 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1147 integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),& 1148 asp(ji,jj,jk-1),bsp(ji,jj,jk-1),& 1149 csp(ji,jj,jk-1),dsp(ji,jj,jk-1)) 1150 END DO 1151 END DO 1152 END DO 1153 1154 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1155 1156 DO jj = 2, jpjm1 1157 DO ji = 2, jpim1 1158 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 1159 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 1160 END DO 1161 END DO 1162 1163 DO jk = 2, jpkm1 1164 DO jj = 2, jpjm1 1165 DO ji = 2, jpim1 1166 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 1167 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 1168 END DO 1169 END DO 1170 END DO 1171 1172 DO jk = 1, jpkm1 1173 DO jj = 2, jpjm1 1174 DO ji = 2, jpim1 1175 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 1176 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 1177 END DO 1178 END DO 1179 END DO 1180 1181 DO jk = 1, jpkm1 1182 DO jj = 2, jpjm1 1183 DO ji = 2, jpim1 1184 1185 pe = 0.; pw = 0.; dze = 0.; dzw = 0. 1186 pn = 0.; ps = 0.; dzn = 0.; dzs = 0. 1187 bhitwe = 0; bhitns = 0 1188 zuijk = zu(ji,jj,jk) 1189 zvijk = zv(ji,jj,jk) 1190 1191 !! for u equation 1192 IF(-fsde3w(ji+1,jj,jk) > zuijk) THEN 1193 DO jk1 = jk+1, mbathy(ji+1,jj) 1194 IF(-fsde3w(ji+1,jj,jk1) > zuijk) THEN 1195 pe = pe + & 1196 integ2(fsde3w(ji+1,jj,jk1-1),fsde3w(ji+1,jj,jk1),& 1197 asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),& 1198 csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1)) 1199 IF(jk1 == mbathy(ji+1,jj)) bhitwe = 1 1200 ELSE 1201 zze = -fsde3w(ji+1,jj,jk1-1) - zuijk 1202 pe = pe + & 1203 integ2(fsde3w(ji+1,jj,jk1-1),-zuijk,& 1204 asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),& 1205 csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1)) 1206 EXIT 1207 ENDIF 1208 END DO 1209 1210 IF(jk == mbathy(ji+1,jj)) bhitwe = 1 1211 1212 IF(bhitwe == 1) zuijk = -fsde3w(ji+1,jj,mbathy(ji+1,jj)) 1213 1214 DO jk1 = jk-1, 1, -1 1215 IF(-fsde3w(ji,jj,jk1) < zuijk) THEN 1216 pw = pw + & 1217 integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),& 1218 asp(ji,jj,jk1),bsp(ji,jj,jk1),& 1219 csp(ji,jj,jk1),dsp(ji,jj,jk1)) 1220 IF(jk1 == 1) THEN 1221 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1222 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1223 * (zuijk + fsde3w(ji,jj,1)) 1224 rhdt1 = max(rhdt1,0.0_wp) 1225 rrw = rhdt1 1226 zzw = zuijk + fsde3w(ji,jj,1) 1227 pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1)) 1228 END IF 1229 ELSE 1230 zzw = zuijk + fsde3w(ji,jj,jk1+1) 1231 pw = pw + & 1232 integ2(-zuijk,fsde3w(ji,jj,jk1+1),& 1233 asp(ji,jj,jk1),bsp(ji,jj,jk1),& 1234 csp(ji,jj,jk1),dsp(ji,jj,jk1)) 1235 EXIT 1236 ENDIF 1237 END DO 1238 1239 IF(jk == 1) THEN 1240 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1241 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1242 * (zuijk + fsde3w(ji,jj,1)) 1243 rhdt1 = max(rhdt1,0.0_wp) 1244 rrw = rhdt1 1245 zzw = zuijk + fsde3w(ji,jj,1) 1246 pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1)) 1247 END IF 1248 1249 ELSE 1250 1251 DO jk1 = jk+1, mbathy(ji,jj) 1252 IF(-fsde3w(ji,jj,jk1) > zuijk) THEN 1253 pw = pw - & 1254 integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),& 1255 asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 1256 csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 1257 IF(jk1 == mbathy(ji,jj)) bhitwe = 1 1258 ELSE 1259 zzw = -fsde3w(ji,jj,jk1-1) - zuijk 1260 pw = pw - & 1261 integ2(fsde3w(ji,jj,jk1-1),-zuijk,& 1262 asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 1263 csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 1264 EXIT 1265 ENDIF 1266 END DO 1267 1268 IF(jk == mbathy(ji,jj)) bhitwe = 1 1269 1270 IF(bhitwe == 1) zuijk = -fsde3w(ji,jj,mbathy(ji,jj)) 1271 1272 DO jk1 = jk-1, 1, -1 1273 IF(-fsde3w(ji+1,jj,jk1) < zuijk) THEN 1274 pe = pe - & 1275 integ2(fsde3w(ji+1,jj,jk1),fsde3w(ji+1,jj,jk1+1),& 1276 asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),& 1277 csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1)) 1278 IF(jk1 == 1) THEN 1279 rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), & 1280 bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) & 1281 * (zuijk + fsde3w(ji+1,jj,1)) 1282 rhdt1 = max(rhdt1,0.0_wp) 1283 rre = rhdt1 1284 zze = zuijk + fsde3w(ji+1,jj,1) 1285 pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1)) 1286 END IF 1287 ELSE 1288 zze = zuijk + fsde3w(ji+1,jj,jk1+1) 1289 pe = pe - & 1290 integ2(-zuijk,fsde3w(ji+1,jj,jk1+1),& 1291 asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),& 1292 csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1)) 1293 EXIT 1294 ENDIF 1295 END DO 1296 1297 IF(jk == 1) THEN 1298 rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), & 1299 bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) & 1300 * (zuijk + fsde3w(ji+1,jj,1)) 1301 rhdt1 = max(rhdt1,0.0_wp) 1302 rre = rhdt1 1303 zze = zuijk + fsde3w(ji+1,jj,1) 1304 pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1)) 1305 END IF 1306 1307 ENDIF 1308 1309 1310 dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1311 IF(lk_vvl) THEN 1312 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj))) 1313 ELSE 1314 dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe) 1315 ENDIF 1316 1317 ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2)*& 1318 & umask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji+1,jj,jk) 1319 1320 !! for v equation 1321 1322 IF(-fsde3w(ji,jj+1,jk) > zvijk) THEN 1323 DO jk1 = jk+1, mbathy(ji,jj+1) 1324 IF(-fsde3w(ji,jj+1,jk1) > zvijk) THEN 1325 pn = pn + & 1326 integ2(fsde3w(ji,jj+1,jk1-1),fsde3w(ji,jj+1,jk1),& 1327 asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),& 1328 csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1)) 1329 IF(jk1 == mbathy(ji,jj+1)) bhitns = 1 1330 ELSE 1331 zzn = -fsde3w(ji,jj+1,jk1-1) - zvijk 1332 pn = pn + & 1333 integ2(fsde3w(ji,jj+1,jk1-1),-zvijk,& 1334 asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),& 1335 csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1)) 1336 EXIT 1337 ENDIF 1338 END DO 1339 1340 IF(jk == mbathy(ji,jj+1)) bhitns = 1 1341 1342 1343 IF(bhitns == 1) zvijk = -fsde3w(ji,jj+1,mbathy(ji,jj+1)) 1344 1345 1346 DO jk1 = jk-1, 1, -1 1347 IF(-fsde3w(ji,jj,jk1) < zvijk) THEN 1348 ps = ps + & 1349 integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),& 1350 asp(ji,jj,jk1),bsp(ji,jj,jk1),& 1351 csp(ji,jj,jk1),dsp(ji,jj,jk1)) 1352 IF(jk1 == 1) THEN 1353 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1354 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1355 * (zvijk + fsde3w(ji,jj,1)) 1356 rhdt1 = max(rhdt1,0.0_wp) 1357 rrs = rhdt1 1358 zzs = zvijk + fsde3w(ji,jj,1) 1359 ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs ) * zzs, zhpi(ji,jj,1)) 1360 END IF 1361 ELSE 1362 zzs = zvijk + fsde3w(ji,jj,jk1+1) 1363 ps = ps + & 1364 integ2(-zvijk,fsde3w(ji,jj,jk1+1),& 1365 asp(ji,jj,jk1),bsp(ji,jj,jk1),& 1366 csp(ji,jj,jk1),dsp(ji,jj,jk1)) 1367 EXIT 1368 ENDIF 1369 END DO 1370 1371 IF(jk == 1) THEN 1372 rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 1373 bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 1374 * (zvijk + fsde3w(ji,jj,1)) 1375 rhdt1 = max(rhdt1,0.0_wp) 1376 rrs = rhdt1 1377 zzs = zvijk + fsde3w(ji,jj,1) 1378 ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs) * zzs, zhpi(ji,jj,1)) 1379 END IF 1380 1381 ELSE 1382 1383 DO jk1 = jk+1, mbathy(ji,jj) 1384 IF(-fsde3w(ji,jj,jk1) > zvijk) THEN 1385 ps = ps - & 1386 integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),& 1387 asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 1388 csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 1389 IF(jk1 == mbathy(ji,jj)) bhitns = 1 1390 ELSE 1391 zzs = -fsde3w(ji,jj,jk1-1) - zvijk 1392 ps = ps - & 1393 integ2(fsde3w(ji,jj,jk1-1),-zvijk,& 1394 asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 1395 csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 1396 EXIT 1397 ENDIF 1398 END DO 1399 1400 IF(jk == mbathy(ji,jj)) bhitns = 1 1401 1402 1403 IF(bhitns == 1) zvijk = -fsde3w(ji,jj,mbathy(ji,jj)) 1404 1405 DO jk1 = jk-1, 1, -1 1406 IF(-fsde3w(ji,jj+1,jk1) < zvijk) THEN 1407 pn = pn - & 1408 integ2(fsde3w(ji,jj+1,jk1),fsde3w(ji,jj+1,jk1+1),& 1409 asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),& 1410 csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1)) 1411 IF(jk1 == 1) THEN 1412 rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), & 1413 bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) & 1414 * (zvijk + fsde3w(ji,jj+1,1)) 1415 rhdt1 = max(rhdt1,0.0_wp) 1416 rrn = rhdt1 1417 zzn = zvijk + fsde3w(ji,jj+1,1) 1418 pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1)) 1419 END IF 1420 1421 ELSE 1422 zzn = zvijk + fsde3w(ji,jj+1,jk1+1) 1423 pn = pn - & 1424 integ2(-zvijk,fsde3w(ji,jj+1,jk1+1),& 1425 asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),& 1426 csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1)) 1427 EXIT 1428 ENDIF 1429 END DO 1430 1431 IF(jk == 1) THEN 1432 rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), & 1433 bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) & 1434 * (zvijk + fsde3w(ji,jj+1,1)) 1435 rhdt1 = max(rhdt1,0.0_wp) 1436 rrn = rhdt1 1437 zzn = zvijk + fsde3w(ji,jj+1,1) 1438 pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1)) 1439 END IF 1440 1441 ENDIF 1442 1443 1444 dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1445 if(lk_vvl) then 1446 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj))) 1447 else 1448 dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn ) 1449 end if 1450 1451 va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*& 1452 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1453 1454 1455 END DO 1456 END DO 1457 END DO 1458 1459 ! Restore fsde3w and rhd 1460 1461 fsde3w(:,:,:) = fsde3w_tmp(:,:,:) 1462 rhd(:,:,:) = rhd_tmp(:,:,:) 1463 1464 ! 1465 IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) ) & 1466 CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 1467 ! 1468 1469 END SUBROUTINE hpg_prj 1470 1471 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 1472 !!---------------------------------------------------------------------- 1473 !! *** ROUTINE cspline *** 1474 !! 1475 !! ** Purpose : constrained cubic spline interpolation 1476 !! 1477 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1478 !! Reference: K.W. Brodlie, A review of mehtods for curve and function 1479 !! drawing, 1980 1480 !! 1481 !!---------------------------------------------------------------------- 1482 IMPLICIT NONE 1483 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1484 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1485 ! the interpoated function 1486 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1487 ! 2: Linear 1488 1489 ! Local Variables 1490 INTEGER :: ji, jj, jk ! dummy loop indices 1491 INTEGER :: jpi, jpj, jpkm1 1492 REAL(wp) :: df1, df2, ddf1, ddf2, tmp1, tmp2, dxtmp 1493 REAL(wp) :: dxtmp1, dxtmp2, alpha 1494 REAL(wp) :: df(size(fsp,3)) 1495 !!---------------------------------------------------------------------- 1496 1497 jpi = size(fsp,1) 1498 jpj = size(fsp,2) 1499 jpkm1 = size(fsp,3) - 1 1500 1501 ! Clean output arrays 1502 asp = 0.0_wp 1503 bsp = 0.0_wp 1504 csp = 0.0_wp 1505 dsp = 0.0_wp 1506 1507 1508 Do ji = 1, jpi 1509 Do jj = 1, jpj 1510 1511 If (polynomial_type == 1) Then !Constrained Cubic Spline 1512 Do jk = 2, jpkm1-1 1513 dxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1514 dxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1515 df1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / dxtmp1 1516 df2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / dxtmp2 1517 1518 alpha = ( dxtmp1 + 2._wp * dxtmp2 ) / ( dxtmp1 + dxtmp2 ) / 3._wp 1519 1520 If(df1 * df2 <= 0._wp) Then 1521 df(jk) = 0._wp 1522 Else 1523 df(jk) = df1 * df2 / ( ( 1._wp - alpha ) * df1 + alpha * df2 ) 1524 End If 1525 End Do 1526 1527 df(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - & 1528 & 0.5_wp * df(2) 1529 df(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1530 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1531 & 0.5_wp * df(jpkm1 - 1) 1532 1533 Do jk = 1, jpkm1 - 1 1534 dxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1535 tmp1 = (df(jk+1) + 2._wp * df(jk)) / dxtmp 1536 tmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / dxtmp / dxtmp 1537 ddf1 = -2._wp * tmp1 + tmp2 1538 tmp1 = (2._wp * df(jk+1) + df(jk)) / dxtmp 1539 ddf2 = 2._wp * tmp1 - tmp2 1540 1541 dsp(ji,jj,jk) = (ddf2 - ddf1) / 6._wp / dxtmp 1542 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * ddf1 - xsp(ji,jj,jk)*ddf2 ) / 2._wp / dxtmp 1543 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / dxtmp - & 1544 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1545 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1546 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1547 & xsp(ji,jj,jk)**2 ) 1548 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1549 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1550 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1551 End Do 1552 1553 Else If (polynomial_type == 2) Then !Linear 1554 1555 Do jk = 1, jpkm1-1 1556 dxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1557 tmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1558 1559 dsp(ji,jj,jk) = 0._wp 1560 csp(ji,jj,jk) = 0._wp 1561 bsp(ji,jj,jk) = tmp1 / dxtmp 1562 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1563 End Do 1564 1565 Else 1566 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1567 End If 1568 1569 End Do 1570 End Do 1571 1572 End Subroutine cspline 1573 1574 1575 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1576 !!---------------------------------------------------------------------- 1577 !! *** ROUTINE interp1 *** 1578 !! 1579 !! ** Purpose : 1-d linear interpolation 1580 !! 1581 !! ** Method : 1582 !! interpolation is straigt forward 1583 !! extrapolation is also permitted (no value limit) 1584 !! 1585 !! H.Liu, Jan 2009, POL 1586 !!---------------------------------------------------------------------- 1587 IMPLICIT NONE 1588 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1589 REAL(wp) :: f ! result of the interpolation (extrapolation) 1590 REAL(wp) :: deltx 1591 !!---------------------------------------------------------------------- 1592 1593 deltx = xr - xl 1594 IF(abs(deltx) <= 10._wp * EPSILON(x)) THEN 1595 f = 0.5_wp * (fl + fr) 1596 ELSE 1597 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / deltx 1598 END IF 1599 1600 END FUNCTION interp1 1601 1602 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1603 !!---------------------------------------------------------------------- 1604 !! *** ROUTINE interp1 *** 1605 !! 1606 !! ** Purpose : 1-d constrained cubic spline interpolation 1607 !! 1608 !! ** Method : cubic spline interpolation 1609 !! 1610 !!---------------------------------------------------------------------- 1611 IMPLICIT NONE 1612 REAL(wp), INTENT(in) :: x, a, b, c, d 1613 REAL(wp) :: f ! value from the interpolation 1614 !!---------------------------------------------------------------------- 1615 1616 f = a + x* ( b + x * ( c + d * x ) ) 1617 1618 END FUNCTION interp2 1619 1620 1621 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1622 !!---------------------------------------------------------------------- 1623 !! *** ROUTINE interp1 *** 1624 !! 1625 !! ** Purpose : deriavtive of a cubic spline function 1626 !! 1627 !! ** Method : 1628 !! 1629 !!---------------------------------------------------------------------- 1630 IMPLICIT NONE 1631 REAL(wp), INTENT(in) :: x, a, b, c, d 1632 REAL(wp) :: f ! value from the interpolation 1633 !!---------------------------------------------------------------------- 1634 1635 f = b + x * ( 2._wp * c + 3._wp * d * x) 1636 1637 END FUNCTION interp3 1638 1639 1640 FUNCTION integ2(xl, xr, a, b, c, d) RESULT(f) 1641 !!---------------------------------------------------------------------- 1642 !! *** ROUTINE interp1 *** 1643 !! 1644 !! ** Purpose : 1-d constrained cubic spline integration 1645 !! 1646 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1647 !! 1648 !!---------------------------------------------------------------------- 1649 IMPLICIT NONE 1650 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1651 REAL(wp) :: a1, a2,a3 1652 REAL(wp) :: f ! integration result 1653 !!---------------------------------------------------------------------- 1654 1655 a1 = 0.5_wp * b 1656 a2 = c / 3.0_wp 1657 a3 = 0.25_wp * d 1658 1659 f = xr * ( a + xr * ( a1 + xr * ( a2 + a3 * xr ) ) ) - & 1660 & xl * ( a + xl * ( a1 + xl * ( a2 + a3 * xl ) ) ) 1661 1662 END FUNCTION integ2 1663 1664 1007 1665 !!====================================================================== 1008 1666 END MODULE dynhpg
Note: See TracChangeset
for help on using the changeset viewer.
