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

Changeset 11083


Ignore:
Timestamp:
2019-06-06T16:29:54+02:00 (6 years ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0_GO6_mixing : Update to be relative to rev 11081 of NEMO_4.0_mirror.

Location:
NEMO/branches/UKMO/NEMO_4.0_GO6_mixing
Files:
1 deleted
40 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/INSTALL.rst

    r10650 r11083  
    104104   .. code:: console 
    105105 
    106       $ svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5 
     106      $ svn co https://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5 
    107107 
    108108Download the NEMO source code 
    … …  
    111111.. code:: console 
    112112 
    113    $ svn co http://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-4.0 
     113   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-4.0 
    114114 
    115115Description of directory tree 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/README.rst

    r10606 r11083  
    1 :Release: |version| 
    2 :Date:    |today| 
     1:Release:  |release| 
     2:Date:     |today| 
     3:SVN rev.: |revision| 
    34 
    45NEMO_ for **Nucleus for European Modelling of the Ocean** is a state-of-the-art modelling framework for 
    … …  
    9091.. Substitutions / Links 
    9192 
    92 .. |NEMO manual| image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1464816.svg 
    93 .. |NEMO guide|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1475325.svg 
    94 .. |SI3 manual|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1471689.svg 
    95 .. |TOP manual|  image:: http://zenodo.org/badge/DOI/10.5281/zenodo.1471700.svg 
     93.. |NEMO manual| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1464816.svg 
     94.. |NEMO guide|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1475325.svg 
     95.. |SI3 manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471689.svg 
     96.. |TOP manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471700.svg 
    9697 
    9798.. |NEMO strategy| replace:: multi-year development strategy 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml

    r9943 r11083  
    7979       <field field_ref="vfxsnw"           name="vfxsnw" /> 
    8080        
    81        <!-- diag error for negative ice volume after advection --> 
    82        <field field_ref="iceneg_pres"      name="sineg_pres" /> 
    83        <field field_ref="iceneg_volu"      name="sineg_volu" /> 
    84        <field field_ref="iceneg_hfx"       name="sineg_hfx"  /> 
    85  
    8681       <!-- categories --> 
    8782       <field field_ref="icemask_cat"      name="simskcat"/> 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/README.rst

    r10606 r11083  
    235235.. literalinclude:: ../../../cfgs/GYRE_PISCES/EXPREF/namelist_ref 
    236236   :language: fortran 
    237    :lines: 306-333 
     237   :lines: 935-960 
    238238 
    239239Input dynamical fields for this configuration (``ORCA2_OFF_v4.0.tar``) comes from 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/SHARED/field_def_nemo-ice.xml

    r10578 r11083  
    160160          <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" /> 
    161161 
    162      <!-- diagnostics of the negative values resulting from the advection scheme --> 
    163      <field id="iceneg_pres"  long_name="Fraction of time steps with negative sea ice volume"                    unit=""     /> 
    164      <field id="iceneg_volu"  long_name="Negative sea ice volume per area arising from advection"                unit="m"    /> 
    165           <field id="iceneg_hfx"   long_name="Negative sea ice heat content (eq. heat flux) arising from advection"   unit="W/m2" /> 
    166  
    167162     <!-- sbcssm variables --> 
    168163          <field id="sst_m"    unit="degC" /> 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/SHARED/namelist_ice_ref

    r10612 r11083  
    154154&namthd_do      !   Ice growth in open water 
    155155!------------------------------------------------------------------------------ 
    156    rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice 
     156   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_himin 
    157157   ln_frazil        = .false.         !  Frazil ice parameterization (ice collection as a function of wind) 
    158158      rn_maxfraz    =   1.0           !     maximum fraction of frazil ice collecting at the ice base 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml

    r9943 r11083  
    7979       <field field_ref="vfxsnw"           name="vfxsnw" /> 
    8080 
    81        <!-- diag error for negative ice volume after advection --> 
    82        <field field_ref="iceneg_pres"      name="sineg_pres" /> 
    83        <field field_ref="iceneg_volu"      name="sineg_volu" /> 
    84        <field field_ref="iceneg_hfx"       name="sineg_hfx"  /> 
    85         
    8681       <!-- categories --> 
    8782       <field field_ref="icemask_cat"      name="simskcat"/> 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/SPITZ12/EXPREF/namelist_cfg

    r10075 r11083  
    8080                     ! Sea-ice : 
    8181   nn_ice      = 2         !  SI3 
     82   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges) 
     83      !                    !  =F levitating ice (no pressure, mass and salt exchanges) 
    8284                     ! Misc. options of sbc : 
    8385   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr ) 
    … …  
    9092      ! 
    9193      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    92       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
     94      ln_Cd_L15   = .true.    !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    9395      ! 
    9496   cn_dir = './'  !  root directory for the bulk data location 
    … …  
    9698   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask ! 
    9799   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   ! 
    98 !! ERAI 
    99 !   sn_wndi     = 'u10_era_spitz'               ,         6         , 'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , '' 
    100 !   sn_wndj     = 'v10_era_spitz'               ,         6         , 'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , '' 
    101 !   sn_qsr      = 'ssrd_era_spitz'              ,         6         , 'ssrd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    102 !   sn_qlw      = 'strd_era_spitz'              ,         6         , 'strd'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    103 !   sn_tair     = 't2_era_spitz'                ,         6         , 't2'      ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    104 !   sn_humi     = 'humi_era_spitz'              ,         6         , 'humi'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    105 !   sn_prec     = 'precip_era_spitz'            ,         6         , 'precip'  ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    106 !   sn_snow     = 'snow_era_spitz'              ,         6         , 'snow'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    107 !   sn_tdif     = 'taudif'                      ,         6         , 'taudif'  ,   .true.    , .false. , 'yearly'  , ''             , '' , '' 
    108 !      rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    109 !! MAR 
    110100   sn_wndi     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'u10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Uwnd' , '' 
    111101   sn_wndj     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,   2 ,  'v10'     ,   .true.    , .false. , 'yearly'  , 'weights_bicub', 'Vwnd' , '' 
    … …  
    129119   !                       ! type of penetration                        (default: NO selection) 
    130120   ln_qsr_rgb  = .true.       !  RGB light penetration (Red-Green-Blue) 
    131 / 
    132 !----------------------------------------------------------------------- 
    133 &namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T) 
    134 !----------------------------------------------------------------------- 
    135 / 
    136 !----------------------------------------------------------------------- 
    137 &namsbc_rnf    !   runoffs                                              (ln_rnf =T) 
    138 !----------------------------------------------------------------------- 
    139121/ 
    140122!!====================================================================== 
    … …  
    236218!----------------------------------------------------------------------- 
    237219   ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
     220   ln_drgimp   = .true.   !  implicit top/bottom friction flag 
    238221/ 
    239222!----------------------------------------------------------------------- 
    240223&namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
    241224!----------------------------------------------------------------------- 
    242    rn_Cd0      =  2.5e-3   !!1.e-3    !  drag coefficient [-] 
    243    rn_Cdmax    =  0.01     !!0.1      !  drag value maximum [-] (logarithmic drag) 
    244    rn_ke0      =  0.       !!2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases) 
     225   rn_Cd0      =  2.5e-3   !  drag coefficient [-] 
     226   rn_Cdmax    =  0.1      !  drag value maximum [-] (logarithmic drag) 
     227   rn_ke0      =  0.       !  background kinetic energy  [m2/s2] (non-linear cases) 
    245228   rn_z0       =  3.e-3    !  roughness [m] (ln_loglayer=T) 
    246    ln_boost    = .false.   !  =T regional boost of Cd0 ; =F constant 
    247       rn_boost =  50.         !  local boost factor  [-] 
    248 / 
    249 !----------------------------------------------------------------------- 
    250 &nambbc        !   bottom temperature boundary condition                (default: OFF) 
    251 !----------------------------------------------------------------------- 
    252    ln_trabbc   = .false.    !  Apply a geothermal heating at the ocean bottom 
    253229/ 
    254230!----------------------------------------------------------------------- 
    … …  
    258234      nn_bbl_ldf  =  1        !  diffusive bbl (=1)   or not (=0) 
    259235      nn_bbl_adv  =  0        !  advective bbl (=1/2) or not (=0) 
    260       rn_ahtbbl   =  1000.    !  lateral mixing coefficient in the bbl  [m2/s] 
    261       rn_gambbl   =  10.      !  advective bbl coefficient                 [s] 
    262236/ 
    263237!!====================================================================== 
    … …  
    293267   !                       !  Coefficients: 
    294268   nn_aht_ijk_t    = 31        !  space/time variation of eddy coefficient: 
    295       !                             !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
     269   !                                !   = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 
    296270/ 
    297271!!====================================================================== 
    … …  
    320294&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    321295!----------------------------------------------------------------------- 
    322    ln_dynvor_enT = .true. !  energy conserving scheme (T-point) 
     296   ln_dynvor_eeT = .true.  !  energy conserving scheme (een using e3t) 
    323297/ 
    324298!----------------------------------------------------------------------- 
    … …  
    352326&namzdf        !   vertical physics manager                             (default: NO selection) 
    353327!----------------------------------------------------------------------- 
     328   !                       ! adaptive-implicit vertical advection 
     329   ln_zad_Aimp = .true.      !  Courant number dependent scheme (Shchepetkin 2015) 
    354330   !                       ! type of vertical closure (required) 
    355331   ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke) 
    356332   ln_zdfgls   = .false.      !  Generic Length Scale closure           (T =>   fill namzdf_gls) 
     333   !                       ! convection 
     334   ln_zdfevd   = .true.       !  enhanced vertical diffusion 
     335   ! 
    357336   ln_zdfddm   = .true.    ! double diffusive mixing 
     337   ! 
    358338   !                       !  Coefficients 
    359339   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F) 
    … …  
    384364/ 
    385365!----------------------------------------------------------------------- 
    386 &nam_diaharm   !   Harmonic analysis of tidal constituents              ('key_diaharm') 
    387 !----------------------------------------------------------------------- 
    388 / 
    389 !----------------------------------------------------------------------- 
    390366&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4") 
    391367!----------------------------------------------------------------------- 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r10535 r11083  
    3535&namdyn         !   Ice dynamics 
    3636!------------------------------------------------------------------------------ 
     37   ln_landfast_L16  = .true.          !  landfast: parameterization from Lemieux 2016 
    3738/ 
    3839!------------------------------------------------------------------------------ 
    … …  
    4344&namdyn_rhg     !   Ice rheology 
    4445!------------------------------------------------------------------------------ 
     46   ln_rhg_EVP       = .true.          !  EVP rheology 
     47      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    4548/ 
    4649!------------------------------------------------------------------------------ 
    … …  
    6770&namthd_do      !   Ice growth in open water 
    6871!------------------------------------------------------------------------------ 
    69    rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_hnewice 
     72   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_himin 
    7073   ln_frazil        = .true.          !  Frazil ice parameterization (ice collection as a function of wind) 
    7174/ 
    … …  
    7780&namthd_pnd     !   Melt ponds 
    7881!------------------------------------------------------------------------------ 
    79    ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
    80    ln_pnd_CST       = .false.         !  activate constant melt ponds 
    81    ln_pnd_alb       = .false.         !  melt ponds affect albedo or not 
     82   ln_pnd_H12       = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
     83   ln_pnd_alb       = .true.          !  melt ponds affect albedo or not 
    8284/ 
    8385 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/Makefile

    r10606 r11083  
    2323# Browse to 127.0.0.1:8000/NEMO_guide.html 
    2424htmllive: 
    25    sphinx-autobuild $(SPHINXOPTS) $(SOURCEDIR) $(BUILDDIR)/html 
     25   sphinx-autobuild $(SPHINXOPTS) $(SOURCEDIR) $(BUILDDIR)/htmllive 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/_templates/layout.html

    r10279 r11083  
    1212   <p>Community ocean model for multifarious space and time scales</p> 
    1313 
    14    <div class="version">{{ release }}</div> 
     14   <div class="version">{{ version }}</div> 
    1515 
    1616   {% include "searchbox.html" %} 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/doc/rst/source/conf.py

    r10600 r11083  
    9292# -- Customisation ----------------------------------------------------------- 
    9393 
     94# Timestamping 
    9495import datetime 
    9596year = datetime.date.today().year 
    … …  
    108109# Include common directives for every rst file 
    109110rst_epilog = open('global.rst', 'r').read() 
     111 
     112# SVN revision 
     113import subprocess 
     114revision = subprocess.check_output("svnversion").decode("utf-8") 
     115rst_prolog = '.. |revision| replace:: %s' % revision 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icecor.F90

    r10888 r11083  
    6666         WRITE(numout,*) '~~~~~~~' 
    6767      ENDIF 
    68       ! 
    6968      !                             !----------------------------------------------------- 
    70       !                             !  ice thickness must exceed himin (for ice diff)    ! 
     69      !                             !  ice thickness must exceed himin (for temp. diff.) ! 
    7170      !                             !----------------------------------------------------- 
    7271      WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 
    … …  
    9796         END DO 
    9897      ENDIF 
    99  
    10098      !                             !----------------------------------------------------- 
    10199      !                             !  Rebin categories with thickness out of bounds     ! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icectl.F90

    r10888 r11083  
    6969      !! 
    7070      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax  
     71      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
    7272      REAL(wp) ::   zvtrp, zetrp 
    7373      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
    … …  
    141141         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    142142 
    143          zvmin = glob_min( 'icectl', v_i ) 
    144          zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    145          zamin = glob_min( 'icectl', a_i ) 
     143         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     144         zvmin  = glob_min( 'icectl', v_i ) 
     145         zamin  = glob_min( 'icectl', a_i ) 
     146         zsmin  = glob_min( 'icectl', sv_i ) 
     147         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     148         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    146149 
    147150         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    152155 
    153156         IF(lwp) THEN 
    154             IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
    155             IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
    156             IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
    157             IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    158             IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10   & 
    159                & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) & 
    160                &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    161             IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    162 !clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90) 
     157            ! check conservation issues 
     158            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
     159            IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
     160            IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
     161            ! check maximum ice concentration 
     162            IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  & 
     163               &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     164            ! check negative values 
     165            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     166            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     167            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin 
     168            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin 
     169            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin 
     170!clem: the following check fails (I think...) 
    163171!            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 
    164172!                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn.F90

    r10888 r11083  
    7474      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7575      !! 
    76       INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     76      INTEGER  ::   ji, jj        ! dummy loop indices 
    7777      REAL(wp) ::   zcoefu, zcoefv 
    78       REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
     78      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i 
    8079      !!-------------------------------------------------------------------- 
    8180      ! 
    … …  
    8988      ENDIF 
    9089      !                       
    91       IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    92          tau_icebfr(:,:) = 0._wp 
    93          DO jl = 1, jpl 
    94             WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    95          END DO 
    96       ENDIF 
    97       ! 
    98       !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
    99       DO jl = 1, jpl 
    100          DO jj = 2, jpjm1 
    101             DO ji = fs_2, fs_jpim1 
    102                zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), & 
    103                   &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
    104                   &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
    105                   &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
    106                zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), & 
    107                   &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
    108                   &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
    109                   &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
    110                zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), & 
    111                   &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
    112                   &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
    113                   &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
    114             END DO 
    115          END DO 
    116       END DO 
    117       CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
    118       ! 
    119       ! 
    120       SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     90      ! retrieve thickness from volume for landfast param. and UMx advection scheme 
     91      WHERE( a_i(:,:,:) >= epsi20 ) 
     92         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 
     93         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 
     94      ELSEWHERE 
     95         h_i(:,:,:) = 0._wp 
     96         h_s(:,:,:) = 0._wp 
     97      END WHERE 
     98      ! 
     99      WHERE( a_ip(:,:,:) >= epsi20 ) 
     100         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     101      ELSEWHERE 
     102         h_ip(:,:,:) = 0._wp 
     103      END WHERE 
     104      ! 
     105      ! 
     106      SELECT CASE( nice_dyn )          !-- Set which dynamics is running 
    121107 
    122108      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
    123          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    124          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    125          CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
    126          CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
    127  
     109         ! 
     110         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     111         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     112         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting  
     113         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections 
     114         ! 
    128115      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    129          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    130          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    131          CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
    132  
     116         ! 
     117         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     118         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     119         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting) 
     120         CALL ice_var_zapsmall                                              ! -- zap small areas 
     121         ! 
    133122      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
    134          ALLOCATE( zdivu_i(jpi,jpj) ) 
     123         ! 
    135124         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
    136125         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
    … …  
    145134         END DO 
    146135         ! --- 
    147          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    148  
    149          ! diagnostics: divergence at T points  
    150          DO jj = 2, jpjm1 
    151             DO ji = 2, jpim1 
    152                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    153                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
    154             END DO 
    155          END DO 
    156          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    157          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    158  
    159          DEALLOCATE( zdivu_i ) 
    160  
     136         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     137         ! 
    161138      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
    162          ALLOCATE( zdivu_i(jpi,jpj) ) 
     139         ! 
    163140         u_ice(:,:) = rn_uice * umask(:,:,1) 
    164141         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    … …  
    166143         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
    167144         ! --- 
    168          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    169  
    170          ! diagnostics: divergence at T points  
    171          DO jj = 2, jpjm1 
    172             DO ji = 2, jpim1 
    173                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    174                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     145         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     146 
     147      END SELECT 
     148      ! 
     149      ! 
     150      ! diagnostics: divergence at T points  
     151      IF( iom_use('icediv') ) THEN 
     152         ! 
     153         SELECT CASE( nice_dyn ) 
     154 
     155         CASE ( np_dynADV1D , np_dynADV2D ) 
     156 
     157            ALLOCATE( zdivu_i(jpi,jpj) ) 
     158            DO jj = 2, jpjm1 
     159               DO ji = 2, jpim1 
     160                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     161                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     162               END DO 
    175163            END DO 
    176          END DO 
    177          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    178          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    179  
    180          DEALLOCATE( zdivu_i ) 
    181  
    182       END SELECT 
     164            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     165            CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     166            DEALLOCATE( zdivu_i ) 
     167 
     168         END SELECT 
     169         ! 
     170      ENDIF 
    183171      ! 
    184172      ! controls 
    … …  
    188176 
    189177 
    190    SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    191       !!------------------------------------------------------------------- 
    192       !!                  ***  ROUTINE Hbig  *** 
    193       !! 
    194       !! ** Purpose : Thickness correction in case advection scheme creates 
    195       !!              abnormally tick ice or snow 
    196       !! 
    197       !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
    198       !!                 (before advection) and reduce it by adapting ice concentration 
    199       !!              2- check whether snow thickness is larger than the surrounding 9-points 
    200       !!                 (before advection) and reduce it by sending the excess in the ocean 
    201       !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
    202       !!                 and reduce it by sending the excess in the ocean 
    203       !!              4- correct pond fraction to avoid a_ip > a_i 
    204       !! 
    205       !! ** input   : Max thickness of the surrounding 9-points 
    206       !!------------------------------------------------------------------- 
    207       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    208       ! 
    209       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    210       REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    211       !!------------------------------------------------------------------- 
    212       ! controls 
    213       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    214       ! 
    215       CALL ice_var_zapsmall                       !-- zap small areas 
    216       ! 
    217       DO jl = 1, jpl 
    218          DO jj = 1, jpj 
    219             DO ji = 1, jpi 
    220                IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    221                   ! 
    222                   !                               ! -- check h_ip -- ! 
    223                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    224                   IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 
    225                      zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 
    226                      IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 
    227                         a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
    228                      ENDIF 
    229                   ENDIF 
    230                   ! 
    231                   !                               ! -- check h_i -- ! 
    232                   ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
    233                   zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    234                   IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
    235                      a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
    236                   ENDIF 
    237                   ! 
    238                   !                               ! -- check h_s -- ! 
    239                   ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
    240                   zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
    241                   IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
    242                      zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    243                      ! 
    244                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
    245                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    246                      ! 
    247                      e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
    248                      v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    249                   ENDIF            
    250                   ! 
    251                   !                               ! -- check snow load -- ! 
    252                   ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
    253                   !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
    254                   !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
    255                   zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    256                   IF( zvs_excess > 0._wp ) THEN 
    257                      zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 ) 
    258                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
    259                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    260                      ! 
    261                      e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
    262                      v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess 
    263                   ENDIF 
    264                    
    265                ENDIF 
    266             END DO 
    267          END DO 
    268       END DO  
    269       !                                           !-- correct pond fraction to avoid a_ip > a_i 
    270       WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    271       ! 
    272       ! controls 
    273       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    274       ! 
    275    END SUBROUTINE Hbig 
    276  
    277  
    278178   SUBROUTINE Hpiling 
    279179      !!------------------------------------------------------------------- 
    … …  
    290190      ! controls 
    291191      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    292       ! 
    293       CALL ice_var_zapsmall                       !-- zap small areas 
    294192      ! 
    295193      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    … …  
    337235         WRITE(numout,*) '~~~~~~~~~~~~' 
    338236         WRITE(numout,*) '   Namelist namdyn:' 
    339          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
    340          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
    341          WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
    342          WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
    343          WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
    344          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
    345          WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
    346          WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
    347          WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
    348          WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
    349          WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
    350          WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
     237         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL 
     238         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV 
     239         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D 
     240         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D 
     241         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')' 
     242         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
     243         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
     244         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home 
     245         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
     246         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
     247         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax 
     248         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile 
    351249         WRITE(numout,*) 
    352250      ENDIF 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv.F90

    r10888 r11083  
    6464      !!---------------------------------------------------------------------- 
    6565      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    66       ! 
    67       INTEGER ::   jl   ! dummy loop indice 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0 
    6966      !!--------------------------------------------------------------------- 
    7067      ! 
    71       IF( ln_timing )   CALL timing_start('icedyn_adv') 
     68      ! controls 
     69      IF( ln_timing    )   CALL timing_start('icedyn_adv')                                                             ! timing 
     70      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    7271      ! 
    7372      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    7675         WRITE(numout,*) '~~~~~~~~~~~' 
    7776      ENDIF 
    78        
    79       ! conservation test 
    80       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    81                       
    82       !---------- 
    83       ! Advection 
    84       !---------- 
     77      ! 
     78      !---------------! 
     79      !== Advection ==! 
     80      !---------------! 
    8581      SELECT CASE( nice_adv ) 
    8682      !                                !-----------------------! 
    8783      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
    8884         !                             !-----------------------! 
    89          CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    90       !                                !-----------------------! 
     85         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     87         !                             !-----------------------! 
    9188      CASE( np_advPRA )                ! PRATHER scheme        ! 
    9289         !                             !-----------------------! 
    93          CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     90         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9492      END SELECT 
    95  
    96       !---------------------------- 
    97       ! Debug the advection schemes 
    98       !---------------------------- 
    99       ! clem: At least one advection scheme above is not strictly positive => UMx 
    100       !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes) 
    101       !       In UMx    , advected fields are not perfectly bounded and negative values can appear. 
    102       !                   These values are usually very small but in some occasions they can also be non-negligible 
    103       !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 
    104       ! 
    105       ! record the negative values resulting from UMx 
    106       zmask(:,:) = 0._wp ! keep the init to 0 here 
    107       DO jl = 1, jpl 
    108          WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp 
    109       END DO 
    110       IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0 
    111       IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only) 
    112       IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only) 
    113          &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2] 
    114       ! 
    115       ! ==> conservation is ensured by calling this routine below, 
    116       !     however the global ice volume is then changed by advection (but errors are small)  
    117       CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    11893 
    11994      !------------ 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_adv_umx.F90

    r10888 r11083  
    1111   !!   'key_si3'                                       SI3 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme 
     13   !!   ice_dyn_adv_umx   : update the tracer fields 
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    15    !!   macho             : ??? 
    16    !!   nonosc_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     15   !!   macho             : compute the fluxes 
     16   !!   nonosc_ice        : limit the fluxes using a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
    … …  
    3232 
    3333   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90 
    34        
    35    REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
    36    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    37     
    38    ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 
    39    INTEGER ::   kn_limiter = 1 
    40  
    41    ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
    42    !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why) 
    43    !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration  
    44    LOGICAL ::   ll_neg = .TRUE. 
    45     
    46    ! alternate directions for upstream 
    47    LOGICAL ::   ll_upsxy = .TRUE. 
    48  
    49    ! alternate directions for high order 
    50    LOGICAL ::   ll_hoxy = .TRUE. 
    51     
    52    ! prelimiter: use it to avoid overshoot in H 
    53    !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why) 
    54    LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
    55  
    56  
     34   ! 
     35   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1 
     36   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2 
     37   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3 
     38   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc 
     39   !                                                      2 = superbee 
     40   !                                                      3 = h3 
     41   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream 
     42   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order 
     43   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6 
     44   !                                                 then interpolate T at u/v points using the upstream scheme 
     45   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D 
     46   ! 
     47   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6 
     48   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120 
     49   ! 
     50   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small 
     51   ! 
    5752   !! * Substitutions 
    5853#  include "vectopt_loop_substitute.h90" 
    … …  
    6459CONTAINS 
    6560 
    66    SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & 
     61   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 
    6762      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    6863      !!---------------------------------------------------------------------- 
    … …  
    7974      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    8075      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     76      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     77      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     78      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    8179      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    8280      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
    … …  
    9391      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    9492      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    95       REAL(wp) ::   zdt 
    96       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
    97       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
     93      REAL(wp) ::   zdt, zvi_cen 
     94      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
     95      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    9896      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
    101       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     97      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
     98      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
     100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     101      ! 
     102      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
    102103      !!---------------------------------------------------------------------- 
    103104      ! 
    104105      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    105106      ! 
    106       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
    107       !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 
    108       !     ...this should not affect too much the stability... Was ok on the tests we did... 
     107      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     108      DO jl = 1, jpl 
     109         DO jj = 2, jpjm1 
     110            DO ji = fs_2, fs_jpim1 
     111               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     112                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     113                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     114                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     115               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     116                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     117                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     118                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     119               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     120                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     121                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     122                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     123            END DO 
     124         END DO 
     125      END DO 
     126      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
     127      ! 
     128      ! 
     129      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     130      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     131      !              this should not affect too much the stability 
    109132      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
    110133      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    … …  
    116139      ELSE                         ;   icycle = 1 
    117140      ENDIF 
    118        
    119141      zdt = rdt_ice / REAL(icycle) 
    120142 
    … …  
    122144      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
    123145      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
    124  
     146      ! 
     147      ! setup transport for each ice cat  
     148      DO jl = 1, jpl 
     149         zu_cat(:,:,jl) = zudy(:,:) 
     150         zv_cat(:,:,jl) = zvdx(:,:) 
     151      END DO 
     152      ! 
    125153      ! --- define velocity for advection: u*grad(H) --- ! 
    126154      DO jj = 2, jpjm1 
    … …  
    154182         END WHERE 
    155183         ! 
    156          ! set u*a=u for advection of A only  
    157          DO jl = 1, jpl 
    158             zua_ho(:,:,jl) = zudy(:,:) 
    159             zva_ho(:,:,jl) = zvdx(:,:) 
    160          END DO 
    161           
     184         ! setup a mask where advection will be upstream 
     185         IF( ll_neg ) THEN 
     186            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) )  
     187            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) )  
     188            DO jl = 1, jpl 
     189               DO jj = 1, jpjm1 
     190                  DO ji = 1, jpim1 
     191                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 
     192                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0 
     193                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF 
     194                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 
     195                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0 
     196                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF 
     197                  END DO 
     198               END DO 
     199            END DO 
     200         ENDIF 
     201         ! 
     202         ! ----------------------- ! 
     203         ! ==> start advection <== ! 
     204         ! ----------------------- ! 
     205         ! 
     206         !== Ice area ==! 
    162207         zamsk = 1._wp 
    163          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area 
    164          zamsk = 0._wp 
    165          ! 
    166          zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
    167          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                !-- Ice volume 
    168          ! 
    169          zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
    170          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                !-- Snw volume 
    171          ! 
    172          zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
    173          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               !-- Salt content 
    174          ! 
    175          DO jk = 1, nlay_i 
    176             zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
    177             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   !-- Ice heat content 
    178          END DO 
    179          ! 
    180          DO jk = 1, nlay_s 
    181             zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
    182             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   !-- Snw heat content 
    183          END DO 
    184          ! 
    185          IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age 
    186             ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 
    187             !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that 
    188             !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration 
    189             !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
    190             !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 
    191             ! set u*a=u for advection of ice age 
    192             DO jl = 1, jpl 
    193                zua_ho(:,:,jl) = zudy(:,:) 
    194                zva_ho(:,:,jl) = zvdx(:,:) 
    195             END DO 
     208         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, & 
     209            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 
     210         ! 
     211         !                             ! --------------------------------- ! 
     212         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- ! 
     213            !                          ! --------------------------------- ! 
     214            zamsk = 0._wp 
     215            !== Ice volume ==! 
     216            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     217            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     218               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     219            !== Snw volume ==!          
     220            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     221            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     222               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
     223            ! 
    196224            zamsk = 1._wp 
    197             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i ) 
     225            !== Salt content ==! 
     226            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     227               &                                      psv_i, psv_i ) 
     228            !== Ice heat content ==! 
     229            DO jk = 1, nlay_i 
     230               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     231                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) ) 
     232            END DO 
     233            !== Snw heat content ==! 
     234            DO jk = 1, nlay_s 
     235               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     236                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) ) 
     237            END DO 
     238            ! 
     239            !                          ! ------------------------------------------ ! 
     240         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- ! 
     241            !                          ! ------------------------------------------ ! 
    198242            zamsk = 0._wp 
     243            !== Ice volume ==! 
     244            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     245            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     246               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     247            !== Snw volume ==!          
     248            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     249            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     250               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
     251            !== Salt content ==! 
     252            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
     253            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     254               &                                      zhvar, psv_i, zua_ups, zva_ups ) 
     255            !== Ice heat content ==! 
     256            DO jk = 1, nlay_i 
     257               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
     258               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     259                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 
     260            END DO 
     261            !== Snw heat content ==! 
     262            DO jk = 1, nlay_s 
     263               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
     264               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     265                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 
     266            END DO 
     267            ! 
     268            !                          ! ----------------------------------------- ! 
     269         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- ! 
     270            !                          ! ----------------------------------------- ! 
     271            zamsk = 0._wp 
     272            ! 
     273            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  & 
     274               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) ) 
     275            ! 
     276            ! inverse of Vi 
     277            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 
     278            ELSEWHERE                        ;   z1_vi(:,:,:) = 0. 
     279            END WHERE 
     280            ! inverse of Vs 
     281            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 
     282            ELSEWHERE                        ;   z1_vs(:,:,:) = 0. 
     283            END WHERE 
     284            ! 
     285            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 
     286            ! 
     287            !== Ice volume ==! 
     288            zuv_ups = zua_ups 
     289            zvv_ups = zva_ups 
     290            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     291            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     292               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     293            !== Salt content ==! 
     294            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 
     295            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 
     296               &                                      zhvar, psv_i, zuv_ups, zvv_ups ) 
     297            !== Ice heat content ==! 
     298            DO jk = 1, nlay_i 
     299               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 
     300               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     301                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 
     302            END DO 
     303            !== Snow volume ==!          
     304            zuv_ups = zua_ups 
     305            zvv_ups = zva_ups 
     306            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     307            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     308               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     309            !== Snw heat content ==! 
     310            DO jk = 1, nlay_s 
     311               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 
     312               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     313                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 
     314            END DO 
     315            ! 
     316            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs ) 
     317            ! 
    199318         ENDIF 
    200319         ! 
    201          IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds 
    202             ! set u*a=u for advection of Ap only  
    203             DO jl = 1, jpl 
    204                zua_ho(:,:,jl) = zudy(:,:) 
    205                zva_ho(:,:,jl) = zvdx(:,:) 
    206             END DO 
    207             ! 
     320         !== Ice age ==! 
     321         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 
    208322            zamsk = 1._wp 
    209             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction 
     323            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     324               &                                      poa_i, poa_i ) 
     325         ENDIF 
     326         ! 
     327         !== melt ponds ==! 
     328         IF ( ln_pnd_H12 ) THEN 
     329            ! fraction 
     330            zamsk = 1._wp 
     331            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & 
     332               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 
     333            ! volume 
    210334            zamsk = 0._wp 
    211             ! 
    212335            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 
    213             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! volume 
     336            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     337               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
    214338         ENDIF 
    215339         ! 
     340         !== Open water area ==! 
    216341         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    217342         DO jj = 2, jpjm1 
    218343            DO ji = fs_2, fs_jpim1 
    219                pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                              !-- Open water area 
     344               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &  
    220345                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    221346            END DO 
    222347         END DO 
    223          CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
    224          ! 
     348         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. ) 
     349         ! 
     350         ! 
     351         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
     352         ! Remove negative values (conservation is ensured) 
     353         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     354         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     355         ! 
     356         ! Make sure ice thickness is not too big 
     357         !    (because ice thickness can be too large where ice concentration is very small) 
     358         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     359 
    225360      END DO 
    226361      ! 
    … …  
    228363 
    229364    
    230    SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
     365   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  & 
     366      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 
    231367      !!---------------------------------------------------------------------- 
    232368      !!                  ***  ROUTINE adv_umx  *** 
    … …  
    235371      !!                 tracers and add it to the general trend of tracer equations 
    236372      !! 
    237       !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H 
     373      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 
    238374      !!                 - calculate tracer H at u and v points (Ultimate) 
    239       !!                 - calculate the high order fluxes using alterning directions (Macho?) 
     375      !!                 - calculate the high order fluxes using alterning directions (Macho) 
    240376      !!                 - apply a limiter on the fluxes (nonosc_ice) 
    241       !!                 - convert this tracer flux to a tracer content flux (uH -> uV) 
    242       !!                 - calculate the high order solution for tracer content V 
     377      !!                 - convert this tracer flux to a "volume" flux (uH -> uV) 
     378      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice) 
     379      !!                 - calculate the high order solution for V 
    243380      !! 
    244       !! ** Action : solve 2 equations => a) da/dt = -div(ua) 
    245       !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 
    246       !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H. 
    247       !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u 
    248       !!                             where ua is the flux from eq. a) 
    249       !!                        - at last we estimate dV/dt = -div(uH*ua/u) 
     381      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA) 
     382      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H) 
     383      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 
    250384      !! 
    251       !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 
    252       !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 
    253       !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
    254       !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     385      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 
     386      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u 
     387      !!                             where uA is the flux from eq. a) 
     388      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 
     389      !!                        - at last we estimate dV/dt = -div(uH * uA / u) 
     390      !! 
     391      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u) 
     392      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)  
     393      !! 
     394      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 
     395      !!           - At the ice edge, Ultimate scheme can lead to: 
     396      !!                              1) negative interpolated tracers at u-v points 
     397      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 
     398      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
     399      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     400      !!                              Solution for 2): we set it to 0 in this case 
    255401      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 
    256402      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 
    257       !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D. 
     403      !!             work on H (and not V). It is partly related to the multi-category approach 
    258404      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    259       !!             concentration is small). 
    260       !! To-do: expand the prelimiter from zalesak to make it work in 2D 
    261       !!---------------------------------------------------------------------- 
    262       REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
    263       INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
    264       INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration 
    265       INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration 
    266       REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step 
    267       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
    268       REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
    269       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    270       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field 
    271       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field 
    272       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
     405      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
     406      !!             since sv_i and e_i are still good. 
     407      !!---------------------------------------------------------------------- 
     408      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     409      INTEGER                         , INTENT(in   )           ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     410      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration 
     411      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration 
     412      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step 
     413      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2 
     414      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u 
     415      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity 
     416      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field 
     417      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field 
     418      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes 
     419      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes 
    273420      ! 
    274421      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
    … …  
    303450            DO jj = 1, jpjm1 
    304451               DO ji = 1, fs_jpim1 
    305                   IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
    306                      zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
    307                      zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     452                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
     453                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     454                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 
    308455                  ELSE 
    309456                     zfu_ho (ji,jj,jl) = 0._wp 
    … …  
    311458                  ENDIF 
    312459                  ! 
    313                   IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
    314                      zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
    315                      zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     460                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
     461                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     462                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 
    316463                  ELSE 
    317464                     zfv_ho (ji,jj,jl) = 0._wp   
    … …  
    321468            END DO 
    322469         END DO 
     470 
     471         ! the new "volume" fluxes must also be "flux corrected" 
     472         ! thus we calculate the upstream solution and apply a limiter again 
     473         DO jl = 1, jpl 
     474            DO jj = 2, jpjm1 
     475               DO ji = fs_2, fs_jpim1 
     476                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 
     477                  ! 
     478                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 
     479               END DO 
     480            END DO 
     481         END DO 
     482         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. ) 
     483         ! 
     484         IF    ( np_limiter == 1 ) THEN 
     485            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
     486         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
     487            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 
     488            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 
     489         ENDIF 
     490         ! 
    323491      ENDIF 
    324       !                                   --ho 
    325       ! in case of advection of A: output u*a 
    326       ! ------------------------------------- 
     492      !                                   --ho    --ups 
     493      ! in case of advection of A: output u*a and u*a 
     494      ! ----------------------------------------------- 
    327495      IF( PRESENT( pua_ho ) ) THEN 
    328496         DO jl = 1, jpl 
    329497            DO jj = 1, jpjm1 
    330498               DO ji = 1, fs_jpim1 
    331                   pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 
    332                   pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 
    333                END DO 
     499                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     500                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     501              END DO 
    334502            END DO 
    335503         END DO 
    … …  
    499667         END DO 
    500668         ! 
    501          IF    ( kn_limiter == 1 ) THEN 
     669         IF    ( np_limiter == 1 ) THEN 
    502670            CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    503          ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     671         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
    504672            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    505673            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    … …  
    517685               END DO 
    518686            END DO 
    519             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     687            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    520688 
    521689            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    … …  
    538706               END DO 
    539707            END DO 
    540             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     708            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    541709 
    542710         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    … …  
    549717               END DO 
    550718            END DO 
    551             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     719            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    552720            ! 
    553721            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    … …  
    570738               END DO 
    571739            END DO 
    572             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     740            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    573741 
    574742         ENDIF 
    575          IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     743         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    576744          
    577745      ENDIF 
    … …  
    609777         ! 
    610778         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    611          CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
     779         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
    612780         !                                                        !--  limiter in x --! 
    613          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     781         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    614782         !                                                        !--  advective form update in zpt  --! 
    615783         DO jl = 1, jpl 
    … …  
    619787                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
    620788                     &                                                                                        * pamsk           & 
    621                      &                             ) * pdt ) * tmask(ji,jj,1)  
     789                     &                             ) * pdt ) * tmask(ji,jj,1) 
    622790               END DO 
    623791            END DO 
    … …  
    627795         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    628796         IF( ll_hoxy ) THEN 
    629             CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
     797            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
    630798         ELSE 
    631             CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
     799            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
    632800         ENDIF 
    633801         !                                                        !--  limiter in y --! 
    634          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     802         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    635803         !          
    636804         ! 
    … …  
    638806         ! 
    639807         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    640          CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
     808         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
    641809         !                                                        !--  limiter in y --! 
    642          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     810         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    643811         !                                                        !--  advective form update in zpt  --! 
    644812         DO jl = 1, jpl 
    … …  
    656824         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    657825         IF( ll_hoxy ) THEN 
    658             CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
     826            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
    659827         ELSE 
    660             CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
     828            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
    661829         ENDIF 
    662830         !                                                        !--  limiter in x --! 
    663          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     831         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    664832         ! 
    665833      ENDIF 
    666834 
    667       IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     835      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    668836      ! 
    669837   END SUBROUTINE macho 
    670838 
    671839 
    672    SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
     840   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    673841      !!--------------------------------------------------------------------- 
    674842      !!                    ***  ROUTINE ultimate_x  *** 
    … …  
    680848      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    681849      !!---------------------------------------------------------------------- 
     850      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    682851      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    683852      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    … …  
    806975            DO jj = 1, jpjm1 
    807976               DO ji = 1, fs_jpim1 
    808                   IF( pt_u(ji,jj,jl) < 0._wp ) THEN 
     977                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    809978                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    810979                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    … …  
    826995    
    827996  
    828    SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
     997   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    829998      !!--------------------------------------------------------------------- 
    830999      !!                    ***  ROUTINE ultimate_y  *** 
    … …  
    8361005      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    8371006      !!---------------------------------------------------------------------- 
     1007      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    8381008      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    8391009      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    … …  
    9591129            DO jj = 1, jpjm1 
    9601130               DO ji = 1, fs_jpim1 
    961                   IF( pt_v(ji,jj,jl) < 0._wp ) THEN 
     1131                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    9621132                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    9631133                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    10231193      !                        |      |      |        |    * 
    10241194      !            t_ups :       i-1     i       i+1       i+2    
    1025       IF( ll_prelimiter_zalesak ) THEN 
     1195      IF( ll_prelim ) THEN 
    10261196          
    10271197         DO jl = 1, jpl 
    … …  
    11021272               ! 
    11031273               !                                  ! up & down beta terms 
    1104                IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
    1105                ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig 
     1274               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 
     1275               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1276               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig 
    11061277               ENDIF 
    11071278               ! 
    1108                IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    1109                ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
     1279               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1280               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
    11101281               ENDIF 
    11111282               ! 
    … …  
    11491320         END DO 
    11501321 
    1151          ! clem test 
    1152 !!         DO jj = 2, jpjm1 
    1153 !!            DO ji = 2, fs_jpim1   ! vector opt. 
    1154 !!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1155 !!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1156 !!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1157 !!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1158 !!                  &         ) * tmask(ji,jj,1) 
    1159 !!               IF( zzt < -epsi20 ) THEN 
    1160 !!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt 
    1161 !!               ENDIF 
    1162 !!            END DO 
    1163 !!         END DO 
    1164  
    11651322      END DO 
    11661323      ! 
    … …  
    12031360               Rjp = zslpx(ji+1,jj,jl) 
    12041361 
    1205                IF( kn_limiter == 3 ) THEN 
     1362               IF( np_limiter == 3 ) THEN 
    12061363 
    12071364                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    … …  
    12191376                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
    12201377 
    1221                ELSEIF( kn_limiter == 2 ) THEN 
     1378               ELSEIF( np_limiter == 2 ) THEN 
    12221379                  IF( Rj /= 0. ) THEN 
    12231380                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
    … …  
    12981455               Rjp = zslpy(ji,jj+1,jl) 
    12991456 
    1300                IF( kn_limiter == 3 ) THEN 
     1457               IF( np_limiter == 3 ) THEN 
    13011458 
    13021459                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    … …  
    13141471                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
    13151472 
    1316                ELSEIF( kn_limiter == 2 ) THEN 
     1473               ELSEIF( np_limiter == 2 ) THEN 
    13171474 
    13181475                  IF( Rj /= 0. ) THEN 
    … …  
    13581515   END SUBROUTINE limiter_y 
    13591516 
     1517 
     1518   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     1519      !!------------------------------------------------------------------- 
     1520      !!                  ***  ROUTINE Hbig  *** 
     1521      !! 
     1522      !! ** Purpose : Thickness correction in case advection scheme creates 
     1523      !!              abnormally tick ice or snow 
     1524      !! 
     1525      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     1526      !!                 (before advection) and reduce it by adapting ice concentration 
     1527      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     1528      !!                 (before advection) and reduce it by sending the excess in the ocean 
     1529      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     1530      !!                 and reduce it by sending the excess in the ocean 
     1531      !!              4- correct pond fraction to avoid a_ip > a_i 
     1532      !! 
     1533      !! ** input   : Max thickness of the surrounding 9-points 
     1534      !!------------------------------------------------------------------- 
     1535      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
     1536      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
     1537      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 
     1538      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     1539      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1540      ! 
     1541      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1542      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra 
     1543      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
     1544      !!------------------------------------------------------------------- 
     1545      ! 
     1546      z1_dt = 1._wp / pdt 
     1547      ! 
     1548      DO jl = 1, jpl 
     1549 
     1550         DO jj = 1, jpj 
     1551            DO ji = 1, jpi 
     1552               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1553                  ! 
     1554                  !                               ! -- check h_ip -- ! 
     1555                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     1556                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1557                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     1558                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     1559                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     1560                     ENDIF 
     1561                  ENDIF 
     1562                  ! 
     1563                  !                               ! -- check h_i -- ! 
     1564                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     1565                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     1566                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1567                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     1568                  ENDIF 
     1569                  ! 
     1570                  !                               ! -- check h_s -- ! 
     1571                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     1572                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     1573                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1574                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     1575                     ! 
     1576                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
     1577                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1578                     ! 
     1579                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1580                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     1581                  ENDIF            
     1582                  ! 
     1583                  !                               ! -- check snow load -- ! 
     1584                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     1585                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     1586                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
     1587                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     1588                  IF( zvs_excess > 0._wp ) THEN 
     1589                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     1590                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     1591                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1592                     ! 
     1593                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1594                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
     1595                  ENDIF 
     1596                   
     1597               ENDIF 
     1598            END DO 
     1599         END DO 
     1600      END DO  
     1601      !                                           !-- correct pond fraction to avoid a_ip > a_i 
     1602      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
     1603      ! 
     1604      ! 
     1605   END SUBROUTINE Hbig 
     1606    
    13601607#else 
    13611608   !!---------------------------------------------------------------------- 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rdgrft.F90

    r10888 r11083  
    142142      INTEGER, PARAMETER ::   jp_itermax = 20     
    143143      !!------------------------------------------------------------------- 
    144       ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 
    145       !       likely due to truncation error ( i.e. 1. - 1. /= 0 ) 
    146       !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    147        
    148144      ! controls 
    149145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
    … …  
    156152      ENDIF       
    157153 
    158       CALL ice_var_zapsmall   ! Zero out categories with very small areas 
    159  
    160154      !-------------------------------- 
    161155      ! 0) Identify grid cells with ice 
    162156      !-------------------------------- 
     157      at_i(:,:) = SUM( a_i, dim=3 ) 
     158      ! 
    163159      npti = 0   ;   nptidx(:) = 0 
    164160      ipti = 0   ;   iptidx(:) = 0 
    165161      DO jj = 1, jpj 
    166162         DO ji = 1, jpi 
    167             IF ( at_i(ji,jj) > 0._wp ) THEN 
     163            IF ( at_i(ji,jj) > epsi10 ) THEN 
    168164               npti           = npti + 1 
    169165               nptidx( npti ) = (jj - 1) * jpi + ji 
    … …  
    178174         
    179175         ! just needed here 
    180          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 
    181          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    182178         ! needed here and in the iteration loop 
    183          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) ) 
    184          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) ) 
    185          CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) ) 
     179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i  ) 
     180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  ) 
     181         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
    186182 
    187183         DO ji = 1, npti 
    … …  
    310306 
    311307      !                       ! Ice thickness needed for rafting 
    312       WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    313       ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp 
     308      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     309      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    314310      END WHERE 
    315311 
    … …  
    329325      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    330326      ! 
    331       WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    332       ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
     327      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     328      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    333329      END WHERE 
    334330      ! 
    … …  
    455451      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    456452      ! NOTE: 0 < aksum <= 1 
    457       WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    458       ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp 
     453      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     454      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    459455      END WHERE 
    460456       
    … …  
    466462         DO ji = 1, npti 
    467463            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    468             IF( zfac > pa_i(ji,jl) ) THEN 
     464            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 
    469465               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
    470466            ENDIF 
    … …  
    510506      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    511507      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a 
     508      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i) 
    512509      ! 
    513510      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice 
    … …  
    518515      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
    519516      !!------------------------------------------------------------------- 
    520  
     517      ! 
     518      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume 
     519      ! 
    521520      ! 1) Change in open water area due to closing and opening 
    522521      !-------------------------------------------------------- 
    … …  
    535534            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    536535 
    537                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    538  
     536               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     537               ELSE                                 ;   z1_ai(ji) = 0._wp 
     538               ENDIF 
     539                
    539540               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    540541               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
    … …  
    549550 
    550551               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    551                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     552               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg 
     553               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0 
     554               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 
     555               ENDIF 
    552556               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    553557 
    554558               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    555559               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    556                virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     560               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw 
    557561               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg 
    558562               sirdg1     = sv_i_2d(ji,jl1)   * afrdg 
    … …  
    726730      END DO ! jl1 
    727731      ! 
     732      ! roundoff errors 
     733      !---------------- 
    728734      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    729       WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp 
    730       WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp 
     735      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
    731736      ! 
    732737   END SUBROUTINE rdgrft_shift 
    … …  
    854859         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    855860         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
    856  
     861         ! 
    857862         !                 !---------------------! 
    858863      CASE( 2 )            !==  from 1D to 2D  ==! 
    … …  
    945950         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    946951      ENDIF 
    947       !                              ! allocate tke arrays 
     952      ! 
     953      IF( .NOT. ln_icethd ) THEN 
     954         rn_porordg = 0._wp 
     955         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 
     956         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 
     957         IF( lwp ) THEN 
     958            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed' 
     959            WRITE(numout,*) '            rn_porordg   = ', rn_porordg 
     960            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg  
     961            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg  
     962            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft  
     963            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft  
     964         ENDIF 
     965      ENDIF 
     966      !                              ! allocate arrays 
    948967      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 
    949968      ! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg.F90

    r10888 r11083  
    5858      !!-------------------------------------------------------------------- 
    5959      INTEGER, INTENT(in) ::   kt     ! ice time step 
     60      ! 
     61      INTEGER  ::   jl   ! dummy loop indices 
    6062      !!-------------------------------------------------------------------- 
    6163      ! controls 
    … …  
    6870         WRITE(numout,*)'~~~~~~~~~~~' 
    6971      ENDIF 
    70  
    71       ! -------- 
    72       ! Rheology 
    73       ! --------    
     72      ! 
     73      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
     74         tau_icebfr(:,:) = 0._wp 
     75         DO jl = 1, jpl 
     76            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     77         END DO 
     78      ENDIF 
     79      ! 
     80      !--------------! 
     81      !== Rheology ==! 
     82      !--------------!    
    7483      SELECT CASE( nice_rhg ) 
    7584      !                                !------------------------! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rhg_evp.F90

    r10888 r11083  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       LOGICAL, PARAMETER ::   ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T) 
     114      LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 
    115115      !                                                                              or only at the main time step (F) 
    116116      INTEGER ::   ji, jj       ! dummy loop indices 
    … …  
    160160 
    161161      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    162       REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     162      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
     163      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    163164      !! --- diags 
    164165      REAL(wp), DIMENSION(jpi,jpj) ::   zswi 
    … …  
    319320 
    320321            ! switches 
    321             zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
    322             zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
     322            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp 
     323            ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF 
     324            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp 
     325            ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF 
    323326 
    324327         END DO 
    … …  
    519522                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    520523                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    521                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     524                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    522525                     &           ) * zmaskV(ji,jj) 
    523526                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    526529                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    527530                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    528                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     531                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    529532                     &            ) * zmaskV(ji,jj) 
    530533                  ENDIF 
    … …  
    567570                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    568571                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    569                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     572                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    570573                     &            ) * zmaskU(ji,jj) 
    571574                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    574577                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    575578                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    576                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     579                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    577580                     &            ) * zmaskU(ji,jj) 
    578581                  ENDIF 
    … …  
    617620                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    618621                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    619                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     622                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    620623                     &            ) * zmaskU(ji,jj) 
    621624                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    624627                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    625628                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    626                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     629                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    627630                     &            ) * zmaskU(ji,jj) 
    628631                  ENDIF 
    … …  
    665668                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    666669                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    667                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     670                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    668671                     &           ) * zmaskV(ji,jj) 
    669672                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    672675                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    673676                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    674                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     677                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    675678                     &            ) * zmaskV(ji,jj) 
    676679                  ENDIF 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/iceitd.F90

    r10888 r11083  
    2121   USE ice1D          ! sea-ice: thermodynamic variables 
    2222   USE ice            ! sea-ice: variables 
     23   USE icevar         ! sea-ice: operations 
    2324   USE icectl         ! sea-ice: conservation tests 
    2425   USE icetab         ! sea-ice: convert 1D<=>2D 
    … …  
    9192      !  1) Identify grid cells with ice 
    9293      !----------------------------------------------------------------------------------------------- 
     94      at_i(:,:) = SUM( a_i, dim=3 ) 
     95      ! 
    9396      npti = 0   ;   nptidx(:) = 0 
    9497      DO jj = 1, jpj 
    … …  
    249252            ! --- g(h) for each thickness category --- !   
    250253            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
    251                &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR   (1:npti,jl)   )   ! out 
     254               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out 
    252255            ! 
    253256         END DO 
    … …  
    389392      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary 
    390393      ! 
    391       INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices 
    392       INTEGER  ::   ii, ij, jl2, jl1   ! local integers 
     394      INTEGER  ::   ji, jl, jk         ! dummy loop indices 
     395      INTEGER  ::   jl2, jl1           ! local integers 
    393396      REAL(wp) ::   ztrans             ! ice/snow transferred 
    394       REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace 
    395       REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    - 
     397      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace 
     398      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    - 
     399      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d 
     400      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d 
    396401      !!------------------------------------------------------------------ 
    397402          
    … …  
    405410      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    406411      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     412      DO jl = 1, jpl 
     413         DO jk = 1, nlay_s 
     414            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     415         END DO 
     416         DO jk = 1, nlay_i 
     417            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     418         END DO 
     419      END DO 
     420      ! to correct roundoff errors on a_i 
     421      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d ) 
    407422 
    408423      !---------------------------------------------------------------------------------------------- 
    … …  
    435450               ELSE                                  ;   zworka(ji) = 0._wp 
    436451               ENDIF 
    437                ! 
    438                ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
    439                !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
    440                !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    441452               ! 
    442453               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
    … …  
    476487         ! 
    477488         DO jk = 1, nlay_s         !--- Snow heat content 
    478             ! 
    479489            DO ji = 1, npti 
    480                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    481                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    482490               ! 
    483491               jl1 = kdonor(ji,jl) 
    … …  
    487495                  ELSE                ;  jl2 = jl 
    488496                  ENDIF 
    489                   ! 
    490                   ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji) 
    491                   e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 
    492                   e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
     497                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji) 
     498                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans 
     499                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans 
    493500               ENDIF 
    494501            END DO 
    … …  
    497504         DO jk = 1, nlay_i         !--- Ice heat content 
    498505            DO ji = 1, npti 
    499                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    500                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    501506               ! 
    502507               jl1 = kdonor(ji,jl) 
    … …  
    506511                  ELSE                ;  jl2 = jl 
    507512                  ENDIF 
    508                   ! 
    509                   ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji) 
    510                   e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 
    511                   e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 
     513                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji) 
     514                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans 
     515                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans 
    512516               ENDIF 
    513517            END DO 
    … …  
    515519         ! 
    516520      END DO                   ! boundaries, 1 to jpl-1 
     521 
     522      !------------------- 
     523      ! 3) roundoff errors 
     524      !------------------- 
     525      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
     526      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
     527      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     528 
     529      ! at_i must be <= rn_amax 
     530      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 
     531      DO jl  = 1, jpl 
     532         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   & 
     533            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 
     534      END DO 
    517535       
    518536      !------------------------------------------------------------------------------- 
    519       ! 3) Update ice thickness and temperature 
     537      ! 4) Update ice thickness and temperature 
    520538      !------------------------------------------------------------------------------- 
    521539      WHERE( a_i_2d(1:npti,:) >= epsi20 ) 
    … …  
    536554      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    537555      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     556      DO jl = 1, jpl 
     557         DO jk = 1, nlay_s 
     558            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     559         END DO 
     560         DO jk = 1, nlay_i 
     561            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     562         END DO 
     563      END DO 
    538564      ! 
    539565   END SUBROUTINE itd_shiftice 
    … …  
    558584      ! 
    559585      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
     586      ! 
     587      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    560588      ! 
    561589      jdonor(:,:) = 0 
    … …  
    635663      END DO 
    636664      ! 
     665      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     666      ! 
    637667   END SUBROUTINE ice_itd_reb 
    638668 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icestp.F90

    r10888 r11083  
    189189         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    190190         ! 
    191          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections 
     191                                        CALL ice_cor( kt , 2 )        ! -- Corrections 
    192192         ! 
    193193                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling) 
    … …  
    323323         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s 
    324324      ENDIF 
     325      !                                        !--- change max ice concentration for roundoff errors 
     326      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 ) 
     327      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 ) 
    325328      !                                        !--- check consistency 
    326329      IF ( jpl > 1 .AND. ln_virtual_itd ) THEN 
    … …  
    431434      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface 
    432435 
    433       tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    434       cnd_ice   (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
    435       qtr_ice_bot(:,:,:) = 0._wp  ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
     436      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     437      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
     438      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 
     439      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
     440      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F) 
    436441      ! 
    437442      ! for control checks (ln_icediachk) 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd.F90

    r10888 r11083  
    102102      ENDIF 
    103103       
    104       CALL ice_var_glo2eqv 
    105  
    106104      !---------------------------------------------! 
    107105      ! computation of friction velocity at T points 
    … …  
    162160            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    163161 
    164             ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    165             IF( zqld > 0._wp ) THEN 
     162            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
     163            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
     164            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    166165               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    167166               qlead(ji,jj) = 0._wp 
    … …  
    178177      ! In case we bypass open-water ice formation 
    179178      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp 
    180       ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere 
     179      ! In case we bypass growing/melting from top and bottom 
    181180      IF( .NOT. ln_icedH ) THEN 
    182          qt_atm_oi  (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) 
    183181         qsb_ice_bot(:,:) = 0._wp 
    184182         fhld       (:,:) = 0._wp 
    … …  
    221219            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
    222220            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    223             ! 
    224             IF( ln_icedH ) THEN                                     ! --- growing/melting --- ! 
    225                               CALL ice_thd_zdf                             ! Ice/Snow Temperature profile 
    226                               CALL ice_thd_dh                              ! Ice/Snow thickness    
    227                               CALL ice_thd_pnd                             ! Melt ponds formation 
    228                               CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping 
     221            !                                       
     222                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
     223            ! 
     224            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
     225                              CALL ice_thd_dh                           ! Ice-Snow thickness    
     226                              CALL ice_thd_pnd                          ! Melt ponds formation 
     227                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    229228            ENDIF 
    230             ! 
    231229                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !     
    232230            ! 
    233                               CALL ice_thd_temp                     ! --- temperature update --- ! 
     231                              CALL ice_thd_temp                     ! --- Temperature update --- ! 
    234232            ! 
    235233            IF( ln_icedH .AND. ln_virtual_itd ) & 
    236                &              CALL ice_thd_mono                     ! --- extra lateral melting if virtual_itd --- ! 
    237             ! 
    238             IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- ! 
     234               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- ! 
     235            ! 
     236            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- ! 
    239237            ! 
    240238                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 
    241239            !                                                       ! --- & Move to 2D arrays --- ! 
    242             ! 
    243240         ENDIF 
    244241         ! 
    245242      END DO 
    246  
     243      ! 
    247244      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    248       ! 
    249                            CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- ! 
    250       !                                                    !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0 
    251245      !                    
    252       IF( jpl > 1      )   CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- ! 
    253       ! 
    254       IF( ln_icedO     )   CALL ice_thd_do                 ! --- frazil ice growing in leads --- ! 
     246      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
     247      ! 
     248      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
    255249      ! 
    256250      ! controls 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_do.F90

    r10888 r11083  
    114114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 
    115115 
    116       CALL ice_var_agg(1) 
    117       CALL ice_var_glo2eqv 
    118  
     116      at_i(:,:) = SUM( a_i, dim=3 ) 
    119117      !------------------------------------------------------------------------------! 
    120118      ! 1) Collection thickness of ice formed in leads and polynyas 
    … …  
    319317 
    320318         ! --- lateral ice growth --- ! 
    321          ! If lateral ice growth gives an ice concentration gt 1, then 
     319         ! If lateral ice growth gives an ice concentration > amax, then 
    322320         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    323321         DO ji = 1, npti 
    324             IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 
    325                zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) 
     322            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error 
     323               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) 
    326324               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    327                za_newice(ji) = za_newice(ji) - zda_res  (ji) 
    328                zv_newice(ji) = zv_newice(ji) - zdv_res  (ji) 
     325               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) ) 
     326               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) ) 
    329327            ELSE 
    330328               zda_res(ji) = 0._wp 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icethd_zdf_bl99.F90

    r10888 r11083  
    206206      ! 
    207207      l_T_converged(:) = .FALSE. 
    208       !                                                          !============================! 
    209208      ! Convergence calculated until all sub-domain grid points have converged 
    210209      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 
    211210      ! but values are not taken into account (results independant of MPI partitioning) 
    212211      ! 
     212      !                                                                            !============================! 
    213213      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    214          !                                                       !============================! 
     214         !                                                                         !============================! 
    215215         iconv = iconv + 1 
    216216         ! 
    … …  
    742742                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
    743743                  ! t_i 
    744                   DO jk = 0, nlay_i 
     744                  DO jk = 1, nlay_i 
    745745                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    746746                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    … …  
    856856         t_i_1d    (1:npti,:) = ztiold        (1:npti,:) 
    857857         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti) 
     858 
     859         !!clem 
     860         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input 
     861         !clem 
    858862      ENDIF 
    859863      ! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icevar.F90

    r10888 r11083  
    4444   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4545   !!   ice_var_zapsmall  : remove very small area and volume 
    46    !!   ice_var_zapneg    : remove negative ice fields (to debug the advection scheme UM3-5) 
     46   !!   ice_var_zapneg    : remove negative ice fields 
     47   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    4748   !!   ice_var_itd       : convert 1-cat to jpl-cat 
    4849   !!   ice_var_itd2      : convert N-cat to jpl-cat 
    … …  
    7172   PUBLIC   ice_var_zapsmall 
    7273   PUBLIC   ice_var_zapneg 
     74   PUBLIC   ice_var_roundoff 
    7375   PUBLIC   ice_var_itd 
    7476   PUBLIC   ice_var_itd2 
    … …  
    229231                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    230232                     ! 
    231                      ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
     233                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3] 
    232234                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C] 
    233235                     ! Conversion q(S,T) -> T (second order equation) 
    … …  
    236238                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    237239                     ! 
    238                   ELSE                                !--- no ice 
     240                  ELSE                                   !--- no ice 
    239241                     t_i(ji,jj,jk,jl) = rt0 
    240242                  ENDIF 
    … …  
    537539 
    538540 
    539    SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     541   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    540542      !!------------------------------------------------------------------- 
    541543      !!                   ***  ROUTINE ice_var_zapneg *** 
    … …  
    543545      !! ** Purpose :   Remove negative sea ice fields and correct fluxes 
    544546      !!------------------------------------------------------------------- 
    545       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    546       ! 
     547      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    547548      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    548549      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
    … …  
    555556      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    556557      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    557       !!------------------------------------------------------------------- 
    558       ! 
     558      ! 
     559      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
     560      REAL(wp) ::   z1_dt 
     561      !!------------------------------------------------------------------- 
     562      ! 
     563      z1_dt = 1._wp / pdt 
    559564      ! 
    560565      DO jl = 1, jpl       !==  loop over the categories  ==! 
    561566         ! 
     567         ! make sure a_i=0 where v_i<=0 
     568         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp 
     569 
    562570         !---------------------------------------- 
    563571         ! zap ice energy and send it to the ocean 
    … …  
    566574            DO jj = 1 , jpj 
    567575               DO ji = 1 , jpi 
    568                   IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    569                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 >0 
     576                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     577                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 
    570578                     pe_i(ji,jj,jk,jl) = 0._wp 
    571579                  ENDIF 
    … …  
    577585            DO jj = 1 , jpj 
    578586               DO ji = 1 , jpi 
    579                   IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    580                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     587                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     588                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 
    581589                     pe_s(ji,jj,jk,jl) = 0._wp 
    582590                  ENDIF 
    … …  
    590598         DO jj = 1 , jpj 
    591599            DO ji = 1 , jpi 
    592                IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    593                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice 
     600               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     601                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 
    594602                  pv_i   (ji,jj,jl) = 0._wp 
    595603               ENDIF 
    596                IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    597                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice 
     604               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     605                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 
    598606                  pv_s   (ji,jj,jl) = 0._wp 
    599607               ENDIF 
    600                IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    601                   sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice 
     608               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     609                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 
    602610                  psv_i  (ji,jj,jl) = 0._wp 
    603611               ENDIF 
    … …  
    616624   END SUBROUTINE ice_var_zapneg 
    617625 
     626 
     627   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     628      !!------------------------------------------------------------------- 
     629      !!                   ***  ROUTINE ice_var_roundoff *** 
     630      !! 
     631      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors 
     632      !!------------------------------------------------------------------- 
     633      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration 
     634      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     635      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume 
     636      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content 
     637      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content 
     638      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
     639      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     640      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
     641      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     642      !!------------------------------------------------------------------- 
     643      ! 
     644      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0 
     645      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0 
     646      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0 
     647      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0 
     648      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0 
     649      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
     650      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
     651      IF ( ln_pnd_H12 ) THEN 
     652         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
     653         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     654      ENDIF 
     655      ! 
     656   END SUBROUTINE ice_var_roundoff 
     657    
    618658    
    619659   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    … …  
    650690      INTEGER  :: idim, i_fill, jl0   
    651691      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    652       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    653       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables 
     692      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables 
     693      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    654694      INTEGER , DIMENSION(4)                  ::   itest 
    655695      !!------------------------------------------------------------------- 
    … …  
    780820      !! 
    781821      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    782       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     822       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    783823      !!               
    784824      !!               3) Expand the filling to the empty cat between jlmin and jlmax  

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icewri.F90

    r10888 r11083  
    145145 
    146146        ! record presence of fast ice 
    147          WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
     147         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
    148148         ELSEWHERE                                                ; zfast(:,:) = 0._wp 
    149149         END WHERE 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/NST/agrif_top_update.F90

    r10888 r11083  
    138138                     DO ji=i1,i2 
    139139                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    140                            trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
     140                           trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) &  
    141141                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    142142                                 &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90

    r10888 r11083  
    8585   ! 
    8686   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
    87    INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
     87   INTEGER, DIMENSION(jp_bdy) ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
    8888   INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
    8989   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90

    r10888 r11083  
    243243                        IF( ln_full_vel_array(jbdy) ) THEN 
    244244                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    245                                      & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   & 
     245                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy),   & 
    246246                                     & fvl=ln_full_vel_array(jbdy)  ) 
    247247                        ELSE 
    … …  
    313313                     jend = jstart + dta%nread(1) - 1 
    314314                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    315                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   & 
     315                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy),   & 
    316316                                  & fvl=ln_full_vel_array(jbdy) ) 
    317317                  ENDIF 
    … …  
    446446      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s 
    447447#endif 
    448       NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy 
     448      NAMELIST/nambdy_dta/ ln_full_vel 
    449449      !!--------------------------------------------------------------------------- 
    450450      ! 
    … …  
    508508      ! Read namelists 
    509509      ! -------------- 
    510       REWIND(numnam_ref) 
    511510      REWIND(numnam_cfg) 
    512511      jfld = 0  
    513512      DO jbdy = 1, nb_bdy          
    514513         IF( nn_dta(jbdy) == 1 ) THEN 
     514            REWIND(numnam_ref) 
    515515            READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 
    516516901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90

    r10888 r11083  
    187187         ! Use characteristics method instead 
    188188         zflag = ABS(flagu) 
    189          zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 
     189         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij) 
    190190         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
    191191      END DO 
    … …  
    205205         ! Use characteristics method instead 
    206206         zflag = ABS(flagv) 
    207          zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 
     207         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv)) 
    208208         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
    209209      END DO 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90

    r10888 r11083  
    5757      INTEGER ::   jbdy   ! BDY set index 
    5858      !!---------------------------------------------------------------------- 
    59       ! 
    60       IF( ln_timing )   CALL timing_start('bdy_ice_thd') 
     59      ! controls 
     60      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing 
     61      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    6162      ! 
    6263      CALL ice_var_glo2eqv 
    … …  
    7879      CALL ice_var_agg(1) 
    7980      ! 
    80       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
    81       IF( ln_timing )   CALL timing_stop('bdy_ice_thd') 
     81      ! controls 
     82      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     83      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
     84      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
    8285      ! 
    8386   END SUBROUTINE bdy_ice 
    … …  
    148151            jpbound = 0   ;   ib = ji   ;   jb = jj 
    149152            ! 
    150             IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1 ; jb = jj 
    151             IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1 ; jb = jj 
    152             IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; ib = ji  ; jb = jj+1 
    153             IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. )   jpbound = 1 ; ib = ji  ; jb = jj-1 
     153            IF( u_ice(ji  ,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1 
     154            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1 
     155            IF( v_ice(ji  ,jj  ) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; jb = jj+1 
     156            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; jb = jj-1 
    154157            ! 
    155158            IF( nn_ice_dta(jbdy) == 0 )   jpbound = 0 ; ib = ji ; jb = jj   ! case ice boundaries = initial conditions 
    … …  
    306309                     ! one of the two zmsk is always 0 (because of zflag) 
    307310                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
    308                      zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
     311                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )  ! 0 if no ice 
    309312                     !   
    310313                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice) 
    … …  
    329332                     ! one of the two zmsk is always 0 (because of zflag) 
    330333                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 
    331                      zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 
     334                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )  ! 0 if no ice 
    332335                     !   
    333336                     ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice) 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90

    r10888 r11083  
    140140      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    141141      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    142       INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending 
     142      INTEGER :: com_east, com_west, com_south, com_north, jpk_max  ! Flags for boundaries sending 
    143143      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving 
    144144      INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates 
    … …  
    397397          IF(lwp) WRITE(numout,*) 
    398398        ENDIF 
    399         IF( nb_jpk_bdy > 0 ) THEN 
     399        IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN 
    400400           IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***' 
    401401        ELSE 
    … …  
    516516         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    517517            &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    518  
    519          IF( nb_jpk_bdy>0 ) THEN 
    520             ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) ) 
    521             ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) ) 
    522             ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) ) 
    523          ELSE 
    524             ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
    525             ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO 
    526             ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO 
    527          ENDIF 
     518          
     519         jpk_max = MAXVAL(nb_jpk_bdy) 
     520         jpk_max = MAX(jpk_max, jpk) 
     521 
     522         ALLOCATE( dta_global(jpbdtau, 1, jpk_max) ) 
     523         ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO 
     524         ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO 
    528525 
    529526         IF ( icount>0 ) THEN 
    530             IF( nb_jpk_bdy>0 ) THEN 
    531                ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) ) 
    532                ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) ) 
    533                ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) ) 
    534             ELSE 
    535                ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
    536                ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO 
    537                ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO   
    538             ENDIF 
     527            ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) ) 
     528            ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO 
     529            ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO   
    539530         ENDIF 
    540531         !  
    … …  
    960951                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    961952                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    962                        if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then 
     953                       if( ii == (nlcit(nowe+1)-1) ) then 
    963954                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    964955                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    … …  
    974965                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    975966                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    976                        if((com_east_b .ne. 1) .and. (ii == 2)) then 
     967                       if( ii == 2 ) then 
    977968                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    978969                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    … …  
    989980                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    990981                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    991                        if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then 
     982                       if( ii == (nlcit(nowe+1)-1) ) then 
    992983                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    993984                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    … …  
    1004995                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    1005996                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    1006                        if((com_east_b .ne. 1) .and. (ii == 2)) then 
     997                       if( ii == 2 ) then 
    1007998                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    1008999                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90

    r10888 r11083  
    2929   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number  
    3030 
    31 !!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc ! 
    32 !!gm          I don't understand why. 
    33    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
    34 !!gm end 
    35  
    3631   PUBLIC   dia_cfl       ! routine called by step.F90 
    3732   PUBLIC   dia_cfl_init  ! routine called by nemogcm 
    … …  
    5550      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    5651      ! 
    57       INTEGER                ::   ji, jj, jk                            ! dummy loop indices 
    58       REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
    59       INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace 
    60 !!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     52      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices 
     53      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars 
     54      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace 
     55      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace 
    6156      !!---------------------------------------------------------------------- 
    6257      ! 
    … …  
    7166      DO jk = 1, jpk       ! calculate Courant numbers 
    7267         DO jj = 1, jpj 
    73             DO ji = 1, fs_jpim1   ! vector opt. 
     68            DO ji = 1, jpi 
    7469               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    7570               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
    … …  
    111106      !                    ! write out to file 
    112107      IF( lwp ) THEN 
    113          WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
     108         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
    114109         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 
    115110         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) 
    … …  
    172167      rCw_max = 0._wp 
    173168      ! 
    174 !!gm required to work 
    175       ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) ) 
    176 !!gm end 
    177       !       
    178169   END SUBROUTINE dia_cfl_init 
    179170 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90

    r10888 r11083  
    7474      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme  
    7575      ! 
    76       INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices 
    77       INTEGER  ::   ii, ifu, ib_bdy   ! local integers 
    78       INTEGER  ::   ij, ifv, igrd     !   -       - 
    79       REAL(wp) ::   zu, zv            ! local scalars 
     76      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices 
     77      INTEGER  ::   ifu, ifv, igrd, ib_bdy   ! local integers 
     78      REAL(wp) ::   zu, zv                   ! local scalars 
    8079      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke 
    8180      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv  
     81      REAL(wp)  :: zweightu, zweightv 
    8282      !!---------------------------------------------------------------------- 
    8383      ! 
    … …  
    9797       
    9898      zhke(:,:,jpk) = 0._wp 
    99        
    100       IF (ln_bdy) THEN 
    101          ! Maria Luneva & Fred Wobus: July-2016 
    102          ! compensate for lack of turbulent kinetic energy on liquid bdy points 
    103          DO ib_bdy = 1, nb_bdy 
    104             IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    105                igrd = 2           ! Copying normal velocity into points outside bdy 
    106                DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    107                   DO jk = 1, jpkm1 
    108                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    109                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    110                      ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
    111                      un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
    112                   END DO 
    113                END DO 
    114                ! 
    115                igrd = 3           ! Copying normal velocity into points outside bdy 
    116                DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    117                   DO jk = 1, jpkm1 
    118                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    119                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    120                      ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
    121                      vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
    122                   END DO 
    123                END DO 
    124             ENDIF 
    125          ENDDO   
    126       ENDIF  
    12799 
    128100      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
    … …  
    140112            END DO 
    141113         END DO 
     114         ! 
     115         IF (ln_bdy) THEN 
     116            ! Maria Luneva & Fred Wobus: July-2016 
     117            ! compensate for lack of turbulent kinetic energy on liquid bdy points 
     118            DO ib_bdy = 1, nb_bdy 
     119               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
     120                  igrd = 1           ! compensating null velocity on the bdy 
     121                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     122                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1 
     123                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1 
     124                     DO jk = 1, jpkm1 
     125                        zhke(ji,jj,jk) = 0._wp 
     126                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk) 
     127                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk) 
     128                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
     129                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
     130                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu)  
     131                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv)  
     132                     END DO 
     133                  END DO 
     134               END IF 
     135               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies 
     136            END DO 
     137         END IF 
    142138         ! 
    143139      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
    … …  
    158154            END DO 
    159155         END DO 
     156         IF (ln_bdy) THEN 
     157            ! Maria Luneva & Fred Wobus: July-2016 
     158            ! compensate for lack of turbulent kinetic energy on liquid bdy points 
     159            DO ib_bdy = 1, nb_bdy 
     160               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
     161                  igrd = 1           ! compensation null velocity on land at the bdy 
     162                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     163                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1 
     164                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1 
     165                     DO jk = 1, jpkm1 
     166                        zhke(ji,jj,jk) = 0._wp 
     167                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) & 
     168                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) )  
     169                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) & 
     170                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) ) 
     171                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
     172                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
     173                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
     174                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
     175                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
     176                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
     177                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
     178                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
     179                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu ) 
     180                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv ) 
     181                     END DO 
     182                  END DO 
     183               END IF 
     184            END DO 
     185         END IF 
    160186         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 
    161187         ! 
    162       END SELECT 
    163  
    164       IF (ln_bdy) THEN 
    165          ! restore velocity masks at points outside boundary 
    166          un(:,:,:) = un(:,:,:) * umask(:,:,:) 
    167          vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
    168       ENDIF       
    169  
     188      END SELECT  
    170189      ! 
    171190      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90

    r10888 r11083  
    297297         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
    298298         IF(lwp) WRITE(numout,*) '~~~~~ ' 
    299          ! 
    300          Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all) 
    301       ENDIF 
     299      ENDIF 
     300      ! 
    302301      ! 
    303302      DO jk = 1, jpkm1            ! calculate Courant numbers 
    … …  
    305304            DO ji = 2, fs_jpim1   ! vector opt. 
    306305               z1_e3w = 1._wp / e3w_n(ji,jj,jk) 
    307                Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    & 
    308                   &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
    309                   &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
    310                   &                        * r1_e1e2t(ji,jj)                                                  & 
    311                   &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
    312                   &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
    313                   &                        * r1_e1e2t(ji,jj)                                                  & 
    314                   &                      ) * z1_e3w 
     306               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &  ! 2*rdt and not r2dt (for restartability) 
     307                  &                             + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     308                  &                                 MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     309                  &                               * r1_e1e2t(ji,jj)                                                  & 
     310                  &                             + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     311                  &                                 MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     312                  &                               * r1_e1e2t(ji,jj)                                                  & 
     313                  &                             ) * z1_e3w 
    315314            END DO 
    316315         END DO 
    317316      END DO 
     317      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 
    318318      ! 
    319319      CALL iom_put("Courant",Cu_adv) 
    320320      ! 
    321       wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero 
    322321      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
    323322         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
    324             DO jj = 2, jpjm1                            ! w where necessary 
    325                DO ji = 2, fs_jpim1   ! vector opt. 
     323            DO jj = 1, jpj                              ! w where necessary 
     324               DO ji = 1, jpi 
    326325                  ! 
    327326                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 
    328327                  ! 
    329                   IF( zCu < Cu_min ) THEN               !<-- Fully explicit 
     328                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
    330329                     zcff = 0._wp 
    331330                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
    … …  
    346345      ELSE 
    347346         ! Fully explicit everywhere 
    348          Cu_adv = 0.0_wp                                ! Reuse array to output coefficient 
     347         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient 
     348         wi    (:,:,:) = 0._wp 
    349349      ENDIF 
    350350      CALL iom_put("wimp",wi)  

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90

    r10888 r11083  
    14801480      LOGICAL         , OPTIONAL, INTENT(in   ) ::   ld_lbc, ld_glb, ld_dlg 
    14811481      !! 
     1482      CHARACTER(len=128)                        ::   ccountname  ! name of a subroutine to count communications 
    14821483      LOGICAL ::   ll_lbc, ll_glb, ll_dlg 
    1483       INTEGER ::    ji,  jj,  jk,  jh, jf   ! dummy loop indices 
     1484      INTEGER ::    ji,  jj,  jk,  jh, jf, jcount   ! dummy loop indices 
    14841485      !!---------------------------------------------------------------------- 
    14851486      ! 
    … …  
    15381539         WRITE(numcom,*) ' ' 
    15391540         WRITE(numcom,*) ' lbc_lnk called' 
    1540          jj = 1 
    1541          DO ji = 2, n_sequence_lbc 
    1542             IF( crname_lbc(ji-1) /= crname_lbc(ji) ) THEN 
    1543                WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(ji-1)) 
    1544                jj = 0 
     1541         DO ji = 1, n_sequence_lbc - 1 
     1542            IF ( crname_lbc(ji) /= 'already counted' ) THEN 
     1543               ccountname = crname_lbc(ji) 
     1544               crname_lbc(ji) = 'already counted' 
     1545               jcount = 1 
     1546               DO jj = ji + 1, n_sequence_lbc 
     1547                  IF ( ccountname ==  crname_lbc(jj) ) THEN 
     1548                     jcount = jcount + 1 
     1549                     crname_lbc(jj) = 'already counted' 
     1550                  END IF 
     1551               END DO 
     1552               WRITE(numcom,'(A, I4, A, A)') ' - ', jcount,' times by subroutine ', TRIM(ccountname) 
    15451553            END IF 
    1546             jj = jj + 1  
    15471554         END DO 
    1548          WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(n_sequence_lbc)) 
     1555         IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN 
     1556            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max)) 
     1557         END IF 
    15491558         WRITE(numcom,*) ' ' 
    15501559         IF ( n_sequence_glb > 0 ) THEN 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90

    r10888 r11083  
    132132      IF( ln_zad_Aimp ) THEN 
    133133         IF( zdf_phy_alloc() /= 0 )   & 
    134         &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' ) 
    135          wi(:,:,:) = 0._wp 
     134            &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' ) 
     135         Cu_adv(:,:,:) = 0._wp 
     136         wi    (:,:,:) = 0._wp 
    136137      ENDIF 
    137138      !                          !==  Background eddy viscosity and diffusivity  ==! 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/TOP/PISCES/P4Z/p4zpoc.F90

    r10888 r11083  
    166166                       &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           & 
    167167                       &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) & 
    168                        &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn)  
     168                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn)  
    169169                       alphat = alphat + alphag(ji,jj,jk,jn) 
    170170                       remint = remint + alphag(ji,jj,jk,jn) * reminp(jn) 

  • NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90

    r10888 r11083  
    9696      !                                ! ==>>>  set by hand non-zero value on first/last columns & rows  
    9797      DO ji = mi0(1), mi1(1)              ! first row of global domain only 
    98          zhu(ji,2) = zht(1,2) 
    99       END DO 
    100        DO ji = mi0(jpi), mi1(jpi)         ! last  row of global domain only 
    101          zhu(ji,2) = zht(jpi,2) 
     98         zhu(ji,2) = zht(ji,2) 
     99      END DO 
     100       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only 
     101         zhu(ji,2) = zht(ji,2) 
    102102      END DO 
    103103      zhu(:,1) = zhu(:,2) 
Note: See TracChangeset for help on using the changeset viewer.