Changeset 1601 for trunk/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2009-08-11T12:09:19+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r1577 r1601 42 42 43 43 !!gm DOCTOR name for the namelist parameter... 44 ! !!! ** Namelist nam_zgr_sco ** 45 REAL(wp) :: sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m) 46 REAL(wp) :: sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 47 REAL(wp) :: theta = 6.0 ! surface control parameter (0<=theta<=20) 48 REAL(wp) :: thetb = 0.75 ! bottom control parameter (0<=thetb<= 1) 49 REAL(wp) :: r_max = 0.15 ! maximum cut-off r-value allowed (0<r_max<1) 44 ! !!! ** Namelist namzgr_sco ** 45 REAL(wp) :: rn_sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m) 46 REAL(wp) :: rn_sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 47 REAL(wp) :: rn_theta = 6.0 ! surface control parameter (0<=rn_theta<=20) 48 REAL(wp) :: rn_thetb = 0.75 ! bottom control parameter (0<=rn_thetb<= 1) 49 REAL(wp) :: rn_rmax = 0.15 ! maximum cut-off r-value allowed (0<rn_rmax<1) 50 LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 51 REAL(wp) :: rn_bb = 0.8 ! stretching parameter for song and haidvogel stretching 52 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 53 REAL(wp) :: rn_hc = 150. ! Critical depth for s-sigma coordinates 50 54 51 55 !! * Substitutions … … 79 83 INTEGER :: ioptio = 0 ! temporary integer 80 84 !! 81 NAMELIST/nam _zgr/ ln_zco, ln_zps, ln_sco82 !!---------------------------------------------------------------------- 83 84 REWIND ( numnam ) ! Read Namelist nam _zgr : vertical coordinate'85 READ ( numnam, nam _zgr )85 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 86 !!---------------------------------------------------------------------- 87 88 REWIND ( numnam ) ! Read Namelist namzgr : vertical coordinate' 89 READ ( numnam, namzgr ) 86 90 87 91 IF(lwp) THEN ! Control print … … 89 93 WRITE(numout,*) 'dom_zgr : vertical coordinate' 90 94 WRITE(numout,*) '~~~~~~~' 91 WRITE(numout,*) ' Namelist nam _zgr : set vertical coordinate'95 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 92 96 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 93 97 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps … … 232 236 ENDIF 233 237 238 !!gm BUG in s-coordinate this does not work! 234 239 ! deepest/shallowest W level Above/Bellow ~10m 235 240 zrefdep = 10. - ( 0.1*MINVAL(e3w_0) ) ! ref. depth with tolerance (10% of minimum layer thickness) 236 241 nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Bellow ~10m 237 242 nla10 = nlb10 - 1 ! deepest W level Above ~10m 243 !!gm end bug 238 244 239 245 IF(lwp) THEN ! control print … … 1001 1007 !!---------------------------------------------------------------------- 1002 1008 ! 1003 pf = ( TANH( theta * ( -(pk-0.5) / REAL(jpkm1) +thetb ) ) &1004 & - TANH( thetb *theta ) ) &1005 & * ( COSH( theta ) &1006 & + COSH( theta * ( 2.e0 *thetb - 1.e0 ) ) ) &1007 & / ( 2.e0 * SINH( theta ) )1009 pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) & 1010 & - TANH( rn_thetb * rn_theta ) ) & 1011 & * ( COSH( rn_theta ) & 1012 & + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) ) ) & 1013 & / ( 2.e0 * SINH( rn_theta ) ) 1008 1014 ! 1009 1015 END FUNCTION fssig 1010 1016 1011 1017 1012 FUNCTION fssig1( pk1, bb ) RESULT( pf1 )1018 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 1013 1019 !!---------------------------------------------------------------------- 1014 1020 !! *** ROUTINE eos_init *** … … 1024 1030 !!---------------------------------------------------------------------- 1025 1031 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 1026 REAL(wp), INTENT(in ) :: bb! Stretching coefficient1032 REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient 1027 1033 REAL(wp) :: pf1 ! sigma value 1028 1034 !!---------------------------------------------------------------------- 1029 1035 ! 1030 IF ( theta == 0 ) then ! uniform sigma1036 IF ( rn_theta == 0 ) then ! uniform sigma 1031 1037 pf1 = -(pk1-0.5) / REAL( jpkm1 ) 1032 1038 ELSE ! stretched sigma 1033 pf1 = (1.0- bb) * (sinh( theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(theta) + &1034 & bb * ( (tanh( theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*theta) ) / &1035 & (2*tanh(0.5* theta) ) )1039 pf1 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + & 1040 & pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) / & 1041 & (2*tanh(0.5*rn_theta) ) ) 1036 1042 ENDIF 1037 1043 ! … … 1077 1083 REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - - 1078 1084 !! 1079 LOGICAL :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true.1080 REAL(wp) :: bb = 0.8 ! stretching parameter for song and haidvogel stretching, bb=0; top only, bb =1; top and bottom1081 REAL(wp) :: hc = 150 ! Critical depth for s-sigma coordinates1082 1085 !!gm never do that !!!! ==> Pb at compilation phase on several computer 1083 1086 REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3 = 0.0d0 … … 1093 1096 !!gm end 1094 1097 !! 1095 NAMELIST/nam _zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max, ln_s_sigma, bb,hc1096 !!---------------------------------------------------------------------- 1097 1098 REWIND( numnam ) ! Read Namelist nam _zgr_sco : sigma-stretching parameters1099 READ ( numnam, nam _zgr_sco )1098 NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 1099 !!---------------------------------------------------------------------- 1100 1101 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters 1102 READ ( numnam, namzgr_sco ) 1100 1103 1101 1104 IF(lwp) THEN ! control print … … 1103 1106 WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate' 1104 1107 WRITE(numout,*) '~~~~~~~~~~~' 1105 WRITE(numout,*) ' Namelist nam_zgr_sco' 1106 WRITE(numout,*) ' sigma-stretching coeffs ' 1107 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max 1108 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min 1109 WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta 1110 WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb 1111 WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max 1112 WRITE(numout,*) ' Critical depth hc = ', hc 1113 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1114 ENDIF 1115 1116 hift(:,:) = sbot_min ! set the minimum depth for the s-coordinate 1117 hifu(:,:) = sbot_min 1118 hifv(:,:) = sbot_min 1119 hiff(:,:) = sbot_min 1108 WRITE(numout,*) ' Namelist namzgr_sco' 1109 WRITE(numout,*) ' sigma-stretching coeffs ' 1110 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max 1111 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min 1112 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta 1113 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb 1114 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax 1115 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1116 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb 1117 WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc 1118 ENDIF 1119 1120 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate 1121 hifu(:,:) = rn_sbot_min 1122 hifv(:,:) = rn_sbot_min 1123 hiff(:,:) = rn_sbot_min 1120 1124 1121 1125 ! ! set maximum ocean depth 1122 bathy(:,:) = MIN( sbot_max, bathy(:,:) )1126 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1123 1127 1124 1128 DO jj = 1, jpj 1125 1129 DO ji = 1, jpi 1126 1130 IF (bathy(ji,jj) .gt. 0.0) THEN 1127 bathy(ji,jj) = MAX( sbot_min, bathy(ji,jj) )1131 bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1128 1132 ENDIF 1129 1133 END DO … … 1139 1143 DO jj = 1, jpj 1140 1144 DO ji = 1, jpi 1141 zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )1145 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1142 1146 END DO 1143 1147 END DO … … 1145 1149 zrmax = 1.e0 1146 1150 ! ! ================ ! 1147 DO WHILE ( jl <= 10000 .AND. zrmax > r _max ) ! Iterative loop !1151 DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1148 1152 ! ! ================ ! 1149 1153 jl = jl + 1 … … 1157 1161 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1158 1162 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1159 IF( zri(ji,jj) > r _max ) zmsk(ji ,jj ) = 1.01160 IF( zri(ji,jj) > r _max ) zmsk(iip1,jj ) = 1.01161 IF( zrj(ji,jj) > r _max ) zmsk(ji ,jj ) = 1.01162 IF( zrj(ji,jj) > r _max ) zmsk(ji ,ijp1) = 1.01163 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0 1164 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1.0 1165 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0 1166 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1.0 1163 1167 END DO 1164 1168 END DO … … 1218 1222 DO ji = 1, jpi 1219 1223 ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 1220 hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)1224 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 1221 1225 END DO 1222 1226 END DO … … 1239 1243 IF(lwp) THEN 1240 1244 WRITE(numout,*) 1241 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min1242 ENDIF 1243 hbatu(:,:) = sbot_min1244 hbatv(:,:) = sbot_min1245 hbatf(:,:) = sbot_min1245 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1246 ENDIF 1247 hbatu(:,:) = rn_sbot_min 1248 hbatv(:,:) = rn_sbot_min 1249 hbatf(:,:) = rn_sbot_min 1246 1250 DO jj = 1, jpjm1 1247 1251 DO ji = 1, jpim1 … … 1259 1263 DO ji = 1, jpi 1260 1264 IF( hbatu(ji,jj) == 0.e0 ) THEN 1261 IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min1265 IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = rn_sbot_min 1262 1266 IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj) 1263 1267 ENDIF … … 1268 1272 DO ji = 1, jpi 1269 1273 IF( hbatv(ji,jj) == 0.e0 ) THEN 1270 IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min1274 IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = rn_sbot_min 1271 1275 IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj) 1272 1276 ENDIF … … 1277 1281 DO ji = 1, jpi 1278 1282 IF( hbatf(ji,jj) == 0.e0 ) THEN 1279 IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min1283 IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = rn_sbot_min 1280 1284 IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj) 1281 1285 ENDIF … … 1307 1311 ! non-dimensional "sigma" for model level depth at w- and t-levels 1308 1312 1309 IF ( ln_s_sigma ) THEN !Song and Haidvogel style stretched sigma for depths below hc, with uniform sigma in shallower waters1310 1311 DO ji =1,jpi1312 DO jj=1,jpj1313 1314 IF (hbatt(ji,jj).GT. hc) THEN !deep water, stretched sigma1313 IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths 1314 ! ! below rn_hc, with uniform sigma in shallower waters 1315 DO ji = 1, jpi 1316 DO jj = 1, jpj 1317 1318 IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma 1315 1319 DO jk = 1, jpk 1316 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, bb )1317 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , bb )1320 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1321 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1318 1322 END DO 1319 1323 ELSE ! shallow water, uniform sigma … … 1342 1346 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp) 1343 1347 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp) 1344 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)- hc)*gsigt3(ji,jj,jk)+hc*zcoeft)1345 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)- hc)*gsigw3(ji,jj,jk)+hc*zcoefw)1346 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)- hc)*gsi3w3(ji,jj,jk)+hc*zcoefw)1348 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft) 1349 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw) 1350 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoefw) 1347 1351 END DO 1348 1352 … … 1367 1371 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1368 1372 1369 e3t(ji,jj,jk)=((hbatt(ji,jj)- hc)*esigt3(ji,jj,jk) +hc/FLOAT(jpkm1))1370 e3u(ji,jj,jk)=((hbatu(ji,jj)- hc)*esigtu3(ji,jj,jk) +hc/FLOAT(jpkm1))1371 e3v(ji,jj,jk)=((hbatv(ji,jj)- hc)*esigtv3(ji,jj,jk) +hc/FLOAT(jpkm1))1372 e3f(ji,jj,jk)=((hbatf(ji,jj)- hc)*esigtf3(ji,jj,jk) +hc/FLOAT(jpkm1))1373 e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1374 e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1375 e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1376 e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1373 1377 ! 1374 e3w (ji,jj,jk)=((hbatt(ji,jj)- hc)*esigw3(ji,jj,jk) +hc/FLOAT(jpkm1))1375 e3uw(ji,jj,jk)=((hbatu(ji,jj)- hc)*esigwu3(ji,jj,jk) +hc/FLOAT(jpkm1))1376 e3vw(ji,jj,jk)=((hbatv(ji,jj)- hc)*esigwv3(ji,jj,jk) +hc/FLOAT(jpkm1))1378 e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1379 e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1380 e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1377 1381 END DO 1378 1382
Note: See TracChangeset
for help on using the changeset viewer.
