Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4687 r5208 35 35 USE oce ! ocean variables 36 36 USE dom_oce ! ocean domain 37 USE sbc_oce ! surface variable (isf) 37 38 USE closea ! closed seas 38 39 USE c1d ! 1D vertical configuration … … 101 102 INTEGER :: ios 102 103 ! 103 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 105 !!---------------------------------------------------------------------- 105 106 ! … … 120 121 WRITE(numout,*) '~~~~~~~' 121 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 122 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 123 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 124 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 125 127 ENDIF 126 128 … … 145 147 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 146 148 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 149 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 147 150 ! 148 151 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain … … 294 297 gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero 295 298 ENDIF 299 300 ! 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 ! 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)) 296 311 297 312 !!gm BUG in s-coordinate this does not work! … … 450 465 END DO 451 466 END DO 467 risfdep(:,:)=0.e0 468 misfdep(:,:)=1 469 ! 470 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 471 IF( cp_cfg == "isomip" ) THEN 472 ! 473 risfdep(:,:)=200.e0 474 misfdep(:,:)=1 475 ij0 = 1 ; ij1 = 40 476 DO jj = mj0(ij0), mj1(ij1) 477 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 478 END DO 479 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 480 ! 481 ELSEIF ( cp_cfg == "isomip2" ) THEN 482 ! 483 risfdep(:,:)=0.e0 484 misfdep(:,:)=1 485 ij0 = 1 ; ij1 = 40 486 DO jj = mj0(ij0), mj1(ij1) 487 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 488 END DO 489 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 490 END IF 452 491 ! 453 492 ! ! ================ ! … … 492 531 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 493 532 CALL iom_close( inum ) 494 ! 533 ! 534 risfdep(:,:)=0._wp 535 misfdep(:,:)=1 536 IF ( ln_isfcav ) THEN 537 CALL iom_open ( 'isf_draft_meter.nc', inum ) 538 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 539 CALL iom_close( inum ) 540 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 541 END IF 542 ! 495 543 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 496 544 ! … … 530 578 ! 531 579 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 580 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 581 WHERE (bathy == risfdep) 582 bathy = 0.0_wp ; risfdep = 0.0_wp 583 END WHERE 584 ! end patch 532 585 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 533 586 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 783 836 END SUBROUTINE zgr_bot_level 784 837 838 SUBROUTINE zgr_top_level 839 !!---------------------------------------------------------------------- 840 !! *** ROUTINE zgr_bot_level *** 841 !! 842 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 843 !! 844 !! ** Method : computes from misfdep with a minimum value of 1 845 !! 846 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 847 !! ocean level at t-, u- & v-points 848 !! (min value = 1) 849 !!---------------------------------------------------------------------- 850 !! 851 INTEGER :: ji, jj ! dummy loop indices 852 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 853 !!---------------------------------------------------------------------- 854 ! 855 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 856 ! 857 CALL wrk_alloc( jpi, jpj, zmik ) 858 ! 859 IF(lwp) WRITE(numout,*) 860 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 861 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 862 ! 863 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 864 ! ! top k-index of W-level (=mikt) 865 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 866 DO ji = 1, jpim1 867 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 868 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 869 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 870 END DO 871 END DO 872 873 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 874 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 875 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 876 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 877 ! 878 CALL wrk_dealloc( jpi, jpj, zmik ) 879 ! 880 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 881 ! 882 END SUBROUTINE zgr_top_level 785 883 786 884 SUBROUTINE zgr_zco … … 861 959 !!---------------------------------------------------------------------- 862 960 !! 863 INTEGER :: ji, jj, jk 961 INTEGER :: ji, jj, jk, jl ! dummy loop indices 864 962 INTEGER :: ik, it ! temporary integers 963 INTEGER :: id, jd, nprocd 964 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 865 965 LOGICAL :: ll_print ! Allow control print for debugging 866 966 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 867 967 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 868 REAL(wp) :: zmax ! Maximum depth968 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 869 969 REAL(wp) :: zdiff ! temporary scalar 870 970 REAL(wp) :: zrefdep ! temporary scalar 971 REAL(wp) :: zbathydiff, zrisfdepdiff 972 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 3D workspace (ISH) 973 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 3D workspace (ISH) 871 974 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 872 975 !!--------------------------------------------------------------------- … … 875 978 ! 876 979 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 980 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 981 CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 877 982 ! 878 983 IF(lwp) WRITE(numout,*) … … 889 994 ENDIF 890 995 891 892 996 ! bathymetry in level (from bathy_meter) 893 997 ! =================== … … 906 1010 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 907 1011 END DO 1012 ! (ISF) compute misfdep 1013 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1014 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1015 END WHERE 1016 1017 ! Compute misfdep for ocean points (i.e. first wet level) 1018 ! find the first ocean level such that the first level thickness 1019 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1020 ! e3t_0 is the reference level thickness 1021 DO jk = 2, jpkm1 1022 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1023 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1024 END DO 1025 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 1026 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1027 END WHERE 1028 1029 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1030 icompt = 0 1031 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1032 DO jl = 1, 10 1033 WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 1034 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1035 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1036 END WHERE 1037 WHERE (mbathy(:,:) <= 0) 1038 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1039 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1040 ENDWHERE 1041 IF( lk_mpp ) THEN 1042 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1043 CALL lbc_lnk( zbathy, 'T', 1. ) 1044 misfdep(:,:) = INT( zbathy(:,:) ) 1045 CALL lbc_lnk( risfdep, 'T', 1. ) 1046 CALL lbc_lnk( bathy, 'T', 1. ) 1047 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1048 CALL lbc_lnk( zbathy, 'T', 1. ) 1049 mbathy(:,:) = INT( zbathy(:,:) ) 1050 ENDIF 1051 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1052 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1053 misfdep(jpi,:) = misfdep( 2 ,:) 1054 ENDIF 1055 1056 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1057 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1058 mbathy(jpi,:) = mbathy( 2 ,:) 1059 ENDIF 1060 1061 ! split last cell if possible (only where water column is 2 cell or less) 1062 DO jk = jpkm1, 1, -1 1063 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1064 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1065 mbathy(:,:) = jk 1066 bathy(:,:) = zmax 1067 END WHERE 1068 END DO 1069 1070 ! split top cell if possible (only where water column is 2 cell or less) 1071 DO jk = 2, jpkm1 1072 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1073 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1074 misfdep(:,:) = jk 1075 risfdep(:,:) = zmax 1076 END WHERE 1077 END DO 1078 1079 1080 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1081 DO jj = 1, jpj 1082 DO ji = 1, jpi 1083 ! find the minimum change option: 1084 ! test bathy 1085 IF (risfdep(ji,jj) .GT. 1) THEN 1086 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1087 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1088 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1089 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1090 1091 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1092 IF (zbathydiff .LE. zrisfdepdiff) THEN 1093 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1094 mbathy(ji,jj)= mbathy(ji,jj) + 1 1095 ELSE 1096 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1097 misfdep(ji,jj) = misfdep(ji,jj) - 1 1098 END IF 1099 END IF 1100 END IF 1101 END DO 1102 END DO 1103 1104 ! At least 2 levels for water thickness at T, U, and V point. 1105 DO jj = 1, jpj 1106 DO ji = 1, jpi 1107 ! find the minimum change option: 1108 ! test bathy 1109 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1110 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1111 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1112 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1113 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1114 IF (zbathydiff .LE. zrisfdepdiff) THEN 1115 mbathy(ji,jj) = mbathy(ji,jj) + 1 1116 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1117 ELSE 1118 misfdep(ji,jj)= misfdep(ji,jj) - 1 1119 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1120 END IF 1121 ENDIF 1122 END DO 1123 END DO 1124 1125 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 1126 DO jj = 1, jpjm1 1127 DO ji = 1, jpim1 1128 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1129 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1130 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1131 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1132 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1133 IF (zbathydiff .LE. zrisfdepdiff) THEN 1134 mbathy(ji,jj) = mbathy(ji,jj) + 1 1135 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & 1136 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1137 ELSE 1138 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1139 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 1140 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1141 END IF 1142 ENDIF 1143 END DO 1144 END DO 1145 1146 IF( lk_mpp ) THEN 1147 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1148 CALL lbc_lnk( zbathy, 'T', 1. ) 1149 misfdep(:,:) = INT( zbathy(:,:) ) 1150 CALL lbc_lnk( risfdep, 'T', 1. ) 1151 CALL lbc_lnk( bathy, 'T', 1. ) 1152 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1153 CALL lbc_lnk( zbathy, 'T', 1. ) 1154 mbathy(:,:) = INT( zbathy(:,:) ) 1155 ENDIF 1156 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 1157 DO jj = 1, jpjm1 1158 DO ji = 1, jpim1 1159 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 1160 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1161 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1162 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1163 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1164 IF (zbathydiff .LE. zrisfdepdiff) THEN 1165 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1166 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1167 ELSE 1168 misfdep(ji,jj) = misfdep(ji,jj) - 1 1169 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1170 END IF 1171 ENDIF 1172 END DO 1173 END DO 1174 1175 1176 IF( lk_mpp ) THEN 1177 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1178 CALL lbc_lnk( zbathy, 'T', 1. ) 1179 misfdep(:,:) = INT( zbathy(:,:) ) 1180 CALL lbc_lnk( risfdep, 'T', 1. ) 1181 CALL lbc_lnk( bathy, 'T', 1. ) 1182 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1183 CALL lbc_lnk( zbathy, 'T', 1. ) 1184 mbathy(:,:) = INT( zbathy(:,:) ) 1185 ENDIF 1186 1187 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 1188 DO jj = 1, jpjm1 1189 DO ji = 1, jpim1 1190 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1191 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1192 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1193 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1194 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1195 IF (zbathydiff .LE. zrisfdepdiff) THEN 1196 mbathy(ji,jj) = mbathy(ji,jj) + 1 1197 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1198 ELSE 1199 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1200 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1201 END IF 1202 ENDIF 1203 ENDDO 1204 ENDDO 1205 1206 IF( lk_mpp ) THEN 1207 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1208 CALL lbc_lnk( zbathy, 'T', 1. ) 1209 misfdep(:,:) = INT( zbathy(:,:) ) 1210 CALL lbc_lnk( risfdep, 'T', 1. ) 1211 CALL lbc_lnk( bathy, 'T', 1. ) 1212 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1213 CALL lbc_lnk( zbathy, 'T', 1. ) 1214 mbathy(:,:) = INT( zbathy(:,:) ) 1215 ENDIF 1216 1217 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 1218 DO jj = 1, jpjm1 1219 DO ji = 1, jpim1 1220 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1221 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1222 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1223 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1224 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1225 IF (zbathydiff .LE. zrisfdepdiff) THEN 1226 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1227 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & 1228 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1229 ELSE 1230 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1231 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) & 1232 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1233 END IF 1234 ENDIF 1235 ENDDO 1236 ENDDO 1237 1238 IF( lk_mpp ) THEN 1239 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1240 CALL lbc_lnk( zbathy, 'T', 1. ) 1241 misfdep(:,:) = INT( zbathy(:,:) ) 1242 CALL lbc_lnk( risfdep, 'T', 1. ) 1243 CALL lbc_lnk( bathy, 'T', 1. ) 1244 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1245 CALL lbc_lnk( zbathy, 'T', 1. ) 1246 mbathy(:,:) = INT( zbathy(:,:) ) 1247 ENDIF 1248 END DO 1249 ! end dig bathy/ice shelf to be compatible 1250 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1251 DO jl = 1,20 1252 1253 ! remove single point "bay" on isf coast line in the ice shelf draft' 1254 DO jk = 1, jpk 1255 WHERE (misfdep==0) misfdep=jpk 1256 zmask=0 1257 WHERE (misfdep .LE. jk) zmask=1 1258 DO jj = 2, jpjm1 1259 DO ji = 2, jpim1 1260 IF (misfdep(ji,jj) .EQ. jk) THEN 1261 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1262 IF (ibtest .LE. 1) THEN 1263 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1264 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 1265 END IF 1266 END IF 1267 END DO 1268 END DO 1269 END DO 1270 WHERE (misfdep==jpk) 1271 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1272 END WHERE 1273 IF( lk_mpp ) THEN 1274 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1275 CALL lbc_lnk( zbathy, 'T', 1. ) 1276 misfdep(:,:) = INT( zbathy(:,:) ) 1277 CALL lbc_lnk( risfdep, 'T', 1. ) 1278 CALL lbc_lnk( bathy, 'T', 1. ) 1279 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1280 CALL lbc_lnk( zbathy, 'T', 1. ) 1281 mbathy(:,:) = INT( zbathy(:,:) ) 1282 ENDIF 1283 1284 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1285 DO jk = jpk,1,-1 1286 zmask=0 1287 WHERE (mbathy .GE. jk ) zmask=1 1288 DO jj = 2, jpjm1 1289 DO ji = 2, jpim1 1290 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 1291 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1292 IF (ibtest .LE. 1) THEN 1293 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1294 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 1295 END IF 1296 END IF 1297 END DO 1298 END DO 1299 END DO 1300 WHERE (mbathy==0) 1301 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1302 END WHERE 1303 IF( lk_mpp ) THEN 1304 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1305 CALL lbc_lnk( zbathy, 'T', 1. ) 1306 misfdep(:,:) = INT( zbathy(:,:) ) 1307 CALL lbc_lnk( risfdep, 'T', 1. ) 1308 CALL lbc_lnk( bathy, 'T', 1. ) 1309 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1310 CALL lbc_lnk( zbathy, 'T', 1. ) 1311 mbathy(:,:) = INT( zbathy(:,:) ) 1312 ENDIF 1313 1314 ! fill hole in ice shelf 1315 zmisfdep = misfdep 1316 zrisfdep = risfdep 1317 WHERE (zmisfdep .LE. 1) zmisfdep=jpk 1318 DO jj = 2, jpjm1 1319 DO ji = 2, jpim1 1320 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1321 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1322 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1) 1323 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1) 1324 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1) 1325 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1) 1326 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1327 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 1328 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1329 END IF 1330 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 1331 misfdep(ji,jj) = ibtest 1332 risfdep(ji,jj) = gdepw_1d(ibtest) 1333 ENDIF 1334 ENDDO 1335 ENDDO 1336 1337 IF( lk_mpp ) THEN 1338 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1339 CALL lbc_lnk( zbathy, 'T', 1. ) 1340 misfdep(:,:) = INT( zbathy(:,:) ) 1341 CALL lbc_lnk( risfdep, 'T', 1. ) 1342 CALL lbc_lnk( bathy, 'T', 1. ) 1343 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1344 CALL lbc_lnk( zbathy, 'T', 1. ) 1345 mbathy(:,:) = INT( zbathy(:,:) ) 1346 ENDIF 1347 ! 1348 !! fill hole in bathymetry 1349 zmbathy (:,:)=mbathy (:,:) 1350 DO jj = 2, jpjm1 1351 DO ji = 2, jpim1 1352 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1353 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1354 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1) 1355 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0 1356 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0 1357 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1358 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1359 IF( ibtest == 0 ) THEN 1360 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1361 END IF 1362 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 1363 mbathy(ji,jj) = ibtest 1364 bathy(ji,jj) = gdepw_1d(ibtest+1) 1365 ENDIF 1366 END DO 1367 END DO 1368 IF( lk_mpp ) THEN 1369 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1370 CALL lbc_lnk( zbathy, 'T', 1. ) 1371 misfdep(:,:) = INT( zbathy(:,:) ) 1372 CALL lbc_lnk( risfdep, 'T', 1. ) 1373 CALL lbc_lnk( bathy, 'T', 1. ) 1374 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1375 CALL lbc_lnk( zbathy, 'T', 1. ) 1376 mbathy(:,:) = INT( zbathy(:,:) ) 1377 ENDIF 1378 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1379 DO jj = 1, jpjm1 1380 DO ji = 1, jpim1 1381 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1382 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1383 END IF 1384 END DO 1385 END DO 1386 IF( lk_mpp ) THEN 1387 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1388 CALL lbc_lnk( zbathy, 'T', 1. ) 1389 misfdep(:,:) = INT( zbathy(:,:) ) 1390 CALL lbc_lnk( risfdep, 'T', 1. ) 1391 CALL lbc_lnk( bathy, 'T', 1. ) 1392 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1393 CALL lbc_lnk( zbathy, 'T', 1. ) 1394 mbathy(:,:) = INT( zbathy(:,:) ) 1395 ENDIF 1396 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1397 DO jj = 1, jpjm1 1398 DO ji = 1, jpim1 1399 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1400 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1401 END IF 1402 END DO 1403 END DO 1404 IF( lk_mpp ) THEN 1405 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1406 CALL lbc_lnk( zbathy, 'T', 1. ) 1407 misfdep(:,:) = INT( zbathy(:,:) ) 1408 CALL lbc_lnk( risfdep, 'T', 1. ) 1409 CALL lbc_lnk( bathy, 'T', 1. ) 1410 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1411 CALL lbc_lnk( zbathy, 'T', 1. ) 1412 mbathy(:,:) = INT( zbathy(:,:) ) 1413 ENDIF 1414 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1415 DO jj = 1, jpjm1 1416 DO ji = 1, jpi 1417 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1418 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1419 END IF 1420 END DO 1421 END DO 1422 IF( lk_mpp ) THEN 1423 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1424 CALL lbc_lnk( zbathy, 'T', 1. ) 1425 misfdep(:,:) = INT( zbathy(:,:) ) 1426 CALL lbc_lnk( risfdep, 'T', 1. ) 1427 CALL lbc_lnk( bathy, 'T', 1. ) 1428 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1429 CALL lbc_lnk( zbathy, 'T', 1. ) 1430 mbathy(:,:) = INT( zbathy(:,:) ) 1431 ENDIF 1432 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1433 DO jj = 1, jpjm1 1434 DO ji = 1, jpi 1435 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1436 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1437 END IF 1438 END DO 1439 END DO 1440 IF( lk_mpp ) THEN 1441 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1442 CALL lbc_lnk( zbathy, 'T', 1. ) 1443 misfdep(:,:) = INT( zbathy(:,:) ) 1444 CALL lbc_lnk( risfdep, 'T', 1. ) 1445 CALL lbc_lnk( bathy, 'T', 1. ) 1446 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1447 CALL lbc_lnk( zbathy, 'T', 1. ) 1448 mbathy(:,:) = INT( zbathy(:,:) ) 1449 ENDIF 1450 ! if not compatible after all check, mask T 1451 DO jj = 1, jpj 1452 DO ji = 1, jpi 1453 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1454 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1455 END IF 1456 END DO 1457 END DO 1458 1459 WHERE (mbathy(:,:) == 1) 1460 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1461 END WHERE 1462 END DO 1463 ! end check compatibility ice shelf/bathy 1464 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1465 WHERE (misfdep(:,:) <= 5) 1466 misfdep = 1; risfdep = 0.0_wp; 1467 END WHERE 1468 1469 IF( icompt == 0 ) THEN 1470 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1471 ELSE 1472 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1473 ENDIF 908 1474 909 1475 ! Scale factors and depth at T- and W-points … … 938 1504 !gm Bug? check the gdepw_1d 939 1505 ! ... on ik 940 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 941 & * ((gdept_1d(ik ) - gdepw_1d(ik) ) &942 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ))1506 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1507 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1508 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 943 1509 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 944 1510 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) … … 973 1539 END DO 974 1540 END DO 1541 ! 1542 ! (ISF) Definition of e3t, u, v, w for ISF case 1543 DO jj = 1, jpj 1544 DO ji = 1, jpi 1545 ik = misfdep(ji,jj) 1546 IF( ik > 1 ) THEN ! ice shelf point only 1547 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1548 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1549 !gm Bug? check the gdepw_0 1550 ! ... on ik 1551 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1552 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1553 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1554 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1555 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1556 1557 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1558 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1559 ENDIF 1560 ! ... on ik / ik-1 1561 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1562 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1563 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1564 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1565 ENDIF 1566 END DO 1567 END DO 1568 ! 1569 it = 0 1570 DO jj = 1, jpj 1571 DO ji = 1, jpi 1572 ik = misfdep(ji,jj) 1573 IF( ik > 1 ) THEN ! ice shelf point only 1574 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1575 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1576 ! test 1577 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1578 IF( zdiff <= 0. .AND. lwp ) THEN 1579 it = it + 1 1580 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1581 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1582 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1583 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1584 ENDIF 1585 ENDIF 1586 END DO 1587 END DO 1588 ! END (ISF) 975 1589 976 1590 ! Scale factors and depth at U-, V-, UW and VW-points … … 991 1605 END DO 992 1606 END DO 1607 ! (ISF) define e3uw 1608 DO jk = 2,jpk 1609 DO jj = 1, jpjm1 1610 DO ji = 1, fs_jpim1 ! vector opt. 1611 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1612 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1613 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1614 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1615 END DO 1616 END DO 1617 END DO 1618 !End (ISF) 1619 993 1620 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 994 1621 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1032 1659 1033 1660 ! Compute gdep3w_0 (vertical sum of e3w) 1034 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1035 DO jk = 2, jpk 1036 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1037 END DO 1038 1661 WHERE (misfdep == 0) misfdep = 1 1662 DO jj = 1,jpj 1663 DO ji = 1,jpi 1664 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1665 DO jk = 2, misfdep(ji,jj) 1666 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1667 END DO 1668 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)) 1669 DO jk = misfdep(ji,jj) + 1, jpk 1670 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1671 END DO 1672 END DO 1673 END DO 1039 1674 ! ! ================= ! 1040 1675 IF(lwp .AND. ll_print) THEN ! Control print ! … … 1066 1701 ! 1067 1702 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1703 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1704 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1068 1705 ! 1069 1706 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
Note: See TracChangeset
for help on using the changeset viewer.
