New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4747 for branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 – NEMO

Ignore:
Timestamp:
2014-08-15T16:32:52+02:00 (10 years ago)
Author:
mathiot
Message:

ISF branch: change to deal with non mask bathymetry (land processor definition, building bathy and ice shelf draft variable), update of hpg (definition of ze3wu in case of zps and vvl) and bfr (in case of 2 cell water column thickness, each cell feels top and bottom friction).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4726 r4747  
    102102      INTEGER ::   ios 
    103103      ! 
    104       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
    105       NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   & 
    106          &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr, nn_fwb, ln_cdgw, nn_isf 
     104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
    107105      !!---------------------------------------------------------------------- 
    108106      ! 
    … …  
    123121         WRITE(numout,*) '~~~~~~~' 
    124122         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate' 
    125          WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco 
    126          WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps 
    127          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 
    128127      ENDIF 
    129128 
    … …  
    533532            risfdep(:,:)=0._wp          
    534533            misfdep(:,:)=1              
    535             IF ( nn_isf == 1 .OR. nn_isf == 4 ) THEN 
     534            IF ( ln_isfcav ) THEN 
    536535               CALL iom_open ( 'isf_draft_meter.nc', inum )  
    537536               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep ) 
    … …  
    577576      !                        
    578577      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
     578         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
     579         WHERE (bathy == risfdep) 
     580            bathy   = 0.0_wp ; risfdep = 0.0_wp 
     581         END WHERE 
     582         ! end patch 
    579583         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    580584         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth 
    … …  
    956960      INTEGER  ::   ik, it           ! temporary integers 
    957961      INTEGER  ::   id, jd, nprocd 
    958       INTEGER  ::   icompt, ibtest   ! (ISF) 
     962      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
    959963      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    960964      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    … …  
    963967      REAL(wp) ::   zdiff            ! temporary scalar 
    964968      REAL(wp) ::   zrefdep          ! temporary scalar 
    965       REAL(wp) ::   eps=0.99            ! small offset to avoid large pool in case bathy slightly greater than risfdep 
    966       REAL(wp), POINTER, DIMENSION(:,:)   ::   zbathy, zmask   ! 3D workspace (ISH) 
     969      REAL(wp) ::   zbathydiff, zrisfdepdiff  
     970      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    967971      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    968972      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
    … …  
    972976      ! 
    973977      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
    974       CALL wrk_alloc( jpi, jpj, zbathy, zmask) 
     978      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    975979      CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    976980      ! 
    … …  
    988992      ENDIF 
    989993 
    990  
    991994      ! bathymetry in level (from bathy_meter) 
    992995      ! =================== 
    … …  
    10061009      END DO 
    10071010      ! (ISF) compute misfdep 
    1008       WHERE( risfdep(:,:) == 0._wp ) ;   misfdep(:,:) = 1   ! no ice shelf : set misfdep to 1   
    1009       ELSEWHERE                     ;   misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1011      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     1012      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    10101013      END WHERE   
    10111014 
    … …  
    10181021         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    10191022      END DO  
    1020       WHERE (risfdep <= e3t_1d(1) .AND. risfdep .GT. 0._wp) 
    1021          risfdep = 0. 
    1022          misfdep= 1 
     1023      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 
     1024         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    10231025      END WHERE 
    10241026  
    … …  
    10261028      icompt = 0  
    10271029! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    1028       DO jl = 1, 10  
     1030      DO jl = 1, 10      
     1031         WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
     1032            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     1033            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     1034         END WHERE 
     1035         WHERE (mbathy(:,:) <= 0)  
     1036            misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
     1037            mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
     1038         ENDWHERE 
    10291039         IF( lk_mpp ) THEN 
    10301040            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    … …  
    10411051            misfdep(jpi,:) = misfdep(  2  ,:)  
    10421052         ENDIF 
    1043  
     1053  
    10441054         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10451055            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10461056            mbathy(jpi,:) = mbathy(  2  ,:) 
    10471057         ENDIF 
    1048  
    1049          WHERE (mbathy == 0)  
    1050             risfdep = 0._wp 
    1051             misfdep= 0 
    1052             bathy  = 0._wp 
    1053          ENDWHERE 
    1054  
    1055          WHERE (bathy(:,:) < risfdep(:,:)+eps) 
    1056             misfdep(:,:) = 0  
    1057             risfdep(:,:)  = 0._wp 
    1058             mbathy(:,:)  = 0 
    1059             bathy(:,:)   = 0._wp 
    1060          END WHERE 
    1061 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     1058  
     1059         ! split last cell if possible (only where water column is 2 cell or less) 
     1060         DO jk = jpkm1, 1, -1 
     1061            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1062            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
     1063               mbathy(:,:) = jk 
     1064               bathy(:,:)  = zmax 
     1065            END WHERE 
     1066         END DO 
     1067  
     1068         ! split top cell if possible (only where water column is 2 cell or less) 
     1069         DO jk = 2, jpkm1 
     1070            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1071            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
     1072               misfdep(:,:) = jk 
     1073               risfdep(:,:) = zmax 
     1074            END WHERE 
     1075         END DO 
     1076  
     1077  
     1078 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    10621079         DO jj = 1, jpj 
    10631080            DO ji = 1, jpi 
    1064                IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1065                   bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1066                   mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1081               ! find the minimum change option: 
     1082               ! test bathy 
     1083               IF (risfdep(ji,jj) .GT. 1) THEN 
     1084               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1085               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1086  
     1087                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
     1088                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1089                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1090                        mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1091                     ELSE 
     1092                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
     1093                        misfdep(ji,jj) = misfdep(ji,jj) - 1 
     1094                     END IF 
     1095                  END IF 
    10671096               END IF 
    10681097            END DO 
    10691098         END DO 
    1070  
    1071  
    1072         ! At least 2 levels for water thickness at T, U, and V point. 
    1073          zmisfdep(:,:)=misfdep(:,:) 
    1074          zmbathy (:,:)=mbathy (:,:) 
    10751099  
    1076          DO jj = 2, jpjm1 
    1077             DO ji = 2, jpim1 
    1078                ! T point 
    1079                IF( zmisfdep(ji,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
    1080                   mbathy(ji,jj) = zmbathy(ji,jj) + 1  
    1081                   bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1100          ! At least 2 levels for water thickness at T, U, and V point. 
     1101         DO jj = 1, jpj 
     1102            DO ji = 1, jpi 
     1103               ! find the minimum change option: 
     1104               ! test bathy 
     1105               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1106                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1107                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1108                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1109                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1110                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1111                  ELSE 
     1112                     misfdep(ji,jj)= misfdep(ji,jj) - 1 
     1113                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1114                  END IF 
    10821115               ENDIF 
    1083                ! V point 
    1084                IF( zmisfdep(ji,jj+1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
    1085                   mbathy(ji,jj) = zmbathy(ji,jj) + 1  
    1086                   bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1116            END DO 
     1117         END DO 
     1118  
     1119 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1120         DO jj = 1, jpjm1 
     1121            DO ji = 1, jpim1 
     1122               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1123                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
     1124                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
     1125                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1126                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1127                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
     1128                  ELSE 
     1129                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
     1130                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1131                  END IF 
    10871132               ENDIF 
    1088                ! V point -1 
    1089                IF( zmisfdep(ji,jj-1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
    1090                   mbathy(ji,jj) = zmbathy(ji,jj) + 1  
    1091                   bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1133            END DO 
     1134         END DO 
     1135  
     1136         IF( lk_mpp ) THEN 
     1137            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1138            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1139            misfdep(:,:) = INT( zbathy(:,:) ) 
     1140            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1141            CALL lbc_lnk( bathy, 'T', 1. ) 
     1142            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1143            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1144            mbathy(:,:) = INT( zbathy(:,:) ) 
     1145         ENDIF 
     1146 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
     1147         DO jj = 1, jpjm1 
     1148            DO ji = 1, jpim1 
     1149               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1150                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
     1151                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
     1152                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1153                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
     1154                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     1155                  ELSE 
     1156                     misfdep(ji,jj)   = misfdep(ji,jj) - 1 
     1157                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
     1158                  END IF 
    10921159               ENDIF 
    1093                ! U point 
    1094                IF( zmisfdep(ji+1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
    1095                   mbathy(ji,jj) = zmbathy(ji,jj) + 1  
    1096                   bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
    1097                ENDIF 
    1098                ! U point -1 
    1099                IF( zmisfdep(ji-1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
    1100                   mbathy(ji,jj) = zmbathy(ji,jj) + 1  
    1101                   bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
    1102                ENDIF 
    1103             END DO 
    1104          END DO 
     1160            END DO 
     1161         END DO 
     1162  
     1163  
    11051164         IF( lk_mpp ) THEN  
    11061165            zbathy(:,:) = FLOAT( misfdep(:,:) )  
    … …  
    11131172            mbathy(:,:) = INT( zbathy(:,:) ) 
    11141173         ENDIF  
    1115  
    1116         ! if single ocean point put as land 
    1117          zmask=1 
    1118          WHERE (mbathy .EQ. 0) zmask=0 
    1119          DO jj = 2, jpjm1 
    1120             DO ji = 2, jpim1 
    1121                ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1122                IF (ibtest .LE. 1) THEN 
    1123                   bathy(ji,jj)=0._wp 
    1124                   mbathy(ji,jj)=0 
    1125                   risfdep(ji,jj)=0._wp 
    1126                   misfdep(ji,jj)=0 
    1127                END IF 
    1128             END DO 
    1129          END DO 
     1174  
     1175 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1176         DO jj = 1, jpjm1 
     1177            DO ji = 1, jpim1 
     1178               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1179                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
     1180                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
     1181                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1182                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1183                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1184                  ELSE 
     1185                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
     1186                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
     1187                  END IF 
     1188               ENDIF 
     1189            ENDDO 
     1190         ENDDO 
     1191  
    11301192         IF( lk_mpp ) THEN  
    11311193            zbathy(:,:) = FLOAT( misfdep(:,:) )  
    … …  
    11371199            CALL lbc_lnk( zbathy, 'T', 1. ) 
    11381200            mbathy(:,:) = INT( zbathy(:,:) ) 
    1139          END IF  
    1140  
    1141          ! if single point on isf coast line 
     1201         ENDIF  
     1202  
     1203 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
     1204         DO jj = 1, jpjm1 
     1205            DO ji = 1, jpim1 
     1206               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1207                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
     1208                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
     1209                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1210                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
     1211                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
     1212                  ELSE 
     1213                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
     1214                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1215                  END IF 
     1216               ENDIF 
     1217            ENDDO 
     1218         ENDDO 
     1219  
     1220         IF( lk_mpp ) THEN 
     1221            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1222            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1223            misfdep(:,:) = INT( zbathy(:,:) ) 
     1224            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1225            CALL lbc_lnk( bathy, 'T', 1. ) 
     1226            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1227            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1228            mbathy(:,:) = INT( zbathy(:,:) ) 
     1229         ENDIF 
     1230      END DO 
     1231      ! end dig bathy/ice shelf to be compatible 
     1232      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
     1233      DO jl = 1,20 
     1234  
     1235 ! remove single point "bay" on isf coast line in the ice shelf draft' 
    11421236         DO jk = 1, jpk 
    11431237            WHERE (misfdep==0) misfdep=jpk 
    … …  
    11551249               END DO 
    11561250            END DO 
    1157             WHERE (misfdep==jpk)  
    1158                misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    1159             END WHERE 
    1160    
    1161             IF( lk_mpp ) THEN  
    1162                zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1163                CALL lbc_lnk( zbathy, 'T', 1. )  
    1164                misfdep(:,:) = INT( zbathy(:,:) )  
    1165                CALL lbc_lnk( risfdep, 'T', 1. )  
    1166                CALL lbc_lnk( bathy, 'T', 1. ) 
    1167                zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1168                CALL lbc_lnk( zbathy, 'T', 1. ) 
    1169                mbathy(:,:) = INT( zbathy(:,:) ) 
    1170             ENDIF  
    1171          END DO 
    1172  
    1173 ! fill hole in ice shelf 
    1174          WHERE (misfdep==0) misfdep=jpk 
    1175          zmisfdep(:,:)=misfdep(:,:) 
    1176          zmbathy (:,:)=mbathy (:,:) 
     1251         END DO 
     1252         WHERE (misfdep==jpk) 
     1253             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1254         END WHERE 
     1255         IF( lk_mpp ) THEN 
     1256            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1257            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1258            misfdep(:,:) = INT( zbathy(:,:) ) 
     1259            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1260            CALL lbc_lnk( bathy, 'T', 1. ) 
     1261            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1262            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1263            mbathy(:,:) = INT( zbathy(:,:) ) 
     1264         ENDIF 
     1265  
     1266 ! remove single point "bay" on bathy coast line beneath an ice shelf' 
     1267         DO jk = jpk,1,-1 
     1268            zmask=0 
     1269            WHERE (mbathy .GE. jk ) zmask=1 
     1270            DO jj = 2, jpjm1 
     1271               DO ji = 2, jpim1 
     1272                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1273                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1274                     IF (ibtest .LE. 1) THEN 
     1275                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
     1276                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 
     1277                     END IF 
     1278                  END IF 
     1279               END DO 
     1280            END DO 
     1281         END DO 
     1282         WHERE (mbathy==0) 
     1283             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1284         END WHERE 
     1285         IF( lk_mpp ) THEN 
     1286            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1287            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1288            misfdep(:,:) = INT( zbathy(:,:) ) 
     1289            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1290            CALL lbc_lnk( bathy, 'T', 1. ) 
     1291            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1292            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1293            mbathy(:,:) = INT( zbathy(:,:) ) 
     1294         ENDIF 
     1295  
     1296 ! fill hole in ice shelf 
     1297         zmisfdep = misfdep 
     1298         zrisfdep = risfdep 
     1299         WHERE (zmisfdep .LE. 1) zmisfdep=jpk 
    11771300         DO jj = 2, jpjm1 
    11781301            DO ji = 2, jpim1 
    1179                ibtest=MIN(zmisfdep(ji-1,jj), zmisfdep(ji+1,jj), zmisfdep(ji,jj-1), zmisfdep(ji,jj+1)) 
    1180                IF( ibtest > zmisfdep(ji,jj)) THEN 
     1302               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
     1303               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
     1304               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
     1305               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
     1306               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
     1307               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
     1308               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1309               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1310                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
     1311               END IF 
     1312               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 
    11811313                  misfdep(ji,jj) = ibtest 
    1182                   risfdep(ji,jj)  = gdepw_1d(ibtest) 
     1314                  risfdep(ji,jj) = gdepw_1d(ibtest) 
    11831315               ENDIF 
    11841316            ENDDO 
    11851317         ENDDO 
    1186          WHERE (misfdep==jpk)  
    1187             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0 
    1188          END WHERE 
     1318  
    11891319         IF( lk_mpp ) THEN  
    11901320            zbathy(:,:) = FLOAT( misfdep(:,:) )  
    … …  
    11971327            mbathy(:,:) = INT( zbathy(:,:) ) 
    11981328         ENDIF  
    1199  
    1200 ! fill hole in bathymetry 
    1201          zmisfdep(:,:)=misfdep(:,:) 
     1329 ! 
     1330 !! fill hole in bathymetry 
    12021331         zmbathy (:,:)=mbathy (:,:) 
    12031332         DO jj = 2, jpjm1 
    12041333            DO ji = 2, jpim1 
    1205                ibtest = MAX(  zmbathy(ji-1,jj), zmbathy(ji+1,jj),   & 
    1206                   &           zmbathy(ji,jj-1), zmbathy(ji,jj+1)  ) 
    1207                IF( ibtest < zmbathy(ji,jj) ) THEN 
     1334               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
     1335               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
     1336               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
     1337               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
     1338               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
     1339               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
     1340               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1341               IF( ibtest == 0 ) THEN 
     1342                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
     1343               END IF 
     1344               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 
    12081345                  mbathy(ji,jj) = ibtest 
    12091346                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
    … …  
    12211358            mbathy(:,:) = INT( zbathy(:,:) ) 
    12221359         ENDIF  
    1223         ! remove 1 cell pool of water stuck between ice shelf and bathymetry (need a 3D flood filling tools to do this properly) 
    1224          DO jk = 1, jpk 
    1225             WHERE (misfdep==0) misfdep=jpk 
    1226             zmisfdep(:,:)=misfdep(:,:) 
    1227             zmbathy (:,:)=mbathy (:,:) 
    1228             DO jj = 2, jpjm1 
    1229                DO ji = 2, jpim1 
    1230                   IF ( jk .GE. zmisfdep(ji,jj) .AND. jk .LE. zmbathy(ji,jj) ) THEN 
    1231                      IF ( (jk > zmbathy(ji,jj+1) .OR. jk < zmisfdep(ji,jj+1)) .AND. & 
    1232                         & (jk > zmbathy(ji,jj-1) .OR. jk < zmisfdep(ji,jj-1)) .AND. & 
    1233                         & (jk > zmbathy(ji+1,jj) .OR. jk < zmisfdep(ji+1,jj)) .AND. & 
    1234                         & (jk > zmbathy(ji-1,jj) .OR. jk < zmisfdep(ji-1,jj))       ) THEN 
    1235                         mbathy(ji,jj) = 0 ; misfdep(ji,jj) = jpk ; risfdep(ji,jj) = 0._wp ; bathy(ji,jj) = 0._wp 
    1236                      ENDIF 
    1237                   ENDIF 
    1238                ENDDO 
    1239             ENDDO 
    1240             WHERE (misfdep==jpk) misfdep=0  
    1241             IF( lk_mpp ) THEN  
    1242                zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1243                CALL lbc_lnk( zbathy, 'T', 1. )  
    1244                misfdep(:,:) = INT( zbathy(:,:) )  
    1245                CALL lbc_lnk( risfdep, 'T', 1. )  
    1246                CALL lbc_lnk( bathy, 'T', 1. ) 
    1247                zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1248                CALL lbc_lnk( zbathy, 'T', 1. ) 
    1249                mbathy(:,:) = INT( zbathy(:,:) ) 
    1250             ENDIF  
    1251          END DO 
     1360 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
     1361         DO jj = 1, jpjm1 
     1362            DO ji = 1, jpim1 
     1363               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1364                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
     1365               END IF 
     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 (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1382                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,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 V point water column less than 2 cells), mask V 
     1397         DO jj = 1, jpjm1 
     1398            DO ji = 1, jpi 
     1399               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
     1400                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,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 (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
     1418                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+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, mask T 
     1433         DO jj = 1, jpj 
     1434            DO ji = 1, jpi 
     1435               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
     1436                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
     1437               END IF 
     1438            END DO 
     1439         END DO 
     1440  
     1441         WHERE (mbathy(:,:) == 1) 
     1442            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
     1443         END WHERE 
    12521444      END DO  
    1253  
    1254       WHERE (mbathy(:,:) < misfdep(:,:)) 
    1255          misfdep(:,:) = 0 
    1256          risfdep(:,:)  = 0._wp 
    1257          mbathy(:,:)  = 0 
    1258          bathy(:,:)   = 0._wp 
     1445! end check compatibility ice shelf/bathy 
     1446      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1447      WHERE (misfdep(:,:) <= 5) 
     1448         misfdep = 1; risfdep = 0.0_wp; 
    12591449      END WHERE 
    1260  
    12611450 
    12621451      IF( icompt == 0 ) THEN  
    … …  
    13381527            ik = misfdep(ji,jj)  
    13391528            IF( ik > 1 ) THEN               ! ice shelf point only  
    1340                 IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1341                 gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1529               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1530               gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    13421531!gm Bug?  check the gdepw_0  
    1343                 !       ... on ik  
    1344                 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1345                    &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1346                    &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1347                 e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1348                 e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1349  
    1350                 IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1351                    e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1532               !       ... on ik  
     1533               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1534                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1535                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1536               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1537               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1538 
     1539               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1540                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    13521541                ENDIF  
    1353                 !       ... on ik / ik-1  
    1354                 e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1355                 e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1542               !       ... on ik / ik-1  
     1543               e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1544               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    13561545! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1357                 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1546               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    13581547            ENDIF  
    13591548         END DO  
Note: See TracChangeset for help on using the changeset viewer.