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 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 – NEMO

Ignore:
Timestamp:
2015-12-01T16:35:30+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4624 r5965  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            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   
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    101102      INTEGER ::   ios 
    102103      ! 
    103       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
     104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
    104105      !!---------------------------------------------------------------------- 
    105106      ! 
     
    120121         WRITE(numout,*) '~~~~~~~' 
    121122         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 
    125127      ENDIF 
    126128 
     
    145147      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    146148                          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 
    147150      ! 
    148151      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain 
     
    295298      ENDIF 
    296299 
     300      IF ( ln_isfcav ) THEN 
     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) 
     302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
     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 
     313 
    297314!!gm BUG in s-coordinate this does not work! 
    298315      ! deepest/shallowest W level Above/Below ~10m 
     
    350367      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
    351368      INTEGER  ::   inum                      ! temporary logical unit 
     369      INTEGER  ::   ierror                    ! error flag 
    352370      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
    353371      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices 
    354372      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
    355373      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars 
    356       INTEGER , POINTER, DIMENSION(:,:) ::   idta   ! global domain integer data 
    357       REAL(wp), POINTER, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
     374      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data 
     375      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
    358376      !!---------------------------------------------------------------------- 
    359377      ! 
    360378      IF( nn_timing == 1 )  CALL timing_start('zgr_bat') 
    361       ! 
    362       CALL wrk_alloc( jpidta, jpjdta, idta ) 
    363       CALL wrk_alloc( jpidta, jpjdta, zdta ) 
    364379      ! 
    365380      IF(lwp) WRITE(numout,*) 
     
    370385         !                                            ! ================== ! 
    371386         !                                            ! global domain level and meter bathymetry (idta,zdta) 
     387         ! 
     388         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 
     389         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 
     390         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 
     391         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 
    372392         ! 
    373393         IF( ntopo == 0 ) THEN                        ! flat basin 
     
    450470            END DO 
    451471         END DO 
     472         risfdep(:,:)=0.e0 
     473         misfdep(:,:)=1 
     474         ! 
     475         DEALLOCATE( idta, zdta ) 
    452476         ! 
    453477         !                                            ! ================ ! 
     
    490514         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
    491515            CALL iom_open ( 'bathy_meter.nc', inum )  
    492             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     516            IF ( ln_isfcav ) THEN 
     517               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 
     518            ELSE 
     519               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  ) 
     520            END IF 
    493521            CALL iom_close( inum ) 
    494522            !                                                 
     523            risfdep(:,:)=0._wp          
     524            misfdep(:,:)=1              
     525            IF ( ln_isfcav ) THEN 
     526               CALL iom_open ( 'isf_draft_meter.nc', inum )  
     527               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep ) 
     528               CALL iom_close( inum ) 
     529               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
     530            END IF 
     531            !        
    495532            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    496533               ! 
     
    530567      !                        
    531568      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
     569         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
     570         IF ( ln_isfcav ) THEN 
     571            WHERE (bathy == risfdep) 
     572               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     573            END WHERE 
     574         END IF 
     575         ! end patch 
    532576         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    533577         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth 
     
    539583         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 
    540584      ENDIF 
    541       ! 
    542       CALL wrk_dealloc( jpidta, jpjdta, idta ) 
    543       CALL wrk_dealloc( jpidta, jpjdta, zdta ) 
    544585      ! 
    545586      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat') 
     
    783824   END SUBROUTINE zgr_bot_level 
    784825 
     826      SUBROUTINE zgr_top_level 
     827      !!---------------------------------------------------------------------- 
     828      !!                    ***  ROUTINE zgr_bot_level  *** 
     829      !! 
     830      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays) 
     831      !! 
     832      !! ** Method  :   computes from misfdep with a minimum value of 1 
     833      !! 
     834      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest  
     835      !!                                     ocean level at t-, u- & v-points 
     836      !!                                     (min value = 1) 
     837      !!---------------------------------------------------------------------- 
     838      !! 
     839      INTEGER ::   ji, jj   ! dummy loop indices 
     840      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik 
     841      !!---------------------------------------------------------------------- 
     842      ! 
     843      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level') 
     844      ! 
     845      CALL wrk_alloc( jpi, jpj, zmik ) 
     846      ! 
     847      IF(lwp) WRITE(numout,*) 
     848      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 
     849      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     850      ! 
     851      mikt(:,:) = MAX( misfdep(:,:) , 1 )    ! top k-index of T-level (=1) 
     852      !                                      ! top k-index of W-level (=mikt) 
     853      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level 
     854         DO ji = 1, jpim1 
     855            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     856            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     857            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     858         END DO 
     859      END DO 
     860 
     861      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     862      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     863      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     864      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     865      ! 
     866      CALL wrk_dealloc( jpi, jpj, zmik ) 
     867      ! 
     868      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level') 
     869      ! 
     870   END SUBROUTINE zgr_top_level 
    785871 
    786872   SUBROUTINE zgr_zco 
     
    862948      !! 
    863949      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    864       INTEGER  ::   ik, it          ! temporary integers 
     950      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
    865951      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    866952      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     
    906992         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    907993      END DO 
     994 
     995      IF ( ln_isfcav ) CALL zgr_isf 
    908996 
    909997      ! Scale factors and depth at T- and W-points 
     
    9381026!gm Bug?  check the gdepw_1d 
    9391027                  !       ... on ik 
    940                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0   (ji,jj,ik+1) - gdepw_1d(ik) )   & 
    941                      &                           * ((gdept_1d(      ik  ) - gdepw_1d(ik) )   & 
    942                      &                           / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )) 
    943                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    944                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1028                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1029                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1030                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1031                  e3t_0  (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1032                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    9451033                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    9461034                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     
    9731061         END DO 
    9741062      END DO 
     1063      ! 
     1064      IF ( ln_isfcav ) THEN 
     1065      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1066         DO jj = 1, jpj  
     1067            DO ji = 1, jpi  
     1068               ik = misfdep(ji,jj)  
     1069               IF( ik > 1 ) THEN               ! ice shelf point only  
     1070                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1071                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1072!gm Bug?  check the gdepw_0  
     1073               !       ... on ik  
     1074                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1075                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1076                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1077                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1078                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1079 
     1080                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1081                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1082                  ENDIF  
     1083               !       ... on ik / ik-1  
     1084                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1085                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1086! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1087                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1088               ENDIF  
     1089            END DO  
     1090         END DO  
     1091      !  
     1092         it = 0  
     1093         DO jj = 1, jpj  
     1094            DO ji = 1, jpi  
     1095               ik = misfdep(ji,jj)  
     1096               IF( ik > 1 ) THEN               ! ice shelf point only  
     1097                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1098                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1099               ! test  
     1100                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1101                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1102                     it = it + 1  
     1103                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1104                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1105                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1106                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1107                  ENDIF  
     1108               ENDIF  
     1109            END DO  
     1110         END DO  
     1111      END IF 
     1112      ! END (ISF) 
    9751113 
    9761114      ! Scale factors and depth at U-, V-, UW and VW-points 
     
    9911129         END DO 
    9921130      END DO 
     1131      IF ( ln_isfcav ) THEN 
     1132      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1133         DO jj = 2, jpjm1  
     1134            DO ji = 2, fs_jpim1   ! vector opt.  
     1135               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1136               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1137               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1138                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1139               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1140               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1141               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1142                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1143            END DO 
     1144         END DO 
     1145      END IF 
     1146 
    9931147      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    9941148      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     
    10321186      
    10331187      ! 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          
     1188      IF ( ln_isfcav ) THEN ! if cavity 
     1189         WHERE (misfdep == 0) misfdep = 1 
     1190         DO jj = 1,jpj 
     1191            DO ji = 1,jpi 
     1192               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1193               DO jk = 2, misfdep(ji,jj) 
     1194                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1195               END DO 
     1196               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)) 
     1197               DO jk = misfdep(ji,jj) + 1, jpk 
     1198                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1199               END DO 
     1200            END DO 
     1201         END DO 
     1202      ELSE ! no cavity 
     1203         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1204         DO jk = 2, jpk 
     1205            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1206         END DO 
     1207      END IF 
    10391208      !                                               ! ================= ! 
    10401209      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     
    10701239      ! 
    10711240   END SUBROUTINE zgr_zps 
     1241 
     1242   SUBROUTINE zgr_isf 
     1243      !!---------------------------------------------------------------------- 
     1244      !!                    ***  ROUTINE zgr_isf  *** 
     1245      !!    
     1246      !! ** Purpose :   check the bathymetry in levels 
     1247      !!    
     1248      !! ** Method  :   THe water column have to contained at least 2 cells 
     1249      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1250      !!                this criterion. 
     1251      !!                  
     1252      !!    
     1253      !! ** Action  : - test compatibility between isfdraft and bathy  
     1254      !!              - bathy and isfdraft are modified 
     1255      !!---------------------------------------------------------------------- 
     1256      !!    
     1257      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     1258      INTEGER  ::   ik, it           ! temporary integers 
     1259      INTEGER  ::   id, jd, nprocd 
     1260      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
     1261      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     1262      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     1263      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     1264      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
     1265      REAL(wp) ::   zdiff            ! temporary scalar 
     1266      REAL(wp) ::   zrefdep          ! temporary scalar 
     1267      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1268      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1269      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
     1270      !!--------------------------------------------------------------------- 
     1271      ! 
     1272      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1273      ! 
     1274      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
     1275      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1276 
     1277 
     1278      ! (ISF) compute misfdep 
     1279      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     1280      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1281      END WHERE   
     1282 
     1283      ! Compute misfdep for ocean points (i.e. first wet level)  
     1284      ! find the first ocean level such that the first level thickness  
     1285      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
     1286      ! e3t_0 is the reference level thickness  
     1287      DO jk = 2, jpkm1  
     1288         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
     1289         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
     1290      END DO  
     1291      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 
     1292         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
     1293      END WHERE 
     1294  
     1295! 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 
     1296      icompt = 0  
     1297! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
     1298      DO jl = 1, 10      
     1299         WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
     1300            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     1301            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     1302         END WHERE 
     1303         WHERE (mbathy(:,:) <= 0)  
     1304            misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
     1305            mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
     1306         ENDWHERE 
     1307         IF( lk_mpp ) THEN 
     1308            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1309            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1310            misfdep(:,:) = INT( zbathy(:,:) ) 
     1311            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1312            CALL lbc_lnk( bathy, 'T', 1. ) 
     1313            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1314            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1315            mbathy(:,:) = INT( zbathy(:,:) ) 
     1316         ENDIF 
     1317         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
     1318            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west  
     1319            misfdep(jpi,:) = misfdep(  2  ,:)  
     1320         ENDIF 
     1321 
     1322         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
     1323            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
     1324            mbathy(jpi,:) = mbathy(  2  ,:) 
     1325         ENDIF 
     1326 
     1327         ! split last cell if possible (only where water column is 2 cell or less) 
     1328         DO jk = jpkm1, 1, -1 
     1329            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1330            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
     1331               mbathy(:,:) = jk 
     1332               bathy(:,:)  = zmax 
     1333            END WHERE 
     1334         END DO 
     1335  
     1336         ! split top cell if possible (only where water column is 2 cell or less) 
     1337         DO jk = 2, jpkm1 
     1338            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1339            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
     1340               misfdep(:,:) = jk 
     1341               risfdep(:,:) = zmax 
     1342            END WHERE 
     1343         END DO 
     1344 
     1345  
     1346 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     1347         DO jj = 1, jpj 
     1348            DO ji = 1, jpi 
     1349               ! find the minimum change option: 
     1350               ! test bathy 
     1351               IF (risfdep(ji,jj) .GT. 1) THEN 
     1352               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
     1353                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1354               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
     1355                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1356  
     1357                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
     1358                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1359                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1360                        mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1361                     ELSE 
     1362                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
     1363                        misfdep(ji,jj) = misfdep(ji,jj) - 1 
     1364                     END IF 
     1365                  END IF 
     1366               END IF 
     1367            END DO 
     1368         END DO 
     1369  
     1370          ! At least 2 levels for water thickness at T, U, and V point. 
     1371         DO jj = 1, jpj 
     1372            DO ji = 1, jpi 
     1373               ! find the minimum change option: 
     1374               ! test bathy 
     1375               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1376                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)& 
     1377                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1378                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  & 
     1379                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1380                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1381                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1382                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1383                  ELSE 
     1384                     misfdep(ji,jj)= misfdep(ji,jj) - 1 
     1385                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1386                  END IF 
     1387               ENDIF 
     1388            END DO 
     1389         END DO 
     1390  
     1391 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1392         DO jj = 1, jpjm1 
     1393            DO ji = 1, jpim1 
     1394               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1395                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   & 
     1396                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
     1397                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  & 
     1398                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
     1399                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1400                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1401                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) & 
     1402                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
     1403                  ELSE 
     1404                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
     1405                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 
     1406                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1407                  END IF 
     1408               ENDIF 
     1409            END DO 
     1410         END DO 
     1411  
     1412         IF( lk_mpp ) THEN 
     1413            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1414            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1415            misfdep(:,:) = INT( zbathy(:,:) ) 
     1416            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1417            CALL lbc_lnk( bathy, 'T', 1. ) 
     1418            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1419            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1420            mbathy(:,:) = INT( zbathy(:,:) ) 
     1421         ENDIF 
     1422 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
     1423         DO jj = 1, jpjm1 
     1424            DO ji = 1, jpim1 
     1425               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1426                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 
     1427                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
     1428                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  & 
     1429                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
     1430                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1431                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
     1432                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     1433                  ELSE 
     1434                     misfdep(ji,jj)   = misfdep(ji,jj) - 1 
     1435                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
     1436                  END IF 
     1437               ENDIF 
     1438            END DO 
     1439         END DO 
     1440  
     1441  
     1442         IF( lk_mpp ) THEN  
     1443            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1444            CALL lbc_lnk( zbathy, 'T', 1. )  
     1445            misfdep(:,:) = INT( zbathy(:,:) )  
     1446            CALL lbc_lnk( risfdep, 'T', 1. )  
     1447            CALL lbc_lnk( bathy, 'T', 1. ) 
     1448            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1449            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1450            mbathy(:,:) = INT( zbathy(:,:) ) 
     1451         ENDIF  
     1452  
     1453 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1454         DO jj = 1, jpjm1 
     1455            DO ji = 1, jpim1 
     1456               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1457                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 
     1458                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
     1459                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 
     1460                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
     1461                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1462                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1463                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1464                  ELSE 
     1465                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
     1466                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
     1467                  END IF 
     1468               ENDIF 
     1469            ENDDO 
     1470         ENDDO 
     1471  
     1472         IF( lk_mpp ) THEN  
     1473            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1474            CALL lbc_lnk( zbathy, 'T', 1. )  
     1475            misfdep(:,:) = INT( zbathy(:,:) )  
     1476            CALL lbc_lnk( risfdep, 'T', 1. )  
     1477            CALL lbc_lnk( bathy, 'T', 1. ) 
     1478            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1479            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1480            mbathy(:,:) = INT( zbathy(:,:) ) 
     1481         ENDIF  
     1482  
     1483 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
     1484         DO jj = 1, jpjm1 
     1485            DO ji = 1, jpim1 
     1486               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1487                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 
     1488                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
     1489                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  & 
     1490                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
     1491                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1492                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
     1493                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  & 
     1494                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
     1495                  ELSE 
     1496                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
     1497                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) & 
     1498                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1499                  END IF 
     1500               ENDIF 
     1501            ENDDO 
     1502         ENDDO 
     1503  
     1504         IF( lk_mpp ) THEN 
     1505            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1506            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1507            misfdep(:,:) = INT( zbathy(:,:) ) 
     1508            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1509            CALL lbc_lnk( bathy, 'T', 1. ) 
     1510            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1511            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1512            mbathy(:,:) = INT( zbathy(:,:) ) 
     1513         ENDIF 
     1514      END DO 
     1515      ! end dig bathy/ice shelf to be compatible 
     1516      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
     1517      DO jl = 1,20 
     1518  
     1519 ! remove single point "bay" on isf coast line in the ice shelf draft' 
     1520         DO jk = 2, jpk 
     1521            WHERE (misfdep==0) misfdep=jpk 
     1522            zmask=0 
     1523            WHERE (misfdep .LE. jk) zmask=1 
     1524            DO jj = 2, jpjm1 
     1525               DO ji = 2, jpim1 
     1526                  IF (misfdep(ji,jj) .EQ. jk) THEN 
     1527                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1528                     IF (ibtest .LE. 1) THEN 
     1529                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
     1530                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 
     1531                     END IF 
     1532                  END IF 
     1533               END DO 
     1534            END DO 
     1535         END DO 
     1536         WHERE (misfdep==jpk) 
     1537             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1538         END WHERE 
     1539         IF( lk_mpp ) THEN 
     1540            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1541            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1542            misfdep(:,:) = INT( zbathy(:,:) ) 
     1543            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1544            CALL lbc_lnk( bathy, 'T', 1. ) 
     1545            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1546            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1547            mbathy(:,:) = INT( zbathy(:,:) ) 
     1548         ENDIF 
     1549  
     1550 ! remove single point "bay" on bathy coast line beneath an ice shelf' 
     1551         DO jk = jpk,1,-1 
     1552            zmask=0 
     1553            WHERE (mbathy .GE. jk ) zmask=1 
     1554            DO jj = 2, jpjm1 
     1555               DO ji = 2, jpim1 
     1556                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1557                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1558                     IF (ibtest .LE. 1) THEN 
     1559                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
     1560                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 
     1561                     END IF 
     1562                  END IF 
     1563               END DO 
     1564            END DO 
     1565         END DO 
     1566         WHERE (mbathy==0) 
     1567             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1568         END WHERE 
     1569         IF( lk_mpp ) THEN 
     1570            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1571            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1572            misfdep(:,:) = INT( zbathy(:,:) ) 
     1573            CALL lbc_lnk( risfdep, 'T', 1. ) 
     1574            CALL lbc_lnk( bathy, 'T', 1. ) 
     1575            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1576            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1577            mbathy(:,:) = INT( zbathy(:,:) ) 
     1578         ENDIF 
     1579  
     1580 ! fill hole in ice shelf 
     1581         zmisfdep = misfdep 
     1582         zrisfdep = risfdep 
     1583         WHERE (zmisfdep .LE. 1) zmisfdep=jpk 
     1584         DO jj = 2, jpjm1 
     1585            DO ji = 2, jpim1 
     1586               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
     1587               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
     1588               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
     1589               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
     1590               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
     1591               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
     1592               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1593               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1594                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
     1595               END IF 
     1596               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 
     1597                  misfdep(ji,jj) = ibtest 
     1598                  risfdep(ji,jj) = gdepw_1d(ibtest) 
     1599               ENDIF 
     1600            ENDDO 
     1601         ENDDO 
     1602  
     1603         IF( lk_mpp ) THEN  
     1604            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1605            CALL lbc_lnk( zbathy, 'T', 1. )  
     1606            misfdep(:,:) = INT( zbathy(:,:) )  
     1607            CALL lbc_lnk( risfdep, 'T', 1. )  
     1608            CALL lbc_lnk( bathy, 'T', 1. ) 
     1609            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1610            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1611            mbathy(:,:) = INT( zbathy(:,:) ) 
     1612         ENDIF  
     1613 ! 
     1614 !! fill hole in bathymetry 
     1615         zmbathy (:,:)=mbathy (:,:) 
     1616         DO jj = 2, jpjm1 
     1617            DO ji = 2, jpim1 
     1618               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
     1619               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
     1620               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
     1621               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
     1622               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
     1623               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
     1624               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1625               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
     1626                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
     1627               END IF 
     1628               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 
     1629                  mbathy(ji,jj) = ibtest 
     1630                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
     1631               ENDIF 
     1632            END DO 
     1633         END DO 
     1634         IF( lk_mpp ) THEN  
     1635            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1636            CALL lbc_lnk( zbathy, 'T', 1. )  
     1637            misfdep(:,:) = INT( zbathy(:,:) )  
     1638            CALL lbc_lnk( risfdep, 'T', 1. )  
     1639            CALL lbc_lnk( bathy, 'T', 1. ) 
     1640            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1641            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1642            mbathy(:,:) = INT( zbathy(:,:) ) 
     1643         ENDIF  
     1644 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
     1645         DO jj = 1, jpjm1 
     1646            DO ji = 1, jpim1 
     1647               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1648                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
     1649               END IF 
     1650            END DO 
     1651         END DO 
     1652         IF( lk_mpp ) THEN  
     1653            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1654            CALL lbc_lnk( zbathy, 'T', 1. )  
     1655            misfdep(:,:) = INT( zbathy(:,:) )  
     1656            CALL lbc_lnk( risfdep, 'T', 1. )  
     1657            CALL lbc_lnk( bathy, 'T', 1. ) 
     1658            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1659            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1660            mbathy(:,:) = INT( zbathy(:,:) ) 
     1661         ENDIF  
     1662 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
     1663         DO jj = 1, jpjm1 
     1664            DO ji = 1, jpim1 
     1665               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1666                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
     1667               END IF 
     1668            END DO 
     1669         END DO 
     1670         IF( lk_mpp ) THEN  
     1671            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1672            CALL lbc_lnk( zbathy, 'T', 1. )  
     1673            misfdep(:,:) = INT( zbathy(:,:) )  
     1674            CALL lbc_lnk( risfdep, 'T', 1. )  
     1675            CALL lbc_lnk( bathy, 'T', 1. ) 
     1676            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1677            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1678            mbathy(:,:) = INT( zbathy(:,:) ) 
     1679         ENDIF  
     1680 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
     1681         DO jj = 1, jpjm1 
     1682            DO ji = 1, jpi 
     1683               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
     1684                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
     1685               END IF 
     1686            END DO 
     1687         END DO 
     1688         IF( lk_mpp ) THEN  
     1689            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1690            CALL lbc_lnk( zbathy, 'T', 1. )  
     1691            misfdep(:,:) = INT( zbathy(:,:) )  
     1692            CALL lbc_lnk( risfdep, 'T', 1. )  
     1693            CALL lbc_lnk( bathy, 'T', 1. ) 
     1694            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1695            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1696            mbathy(:,:) = INT( zbathy(:,:) ) 
     1697         ENDIF  
     1698 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
     1699         DO jj = 1, jpjm1 
     1700            DO ji = 1, jpi 
     1701               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
     1702                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ; 
     1703               END IF 
     1704            END DO 
     1705         END DO 
     1706         IF( lk_mpp ) THEN  
     1707            zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1708            CALL lbc_lnk( zbathy, 'T', 1. )  
     1709            misfdep(:,:) = INT( zbathy(:,:) )  
     1710            CALL lbc_lnk( risfdep, 'T', 1. )  
     1711            CALL lbc_lnk( bathy, 'T', 1. ) 
     1712            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1713            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1714            mbathy(:,:) = INT( zbathy(:,:) ) 
     1715         ENDIF  
     1716 ! if not compatible after all check, mask T 
     1717         DO jj = 1, jpj 
     1718            DO ji = 1, jpi 
     1719               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
     1720                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
     1721               END IF 
     1722            END DO 
     1723         END DO 
     1724  
     1725         WHERE (mbathy(:,:) == 1) 
     1726            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
     1727         END WHERE 
     1728      END DO  
     1729! end check compatibility ice shelf/bathy 
     1730      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1731      WHERE (misfdep(:,:) <= 5) 
     1732         misfdep = 1; risfdep = 0.0_wp; 
     1733      END WHERE 
     1734 
     1735      IF( icompt == 0 ) THEN  
     1736         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
     1737      ELSE  
     1738         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
     1739      ENDIF  
     1740 
     1741      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
     1742      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
     1743 
     1744      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1745       
     1746   END SUBROUTINE 
    10721747 
    10731748   SUBROUTINE zgr_sco 
     
    14452120            DO jk = 1, jpkm1 
    14462121               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    1447                IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
    1448             END DO 
     2122            END DO 
     2123            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
    14492124         END DO 
    14502125      END DO 
Note: See TracChangeset for help on using the changeset viewer.