Changeset 1348 for trunk/NEMO/OPA_SRC/DOM
- Timestamp:
- 2009-03-30T18:42:51+02:00 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r1273 r1348 112 112 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 113 113 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 114 ! 114 115 !!bug gm control print: 116 IF( nprint == 1 .AND. lwp ) THEN 117 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 118 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 119 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) ) 120 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), & 121 & ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), & 122 & ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), & 123 & ' w ', MINVAL( fse3w(:,:,:) ) 124 125 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 126 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) ) 127 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), & 128 & ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), & 129 & ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), & 130 & ' w ', MAXVAL( fse3w(:,:,:) ) 131 ENDIF 132 !!!bug gm 133 115 134 END SUBROUTINE dom_zgr 116 135 … … 264 283 INTEGER :: ji, jj, jl, jk ! dummy loop indices 265 284 INTEGER :: inum ! temporary logical unit 285 INTEGER :: ii_bump, ij_bump, ih ! bump center position 266 286 INTEGER :: ii0, ii1, ij0, ij1 ! indices 267 INTEGER :: ii_bump, ij_bump, ih ! bump center position268 287 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 269 288 REAL(wp) :: zi , zj , zh ! temporary scalars … … 394 413 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 395 414 CALL iom_close (inum) 396 ! ! =====================415 ! ! ===================== 397 416 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 398 417 ! ! ===================== … … 422 441 ! 423 442 ENDIF 424 425 443 ! 444 ENDIF 426 445 ! ! =============== ! 427 446 ELSE ! error ! … … 996 1015 997 1016 1017 FUNCTION fssig1( pk1, bb ) RESULT( pf1 ) 1018 !!---------------------------------------------------------------------- 1019 !! *** ROUTINE eos_init *** 1020 !! 1021 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1022 !! 1023 !! ** Method : the function provides the non-dimensional position of 1024 !! T and W (i.e. between 0 and 1) 1025 !! T-points at integer values (between 1 and jpk) 1026 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1027 !! 1028 !! Reference : ??? 1029 !!---------------------------------------------------------------------- 1030 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 1031 REAL(wp), INTENT(in ) :: bb ! Stretching coefficient 1032 REAL(wp) :: pf1 ! sigma value 1033 1034 !!---------------------------------------------------------------------- 1035 ! 1036 IF ( theta == 0 ) then !uniform sigma 1037 pf1 = -(pk1-0.5)/REAL(jpkm1) 1038 ELSE ! stretched sigma 1039 pf1 = (1.0-bb) * (sinh( theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(theta) + & 1040 bb * ( (tanh( theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*theta) ) / & 1041 (2*tanh(0.5*theta) ) ) 1042 ENDIF 1043 1044 END FUNCTION fssig1 1045 1046 998 1047 SUBROUTINE zgr_sco 999 1048 !!---------------------------------------------------------------------- … … 1033 1082 REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace 1034 1083 REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - - 1035 !! 1036 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 1084 1085 LOGICAL :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true. 1086 REAL(wp) :: bb = 0.8 ! stretching parameter for song and haidvogel stretching, bb=0; top only, bb =1; top and bottom 1087 REAL(wp) :: hc = 150 ! Critical depth for s-sigma coordinates 1088 1089 REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3, gsigt3, gsi3w3, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1090 !! 1091 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max, ln_s_sigma, bb, hc 1037 1092 !!---------------------------------------------------------------------- 1038 1093 … … 1046 1101 WRITE(numout,*) ' Namelist nam_zgr_sco' 1047 1102 WRITE(numout,*) ' sigma-stretching coeffs ' 1048 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max 1049 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min 1050 WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta 1051 WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb 1052 WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max 1103 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max 1104 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min 1105 WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta 1106 WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb 1107 WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max 1108 WRITE(numout,*) ' Critical depth hc = ', hc 1109 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1053 1110 ENDIF 1054 1111 … … 1057 1114 hifv(:,:) = sbot_min 1058 1115 hiff(:,:) = sbot_min 1059 ! ! set maximum ocean depth 1116 1117 ! ! set maximum ocean depth 1060 1118 bathy(:,:) = MIN( sbot_max, bathy(:,:) ) 1061 1119 … … 1063 1121 ! ! Define the envelop bathymetry (hbatt) 1064 1122 ! ! ============================= 1065 ! Smooth the bathymetry (if required) 1123 ! Smooth the bathymetry (if required) 1066 1124 scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea) 1067 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1125 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1068 1126 ! 1069 1127 ! use r-value to create hybrid coordinates … … 1099 1157 DO jj = 1, nlcj 1100 1158 DO ji = 1, nlci 1101 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 1102 END DO 1103 END DO 1104 ! 1159 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 1160 END DO 1161 END DO 1162 ! 1105 1163 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1106 1164 ! … … 1124 1182 END DO 1125 1183 END DO 1126 ! 1184 ! 1127 1185 DO jj = 1, nlcj 1128 1186 DO ji = 1, nlci … … 1237 1295 ! 1238 1296 ! non-dimensional "sigma" for model level depth at w- and t-levels 1239 DO jk = 1, jpk 1240 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1241 gsigt(jk) = -fssig( REAL(jk,wp) ) 1242 END DO 1243 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1244 ! 1245 ! Coefficients for vertical scale factors at w-, t- levels 1297 1298 IF ( ln_s_sigma ) THEN !Song and Haidvogel style stretched sigma for depths below hc, with uniform sigma in shallower waters 1299 1300 DO ji=1,jpi 1301 DO jj=1,jpj 1302 1303 IF (hbatt(ji,jj).GT.hc) THEN !deep water, stretched sigma 1304 DO jk = 1, jpk 1305 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, bb ) 1306 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , bb ) 1307 END DO 1308 ELSE ! shallow water, uniform sigma 1309 DO jk = 1, jpk 1310 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1311 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1312 END DO 1313 ENDIF 1314 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1315 1316 1317 DO jk = 1, jpkm1 1318 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1319 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1320 END DO 1321 esigw3(ji,jj,1) = esigw3(ji,jj,2) 1322 esigt3(ji,jj,jpk) = esigt3(ji,jj,jpkm1) 1323 1324 ! Coefficients for vertical depth as the sum of e3w scale factors 1325 gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 1326 DO jk = 2, jpk 1327 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1328 END DO 1329 1330 DO jk = 1, jpk 1331 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp) 1332 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp) 1333 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigt3(ji,jj,jk)+hc*zcoeft) 1334 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsigw3(ji,jj,jk)+hc*zcoefw) 1335 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-hc)*gsi3w3(ji,jj,jk)+hc*zcoefw) 1336 END DO 1337 1338 ENDDO ! for all jj's 1339 ENDDO ! for all ji's 1340 1341 1342 DO ji=1,jpi 1343 DO jj=1,jpj 1344 1345 DO jk = 1, jpk 1346 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 1347 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1348 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 1349 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1350 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) + & 1351 hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 1352 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1353 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 1354 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1355 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 1356 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1357 1358 e3t(ji,jj,jk)=((hbatt(ji,jj)-hc)*esigt3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1359 e3u(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigtu3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1360 e3v(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigtv3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1361 e3f(ji,jj,jk)=((hbatf(ji,jj)-hc)*esigtf3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1362 ! 1363 e3w (ji,jj,jk)=((hbatt(ji,jj)-hc)*esigw3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1364 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hc)*esigwu3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1365 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hc)*esigwv3(ji,jj,jk) + hc/FLOAT(jpkm1)) 1366 END DO 1367 1368 ENDDO 1369 ENDDO 1370 1371 ELSE ! not ln_s_sigma 1372 1373 DO jk = 1, jpk 1374 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1375 gsigt(jk) = -fssig( REAL(jk,wp) ) 1376 END DO 1377 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1378 ! 1379 ! Coefficients for vertical scale factors at w-, t- levels 1246 1380 !!gm bug : define it from analytical function, not like juste bellow.... 1247 1381 !!gm or betteroffer the 2 possibilities.... 1248 DO jk = 1, jpkm11249 esigt(jk ) = gsigw(jk+1) - gsigw(jk)1250 esigw(jk+1) = gsigt(jk+1) - gsigt(jk)1251 END DO1252 esigw( 1 ) = esigw( 2 )1253 esigt(jpk) = esigt(jpkm1)1382 DO jk = 1, jpkm1 1383 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1384 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1385 END DO 1386 esigw( 1 ) = esigw( 2 ) 1387 esigt(jpk) = esigt(jpkm1) 1254 1388 1255 1389 !!gm original form … … 1259 1393 !!org END DO 1260 1394 !!gm 1261 !1262 ! Coefficients for vertical depth as the sum of e3w scale factors1263 gsi3w(1) = 0.5 * esigw(1)1264 DO jk = 2, jpk1265 gsi3w(jk) = gsi3w(jk-1) + esigw(jk)1266 END DO1395 ! 1396 ! Coefficients for vertical depth as the sum of e3w scale factors 1397 gsi3w(1) = 0.5 * esigw(1) 1398 DO jk = 2, jpk 1399 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1400 END DO 1267 1401 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1268 DO jk = 1, jpk1269 zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)1270 zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)1271 gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)1272 gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)1273 gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw)1274 END DO1402 DO jk = 1, jpk 1403 zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 1404 zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 1405 gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 1406 gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 1407 gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 1408 END DO 1275 1409 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1276 DO jj = 1, jpj 1277 DO ji = 1, jpi 1278 DO jk = 1, jpk 1279 e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1280 e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1281 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1282 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 1283 ! 1284 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1285 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1286 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1287 END DO 1288 END DO 1289 END DO 1410 DO jj = 1, jpj 1411 DO ji = 1, jpi 1412 DO jk = 1, jpk 1413 e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1414 e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1415 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1416 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 1417 ! 1418 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1419 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1420 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1421 END DO 1422 END DO 1423 END DO 1424 1425 ENDIF ! ln_s_sigma 1426 1427 1290 1428 ! 1291 1429 ! HYBRID :
Note: See TracChangeset
for help on using the changeset viewer.