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 – NEMO

Changeset 4747


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).

Location:
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/CONFIG/ISOMIP/EXP00/namelist_cfg

    r4666 r4747  
    2424&namrun        !   parameters of the run 
    2525!----------------------------------------------------------------------- 
    26    nn_no       =  1   !  job number (no more used...) 
     26   nn_no       =  1        !  job number (no more used...) 
    2727   cn_exp      =  "ISOMIP" !  experience name 
    28    nn_it000    =  1   !  first time step 
     28   nn_it000    =  1        !  first time step 
    2929   nn_itend    =  490560   !  last  time step (std 5475) 
    3030   nn_date0    =  010101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    3131   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    32    ln_rstart   = .false. !  start from rest (F) or from a restart file (T) 
     32   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    3333   nn_euler    =       1   !  = 0 : start with forward time step if ln_rstart=.true. 
    3434   nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T 
    … …  
    8282   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F) 
    8383   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F) 
     84   ln_isfcav   = .true.    !  ice shelf cavity 
    8485/ 
    8586!----------------------------------------------------------------------- 
    … …  
    165166!!   namtra_qsr      penetrative solar radiation 
    166167!!   namsbc_rnf      river runoffs 
     168!!   namsbc_isf      ice shelf melting/freezing 
    167169!!   namsbc_apr      Atmospheric Pressure 
    168170!!   namsbc_ssr      sea surface restoring term (for T and/or S) 
    … …  
    184186                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    185187   ln_rnf      = .false.   !  runoffs                                   (T => fill namsbc_rnf) 
    186    nn_isf      = 1         !  0 = no isf               / 1 = presence of ISF  
     188   nn_isf      = 1         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
     189                           !  0 = no isf               / 1 = presence of ISF  
    187190                           !  2 = bg03 parametrisation / 3 = rnf file for isf 
    188191                           !  4 = ISF are prescribed 
     192                           !  options 1 and 4 need ln_isfcav = .true. (domzgr) 
    189193   ln_ssr      = .false.   !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    190194   nn_fwb      = 1         !  FreshWater Budget: =0 unchecked 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4666 r4747  
    8282   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F) 
    8383   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F) 
     84   ln_isfcav   = .false.   !  ice shelf cavity 
    8485/ 
    8586!----------------------------------------------------------------------- 
    … …  
    208209!!   namtra_qsr      penetrative solar radiation 
    209210!!   namsbc_rnf      river runoffs 
     211!!   namsbc_isf      ice shelf melting/freezing 
    210212!!   namsbc_apr      Atmospheric Pressure 
    211213!!   namsbc_ssr      sea surface restoring term (for T and/or S) 
    … …  
    232234                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    233235   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    234    ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
    235    nn_isf      = 1         !  0=no isf                   1 = presence of ISF  
     236   ln_rnf      = .true.    !  runoffs                                   (T   => fill namsbc_rnf) 
     237   nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
     238                           !  0 =no isf                  1 = presence of ISF  
    236239                           !  2 = bg03 parametrisation   3 = rnf file for isf    
    237240                           !  4 = ISF fwf specified 
     241                           !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    238242   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    239243   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
    … …  
    640644   rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2) 
    641645   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T  
    642    ln_loglayer = .false.   !  logarithmic formulation (non linear case) 
    643646   ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file ) 
    644647   rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T) 
     648   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case) 
     649   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T 
     650   rn_tfri2_max =   1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T) 
     651   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2) 
     652   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T 
     653   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file ) 
     654   rn_tfrien   =    50.    !  local multiplying factor of tfr (ln_tfr2d=T) 
     655 
    645656   ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true) 
     657   ln_loglayer = .false.   !  logarithmic formulation (non linear case) 
    646658/ 
    647659!----------------------------------------------------------------------- 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r4666 r4747  
    193193         CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity 
    194194      END IF 
     195      ! multiply by umask to prevent not numerical value error in the ioserver sometimes 
    195196      IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 
    196          CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport 
    197          CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport 
     197         CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport 
     198         CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport 
    198199      ELSE 
    199          CALL iom_put( "uoce" , un                         )    ! i-current 
    200          CALL iom_put( "voce" , vn                         )    ! j-current 
     200         CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:)                  )    ! i-current 
     201         CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:)                  )    ! j-current 
    201202         DO jj = 1, jpj 
    202203            DO ji = 1, jpi 

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

    r4726 r4747  
    175175   LOGICAL, PUBLIC ::   ln_zps        !: z-coordinate - partial step 
    176176   LOGICAL, PUBLIC ::   ln_sco        !: s-coordinate or hybrid z-s coordinate 
    177    LOGICAL, PUBLIC ::   ln_isf        !: presence of ISF  
     177   LOGICAL, PUBLIC ::   ln_isfcav     !: presence of ISF  
    178178 
    179179   !! All coordinates 

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

    r4726 r4747  
    849849               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    850850               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 
     851               ! needed to restart if land processor not computed  
     852               WHERE ( tmask(:,:,:) == 0.0_wp )  
     853                  fse3t_n(:,:,:) = e3t_0(:,:,:) 
     854                  fse3t_b(:,:,:) = e3t_0(:,:,:) 
     855               END WHERE 
    851856               IF( neuler == 0 ) THEN 
    852857                  fse3t_b(:,:,:) = fse3t_n(:,:,:) 

  • 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  

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4726 r4747  
    137137      READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 
    138138902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 
    139       WRITE ( numond, namdyn_hpg ) 
     139      IF(lwm) WRITE ( numond, namdyn_hpg ) 
    140140      ! 
    141141      IF(lwp) THEN                   ! Control print 
    … …  
    397397      ! iniitialised to 0. zhpi zhpi  
    398398      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 
    399        
    400       ! assume density of water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
     399 
     400!==================================================================================      
     401!=====Compute iceload and contribution of the half first wet layer =================  
     402!=================================================================================== 
     403 
     404      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
    401405      ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
     406 
     407      ! compute density of the water displaced by the ice shelf  
    402408      zrhd = rhd ! save rhd 
    403409      DO jk = 1, jpk 
    … …  
    434440            END DO 
    435441            IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 
    436                                &                              * ( fsdepw(ji,jj,ikt) - gdept_1d(ikt-1) ) 
     442                               &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
    437443         END DO 
    438444      END DO 
    … …  
    464470         END DO 
    465471      END DO 
    466       ! 
     472!==================================================================================      
     473!===== Compute partial cell contribution for the top cell =========================  
     474!================================================================================== 
    467475      DO jj = 2, jpjm1 
    468476         DO ji = fs_2, fs_jpim1   ! vector opt. 
    469477            iku = miku(ji,jj) ;  
    470478            zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 
    471             ze3wu  = fse3w(ji+1,jj,iku+1)-fse3w(ji,jj,iku+1)  
     479            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    472480            ! u direction 
    473481            IF ( iku .GT. 1 ) THEN 
    … …  
    478486               ! corrective term ( = 0 if z coordinate ) 
    479487               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
    480                ! zhpi will be added in interior loop and zuapint will be removed in the interior loop 
     488               ! zhpi will be added in interior loop 
    481489               ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap 
    482490               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    483491               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku) 
    484492 
    485                ! zhpi will be added in interior loop and zuapint will be removed in the interior loop 
    486                ! case iku + 1 (remove the zphi term added in the interior loop) 
    487                zhpiint        =  zcoef0 / e1u(ji,jj)                                                             &     
     493               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 
     494               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &     
    488495                  &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
    489496                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
    … …  
    492499               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    493500            END IF 
     501                
    494502            ! v direction 
    495503            ikv = mikv(ji,jj) 
    496             ze3wv  = fse3w(ji,jj+1,ikv+1)-fse3w(ji,jj,ikv+1)  
     504            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    497505            IF ( ikv .GT. 1 ) THEN 
    498506               ! case ikv 
    … …  
    502510               ! corrective term ( = 0 if z coordinate ) 
    503511               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
    504                ! zhpi will be added in interior loop and zvapint will be removed in the interior loop 
     512               ! zhpi will be added in interior loop 
    505513               va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap 
    506514               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    507515               IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv)  
    508                ! zhpi will be added in interior loop and zvapint will be removed in the interior loop 
    509                ! case ikv + 1 (remove the zphj term added in the interior loop) 
    510                zhpjint        =  zcoef0 / e2v(ji,jj)                                                            & 
     516                
     517               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 
     518               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              & 
    511519                  &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
    512                   &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)  & 
     520                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    & 
    513521                  &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
    514522                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  ) 
    … …  
    518526      END DO 
    519527 
    520       ! interior value (2=<jk=<jpkm1) 
     528!==================================================================================      
     529!===== Compute interior value =====================================================  
     530!================================================================================== 
     531 
    521532      DO jj = 2, jpjm1 
    522533         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    524535            DO jk = 2, jpkm1 
    525536               ! hydrostatic pressure gradient along s-surfaces 
    526                ! s-coordinate pressure gradient correction 
     537               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    527538               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              & 
    528539                  &                            + zcoef0 / e1u(ji,jj)                                                           & 
    … …  
    531542                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 & 
    532543                  &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
    533              !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
    534                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     544               ! s-coordinate pressure gradient correction 
     545               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
     546               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    & 
    535547                  &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
    536548               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    537549 
     550               ! hydrostatic pressure gradient along s-surfaces 
     551               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    538552               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              & 
    539553                  &                            + zcoef0 / e2v(ji,jj)                                                           & 
    … …  
    542556                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 & 
    543557                  &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
    544              !corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
    545                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     558               ! s-coordinate pressure gradient correction 
     559               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
     560               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     & 
    546561                  &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
    547562               ! add to the general momentum trend 
    … …  
    550565         END DO 
    551566      END DO 
    552       
    553       ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90) 
     567 
     568!==================================================================================      
     569!===== Compute bottom cell contribution (partial cell) ============================  
     570!================================================================================== 
     571 
    554572# if defined key_vectopt_loop 
    555573         jj = 1 
    … …  
    589607# endif 
    590608      END DO 
    591  
    592  
     609      
     610      ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 
    593611      rhd = zrhd 
    594612      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r4666 r4747  
    220220      ENDIF 
    221221      ! 
    222       IF( lk_vvl ) THEN     ! need it for (ISF) 
    223                      CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    224                      CALL iom_get( numror, jpdom_autoglo, 'fse3t  ', fse3t  (:,:,:) ) 
    225                      CALL iom_get( numror, jpdom_autoglo, 'fse3w  ', fse3w  (:,:,:) ) 
    226                      CALL iom_get( numror, jpdom_autoglo, 'fsdepw ', fsdepw (:,:,:) ) 
    227       END IF 
    228222      CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields 
    229223      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn      ) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

    r4647 r4747  
    6767         imask                                ! temporary global workspace 
    6868      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   & 
    69          zdta                   ! temporary data workspace 
     69         zdta, zdtaisf                     ! temporary data workspace 
    7070      REAL(wp) ::   zidom , zjdom          ! temporary scalars 
    7171 
    7272      ! read namelist for ln_zco 
    73       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
     73      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
    7474 
    7575      !!---------------------------------------------------------------------- 
    … …  
    109109      ENDIF 
    110110      CALL iom_close (inum) 
     111       
     112      ! used to compute the land processor in case of not masked bathy file. 
     113      zdtaisf(:,:) = 0.0_wp 
     114      IF ( ln_isfcav ) THEN 
     115         CALL iom_open ( 'bathy_meter.nc', inum )   ! Meter bathy in case of partial steps 
     116         CALL iom_get ( inum, jpdom_unknown, 'isf_draft' , zdtaisf, kstart=(/jpizoom,jpjzoom/), kcount=(/jpiglo,jpjglo/) ) 
     117      END IF 
     118      CALL iom_close (inum) 
    111119 
    112120      ! land/sea mask over the global/zoom domain 
    113121 
    114122      imask(:,:)=1 
    115       WHERE ( zdta(:,:) <= 0. ) imask = 0 
     123      WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0 
    116124 
    117125      !  1. Dimension arrays for subdomains 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r4724 r4747  
    130130               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    131131               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    132                ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
    133                ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
     132               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
     133               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
     134               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
     135               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
     136               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    134137               ! 
    135138               ! i- direction 
    … …  
    189192               iku = mbku(ji,jj) 
    190193               ikv = mbkv(ji,jj) 
    191                ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
    192                ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
    193  
    194                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    195                ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    196                ENDIF 
    197                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
    198                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     194               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
     195               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
     196 
     197               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
     198               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
     199               ENDIF 
     200               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1 
     201               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2 
    199202               ENDIF 
    200203# if ! defined key_vectopt_loop 
    … …  
    218221               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    219222               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    220                ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
    221                ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
     223               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
     224               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    222225               IF( ze3wu >= 0._wp ) THEN  
    223226                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    224                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i:  
     227                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1  
    225228                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    226229                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
    … …  
    267270               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
    268271               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
    269                ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
    270                ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
    271272               ! 
     273               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
     274               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
     275               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
     276               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
     277               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))  
     278               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    272279               ! i- direction 
    273280               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    … …  
    322329               iku = miku(ji,jj) 
    323330               ikv = mikv(ji,jj) 
    324                ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
    325                ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
    326  
    327                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    328                ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    329                ENDIF 
    330                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
    331                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     331               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
     332               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
     333 
     334               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
     335               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
     336               ENDIF 
     337               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
     338               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
    332339               ENDIF 
    333340# if ! defined key_vectopt_loop 
    … …  
    351358               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    352359               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
    353                ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
    354                ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
     360               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
     361               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    355362               IF( ze3wu >= 0._wp ) THEN 
    356363                 sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4724 r4747  
    167167               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv 
    168168               ! 
     169               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     170               IF ( miku(ji,jj) + 2 .LE. mbku(ji,jj) ) THEN 
     171                  bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
     172                               &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
     173                               &          * zecu * (1._wp - umask(ji,jj,1)) 
     174               END IF 
     175               IF ( mikv(ji,jj) + 2 .LE. mbkv(ji,jj) ) THEN 
     176                  bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
     177                               &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
     178                               &          * zecv 
     179               END IF 
    169180               ! (ISF) ======================================================================== 
    170181               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
    … …  
    182193               tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
    183194               ! (ISF) END ==================================================================== 
    184  
     195               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     196               IF ( miku(ji,jj) + 2 .LE. mbku(ji,jj) ) THEN 
     197                  tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
     198                               &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
     199                               &          * zecu * (1._wp - umask(ji,jj,1)) 
     200               END IF 
     201               IF ( mikv(ji,jj) + 2 .LE. mbkv(ji,jj) ) THEN 
     202                  tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
     203                               &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
     204                               &          * zecv * (1._wp - vmask(ji,jj,1)) 
     205               END IF 
    185206            END DO 
    186207         END DO 
Note: See TracChangeset for help on using the changeset viewer.