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 5619 for branches/NERC/dev_r5589_is_oce_cpl – NEMO

Ignore:
Timestamp:
2015-07-20T19:43:15+02:00 (9 years ago)
Author:
mathiot
Message:

ocean/ice sheet coupling: initial commit

Location:
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM
Files:
1 added
25 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/field_def.xml

    r5517 r5619  
    193193 
    194194         <!-- * variable related to ice shelf forcing * --> 
    195          <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
    196          <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     /> 
    197          <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      /> 
    198          <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      /> 
    199          <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="1e-3"     /> 
    200          <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degC"     /> 
     195         <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  /> 
     196         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     /> 
     197         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "                 unit="m/s"      /> 
     198         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "                 unit="m/s"      /> 
     199    <field id="stbl"         long_name="salinity in the Losh tbl                    "                 unit="PSU"      /> 
     200    <field id="ttbl"         long_name="temperature in the Losh tbl                 "                 unit="C"        /> 
     201    <field id="utbl"         long_name="zonal current in the Losh tbl at T point    "                 unit="m/s"      /> 
     202    <field id="vtbl"         long_name="merid current in the Losh tbl at T point    "                 unit="m/s"      /> 
     203    <field id="thermald"     long_name="thermal driving of ice shelf melting        "                 unit="C"        /> 
     204    <field id="tfrz"         long_name="top freezing point (used to compute melt)   "                 unit="C"        /> 
     205    <field id="tinsitu"      long_name="top insitu temperature (used to cmpt melt)  "                 unit="C"        /> 
     206    <field id="ustar"        long_name="ustar at T point used in ice shelf melting  "                 unit="m/s"      /> 
    201207 
    202208         <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5577 r5619  
    4040   cn_ocerst_out = "restart"   !  suffix of ocean restart name (output) 
    4141   cn_ocerst_outdir = "."      !  directory in which to write output ocean restarts 
     42   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4243   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4344   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
     
    225226!!   namsbc_rnf      river runoffs 
    226227!!   namsbc_isf      ice shelf melting/freezing 
     228!!   namsbc_iscpl    coupling option between land ice model and ocean 
    227229!!   namsbc_apr      Atmospheric Pressure 
    228230!!   namsbc_ssr      sea surface restoring term (for T and/or S) 
     
    469471/ 
    470472!----------------------------------------------------------------------- 
     473&namsbc_iscpl  !   land ice / ocean coupling option                      
     474!----------------------------------------------------------------------- 
     475   rn_fiscpl = 43800    ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 
     476   ln_hfb    = .false.  ! activate conservation module (conservation exact after a time of rn_fiscpl) 
     477/ 
     478!----------------------------------------------------------------------- 
    471479&namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk 
    472480!----------------------------------------------------------------------- 
     
    802810   rn_smsh          =     1.    !  Smagorinsky diffusivity: = 0 - use only sheer 
    803811   rn_aht_m         =  2000.    !  upper limit or stability criteria for lateral eddy diffusivity (m2/s) 
     812/ 
     813!----------------------------------------------------------------------- 
     814&namtra_dmpfile    !   tracer: T & S newtonian damping 
     815!----------------------------------------------------------------------- 
     816!          !  file name                            ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     817!          !                                       !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
     818   sn_dmpt  = 'resto',         -1        ,'Tinit' ,    .true.    , .true. , 'yearly'   , ''       ,   ''    ,    '' 
     819   sn_dmps  = 'resto',         -1        ,'Sinit' ,    .true.    , .true. , 'yearly'   , ''       ,   ''    ,    '' 
    804820/ 
    805821!----------------------------------------------------------------------- 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r5586 r5619  
    197197                     ! Elevation 
    198198                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)         
    199                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj) 
    200                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj) 
     199                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
     200                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
    201201                  END DO 
    202202               END DO 
     
    326326               X1= ana_amp(ji,jj,jh,1) 
    327327               X2=-ana_amp(ji,jj,jh,2) 
    328                out_u(ji,jj,       jh) = X1 * umask_i(ji,jj) 
    329                out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj) 
     328               out_u(ji,jj,       jh) = X1 * ssumask(ji,jj) 
     329               out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 
    330330            ENDDO 
    331331         ENDDO 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5120 r5619  
    165165      DO jk = 1, jpkm1 
    166166         ! volume variation (calculated with scale factors) 
    167          zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) & 
    168             &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
     167         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * ( tmask(:,:,jk) & 
     168            &                           * fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
    169169         ! heat content variation 
    170          zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) &  
    171             &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) ) 
     170         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * ( tmask(:,:,jk) &  
     171            &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) ) 
    172172         ! salt content variation 
    173          zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   & 
    174             &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) ) 
     173         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * ( tmask(:,:,jk)   & 
     174            &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) ) 
    175175      ENDDO 
    176176 
     
    200200!      ENDIF 
    201201!!gm end 
    202  
     202       
     203      IF (lwp) PRINT *, 'ISCPL CONS HEAT ', kt, zdiff_hc / zvol_tot, zdiff_sc / zvol_tot 
     204      IF (lwp) PRINT *, 'ISCPL CONS VOL  ', kt, zdiff_v1 * 1.e-9, zdiff_v2 * 1.e-9 
    203205 
    204206      IF( lk_vvl ) THEN 
     
    217219        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J)  
    218220        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)   
     221        CALL iom_put( 'bgvolssh' , (zdiff_v1+zdiff_v2) * 1.e-9    )   ! volume ssh variation (km3)   
    220222        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    221223        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)  
     
    277279          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh 
    278280          DO jk = 1, jpk 
    279              e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors 
    280              hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content 
    281              sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content 
     281             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)    * tmask(:,:,jk)                    ! initial vertical scale factors 
     282             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content 
     283             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content 
    282284          END DO 
    283285          frc_v = 0._wp                                           ! volume       trend due to forcing 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5123 r5619  
    4343   INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    4444   INTEGER , PUBLIC ::   nn_euler        !: =0 start with forward time step or not (=1) 
     45   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    4546   LOGICAL , PUBLIC ::   ln_crs          !: Apply grid coarsening to dynamical model output or online passive tracers 
    4647 
     
    252253   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level 
    253254   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level 
    254    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters) 
    255    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 
    256    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(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function 
    257258 
    258259   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF) 
    259260   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF) 
    260261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                       (ISF) 
    261    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask                   !: surface domain T-point mask  
    262  
     262 
     263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts 
    263264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    264265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
     
    389390         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
    390391 
    391       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
    392          &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
    393          &     bmask(jpi,jpj)   ,                                                       & 
    394          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
     392      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       & 
     393         &     tmask_i(jpi,jpj) ,                                                          &  
     394         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 
     395         &     bmask(jpi,jpj)   ,                                                          & 
     396         &     mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
    395397 
    396398! (ISF) Allocation of basic array    
    397399      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     & 
    398400         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
    399          &     mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 
     401         &     mikf(jpi,jpj), STAT=ierr(10) ) 
    400402 
    401403      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r5619  
    111111      END DO 
    112112      !                                        ! Inverse of the local depth 
    113       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
    114       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     113      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 
     114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 
    115115 
    116116                             CALL dom_stp      ! time step 
    117       IF( nmsh /= 0     )   CALL dom_wri      ! Create a domain file 
     117      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )   CALL dom_wri      ! Create a domain file 
    118118      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control 
    119119      ! 
     
    138138         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
    139139         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140          &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
     140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler, ln_iscpl 
    141141      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   & 
    142142         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  & 
     
    192192         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    193193         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     194         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl 
    194195      ENDIF 
    195196 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5552 r5619  
    264264         END DO 
    265265      END DO 
    266       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 
     266      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 
    267267      DO jj = 1, jpjm1 
    268268         DO ji = 1, fs_jpim1   ! vector loop 
    269             umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    270             vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     269            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     270            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    271271         END DO 
    272272         DO ji = 1, jpim1      ! NO vector opt. 
    273             fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
     273            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    274274               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    275275         END DO 
    276276      END DO 
    277       CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions 
    278       CALL lbc_lnk( vmask, 'V', 1._wp ) 
    279       CALL lbc_lnk( fmask, 'F', 1._wp ) 
    280       CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
    281       CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
    282       CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
     277      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions 
     278      CALL lbc_lnk( vmask  , 'V', 1._wp ) 
     279      CALL lbc_lnk( fmask  , 'F', 1._wp ) 
     280      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions 
     281      CALL lbc_lnk( ssvmask, 'V', 1._wp ) 
     282      CALL lbc_lnk( ssfmask, 'F', 1._wp ) 
    283283 
    284284      ! 3. Ocean/land mask at wu-, wv- and w points  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90

    r3294 r5619  
    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  *** 
     
    4040      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
    4141      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
     42      INTEGER         , INTENT(in   ), OPTIONAL :: kkk 
    4243      CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W' 
    4344      ! 
    4445      INTEGER , DIMENSION(2) ::   iloc 
     46      INTEGER :: jk 
    4547      REAL(wp)               ::   zlon, zmini 
    4648      REAL(wp), POINTER, DIMENSION(:,:) ::  zglam, zgphi, zmask, zdist 
     
    5254      ! 
    5355      zmask(:,:) = 0._wp 
     56      jk = 1 
     57      IF (PRESENT(kkk)) jk=kkk 
    5458      SELECT CASE( cdgrid ) 
    55       CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1) 
    56       CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1) 
    57       CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1) 
    58       CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1) 
     59      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,jk) 
     60      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,jk) 
     61      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,jk) 
     62      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,jk) 
    5963      END SELECT 
    6064 
    61       zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
    62       zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
    63       IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
    64       IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
    65  
    66       zglam(:,:) = zglam(:,:) - zlon 
     65      IF (jphgr_msh .NE. 2 .AND. jphgr_msh .NE. 3) THEN 
     66         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
     67         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
     68         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
     69         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
     70         zglam(:,:) = zglam(:,:) - zlon 
     71      ELSE 
     72         zglam(:,:) = zglam(:,:) - plon 
     73      END IF 
     74! 
    6775      zgphi(:,:) = zgphi(:,:) - plat 
    6876      zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5506 r5619  
    2828   USE in_out_manager  ! I/O manager 
    2929   USE iom             ! I/O manager library 
    30    USE restart         ! ocean restart 
     30   USE restart  , ONLY : rst_read_open    ! ocean restart 
    3131   USE lib_mpp         ! distributed memory computing library 
    3232   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    199199         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    200200      END DO 
    201       hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 
    202       hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 
     201      hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) ) 
     202      hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) ) 
    203203 
    204204      ! Restoring frequencies for z_tilde coordinate 
     
    545545      END DO 
    546546      !                                        ! Inverse of the local depth 
    547       hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
    548       hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     547      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 
     548      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 
    549549 
    550550      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5506 r5619  
    365365      !!              - bathy : meter bathymetry (in meters) 
    366366      !!---------------------------------------------------------------------- 
    367       INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
     367      INTEGER  ::   ji, jj, jk            ! dummy loop indices 
    368368      INTEGER  ::   inum                      ! temporary logical unit 
    369369      INTEGER  ::   ierror                    ! error flag 
     
    472472         risfdep(:,:)=0.e0 
    473473         misfdep(:,:)=1 
     474         ! 
     475         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
     476         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN  
     477            risfdep(:,:)=200.e0  
     478            misfdep(:,:)=1  
     479            ij0 = 1 ; ij1 = 40  
     480            DO jj = mj0(ij0), mj1(ij1)  
     481               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
     482            END DO  
     483            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
     484         !  
     485         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 
     486         !  
     487            risfdep(:,:)=0.e0 
     488            misfdep(:,:)=1 
     489            ij0 = 1 ; ij1 = 40 
     490            DO jj = mj0(ij0), mj1(ij1) 
     491               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
     492            END DO 
     493            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
     494         END IF 
    474495         ! 
    475496         DEALLOCATE( idta, zdta ) 
     
    529550               WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    530551            END IF 
     552            ! set grounded point to 0 
     553            WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 
     554               misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     555               mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     556            END WHERE 
    531557            !        
    532558            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
     
    952978      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    953979      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    954       REAL(wp) ::   zmax             ! Maximum depth 
    955980      REAL(wp) ::   zdiff            ! temporary scalar 
    956       REAL(wp) ::   zrefdep          ! temporary scalar 
     981      REAL(wp) ::   zmax             ! temporary scalar 
    957982      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
    958983      !!--------------------------------------------------------------------- 
     
    9931018      END DO 
    9941019 
    995       IF ( ln_isfcav ) CALL zgr_isf 
    996  
    9971020      ! Scale factors and depth at T- and W-points 
    9981021      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     
    10021025         e3w_0  (:,:,jk) = e3w_1d  (jk) 
    10031026      END DO 
    1004       !  
    1005       DO jj = 1, jpj 
    1006          DO ji = 1, jpi 
    1007             ik = mbathy(ji,jj) 
    1008             IF( ik > 0 ) THEN               ! ocean point only 
    1009                ! max ocean level case 
    1010                IF( ik == jpkm1 ) THEN 
    1011                   zdepwp = bathy(ji,jj) 
    1012                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1013                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1014                   e3t_0(ji,jj,ik  ) = ze3tp 
    1015                   e3t_0(ji,jj,ik+1) = ze3tp 
    1016                   e3w_0(ji,jj,ik  ) = ze3wp 
    1017                   e3w_0(ji,jj,ik+1) = ze3tp 
    1018                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1019                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1020                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1021                   ! 
    1022                ELSE                         ! standard case 
    1023                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1024                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1027       
     1028      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 
     1029      IF ( ln_isfcav ) CALL zgr_isf 
     1030 
     1031      ! Scale factors and depth at T- and W-points 
     1032      IF ( .NOT. ln_isfcav ) THEN 
     1033         DO jj = 1, jpj 
     1034            DO ji = 1, jpi 
     1035               ik = mbathy(ji,jj) 
     1036               IF( ik > 0 ) THEN               ! ocean point only 
     1037                  ! max ocean level case 
     1038                  IF( ik == jpkm1 ) THEN 
     1039                     zdepwp = bathy(ji,jj) 
     1040                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1041                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1042                     e3t_0(ji,jj,ik  ) = ze3tp 
     1043                     e3t_0(ji,jj,ik+1) = ze3tp 
     1044                     e3w_0(ji,jj,ik  ) = ze3wp 
     1045                     e3w_0(ji,jj,ik+1) = ze3tp 
     1046                     gdepw_0(ji,jj,ik+1) = zdepwp 
     1047                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1048                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1049                     ! 
     1050                  ELSE                         ! standard case 
     1051                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1052                     ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1053                     ENDIF 
     1054   !gm Bug?  check the gdepw_1d 
     1055                     !       ... on ik 
     1056                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1057                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1058                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1059                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1060                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1061                     e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1062                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1063                     !       ... on ik+1 
     1064                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1065                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1066                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    10251067                  ENDIF 
    1026 !gm Bug?  check the gdepw_1d 
    1027                   !       ... on ik 
    1028                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1029                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1030                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1031                   e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1032                      &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1033                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1034                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1035                   !       ... on ik+1 
    1036                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1037                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1038                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    10391068               ENDIF 
    1040             ENDIF 
    1041          END DO 
    1042       END DO 
    1043       ! 
    1044       it = 0 
    1045       DO jj = 1, jpj 
    1046          DO ji = 1, jpi 
    1047             ik = mbathy(ji,jj) 
    1048             IF( ik > 0 ) THEN               ! ocean point only 
    1049                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1050                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1051                ! test 
    1052                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1053                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1054                   it = it + 1 
    1055                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1056                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1057                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1058                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1069            END DO 
     1070         END DO 
     1071         ! 
     1072         it = 0 
     1073         DO jj = 1, jpj 
     1074            DO ji = 1, jpi 
     1075               ik = mbathy(ji,jj) 
     1076               IF( ik > 0 ) THEN               ! ocean point only 
     1077                  e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1078                  e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1079                  ! test 
     1080                  zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1081                  IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1082                     it = it + 1 
     1083                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1084                     WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1085                     WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1086                     WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1087                  ENDIF 
    10591088               ENDIF 
    1060             ENDIF 
    1061          END DO 
    1062       END DO 
    1063       ! 
    1064       IF ( ln_isfcav ) THEN 
    1065       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1066          DO jj = 1, jpj  
    1067             DO ji = 1, jpi  
    1068                ik = misfdep(ji,jj)  
    1069                IF( ik > 1 ) THEN               ! ice shelf point only  
    1070                   IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1071                   gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1072 !gm Bug?  check the gdepw_0  
    1073                !       ... on ik  
    1074                   gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1075                      &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1076                      &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1077                   e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1078                   e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1079  
    1080                   IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1081                      e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1082                   ENDIF  
    1083                !       ... on ik / ik-1  
    1084                   e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1085                   e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1086 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1087                   gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1088                ENDIF  
    1089             END DO  
    1090          END DO  
    1091       !  
    1092          it = 0  
    1093          DO jj = 1, jpj  
    1094             DO ji = 1, jpi  
    1095                ik = misfdep(ji,jj)  
    1096                IF( ik > 1 ) THEN               ! ice shelf point only  
    1097                   e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1098                   e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1099                ! test  
    1100                   zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1101                   IF( zdiff <= 0. .AND. lwp ) THEN   
    1102                      it = it + 1  
    1103                      WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1104                      WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1105                      WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1106                      WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1107                   ENDIF  
    1108                ENDIF  
    1109             END DO  
    1110          END DO  
     1089            END DO 
     1090         END DO 
    11111091      END IF 
    1112       ! END (ISF) 
    1113  
     1092      ! 
    11141093      ! Scale factors and depth at U-, V-, UW and VW-points 
    11151094      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     
    11191098         e3vw_0(:,:,jk) = e3w_1d(jk) 
    11201099      END DO 
     1100 
    11211101      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    11221102         DO jj = 1, jpjm1 
     
    11481128      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    11491129      ! 
     1130 
    11501131      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    11511132         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     
    12551236      !!---------------------------------------------------------------------- 
    12561237      !!    
    1257       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    1258       INTEGER  ::   ik, it           ! temporary integers 
    1259       INTEGER  ::   id, jd, nprocd 
     1238      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices 
     1239      INTEGER  ::   ik, it               ! temporary integers 
    12601240      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
    1261       LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     1241      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t 
     1242      REAL(wp) ::   zmax             ! Maximum and minimum depth 
     1243      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
    12621244      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    1263       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    1264       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
     1245      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
    12651246      REAL(wp) ::   zdiff            ! temporary scalar 
    1266       REAL(wp) ::   zrefdep          ! temporary scalar 
    1267       REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
    12681247      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
    12691248      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
     
    12801259      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    12811260      END WHERE   
     1261 
     1262      ! set grounded point to 0 
     1263      WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 
     1264         misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
     1265         mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     1266      END WHERE 
    12821267 
    12831268      ! Compute misfdep for ocean points (i.e. first wet level)  
     
    12921277         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    12931278      END WHERE 
     1279 
     1280      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1281      IF ( cp_cfg .NE. "isomip" ) THEN 
     1282         WHERE (risfdep(:,:) < 100 ) 
     1283            misfdep = 1; risfdep = 0.0_wp; 
     1284         END WHERE 
     1285      END IF 
    12941286  
    12951287! 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 
     
    12971289! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    12981290      DO jl = 1, 10      
    1299          WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
     1291         WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 ) 
    13001292            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    13011293            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     
    13261318 
    13271319         ! split last cell if possible (only where water column is 2 cell or less) 
    1328          DO jk = jpkm1, 1, -1 
    1329             zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1330             WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1331                mbathy(:,:) = jk 
    1332                bathy(:,:)  = zmax 
    1333             END WHERE 
    1334          END DO 
     1320         !DO jk = jpkm1, 1, -1 
     1321         !   zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1322         !   WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
     1323         !      mbathy(:,:) = jk 
     1324         !      bathy(:,:)  = zmax 
     1325         !   END WHERE 
     1326         !END DO 
    13351327  
    13361328         ! split top cell if possible (only where water column is 2 cell or less) 
     
    13501342               ! test bathy 
    13511343               IF (risfdep(ji,jj) .GT. 1) THEN 
    1352                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    1353                  &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1354                zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    1355                  &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1344               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1345               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    13561346  
    13571347                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1358                      IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1359                         bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1360                         mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1361                      ELSE 
     1348!                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1349!                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1350!                        mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1351!                     ELSE 
    13621352                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    13631353                        misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1364                      END IF 
     1354!                     END IF 
    13651355                  END IF 
    13661356               END IF 
     
    13741364               ! test bathy 
    13751365               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1376                   zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)& 
    1377                     & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1378                   zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  & 
    1379                     & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1380                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1381                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1382                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1383                   ELSE 
     1366!                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1367!                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1368!                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1369!                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1370!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1371!                  ELSE 
    13841372                     misfdep(ji,jj)= misfdep(ji,jj) - 1 
    13851373                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1386                   END IF 
     1374!                  END IF 
    13871375               ENDIF 
    13881376            END DO 
     
    13931381            DO ji = 1, jpim1 
    13941382               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1395                   zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   & 
    1396                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1397                   zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  & 
    1398                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1399                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1400                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1401                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) & 
    1402                    &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1403                   ELSE 
     1383!                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
     1384!                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
     1385!                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1386!                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1387!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
     1388!                  ELSE 
    14041389                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1405                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 
    1406                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1407                   END IF 
     1390                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1391!                  END IF 
    14081392               ENDIF 
    14091393            END DO 
     
    14241408            DO ji = 1, jpim1 
    14251409               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1426                   zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 
    1427                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1428                   zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  & 
    1429                    &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1430                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1431                      mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1432                      bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1433                   ELSE 
     1410!                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
     1411!                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
     1412!                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1413!                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
     1414!                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     1415!                  ELSE 
    14341416                     misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    14351417                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1436                   END IF 
     1418!                  END IF 
    14371419               ENDIF 
    14381420            END DO 
     
    14551437            DO ji = 1, jpim1 
    14561438               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1457                   zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 
    1458                     &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1459                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 
    1460                     &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1461                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1462                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1463                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1464                   ELSE 
     1439!                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
     1440!                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
     1441!                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1442!                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1443!                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1444!                  ELSE 
    14651445                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    14661446                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1467                   END IF 
     1447!                  END IF 
    14681448               ENDIF 
    14691449            ENDDO 
     
    14851465            DO ji = 1, jpim1 
    14861466               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1487                   zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 
    1488                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1489                   zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  & 
    1490                       &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1491                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1492                      mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1493                      bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  & 
    1494                       &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1495                   ELSE 
     1467!                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
     1468!                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
     1469!                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1470!                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
     1471!                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
     1472!                  ELSE 
    14961473                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1497                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) & 
    1498                       &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1499                   END IF 
     1474                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1475!                  END IF 
    15001476               ENDIF 
    15011477            ENDDO 
     
    17291705! end check compatibility ice shelf/bathy 
    17301706      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1731       WHERE (misfdep(:,:) <= 5) 
    1732          misfdep = 1; risfdep = 0.0_wp; 
    1733       END WHERE 
    1734  
    17351707      IF( icompt == 0 ) THEN  
    17361708         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
     
    17381710         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
    17391711      ENDIF  
     1712 
     1713      ! compute scale factor and depth at T- and W- points 
     1714      DO jj = 1, jpj 
     1715         DO ji = 1, jpi 
     1716            ik = mbathy(ji,jj) 
     1717            IF( ik > 0 ) THEN               ! ocean point only 
     1718               ! max ocean level case 
     1719               IF( ik == jpkm1 ) THEN 
     1720                  zdepwp = bathy(ji,jj) 
     1721                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1722                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1723                  e3t_0(ji,jj,ik  ) = ze3tp 
     1724                  e3t_0(ji,jj,ik+1) = ze3tp 
     1725                  e3w_0(ji,jj,ik  ) = ze3wp 
     1726                  e3w_0(ji,jj,ik+1) = ze3tp 
     1727                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1728                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1729                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1730                  ! 
     1731               ELSE                         ! standard case 
     1732                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1733                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1734                  ENDIF 
     1735      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1736!gm Bug?  check the gdepw_1d 
     1737                  !       ... on ik 
     1738                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1739                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1740                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1741                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  ) 
     1742                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1) 
     1743                  !       ... on ik+1 
     1744                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1745                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1746               ENDIF 
     1747            ENDIF 
     1748         END DO 
     1749      END DO 
     1750      ! 
     1751      it = 0 
     1752      DO jj = 1, jpj 
     1753         DO ji = 1, jpi 
     1754            ik = mbathy(ji,jj) 
     1755            IF( ik > 0 ) THEN               ! ocean point only 
     1756               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1757               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1758               ! test 
     1759               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1760               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1761                  it = it + 1 
     1762                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1763                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1764                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1765                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1766               ENDIF 
     1767            ENDIF 
     1768         END DO 
     1769      END DO 
     1770      ! 
     1771      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1772      DO jj = 1, jpj  
     1773         DO ji = 1, jpi  
     1774            ik = misfdep(ji,jj)  
     1775            IF( ik > 1 ) THEN               ! ice shelf point only  
     1776               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1777               gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1778!gm Bug?  check the gdepw_0  
     1779            !       ... on ik  
     1780               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1781                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1782                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1783               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1784               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1785 
     1786               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1787                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1788               ENDIF  
     1789            !       ... on ik / ik-1  
     1790               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1791               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1792! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1793               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1794            ENDIF  
     1795         END DO  
     1796      END DO  
     1797    
     1798      it = 0  
     1799      DO jj = 1, jpj  
     1800         DO ji = 1, jpi  
     1801            ik = misfdep(ji,jj)  
     1802            IF( ik > 1 ) THEN               ! ice shelf point only  
     1803               e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1804               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1805            ! test  
     1806               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1807               IF( zdiff <= 0. .AND. lwp ) THEN   
     1808                  it = it + 1  
     1809                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1810                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1811                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1812                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1813               ENDIF  
     1814            ENDIF  
     1815         END DO  
     1816      END DO  
    17401817 
    17411818      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r5215 r5619  
    3434 
    3535   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     36   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsddmp   ! structure of input SST (file informations, fields read) 
    3637 
    3738   !! * Substitutions 
     
    6061      TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read 
    6162      TYPE(FLD_N)                   ::   sn_tem, sn_sal 
     63      TYPE(FLD_N)                   ::   sn_dmpt, sn_dmps 
    6264      !! 
    6365      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal 
     66      NAMELIST/namtra_dmpfile/ sn_dmpt, sn_dmps 
    6467      INTEGER  ::   ios 
    6568      !!---------------------------------------------------------------------- 
     
    7881902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp ) 
    7982      IF(lwm) WRITE ( numond, namtsd ) 
     83 
     84      REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 
     85      READ  ( numnam_ref, namtra_dmpfile, IOSTAT = ios, ERR = 903) 
     86903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 
     87 
     88      REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term 
     89      READ  ( numnam_cfg, namtra_dmpfile, IOSTAT = ios, ERR = 904 ) 
     90904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 
    8091 
    8192      IF( PRESENT( ld_tradmp ) )   ln_tsd_tradmp = .TRUE.     ! forces the initialization when tradmp is used 
     
    105116         ! 
    106117         ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     118         ALLOCATE( sf_tsddmp(jpts), STAT=ierr0 ) 
    107119         IF( ierr0 > 0 ) THEN 
    108120            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN 
     
    113125                                ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
    114126         IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
     127         ! dmp file 
     128                                 ALLOCATE( sf_tsddmp(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
     129         IF( sn_dmpt%ln_tint )   ALLOCATE( sf_tsddmp(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
     130                                 ALLOCATE( sf_tsddmp(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
     131         IF( sn_dmps%ln_tint )   ALLOCATE( sf_tsddmp(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    115132         ! 
    116133         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 
     
    120137         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal 
    121138         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' ) 
     139         slf_i(jp_tem) = sn_dmpt   ;   slf_i(jp_sal) = sn_dmps 
     140         CALL fld_fill( sf_tsddmp, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' ) 
    122141         ! 
    123142      ENDIF 
     
    128147 
    129148 
    130    SUBROUTINE dta_tsd( kt, ptsd ) 
     149   SUBROUTINE dta_tsd( kt, ptsd, ptsddmp ) 
    131150      !!---------------------------------------------------------------------- 
    132151      !!                   ***  ROUTINE dta_tsd  *** 
     
    145164      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step 
    146165      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data 
     166      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), OPTIONAL, INTENT(  out) ::   ptsddmp   ! T & S data 
    147167      ! 
    148168      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
     
    155175      ! 
    156176      CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==! 
     177      IF ( PRESENT(ptsddmp) ) THEN 
     178         CALL fld_read( kt, 1, sf_tsddmp )      !==   read T & S data at kt time step   ==! 
     179         ptsddmp(:,:,:,jp_tem) = sf_tsddmp(jp_tem)%fnow(:,:,:)    ! NO mask 
     180         ptsddmp(:,:,:,jp_sal) = sf_tsddmp(jp_sal)%fnow(:,:,:)  
     181      END IF 
    157182      ! 
    158183      ! 
     
    304329         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta ) 
    305330                                        DEALLOCATE( sf_tsd              )     ! the structure itself 
     331         IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run' 
     332                                        DEALLOCATE( sf_tsddmp(jp_tem)%fnow )     ! T arrays in the structure 
     333         IF( sf_tsddmp(jp_tem)%ln_tint )   DEALLOCATE( sf_tsddmp(jp_tem)%fdta ) 
     334                                        DEALLOCATE( sf_tsddmp(jp_sal)%fnow )     ! S arrays in the structure 
     335         IF( sf_tsddmp(jp_sal)%ln_tint )   DEALLOCATE( sf_tsddmp(jp_sal)%fdta ) 
     336                                        DEALLOCATE( sf_tsddmp              )     ! the structure itself 
    306337      ENDIF 
    307338      ! 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5332 r5619  
    4747   USE wrk_nemo        ! Memory allocation 
    4848   USE timing          ! Timing 
     49   USE sbc_iscpl 
    4950 
    5051   IMPLICIT NONE 
     
    9192         !                                    ! ------------------- 
    9293         CALL rst_read                           ! Read the restart file 
     94         IF (ln_iscpl) CALL rst_iscpl            ! extraloate restart to wet and dry 
    9395         CALL day_init                           ! model calendar (using both namelist and restart infos) 
    9496      ELSE 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r5147 r5619  
    7373   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
    7474#else 
    75    REAL(wp), PUBLIC ::   rhoic    =  900._wp         !: volumic mass of sea ice                               [kg/m3] 
     75   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    7676   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: conductivity of the ice                               [W/m/K] 
    7777   REAL(wp), PUBLIC ::   rcpic    =    1.8837e+6_wp  !: volumetric specific heat for ice                      [J/m3/K] 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r5516 r5619  
    2929   USE sbcrnf          ! river runoff  
    3030   USE sbcisf          ! ice shelf  
     31   USE sbc_iscpl        ! ice shelf  
    3132   USE cla             ! cross land advection             (cla_div routine) 
    3233   USE in_out_manager  ! I/O manager 
     
    328329 
    329330      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field) 
    330       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
     331      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div  ( hdivn )      ! ice shelf (update hdivn field) 
     332      IF( ln_iscpl  .AND. ln_hfb )            CALL sbc_iscpl_div( hdivn )      ! ice shelf (update hdivn field) 
    331333      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    332334      ! 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r5467 r5619  
    353353            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    354354         END DO 
    355          hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 
    356          hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 
     355         hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     356         hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    357357      ENDIF 
    358358      ! 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r5407 r5619  
    2424   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
    2525   USE divcur          ! hor. divergence and curl      (div & cur routines) 
     26   USE wrk_nemo 
    2627 
    2728   IMPLICIT NONE 
     
    145146                     CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       ) 
    146147#endif 
     148 
     149 
     150 
     151                  IF ( ln_iscpl ) THEN  
     152                     CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask     ) ! need to extrapolate T/S 
     153                     CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask     ) ! need to correct barotropic velocity 
     154                     CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask     ) ! need to correct barotropic velocity 
     155                     CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask    ) ! need to correct barotropic velocity 
     156                     CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) ! need to compute temperature correction 
     157                     CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) ) ! need to compute volume      correction  ???? 
     158                     CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) ) ! need to compute volume      correction  ???? 
     159                     CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute volume      correction  ???? 
     160                  END IF 
    147161      IF( kt == nitrst ) THEN 
    148162         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90

    r5429 r5619  
    3131   END INTERFACE 
    3232 
     33   INTERFACE lbc_sum 
     34      MODULE PROCEDURE mpp_sum_3d, mpp_sum_2d 
     35   END INTERFACE 
     36 
    3337   INTERFACE lbc_bdy_lnk 
    3438      MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d 
     
    4549   PUBLIC lbc_lnk       ! ocean lateral boundary conditions 
    4650   PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions 
     51   PUBLIC lbc_sum 
    4752   PUBLIC lbc_lnk_e 
    4853   PUBLIC lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

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

    r5130 r5619  
    136136 
    137137      imask(:,:)=1 
    138       WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0 
     138      WHERE ( zdta(:,:) - zdtaisf(:,:) <= 1e-2 ) imask = 0 
    139139 
    140140      !  1. Dimension arrays for subdomains 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r5120 r5619  
    111111            zcoef = z_fwf * rcp 
    112112            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1) 
     113            sfx(:,:) = sfx(:,:) + z_fwf * sss_m      * tmask(:,:,1) 
    113114            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    114115         ENDIF 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r5541 r5619  
    6363 
    6464   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
    65    REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ? 
     65   REAL(wp), PUBLIC, SAVE ::   kappa  =    0.0_wp    ! phycst ? 
    6666   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ? 
    6767   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ? 
     
    125125         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf  
    126126         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
     127         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
     128         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
    127129         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
    128130         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 
     
    150152            !: read effective lenght (BG03) 
    151153            IF (nn_isf == 2) THEN 
    152                ! Read Data and save some integral values 
     154               cvarLeff = 'soLeff'  
    153155               CALL iom_open( sn_Leff_isf%clname, inum ) 
    154                cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale 
    155156               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
    156157               CALL iom_close(inum) 
     
    237238            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    238239            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 
    239             CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U') 
    240             CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V') 
     240            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
     241            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
    241242            ! iom print 
    242243            CALL iom_put('ttbl',ttbl(:,:)) 
    243244            CALL iom_put('stbl',stbl(:,:)) 
    244             CALL iom_put('utbl',utbl(:,:)) 
    245             CALL iom_put('vtbl',vtbl(:,:)) 
     245            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
     246            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
    246247            ! compute fwf and heat flux 
    247248            CALL sbc_isf_cav (kt) 
     
    296297         !  
    297298         ! output 
    298          CALL iom_put('qisf'  , qisf) 
    299          IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce ) 
     299         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf) 
     300         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf) 
    300301      END IF 
    301302   
     
    416417      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b 
    417418      REAL(wp) ::   zgammat, zgammas 
    418       REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==! 
     419      REAL(wp) ::   zeps = 1.e-20_wp        !==   Local constant initialization   ==! 
    419420      INTEGER  ::   ji, jj     ! dummy loop indices 
    420421      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
     
    506507                     zeps1=rcp*rau0*zgammat 
    507508                     zeps2=lfusisf*rau0*zgammas 
    508                      zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj) 
     509                     zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 
    509510                     zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    510511                     zeps6=zeps4-zti(ji,jj) 
    511512                     zeps7=zeps4-tsurf 
    512513                     zaqe=zlamb1 * (zeps1 + zeps3) 
    513                      zaqer=0.5/zaqe 
     514                     zaqer=0.5/MIN(zaqe,-zeps) 
    514515                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 
    515516                     zcqe=zeps2*stbl(ji,jj) 
     
    520521                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    521522   
    522 ! zfwflx is upward water flux 
    523                      zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz ) 
     523! zfwflx is upward water flux (MAX usefull if kappa = 0 
     524                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 
    524525! zhtflx is upward heat flux (out of ocean) 
    525526! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 
     
    527528! zwflx is upward water flux 
    528529! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 
    529                      zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 
     530!!!!!!!!                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 
    530531! test convergence and compute gammat 
    531532                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r5541 r5619  
    921921               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    922922                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    923                   &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     923                  &            / fse3w(ji,jj,jk) * wmask(ji,jj,jk) 
    924924            END DO 
    925925         END DO 
     
    12421242      ! 
    12431243      rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
    1244       rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
     1244      rcp         = 3974._wp      !: heat capacity     [J/K] 
    12451245      ! 
    12461246      IF(lwp) THEN                ! Control print 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5102 r5619  
    102102      REAL(wp) ::   zta, zsa             ! local scalars 
    103103      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta  
     104      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dtadmp  
    104105      !!---------------------------------------------------------------------- 
    105106      ! 
    106107      IF( nn_timing == 1 )  CALL timing_start( 'tra_dmp') 
    107108      ! 
    108       CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta ) 
    109       ! 
     109      CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta, zts_dtadmp ) 
    110110      !                           !==   input T-S data at kt   ==! 
    111       CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt 
     111      CALL dta_tsd( kt, zts_dta, zts_dtadmp )            ! read and interpolates T-S data at kt 
     112      zts_dta=zts_dtadmp 
    112113      ! 
    113114      SELECT CASE ( nn_zdmp )     !==    type of damping   ==! 
     
    175176         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    176177      ! 
    177       CALL wrk_dealloc( jpi, jpj, jpk, jpts,  zts_dta ) 
     178      CALL wrk_dealloc( jpi, jpj, jpk, jpts,  zts_dta, zts_dtadmp ) 
    178179      ! 
    179180      IF( nn_timing == 1 )  CALL timing_stop( 'tra_dmp') 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r5431 r5619  
    2222   USE sbcrnf          ! River runoff   
    2323   USE sbcisf          ! Ice shelf    
     24   USE sbc_iscpl       ! 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 
    121       REAL(wp) ::   zalpha, zhk 
    122       REAL(wp) ::  zt_frz, zpress 
     121      REAL(wp) ::   zt_frz, zpress 
    123122      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    124123      !!---------------------------------------------------------------------- 
     
    220219      ! 
    221220      IF( nn_isf > 0 ) THEN 
    222          zfact = 0.5e0 
     221         zfact = 0.5_wp 
    223222         DO jj = 2, jpj 
    224223            DO ji = fs_2, fs_jpim1 
     
    233232               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    234233!                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
    235                   zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
     234                  zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
    236235               ! compute trend 
    237236                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
     
    246245               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    247246!               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04 
    248                zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 
     247               zt_frz = -1.9_wp !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 
    249248               ! compute trend 
    250249               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
     
    288287      ENDIF 
    289288  
    290       IF( l_trdtra )   THEN                      ! send trends for further diagnostics 
     289      !---------------------------------------- 
     290      !        Ice Sheet coupling imbalance correction to have conservation 
     291      !---------------------------------------- 
     292      ! 
     293      IF( ln_iscpl .AND. ln_hfb) THEN         ! input of heat and salt due to river runoff  
     294         DO jk = 1,jpk 
     295            DO jj = 2, jpj  
     296               DO ji = fs_2, fs_jpim1 
     297                  zdep = 1._wp / fse3t_n(ji,jj,jk)  
     298                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem)                       & 
     299                      &                                         * zdep 
     300                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal)                       & 
     301                      &                                         * zdep   
     302               END DO   
     303            END DO   
     304         END DO 
     305      ENDIF 
     306 
     307      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    291308         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    292309         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90

    r4161 r5619  
    2525 
    2626   PUBLIC   glob_sum   ! used in many places 
     27   PUBLIC   glob_sum_full   ! used in many places 
    2728   PUBLIC   DDPDD      ! also used in closea module 
    2829   PUBLIC   glob_min, glob_max 
     
    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_2d *** 
     166      !! 
     167      !! ** Purpose : perform a sum in calling DDPDD routine 
     168      !!---------------------------------------------------------------------- 
     169      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     170      REAL(wp)                             ::   glob_sum_full_2d   ! global masked sum 
     171      !! 
     172      !!----------------------------------------------------------------------- 
     173      ! 
     174      glob_sum_full_2d = SUM( ptab(:,:) ) 
     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_3d *** 
     182      !! 
     183      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 
     184      !!---------------------------------------------------------------------- 
     185      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     186      REAL(wp)                               ::   glob_sum_full_3d   ! global masked 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_3d = 0.e0 
     195      DO jk = 1, ijpk 
     196         glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk) ) 
     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_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 masked sum 
     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) 
     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_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 masked sum 
     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) 
     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.