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 5462 for branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 – NEMO

Ignore:
Timestamp:
2015-06-22T16:40:58+02:00 (9 years ago)
Author:
davestorkey
Message:

Update UKMO/restart_datestamp branch to revision 5442 of trunk. Also clear the SVN keywords again - necessary after the update step.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5420 r5462  
    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 
    … …  
    3536   USE oce               ! ocean variables 
    3637   USE dom_oce           ! ocean domain 
    37    USE sbc_oce           ! surface variable (isf) 
    3838   USE closea            ! closed seas 
    3939   USE c1d               ! 1D vertical configuration 
    … …  
    298298      ENDIF 
    299299 
     300      IF ( ln_isfcav ) THEN 
    300301! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    301302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    302       DO jk = 1, jpkm1 
    303          e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    304       END DO 
    305       e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
    306  
    307       DO jk = 2, jpk 
    308          e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
    309       END DO 
    310       e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     303         DO jk = 1, jpkm1 
     304            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305         END DO 
     306         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308         DO jk = 2, jpk 
     309            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310         END DO 
     311         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     312      END IF 
    311313 
    312314!!gm BUG in s-coordinate this does not work! 
    … …  
    471473         misfdep(:,:)=1 
    472474         ! 
    473          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    474          IF( cp_cfg == "isomip" ) THEN  
    475            !  
    476            risfdep(:,:)=200.e0  
    477            misfdep(:,:)=1  
    478            ij0 = 1 ; ij1 = 40  
    479            DO jj = mj0(ij0), mj1(ij1)  
    480               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    481                 END DO  
    482             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    483            !  
    484          ELSEIF ( cp_cfg == "isomip2" ) THEN 
    485          !  
    486             risfdep(:,:)=0.e0 
    487             misfdep(:,:)=1 
    488             ij0 = 1 ; ij1 = 40 
    489             DO jj = mj0(ij0), mj1(ij1) 
    490                risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
    491             END DO 
    492             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    493          END IF 
    494          ! 
    495475         DEALLOCATE( idta, zdta ) 
    496476         ! 
    … …  
    504484            CALL iom_close( inum ) 
    505485            mbathy(:,:) = INT( bathy(:,:) ) 
     486            ! 
     487            ! CL :  add Amazon deeper 
     488            IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration  
     489                 ii0 = 230   ;   ii1 = 245                   ! Amazon area  
     490                 ij0 = 140   ;   ij1 = 155                   ! no ocean shallower than 30 meters 
     491                 DO ji = mi0(ii0), mi1(ii1) 
     492                    DO jj = mj0(ij0), mj1(ij1) 
     493                       IF( bathy(ji,jj) .LE. 30. .AND. bathy(ji,jj) .GT. 0.0 ) bathy(ji,jj) = 30._wp 
     494                    END DO 
     495                 END DO 
     496                 IF(lwp) WRITE(numout,*) 
     497                 IF(lwp) WRITE(numout,*) '             orca_r1: Amazon area not shallower than 30 meters for: ' 
     498                 IF(lwp) WRITE(numout,*) '                    Longitude index ',ii0, ii0 
     499                 IF(lwp) WRITE(numout,*) '                    Latitude index  ',ij0, ij0 
     500            ENDIF 
    506501            !                                                ! ===================== 
    507502            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    … …  
    534529         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
    535530            CALL iom_open ( 'bathy_meter.nc', inum )  
    536             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     531            IF ( ln_isfcav ) THEN 
     532               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 
     533            ELSE 
     534               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  ) 
     535            END IF 
    537536            CALL iom_close( inum ) 
    538             !   
     537            !                                                 
    539538            risfdep(:,:)=0._wp          
    540539            misfdep(:,:)=1              
    … …  
    584583      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    585584         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    586          WHERE (bathy == risfdep) 
    587             bathy   = 0.0_wp ; risfdep = 0.0_wp 
    588          END WHERE 
     585         IF ( ln_isfcav ) THEN 
     586            WHERE (bathy == risfdep) 
     587               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     588            END WHERE 
     589         END IF 
    589590         ! end patch 
    590591         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    … …  
    961962      !!---------------------------------------------------------------------- 
    962963      !! 
     964      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     965      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
     966      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     967      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     968      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     969      REAL(wp) ::   zmax             ! Maximum depth 
     970      REAL(wp) ::   zdiff            ! temporary scalar 
     971      REAL(wp) ::   zrefdep          ! temporary scalar 
     972      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     973      !!--------------------------------------------------------------------- 
     974      ! 
     975      IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
     976      ! 
     977      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     978      ! 
     979      IF(lwp) WRITE(numout,*) 
     980      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     981      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     982      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     983 
     984      ll_print = .FALSE.                   ! Local variable for debugging 
     985       
     986      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     987         WRITE(numout,*) 
     988         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     989         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     990      ENDIF 
     991 
     992 
     993      ! bathymetry in level (from bathy_meter) 
     994      ! =================== 
     995      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     996      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     997      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     998      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     999      END WHERE 
     1000 
     1001      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     1002      ! find the number of ocean levels such that the last level thickness 
     1003      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     1004      ! e3t_1d is the reference level thickness 
     1005      DO jk = jpkm1, 1, -1 
     1006         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1007         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     1008      END DO 
     1009 
     1010      IF ( ln_isfcav ) CALL zgr_isf 
     1011 
     1012      ! Scale factors and depth at T- and W-points 
     1013      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     1014         gdept_0(:,:,jk) = gdept_1d(jk) 
     1015         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     1016         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     1017         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     1018      END DO 
     1019      !  
     1020      DO jj = 1, jpj 
     1021         DO ji = 1, jpi 
     1022            ik = mbathy(ji,jj) 
     1023            IF( ik > 0 ) THEN               ! ocean point only 
     1024               ! max ocean level case 
     1025               IF( ik == jpkm1 ) THEN 
     1026                  zdepwp = bathy(ji,jj) 
     1027                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1028                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1029                  e3t_0(ji,jj,ik  ) = ze3tp 
     1030                  e3t_0(ji,jj,ik+1) = ze3tp 
     1031                  e3w_0(ji,jj,ik  ) = ze3wp 
     1032                  e3w_0(ji,jj,ik+1) = ze3tp 
     1033                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1034                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1035                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1036                  ! 
     1037               ELSE                         ! standard case 
     1038                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1039                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1040                  ENDIF 
     1041!gm Bug?  check the gdepw_1d 
     1042                  !       ... on ik 
     1043                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1044                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1045                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1046                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1047                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1048                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1049                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1050                  !       ... on ik+1 
     1051                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1052                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1053                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1054               ENDIF 
     1055            ENDIF 
     1056         END DO 
     1057      END DO 
     1058      ! 
     1059      it = 0 
     1060      DO jj = 1, jpj 
     1061         DO ji = 1, jpi 
     1062            ik = mbathy(ji,jj) 
     1063            IF( ik > 0 ) THEN               ! ocean point only 
     1064               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1065               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1066               ! test 
     1067               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1068               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1069                  it = it + 1 
     1070                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1071                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1072                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1073                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1074               ENDIF 
     1075            ENDIF 
     1076         END DO 
     1077      END DO 
     1078      ! 
     1079      IF ( ln_isfcav ) THEN 
     1080      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1081         DO jj = 1, jpj  
     1082            DO ji = 1, jpi  
     1083               ik = misfdep(ji,jj)  
     1084               IF( ik > 1 ) THEN               ! ice shelf point only  
     1085                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1086                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1087!gm Bug?  check the gdepw_0  
     1088               !       ... on ik  
     1089                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1090                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1091                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1092                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1093                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1094 
     1095                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1096                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1097                  ENDIF  
     1098               !       ... on ik / ik-1  
     1099                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1100                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1101! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1102                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1103               ENDIF  
     1104            END DO  
     1105         END DO  
     1106      !  
     1107         it = 0  
     1108         DO jj = 1, jpj  
     1109            DO ji = 1, jpi  
     1110               ik = misfdep(ji,jj)  
     1111               IF( ik > 1 ) THEN               ! ice shelf point only  
     1112                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1113                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1114               ! test  
     1115                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1116                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1117                     it = it + 1  
     1118                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1119                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1120                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1121                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1122                  ENDIF  
     1123               ENDIF  
     1124            END DO  
     1125         END DO  
     1126      END IF 
     1127      ! END (ISF) 
     1128 
     1129      ! Scale factors and depth at U-, V-, UW and VW-points 
     1130      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1131         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1132         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1133         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1134         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1135      END DO 
     1136      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1137         DO jj = 1, jpjm1 
     1138            DO ji = 1, fs_jpim1   ! vector opt. 
     1139               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1140               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1141               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1142               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1143            END DO 
     1144         END DO 
     1145      END DO 
     1146      IF ( ln_isfcav ) THEN 
     1147      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1148         DO jj = 2, jpjm1  
     1149            DO ji = 2, fs_jpim1   ! vector opt.  
     1150               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1151               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1152               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1153                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1154               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1155               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1156               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1157                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1158            END DO 
     1159         END DO 
     1160      END IF 
     1161 
     1162      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1163      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1164      ! 
     1165      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1166         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1167         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1168         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1169         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1170      END DO 
     1171       
     1172      ! Scale factor at F-point 
     1173      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1174         e3f_0(:,:,jk) = e3t_1d(jk) 
     1175      END DO 
     1176      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1177         DO jj = 1, jpjm1 
     1178            DO ji = 1, fs_jpim1   ! vector opt. 
     1179               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1180            END DO 
     1181         END DO 
     1182      END DO 
     1183      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1184      ! 
     1185      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1186         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1187      END DO 
     1188!!gm  bug ? :  must be a do loop with mj0,mj1 
     1189      !  
     1190      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1191      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1192      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1193      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1194      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1195 
     1196      ! Control of the sign 
     1197      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1198      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1199      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1200      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1201      
     1202      ! Compute gdep3w_0 (vertical sum of e3w) 
     1203      IF ( ln_isfcav ) THEN ! if cavity 
     1204         WHERE (misfdep == 0) misfdep = 1 
     1205         DO jj = 1,jpj 
     1206            DO ji = 1,jpi 
     1207               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1208               DO jk = 2, misfdep(ji,jj) 
     1209                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1210               END DO 
     1211               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)) 
     1212               DO jk = misfdep(ji,jj) + 1, jpk 
     1213                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1214               END DO 
     1215            END DO 
     1216         END DO 
     1217      ELSE ! no cavity 
     1218         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1219         DO jk = 2, jpk 
     1220            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1221         END DO 
     1222      END IF 
     1223      !                                               ! ================= ! 
     1224      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1225         !                                            ! ================= ! 
     1226         DO jj = 1,jpj 
     1227            DO ji = 1, jpi 
     1228               ik = MAX( mbathy(ji,jj), 1 ) 
     1229               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1230               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1231               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1232               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1233               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1234               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1235            END DO 
     1236         END DO 
     1237         WRITE(numout,*) 
     1238         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1239         WRITE(numout,*) 
     1240         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1241         WRITE(numout,*) 
     1242         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1243         WRITE(numout,*) 
     1244         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1245         WRITE(numout,*) 
     1246         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1247         WRITE(numout,*) 
     1248         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1249      ENDIF   
     1250      ! 
     1251      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1252      ! 
     1253      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1254      ! 
     1255   END SUBROUTINE zgr_zps 
     1256 
     1257   SUBROUTINE zgr_isf 
     1258      !!---------------------------------------------------------------------- 
     1259      !!                    ***  ROUTINE zgr_isf  *** 
     1260      !!    
     1261      !! ** Purpose :   check the bathymetry in levels 
     1262      !!    
     1263      !! ** Method  :   THe water column have to contained at least 2 cells 
     1264      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1265      !!                this criterion. 
     1266      !!                  
     1267      !!    
     1268      !! ** Action  : - test compatibility between isfdraft and bathy  
     1269      !!              - bathy and isfdraft are modified 
     1270      !!---------------------------------------------------------------------- 
     1271      !!    
    9631272      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    9641273      INTEGER  ::   ik, it           ! temporary integers 
    … …  
    9711280      REAL(wp) ::   zdiff            ! temporary scalar 
    9721281      REAL(wp) ::   zrefdep          ! temporary scalar 
    973       REAL(wp) ::   zbathydiff, zrisfdepdiff  
    974       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    975       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    976       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     1282      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1283      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1284      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    9771285      !!--------------------------------------------------------------------- 
    9781286      ! 
    979       IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    980       ! 
    981       CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     1287      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1288      ! 
    9821289      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    983       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    984       ! 
    985       IF(lwp) WRITE(numout,*) 
    986       IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
    987       IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    988       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    989  
    990       ll_print = .FALSE.                   ! Local variable for debugging 
    991        
    992       IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    993          WRITE(numout,*) 
    994          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    995          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    996       ENDIF 
    997  
    998       ! bathymetry in level (from bathy_meter) 
    999       ! =================== 
    1000       zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
    1001       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    1002       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
    1003       ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
    1004       END WHERE 
    1005  
    1006       ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    1007       ! find the number of ocean levels such that the last level thickness 
    1008       ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
    1009       ! e3t_1d is the reference level thickness 
    1010       DO jk = jpkm1, 1, -1 
    1011          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1012          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    1013       END DO 
     1290      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1291 
     1292 
    10141293      ! (ISF) compute misfdep 
    10151294      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    … …  
    10551334            misfdep(jpi,:) = misfdep(  2  ,:)  
    10561335         ENDIF 
    1057   
     1336 
    10581337         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10591338            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10601339            mbathy(jpi,:) = mbathy(  2  ,:) 
    10611340         ENDIF 
    1062   
     1341 
    10631342         ! split last cell if possible (only where water column is 2 cell or less) 
    10641343         DO jk = jpkm1, 1, -1 
    … …  
    10781357            END WHERE 
    10791358         END DO 
    1080   
     1359 
    10811360  
    10821361 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    … …  
    12541533  
    12551534 ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1256          DO jk = 1, jpk 
     1535         DO jk = 2, jpk 
    12571536            WHERE (misfdep==0) misfdep=jpk 
    12581537            zmask=0 
    … …  
    13591638               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    13601639               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1361                IF( ibtest == 0 ) THEN 
     1640               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    13621641                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    13631642               END IF 
    … …  
    14751754      ENDIF  
    14761755 
    1477       ! Scale factors and depth at T- and W-points 
    1478       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    1479          gdept_0(:,:,jk) = gdept_1d(jk) 
    1480          gdepw_0(:,:,jk) = gdepw_1d(jk) 
    1481          e3t_0  (:,:,jk) = e3t_1d  (jk) 
    1482          e3w_0  (:,:,jk) = e3w_1d  (jk) 
    1483       END DO 
    1484       !  
    1485       DO jj = 1, jpj 
    1486          DO ji = 1, jpi 
    1487             ik = mbathy(ji,jj) 
    1488             IF( ik > 0 ) THEN               ! ocean point only 
    1489                ! max ocean level case 
    1490                IF( ik == jpkm1 ) THEN 
    1491                   zdepwp = bathy(ji,jj) 
    1492                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1493                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1494                   e3t_0(ji,jj,ik  ) = ze3tp 
    1495                   e3t_0(ji,jj,ik+1) = ze3tp 
    1496                   e3w_0(ji,jj,ik  ) = ze3wp 
    1497                   e3w_0(ji,jj,ik+1) = ze3tp 
    1498                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1499                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1500                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1501                   ! 
    1502                ELSE                         ! standard case 
    1503                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1504                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1505                   ENDIF 
    1506 !gm Bug?  check the gdepw_1d 
    1507                   !       ... on ik 
    1508                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1509                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1510                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1511                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1512                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1513                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1514                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1515                   !       ... on ik+1 
    1516                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1517                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1518                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1519                ENDIF 
    1520             ENDIF 
    1521          END DO 
    1522       END DO 
    1523       ! 
    1524       it = 0 
    1525       DO jj = 1, jpj 
    1526          DO ji = 1, jpi 
    1527             ik = mbathy(ji,jj) 
    1528             IF( ik > 0 ) THEN               ! ocean point only 
    1529                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1530                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1531                ! test 
    1532                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1533                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1534                   it = it + 1 
    1535                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1536                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1537                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1538                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1539                ENDIF 
    1540             ENDIF 
    1541          END DO 
    1542       END DO 
    1543       ! 
    1544       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1545       DO jj = 1, jpj  
    1546          DO ji = 1, jpi  
    1547             ik = misfdep(ji,jj)  
    1548             IF( ik > 1 ) THEN               ! ice shelf point only  
    1549                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1550                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1551 !gm Bug?  check the gdepw_0  
    1552                !       ... on ik  
    1553                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1554                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1555                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1556                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1557                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1558  
    1559                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1560                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1561                 ENDIF  
    1562                !       ... on ik / ik-1  
    1563                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1564                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1565 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1566                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1567             ENDIF  
    1568          END DO  
    1569       END DO  
    1570       !  
    1571       it = 0  
    1572       DO jj = 1, jpj  
    1573          DO ji = 1, jpi  
    1574             ik = misfdep(ji,jj)  
    1575             IF( ik > 1 ) THEN               ! ice shelf point only  
    1576                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1577                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1578                ! test  
    1579                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1580                IF( zdiff <= 0. .AND. lwp ) THEN   
    1581                   it = it + 1  
    1582                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1583                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1584                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1585                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1586                ENDIF  
    1587             ENDIF  
    1588          END DO  
    1589       END DO  
    1590       ! END (ISF) 
    1591  
    1592       ! Scale factors and depth at U-, V-, UW and VW-points 
    1593       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1594          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1595          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1596          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1597          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1598       END DO 
    1599       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1600          DO jj = 1, jpjm1 
    1601             DO ji = 1, fs_jpim1   ! vector opt. 
    1602                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1603                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1604                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1605                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1606             END DO 
    1607          END DO 
    1608       END DO 
    1609       ! (ISF) define e3uw 
    1610       DO jk = 2,jpk                           
    1611          DO jj = 1, jpjm1  
    1612             DO ji = 1, fs_jpim1   ! vector opt.  
    1613                e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1614                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1615                e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1616                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1617             END DO  
    1618          END DO  
    1619       END DO 
    1620       !End (ISF) 
    1621        
    1622       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1623       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1624       ! 
    1625       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1626          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1627          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1628          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1629          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1630       END DO 
    1631        
    1632       ! Scale factor at F-point 
    1633       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1634          e3f_0(:,:,jk) = e3t_1d(jk) 
    1635       END DO 
    1636       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1637          DO jj = 1, jpjm1 
    1638             DO ji = 1, fs_jpim1   ! vector opt. 
    1639                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1640             END DO 
    1641          END DO 
    1642       END DO 
    1643       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1644       ! 
    1645       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1646          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1647       END DO 
    1648 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1649       !  
    1650       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1651       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1652       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1653       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1654       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1655  
    1656       ! Control of the sign 
    1657       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1658       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1659       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1660       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1661       
    1662       ! Compute gdep3w_0 (vertical sum of e3w) 
    1663       WHERE (misfdep == 0) misfdep = 1 
    1664       DO jj = 1,jpj 
    1665          DO ji = 1,jpi 
    1666             gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1667             DO jk = 2, misfdep(ji,jj) 
    1668                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1669             END DO 
    1670             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)) 
    1671             DO jk = misfdep(ji,jj) + 1, jpk 
    1672                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1673             END DO 
    1674         END DO 
    1675       END DO 
    1676       !                                               ! ================= ! 
    1677       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1678          !                                            ! ================= ! 
    1679          DO jj = 1,jpj 
    1680             DO ji = 1, jpi 
    1681                ik = MAX( mbathy(ji,jj), 1 ) 
    1682                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1683                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1684                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1685                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1686                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1687                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1688             END DO 
    1689          END DO 
    1690          WRITE(numout,*) 
    1691          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1692          WRITE(numout,*) 
    1693          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1694          WRITE(numout,*) 
    1695          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1696          WRITE(numout,*) 
    1697          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1698          WRITE(numout,*) 
    1699          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1700          WRITE(numout,*) 
    1701          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1702       ENDIF   
    1703       ! 
    1704       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    17051756      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    17061757      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1707       ! 
    1708       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1709       ! 
    1710    END SUBROUTINE zgr_zps 
     1758 
     1759      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1760       
     1761   END SUBROUTINE 
    17111762 
    17121763   SUBROUTINE zgr_sco 
Note: See TracChangeset for help on using the changeset viewer.