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

Changeset 6006


Ignore:
Timestamp:
2015-12-04T17:56:07+01:00 (8 years ago)
Author:
mathiot
Message:

Merged ice sheet coupling branch

Location:
branches/2015/dev_MetOffice_merge_2015
Files:
27 edited
4 copied

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DOM.tex

    r5120 r6006  
    498498Contrary to the horizontal grid, the vertical grid is computed in the code and no  
    499499provision is made for reading it from a file. The only input file is the bathymetry  
    500 (in meters) (\ifile{bathy\_meter})  
     500(in meters) (\ifile{bathy\_meter}).  
    501501\footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the  
    502502\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point  
    503503in each water column is by-passed}.  
     504If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft  
     505(in meters) (\ifile{isf\_draft\_meter}) is needed and all the location where the isf cavity thinnest 
     506 than \np{rn\_isfhmin} meters are grounded (ie masked).  
     507 
    504508After reading the bathymetry, the algorithm for vertical grid definition differs  
    505509between the different options: 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex

    r5120 r6006  
    996996 at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}  
    997997(called from \mdl{divcur}).  
     998 
     999\section{ Ice sheet coupling} 
     1000\label{SBC_iscpl} 
     1001%------------------------------------------namsbc_iscpl---------------------------------------------------- 
     1002\namdisplay{namsbc_iscpl} 
     1003%-------------------------------------------------------------------------------------------------------- 
     1004Ice sheet/ocean coupling is done through file exchange at the restart step. NEMO, at each restart step,  
     1005read the bathymetry and ice shelf draft variable in a netcdf file.  
     1006If \np{ln\_iscpl = ~true}, the isf draft is assume to be different at each restart step  
     1007with potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics. 
     1008The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases: 
     1009\begin{description} 
     1010\item[Thin a cell down:] 
     1011   T/S/ssh are unchanged and U/V in the top cell are corrected to keep the barotropic transport (bt) constant ($bt_b=bt_n$). 
     1012\item[Enlarge  a cell:] 
     1013   See case "Thin a cell down" 
     1014\item[Dry a cell:] 
     1015   mask, T/S, U/V and ssh are set to 0. Furthermore, U/V into the water column are modified to satisfy ($bt_b=bt_n$). 
     1016\item[Wet a cell:]  
     1017   mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$ and U/V set to 0. If no neighbours along i,j and k, T/S/U/V and mask are set to 0. 
     1018\item[Dry a column:] 
     1019   mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0. 
     1020\item[Wet a column:] 
     1021   set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0. If no neighbour, T/S/U/V and mask set to 0. 
     1022\end{description} 
     1023The extrapolation is call \np{nn\_drown} times. It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps, 
     1024 the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\ 
     1025This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe.  
     1026This could lead to a trend in heat/salt content and volume. In order to remove the trend and keep the conservation level as close to 0 as possible, 
     1027 a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location.  
     1028Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added, 
     1029 over a period of \np{rn\_fiscpl} time step, into the system.  
     1030So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\ 
     1031 
     1032As the before and now fields are not compatible (modification of the geometry), the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$. 
    9981033% 
    9991034% ================================================================ 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namdom

    r4560 r6006  
    66   nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
    77   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
     8   rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice 
    89   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
    910   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namrun

    r4147 r6006  
    1515   cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input) 
    1616   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
     17   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model (ln_iscpl = T => namsbc_iscpl) 
    1718   nn_istate   =       0   !  output the initial state (1) or not (0) 
    1819   nn_stock    =    5475   !  frequency of creation of a restart file (modulo referenced to 1) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml

    r6005 r6006  
    200200 
    201201         <!-- * variable related to ice shelf forcing * --> 
    202          <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
    203          <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     /> 
    204          <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      /> 
    205          <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      /> 
    206          <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="1e-3"     /> 
    207          <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degC"     /> 
     202         <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  /> 
     203         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     /> 
     204         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "                 unit="m/s"      /> 
     205         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "                 unit="m/s"      /> 
     206    <field id="stbl"         long_name="salinity in the Losh tbl                    "                 unit="PSU"      /> 
     207    <field id="ttbl"         long_name="temperature in the Losh tbl                 "                 unit="C"        /> 
     208    <field id="utbl"         long_name="zonal current in the Losh tbl at T point    "                 unit="m/s"      /> 
     209    <field id="vtbl"         long_name="merid current in the Losh tbl at T point    "                 unit="m/s"      /> 
     210    <field id="thermald"     long_name="thermal driving of ice shelf melting        "                 unit="C"        /> 
     211    <field id="tfrz"         long_name="top freezing point (used to compute melt)   "                 unit="C"        /> 
     212    <field id="tinsitu"      long_name="top insitu temperature (used to cmpt melt)  "                 unit="C"        /> 
     213    <field id="ustar"        long_name="ustar at T point used in ice shelf melting  "                 unit="m/s"      /> 
    208214 
    209215         <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6005 r6006  
    4343   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
    4444   cn_ocerst_outdir = "."      !  directory in which to write output ocean restarts 
     45   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4546   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4647   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
    … …  
    128129   nn_msh      =    1      !  create (=1) a mesh file or not (=0) 
    129130   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
     131   rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice 
    130132   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
    131133   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 
    … …  
    212214!!   namsbc_rnf      river runoffs 
    213215!!   namsbc_isf      ice shelf melting/freezing 
     216!!   namsbc_iscpl    coupling option between land ice model and ocean 
    214217!!   namsbc_apr      Atmospheric Pressure 
    215218!!   namsbc_ssr      sea surface restoring term (for T and/or S) 
    … …  
    456459/ 
    457460!----------------------------------------------------------------------- 
     461&namsbc_iscpl  !   land ice / ocean coupling option                      
     462!----------------------------------------------------------------------- 
     463   nn_drown  = 10       ! number of iteration of the extrapolation loop (fill the new wet cells) 
     464   ln_hsb    = .false.  ! activate conservation module (conservation exact after a time of rn_fiscpl) 
     465   nn_fiscpl = 43800    ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 
     466/ 
     467!----------------------------------------------------------------------- 
    458468&namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk 
    459469!----------------------------------------------------------------------- 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90

    r5836 r6006  
    116116         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
    117117         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    118          &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
    119       NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   & 
     118         &             nn_write, ln_iscpl, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
     119      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin, rn_isfhmin,  & 
    120120         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            & 
    121121         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, & 
    … …  
    813813      DO jj = 1, jpjm1 
    814814         DO ji = 1, fs_jpim1   ! vector loop 
    815             umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    816             vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     815            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     816            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    817817         END DO 
    818818         DO ji = 1, jpim1      ! NO vector opt. 
    819             fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
     819            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    820820               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    821821         END DO 
    822822      END DO 
    823       CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
    824       CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
    825       CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
     823      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions 
     824      CALL lbc_lnk( ssvmask, 'V', 1._wp ) 
     825      CALL lbc_lnk( ssfmask, 'F', 1._wp ) 
    826826 
    827827      ! 3. Ocean/land mask at wu-, wv- and w points  

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90

    r5836 r6006  
    460460      ENDIF 
    461461 
    462       IF( nn_timing == 1 )   CALL timing_start('dia_fwb') 
     462      IF( nn_timing == 1 )   CALL timing_stop('dia_fwb') 
    463463 
    464464 9005 FORMAT(1X,A,ES24.16) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r5930 r6006  
    135135      DO jk=1,nb_ana 
    136136       DO ji=1,jpmax_harmo 
    137           IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN 
     137          IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 
    138138             name(jk) = ji 
    139139             EXIT 
    … …  
    195195                     ! Elevation 
    196196                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)         
    197                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj) 
    198                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj) 
     197                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
     198                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
    199199                  END DO 
    200200               END DO 
    … …  
    324324               X1= ana_amp(ji,jj,jh,1) 
    325325               X2=-ana_amp(ji,jj,jh,2) 
    326                out_u(ji,jj,       jh) = X1 * umask_i(ji,jj) 
    327                out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj) 
     326               out_u(ji,jj,       jh) = X1 * ssumask(ji,jj) 
     327               out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 
    328328            ENDDO 
    329329         ENDDO 
    … …  
    488488            DO jj_sd = ji_sd, ninco 
    489489               zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
    490                IF( zval2.GE.zval1 )THEN 
     490               IF( zval2 >= zval1 )THEN 
    491491                  ipivot(ji_sd) = jj_sd 
    492492                  zval1         = zval2 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5643 r6006  
    4646   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends 
    4747   ! 
    48    REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          ! 
     48   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf  
     49   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          ! 
    4950   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   ! 
    5051   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
    … …  
    138139      ! 2 -  Content variations ! 
    139140      ! ------------------------ ! 
     141      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl) 
    140142      zdiff_v2 = 0._wp 
    141143      zdiff_hc = 0._wp 
    … …  
    143145 
    144146      ! volume variation (calculated with ssh) 
    145       zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) ) 
     147      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) ) 
    146148 
    147149      ! heat & salt content variation (associated with ssh) 
    … …  
    158160            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) )  
    159161         END IF 
    160          z_ssh_hc = glob_sum( z2d0 )  
    161          z_ssh_sc = glob_sum( z2d1 )  
     162         z_ssh_hc = glob_sum_full( z2d0 )  
     163         z_ssh_sc = glob_sum_full( z2d1 )  
    162164      ENDIF 
    163165 
    164166      DO jk = 1, jpkm1 
    165167         ! volume variation (calculated with scale factors) 
    166          zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) & 
    167             &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
    168          ! heat content variation 
    169          zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) &  
    170             &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) ) 
     168         zdiff_v2 = zdiff_v2 + glob_sum_full( surf    (:,:) * fse3t_n  (:,:,jk) * tmask(:,:,jk) & 
     169              &                             - surf_ini(:,:) *   e3t_ini(:,:,jk) ) 
     170         ! heat content variation  
     171         zdiff_hc = zdiff_hc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) & 
     172              &                             - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 
    171173         ! salt content variation 
    172          zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   & 
    173             &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) ) 
     174         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) &  
     175              &                             - surf_ini(:,:) * sc_loc_ini(:,:,jk) ) 
    174176      ENDDO 
    175177 
    … …  
    191193      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors) 
    192194      DO jk = 1, jpkm1 
    193          zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
     195         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
    194196      END DO 
    195197 
    … …  
    214216        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu) 
    215217        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J)  
    216         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3) 
    217         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)   
     218        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9   )              ! Salt content variation (psu*km3) 
     219        CALL iom_put( 'bgvolssh' , zdiff_v1  * 1.e-9   )              ! volume ssh variation (km3)   
    218220        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    219221        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)  
    … …  
    261263              CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 
    262264           ENDIF 
     265           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling 
    263266           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 
    264267           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) 
    … …  
    273276          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state ' 
    274277          IF(lwp) WRITE(numout,*) '~~~~~~~' 
    275           ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh 
     278          surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface 
     279          ssh_ini(:,:) = sshn(:,:)                          ! initial ssh 
    276280          DO jk = 1, jpk 
    277              e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors 
    278              hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content 
    279              sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content 
     281             ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance). 
     282             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)    * tmask(:,:,jk)                    ! initial vertical scale factors 
     283             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content 
     284             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content 
    280285          END DO 
    281286          frc_v = 0._wp                                           ! volume       trend due to forcing 
    … …  
    312317           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 
    313318        ENDIF 
     319        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )      ! ice sheet coupling 
    314320        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 
    315321        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) 
    … …  
    379385      ! 1 - Allocate memory ! 
    380386      ! ------------------- ! 
    381       ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 
    382          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
     387      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), & 
     388         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  ) 
    383389      IF( ierror > 0 ) THEN 
    384390         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5997 r6006  
    3232   REAL(wp), PUBLIC ::   rn_bathy        !: depth of flat bottom (active if nn_bathy=0; if =0 depth=jpkm1) 
    3333   REAL(wp), PUBLIC ::   rn_hmin         !: minimum ocean depth (>0) or minimum number of ocean levels (<0) 
     34   REAL(wp), PUBLIC ::   rn_isfhmin      !: threshold to discriminate grounded ice to floating ice 
    3435   REAL(wp), PUBLIC ::   rn_e3zps_min    !: miminum thickness for partial steps (meters) 
    3536   REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps 
    … …  
    4344   INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    4445   INTEGER , PUBLIC ::   nn_euler        !: =0 start with forward time step or not (=1) 
     46   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    4547   LOGICAL , PUBLIC ::   ln_crs          !: Apply grid coarsening to dynamical model output or online passive tracers 
    4648 
    … …  
    251253   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level 
    252254   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level 
    253    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters) 
    254    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 
    255    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters) 
     256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask 
     257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book) 
     258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function 
    256259 
    257260   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF) 
    258261   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF) 
    259262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                       (ISF) 
    260    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask                   !: surface domain T-point mask  
    261  
     263 
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts 
    262265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    263266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    … …  
    386389         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 
    387390 
    388       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
    389          &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
    390          &     bmask  (jpi,jpj) ,                                                       & 
    391          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
     391      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       & 
     392         &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       &  
     393         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 
     394         &     bmask(jpi,jpj)   ,                                                          & 
     395         &     mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
    392396 
    393397! (ISF) Allocation of basic array    
    394       ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),                   & 
    395          &      mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
    396          &      mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
     398      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     & 
     399         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
     400         &     mikf(jpi,jpj), STAT=ierr(10) ) 
    397401 
    398402      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5997 r6006  
    112112      END DO 
    113113      !                                        ! Inverse of the local depth 
    114       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
    115       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     114      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 
     115      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 
    116116 
    117117                             CALL dom_stp      ! time step 
    118       IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file 
     118      ! 
     119      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file 
     120      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file 
     121      ! 
    119122      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control 
    120123      ! 
    … …  
    135138      !!---------------------------------------------------------------------- 
    136139      USE ioipsl 
    137       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
    138                        nn_no   , cn_exp  , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,  & 
    139          &             nn_it000, nn_itend, nn_date0    , nn_time0     , nn_leapy  , nn_istate ,  & 
    140          &             nn_stock, nn_write, ln_dimgnnn  , ln_mskland   , ln_clobber, nn_chunksz,  & 
    141          &             nn_euler, ln_cfmeta 
    142       NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   & 
    143          &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  & 
    144          &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    & 
    145          &             jphgr_msh, & 
    146          &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 
    147          &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 
     140      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
     141                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     & 
     142         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
     143         &             nn_stock, nn_write , ln_dimgnnn  , ln_mskland   , ln_clobber, nn_chunksz,     & 
     144         &             nn_euler, ln_cfmeta, ln_iscpl 
     145      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, & 
     146         &             nn_acc   , rn_atfp , rn_rdt      , rn_rdtmin   ,                              & 
     147         &             rn_rdtmax, rn_rdth , nn_closea   , ln_crs      ,                              & 
     148         &             jphgr_msh,                                                                    & 
     149         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
     150         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
    148151         &             ppa2, ppkth2, ppacr2 
    149152#if defined key_netcdf4 
    … …  
    193196         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    194197         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     198         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl 
    195199      ENDIF 
    196200 
    … …  
    260264         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin 
    261265         WRITE(numout,*) '      min number of ocean level (<0)       ' 
     266         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)' 
    262267         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)' 
    263268         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5930 r6006  
    183183      ! -------------------- 
    184184      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf 
     185 
     186      tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere 
    185187      iif = jpreci                         ! ??? 
    186188      iil = nlci - jpreci + 1 
    … …  
    188190      ijl = nlcj - jprecj + 1 
    189191 
    190       tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns 
    191       tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
    192       tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows 
    193       tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
     192      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns 
     193      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
     194      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows 
     195      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
    194196 
    195197      ! north fold mask 
    … …  
    202204         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row 
    203205            DO ji = iif+1, iil-1 
    204                tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji)) 
     206               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 
    205207            END DO 
    206208         ENDIF 
    207209      ENDIF 
     210      
     211      tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 
     212 
    208213      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot 
    209214         tpol(     1    :jpiglo) = 0._wp 
    … …  
    225230         END DO 
    226231      END DO 
    227       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 
     232      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 
    228233      DO jj = 1, jpjm1 
    229234         DO ji = 1, fs_jpim1   ! vector loop 
    230             umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    231             vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     235            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     236            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    232237         END DO 
    233238         DO ji = 1, jpim1      ! NO vector opt. 
    234             fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
     239            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    235240               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    236241         END DO 
    237242      END DO 
    238       CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions 
    239       CALL lbc_lnk( vmask, 'V', 1._wp ) 
    240       CALL lbc_lnk( fmask, 'F', 1._wp ) 
    241       CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
    242       CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
    243       CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
     243      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions 
     244      CALL lbc_lnk( vmask  , 'V', 1._wp ) 
     245      CALL lbc_lnk( fmask  , 'F', 1._wp ) 
     246      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions 
     247      CALL lbc_lnk( ssvmask, 'V', 1._wp ) 
     248      CALL lbc_lnk( ssfmask, 'F', 1._wp ) 
    244249 
    245250      ! 3. Ocean/land mask at wu-, wv- and w points  

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

    r5836 r6006  
    2828CONTAINS 
    2929 
    30    SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid ) 
     30   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk ) 
    3131      !!---------------------------------------------------------------------- 
    3232      !!                    ***  ROUTINE dom_ngb  *** 
    … …  
    3939      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
    4040      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
     41      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used 
    4142      CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W' 
    4243      ! 
     44      INTEGER :: ik         ! working level 
    4345      INTEGER , DIMENSION(2) ::   iloc 
    4446      REAL(wp)               ::   zlon, zmini 
    … …  
    5153      ! 
    5254      zmask(:,:) = 0._wp 
     55      ik = 1 
     56      IF ( PRESENT(kkk) ) ik=kkk 
    5357      SELECT CASE( cdgrid ) 
    54       CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1) 
    55       CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1) 
    56       CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1) 
    57       CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1) 
     58      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik) 
     59      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik) 
     60      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik) 
     61      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik) 
    5862      END SELECT 
    5963 
    60       zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
    61       zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
    62       IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
    63       IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
     64      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN 
     65         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
     66         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
     67         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
     68         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
     69         zglam(:,:) = zglam(:,:) - zlon 
     70      ELSE 
     71         zglam(:,:) = zglam(:,:) - plon 
     72      END IF 
    6473 
    65       zglam(:,:) = zglam(:,:) - zlon 
    6674      zgphi(:,:) = zgphi(:,:) - plat 
    6775      zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5836 r6006  
    2424   USE in_out_manager  ! I/O manager 
    2525   USE iom             ! I/O manager library 
    26    USE restart         ! ocean restart 
     26   USE restart, ONLY : rst_read_open    ! ocean restart 
    2727   USE lib_mpp         ! distributed memory computing library 
    2828   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    192192         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    193193      END DO 
    194       hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 
    195       hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 
     194      hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) ) 
     195      hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) ) 
    196196 
    197197      ! Restoring frequencies for z_tilde coordinate 
    … …  
    325325         ! -------------------------------------------------- 
    326326         IF( ln_vvl_ztilde ) THEN 
    327             IF( kt .GT. nit000 ) THEN 
     327            IF( kt > nit000 ) THEN 
    328328               DO jk = 1, jpkm1 
    329329                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    … …  
    419419         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    420420         ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    421          IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN 
     421         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
    422422            IF( lk_mpp ) THEN 
    423423               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 
    … …  
    538538      END DO 
    539539      !                                        ! Inverse of the local depth 
    540       hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
    541       hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     540      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 
     541      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 
    542542 
    543543      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5836 r6006  
    379379      !!              - bathy : meter bathymetry (in meters) 
    380380      !!---------------------------------------------------------------------- 
    381       INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
     381      INTEGER  ::   ji, jj, jk            ! dummy loop indices 
    382382      INTEGER  ::   inum                      ! temporary logical unit 
    383383      INTEGER  ::   ierror                    ! error flag 
    … …  
    541541               CALL iom_close( inum ) 
    542542               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
     543 
     544               ! set grounded point to 0  
     545               ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 
     546               WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 
     547                  misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     548                  mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     549               END WHERE 
    543550            END IF 
    544551            !        
    … …  
    578585      !                        
    579586      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    580          ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    581          IF ( ln_isfcav ) THEN 
    582             WHERE (bathy == risfdep) 
    583                bathy   = 0.0_wp ; risfdep = 0.0_wp 
    584             END WHERE 
    585          END IF 
    586          ! end patch 
    587587         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    588588         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth 
    … …  
    835835   END SUBROUTINE zgr_bot_level 
    836836 
    837       SUBROUTINE zgr_top_level 
    838       !!---------------------------------------------------------------------- 
    839       !!                    ***  ROUTINE zgr_bot_level  *** 
     837   SUBROUTINE zgr_top_level 
     838      !!---------------------------------------------------------------------- 
     839      !!                    ***  ROUTINE zgr_top_level  *** 
    840840      !! 
    841841      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays) 
    … …  
    963963      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    964964      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    965       REAL(wp) ::   zmax             ! Maximum depth 
    966965      REAL(wp) ::   zdiff            ! temporary scalar 
    967       REAL(wp) ::   zrefdep          ! temporary scalar 
     966      REAL(wp) ::   zmax             ! temporary scalar 
    968967      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
    969968      !!--------------------------------------------------------------------- 
    … …  
    10041003      END DO 
    10051004 
    1006       IF ( ln_isfcav ) CALL zgr_isf 
    1007  
    10081005      ! Scale factors and depth at T- and W-points 
    10091006      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    … …  
    10131010         e3w_0  (:,:,jk) = e3w_1d  (jk) 
    10141011      END DO 
     1012       
     1013      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 
     1014      IF ( ln_isfcav ) CALL zgr_isf 
     1015 
     1016      ! Scale factors and depth at T- and W-points 
     1017      IF ( .NOT. ln_isfcav ) THEN 
     1018         DO jj = 1, jpj 
     1019            DO ji = 1, jpi 
     1020               ik = mbathy(ji,jj) 
     1021               IF( ik > 0 ) THEN               ! ocean point only 
     1022                  ! max ocean level case 
     1023                  IF( ik == jpkm1 ) THEN 
     1024                     zdepwp = bathy(ji,jj) 
     1025                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1026                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1027                     e3t_0(ji,jj,ik  ) = ze3tp 
     1028                     e3t_0(ji,jj,ik+1) = ze3tp 
     1029                     e3w_0(ji,jj,ik  ) = ze3wp 
     1030                     e3w_0(ji,jj,ik+1) = ze3tp 
     1031                     gdepw_0(ji,jj,ik+1) = zdepwp 
     1032                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1033                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1034                     ! 
     1035                  ELSE                         ! standard case 
     1036                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1037                     ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1038                     ENDIF 
     1039   !gm Bug?  check the gdepw_1d 
     1040                     !       ... on ik 
     1041                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1042                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1043                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1044                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1045                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1046                     e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1047                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1048                     !       ... on ik+1 
     1049                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1050                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1051                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1052                  ENDIF 
     1053               ENDIF 
     1054            END DO 
     1055         END DO 
     1056         ! 
     1057         it = 0 
     1058         DO jj = 1, jpj 
     1059            DO ji = 1, jpi 
     1060               ik = mbathy(ji,jj) 
     1061               IF( ik > 0 ) THEN               ! ocean point only 
     1062                  e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1063                  e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1064                  ! test 
     1065                  zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1066                  IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1067                     it = it + 1 
     1068                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1069                     WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1070                     WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1071                     WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1072                  ENDIF 
     1073               ENDIF 
     1074            END DO 
     1075         END DO 
     1076      END IF 
     1077      ! 
     1078      ! Scale factors and depth at U-, V-, UW and VW-points 
     1079      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1080         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1081         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1082         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1083         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1084      END DO 
     1085 
     1086      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1087         DO jj = 1, jpjm1 
     1088            DO ji = 1, fs_jpim1   ! vector opt. 
     1089               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1090               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1091               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1092               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1093            END DO 
     1094         END DO 
     1095      END DO 
     1096      IF ( ln_isfcav ) THEN 
     1097      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1098         DO jj = 2, jpjm1  
     1099            DO ji = 2, fs_jpim1   ! vector opt.  
     1100               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1101               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1102               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1103                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1104               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1105               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1106               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1107                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1108            END DO 
     1109         END DO 
     1110      END IF 
     1111 
     1112      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1113      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1114      ! 
     1115 
     1116      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1117         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1118         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1119         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1120         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1121      END DO 
     1122       
     1123      ! Scale factor at F-point 
     1124      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1125         e3f_0(:,:,jk) = e3t_1d(jk) 
     1126      END DO 
     1127      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1128         DO jj = 1, jpjm1 
     1129            DO ji = 1, fs_jpim1   ! vector opt. 
     1130               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1131            END DO 
     1132         END DO 
     1133      END DO 
     1134      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1135      ! 
     1136      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1137         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1138      END DO 
     1139!!gm  bug ? :  must be a do loop with mj0,mj1 
    10151140      !  
     1141      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1142      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1143      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1144      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1145      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1146 
     1147      ! Control of the sign 
     1148      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1149      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1150      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1151      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1152      
     1153      ! Compute gdep3w_0 (vertical sum of e3w) 
     1154      IF ( ln_isfcav ) THEN ! if cavity 
     1155         WHERE (misfdep == 0) misfdep = 1 
     1156         DO jj = 1,jpj 
     1157            DO ji = 1,jpi 
     1158               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1159               DO jk = 2, misfdep(ji,jj) 
     1160                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1161               END DO 
     1162               IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1163               DO jk = misfdep(ji,jj) + 1, jpk 
     1164                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1165               END DO 
     1166            END DO 
     1167         END DO 
     1168      ELSE ! no cavity 
     1169         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1170         DO jk = 2, jpk 
     1171            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1172         END DO 
     1173      END IF 
     1174      !                                               ! ================= ! 
     1175      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1176         !                                            ! ================= ! 
     1177         DO jj = 1,jpj 
     1178            DO ji = 1, jpi 
     1179               ik = MAX( mbathy(ji,jj), 1 ) 
     1180               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1181               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1182               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1183               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1184               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1185               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1186            END DO 
     1187         END DO 
     1188         WRITE(numout,*) 
     1189         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1190         WRITE(numout,*) 
     1191         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1192         WRITE(numout,*) 
     1193         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1194         WRITE(numout,*) 
     1195         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1196         WRITE(numout,*) 
     1197         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1198         WRITE(numout,*) 
     1199         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1200      ENDIF   
     1201      ! 
     1202      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1203      ! 
     1204      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1205      ! 
     1206   END SUBROUTINE zgr_zps 
     1207 
     1208   SUBROUTINE zgr_isf 
     1209      !!---------------------------------------------------------------------- 
     1210      !!                    ***  ROUTINE zgr_isf  *** 
     1211      !!    
     1212      !! ** Purpose :   check the bathymetry in levels 
     1213      !!    
     1214      !! ** Method  :   THe water column have to contained at least 2 cells 
     1215      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1216      !!                this criterion. 
     1217      !!    
     1218      !! ** Action  : - test compatibility between isfdraft and bathy  
     1219      !!              - bathy and isfdraft are modified 
     1220      !!---------------------------------------------------------------------- 
     1221      !!    
     1222      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices 
     1223      INTEGER  ::   ik, it               ! temporary integers 
     1224      INTEGER  ::   icompt, ibtest       ! (ISF) 
     1225      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF) 
     1226      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF) 
     1227      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t 
     1228      REAL(wp) ::   zmax             ! Maximum and minimum depth 
     1229      REAL(wp) ::   zbathydiff       ! isf temporary scalar 
     1230      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar 
     1231      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     1232      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
     1233      REAL(wp) ::   zdiff            ! temporary scalar 
     1234      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1235      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
     1236      !!--------------------------------------------------------------------- 
     1237      ! 
     1238      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1239      ! 
     1240      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
     1241      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1242 
     1243      ! (ISF) compute misfdep 
     1244      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     1245      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1246      END WHERE   
     1247 
     1248      ! Compute misfdep for ocean points (i.e. first wet level)  
     1249      ! find the first ocean level such that the first level thickness  
     1250      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
     1251      ! e3t_0 is the reference level thickness  
     1252      DO jk = 2, jpkm1  
     1253         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
     1254         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
     1255      END DO  
     1256      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 
     1257         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
     1258      END WHERE 
     1259 
     1260! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
     1261      icompt = 0  
     1262! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
     1263      DO jl = 1, 10      
     1264         ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 
     1265         WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 
     1266            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     1267            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     1268         END WHERE 
     1269         WHERE (mbathy(:,:) <= 0)  
     1270            misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
     1271            mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
     1272         ENDWHERE 
     1273         IF( lk_mpp ) THEN 
     1274            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1275            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1276            misfdep(:,:) = INT( zbathy(:,:) ) 
     1277 
     1278            CALL lbc_lnk( risfdep,'T', 1. ) 
     1279            CALL lbc_lnk( bathy,  'T', 1. ) 
     1280 
     1281            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1282            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1283            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1284         ENDIF 
     1285         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
     1286            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west  
     1287            misfdep(jpi,:) = misfdep(  2  ,:)  
     1288            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west 
     1289            mbathy(jpi,:)  = mbathy(  2  ,:) 
     1290         ENDIF 
     1291 
     1292         ! split last cell if possible (only where water column is 2 cell or less) 
     1293         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 
     1294         IF ( .NOT. ln_iscpl) THEN 
     1295            DO jk = jpkm1, 1, -1 
     1296               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1297               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
     1298                  mbathy(:,:) = jk 
     1299                  bathy(:,:)  = zmax 
     1300               END WHERE 
     1301            END DO 
     1302         END IF 
     1303  
     1304         ! split top cell if possible (only where water column is 2 cell or less) 
     1305         DO jk = 2, jpkm1 
     1306            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1307            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
     1308               misfdep(:,:) = jk 
     1309               risfdep(:,:) = zmax 
     1310            END WHERE 
     1311         END DO 
     1312 
     1313  
     1314 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     1315         DO jj = 1, jpj 
     1316            DO ji = 1, jpi 
     1317               ! find the minimum change option: 
     1318               ! test bathy 
     1319               IF (risfdep(ji,jj) > 1) THEN 
     1320                  IF ( .NOT. ln_iscpl ) THEN 
     1321                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
     1322                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1323                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
     1324                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1325  
     1326                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
     1327                        IF (zbathydiff <= zrisfdepdiff) THEN 
     1328                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1329                           mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1330                        ELSE 
     1331                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
     1332                           misfdep(ji,jj) = misfdep(ji,jj) - 1 
     1333                        END IF 
     1334                     ENDIF 
     1335                  ELSE 
     1336                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
     1337                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
     1338                        misfdep(ji,jj) = misfdep(ji,jj) - 1 
     1339                     END IF 
     1340                  END IF 
     1341               END IF 
     1342            END DO 
     1343         END DO 
     1344  
     1345         ! At least 2 levels for water thickness at T, U, and V point. 
     1346         DO jj = 1, jpj 
     1347            DO ji = 1, jpi 
     1348               ! find the minimum change option: 
     1349               ! test bathy 
     1350               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
     1351                  IF ( .NOT. ln_iscpl ) THEN  
     1352                     zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1353                         &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1354                     zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) &  
     1355                         &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1356                     IF (zbathydiff <= zrisfdepdiff) THEN 
     1357                        mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1358                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1359                     ELSE 
     1360                        misfdep(ji,jj)= misfdep(ji,jj) - 1 
     1361                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1362                     END IF 
     1363                  ELSE 
     1364                     misfdep(ji,jj)= misfdep(ji,jj) - 1 
     1365                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1366                  END IF 
     1367               ENDIF 
     1368            END DO 
     1369         END DO 
     1370  
     1371 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1372         DO jj = 1, jpjm1 
     1373            DO ji = 1, jpim1 
     1374               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
     1375                  IF ( .NOT. ln_iscpl ) THEN  
     1376                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1377                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
     1378                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 
     1379                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
     1380                     IF (zbathydiff <= zrisfdepdiff) THEN 
     1381                        mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1382                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
     1383                     ELSE 
     1384                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
     1385                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1386                     END IF 
     1387                  ELSE 
     1388                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
     1389                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1390                  END IF 
     1391               ENDIF 
     1392            END DO 
     1393         END DO 
     1394  
     1395         IF( lk_mpp ) THEN 
     1396            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1397            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1398            misfdep(:,:) = INT( zbathy(:,:) ) 
     1399 
     1400            CALL lbc_lnk( risfdep,'T', 1. ) 
     1401            CALL lbc_lnk( bathy,  'T', 1. ) 
     1402 
     1403            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1404            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1405            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1406         ENDIF 
     1407 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
     1408         DO jj = 1, jpjm1 
     1409            DO ji = 1, jpim1 
     1410               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 
     1411                  IF ( .NOT. ln_iscpl ) THEN  
     1412                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 
     1413                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
     1414                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) & 
     1415                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
     1416                     IF (zbathydiff <= zrisfdepdiff) THEN 
     1417                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
     1418                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     1419                     ELSE 
     1420                        misfdep(ji,jj)   = misfdep(ji,jj) - 1 
     1421                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
     1422                     END IF 
     1423                  ELSE 
     1424                     misfdep(ji,jj)   = misfdep(ji,jj) - 1 
     1425                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
     1426                  END IF 
     1427               ENDIF 
     1428            END DO 
     1429         END DO 
     1430  
     1431  
     1432         IF( lk_mpp ) THEN  
     1433            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1434            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1435            misfdep(:,:) = INT( zbathy(:,:) ) 
     1436 
     1437            CALL lbc_lnk( risfdep,'T', 1. ) 
     1438            CALL lbc_lnk( bathy,  'T', 1. ) 
     1439 
     1440            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1441            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1442            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1443         ENDIF  
     1444  
     1445 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1446         DO jj = 1, jpjm1 
     1447            DO ji = 1, jpim1 
     1448               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
     1449                  IF ( .NOT. ln_iscpl ) THEN  
     1450                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1451                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
     1452                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 
     1453                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
     1454                  IF (zbathydiff <= zrisfdepdiff) THEN 
     1455                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1456                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1457                  ELSE 
     1458                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
     1459                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
     1460                  END IF 
     1461                  ELSE 
     1462                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
     1463                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
     1464                  ENDIF 
     1465               ENDIF 
     1466            ENDDO 
     1467         ENDDO 
     1468  
     1469         IF( lk_mpp ) THEN  
     1470            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1471            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1472            misfdep(:,:) = INT( zbathy(:,:) ) 
     1473 
     1474            CALL lbc_lnk( risfdep,'T', 1. ) 
     1475            CALL lbc_lnk( bathy,  'T', 1. ) 
     1476 
     1477            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1478            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1479            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1480         ENDIF  
     1481  
     1482 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
     1483         DO jj = 1, jpjm1 
     1484            DO ji = 1, jpim1 
     1485               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 
     1486                  IF ( .NOT. ln_iscpl ) THEN  
     1487                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 
     1488                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
     1489                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) & 
     1490                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
     1491                     IF (zbathydiff <= zrisfdepdiff) THEN 
     1492                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
     1493                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
     1494                     ELSE 
     1495                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
     1496                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1497                     END IF 
     1498                  ELSE 
     1499                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
     1500                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1501                  ENDIF 
     1502               ENDIF 
     1503            ENDDO 
     1504         ENDDO 
     1505  
     1506         IF( lk_mpp ) THEN 
     1507            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1508            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1509            misfdep(:,:) = INT( zbathy(:,:) ) 
     1510 
     1511            CALL lbc_lnk( risfdep,'T', 1. ) 
     1512            CALL lbc_lnk( bathy,  'T', 1. ) 
     1513 
     1514            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1515            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1516            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1517         ENDIF 
     1518      END DO 
     1519      ! end dig bathy/ice shelf to be compatible 
     1520      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
     1521      DO jl = 1,20 
     1522  
     1523 ! remove single point "bay" on isf coast line in the ice shelf draft' 
     1524         DO jk = 2, jpk 
     1525            WHERE (misfdep==0) misfdep=jpk 
     1526            zmask=0._wp 
     1527            WHERE (misfdep <= jk) zmask=1 
     1528            DO jj = 2, jpjm1 
     1529               DO ji = 2, jpim1 
     1530                  IF (misfdep(ji,jj) == jk) THEN 
     1531                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1532                     IF (ibtest <= 1) THEN 
     1533                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
     1534                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 
     1535                     END IF 
     1536                  END IF 
     1537               END DO 
     1538            END DO 
     1539         END DO 
     1540         WHERE (misfdep==jpk) 
     1541             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
     1542         END WHERE 
     1543         IF( lk_mpp ) THEN 
     1544            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1545            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1546            misfdep(:,:) = INT( zbathy(:,:) ) 
     1547 
     1548            CALL lbc_lnk( risfdep,'T', 1. ) 
     1549            CALL lbc_lnk( bathy,  'T', 1. ) 
     1550 
     1551            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1552            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1553            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1554         ENDIF 
     1555  
     1556 ! remove single point "bay" on bathy coast line beneath an ice shelf' 
     1557         DO jk = jpk,1,-1 
     1558            zmask=0._wp 
     1559            WHERE (mbathy >= jk ) zmask=1 
     1560            DO jj = 2, jpjm1 
     1561               DO ji = 2, jpim1 
     1562                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 
     1563                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1564                     IF (ibtest <= 1) THEN 
     1565                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
     1566                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 
     1567                     END IF 
     1568                  END IF 
     1569               END DO 
     1570            END DO 
     1571         END DO 
     1572         WHERE (mbathy==0) 
     1573             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
     1574         END WHERE 
     1575         IF( lk_mpp ) THEN 
     1576            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
     1577            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1578            misfdep(:,:) = INT( zbathy(:,:) ) 
     1579 
     1580            CALL lbc_lnk( risfdep,'T', 1. ) 
     1581            CALL lbc_lnk( bathy,  'T', 1. ) 
     1582 
     1583            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1584            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1585            mbathy(:,:)  = INT( zbathy(:,:) ) 
     1586         ENDIF 
     1587  
     1588 ! fill hole in ice shelf 
     1589         zmisfdep = misfdep 
     1590         zrisfdep = risfdep 
     1591         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 
     1592         DO jj = 2, jpjm1 
     1593            DO ji = 2, jpim1 
     1594               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
     1595               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
     1596               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk 
     1597               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk 
     1598               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk 
     1599               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk 
     1600               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1601               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 
     1602                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
     1603               END IF 
     1604               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 
     1605                  misfdep(ji,jj) = ibtest 
     1606                  risfdep(ji,jj) = gdepw_1d(ibtest) 
     1607               ENDIF 
     1608            ENDDO 
     1609         ENDDO 
     1610  
     1611         IF( lk_mpp ) THEN  
     1612            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1613            CALL lbc_lnk( zbathy,  'T', 1. )  
     1614            misfdep(:,:) = INT( zbathy(:,:) )  
     1615 
     1616            CALL lbc_lnk( risfdep, 'T', 1. )  
     1617            CALL lbc_lnk( bathy,   'T', 1. ) 
     1618 
     1619            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1620            CALL lbc_lnk( zbathy,  'T', 1. ) 
     1621            mbathy(:,:) = INT( zbathy(:,:) ) 
     1622         ENDIF  
     1623 ! 
     1624 !! fill hole in bathymetry 
     1625         zmbathy (:,:)=mbathy (:,:) 
     1626         DO jj = 2, jpjm1 
     1627            DO ji = 2, jpim1 
     1628               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
     1629               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
     1630               IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0 
     1631               IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0 
     1632               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
     1633               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
     1634               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
     1635               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 
     1636                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
     1637               END IF 
     1638               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 
     1639                  mbathy(ji,jj) = ibtest 
     1640                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
     1641               ENDIF 
     1642            END DO 
     1643         END DO 
     1644         IF( lk_mpp ) THEN  
     1645            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1646            CALL lbc_lnk( zbathy,  'T', 1. )  
     1647            misfdep(:,:) = INT( zbathy(:,:) )  
     1648 
     1649            CALL lbc_lnk( risfdep, 'T', 1. )  
     1650            CALL lbc_lnk( bathy,   'T', 1. ) 
     1651 
     1652            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1653            CALL lbc_lnk( zbathy,  'T', 1. ) 
     1654            mbathy(:,:) = INT( zbathy(:,:) ) 
     1655         ENDIF  
     1656 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
     1657         DO jj = 1, jpjm1 
     1658            DO ji = 1, jpim1 
     1659               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
     1660                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
     1661               END IF 
     1662            END DO 
     1663         END DO 
     1664         IF( lk_mpp ) THEN  
     1665            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1666            CALL lbc_lnk( zbathy,  'T', 1. )  
     1667            misfdep(:,:) = INT( zbathy(:,:) )  
     1668 
     1669            CALL lbc_lnk( risfdep, 'T', 1. )  
     1670            CALL lbc_lnk( bathy,   'T', 1. ) 
     1671 
     1672            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1673            CALL lbc_lnk( zbathy,  'T', 1. ) 
     1674            mbathy(:,:) = INT( zbathy(:,:) ) 
     1675         ENDIF  
     1676 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
     1677         DO jj = 1, jpjm1 
     1678            DO ji = 1, jpim1 
     1679               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
     1680                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
     1681               END IF 
     1682            END DO 
     1683         END DO 
     1684         IF( lk_mpp ) THEN  
     1685            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1686            CALL lbc_lnk( zbathy, 'T', 1. )  
     1687            misfdep(:,:) = INT( zbathy(:,:) )  
     1688 
     1689            CALL lbc_lnk( risfdep,'T', 1. )  
     1690            CALL lbc_lnk( bathy,  'T', 1. ) 
     1691 
     1692            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1693            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1694            mbathy(:,:) = INT( zbathy(:,:) ) 
     1695         ENDIF  
     1696 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
     1697         DO jj = 1, jpjm1 
     1698            DO ji = 1, jpi 
     1699               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
     1700                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
     1701               END IF 
     1702            END DO 
     1703         END DO 
     1704         IF( lk_mpp ) THEN  
     1705            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1706            CALL lbc_lnk( zbathy, 'T', 1. )  
     1707            misfdep(:,:) = INT( zbathy(:,:) )  
     1708 
     1709            CALL lbc_lnk( risfdep,'T', 1. )  
     1710            CALL lbc_lnk( bathy,  'T', 1. ) 
     1711 
     1712            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1713            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1714            mbathy(:,:) = INT( zbathy(:,:) ) 
     1715         ENDIF  
     1716 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
     1717         DO jj = 1, jpjm1 
     1718            DO ji = 1, jpi 
     1719               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
     1720                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 
     1721               END IF 
     1722            END DO 
     1723         END DO 
     1724         IF( lk_mpp ) THEN  
     1725            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1726            CALL lbc_lnk( zbathy, 'T', 1. )  
     1727            misfdep(:,:) = INT( zbathy(:,:) )  
     1728 
     1729            CALL lbc_lnk( risfdep,'T', 1. )  
     1730            CALL lbc_lnk( bathy,  'T', 1. ) 
     1731 
     1732            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1733            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1734            mbathy(:,:) = INT( zbathy(:,:) ) 
     1735         ENDIF  
     1736 ! if not compatible after all check, mask T 
     1737         DO jj = 1, jpj 
     1738            DO ji = 1, jpi 
     1739               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
     1740                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
     1741               END IF 
     1742            END DO 
     1743         END DO 
     1744  
     1745         WHERE (mbathy(:,:) == 1) 
     1746            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
     1747         END WHERE 
     1748      END DO  
     1749! end check compatibility ice shelf/bathy 
     1750      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1751      IF( icompt == 0 ) THEN  
     1752         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
     1753      ELSE  
     1754         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
     1755      ENDIF  
     1756 
     1757      ! compute scale factor and depth at T- and W- points 
    10161758      DO jj = 1, jpj 
    10171759         DO ji = 1, jpi 
    … …  
    10351777                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    10361778                  ENDIF 
     1779      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    10371780!gm Bug?  check the gdepw_1d 
    10381781                  !       ... on ik 
    … …  
    10401783                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    10411784                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1042                   e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1043                      &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1044                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1045                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1785                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  ) 
     1786                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1) 
    10461787                  !       ... on ik+1 
    10471788                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    10481789                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1049                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    10501790               ENDIF 
    10511791            ENDIF 
    … …  
    10731813      END DO 
    10741814      ! 
    1075       IF ( ln_isfcav ) THEN 
    10761815      ! (ISF) Definition of e3t, u, v, w for ISF case 
    1077          DO jj = 1, jpj  
    1078             DO ji = 1, jpi  
    1079                ik = misfdep(ji,jj)  
    1080                IF( ik > 1 ) THEN               ! ice shelf point only  
    1081                   IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1082                   gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1816      DO jj = 1, jpj  
     1817         DO ji = 1, jpi  
     1818            ik = misfdep(ji,jj)  
     1819            IF( ik > 1 ) THEN               ! ice shelf point only  
     1820               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1821               gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    10831822!gm Bug?  check the gdepw_0  
    1084                !       ... on ik  
    1085                   gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1086                      &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1087                      &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1088                   e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1089                   e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1090  
    1091                   IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1092                      e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1093                   ENDIF  
    1094                !       ... on ik / ik-1  
    1095                   e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1096                   e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1823            !       ... on ik  
     1824               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1825                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1826                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1827               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1828               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1829 
     1830               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1831                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1832               ENDIF  
     1833            !       ... on ik / ik-1  
     1834               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik)  
     1835               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    10971836! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1098                   gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1837               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1838            ENDIF  
     1839         END DO  
     1840      END DO  
     1841    
     1842      it = 0  
     1843      DO jj = 1, jpj  
     1844         DO ji = 1, jpi  
     1845            ik = misfdep(ji,jj)  
     1846            IF( ik > 1 ) THEN               ! ice shelf point only  
     1847               e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1848               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1849            ! test  
     1850               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1851               IF( zdiff <= 0. .AND. lwp ) THEN   
     1852                  it = it + 1  
     1853                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1854                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1855                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1856                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    10991857               ENDIF  
    1100             END DO  
     1858            ENDIF  
    11011859         END DO  
    1102       !  
    1103          it = 0  
    1104          DO jj = 1, jpj  
    1105             DO ji = 1, jpi  
    1106                ik = misfdep(ji,jj)  
    1107                IF( ik > 1 ) THEN               ! ice shelf point only  
    1108                   e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1109                   e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1110                ! test  
    1111                   zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1112                   IF( zdiff <= 0. .AND. lwp ) THEN   
    1113                      it = it + 1  
    1114                      WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1115                      WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1116                      WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1117                      WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1118                   ENDIF  
    1119                ENDIF  
    1120             END DO  
    1121          END DO  
    1122       END IF 
    1123       ! END (ISF) 
    1124  
    1125       ! Scale factors and depth at U-, V-, UW and VW-points 
    1126       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1127          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1128          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1129          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1130          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1131       END DO 
    1132       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1133          DO jj = 1, jpjm1 
    1134             DO ji = 1, fs_jpim1   ! vector opt. 
    1135                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1136                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1137                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1138                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1139             END DO 
    1140          END DO 
    1141       END DO 
    1142       IF ( ln_isfcav ) THEN 
    1143       ! (ISF) define e3uw (adapted for 2 cells in the water column) 
    1144          DO jj = 2, jpjm1  
    1145             DO ji = 2, fs_jpim1   ! vector opt.  
    1146                ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
    1147                ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
    1148                IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
    1149                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
    1150                ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
    1151                ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
    1152                IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
    1153                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
    1154             END DO 
    1155          END DO 
    1156       END IF 
    1157  
    1158       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1159       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1160       ! 
    1161       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1162          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1163          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1164          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1165          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1166       END DO 
    1167        
    1168       ! Scale factor at F-point 
    1169       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1170          e3f_0(:,:,jk) = e3t_1d(jk) 
    1171       END DO 
    1172       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1173          DO jj = 1, jpjm1 
    1174             DO ji = 1, fs_jpim1   ! vector opt. 
    1175                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1176             END DO 
    1177          END DO 
    1178       END DO 
    1179       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1180       ! 
    1181       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1182          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1183       END DO 
    1184 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1185       !  
    1186       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1187       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1188       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1189       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1190       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1191  
    1192       ! Control of the sign 
    1193       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1194       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1195       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1196       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1197       
    1198       ! Compute gdep3w_0 (vertical sum of e3w) 
    1199       IF ( ln_isfcav ) THEN ! if cavity 
    1200          WHERE (misfdep == 0) misfdep = 1 
    1201          DO jj = 1,jpj 
    1202             DO ji = 1,jpi 
    1203                gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1204                DO jk = 2, misfdep(ji,jj) 
    1205                   gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1206                END DO 
    1207                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)) 
    1208                DO jk = misfdep(ji,jj) + 1, jpk 
    1209                   gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1210                END DO 
    1211             END DO 
    1212          END DO 
    1213       ELSE ! no cavity 
    1214          gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
    1215          DO jk = 2, jpk 
    1216             gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
    1217          END DO 
    1218       END IF 
    1219       !                                               ! ================= ! 
    1220       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1221          !                                            ! ================= ! 
    1222          DO jj = 1,jpj 
    1223             DO ji = 1, jpi 
    1224                ik = MAX( mbathy(ji,jj), 1 ) 
    1225                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1226                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1227                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1228                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1229                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1230                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1231             END DO 
    1232          END DO 
    1233          WRITE(numout,*) 
    1234          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1235          WRITE(numout,*) 
    1236          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1237          WRITE(numout,*) 
    1238          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1239          WRITE(numout,*) 
    1240          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1241          WRITE(numout,*) 
    1242          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1243          WRITE(numout,*) 
    1244          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1245       ENDIF   
    1246       ! 
    1247       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    1248       ! 
    1249       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1250       ! 
    1251    END SUBROUTINE zgr_zps 
    1252  
    1253    SUBROUTINE zgr_isf 
    1254       !!---------------------------------------------------------------------- 
    1255       !!                    ***  ROUTINE zgr_isf  *** 
    1256       !!    
    1257       !! ** Purpose :   check the bathymetry in levels 
    1258       !!    
    1259       !! ** Method  :   THe water column have to contained at least 2 cells 
    1260       !!                Bathymetry and isfdraft are modified (dig/close) to respect 
    1261       !!                this criterion. 
    1262       !!                  
    1263       !!    
    1264       !! ** Action  : - test compatibility between isfdraft and bathy  
    1265       !!              - bathy and isfdraft are modified 
    1266       !!---------------------------------------------------------------------- 
    1267       !!    
    1268       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    1269       INTEGER  ::   ik, it           ! temporary integers 
    1270       INTEGER  ::   id, jd, nprocd 
    1271       INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
    1272       LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    1273       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    1274       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    1275       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
    1276       REAL(wp) ::   zdiff            ! temporary scalar 
    1277       REAL(wp) ::   zrefdep          ! temporary scalar 
    1278       REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
    1279       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
    1280       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    1281       !!--------------------------------------------------------------------- 
    1282       ! 
    1283       IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
    1284       ! 
    1285       CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    1286       CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
    1287  
    1288  
    1289       ! (ISF) compute misfdep 
    1290       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    1291       ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    1292       END WHERE   
    1293  
    1294       ! Compute misfdep for ocean points (i.e. first wet level)  
    1295       ! find the first ocean level such that the first level thickness  
    1296       ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
    1297       ! e3t_0 is the reference level thickness  
    1298       DO jk = 2, jpkm1  
    1299          zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
    1300          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    13011860      END DO  
    1302       WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 
    1303          risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    1304       END WHERE 
    1305   
    1306 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
    1307       icompt = 0  
    1308 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    1309       DO jl = 1, 10      
    1310          WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
    1311             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    1312             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
    1313          END WHERE 
    1314          WHERE (mbathy(:,:) <= 0)  
    1315             misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
    1316             mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
    1317          ENDWHERE 
    1318          IF( lk_mpp ) THEN 
    1319             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1320             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1321             misfdep(:,:) = INT( zbathy(:,:) ) 
    1322             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1323             CALL lbc_lnk( bathy, 'T', 1. ) 
    1324             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1325             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1326             mbathy(:,:) = INT( zbathy(:,:) ) 
    1327          ENDIF 
    1328          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
    1329             misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west  
    1330             misfdep(jpi,:) = misfdep(  2  ,:)  
    1331          ENDIF 
    1332  
    1333          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    1334             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    1335             mbathy(jpi,:) = mbathy(  2  ,:) 
    1336          ENDIF 
    1337  
    1338          ! split last cell if possible (only where water column is 2 cell or less) 
    1339          DO jk = jpkm1, 1, -1 
    1340             zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1341             WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1342                mbathy(:,:) = jk 
    1343                bathy(:,:)  = zmax 
    1344             END WHERE 
    1345          END DO 
    1346   
    1347          ! split top cell if possible (only where water column is 2 cell or less) 
    1348          DO jk = 2, jpkm1 
    1349             zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1350             WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
    1351                misfdep(:,:) = jk 
    1352                risfdep(:,:) = zmax 
    1353             END WHERE 
    1354          END DO 
    1355  
    1356   
    1357  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    1358          DO jj = 1, jpj 
    1359             DO ji = 1, jpi 
    1360                ! find the minimum change option: 
    1361                ! test bathy 
    1362                IF (risfdep(ji,jj) .GT. 1) THEN 
    1363                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    1364                  &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1365                zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    1366                  &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1367   
    1368                   IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1369                      IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1370                         bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1371                         mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1372                      ELSE 
    1373                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1374                         misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1375                      END IF 
    1376                   END IF 
    1377                END IF 
    1378             END DO 
    1379          END DO 
    1380   
    1381           ! At least 2 levels for water thickness at T, U, and V point. 
    1382          DO jj = 1, jpj 
    1383             DO ji = 1, jpi 
    1384                ! find the minimum change option: 
    1385                ! test bathy 
    1386                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1387                   zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)& 
    1388                     & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1389                   zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  & 
    1390                     & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1391                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1392                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1393                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1394                   ELSE 
    1395                      misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1396                      risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1397                   END IF 
    1398                ENDIF 
    1399             END DO 
    1400          END DO 
    1401   
    1402  ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
    1403          DO jj = 1, jpjm1 
    1404             DO ji = 1, jpim1 
    1405                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1406                   zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   & 
    1407                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1408                   zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  & 
    1409                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1410                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1411                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1412                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) & 
    1413                    &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1414                   ELSE 
    1415                      misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1416                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 
    1417                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1418                   END IF 
    1419                ENDIF 
    1420             END DO 
    1421          END DO 
    1422   
    1423          IF( lk_mpp ) THEN 
    1424             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1425             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1426             misfdep(:,:) = INT( zbathy(:,:) ) 
    1427             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1428             CALL lbc_lnk( bathy, 'T', 1. ) 
    1429             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1430             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1431             mbathy(:,:) = INT( zbathy(:,:) ) 
    1432          ENDIF 
    1433  ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
    1434          DO jj = 1, jpjm1 
    1435             DO ji = 1, jpim1 
    1436                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1437                   zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 
    1438                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1439                   zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  & 
    1440                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1441                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1442                      mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1443                      bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1444                   ELSE 
    1445                      misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1446                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1447                   END IF 
    1448                ENDIF 
    1449             END DO 
    1450          END DO 
    1451   
    1452   
    1453          IF( lk_mpp ) THEN  
    1454             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1455             CALL lbc_lnk( zbathy, 'T', 1. )  
    1456             misfdep(:,:) = INT( zbathy(:,:) )  
    1457             CALL lbc_lnk( risfdep, 'T', 1. )  
    1458             CALL lbc_lnk( bathy, 'T', 1. ) 
    1459             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1460             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1461             mbathy(:,:) = INT( zbathy(:,:) ) 
    1462          ENDIF  
    1463   
    1464  ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
    1465          DO jj = 1, jpjm1 
    1466             DO ji = 1, jpim1 
    1467                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1468                   zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 
    1469                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1470                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 
    1471                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1472                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1473                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1474                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1475                   ELSE 
    1476                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1477                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1478                   END IF 
    1479                ENDIF 
    1480             ENDDO 
    1481          ENDDO 
    1482   
    1483          IF( lk_mpp ) THEN  
    1484             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1485             CALL lbc_lnk( zbathy, 'T', 1. )  
    1486             misfdep(:,:) = INT( zbathy(:,:) )  
    1487             CALL lbc_lnk( risfdep, 'T', 1. )  
    1488             CALL lbc_lnk( bathy, 'T', 1. ) 
    1489             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1490             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1491             mbathy(:,:) = INT( zbathy(:,:) ) 
    1492          ENDIF  
    1493   
    1494  ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
    1495          DO jj = 1, jpjm1 
    1496             DO ji = 1, jpim1 
    1497                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1498                   zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 
    1499                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1500                   zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  & 
    1501                       &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1502                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1503                      mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1504                      bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  & 
    1505                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1506                   ELSE 
    1507                      misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1508                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) & 
    1509                       &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1510                   END IF 
    1511                ENDIF 
    1512             ENDDO 
    1513          ENDDO 
    1514   
    1515          IF( lk_mpp ) THEN 
    1516             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1517             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1518             misfdep(:,:) = INT( zbathy(:,:) ) 
    1519             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1520             CALL lbc_lnk( bathy, 'T', 1. ) 
    1521             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1522             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1523             mbathy(:,:) = INT( zbathy(:,:) ) 
    1524          ENDIF 
    1525       END DO 
    1526       ! end dig bathy/ice shelf to be compatible 
    1527       ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
    1528       DO jl = 1,20 
    1529   
    1530  ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1531          DO jk = 2, jpk 
    1532             WHERE (misfdep==0) misfdep=jpk 
    1533             zmask=0 
    1534             WHERE (misfdep .LE. jk) zmask=1 
    1535             DO jj = 2, jpjm1 
    1536                DO ji = 2, jpim1 
    1537                   IF (misfdep(ji,jj) .EQ. jk) THEN 
    1538                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1539                      IF (ibtest .LE. 1) THEN 
    1540                         risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
    1541                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 
    1542                      END IF 
    1543                   END IF 
    1544                END DO 
    1545             END DO 
    1546          END DO 
    1547          WHERE (misfdep==jpk) 
    1548              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
    1549          END WHERE 
    1550          IF( lk_mpp ) THEN 
    1551             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1552             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1553             misfdep(:,:) = INT( zbathy(:,:) ) 
    1554             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1555             CALL lbc_lnk( bathy, 'T', 1. ) 
    1556             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1557             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1558             mbathy(:,:) = INT( zbathy(:,:) ) 
    1559          ENDIF 
    1560   
    1561  ! remove single point "bay" on bathy coast line beneath an ice shelf' 
    1562          DO jk = jpk,1,-1 
    1563             zmask=0 
    1564             WHERE (mbathy .GE. jk ) zmask=1 
    1565             DO jj = 2, jpjm1 
    1566                DO ji = 2, jpim1 
    1567                   IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 
    1568                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1569                      IF (ibtest .LE. 1) THEN 
    1570                         bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
    1571                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 
    1572                      END IF 
    1573                   END IF 
    1574                END DO 
    1575             END DO 
    1576          END DO 
    1577          WHERE (mbathy==0) 
    1578              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
    1579          END WHERE 
    1580          IF( lk_mpp ) THEN 
    1581             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
    1582             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1583             misfdep(:,:) = INT( zbathy(:,:) ) 
    1584             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1585             CALL lbc_lnk( bathy, 'T', 1. ) 
    1586             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1587             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1588             mbathy(:,:) = INT( zbathy(:,:) ) 
    1589          ENDIF 
    1590   
    1591  ! fill hole in ice shelf 
    1592          zmisfdep = misfdep 
    1593          zrisfdep = risfdep 
    1594          WHERE (zmisfdep .LE. 1) zmisfdep=jpk 
    1595          DO jj = 2, jpjm1 
    1596             DO ji = 2, jpim1 
    1597                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    1598                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1599                IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
    1600                IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
    1601                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
    1602                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
    1603                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1604                IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
    1605                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
    1606                END IF 
    1607                IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 
    1608                   misfdep(ji,jj) = ibtest 
    1609                   risfdep(ji,jj) = gdepw_1d(ibtest) 
    1610                ENDIF 
    1611             ENDDO 
    1612          ENDDO 
    1613   
    1614          IF( lk_mpp ) THEN  
    1615             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1616             CALL lbc_lnk( zbathy, 'T', 1. )  
    1617             misfdep(:,:) = INT( zbathy(:,:) )  
    1618             CALL lbc_lnk( risfdep, 'T', 1. )  
    1619             CALL lbc_lnk( bathy, 'T', 1. ) 
    1620             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1621             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1622             mbathy(:,:) = INT( zbathy(:,:) ) 
    1623          ENDIF  
    1624  ! 
    1625  !! fill hole in bathymetry 
    1626          zmbathy (:,:)=mbathy (:,:) 
    1627          DO jj = 2, jpjm1 
    1628             DO ji = 2, jpim1 
    1629                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    1630                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1631                IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
    1632                IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    1633                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
    1634                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    1635                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1636                IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    1637                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    1638                END IF 
    1639                IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 
    1640                   mbathy(ji,jj) = ibtest 
    1641                   bathy(ji,jj)  = gdepw_1d(ibtest+1)  
    1642                ENDIF 
    1643             END DO 
    1644          END DO 
    1645          IF( lk_mpp ) THEN  
    1646             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1647             CALL lbc_lnk( zbathy, 'T', 1. )  
    1648             misfdep(:,:) = INT( zbathy(:,:) )  
    1649             CALL lbc_lnk( risfdep, 'T', 1. )  
    1650             CALL lbc_lnk( bathy, 'T', 1. ) 
    1651             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1652             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1653             mbathy(:,:) = INT( zbathy(:,:) ) 
    1654          ENDIF  
    1655  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1656          DO jj = 1, jpjm1 
    1657             DO ji = 1, jpim1 
    1658                IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
    1659                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1660                END IF 
    1661             END DO 
    1662          END DO 
    1663          IF( lk_mpp ) THEN  
    1664             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1665             CALL lbc_lnk( zbathy, 'T', 1. )  
    1666             misfdep(:,:) = INT( zbathy(:,:) )  
    1667             CALL lbc_lnk( risfdep, 'T', 1. )  
    1668             CALL lbc_lnk( bathy, 'T', 1. ) 
    1669             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1670             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1671             mbathy(:,:) = INT( zbathy(:,:) ) 
    1672          ENDIF  
    1673  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1674          DO jj = 1, jpjm1 
    1675             DO ji = 1, jpim1 
    1676                IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
    1677                   mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
    1678                END IF 
    1679             END DO 
    1680          END DO 
    1681          IF( lk_mpp ) THEN  
    1682             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1683             CALL lbc_lnk( zbathy, 'T', 1. )  
    1684             misfdep(:,:) = INT( zbathy(:,:) )  
    1685             CALL lbc_lnk( risfdep, 'T', 1. )  
    1686             CALL lbc_lnk( bathy, 'T', 1. ) 
    1687             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1688             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1689             mbathy(:,:) = INT( zbathy(:,:) ) 
    1690          ENDIF  
    1691  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1692          DO jj = 1, jpjm1 
    1693             DO ji = 1, jpi 
    1694                IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
    1695                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1696                END IF 
    1697             END DO 
    1698          END DO 
    1699          IF( lk_mpp ) THEN  
    1700             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1701             CALL lbc_lnk( zbathy, 'T', 1. )  
    1702             misfdep(:,:) = INT( zbathy(:,:) )  
    1703             CALL lbc_lnk( risfdep, 'T', 1. )  
    1704             CALL lbc_lnk( bathy, 'T', 1. ) 
    1705             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1706             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1707             mbathy(:,:) = INT( zbathy(:,:) ) 
    1708          ENDIF  
    1709  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1710          DO jj = 1, jpjm1 
    1711             DO ji = 1, jpi 
    1712                IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
    1713                   mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ; 
    1714                END IF 
    1715             END DO 
    1716          END DO 
    1717          IF( lk_mpp ) THEN  
    1718             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1719             CALL lbc_lnk( zbathy, 'T', 1. )  
    1720             misfdep(:,:) = INT( zbathy(:,:) )  
    1721             CALL lbc_lnk( risfdep, 'T', 1. )  
    1722             CALL lbc_lnk( bathy, 'T', 1. ) 
    1723             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1724             CALL lbc_lnk( zbathy, 'T', 1. ) 
    1725             mbathy(:,:) = INT( zbathy(:,:) ) 
    1726          ENDIF  
    1727  ! if not compatible after all check, mask T 
    1728          DO jj = 1, jpj 
    1729             DO ji = 1, jpi 
    1730                IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
    1731                   misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
    1732                END IF 
    1733             END DO 
    1734          END DO 
    1735   
    1736          WHERE (mbathy(:,:) == 1) 
    1737             mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
    1738          END WHERE 
    1739       END DO  
    1740 ! end check compatibility ice shelf/bathy 
    1741       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1742       WHERE (misfdep(:,:) <= 5) 
    1743          misfdep = 1; risfdep = 0.0_wp; 
    1744       END WHERE 
    1745  
    1746       IF( icompt == 0 ) THEN  
    1747          IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
    1748       ELSE  
    1749          IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
    1750       ENDIF  
    17511861 
    17521862      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    … …  
    17551865      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
    17561866       
    1757    END SUBROUTINE 
     1867   END SUBROUTINE zgr_isf 
    17581868 
    17591869   SUBROUTINE zgr_sco 
    … …  
    20842194      fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    20852195      ! 
    2086       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    2087       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    2088       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    2089       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    2090       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    2091       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    2092       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2196      where (e3t_0   (:,:,:) == 0.0)  e3t_0(:,:,:)  = 1.0_wp 
     2197      where (e3u_0   (:,:,:) == 0.0)  e3u_0(:,:,:)  = 1.0_wp 
     2198      where (e3v_0   (:,:,:) == 0.0)  e3v_0(:,:,:)  = 1.0_wp 
     2199      where (e3f_0   (:,:,:) == 0.0)  e3f_0(:,:,:)  = 1.0_wp 
     2200      where (e3w_0   (:,:,:) == 0.0)  e3w_0(:,:,:)  = 1.0_wp 
     2201      where (e3uw_0  (:,:,:) == 0.0)  e3uw_0(:,:,:) = 1.0_wp 
     2202      where (e3vw_0  (:,:,:) == 0.0)  e3vw_0(:,:,:) = 1.0_wp 
    20932203 
    20942204#if defined key_agrif 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5930 r6006  
    3535   USE dtauvd          ! data: U & V current             (dta_uvd routine) 
    3636   USE domvvl          ! varying vertical mesh 
     37   USE iscplrst        ! ice sheet coupling 
    3738   ! 
    3839   USE in_out_manager  ! I/O manager 
    … …  
    8586      IF( ln_rstart ) THEN                    ! Restart from a file 
    8687         !                                    ! ------------------- 
    87          CALL rst_read                           ! Read the restart file 
    88          CALL day_init                           ! model calendar (using both namelist and restart infos) 
     88         CALL rst_read                        ! Read the restart file 
     89         IF (ln_iscpl)       CALL iscpl_stp   ! extraloate restart to wet and dry 
     90         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    8991      ELSE 
    9092         !                                    ! Start from rest 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90

    r5836 r6006  
    2323   USE sbcrnf          ! river runoff  
    2424   USE sbcisf          ! ice shelf 
     25   USE iscplhsb        ! ice sheet / ocean coupling 
     26   USE iscplini        ! ice sheet / ocean coupling 
    2527   ! 
    2628   USE in_out_manager  ! I/O manager 
    … …  
    9395      IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    9496      ! 
     97      IF( ln_iscpl  .AND. ln_hsb     )   CALL iscpl_div( hdivn )   !==  ice sheet  ==!   (update hdivn field) 
     98      ! 
    9599      CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change) 
    96100      ! 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r5930 r6006  
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    97       INTEGER  ::   iku, ikv     ! local integers 
    9897      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars 
    9998      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    … …  
    320319            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    321320         END DO 
    322          hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 
    323          hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 
     321         hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     322         hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    324323      ENDIF 
    325324      ! 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r6005 r6006  
    142142                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn      ) 
    143143                     CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop      ) 
     144 
     145                  ! extra variable needed for the ice sheet coupling 
     146                  IF ( ln_iscpl ) THEN  
     147                     CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask     ) ! need to extrapolate T/S 
     148                     CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask     ) ! need to correct barotropic velocity 
     149                     CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask     ) ! need to correct barotropic velocity 
     150                     CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask    ) ! need to correct barotropic velocity 
     151                     CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )   ! need to compute temperature correction 
     152                     CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) )   ! need to compute bt conservation 
     153                     CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) )   ! need to compute bt conservation 
     154                     CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute extrapolation if vvl 
     155                  END IF 
    144156      ENDIF 
    145157       
    146158      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    )   
    147        
     159 
    148160      IF( kt == nitrst ) THEN 
    149161         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
    … …  
    209221      REAL(wp) ::   zrdt, zrdttra1 
    210222      INTEGER  ::   jk 
    211       LOGICAL  ::   llok 
    212223      !!---------------------------------------------------------------------- 
    213224 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90

    r5429 r6006  
    1717   !!---------------------------------------------------------------------- 
    1818   !!   lbc_lnk      : generic interface for mpp_lnk_3d and mpp_lnk_2d routines defined in lib_mpp 
     19   !!   lbc_sum      : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d routines defined in lib_mpp 
    1920   !!   lbc_lnk_e    : generic interface for mpp_lnk_2d_e routine defined in lib_mpp 
    2021   !!   lbc_bdy_lnk  : generic interface for mpp_lnk_bdy_2d and mpp_lnk_bdy_3d routines defined in lib_mpp 
    … …  
    3132   END INTERFACE 
    3233 
     34!JMM interface not defined if not key_mpp_mpi : likely do not compile without this CPP key !!!! 
     35   INTERFACE lbc_sum 
     36      MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 
     37   END INTERFACE 
     38 
    3339   INTERFACE lbc_bdy_lnk 
    3440      MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d 
    … …  
    4551   PUBLIC lbc_lnk       ! ocean lateral boundary conditions 
    4652   PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions 
     53   PUBLIC lbc_sum 
    4754   PUBLIC lbc_lnk_e 
    4855   PUBLIC lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    … …  
    5966   !!   Default option                              shared memory computing 
    6067   !!---------------------------------------------------------------------- 
    61    !!   lbc_lnk      : generic interface for lbc_lnk_3d and lbc_lnk_2d 
    62    !!   lbc_lnk_3d   : set the lateral boundary condition on a 3D variable on ocean mesh 
    63    !!   lbc_lnk_2d   : set the lateral boundary condition on a 2D variable on ocean mesh 
    64    !!   lbc_bdy_lnk  : set the lateral BDY boundary condition 
     68   !!   lbc_sum       : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d  
     69   !!   lbc_lnk_sum_3d: compute sum over the halos on a 3D variable on ocean mesh 
     70   !!   lbc_lnk_sum_3d: compute sum over the halos on a 2D variable on ocean mesh 
     71   !!   lbc_lnk       : generic interface for lbc_lnk_3d and lbc_lnk_2d 
     72   !!   lbc_lnk_3d    : set the lateral boundary condition on a 3D variable on ocean mesh 
     73   !!   lbc_lnk_2d    : set the lateral boundary condition on a 2D variable on ocean mesh 
     74   !!   lbc_bdy_lnk   : set the lateral BDY boundary condition 
    6575   !!---------------------------------------------------------------------- 
    6676   USE oce             ! ocean dynamics and tracers    
    … …  
    7484   INTERFACE lbc_lnk 
    7585      MODULE PROCEDURE lbc_lnk_3d_gather, lbc_lnk_3d, lbc_lnk_2d 
     86   END INTERFACE 
     87 
     88   INTERFACE lbc_sum 
     89      MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 
    7690   END INTERFACE 
    7791 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r5836 r6006  
    7272   PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e 
    7373   PUBLIC   mpp_lnk_2d_9  
     74   PUBLIC   mpp_lnk_sum_3d, mpp_lnk_sum_2d 
    7475   PUBLIC   mppscatter, mppgather 
    7576   PUBLIC   mpp_ini_ice, mpp_ini_znl 
    … …  
    14021403   END SUBROUTINE mpp_lnk_2d_e 
    14031404 
     1405   SUBROUTINE mpp_lnk_sum_3d( ptab, cd_type, psgn, cd_mpp, pval ) 
     1406      !!---------------------------------------------------------------------- 
     1407      !!                  ***  routine mpp_lnk_sum_3d  *** 
     1408      !! 
     1409      !! ** Purpose :   Message passing manadgement (sum the overlap region) 
     1410      !! 
     1411      !! ** Method  :   Use mppsend and mpprecv function for passing mask 
     1412      !!      between processors following neighboring subdomains. 
     1413      !!            domain parameters 
     1414      !!                    nlci   : first dimension of the local subdomain 
     1415      !!                    nlcj   : second dimension of the local subdomain 
     1416      !!                    nbondi : mark for "east-west local boundary" 
     1417      !!                    nbondj : mark for "north-south local boundary" 
     1418      !!                    noea   : number for local neighboring processors 
     1419      !!                    nowe   : number for local neighboring processors 
     1420      !!                    noso   : number for local neighboring processors 
     1421      !!                    nono   : number for local neighboring processors 
     1422      !! 
     1423      !! ** Action  :   ptab with update value at its periphery 
     1424      !! 
     1425      !!---------------------------------------------------------------------- 
     1426      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied 
     1427      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points 
     1428      !                                                             ! = T , U , V , F , W points 
     1429      REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary 
     1430      !                                                             ! =  1. , the sign is kept 
     1431      CHARACTER(len=3), OPTIONAL      , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only 
     1432      REAL(wp)        , OPTIONAL      , INTENT(in   ) ::   pval     ! background value (used at closed boundaries) 
     1433      !! 
     1434      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
     1435      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers 
     1436      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend 
     1437      REAL(wp) ::   zland 
     1438      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend 
     1439      ! 
     1440      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north 
     1441      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east 
     1442 
     1443      !!---------------------------------------------------------------------- 
     1444       
     1445      ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   & 
     1446         &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  ) 
     1447 
     1448      ! 
     1449      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value 
     1450      ELSE                         ;   zland = 0.e0      ! zero by default 
     1451      ENDIF 
     1452 
     1453      ! 1. standard boundary treatment 
     1454      ! ------------------------------ 
     1455      ! 2. East and west directions exchange 
     1456      ! ------------------------------------ 
     1457      ! we play with the neigbours AND the row number because of the periodicity 
     1458      ! 
     1459      SELECT CASE ( nbondi )      ! Read lateral conditions 
     1460      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case) 
     1461      iihom = nlci-jpreci 
     1462         DO jl = 1, jpreci 
     1463            zt3ew(:,jl,:,1) = ptab(jl      ,:,:) ; ptab(jl      ,:,:) = 0.0_wp 
     1464            zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp  
     1465         END DO 
     1466      END SELECT 
     1467      ! 
     1468      !                           ! Migrations 
     1469      imigr = jpreci * jpj * jpk 
     1470      ! 
     1471      SELECT CASE ( nbondi ) 
     1472      CASE ( -1 ) 
     1473         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 ) 
     1474         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 
     1475         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1476      CASE ( 0 ) 
     1477         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 
     1478         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 ) 
     1479         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea ) 
     1480         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 
     1481         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1482         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 
     1483      CASE ( 1 ) 
     1484         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 ) 
     1485         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe ) 
     1486         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1487      END SELECT 
     1488      ! 
     1489      !                           ! Write lateral conditions 
     1490      iihom = nlci-nreci 
     1491      ! 
     1492      SELECT CASE ( nbondi ) 
     1493      CASE ( -1 ) 
     1494         DO jl = 1, jpreci 
     1495            ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2) 
     1496         END DO 
     1497      CASE ( 0 ) 
     1498         DO jl = 1, jpreci 
     1499            ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 
     1500            ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2) 
     1501         END DO 
     1502      CASE ( 1 ) 
     1503         DO jl = 1, jpreci 
     1504            ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2) 
     1505         END DO 
     1506      END SELECT 
     1507 
     1508 
     1509      ! 3. North and south directions 
     1510      ! ----------------------------- 
     1511      ! always closed : we play only with the neigbours 
     1512      ! 
     1513      IF( nbondj /= 2 ) THEN      ! Read lateral conditions 
     1514         ijhom = nlcj-jprecj 
     1515         DO jl = 1, jprecj 
     1516            zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp 
     1517            zt3ns(:,jl,:,1) = ptab(:,jl      ,:) ; ptab(:,jl      ,:) = 0.0_wp 
     1518         END DO 
     1519      ENDIF 
     1520      ! 
     1521      !                           ! Migrations 
     1522      imigr = jprecj * jpi * jpk 
     1523      ! 
     1524      SELECT CASE ( nbondj ) 
     1525      CASE ( -1 ) 
     1526         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 ) 
     1527         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 
     1528         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1529      CASE ( 0 ) 
     1530         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 
     1531         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 ) 
     1532         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono ) 
     1533         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 
     1534         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1535         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 
     1536      CASE ( 1 ) 
     1537         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 ) 
     1538         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso ) 
     1539         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 
     1540      END SELECT 
     1541      ! 
     1542      !                           ! Write lateral conditions 
     1543      ijhom = nlcj-nrecj 
     1544      ! 
     1545      SELECT CASE ( nbondj ) 
     1546      CASE ( -1 ) 
     1547         DO jl = 1, jprecj 
     1548            ptab(:,ijhom+jl,:) = ptab(:,ijhom+jl,:) + zt3ns(:,jl,:,2) 
     1549         END DO 
     1550      CASE ( 0 ) 
     1551         DO jl = 1, jprecj 
     1552            ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl,:,2) 
     1553            ptab(:,ijhom +jl,:) = ptab(:,ijhom +jl,:) + zt3ns(:,jl,:,2) 
     1554         END DO 
     1555      CASE ( 1 ) 
     1556         DO jl = 1, jprecj 
     1557            ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl   ,:,2) 
     1558         END DO 
     1559      END SELECT 
     1560 
     1561 
     1562      ! 4. north fold treatment 
     1563      ! ----------------------- 
     1564      ! 
     1565      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 
     1566         ! 
     1567         SELECT CASE ( jpni ) 
     1568         CASE ( 1 )     ;   CALL lbc_nfd      ( ptab, cd_type, psgn )   ! only 1 northern proc, no mpp 
     1569         CASE DEFAULT   ;   CALL mpp_lbc_north( ptab, cd_type, psgn )   ! for all northern procs. 
     1570         END SELECT 
     1571         ! 
     1572      ENDIF 
     1573      ! 
     1574      DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we ) 
     1575      ! 
     1576   END SUBROUTINE mpp_lnk_sum_3d 
     1577 
     1578   SUBROUTINE mpp_lnk_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 
     1579      !!---------------------------------------------------------------------- 
     1580      !!                  ***  routine mpp_lnk_sum_2d  *** 
     1581      !! 
     1582      !! ** Purpose :   Message passing manadgement for 2d array (sum the overlap region) 
     1583      !! 
     1584      !! ** Method  :   Use mppsend and mpprecv function for passing mask 
     1585      !!      between processors following neighboring subdomains. 
     1586      !!            domain parameters 
     1587      !!                    nlci   : first dimension of the local subdomain 
     1588      !!                    nlcj   : second dimension of the local subdomain 
     1589      !!                    nbondi : mark for "east-west local boundary" 
     1590      !!                    nbondj : mark for "north-south local boundary" 
     1591      !!                    noea   : number for local neighboring processors 
     1592      !!                    nowe   : number for local neighboring processors 
     1593      !!                    noso   : number for local neighboring processors 
     1594      !!                    nono   : number for local neighboring processors 
     1595      !! 
     1596      !!---------------------------------------------------------------------- 
     1597      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d     ! 2D array on which the boundary condition is applied 
     1598      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points 
     1599      !                                                         ! = T , U , V , F , W and I points 
     1600      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary 
     1601      !                                                         ! =  1. , the sign is kept 
     1602      CHARACTER(len=3), OPTIONAL  , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only 
     1603      REAL(wp)        , OPTIONAL  , INTENT(in   ) ::   pval     ! background value (used at closed boundaries) 
     1604      !! 
     1605      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1606      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers 
     1607      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend 
     1608      REAL(wp) ::   zland 
     1609      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend 
     1610      ! 
     1611      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north 
     1612      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east 
     1613 
     1614      !!---------------------------------------------------------------------- 
     1615 
     1616      ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2),  & 
     1617         &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   ) 
     1618 
     1619      ! 
     1620      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value 
     1621      ELSE                         ;   zland = 0.e0      ! zero by default 
     1622      ENDIF 
     1623 
     1624      ! 1. standard boundary treatment 
     1625      ! ------------------------------ 
     1626      ! 2. East and west directions exchange 
     1627      ! ------------------------------------ 
     1628      ! we play with the neigbours AND the row number because of the periodicity 
     1629      ! 
     1630      SELECT CASE ( nbondi )      ! Read lateral conditions 
     1631      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case) 
     1632         iihom = nlci - jpreci 
     1633         DO jl = 1, jpreci 
     1634            zt2ew(:,jl,1) = pt2d(jl       ,:) ; pt2d(jl       ,:) = 0.0_wp 
     1635            zt2we(:,jl,1) = pt2d(iihom +jl,:) ; pt2d(iihom +jl,:) = 0.0_wp 
     1636         END DO 
     1637      END SELECT 
     1638      ! 
     1639      !                           ! Migrations 
     1640      imigr = jpreci * jpj 
     1641      ! 
     1642      SELECT CASE ( nbondi ) 
     1643      CASE ( -1 ) 
     1644         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 ) 
     1645         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 
     1646         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1647      CASE ( 0 ) 
     1648         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 
     1649         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 ) 
     1650         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea ) 
     1651         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 
     1652         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1653         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 
     1654      CASE ( 1 ) 
     1655         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 ) 
     1656         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe ) 
     1657         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1658      END SELECT 
     1659      ! 
     1660      !                           ! Write lateral conditions 
     1661      iihom = nlci-nreci 
     1662      ! 
     1663      SELECT CASE ( nbondi ) 
     1664      CASE ( -1 ) 
     1665         DO jl = 1, jpreci 
     1666            pt2d(iihom+jl,:) = pt2d(iihom+jl,:) + zt2ew(:,jl,2) 
     1667         END DO 
     1668      CASE ( 0 ) 
     1669         DO jl = 1, jpreci 
     1670            pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 
     1671            pt2d(iihom +jl,:) = pt2d(iihom +jl,:) + zt2ew(:,jl,2) 
     1672         END DO 
     1673      CASE ( 1 ) 
     1674         DO jl = 1, jpreci 
     1675            pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2) 
     1676         END DO 
     1677      END SELECT 
     1678 
     1679 
     1680      ! 3. North and south directions 
     1681      ! ----------------------------- 
     1682      ! always closed : we play only with the neigbours 
     1683      ! 
     1684      IF( nbondj /= 2 ) THEN      ! Read lateral conditions 
     1685         ijhom = nlcj - jprecj 
     1686         DO jl = 1, jprecj 
     1687            zt2sn(:,jl,1) = pt2d(:,ijhom +jl) ; pt2d(:,ijhom +jl) = 0.0_wp 
     1688            zt2ns(:,jl,1) = pt2d(:,jl       ) ; pt2d(:,jl       ) = 0.0_wp 
     1689         END DO 
     1690      ENDIF 
     1691      ! 
     1692      !                           ! Migrations 
     1693      imigr = jprecj * jpi 
     1694      ! 
     1695      SELECT CASE ( nbondj ) 
     1696      CASE ( -1 ) 
     1697         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 ) 
     1698         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 
     1699         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1700      CASE ( 0 ) 
     1701         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 
     1702         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 ) 
     1703         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono ) 
     1704         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 
     1705         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1706         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 
     1707      CASE ( 1 ) 
     1708         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 ) 
     1709         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso ) 
     1710         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 
     1711      END SELECT 
     1712      ! 
     1713      !                           ! Write lateral conditions 
     1714      ijhom = nlcj-nrecj 
     1715      ! 
     1716      SELECT CASE ( nbondj ) 
     1717      CASE ( -1 ) 
     1718         DO jl = 1, jprecj 
     1719            pt2d(:,ijhom+jl) = pt2d(:,ijhom+jl) + zt2ns(:,jl,2) 
     1720         END DO 
     1721      CASE ( 0 ) 
     1722         DO jl = 1, jprecj 
     1723            pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 
     1724            pt2d(:,ijhom +jl) = pt2d(:,ijhom +jl) + zt2ns(:,jl,2) 
     1725         END DO 
     1726      CASE ( 1 ) 
     1727         DO jl = 1, jprecj 
     1728            pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2) 
     1729         END DO 
     1730      END SELECT 
     1731 
     1732 
     1733      ! 4. north fold treatment 
     1734      ! ----------------------- 
     1735      ! 
     1736      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 
     1737         ! 
     1738         SELECT CASE ( jpni ) 
     1739         CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d, cd_type, psgn )   ! only 1 northern proc, no mpp 
     1740         CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d, cd_type, psgn )   ! for all northern procs. 
     1741         END SELECT 
     1742         ! 
     1743      ENDIF 
     1744      ! 
     1745      DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we ) 
     1746      ! 
     1747   END SUBROUTINE mpp_lnk_sum_2d 
    14041748 
    14051749   SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

    r5130 r6006  
    136136 
    137137      imask(:,:)=1 
    138       WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0 
     138      WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0 
    139139 
    140140      !  1. Dimension arrays for subdomains 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r5836 r6006  
    119119         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf  
    120120         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
     121         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
     122         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
    121123         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
    122124         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 
    … …  
    230232            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    231233            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 
    232             CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U') 
    233             CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V') 
     234            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
     235            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
    234236            ! iom print 
    235237            CALL iom_put('ttbl',ttbl(:,:)) 
    236238            CALL iom_put('stbl',stbl(:,:)) 
    237             CALL iom_put('utbl',utbl(:,:)) 
    238             CALL iom_put('vtbl',vtbl(:,:)) 
     239            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
     240            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
    239241            ! compute fwf and heat flux 
    240242            CALL sbc_isf_cav (kt) 
    … …  
    356358    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
    357359     ! 
    358  
    359     ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 
    360360    DO ji = 1, jpi 
    361361       DO jj = 1, jpj 
    … …  
    366366    ! 3. -----------the average temperature between 200m and 600m --------------------- 
    367367             DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    368              ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
    369              ! after verif with UNESCO, wrong sign in BG eq. 2 
    370368             ! Calculate freezing temperature 
    371                 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04  
    372                 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress)  
     369                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik))  
    373370                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp 
    374371             ENDDO 
    … …  
    418415      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b 
    419416      REAL(wp) ::   zgammat, zgammas 
    420       REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==! 
     417      REAL(wp) ::   zeps = 1.e-20_wp        !==   Local constant initialization   ==! 
    421418      INTEGER  ::   ji, jj     ! dummy loop indices 
    422419      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
    … …  
    503500                     zeps1=rcp*rau0*zgammat 
    504501                     zeps2=lfusisf*rau0*zgammas 
    505                      zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj) 
     502                     zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 
    506503                     zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    507504                     zeps6=zeps4-zti(ji,jj) 
    508505                     zeps7=zeps4-tsurf 
    509506                     zaqe=zlamb1 * (zeps1 + zeps3) 
    510                      zaqer=0.5/zaqe 
     507                     zaqer=0.5/MIN(zaqe,-zeps) 
    511508                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 
    512509                     zcqe=zeps2*stbl(ji,jj) 
    … …  
    517514                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    518515   
    519 ! zfwflx is upward water flux 
    520                      zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz ) 
     516! zfwflx is upward water flux (MAX usefull if kappa = 0 
     517                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 
    521518! zhtflx is upward heat flux (out of ocean) 
    522519! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5836 r6006  
    100100      ! 
    101101      CALL wrk_alloc( jpi,jpj,jpk,jpts,   zts_dta ) 
    102       ! 
    103102      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    104103         CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts )  

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r5643 r6006  
    2222   USE sbcrnf          ! River runoff   
    2323   USE sbcisf          ! Ice shelf    
     24   USE iscplini        ! Ice sheet coupling 
    2425   USE traqsr          ! solar radiation penetration 
    2526   USE trd_oce         ! trends: ocean variables 
    … …  
    117118      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices   
    118119      INTEGER  ::   ikt, ikb  
    119       INTEGER  ::   nk_isf 
    120120      REAL(wp) ::   zfact, z1_e3t, zdep 
    121121      REAL(wp) ::   zalpha, zhk 
    … …  
    219219      ! 
    220220      IF( nn_isf > 0 ) THEN 
    221          zfact = 0.5e0 
     221         zfact = 0.5_wp 
    222222         DO jj = 2, jpj 
    223223            DO ji = fs_2, fs_jpim1 
    … …  
    230230               ! sign - because fwf sign of evapo (rnf sign of precip) 
    231231               DO jk = ikt, ikb - 1 
    232                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    233232               ! compute trend 
    234233                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    … …  
    239238    
    240239               ! level partially include in ice shelf boundary layer  
    241                ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    242240               ! compute trend 
    243241               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    … …  
    279277      ENDIF 
    280278  
    281       IF( l_trdtra )   THEN                      ! send trends for further diagnostics 
     279      !---------------------------------------- 
     280      !        Ice Sheet coupling imbalance correction to have conservation 
     281      !---------------------------------------- 
     282      ! 
     283      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff  
     284         DO jk = 1,jpk 
     285            DO jj = 2, jpj  
     286               DO ji = fs_2, fs_jpim1 
     287                  zdep = 1._wp / fse3t_n(ji,jj,jk)  
     288                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem)                       & 
     289                      &                                         * zdep 
     290                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal)                       & 
     291                      &                                         * zdep   
     292               END DO   
     293            END DO   
     294         END DO 
     295      ENDIF 
     296 
     297      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    282298         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    283299         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90

    r4161 r6006  
    2424   PRIVATE 
    2525 
    26    PUBLIC   glob_sum   ! used in many places 
    27    PUBLIC   DDPDD      ! also used in closea module 
     26   PUBLIC   glob_sum      ! used in many places (masked with tmask_i) 
     27   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie omly over the halos) 
     28   PUBLIC   DDPDD         ! also used in closea module 
    2829   PUBLIC   glob_min, glob_max 
    2930#if defined key_nosignedzero 
    … …  
    3435      MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, & 
    3536         &             glob_sum_2d_a, glob_sum_3d_a 
     37   END INTERFACE 
     38   INTERFACE glob_sum_full 
     39      MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d 
    3640   END INTERFACE 
    3741   INTERFACE glob_min 
    … …  
    156160      ! 
    157161   END FUNCTION glob_sum_3d_a 
     162 
     163   FUNCTION glob_sum_full_2d( ptab ) 
     164      !!---------------------------------------------------------------------- 
     165      !!                  ***  FUNCTION  glob_sum_full_2d *** 
     166      !! 
     167      !! ** Purpose : perform a sum in calling DDPDD routine (nomask) 
     168      !!---------------------------------------------------------------------- 
     169      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     170      REAL(wp)                             ::   glob_sum_full_2d   ! global sum 
     171      !! 
     172      !!----------------------------------------------------------------------- 
     173      ! 
     174      glob_sum_full_2d = SUM( ptab(:,:) * tmask_h(:,:) ) 
     175      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_2d ) 
     176      ! 
     177   END FUNCTION glob_sum_full_2d 
     178 
     179   FUNCTION glob_sum_full_3d( ptab ) 
     180      !!---------------------------------------------------------------------- 
     181      !!                  ***  FUNCTION  glob_sum_full_3d *** 
     182      !! 
     183      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine (nomask) 
     184      !!---------------------------------------------------------------------- 
     185      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     186      REAL(wp)                               ::   glob_sum_full_3d   ! global sum 
     187      !! 
     188      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     189      INTEGER    ::   ijpk ! local variables: size of ptab 
     190      !!----------------------------------------------------------------------- 
     191      ! 
     192      ijpk = SIZE(ptab,3) 
     193      ! 
     194      glob_sum_full_3d = 0.e0 
     195      DO jk = 1, ijpk 
     196         glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) * tmask_h(:,:) ) 
     197      END DO 
     198      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_3d ) 
     199      ! 
     200   END FUNCTION glob_sum_full_3d 
     201 
    158202 
    159203#else   
    … …  
    314358   END FUNCTION glob_sum_3d_a    
    315359 
     360   FUNCTION glob_sum_full_2d( ptab ) 
     361      !!---------------------------------------------------------------------- 
     362      !!                  ***  FUNCTION  glob_sum_full_2d *** 
     363      !! 
     364      !! ** Purpose : perform a sum in calling DDPDD routine 
     365      !!---------------------------------------------------------------------- 
     366      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     367      REAL(wp)                             ::   glob_sum_full_2d   ! global sum (nomask) 
     368      !! 
     369      COMPLEX(wp)::   ctmp 
     370      REAL(wp)   ::   ztmp 
     371      INTEGER    ::   ji, jj   ! dummy loop indices 
     372      !!----------------------------------------------------------------------- 
     373      ! 
     374      ztmp = 0.e0 
     375      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     376      DO jj = 1, jpj 
     377         DO ji =1, jpi 
     378         ztmp =  ptab(ji,jj) * tmask_h(ji,jj)  
     379         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     380         END DO 
     381      END DO 
     382      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain 
     383      glob_sum_full_2d = REAL(ctmp,wp) 
     384      ! 
     385   END FUNCTION glob_sum_full_2d 
     386 
     387   FUNCTION glob_sum_full_3d( ptab ) 
     388      !!---------------------------------------------------------------------- 
     389      !!                  ***  FUNCTION  glob_sum_full_3d *** 
     390      !! 
     391      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 
     392      !!---------------------------------------------------------------------- 
     393      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     394      REAL(wp)                               ::   glob_sum_full_3d   ! global sum (nomask) 
     395      !! 
     396      COMPLEX(wp)::   ctmp 
     397      REAL(wp)   ::   ztmp 
     398      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     399      INTEGER    ::   ijpk ! local variables: size of ptab 
     400      !!----------------------------------------------------------------------- 
     401      ! 
     402      ijpk = SIZE(ptab,3) 
     403      ! 
     404      ztmp = 0.e0 
     405      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     406      DO jk = 1, ijpk 
     407         DO jj = 1, jpj 
     408            DO ji =1, jpi 
     409            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj) 
     410            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     411            END DO 
     412         END DO 
     413      END DO 
     414      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain 
     415      glob_sum_full_3d = REAL(ctmp,wp) 
     416      ! 
     417   END FUNCTION glob_sum_full_3d 
     418 
     419 
     420 
    316421#endif 
    317422 
Note: See TracChangeset for help on using the changeset viewer.