Changeset 5120 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2015-03-03T17:11:55+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5118 r5120 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 19 20 !!---------------------------------------------------------------------- 20 21 … … 35 36 USE oce ! ocean variables 36 37 USE dom_oce ! ocean domain 37 USE sbc_oce ! surface variable (isf)38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration … … 298 298 ENDIF 299 299 300 IF ( ln_isfcav ) THEN 300 301 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 301 302 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 302 DO jk = 1, jpkm1 303 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 304 END DO 305 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 306 307 DO jk = 2, jpk 308 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 309 END DO 310 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 303 DO jk = 1, jpkm1 304 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 305 END DO 306 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 307 308 DO jk = 2, jpk 309 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 310 END DO 311 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 312 END IF 311 313 312 314 !!gm BUG in s-coordinate this does not work! … … 472 474 ! 473 475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 474 IF( cp_cfg == "isomip" ) THEN 475 ! 476 risfdep(:,:)=200.e0 477 misfdep(:,:)=1 478 ij0 = 1 ; ij1 = 40 479 DO jj = mj0(ij0), mj1(ij1) 480 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 481 END DO 476 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 477 risfdep(:,:)=200.e0 478 misfdep(:,:)=1 479 ij0 = 1 ; ij1 = 40 480 DO jj = mj0(ij0), mj1(ij1) 481 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 482 END DO 482 483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 483 484 ELSEIF ( cp_cfg == "isomip2" ) THEN484 ! 485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 485 486 ! 486 487 risfdep(:,:)=0.e0 … … 540 541 END IF 541 542 CALL iom_close( inum ) 542 ! 543 ! 543 544 risfdep(:,:)=0._wp 544 545 misfdep(:,:)=1 … … 588 589 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 589 590 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 590 WHERE (bathy == risfdep) 591 bathy = 0.0_wp ; risfdep = 0.0_wp 592 END WHERE 591 IF ( ln_isfcav ) THEN 592 WHERE (bathy == risfdep) 593 bathy = 0.0_wp ; risfdep = 0.0_wp 594 END WHERE 595 END IF 593 596 ! end patch 594 597 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level … … 965 968 !!---------------------------------------------------------------------- 966 969 !! 970 INTEGER :: ji, jj, jk ! dummy loop indices 971 INTEGER :: ik, it ! temporary integers 972 LOGICAL :: ll_print ! Allow control print for debugging 973 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 974 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 975 REAL(wp) :: zmax ! Maximum depth 976 REAL(wp) :: zdiff ! temporary scalar 977 REAL(wp) :: zrefdep ! temporary scalar 978 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 979 !!--------------------------------------------------------------------- 980 ! 981 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 982 ! 983 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 984 ! 985 IF(lwp) WRITE(numout,*) 986 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 987 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 988 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 989 990 ll_print = .FALSE. ! Local variable for debugging 991 992 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 993 WRITE(numout,*) 994 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 995 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 996 ENDIF 997 998 999 ! bathymetry in level (from bathy_meter) 1000 ! =================== 1001 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 1002 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 1003 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1004 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1005 END WHERE 1006 1007 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1008 ! find the number of ocean levels such that the last level thickness 1009 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1010 ! e3t_1d is the reference level thickness 1011 DO jk = jpkm1, 1, -1 1012 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1013 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1014 END DO 1015 1016 IF ( ln_isfcav ) CALL zgr_isf 1017 1018 ! Scale factors and depth at T- and W-points 1019 DO jk = 1, jpk ! intitialization to the reference z-coordinate 1020 gdept_0(:,:,jk) = gdept_1d(jk) 1021 gdepw_0(:,:,jk) = gdepw_1d(jk) 1022 e3t_0 (:,:,jk) = e3t_1d (jk) 1023 e3w_0 (:,:,jk) = e3w_1d (jk) 1024 END DO 1025 ! 1026 DO jj = 1, jpj 1027 DO ji = 1, jpi 1028 ik = mbathy(ji,jj) 1029 IF( ik > 0 ) THEN ! ocean point only 1030 ! max ocean level case 1031 IF( ik == jpkm1 ) THEN 1032 zdepwp = bathy(ji,jj) 1033 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1034 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1035 e3t_0(ji,jj,ik ) = ze3tp 1036 e3t_0(ji,jj,ik+1) = ze3tp 1037 e3w_0(ji,jj,ik ) = ze3wp 1038 e3w_0(ji,jj,ik+1) = ze3tp 1039 gdepw_0(ji,jj,ik+1) = zdepwp 1040 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1041 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1042 ! 1043 ELSE ! standard case 1044 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1045 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1046 ENDIF 1047 !gm Bug? check the gdepw_1d 1048 ! ... on ik 1049 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1050 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1051 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1052 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1053 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1054 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1055 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1056 ! ... on ik+1 1057 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1058 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1059 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1060 ENDIF 1061 ENDIF 1062 END DO 1063 END DO 1064 ! 1065 it = 0 1066 DO jj = 1, jpj 1067 DO ji = 1, jpi 1068 ik = mbathy(ji,jj) 1069 IF( ik > 0 ) THEN ! ocean point only 1070 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1071 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1072 ! test 1073 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1074 IF( zdiff <= 0._wp .AND. lwp ) THEN 1075 it = it + 1 1076 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1077 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1078 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1079 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1080 ENDIF 1081 ENDIF 1082 END DO 1083 END DO 1084 ! 1085 IF ( ln_isfcav ) THEN 1086 ! (ISF) Definition of e3t, u, v, w for ISF case 1087 DO jj = 1, jpj 1088 DO ji = 1, jpi 1089 ik = misfdep(ji,jj) 1090 IF( ik > 1 ) THEN ! ice shelf point only 1091 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1092 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1093 !gm Bug? check the gdepw_0 1094 ! ... on ik 1095 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1096 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1097 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1098 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1099 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1100 1101 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1102 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1103 ENDIF 1104 ! ... on ik / ik-1 1105 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1106 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1107 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1108 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1109 ENDIF 1110 END DO 1111 END DO 1112 ! 1113 it = 0 1114 DO jj = 1, jpj 1115 DO ji = 1, jpi 1116 ik = misfdep(ji,jj) 1117 IF( ik > 1 ) THEN ! ice shelf point only 1118 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1119 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1120 ! test 1121 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1122 IF( zdiff <= 0. .AND. lwp ) THEN 1123 it = it + 1 1124 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1125 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1126 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1127 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1128 ENDIF 1129 ENDIF 1130 END DO 1131 END DO 1132 END IF 1133 ! END (ISF) 1134 1135 ! Scale factors and depth at U-, V-, UW and VW-points 1136 DO jk = 1, jpk ! initialisation to z-scale factors 1137 e3u_0 (:,:,jk) = e3t_1d(jk) 1138 e3v_0 (:,:,jk) = e3t_1d(jk) 1139 e3uw_0(:,:,jk) = e3w_1d(jk) 1140 e3vw_0(:,:,jk) = e3w_1d(jk) 1141 END DO 1142 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1143 DO jj = 1, jpjm1 1144 DO ji = 1, fs_jpim1 ! vector opt. 1145 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1146 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1147 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1148 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1149 END DO 1150 END DO 1151 END DO 1152 IF ( ln_isfcav ) THEN 1153 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1154 ! Need to test if the modification of only mikt and mbkt levels is enough 1155 DO jk = 2,jpk 1156 DO jj = 1, jpjm1 1157 DO ji = 1, fs_jpim1 ! vector opt. 1158 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1159 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1160 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1161 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1162 END DO 1163 END DO 1164 END DO 1165 END IF 1166 1167 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1168 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1169 ! 1170 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1171 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1172 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1173 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1174 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1175 END DO 1176 1177 ! Scale factor at F-point 1178 DO jk = 1, jpk ! initialisation to z-scale factors 1179 e3f_0(:,:,jk) = e3t_1d(jk) 1180 END DO 1181 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1182 DO jj = 1, jpjm1 1183 DO ji = 1, fs_jpim1 ! vector opt. 1184 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1185 END DO 1186 END DO 1187 END DO 1188 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1189 ! 1190 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1191 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1192 END DO 1193 !!gm bug ? : must be a do loop with mj0,mj1 1194 ! 1195 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1196 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1197 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1198 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1199 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1200 1201 ! Control of the sign 1202 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1203 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1204 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1205 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1206 1207 ! Compute gdep3w_0 (vertical sum of e3w) 1208 IF ( ln_isfcav ) THEN ! if cavity 1209 WHERE (misfdep == 0) misfdep = 1 1210 DO jj = 1,jpj 1211 DO ji = 1,jpi 1212 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1213 DO jk = 2, misfdep(ji,jj) 1214 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1215 END DO 1216 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1217 DO jk = misfdep(ji,jj) + 1, jpk 1218 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1219 END DO 1220 END DO 1221 END DO 1222 ELSE ! no cavity 1223 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1224 DO jk = 2, jpk 1225 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1226 END DO 1227 END IF 1228 ! ! ================= ! 1229 IF(lwp .AND. ll_print) THEN ! Control print ! 1230 ! ! ================= ! 1231 DO jj = 1,jpj 1232 DO ji = 1, jpi 1233 ik = MAX( mbathy(ji,jj), 1 ) 1234 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1235 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1236 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1237 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1238 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1239 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1240 END DO 1241 END DO 1242 WRITE(numout,*) 1243 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1244 WRITE(numout,*) 1245 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1246 WRITE(numout,*) 1247 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1248 WRITE(numout,*) 1249 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1250 WRITE(numout,*) 1251 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1252 WRITE(numout,*) 1253 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1254 ENDIF 1255 ! 1256 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1257 ! 1258 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1259 ! 1260 END SUBROUTINE zgr_zps 1261 1262 SUBROUTINE zgr_isf 1263 !!---------------------------------------------------------------------- 1264 !! *** ROUTINE zgr_isf *** 1265 !! 1266 !! ** Purpose : check the bathymetry in levels 1267 !! 1268 !! ** Method : THe water column have to contained at least 2 cells 1269 !! Bathymetry and isfdraft are modified (dig/close) to respect 1270 !! this criterion. 1271 !! 1272 !! 1273 !! ** Action : - test compatibility between isfdraft and bathy 1274 !! - bathy and isfdraft are modified 1275 !!---------------------------------------------------------------------- 1276 !! 967 1277 INTEGER :: ji, jj, jk, jl ! dummy loop indices 968 1278 INTEGER :: ik, it ! temporary integers … … 975 1285 REAL(wp) :: zdiff ! temporary scalar 976 1286 REAL(wp) :: zrefdep ! temporary scalar 977 REAL(wp) :: zbathydiff, zrisfdepdiff 978 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 3D workspace (ISH) 979 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 3D workspace (ISH) 980 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 1287 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1288 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1289 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 981 1290 !!--------------------------------------------------------------------- 982 1291 ! 983 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 984 ! 985 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 1292 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1293 ! 986 1294 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 987 CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 988 ! 989 IF(lwp) WRITE(numout,*) 990 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 991 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 992 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 993 994 ll_print = .FALSE. ! Local variable for debugging 995 996 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 997 WRITE(numout,*) 998 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 999 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 1000 ENDIF 1001 1002 ! bathymetry in level (from bathy_meter) 1003 ! =================== 1004 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 1005 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 1006 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1007 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1008 END WHERE 1009 1010 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1011 ! find the number of ocean levels such that the last level thickness 1012 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1013 ! e3t_1d is the reference level thickness 1014 DO jk = jpkm1, 1, -1 1015 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1016 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1017 END DO 1295 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1296 1297 1018 1298 ! (ISF) compute misfdep 1019 1299 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 … … 1059 1339 misfdep(jpi,:) = misfdep( 2 ,:) 1060 1340 ENDIF 1061 1341 1062 1342 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1063 1343 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1064 1344 mbathy(jpi,:) = mbathy( 2 ,:) 1065 1345 ENDIF 1066 1346 1067 1347 ! split last cell if possible (only where water column is 2 cell or less) 1068 1348 DO jk = jpkm1, 1, -1 … … 1082 1362 END WHERE 1083 1363 END DO 1084 1364 1085 1365 1086 1366 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition … … 1363 1643 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1364 1644 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1365 IF( ibtest == 0 ) THEN1645 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 1366 1646 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1367 1647 END IF … … 1479 1759 ENDIF 1480 1760 1481 ! Scale factors and depth at T- and W-points1482 DO jk = 1, jpk ! intitialization to the reference z-coordinate1483 gdept_0(:,:,jk) = gdept_1d(jk)1484 gdepw_0(:,:,jk) = gdepw_1d(jk)1485 e3t_0 (:,:,jk) = e3t_1d (jk)1486 e3w_0 (:,:,jk) = e3w_1d (jk)1487 END DO1488 !1489 DO jj = 1, jpj1490 DO ji = 1, jpi1491 ik = mbathy(ji,jj)1492 IF( ik > 0 ) THEN ! ocean point only1493 ! max ocean level case1494 IF( ik == jpkm1 ) THEN1495 zdepwp = bathy(ji,jj)1496 ze3tp = bathy(ji,jj) - gdepw_1d(ik)1497 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )1498 e3t_0(ji,jj,ik ) = ze3tp1499 e3t_0(ji,jj,ik+1) = ze3tp1500 e3w_0(ji,jj,ik ) = ze3wp1501 e3w_0(ji,jj,ik+1) = ze3tp1502 gdepw_0(ji,jj,ik+1) = zdepwp1503 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp1504 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp1505 !1506 ELSE ! standard case1507 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)1508 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)1509 ENDIF1510 !gm Bug? check the gdepw_1d1511 ! ... on ik1512 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1513 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1514 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1515 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1516 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1517 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1518 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1519 ! ... on ik+11520 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1521 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1522 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1523 ENDIF1524 ENDIF1525 END DO1526 END DO1527 !1528 it = 01529 DO jj = 1, jpj1530 DO ji = 1, jpi1531 ik = mbathy(ji,jj)1532 IF( ik > 0 ) THEN ! ocean point only1533 e3tp (ji,jj) = e3t_0(ji,jj,ik)1534 e3wp (ji,jj) = e3w_0(ji,jj,ik)1535 ! test1536 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1537 IF( zdiff <= 0._wp .AND. lwp ) THEN1538 it = it + 11539 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1540 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1541 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1542 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1543 ENDIF1544 ENDIF1545 END DO1546 END DO1547 !1548 ! (ISF) Definition of e3t, u, v, w for ISF case1549 DO jj = 1, jpj1550 DO ji = 1, jpi1551 ik = misfdep(ji,jj)1552 IF( ik > 1 ) THEN ! ice shelf point only1553 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)1554 gdepw_0(ji,jj,ik) = risfdep(ji,jj)1555 !gm Bug? check the gdepw_01556 ! ... on ik1557 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &1558 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &1559 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) )1560 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)1561 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)1562 1563 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)1564 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)1565 ENDIF1566 ! ... on ik / ik-11567 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1568 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)1569 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code1570 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)1571 ENDIF1572 END DO1573 END DO1574 !1575 it = 01576 DO jj = 1, jpj1577 DO ji = 1, jpi1578 ik = misfdep(ji,jj)1579 IF( ik > 1 ) THEN ! ice shelf point only1580 e3tp (ji,jj) = e3t_0(ji,jj,ik )1581 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1582 ! test1583 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1584 IF( zdiff <= 0. .AND. lwp ) THEN1585 it = it + 11586 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1587 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1588 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1589 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1590 ENDIF1591 ENDIF1592 END DO1593 END DO1594 ! END (ISF)1595 1596 ! Scale factors and depth at U-, V-, UW and VW-points1597 DO jk = 1, jpk ! initialisation to z-scale factors1598 e3u_0 (:,:,jk) = e3t_1d(jk)1599 e3v_0 (:,:,jk) = e3t_1d(jk)1600 e3uw_0(:,:,jk) = e3w_1d(jk)1601 e3vw_0(:,:,jk) = e3w_1d(jk)1602 END DO1603 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1604 DO jj = 1, jpjm11605 DO ji = 1, fs_jpim1 ! vector opt.1606 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1607 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1608 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1609 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1610 END DO1611 END DO1612 END DO1613 ! (ISF) define e3uw1614 DO jk = 2,jpk1615 DO jj = 1, jpjm11616 DO ji = 1, fs_jpim1 ! vector opt.1617 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) &1618 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) )1619 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) &1620 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) )1621 END DO1622 END DO1623 END DO1624 !End (ISF)1625 1626 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1627 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1628 !1629 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1630 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1631 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1632 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1633 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1634 END DO1635 1636 ! Scale factor at F-point1637 DO jk = 1, jpk ! initialisation to z-scale factors1638 e3f_0(:,:,jk) = e3t_1d(jk)1639 END DO1640 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1641 DO jj = 1, jpjm11642 DO ji = 1, fs_jpim1 ! vector opt.1643 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1644 END DO1645 END DO1646 END DO1647 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1648 !1649 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1650 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1651 END DO1652 !!gm bug ? : must be a do loop with mj0,mj11653 !1654 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21655 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1656 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1657 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1658 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1659 1660 ! Control of the sign1661 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1662 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1663 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1664 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1665 1666 ! Compute gdep3w_0 (vertical sum of e3w)1667 WHERE (misfdep == 0) misfdep = 11668 DO jj = 1,jpj1669 DO ji = 1,jpi1670 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1671 DO jk = 2, misfdep(ji,jj)1672 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1673 END DO1674 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1675 DO jk = misfdep(ji,jj) + 1, jpk1676 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1677 END DO1678 END DO1679 END DO1680 ! ! ================= !1681 IF(lwp .AND. ll_print) THEN ! Control print !1682 ! ! ================= !1683 DO jj = 1,jpj1684 DO ji = 1, jpi1685 ik = MAX( mbathy(ji,jj), 1 )1686 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1687 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1688 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1689 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1690 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1691 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1692 END DO1693 END DO1694 WRITE(numout,*)1695 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1696 WRITE(numout,*)1697 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1698 WRITE(numout,*)1699 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1700 WRITE(numout,*)1701 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1702 WRITE(numout,*)1703 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1704 WRITE(numout,*)1705 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1706 ENDIF1707 !1708 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1709 1761 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1710 1762 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1711 ! 1712 IF( nn_timing == 1 ) CALL timing_stop('zgr_ zps')1713 !1714 END SUBROUTINE zgr_zps1763 1764 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1765 1766 END SUBROUTINE 1715 1767 1716 1768 SUBROUTINE zgr_sco
Note: See TracChangeset
for help on using the changeset viewer.
