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

Changeset 13877


Ignore:
Timestamp:
2020-11-25T15:23:00+01:00 (3 years ago)
Author:
mocavero
Message:

Align branch with the trunk@13747

Location:
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3
Files:
70 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3

    • Property svn:externals

      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13292        sette 
         10^/utils/CI/sette@13559        sette 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/INSTALL.rst

    r11734 r13877  
    122122.. code:: console 
    123123 
    124    $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/trunk 
     124   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/r4.0/r4.0.4 
    125125 
    126126Description of 1\ :sup:`st` level tree structure 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml

    r9896 r13877  
    5252       <field field_ref="normstr"          name="normstr" /> 
    5353       <field field_ref="sheastr"          name="sheastr" /> 
    54        <field field_ref="isig1"            name="isig1"   /> 
    55        <field field_ref="isig2"            name="isig2"   /> 
    56        <field field_ref="isig3"            name="isig3"   /> 
     54       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     55       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5756        
    5857       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/C1D_PAPA/EXPREF/namelist_cfg

    r13631 r13877  
    428428/ 
    429429!----------------------------------------------------------------------- 
    430 !----------------------------------------------------------------------- 
    431 / 
    432 !----------------------------------------------------------------------- 
    433430&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    434431!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml

    r11945 r13877  
    5353       <field field_ref="normstr"          name="normstr" /> 
    5454       <field field_ref="sheastr"          name="sheastr" /> 
    55        <field field_ref="isig1"            name="isig1"   /> 
    56        <field field_ref="isig2"            name="isig2"   /> 
    57        <field field_ref="isig3"            name="isig3"   /> 
     55       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     56       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5857        
    5958       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml

    r12377 r13877  
    5353       <field field_ref="normstr"          name="normstr" /> 
    5454       <field field_ref="sheastr"          name="sheastr" /> 
    55        <field field_ref="isig1"            name="isig1"   /> 
    56        <field field_ref="isig2"            name="isig2"   /> 
    57        <field field_ref="isig3"            name="isig3"   /> 
     55       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     56       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5857        
    5958       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg

    r13631 r13877  
    388388/ 
    389389!----------------------------------------------------------------------- 
    390 !----------------------------------------------------------------------- 
    391 / 
    392 !----------------------------------------------------------------------- 
    393390&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    394391!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg

    r13631 r13877  
    384384/ 
    385385!----------------------------------------------------------------------- 
    386 !----------------------------------------------------------------------- 
    387 / 
    388 !----------------------------------------------------------------------- 
    389386&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    390387!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml

    r9896 r13877  
    5252       <field field_ref="normstr"          name="normstr" /> 
    5353       <field field_ref="sheastr"          name="sheastr" /> 
    54        <field field_ref="isig1"            name="isig1"   /> 
    55        <field field_ref="isig2"            name="isig2"   /> 
    56        <field field_ref="isig3"            name="isig3"   /> 
     54       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     55       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5756        
    5857       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/field_def_nemo-ice.xml

    r13631 r13877  
    7575          <field id="utau_bi"      long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress"              unit="N/m2" /> 
    7676          <field id="vtau_bi"      long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress"              unit="N/m2" /> 
    77           <field id="isig1"        long_name="1st principal stress component for EVP rhg"                                                                        unit=""     /> 
    78           <field id="isig2"        long_name="2nd principal stress component for EVP rhg"                                                                        unit=""     /> 
    79           <field id="isig3"        long_name="convergence measure for EVP rheology (must be around 1)"                                                           unit=""     /> 
     77          <field id="sig1_pnorm"   long_name="P-normalized 1st principal stress component"                                                                       unit=""     /> 
     78          <field id="sig2_pnorm"   long_name="P-normalized 2nd principal stress component"                                                                       unit=""     /> 
    8079          <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  /> 
    8180          <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  /> 
    … …  
    165164 
    166165          <!-- diags --> 
    167           <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" /> 
     166          <field id="hfxdhc"        long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" /> 
     167          <!-- available if ln_icediachk=T --> 
     168          <field id="icedrift_mass" long_name="Ice mass drift (conservation check)"   unit="kg/m2/s" /> 
     169          <field id="icedrift_salt" long_name="Ice salt drift (conservation check)"   unit="kg/m2/s" /> 
     170          <field id="icedrift_heat" long_name="Ice heat drift (conservation check)"   unit="W/m2"    /> 
    168171 
    169172     <!-- sbcssm variables --> 
    … …  
    400403     <field field_ref="normstr"          name="normstr" /> 
    401404     <field field_ref="sheastr"          name="sheastr" /> 
    402      <field field_ref="isig1"            name="isig1"   /> 
    403      <field field_ref="isig2"            name="isig2"   /> 
    404      <field field_ref="isig3"            name="isig3"   /> 
     405     <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     406     <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    405407      
    406408     <!-- heat fluxes --> 
    … …  
    459461          <field field_ref="vfxthin"          name="vfxthin" /> 
    460462    
    461      <!-- diag error for negative ice volume after advection --> 
    462      <field field_ref="iceneg_pres"      name="sineg_pres" /> 
    463      <field field_ref="iceneg_volu"      name="sineg_volu" /> 
    464      <field field_ref="iceneg_hfx"       name="sineg_hfx"  /> 
    465463        </field_group> 
    466464 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/namelist_ref

    r13631 r13877  
    12641264!!gm   ln_trdmld_instant = .false.         !  flag to diagnose trends of instantantaneous or mean ML T/S 
    12651265!!gm 
    1266 / 
    12671266!----------------------------------------------------------------------- 
    12681267&namhsb        !  Heat and salt budgets                                 (default: OFF) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml

    r11536 r13877  
    5353       <field field_ref="normstr"          name="normstr" /> 
    5454       <field field_ref="sheastr"          name="sheastr" /> 
    55        <field field_ref="isig1"            name="isig1"   /> 
    56        <field field_ref="isig2"            name="isig2"   /> 
    57        <field field_ref="isig3"            name="isig3"   /> 
     55       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     56       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5857        
    5958       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/file_def_nemo-ice.xml

    r12905 r13877  
    5050       <field field_ref="normstr"          name="normstr" /> 
    5151       <field field_ref="sheastr"          name="sheastr" /> 
    52        <field field_ref="isig1"            name="isig1"   /> 
    53        <field field_ref="isig2"            name="isig2"   /> 
    54        <field field_ref="isig3"            name="isig3"   /> 
     52       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     53       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5554        
    5655       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/namelist_cfg

    r13631 r13877  
    583583/ 
    584584!----------------------------------------------------------------------- 
    585 !----------------------------------------------------------------------- 
    586 / 
    587 !----------------------------------------------------------------------- 
    588585&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    589586!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/NEMO_manual_state.txt

    r10569 r13877  
    3939         namdia: iiceprt jiceprt 
    4040    nam_diaharm: nit000_han nitend_han nstep_han tname(1) tname(2) 
    41          namdrg: ln_OFF 
     41         namdrg: ln_drg_OFF 
    4242     namdrg_bot: rn_Cd0 rn_Uc0 rn_Cdmax 
    4343     namdrg_top: rn_Cd0 rn_Uc0 rn_Cdmax 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_TRA.tex

    r11693 r13877  
    12291229In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, 
    12301230with $\rho_o$ a reference density. 
    1231 Called \textit{rau0} in the code, 
     1231Called \textit{rho0} in the code, 
    12321232$\rho_o$ is set in \mdl{phycst} to a value of \texttt{1,026} $Kg/m^3$. 
    12331233This is a sensible choice for the reference density used in a Boussinesq ocean climate model, 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex

    r11693 r13877  
    11601160\] 
    11611161When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}. 
    1162 Setting \np[=.true.]{ln_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 
     1162Setting \np[=.true.]{ln_drg_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 
    11631163 
    11641164These values are assigned in \mdl{zdfdrg}. 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nambdy_dta

    r11703 r13877  
    2929   bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    3030   bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     31   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    3132   ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds 
    3233   rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice 
    … …  
    3536   rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            -- 
    3637   rn_ice_hpnd = 0.05         !       --   pond depth                          -- 
     38   rn_ice_hlid = 0.0          !       --   pond lid depth                      -- 
    3739/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdia

    r11703 r13877  
    88   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    99   ln_icectl        = .false.         !  ice points output for debug (T or F) 
    10    iiceprt          =  10             !  i-index for debug 
    11    jiceprt          =  10             !  j-index for debug 
     10      iiceprt       =  10             !     i-index for debug 
     11      jiceprt       =  10             !     j-index for debug 
    1212/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg

    r10075 r13877  
    22&namdrg        !   top/bottom drag coefficient                          (default: NO selection) 
    33!----------------------------------------------------------------------- 
    4    ln_OFF      = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
     4   ln_drg_OFF  = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
    55   ln_lin      = .false.   !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top) 
    66   ln_non_lin  = .false.   !  non-linear  drag: Cd = Cd0 |U| 
    … …  
    88   ! 
    99   ln_drgimp   = .true.    !  implicit top/bottom friction flag 
     10      ln_drgice_imp = .false. ! implicit ice-ocean drag 
    1011/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_bot

    r10075 r13877  
    11!----------------------------------------------------------------------- 
    2 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     2&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    33!----------------------------------------------------------------------- 
    44   rn_Cd0      =  1.e-3    !  drag coefficient [-] 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_top

    r10075 r13877  
    11!----------------------------------------------------------------------- 
    2 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
     2&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
    33!----------------------------------------------------------------------- 
    44   rn_Cd0      =  1.e-3    !  drag coefficient [-] 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn

    r11703 r13877  
    1010   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    1111   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
    12       rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
     12      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
    1313                                      !          recommended range: [0.1 ; 0.25] 
    14       rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
    15       rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
    16       rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??] 
     14      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
     15      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
     16      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??] 
    1717/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn_rhg

    r11025 r13877  
    99      rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    1010                                      !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
     11   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg) 
    1112/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nameos

    r10075 r13877  
    77                                 ! 
    88   !                     ! S-EOS coefficients (ln_seos=T): 
    9    !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
     9   !                             !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
    1010   rn_a0       =  1.6550e-1      !  thermal expension coefficient 
    1111   rn_b0       =  7.6554e-1      !  saline  expension coefficient 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namini

    r11703 r13877  
    33!------------------------------------------------------------------------------ 
    44   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
    5    ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F) 
     5   nn_iceini_file   =   0             !     0 = Initialise sea ice based on SSTs 
     6                                      !     1 = Initialise sea ice from single category netcdf file 
     7                                      !     2 = Initialise sea ice from multi category restart file 
    68   rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze) 
    79   rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North 
    … …  
    2325   rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North 
    2426   rn_hpd_ini_s     =   0.05          !        "            "             South 
    25    ! -- for ln_iceini_file = T 
     27   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North 
     28   rn_hld_ini_s     =   0.0           !        "            "             South 
     29   ! -- for nn_iceini_file = 1 
    2630   sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    2731   sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    … …  
    3438   sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    3539   sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     40   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    3641   cn_dir='./' 
    3742/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_blk

    r12377 r13877  
    3535   sn_tair     = 't_10.15JUNE2009_fill'       ,    6.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    3636   sn_humi     = 'q_10.15JUNE2009_fill'       ,    6.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    37    sn_hpgi     = 'NONE'                       ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , '' 
    38    sn_hpgj     = 'NONE'                       ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , '' 
    3937   sn_prec     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    4038   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    4139   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
     40   sn_uoatm    = 'NOT USED'                   ,    6.        , 'UOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', '' 
     41   sn_voatm    = 'NOT USED'                   ,    6.        , 'VOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', '' 
     42   sn_cc       = 'NOT USED'                   ,   24.        , 'CC'      ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
     43   sn_hpgi     = 'NOT USED'                   ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , '' 
     44   sn_hpgj     = 'NOT USED'                   ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , '' 
    4245/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_cpl

    r10075 r13877  
    22&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3") 
    33!----------------------------------------------------------------------- 
    4    nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data 
    5    ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models 
    6    !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
    7    nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1) 
    8  
     4   nn_cplmodel       =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data 
     5   ln_usecplmask     = .false. !  use a coupling mask file to merge data received from several models 
     6   !                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     7   ln_scale_ice_flux = .false. !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 
     8   nn_cats_cpl       =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1) 
    99   !_____________!__________________________!____________!_____________!______________________!________! 
    1010   !             !        description       !  multiple  !    vector   !       vector         ! vector ! 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd

    r11025 r13877  
    66   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F) 
    77   ln_icedS         = .true.          !  activate brine drainage (T) or not (F) 
     8   !                                    
     9   ln_leadhfx       = .true.          !  heat in the leads is used to melt sea-ice before warming the ocean 
    810/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_pnd

    r11536 r13877  
    22&namthd_pnd     !   Melt ponds 
    33!------------------------------------------------------------------------------ 
    4    ln_pnd           = .false.         !  activate melt ponds or not 
    5      ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
    6      ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    7        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
    8        rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC 
    9      ln_pnd_alb     = .false.         !  melt ponds affect albedo or not 
     4   ln_pnd            = .false.        !  activate melt ponds or not 
     5      ln_pnd_LEV     = .false.        !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 
     6         rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 
     7         rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 
     8      ln_pnd_CST     = .false.        !  constant  melt ponds 
     9         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC 
     10         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC 
     11      ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV) 
     12      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo 
    1013/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_zdf

    r11025 r13877  
    77   rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 
    88                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    9    rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m] 
     9   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice                     [1/m] 
     10   rn_kappa_s       =  10.0           !  nn_qtrice = 0: radiation attenuation coefficient in snow         [1/m] 
     11   rn_kappa_smlt    =   7.0           !  nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m] 
     12   rn_kappa_sdry    =  10.0           !                 radiation attenuation coefficient in dry snow     [1/m] 
     13   ln_zdf_chkcvg    = .false.         !  check convergence of heat diffusion scheme (output variable: tice_cvg) 
    1014/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_gls

    r9355 r13877  
    1313   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
    1414   !                             ! =3 requires ln_wave=T 
     15   nn_z0_ice     =   1     !  attenutaion of surface wave breaking under ice 
     16   !                       !           = 0 no impact of ice cover 
     17   !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i) 
     18   !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i 
     19   !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i) 
    1520   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    1621   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_tke

    r10075 r13877  
    1515   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F) 
    1616   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value 
    17    ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE 
    1817   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002) 
    1918      rn_lc       =   0.15    !  coef. associated to Langmuir cells 
    … …  
    2625                              !        = 0  constant 10 m length scale 
    2726                              !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
    28       rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4    
     27   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice 
     28   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking 
     29   !                       !           = 1 weigthed by 1-TANH(10*fr_i) 
     30   !                       !           = 2 weighted by 1-fr_i 
     31   !                       !           = 3 weighted by 1-MIN(1,4*fr_i) 
    2932/ 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice.F90

    r13571 r13877  
    391391   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice         !: ice volume variation   [m/s]  
    392392   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
     393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1]  
     394   ! 
     395   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s) 
     396   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_salt     !: advection of salt (g/m2/s) 
     397   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_heat     !: advection of heat (W/m2) 
    393398   ! 
    394399   !!---------------------------------------------------------------------- 
    … …  
    493498      ! * Ice diagnostics 
    494499      ii = ii + 1 
    495       ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),   &  
    496          &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),   & 
    497          &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) ) 
     500      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      &  
     501         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      & 
     502         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj),  & 
     503         &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 
    498504 
    499505      ! * Ice conservation 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice1d.F90

    r13571 r13877  
    145145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sst_1d 
    146146   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frq_m_1d 
    147148 
    148149   ! convergence check 
    … …  
    225226      ! 
    226227      ii = ii + 1 
    227       ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) ) 
     228      ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , frq_m_1d(jpij) , STAT=ierr(ii) ) 
    228229      ! 
    229230      ii = ii + 1 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icecor.F90

    r13630 r13877  
    5555      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    5656      REAL(wp) ::   zsal, zzc 
    57       REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag 
    5857      !!---------------------------------------------------------------------- 
    5958      ! controls 
    … …  
    9594               zsal = sv_i(ji,jj,jl) 
    9695               sv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  ) 
    97                sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux 
     96               IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 
     97                  &   sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux 
    9898            END_2D 
    9999         END DO 
    100100      ENDIF 
    101       !                             !----------------------------------------------------- 
    102       CALL ice_var_zapsmall         !  Zap small values                                  ! 
    103       !                             !----------------------------------------------------- 
    104101 
     102      IF( kn /= 0 ) THEN   ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 
     103         !                                                              otherwise conservation diags will fail 
     104         !                          !----------------------------------------------------- 
     105         CALL ice_var_zapsmall      !  Zap small values                                  ! 
     106         !                          !----------------------------------------------------- 
     107      ENDIF 
    105108      !                             !----------------------------------------------------- 
    106109      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
    … …  
    119122#endif 
    120123      ENDIF 
    121  
    122       !                             !----------------------------------------------------- 
    123       SELECT CASE( kn )             !  Diagnostics                                       ! 
    124       !                             !----------------------------------------------------- 
    125       CASE( 1 )                        !--- dyn trend diagnostics 
    126          ! 
    127          IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
    128             diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &      ! W.m-2 
    129                &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
    130             diag_sice(:,:) =   SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    131             diag_vice(:,:) =   SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    132             diag_vsnw(:,:) =   SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
    133          ENDIF 
    134          !                       ! concentration tendency (dynamics) 
    135          IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
    136             zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice  
    137             CALL iom_put( 'afxdyn' , zafx ) 
    138          ENDIF 
    139          ! 
    140       CASE( 2 )                        !--- thermo trend diagnostics & ice aging 
    141          ! 
    142          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice   ! ice natural aging incrementation 
    143          ! 
    144          IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
    145             diag_heat(:,:) = diag_heat(:,:) & 
    146                &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 
    147                &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
    148             diag_sice(:,:) = diag_sice(:,:) & 
    149                &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    150             diag_vice(:,:) = diag_vice(:,:) & 
    151                &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    152             diag_vsnw(:,:) = diag_vsnw(:,:) & 
    153                &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
    154             CALL iom_put ( 'hfxdhc' , diag_heat )  
    155          ENDIF 
    156          !                       ! concentration tendency (total + thermo) 
    157          IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
    158             zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
    159             CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) 
    160             CALL iom_put( 'afxtot' , zafx ) 
    161          ENDIF 
    162          ! 
    163       END SELECT 
    164124      ! 
    165125      ! controls 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icectl.F90

    r13571 r13877  
    4343   PUBLIC   ice_prt 
    4444   PUBLIC   ice_prt3D 
     45   PUBLIC   ice_drift_wri 
     46   PUBLIC   ice_drift_init 
    4547 
    4648   ! thresold rates for conservation 
    … …  
    4951   REAL(wp), PARAMETER ::   zchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 
    5052   REAL(wp), PARAMETER ::   zchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 
     53 
     54   ! for drift outputs 
     55   CHARACTER(LEN=50)   ::   clname="icedrift_diagnostics.ascii"   ! ascii filename 
     56   INTEGER             ::   numicedrift                           ! outfile unit 
     57   REAL(wp)            ::   rdiag_icemass, rdiag_icesalt, rdiag_iceheat  
     58   REAL(wp)            ::   rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat  
    5159    
    5260   !! * Substitutions 
    … …  
    132140 
    133141         ! -- advection scheme is conservative? -- ! 
    134          zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather) 
    135          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 (only for Prather) 
     142         zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 
     143         zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
    136144 
    137145         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    156164               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
    157165            ! check if advection scheme is conservative 
    158             !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 
    159             !    so the formulation for conservation is different (and not coded)  
    160             !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 
    161             !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    162             !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 
     166            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     167               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
     168            IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     169               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rdt_ice 
    163170         ENDIF 
    164171         ! 
    … …  
    186193      ! water flux 
    187194      ! -- mass diag -- ! 
    188       zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 
     195      zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
     196         &                              + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 
    189197 
    190198      ! -- salt diag -- ! 
    191       zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 
     199      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 
    192200 
    193201      ! -- heat diag -- ! 
    194       ! clem: not the good formulation 
    195       !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
    196       !!   &                              ) * e1e2t ) 
     202      zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     203      ! equivalent to this: 
     204      !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
     205      !!   &                                          - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 
     206      !!   &                                          ) * e1e2t ) 
    197207 
    198208      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    204214         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
    205215            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice 
    206          !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
     216         IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
     217            &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
    207218      ENDIF 
    208219      ! 
    … …  
    725736       
    726737   END SUBROUTINE ice_prt3D 
     738 
     739 
     740   SUBROUTINE ice_drift_wri( kt ) 
     741      !!------------------------------------------------------------------- 
     742      !!                     ***  ROUTINE ice_drift_wri *** 
     743      !! 
     744      !! ** Purpose : conservation of mass, salt and heat 
     745      !!              write the drift in a ascii file at each time step 
     746      !!              and the total run drifts 
     747      !!------------------------------------------------------------------- 
     748      INTEGER, INTENT(in) ::   kt   ! ice time-step index 
     749      ! 
     750      INTEGER  ::   ji, jj 
     751      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 
     752      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 
     753      !!------------------------------------------------------------------- 
     754      ! 
     755      IF( kt == nit000 .AND. lwp ) THEN 
     756         WRITE(numout,*) 
     757         WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 
     758         WRITE(numout,*) '~~~~~~~~~~~~~' 
     759      ENDIF 
     760      ! 
     761      ! 2D budgets (must be close to 0) 
     762      IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 
     763         DO_2D( 1, 1, 1, 1 ) 
     764            zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 
     765               &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 
     766            zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 
     767            zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 
     768         END_2D 
     769         ! 
     770         ! write outputs 
     771         CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 
     772         CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 
     773         CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 
     774      ENDIF 
     775 
     776      ! -- mass diag -- ! 
     777      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
     778         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 
     779      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 
     780 
     781      ! -- salt diag -- ! 
     782      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 
     783      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 
     784 
     785      ! -- heat diag -- ! 
     786      zdiag_heat     = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     787      zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
     788 
     789      !                    ! write out to file 
     790      IF( lwp ) THEN 
     791         ! check global drift (must be close to 0) 
     792         WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zdiag_mass 
     793         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zdiag_salt 
     794         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zdiag_heat 
     795         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     796         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zdiag_adv_mass 
     797         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zdiag_adv_salt 
     798         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zdiag_adv_heat 
     799      ENDIF 
     800      !                    ! drifts 
     801      rdiag_icemass = rdiag_icemass + zdiag_mass 
     802      rdiag_icesalt = rdiag_icesalt + zdiag_salt 
     803      rdiag_iceheat = rdiag_iceheat + zdiag_heat 
     804      rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 
     805      rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 
     806      rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 
     807      ! 
     808      !                    ! output drifts and close ascii file 
     809      IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 
     810         ! to ascii file 
     811         WRITE(numicedrift,*) '******************************************' 
     812         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift     [kg]', rdiag_icemass 
     813         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 
     814         WRITE(numicedrift,*) '******************************************' 
     815         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift     [kg]', rdiag_icesalt 
     816         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 
     817         WRITE(numicedrift,*) '******************************************' 
     818         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift     [W] ', rdiag_iceheat 
     819         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 
     820         CLOSE( numicedrift ) 
     821         ! 
     822         ! to ocean output 
     823         WRITE(numout,*) 
     824         WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 
     825         WRITE(numout,*) '~~~~~~~~~~~~~' 
     826         ! check global drift (must be close to 0) 
     827         WRITE(numout,*) '   sea-ice mass drift     [kg] = ', rdiag_icemass 
     828         WRITE(numout,*) '   sea-ice salt drift     [kg] = ', rdiag_icesalt 
     829         WRITE(numout,*) '   sea-ice heat drift     [W]  = ', rdiag_iceheat 
     830         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     831         WRITE(numout,*) '   sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 
     832         WRITE(numout,*) '   sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 
     833         WRITE(numout,*) '   sea-ice heat drift adv [W]  = ', rdiag_adv_iceheat 
     834      ENDIF 
     835      ! 
     836   END SUBROUTINE ice_drift_wri 
     837 
     838   SUBROUTINE ice_drift_init 
     839      !!---------------------------------------------------------------------- 
     840      !!                  ***  ROUTINE ice_drift_init  *** 
     841      !!                    
     842      !! ** Purpose :   create output file, initialise arrays 
     843      !!---------------------------------------------------------------------- 
     844      ! 
     845      IF( .NOT.ln_icediachk ) RETURN ! exit 
     846      ! 
     847      IF(lwp) THEN 
     848         WRITE(numout,*) 
     849         WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 
     850         WRITE(numout,*) '~~~~~~~~~~~~~' 
     851         WRITE(numout,*) 
     852         ! 
     853         ! create output ascii file 
     854         CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 
     855         WRITE(numicedrift,*) 'Timestep  Drifts' 
     856         WRITE(numicedrift,*) '******************************************' 
     857      ENDIF 
     858      ! 
     859      rdiag_icemass = 0._wp 
     860      rdiag_icesalt = 0._wp 
     861      rdiag_iceheat = 0._wp 
     862      rdiag_adv_icemass = 0._wp 
     863      rdiag_adv_icesalt = 0._wp 
     864      rdiag_adv_iceheat = 0._wp 
     865      ! 
     866   END SUBROUTINE ice_drift_init 
    727867       
    728868#else 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90

    r13807 r13877  
    8888      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8989      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    90       REAL(wp) ::   zdt                     !   -      - 
     90      REAL(wp) ::   zdt, z1_dt              !   -      - 
    9191      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9292      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    … …  
    100100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    101101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     102      !! diagnostics 
     103      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    102104      !!---------------------------------------------------------------------- 
    103105      ! 
    … …  
    109111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    110112      END WHERE 
    111       DO jl = 1, jpl 
    112          DO_2D( 0, 0, 0, 0 ) 
    113             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    114                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    115                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    116                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    117             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    118                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    119                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    120                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    121             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    122                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    123                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    124                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    125             zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    126                &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    127                &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    128                &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    129          END_2D 
    130       END DO 
     113      CALL icemax3D( ph_i , zhi_max ) 
     114      CALL icemax3D( ph_s , zhs_max ) 
     115      CALL icemax3D( ph_ip, zhip_max) 
     116      CALL icemax3D( zs_i , zsi_max ) 
    131117#if defined key_mpi3 
    132118      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    … …  
    145131         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    146132         END WHERE 
    147       END DO 
    148       DO jl = 1, jpl 
    149          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    150             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
    151                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    152                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    153                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    154          END_3D 
    155       END DO 
    156       DO jl = 1, jpl 
    157          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    158             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
    159                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    160                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    161                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    162          END_3D 
    163       END DO 
     133      END DO    
     134      CALL icemax4D( ze_i , zei_max ) 
     135      CALL icemax4D( ze_s , zes_max ) 
    164136#if defined key_mpi3 
    165       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    166       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     137      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     138      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
    167139#else 
    168       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    169       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     140      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     141      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
    170142#endif 
    171143      ! 
    … …  
    184156      ENDIF 
    185157      zdt = rDt_ice / REAL(icycle) 
     158      z1_dt = 1._wp / zdt 
    186159       
    187160      ! --- transport --- ! 
    … …  
    190163 
    191164      DO jt = 1, icycle 
     165 
     166         ! diagnostics 
     167         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     168         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     169         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     170            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    192171 
    193172         ! record at_i before advection (for open water) 
    … …  
    288267                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
    289268               ENDIF 
    290            ENDIF 
    291             ! 
     269            ENDIF 
     270            ! 
     271         ENDIF 
     272          
     273         ! --- Lateral boundary conditions --- ! 
     274         !     caution: for gradients (sx and sy) the sign changes 
     275         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume 
     276            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  & 
     277            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume 
     278            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  ) 
     279         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity 
     280            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  & 
     281            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration 
     282            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  ) 
     283         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age 
     284            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  ) 
     285         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy 
     286            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  )  
     287         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
     288            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
     289         IF ( ln_pnd_LEV ) THEN 
     290            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
     291               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
     292               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume 
     293               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )  
     294            IF ( ln_pnd_lids ) THEN 
     295               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume 
     296                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )  
     297            ENDIF 
    292298         ENDIF 
    293299 
    … …  
    325331         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp ) 
    326332#endif 
     333         ! 
     334         ! --- diagnostics --- ! 
     335         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     336            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     337         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     338            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     339         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     340            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     341            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    327342         ! 
    328343         ! --- Ensure non-negative fields --- ! 
    … …  
    363378      !!  
    364379      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     380      INTEGER  ::   jj0                                  ! dummy loop indices 
    365381      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    366382      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    367383      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     384      REAL(wp) ::   zpsm, zps0 
     385      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    368386      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    369387      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    370388      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    371389      !----------------------------------------------------------------------- 
     390      ! in order to avoid lbc_lnk (communications): 
     391      !    jj loop must be 1:jpj   if adv_x is called first 
     392      !                and 2:jpj-1 if adv_x is called second 
     393      jj0 = NINT(pcrh) 
    372394      ! 
    373395      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    376398         ! 
    377399         ! Limitation of moments.                                            
    378          DO_2D( 0, 0, 1, 1 ) 
    379             !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    380             psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    381             ! 
    382             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    383             zs1max  = 1.5 * zslpmax 
    384             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    385             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    386                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    387             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    388  
    389             ps0 (ji,jj,jl) = zslpmax   
    390             psx (ji,jj,jl) = zs1new         * rswitch 
    391             psxx(ji,jj,jl) = zs2new         * rswitch 
    392             psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    393             psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    394             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    395          END_2D 
    396  
    397          !  Calculate fluxes and moments between boxes i<-->i+1               
    398          DO_2D( 0, 0, 1, 1 )                   !  Flux from i to i+1 WHEN u GT 0 
    399             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    400             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    401             zalfq        =  zalf * zalf 
    402             zalf1        =  1.0 - zalf 
    403             zalf1q       =  zalf1 * zalf1 
    404             ! 
    405             zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    406             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    407             zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    408             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    409             zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    410             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    411             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    412  
    413             !  Readjust moments remaining in the box. 
    414             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    415             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    416             psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    417             psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    418             psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    419             psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    420             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    421          END_2D 
    422  
    423          DO_2D( 0, 0, 1, 0 )                   !  Flux from i+1 to i when u LT 0. 
    424             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    425             zalg  (ji,jj) = zalf 
    426             zalfq         = zalf * zalf 
    427             zalf1         = 1.0 - zalf 
    428             zalg1 (ji,jj) = zalf1 
    429             zalf1q        = zalf1 * zalf1 
    430             zalg1q(ji,jj) = zalf1q 
    431             ! 
    432             zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
    433             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
    434                &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
    435             zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
    436             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
    437             zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
    438             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
    439             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    440          END_2D 
    441  
    442          DO_2D( 0, 0, 0, 0 )                   !  Readjust moments remaining in the box. 
    443             zbt  =       zbet(ji-1,jj) 
    444             zbt1 = 1.0 - zbet(ji-1,jj) 
    445             ! 
    446             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    447             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    448             psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    449             psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    450             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    451             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    452             psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    453          END_2D 
    454  
    455          !   Put the temporary moments into appropriate neighboring boxes.     
    456          DO_2D( 0, 0, 0, 0 )                   !   Flux from i to i+1 IF u GT 0. 
    457             zbt  =       zbet(ji-1,jj) 
    458             zbt1 = 1.0 - zbet(ji-1,jj) 
    459             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    460             zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    461             zalf1         = 1.0 - zalf 
    462             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    463             ! 
    464             ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    465             psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    466             psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    467                &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
    468                &            + zbt1 * psxx(ji,jj,jl) 
    469             psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    470                &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    471                &            + zbt1 * psxy(ji,jj,jl) 
    472             psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    473             psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    474          END_2D 
    475  
    476          DO_2D( 0, 0, 0, 0 )                   !  Flux from i+1 to i IF u LT 0. 
    477             zbt  =       zbet(ji,jj) 
    478             zbt1 = 1.0 - zbet(ji,jj) 
    479             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    480             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    481             zalf1         = 1.0 - zalf 
    482             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    483             ! 
    484             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    485             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    486             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    487                &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    488                &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    489             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    490                &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    491             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    492             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
    493          END_2D 
    494  
     400         DO jj = Njs0 - jj0, Nje0 + jj0 
     401             
     402            DO ji = Nis0 - 1, Nie0 + 1 
     403 
     404               zpsm  = psm (ji,jj,jl) ! optimization 
     405               zps0  = ps0 (ji,jj,jl) 
     406               zpsx  = psx (ji,jj,jl) 
     407               zpsxx = psxx(ji,jj,jl) 
     408               zpsy  = psy (ji,jj,jl) 
     409               zpsyy = psyy(ji,jj,jl) 
     410               zpsxy = psxy(ji,jj,jl) 
     411 
     412               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     413               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     414               ! 
     415               zslpmax = MAX( 0._wp, zps0 ) 
     416               zs1max  = 1.5 * zslpmax 
     417               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     418               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
     419               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     420 
     421               zps0  = zslpmax   
     422               zpsx  = zs1new  * rswitch 
     423               zpsxx = zs2new  * rswitch 
     424               zpsy  = zpsy    * rswitch 
     425               zpsyy = zpsyy   * rswitch 
     426               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     427 
     428               !  Calculate fluxes and moments between boxes i<-->i+1               
     429               !                                !  Flux from i to i+1 WHEN u GT 0  
     430               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     431               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
     432               zalfq        =  zalf * zalf 
     433               zalf1        =  1.0 - zalf 
     434               zalf1q       =  zalf1 * zalf1 
     435               ! 
     436               zfm (ji,jj)  =  zalf  *   zpsm  
     437               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     438               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     439               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     440               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     441               zfxy(ji,jj)  =  zalfq *   zpsxy 
     442               zfyy(ji,jj)  =  zalf  *   zpsyy 
     443 
     444               !                                !  Readjust moments remaining in the box. 
     445               zpsm  =  zpsm  - zfm(ji,jj) 
     446               zps0  =  zps0  - zf0(ji,jj) 
     447               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     448               zpsxx =  zalf1  * zalf1q * zpsxx 
     449               zpsy  =  zpsy  - zfy (ji,jj) 
     450               zpsyy =  zpsyy - zfyy(ji,jj) 
     451               zpsxy =  zalf1q * zpsxy 
     452               ! 
     453               psm (ji,jj,jl) = zpsm ! optimization 
     454               ps0 (ji,jj,jl) = zps0  
     455               psx (ji,jj,jl) = zpsx  
     456               psxx(ji,jj,jl) = zpsxx 
     457               psy (ji,jj,jl) = zpsy  
     458               psyy(ji,jj,jl) = zpsyy 
     459               psxy(ji,jj,jl) = zpsxy 
     460               ! 
     461            END DO 
     462             
     463            DO ji = Nis0 - 1, Nie0 
     464               !                                !  Flux from i+1 to i when u LT 0. 
     465               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     466               zalg  (ji,jj) = zalf 
     467               zalfq         = zalf * zalf 
     468               zalf1         = 1.0 - zalf 
     469               zalg1 (ji,jj) = zalf1 
     470               zalf1q        = zalf1 * zalf1 
     471               zalg1q(ji,jj) = zalf1q 
     472               ! 
     473               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     474               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     475                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     476               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     477               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     478               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     479               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     480               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     481            END DO 
     482 
     483            DO ji = Nis0, Nie0 
     484               ! 
     485               zpsm  = psm (ji,jj,jl) ! optimization 
     486               zps0  = ps0 (ji,jj,jl) 
     487               zpsx  = psx (ji,jj,jl) 
     488               zpsxx = psxx(ji,jj,jl) 
     489               zpsy  = psy (ji,jj,jl) 
     490               zpsyy = psyy(ji,jj,jl) 
     491               zpsxy = psxy(ji,jj,jl) 
     492               !                                !  Readjust moments remaining in the box. 
     493               zbt  =       zbet(ji-1,jj) 
     494               zbt1 = 1.0 - zbet(ji-1,jj) 
     495               ! 
     496               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     497               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     498               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     499               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     500               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) ) 
     501               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     502               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     503 
     504               !   Put the temporary moments into appropriate neighboring boxes.     
     505               !                                !   Flux from i to i+1 IF u GT 0. 
     506               zbt   =       zbet(ji-1,jj) 
     507               zbt1  = 1.0 - zbet(ji-1,jj) 
     508               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     509               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     510               zalf1 = 1.0 - zalf 
     511               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     512               ! 
     513               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     514               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     515               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     516                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     517                  &            + zbt1 * zpsxx 
     518               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     519                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     520                  &            + zbt1 * zpsxy 
     521               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     522               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     523 
     524               !                                !  Flux from i+1 to i IF u LT 0. 
     525               zbt   =       zbet(ji,jj) 
     526               zbt1  = 1.0 - zbet(ji,jj) 
     527               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     528               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     529               zalf1 = 1.0 - zalf 
     530               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     531               ! 
     532               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     533               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     534               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     535                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     536                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     537               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     538                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     539               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     540               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     541               ! 
     542               psm (ji,jj,jl) = zpsm  ! optimization 
     543               ps0 (ji,jj,jl) = zps0  
     544               psx (ji,jj,jl) = zpsx  
     545               psxx(ji,jj,jl) = zpsxx 
     546               psy (ji,jj,jl) = zpsy  
     547               psyy(ji,jj,jl) = zpsyy 
     548               psxy(ji,jj,jl) = zpsxy 
     549            END DO 
     550            ! 
     551         END DO 
     552         ! 
    495553      END DO 
    496  
    497       !-- Lateral boundary conditions 
    498 #if defined key_mpi3 
    499       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    500          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    501          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    502 #else 
    503       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    504          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    505          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    506 #endif 
    507       ! 
     554      !       
    508555   END SUBROUTINE adv_x 
    509556 
    … …  
    526573      !! 
    527574      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     575      INTEGER  ::   ji0                                  ! dummy loop indices 
    528576      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    529577      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    530578      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     579      REAL(wp) ::   zpsm, zps0 
     580      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    531581      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    532582      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    533583      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    534584      !--------------------------------------------------------------------- 
     585      ! in order to avoid lbc_lnk (communications): 
     586      !    ji loop must be 1:jpi   if adv_y is called first 
     587      !                and 2:jpi-1 if adv_y is called second 
     588      ji0 = NINT(pcrh) 
    535589      ! 
    536590      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    539593         ! 
    540594         ! Limitation of moments. 
    541          DO_2D( 1, 1, 0, 0 ) 
    542             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    543             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    544             ! 
    545             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     595         DO_2D( 1, 1, ji0, ji0 ) 
     596            ! 
     597            zpsm  = psm (ji,jj,jl) ! optimization 
     598            zps0  = ps0 (ji,jj,jl) 
     599            zpsx  = psx (ji,jj,jl) 
     600            zpsxx = psxx(ji,jj,jl) 
     601            zpsy  = psy (ji,jj,jl) 
     602            zpsyy = psyy(ji,jj,jl) 
     603            zpsxy = psxy(ji,jj,jl) 
     604            ! 
     605            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     606            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     607            ! 
     608            zslpmax = MAX( 0._wp, zps0 ) 
    546609            zs1max  = 1.5 * zslpmax 
    547             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    548             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    549                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     610            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     611            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    550612            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    551613            ! 
    552             ps0 (ji,jj,jl) = zslpmax   
    553             psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    554             psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    555             psy (ji,jj,jl) = zs1new         * rswitch 
    556             psyy(ji,jj,jl) = zs2new         * rswitch 
    557             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    558          END_2D 
    559   
    560          !  Calculate fluxes and moments between boxes j<-->j+1               
    561          DO_2D( 1, 1, 0, 0 )                   !  Flux from j to j+1 WHEN v GT 0 
     614            zps0  = zslpmax   
     615            zpsx  = zpsx  * rswitch 
     616            zpsxx = zpsxx * rswitch 
     617            zpsy  = zs1new         * rswitch 
     618            zpsyy = zs2new         * rswitch 
     619            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     620 
     621            !  Calculate fluxes and moments between boxes j<-->j+1               
     622            !                                !  Flux from j to j+1 WHEN v GT 0    
    562623            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    563             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     624            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    564625            zalfq        =  zalf * zalf 
    565626            zalf1        =  1.0 - zalf 
    566627            zalf1q       =  zalf1 * zalf1 
    567628            ! 
    568             zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    569             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    570             zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    571             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    572             zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    573             zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    574             zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
    575             ! 
    576             !  Readjust moments remaining in the box. 
    577             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    578             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    579             psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    580             psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    581             psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    582             psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    583             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     629            zfm (ji,jj)  =  zalf  * zpsm 
     630            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     631            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     632            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     633            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     634            zfxy(ji,jj)  =  zalfq * zpsxy 
     635            zfxx(ji,jj)  =  zalf  * zpsxx 
     636            ! 
     637            !                                !  Readjust moments remaining in the box. 
     638            zpsm   =  zpsm  - zfm(ji,jj) 
     639            zps0   =  zps0  - zf0(ji,jj) 
     640            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     641            zpsyy  =  zalf1 * zalf1q * zpsyy 
     642            zpsx   =  zpsx  - zfx(ji,jj) 
     643            zpsxx  =  zpsxx - zfxx(ji,jj) 
     644            zpsxy  =  zalf1q * zpsxy 
     645            ! 
     646            psm (ji,jj,jl) = zpsm ! optimization 
     647            ps0 (ji,jj,jl) = zps0  
     648            psx (ji,jj,jl) = zpsx  
     649            psxx(ji,jj,jl) = zpsxx 
     650            psy (ji,jj,jl) = zpsy  
     651            psyy(ji,jj,jl) = zpsyy 
     652            psxy(ji,jj,jl) = zpsxy 
    584653         END_2D 
    585654         ! 
    586          DO_2D( 1, 0, 0, 0 )                   !  Flux from j+1 to j when v LT 0. 
     655         DO_2D( 1, 0, ji0, ji0 ) 
     656            !                                !  Flux from j+1 to j when v LT 0. 
    587657            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    588658            zalg  (ji,jj) = zalf 
    … …  
    603673         END_2D 
    604674 
    605          !  Readjust moments remaining in the box.  
    606          DO_2D( 0, 0, 0, 0 ) 
     675         DO_2D( 0, 0, ji0, ji0 ) 
     676            !                                !  Readjust moments remaining in the box. 
    607677            zbt  =         zbet(ji,jj-1) 
    608678            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    609679            ! 
    610             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    611             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    612             psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    613             psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    614             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    615             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    616             psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     680            zpsm  = psm (ji,jj,jl) ! optimization 
     681            zps0  = ps0 (ji,jj,jl) 
     682            zpsx  = psx (ji,jj,jl) 
     683            zpsxx = psxx(ji,jj,jl) 
     684            zpsy  = psy (ji,jj,jl) 
     685            zpsyy = psyy(ji,jj,jl) 
     686            zpsxy = psxy(ji,jj,jl) 
     687            ! 
     688            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     689            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     690            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     691            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     692            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     693            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     694            zpsxy = zalg1q(ji,jj-1) * zpsxy 
     695 
     696            !   Put the temporary moments into appropriate neighboring boxes.     
     697            !                                !   Flux from j to j+1 IF v GT 0. 
     698            zbt   =       zbet(ji,jj-1) 
     699            zbt1  = 1.0 - zbet(ji,jj-1) 
     700            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     701            zalf  = zbt * zfm(ji,jj-1) / zpsm  
     702            zalf1 = 1.0 - zalf 
     703            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     704            ! 
     705            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     706            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     707               &             + zbt1 * zpsy   
     708            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     709               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     710               &             + zbt1 * zpsyy 
     711            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     712               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     713               &             + zbt1 * zpsxy 
     714            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     715            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
     716 
     717            !                                !  Flux from j+1 to j IF v LT 0. 
     718            zbt   =       zbet(ji,jj) 
     719            zbt1  = 1.0 - zbet(ji,jj) 
     720            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     721            zalf  = zbt1 * zfm(ji,jj) / zpsm 
     722            zalf1 = 1.0 - zalf 
     723            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     724            ! 
     725            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     726            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     727            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     728               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     729               &                         + ( zalf1 - zalf ) * ztemp ) ) 
     730            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     731               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     732            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     733            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     734            ! 
     735            psm (ji,jj,jl) = zpsm ! optimization 
     736            ps0 (ji,jj,jl) = zps0  
     737            psx (ji,jj,jl) = zpsx  
     738            psxx(ji,jj,jl) = zpsxx 
     739            psy (ji,jj,jl) = zpsy  
     740            psyy(ji,jj,jl) = zpsyy 
     741            psxy(ji,jj,jl) = zpsxy 
    617742         END_2D 
    618  
    619          !   Put the temporary moments into appropriate neighboring boxes.     
    620          DO_2D( 0, 0, 0, 0 )                   !  Flux from j to j+1 IF v GT 0. 
    621             zbt  =       zbet(ji,jj-1) 
    622             zbt1 = 1.0 - zbet(ji,jj-1) 
    623             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    624             zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    625             zalf1         = 1.0 - zalf 
    626             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    627             ! 
    628             ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    629             psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    630                &             + zbt1 * psy(ji,jj,jl)   
    631             psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    632                &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    633                &             + zbt1 * psyy(ji,jj,jl) 
    634             psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    635                &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    636                &             + zbt1 * psxy(ji,jj,jl) 
    637             psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    638             psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
    639          END_2D 
    640  
    641          DO_2D( 0, 0, 0, 0 )                   !  Flux from j+1 to j IF v LT 0. 
    642             zbt  =       zbet(ji,jj) 
    643             zbt1 = 1.0 - zbet(ji,jj) 
    644             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    645             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    646             zalf1         = 1.0 - zalf 
    647             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    648             ! 
    649             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    650             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    651             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    652                &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    653                &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    654             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    655                &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    656             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    657             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
    658          END_2D 
    659  
     743         ! 
    660744      END DO 
    661  
    662       !-- Lateral boundary conditions 
    663 #if defined key_mpi3 
    664       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    665          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    666          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    667 #else 
    668       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    669          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    670          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    671 #endif 
    672745      ! 
    673746   END SUBROUTINE adv_y 
    … …  
    897970            ! 
    898971            !                                                        ! ice thickness 
    899             CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice ) 
    900             CALL iom_get( numrir, jpdom_auto, 'syice' , syice ) 
     972            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 
     973            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 
    901974            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 
    902975            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 
    903976            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 
    904977            !                                                        ! snow thickness 
    905             CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  ) 
    906             CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  ) 
     978            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp ) 
     979            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp ) 
    907980            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  ) 
    908981            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  ) 
    909982            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  ) 
    910983            !                                                        ! ice concentration 
    911             CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    ) 
    912             CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    ) 
     984            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp ) 
     985            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp ) 
    913986            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   ) 
    914987            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   ) 
    915988            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   ) 
    916989            !                                                        ! ice salinity 
    917             CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal ) 
    918             CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal ) 
     990            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 
     991            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 
    919992            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 
    920993            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 
    921994            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 
    922995            !                                                        ! ice age 
    923             CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage ) 
    924             CALL iom_get( numrir, jpdom_auto, 'syage' , syage ) 
     996            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 
     997            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 
    925998            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 
    926999            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 
    … …  
    9291002            DO jk = 1, nlay_s 
    9301003               WRITE(zchar1,'(I2.2)') jk 
    931                znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
    932                znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
     1004               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
     1005               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
    9331006               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
    9341007               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
    … …  
    9381011            DO jk = 1, nlay_i 
    9391012               WRITE(zchar1,'(I2.2)') jk 
    940                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
    941                znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
     1013               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
     1014               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
    9421015               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
    9431016               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
    … …  
    9461019            ! 
    9471020            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
    948                IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
    949                   CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap ) 
    950                   CALL iom_get( numrir, jpdom_auto, 'syap' , syap ) 
     1021               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     1022                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 
     1023                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 
    9511024                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
    9521025                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
    9531026                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
    9541027                  !                                                     ! melt pond volume 
    955                   CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp ) 
    956                   CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp ) 
     1028                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 
     1029                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 
    9571030                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
    9581031                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
    … …  
    9641037                  ! 
    9651038               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
    966                   IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
    967                      CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl ) 
    968                      CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl ) 
     1039                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     1040                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 
     1041                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 
    9691042                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 
    9701043                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 
    … …  
    10811154   END SUBROUTINE adv_pra_rst 
    10821155 
     1156   SUBROUTINE icemax3D( pice , pmax ) 
     1157      !!--------------------------------------------------------------------- 
     1158      !!                   ***  ROUTINE icemax3D ***                      
     1159      !! ** Purpose :  compute the max of the 9 points around 
     1160      !!---------------------------------------------------------------------- 
     1161      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1162      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1163      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1164      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1165      !!---------------------------------------------------------------------- 
     1166      DO jl = 1, jpl 
     1167         DO jj = Njs0-1, Nje0+1     
     1168            DO ji = Nis0, Nie0 
     1169               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1170            END DO 
     1171         END DO 
     1172         DO jj = Njs0, Nje0     
     1173            DO ji = Nis0, Nie0 
     1174               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1175            END DO 
     1176         END DO 
     1177      END DO 
     1178   END SUBROUTINE icemax3D 
     1179 
     1180   SUBROUTINE icemax4D( pice , pmax ) 
     1181      !!--------------------------------------------------------------------- 
     1182      !!                   ***  ROUTINE icemax4D ***                      
     1183      !! ** Purpose :  compute the max of the 9 points around 
     1184      !!---------------------------------------------------------------------- 
     1185      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1186      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1187      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1188      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1189      !!---------------------------------------------------------------------- 
     1190      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1191      DO jl = 1, jpl 
     1192         DO jk = 1, jlay 
     1193            DO jj = Njs0-1, Nje0+1     
     1194               DO ji = Nis0, Nie0 
     1195                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1196               END DO 
     1197            END DO 
     1198            DO jj = Njs0, Nje0     
     1199               DO ji = Nis0, Nie0 
     1200                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1201               END DO 
     1202            END DO 
     1203         END DO 
     1204      END DO 
     1205   END SUBROUTINE icemax4D 
     1206 
    10831207#else 
    10841208   !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90

    r13807 r13877  
    9292      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    9393      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    94       REAL(wp) ::   zdt, zvi_cen 
     94      REAL(wp) ::   zdt, z1_dt, zvi_cen 
    9595      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9696      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box 
    … …  
    104104      ! 
    105105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
     106      !! diagnostics 
     107      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    106108      !!---------------------------------------------------------------------- 
    107109      ! 
    … …  
    113115      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
    114116      END WHERE 
    115       DO jl = 1, jpl 
    116          DO_2D( 0, 0, 0, 0 ) 
    117             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    118                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    119                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    120                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    121             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    122                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    123                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    124                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    125             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    126                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    127                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    128                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    129             zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
    130                &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
    131                &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
    132                &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    133          END_2D 
    134       END DO 
     117      CALL icemax3D( ph_i , zhi_max ) 
     118      CALL icemax3D( ph_s , zhs_max ) 
     119      CALL icemax3D( ph_ip, zhip_max) 
     120      CALL icemax3D( zs_i , zsi_max ) 
    135121#if defined key_mpi3 
    136122      CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    … …  
    149135         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    150136         END WHERE 
    151       END DO 
    152       DO jl = 1, jpl 
    153          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    154             zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
    155                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    156                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    157                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    158          END_3D 
    159       END DO 
    160       DO jl = 1, jpl 
    161          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    162             zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
    163                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    164                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    165                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    166          END_3D 
    167       END DO 
    168 #if defined key_mpi3 
    169       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    170       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    171 #else 
    172       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    173       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    174 #endif 
     137      END DO    
     138      CALL icemax4D( ze_i , zei_max ) 
     139      CALL icemax4D( ze_s , zes_max ) 
     140      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 
     141      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 
    175142      ! 
    176143      ! 
    … …  
    188155      ENDIF 
    189156      zdt = rDt_ice / REAL(icycle) 
     157      z1_dt = 1._wp / zdt 
    190158 
    191159      ! --- transport --- ! 
    … …  
    216184      !---------------! 
    217185      DO jt = 1, icycle 
     186 
     187         ! diagnostics 
     188         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     189         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     190         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     191            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    218192 
    219193         ! record at_i before advection (for open water) 
    … …  
    386360            ENDIF 
    387361         ENDIF 
     362 
     363         ! --- Lateral boundary conditions --- ! 
     364         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     365            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
     366               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
     367         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     368            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
     369               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
     370         ELSE 
     371            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 
     372         ENDIF 
     373         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 
     374         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 
    388375         ! 
    389376         !== Open water area ==! 
    … …  
    399386#endif 
    400387         ! 
     388         ! --- diagnostics --- ! 
     389         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     390            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     391         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     392            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     393         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     394            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     395            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    401396         ! 
    402397         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
    … …  
    458453      !!             work on H (and not V). It is partly related to the multi-category approach 
    459454      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    460       !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
    461       !!             since sv_i and e_i are still good. 
     455      !!             concentration is small). We also limit S and T. 
    462456      !!---------------------------------------------------------------------- 
    463457      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    … …  
    503497      IF( pamsk == 0._wp ) THEN 
    504498         DO jl = 1, jpl 
    505             DO_2D( 1, 0, 1, 0 ) 
     499            DO_2D( 0, 0, 1, 0 ) 
    506500               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
    507501                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
    … …  
    512506               ENDIF 
    513507               ! 
     508            END_2D 
     509            DO_2D( 1, 0, 0, 0 ) 
    514510               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
    515511                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
    … …  
    550546      IF( PRESENT( pua_ho ) ) THEN 
    551547         DO jl = 1, jpl 
    552             DO_2D( 1, 0, 1, 0 ) 
    553                pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
    554                pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     548            DO_2D( 0, 0, 1, 0 ) 
     549               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 
     550               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 
     551            END_2D 
     552            DO_2D( 1, 0, 0, 0 ) 
     553               pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     554               pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
    555555            END_2D 
    556556         END DO 
    … …  
    566566         END_2D 
    567567      END DO 
    568 #if defined key_mpi3 
    569       CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', ptc, 'T',  1.0_wp ) 
    570 #else 
    571       CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp ) 
    572 #endif 
    573568      ! 
    574569   END SUBROUTINE adv_umx 
    … …  
    609604            ! 
    610605            DO jl = 1, jpl              !-- flux in x-direction 
    611                DO_2D( 1, 0, 1, 0 ) 
     606               DO_2D( 1, 1, 1, 0 ) 
    612607                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
    613608               END_2D 
    … …  
    615610            ! 
    616611            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    617                DO_2D( 0, 0, 0, 0 ) 
     612               DO_2D( 1, 1, 0, 0 ) 
    618613                  ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
    619614                     &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    622617               END_2D 
    623618            END DO 
    624 #if defined key_mpi3 
    625             CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    626 #else 
    627             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    628 #endif 
    629619            ! 
    630620            DO jl = 1, jpl              !-- flux in y-direction 
    631                DO_2D( 1, 0, 1, 0 ) 
     621               DO_2D( 1, 0, 0, 0 ) 
    632622                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
    633623               END_2D 
    … …  
    637627            ! 
    638628            DO jl = 1, jpl              !-- flux in y-direction 
    639                DO_2D( 1, 0, 1, 0 ) 
     629               DO_2D( 1, 0, 1, 1 ) 
    640630                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    641631               END_2D 
    … …  
    643633            ! 
    644634            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    645                DO_2D( 0, 0, 0, 0 ) 
     635               DO_2D( 0, 0, 1, 1 ) 
    646636                  ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
    647637                     &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    650640               END_2D 
    651641            END DO 
    652 #if defined key_mpi3 
    653             CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    654 #else 
    655             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    656 #endif 
    657642            ! 
    658643            DO jl = 1, jpl              !-- flux in x-direction 
    659                DO_2D( 1, 0, 1, 0 ) 
     644               DO_2D( 0, 0, 1, 0 ) 
    660645                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
    661646               END_2D 
    … …  
    710695         ! 
    711696         DO jl = 1, jpl 
    712             DO_2D( 1, 0, 1, 0 ) 
     697            DO_2D( 1, 1, 1, 0 ) 
    713698               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) ) 
     699            END_2D 
     700            DO_2D( 1, 0, 1, 1 ) 
    714701               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
    715702            END_2D 
    … …  
    728715            ! 
    729716            DO jl = 1, jpl              !-- flux in x-direction 
    730                DO_2D( 1, 0, 1, 0 ) 
     717               DO_2D( 1, 1, 1, 0 ) 
    731718                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    732719               END_2D 
    … …  
    735722 
    736723            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    737                DO_2D( 0, 0, 0, 0 ) 
     724               DO_2D( 1, 1, 0, 0 ) 
    738725                  ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
    739726                     &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    742729               END_2D 
    743730            END DO 
    744 #if defined key_mpi3 
    745             CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    746 #else 
    747             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    748 #endif 
    749731 
    750732            DO jl = 1, jpl              !-- flux in y-direction 
    751                DO_2D( 1, 0, 1, 0 ) 
     733               DO_2D( 1, 0, 0, 0 ) 
    752734                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    753735               END_2D 
    … …  
    758740            ! 
    759741            DO jl = 1, jpl              !-- flux in y-direction 
    760                DO_2D( 1, 0, 1, 0 ) 
     742               DO_2D( 1, 0, 1, 1 ) 
    761743                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    762744               END_2D 
    … …  
    765747            ! 
    766748            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    767                DO_2D( 0, 0, 0, 0 ) 
     749               DO_2D( 0, 0, 1, 1 ) 
    768750                  ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
    769751                     &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    772754               END_2D 
    773755            END DO 
    774 #if defined key_mpi3 
    775             CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    776 #else 
    777             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    778 #endif 
    779756            ! 
    780757            DO jl = 1, jpl              !-- flux in x-direction 
    781                DO_2D( 1, 0, 1, 0 ) 
     758               DO_2D( 0, 0, 1, 0 ) 
    782759                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    783760               END_2D 
    … …  
    952929         !         
    953930         DO jl = 1, jpl 
    954             DO_2D( 1, 0, 1, 0 ) 
     931            DO_2D( 0, 0, 1, 0 ) 
    955932               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    956933                  &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    … …  
    961938         ! 
    962939         DO jl = 1, jpl 
    963             DO_2D( 1, 0, 1, 0 ) 
     940            DO_2D( 0, 0, 1, 0 ) 
    964941               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    965942               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    … …  
    971948         ! 
    972949         DO jl = 1, jpl 
    973             DO_2D( 1, 0, 1, 0 ) 
     950            DO_2D( 0, 0, 1, 0 ) 
    974951               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    975952               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    985962         ! 
    986963         DO jl = 1, jpl 
    987             DO_2D( 1, 0, 1, 0 ) 
     964            DO_2D( 0, 0, 1, 0 ) 
    988965               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    989966               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    999976         ! 
    1000977         DO jl = 1, jpl 
    1001             DO_2D( 1, 0, 1, 0 ) 
     978            DO_2D( 0, 0, 1, 0 ) 
    1002979               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    1003980               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    1020997      IF( ll_neg ) THEN 
    1021998         DO jl = 1, jpl 
    1022             DO_2D( 1, 0, 1, 0 ) 
     999            DO_2D( 0, 0, 1, 0 ) 
    10231000               IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    10241001                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    … …  
    10301007      !                                                     !-- High order flux in i-direction  --! 
    10311008      DO jl = 1, jpl 
    1032          DO_2D( 1, 0, 1, 0 ) 
     1009         DO_2D( 0, 0, 1, 0 ) 
    10331010            pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
    10341011         END_2D 
    … …  
    10961073      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    10971074         DO jl = 1, jpl 
    1098             DO_2D( 1, 0, 1, 0 ) 
     1075            DO_2D( 1, 0, 0, 0 ) 
    10991076               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    11001077                  &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    11041081      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    11051082         DO jl = 1, jpl 
    1106             DO_2D( 1, 0, 1, 0 ) 
     1083            DO_2D( 1, 0, 0, 0 ) 
    11071084               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    11081085               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    … …  
    11131090      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    11141091         DO jl = 1, jpl 
    1115             DO_2D( 1, 0, 1, 0 ) 
     1092            DO_2D( 1, 0, 0, 0 ) 
    11161093               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    11171094               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    11261103      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    11271104         DO jl = 1, jpl 
    1128             DO_2D( 1, 0, 1, 0 ) 
     1105            DO_2D( 1, 0, 0, 0 ) 
    11291106               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    11301107               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    11391116      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    11401117         DO jl = 1, jpl 
    1141             DO_2D( 1, 0, 1, 0 ) 
     1118            DO_2D( 1, 0, 0, 0 ) 
    11421119               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    11431120               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    11601137      IF( ll_neg ) THEN 
    11611138         DO jl = 1, jpl 
    1162             DO_2D( 1, 0, 1, 0 ) 
     1139            DO_2D( 1, 0, 0, 0 ) 
    11631140               IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    11641141                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    … …  
    11701147      !                                                     !-- High order flux in j-direction  --! 
    11711148      DO jl = 1, jpl 
    1172          DO_2D( 1, 0, 1, 0 ) 
     1149         DO_2D( 1, 0, 0, 0 ) 
    11731150            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    11741151         END_2D 
    … …  
    12061183      ! -------------------------------------------------- 
    12071184      DO jl = 1, jpl 
    1208          DO_2D( 1, 0, 1, 0 ) 
     1185         DO_2D( 0, 0, 1, 0 ) 
    12091186            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
     1187         END_2D 
     1188         DO_2D( 1, 0, 0, 0 ) 
    12101189            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    12111190         END_2D 
    … …  
    13251304      ! --------------------------------- 
    13261305      DO jl = 1, jpl 
    1327          DO_2D( 1, 0, 1, 0 ) 
     1306         DO_2D( 0, 0, 1, 0 ) 
    13281307            zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    13291308            zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
    … …  
    13361315         END_2D 
    13371316 
    1338          DO_2D( 1, 0, 1, 0 ) 
     1317         DO_2D( 1, 0, 0, 0 ) 
    13391318            zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    13401319            zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
    … …  
    17091688   END SUBROUTINE Hsnow 
    17101689 
     1690   SUBROUTINE icemax3D( pice , pmax ) 
     1691      !!--------------------------------------------------------------------- 
     1692      !!                   ***  ROUTINE icemax3D ***                      
     1693      !! ** Purpose :  compute the max of the 9 points around 
     1694      !!---------------------------------------------------------------------- 
     1695      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1696      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1697      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1698      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1699      !!---------------------------------------------------------------------- 
     1700      DO jl = 1, jpl 
     1701         DO jj = Njs0-1, Nje0+1     
     1702            DO ji = Nis0, Nie0 
     1703               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1704            END DO 
     1705         END DO 
     1706         DO jj = Njs0, Nje0     
     1707            DO ji = Nis0, Nie0 
     1708               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1709            END DO 
     1710         END DO 
     1711      END DO 
     1712   END SUBROUTINE icemax3D 
     1713 
     1714   SUBROUTINE icemax4D( pice , pmax ) 
     1715      !!--------------------------------------------------------------------- 
     1716      !!                   ***  ROUTINE icemax4D ***                      
     1717      !! ** Purpose :  compute the max of the 9 points around 
     1718      !!---------------------------------------------------------------------- 
     1719      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1720      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1721      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1722      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1723      !!---------------------------------------------------------------------- 
     1724      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1725      DO jl = 1, jpl 
     1726         DO jk = 1, jlay 
     1727            DO jj = Njs0-1, Nje0+1     
     1728               DO ji = Nis0, Nie0 
     1729                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1730               END DO 
     1731            END DO 
     1732            DO jj = Njs0, Nje0     
     1733               DO ji = Nis0, Nie0 
     1734                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1735               END DO 
     1736            END DO 
     1737         END DO 
     1738      END DO 
     1739   END SUBROUTINE icemax4D 
    17111740 
    17121741#else 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rdgrft.F90

    r13630 r13877  
    349349               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 
    350350                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  & 
    351                      &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar ) 
     351                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar     ) * z1_gstar ) 
    352352               ELSE 
    353353                  apartf(ji,jl) = 0._wp 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90

    r13807 r13877  
    148148      ! 
    149149      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     150      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension 
    150151      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    151152      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    … …  
    170171      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    171172      !! --- diags 
    172       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
     173      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
     174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p          
    173175      !! --- SIMIP diags 
    174176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    … …  
    760762            &   ) * r1_e1e2t(ji,jj) 
    761763         zdt2 = zdt * zdt 
     764 
     765         zten_i(ji,jj) = zdt 
    762766          
    763767         ! shear**2 at T points (doc eq. A16) 
    … …  
    775779          
    776780         ! delta at T points 
    777          zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    778          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0 
    779          pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch 
     781         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta   
     782         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 
     783         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 
    780784 
    781785      END_2D 
    782786#if defined key_mpi3 
    783       CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 
    784 #else 
    785       CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 
     787      CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 
     788         &                                  zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp ) 
     789#else 
     790      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 
     791         &                                  zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp ) 
    786792#endif 
    787793       
    788794      ! --- Store the stress tensor for the next time step --- ! 
    789 #if defined key_mpi3 
    790       CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp ) 
    791 #else 
    792       CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp ) 
    793 #endif 
    794795      pstress1_i (:,:) = zs1 (:,:) 
    795796      pstress2_i (:,:) = zs2 (:,:) 
    … …  
    825826      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    826827 
    827       ! --- stress tensor --- ! 
    828       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    829          ! 
    830          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     828      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     829      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     830         ! 
     831         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    831832         !          
    832          DO_2D( 0, 0, 0, 0 ) 
    833             zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    834                &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    835                &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    836  
    837             zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    838  
    839             zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    840  
    841 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    842 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    843 !!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    844 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    845             zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    846             zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    847             zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    848          END_2D 
    849 #if defined key_mpi3 
    850          CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 
    851 #else 
    852          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 
    853 #endif 
    854          ! 
    855          CALL iom_put( 'isig1' , zsig1 ) 
    856          CALL iom_put( 'isig2' , zsig2 ) 
    857          CALL iom_put( 'isig3' , zsig3 ) 
    858          ! 
    859          ! Stress tensor invariants (normal and shear stress N/m) 
    860          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    861          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    862  
    863          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     833         DO_2D( 1, 1, 1, 1 ) 
     834             
     835            ! Ice stresses 
     836            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     837            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
     838            ! I know, this can be confusing... 
     839            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     840            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     841            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     842            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     843             
     844            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
     845            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     846            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     847                
     848         END_2D          
     849         ! 
     850         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
     851         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
     852         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     853          
     854         DEALLOCATE ( zsig_I, zsig_II ) 
     855          
     856      ENDIF 
     857 
     858      ! --- Normalized stress tensor principal components --- ! 
     859      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 
     860      ! Recommendation 1 : we use ice strength, not replacement pressure 
     861      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
     862      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     863         ! 
     864         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
     865         !          
     866         DO_2D( 1, 1, 1, 1 ) 
     867             
     868            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     869            !                        and **deformations** at current iterates 
     870            !                        following Lemieux & Dupont (2020) 
     871            zfac             =   zp_delt(ji,jj) 
     872            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
     873            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     874            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     875             
     876            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     877            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                            ! 1st stress invariant, aka average normal stress, aka negative pressure 
     878            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )   ! 2nd  ''       '', aka maximum shear stress 
     879             
     880            ! Normalized  principal stresses (used to display the ellipse) 
     881            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     882            zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     883            zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     884         END_2D               
     885         ! 
     886         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     887         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     888 
     889         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     890          
    864891      ENDIF 
    865892 
    … …  
    10081035         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
    10091036         ! close file 
    1010          IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     1037         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
    10111038      ENDIF 
    10121039       

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceitd.F90

    r13571 r13877  
    627627         END_2D 
    628628         ! 
    629 !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
    630          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
    631          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
    632          ! 
    633          DO ji = 1, npti 
    634             jdonor(ji,jl)  = jl  
    635             ! how much of a_i you send in cat sup is somewhat arbitrary 
    636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    637 !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    638 !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    640 !!          zdaice(ji,jl)  = a_i_1d(ji) 
    641 !!          zdvice(ji,jl)  = v_i_1d(ji) 
    642 !!clem: these are from UCL and work ok 
    643             zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    644             zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
    645          END DO 
    646          ! 
    647          IF( npti > 0 ) THEN 
     629         IF( npti > 0 ) THEN             
     630            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
     631            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     632            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 
     633            ! 
     634            DO ji = 1, npti 
     635               jdonor(ji,jl)  = jl  
     636               ! how much of a_i you send in cat sup is somewhat arbitrary 
     637               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     638               !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
     639               !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
     640               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
     641               !!          zdaice(ji,jl)  = a_i_1d(ji) 
     642               !!          zdvice(ji,jl)  = v_i_1d(ji) 
     643               !!clem: these are from UCL and work ok 
     644               zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
     645               zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     646            END DO 
     647            ! 
    648648            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1 
    649649            ! Reset shift parameters 
    … …  
    666666         END_2D 
    667667         ! 
    668          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
    669          CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
    670          DO ji = 1, npti 
    671             jdonor(ji,jl) = jl + 1 
    672             zdaice(ji,jl) = a_i_1d(ji)  
    673             zdvice(ji,jl) = v_i_1d(ji) 
    674          END DO 
    675          ! 
    676668         IF( npti > 0 ) THEN 
     669            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 
     670            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 
     671            DO ji = 1, npti 
     672               jdonor(ji,jl) = jl + 1 
     673               zdaice(ji,jl) = a_i_1d(ji)  
     674               zdvice(ji,jl) = v_i_1d(ji) 
     675            END DO 
     676            ! 
    677677            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl 
    678678            ! Reset shift parameters 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icestp.F90

    r13571 r13877  
    5555   USE icedyn         ! sea-ice: dynamics 
    5656   USE icethd         ! sea-ice: thermodynamics 
    57    USE icecor         ! sea-ice: corrections 
    5857   USE iceupdate      ! sea-ice: sea surface boundary condition update 
    5958   USE icedia         ! sea-ice: budget diagnostics 
    … …  
    8685   PUBLIC   ice_init   ! called by sbcmod.F90 
    8786 
     87   !! * Substitutions 
     88#  include "do_loop_substitute.h90" 
    8889   !!---------------------------------------------------------------------- 
    8990   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    160161         IF( ln_icedyn .AND. .NOT.lk_c1d )   & 
    161162            &                           CALL ice_dyn( kt, Kmm )       ! -- Ice dynamics 
     163         ! 
     164                                        CALL diag_trends( 1 )         ! record dyn trends 
    162165         ! 
    163166         !                          !==  lateral boundary conditions  ==! 
    … …  
    188191         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    189192         ! 
    190                                         CALL ice_cor( kt , 2 )        ! -- Corrections 
    191          ! 
     193                                        CALL diag_trends( 2 )         ! record thermo trends 
    192194                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling) 
    193195                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling) 
    … …  
    196198         ! 
    197199         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
     200         ! 
     201         IF( ln_icediachk )             CALL ice_drift_wri( kt )      ! -- Diagnostics outputs for conservation  
    198202         ! 
    199203                                        CALL ice_wri( kt )            ! -- Ice outputs  
    … …  
    208212      ! --- Ocean time step --- ! 
    209213      !-------------------------! 
    210       IF( ln_icedyn )                   CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )   ! -- update surface ocean stresses 
     214      CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )         ! -- update surface ocean stresses 
    211215!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    212216      ! 
    … …  
    281285      ! 
    282286      CALL ice_dia_init                ! initialization for diags 
     287      ! 
     288      CALL ice_drift_init              ! initialization for diags of conservation 
    283289      ! 
    284290      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
    … …  
    341347      ENDIF 
    342348      ! 
    343       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    344       ! 
    345349      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse 
    346350      r1_Dt_ice = 1._wp / rDt_ice 
    … …  
    392396      !!               of the time step 
    393397      !!---------------------------------------------------------------------- 
    394       INTEGER  ::   ji, jj      ! dummy loop index 
    395       !!---------------------------------------------------------------------- 
    396       sfx    (:,:) = 0._wp   ; 
    397       sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
    398       sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    399       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    400       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    401       sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
    402       ! 
    403       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    404       wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    405       wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    406       wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    407       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    408       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
    409       wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 
    410       wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 
    411       wfx_snw_sni(:,:) = 0._wp  
    412       wfx_pnd(:,:) = 0._wp 
    413  
    414       hfx_thd(:,:) = 0._wp   ; 
    415       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    416       hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
    417       hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
    418       hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    419       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
    420       hfx_err_dif(:,:) = 0._wp 
    421       wfx_err_sub(:,:) = 0._wp 
    422       ! 
    423       diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp 
    424       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp 
    425  
    426       ! SIMIP diagnostics 
    427       qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 
    428       t_si       (:,:,:) = rt0   ! temp at the ice-snow interface 
    429  
    430       tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    431       cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
    432       qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 
    433       qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
    434       qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F) 
    435       ! 
    436       ! for control checks (ln_icediachk) 
    437       diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp 
    438       diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp 
    439       diag_trp_sv(:,:) = 0._wp 
    440        
     398      INTEGER  ::   ji, jj, jl      ! dummy loop index 
     399      !!---------------------------------------------------------------------- 
     400 
     401      DO_2D( 1, 1, 1, 1 ) 
     402         sfx    (ji,jj) = 0._wp   ; 
     403         sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp 
     404         sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp 
     405         sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp 
     406         sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp 
     407         sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp 
     408         ! 
     409         wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp 
     410         wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp 
     411         wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp 
     412         wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp 
     413         wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp 
     414         wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp   
     415         wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 
     416         wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 
     417         wfx_snw_sni(ji,jj) = 0._wp  
     418         wfx_pnd(ji,jj) = 0._wp 
     419 
     420         hfx_thd(ji,jj) = 0._wp   ; 
     421         hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 
     422         hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp 
     423         hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp 
     424         hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp 
     425         hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp 
     426         hfx_err_dif(ji,jj) = 0._wp 
     427         wfx_err_sub(ji,jj) = 0._wp 
     428         ! 
     429         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
     430         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
     431 
     432         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     433         qsb_ice_bot(ji,jj) = 0._wp   ! (needed if ln_icethd=F) 
     434 
     435         fhld(ji,jj) = 0._wp   ! needed if ln_icethd=F 
     436 
     437         ! for control checks (ln_icediachk) 
     438         diag_trp_vi(ji,jj) = 0._wp   ;   diag_trp_vs(ji,jj) = 0._wp 
     439         diag_trp_ei(ji,jj) = 0._wp   ;   diag_trp_es(ji,jj) = 0._wp 
     440         diag_trp_sv(ji,jj) = 0._wp 
     441         ! 
     442         diag_adv_mass(ji,jj) = 0._wp 
     443         diag_adv_salt(ji,jj) = 0._wp 
     444         diag_adv_heat(ji,jj) = 0._wp 
     445      END_2D 
     446 
     447      DO jl = 1, jpl 
     448         DO_2D( 1, 1, 1, 1 ) 
     449            ! SIMIP diagnostics 
     450            t_si       (ji,jj,jl) = rt0     ! temp at the ice-snow interface 
     451            qcn_ice_bot(ji,jj,jl) = 0._wp 
     452            qcn_ice_top(ji,jj,jl) = 0._wp   ! conductive fluxes 
     453            cnd_ice    (ji,jj,jl) = 0._wp   ! effective conductivity at the top of ice/snow (ln_cndflx=T) 
     454            qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T) 
     455            qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs 
     456         END_2D 
     457      ENDDO 
     458 
    441459   END SUBROUTINE diag_set0 
     460 
     461 
     462   SUBROUTINE diag_trends( kn ) 
     463      !!---------------------------------------------------------------------- 
     464      !!                  ***  ROUTINE diag_trends  *** 
     465      !! 
     466      !! ** purpose : diagnostics of the trends. Used for conservation purposes 
     467      !!              and outputs 
     468      !!---------------------------------------------------------------------- 
     469      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo 
     470      !!---------------------------------------------------------------------- 
     471      ! 
     472      ! --- trends of heat, salt, mass (used for conservation controls) 
     473      IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
     474         ! 
     475         diag_heat(:,:) = diag_heat(:,:) & 
     476            &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 
     477            &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
     478         diag_sice(:,:) = diag_sice(:,:) & 
     479            &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     480         diag_vice(:,:) = diag_vice(:,:) & 
     481            &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
     482         diag_vsnw(:,:) = diag_vsnw(:,:) & 
     483            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
     484         ! 
     485         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 
     486         ! 
     487      ENDIF 
     488      ! 
     489      ! --- trends of concentration (used for simip outputs) 
     490      IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
     491         ! 
     492         diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
     493         ! 
     494         IF( kn == 1 )   CALL iom_put( 'afxdyn' , diag_aice )                                           ! dyn trend 
     495         IF( kn == 2 )   CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 
     496         IF( kn == 2 )   CALL iom_put( 'afxtot' , diag_aice )                                           ! total trend 
     497         ! 
     498      ENDIF 
     499      ! 
     500   END SUBROUTINE diag_trends 
    442501 
    443502#else 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icetab.F90

    r10069 r13877  
    4040      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index 
    4141      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field 
    42       REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field 
     42      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab1d    ! output 1D field 
    4343      ! 
    4444      INTEGER ::   jl, jn, jid, jjd 
    … …  
    6161      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind  ! input index 
    6262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   tab2d    ! input 2D field 
    63       REAL(wp), DIMENSION(ndim1d) , INTENT(  out) ::   tab1d    ! output 1D field 
     63      REAL(wp), DIMENSION(ndim1d) , INTENT(inout) ::   tab1d    ! output 1D field 
    6464      ! 
    6565      INTEGER ::   jn , jid, jjd 
    … …  
    8080      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index 
    8181      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field 
    82       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field 
     82      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab2d     ! output 2D field 
    8383      ! 
    8484      INTEGER ::   jl, jn, jid, jjd 
    … …  
    101101      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind   ! input index 
    102102      REAL(wp), DIMENSION(ndim1d) , INTENT(in   ) ::   tab1d     ! input 1D field 
    103       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   tab2d     ! output 2D field 
     103      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   tab2d     ! output 2D field 
    104104      ! 
    105105      INTEGER ::   jn , jid, jjd 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90

    r13630 r13877  
    1818   USE ice            ! sea-ice: variables 
    1919!!gm list trop longue ==>>> why not passage en argument d'appel ? 
    20    USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl 
     20   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 
    2121   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    2222      &                 qml_ice, qcn_ice, qtr_ice_top 
    … …  
    3030   USE icethd_pnd     ! sea-ice: melt ponds 
    3131   USE iceitd         ! sea-ice: remapping thickness distribution 
     32   USE icecor         ! sea-ice: corrections 
    3233   USE icetab         ! sea-ice: 1D <==> 2D transformation 
    3334   USE icevar         ! sea-ice: operations 
    … …  
    5253   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F) 
    5354   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
    54    LOGICAL ::   ln_leadhfx       !  heat in the leads is used to melt sea-ice before warming the ocean 
     55   LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean 
    5556 
    5657   !! for convergence tests 
    … …  
    9192      ! 
    9293      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    93       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 
    94       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
    95       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
    96       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
     94      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 
     95      REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
     96      REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient 
     97      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
    9798      ! 
    9899      !!------------------------------------------------------------------- 
    … …  
    124125               &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    125126               &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 
     127            zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 
     128               &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
    126129         END_2D 
    127130      ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean 
    … …  
    130133               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    131134               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     135            zvel(ji,jj) = 0._wp 
    132136         END_2D 
    133137      ENDIF 
    134138#if defined key_mpi3 
    135       CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T',  1.0_wp ) 
     139      CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
    136140#else 
    137       CALL lbc_lnk( 'icethd', zfric, 'T', 1.0_wp ) 
     141      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
    138142#endif 
    139143      ! 
    … …  
    144148         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    145149         ! 
    146          !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    147          !           !  practically no "direct lateral ablation" 
    148          !            
    149          !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    150          !           !  temperature and turbulent mixing (McPhee, 1992) 
    151          ! 
    152150         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    153151         zqld =  tmask(ji,jj,1) * rDt_ice *  & 
    … …  
    155153            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    156154 
    157          ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 
     155         ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
     156         !     (mostly<0 but >0 if supercooling) 
    158157         zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 
    159158         zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
    160  
    161          ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 
     159         zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
     160 
     161         ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
     162         !     (mostly>0 but <0 if supercooling) 
    162163         zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin )  
    163          qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
    164  
    165          qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     164         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
     165          
    166166         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    167167         !                              the freezing point, so that we do not have SST < T_freeze 
    168          !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 
    169  
    170          !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 
    171          qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    172  
    173          ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
    174          ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
    175          IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    176             IF( ln_leadhfx ) THEN   ;   fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    177             ELSE                    ;   fhld(ji,jj) = 0._wp 
    178             ENDIF 
     168         !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
     169         !                              The following formulation is ok for both normal conditions and supercooling 
     170         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     171 
     172         ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
     173         !     qlead is the energy received from the atm. in the leads. 
     174         !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
     175         !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
     176         IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     177            ! upper bound for fhld: fhld should be equal to zqld 
     178            !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
     179            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
     180            !                        The following formulation is ok for both normal conditions and supercooling 
     181            fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
     182               &                                 - qsb_ice_bot(ji,jj) ) 
    179183            qlead(ji,jj) = 0._wp 
    180184         ELSE 
    181185            fhld (ji,jj) = 0._wp 
     186            ! upper bound for qlead: qlead should be equal to zqld 
     187            !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
     188            !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
     189            !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
     190            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
     191            !                        The following formulation is ok for both normal conditions and supercooling 
     192            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    182193         ENDIF 
    183194         ! 
    184          ! Net heat flux on top of the ice-ocean [W.m-2] 
    185          ! --------------------------------------------- 
    186          qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
     195         ! If ice is landfast and ice concentration reaches its max 
     196         ! => stop ice formation in open water 
     197         IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
     198         ! 
     199         ! If the grid cell is almost fully covered by ice (no leads) 
     200         ! => stop ice formation in open water 
     201         IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
     202         ! 
     203         ! If ln_leadhfx is false 
     204         ! => do not use energy of the leads to melt sea-ice 
     205         IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
     206         ! 
    187207      END_2D 
    188208       
    … …  
    195215      ENDIF 
    196216 
    197       ! --------------------------------------------------------------------- 
    198       ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    199       ! --------------------------------------------------------------------- 
    200       !     First  step here              :  non solar + precip - qlead - qsensible 
    201       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)  
    202       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    203       qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    204          &             - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    205          &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux 
    206          &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt  
    207       !                                                                           !    (fhld should be 0 while bott growth) 
    208217      !-------------------------------------------------------------------------------------------! 
    209218      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
    … …  
    259268      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
    260269      ! 
     270                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
     271      ! 
     272      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation      
     273      ! 
    261274      ! convergence tests 
    262275      IF( ln_zdf_chkcvg ) THEN 
    … …  
    422435         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       ) 
    423436         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   ) 
    424          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     ) 
    425437         ! 
    426438         ! ocean surface fields 
    427439         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 
    428440         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 
     441         CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 
    429442         ! 
    430443         ! to update ice age 
    … …  
    514527         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     ) 
    515528         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 
    516          CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   ) 
    517529         ! 
    518530         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) ) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_dh.F90

    r13571 r13877  
    139139      ! 
    140140      DO ji = 1, npti 
    141          zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji)  
     141         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji)  
    142142         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 
    143143      END DO 
    … …  
    556556         !     
    557557         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    558          qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 
     558         !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 
    559559 
    560560         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_do.F90

    r13630 r13877  
    131131 
    132132      ! Default new ice thickness 
    133       WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 
    134       ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp 
     133      WHERE( qlead(:,:) < 0._wp ) ! cooling 
     134         ht_i_new(:,:) = rn_hinew 
     135      ELSEWHERE 
     136         ht_i_new(:,:) = 0._wp 
    135137      END WHERE 
    136138 
    … …  
    146148         ! 
    147149         DO_2D( 0, 0, 0, 0 ) 
    148             IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 
     150            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    149151               ! -- Wind stress -- ! 
    150152               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
    … …  
    202204      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 
    203205      !------------------------------------------------------------------------------! 
    204       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 
     206      ! it occurs if cooling 
    205207 
    206208      ! Identify grid points where new ice forms 
    207209      npti = 0   ;   nptidx(:) = 0 
    208210      DO_2D( 1, 1, 1, 1 ) 
    209          IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 
     211         IF ( qlead(ji,jj)  <  0._wp ) THEN 
    210212            npti = npti + 1 
    211213            nptidx( npti ) = (jj - 1) * jpi + ji 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceupdate.F90

    r13630 r13877  
    2424   USE traqsr         ! add penetration of solar flux in the calculation of heat budget 
    2525   USE icectl         ! sea-ice: control prints 
    26    USE bdy_oce , ONLY : ln_bdy 
    2726   USE zdfdrg  , ONLY : ln_drgice_imp 
    2827   ! 
    … …  
    9291      ! 
    9392      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    94       REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    9593      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
    9694      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
    … …  
    103101         WRITE(numout,*)'~~~~~~~~~~~~~~' 
    104102      ENDIF 
     103 
     104      ! Net heat flux on top of the ice-ocean (W.m-2) 
     105      !---------------------------------------------- 
     106      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:)  
    105107 
    106108      ! --- case we bypass ice thermodynamics --- ! 
    … …  
    115117      DO_2D( 1, 1, 1, 1 ) 
    116118 
    117          ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
     119         ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)  
    118120         !--------------------------------------------------- 
    119121         zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) 
    … …  
    121123         ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2)  
    122124         !--------------------------------------------------- 
    123          zqmass           = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
    124          qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 
    125  
    126          ! Add the residual from heat diffusion equation and sublimation (W.m-2) 
    127          !---------------------------------------------------------------------- 
    128          qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) +   & 
    129             &             ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 
    130  
     125         qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 
     126            &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 
     127            &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 
     128            &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)                  
     129          
    131130         ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    132131         !---------------------------------------------------------------------------- 
    133          qsr(ji,jj) = zqsr                                       
     132         ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 
     133         ! else ( cooling or no ice left ), then we suppose that     no    solar flux has been consumed 
     134         ! 
     135         IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN   !-- warming and some ice remains 
     136            !                                        solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 
     137            qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 
     138               !                                   + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice) 
     139               &             + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) ) 
     140            ! 
     141         ELSE                                                       !-- cooling or no ice left 
     142            qsr(ji,jj) = zqsr 
     143         ENDIF 
     144         ! 
     145         ! the non-solar is simply derived from the solar flux 
    134146         qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr               
    135  
     147          
    136148         ! Mass flux at the atm. surface        
    137149         !----------------------------------- 
    … …  
    140152         ! Mass flux at the ocean surface       
    141153         !------------------------------------ 
    142          !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 
    143          !  -------------------------------------------------------------------------------------  
    144          !  The idea of this approach is that the system that we consider is the ICE-OCEAN system 
    145          !  Thus  FW  flux  =  External ( E-P+snow melt) 
    146          !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW 
    147          !                     Associated to Ice formation AND Ice melting 
    148          !                     Even if i see Ice melting as a FW and SALT flux 
    149          !         
    150          ! mass flux from ice/ocean 
     154         ! ice-ocean  mass flux 
    151155         wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    152156            &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
    153  
    154          ! add the snow melt water to snow mass flux to the ocean 
     157          
     158         ! snw-ocean mass flux 
    155159         wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 
    156  
    157          ! mass flux at the ocean/ice interface 
    158          fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model 
    159          emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    160  
     160          
     161         ! total mass flux at the ocean/ice interface 
     162         fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model 
     163         emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux 
    161164 
    162165         ! Salt flux at the ocean surface       
    … …  
    262265      CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change 
    263266      CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt  
    264       CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai) 
     267      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion 
    265268 
    266269      ! heat fluxes associated with mass exchange (freeze/melt/precip...) 
    … …  
    279282      !--------- 
    280283#if ! defined key_agrif 
    281       IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation 
     284      IF( ln_icediachk      )   CALL ice_cons_final('iceupdate')                                       ! conservation 
    282285#endif 
    283       IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    284       IF( sn_cfctl%l_prtctl              )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    285       IF( ln_timing                      )   CALL timing_stop   ('ice_update')                                      ! timing 
     286      IF( ln_icectl         )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
     287      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
     288      IF( ln_timing         )   CALL timing_stop   ('ice_update')                                      ! timing 
    286289      ! 
    287290   END SUBROUTINE ice_update_flx 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/BDY/bdyice.F90

    r13630 r13877  
    6161      !!---------------------------------------------------------------------- 
    6262      ! controls 
    63       IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing 
    64       IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    65       IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
     63      IF( ln_timing )   CALL timing_start('bdy_ice_thd')   ! timing 
    6664      ! 
    6765      CALL ice_var_glo2eqv 
    … …  
    120118      ! 
    121119      ! controls 
    122       IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
    123       IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    124       IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    125       IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
     120      IF( ln_icectl )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )   ! prints 
     121      IF( ln_timing )   CALL timing_stop ('bdy_ice_thd')                                       ! timing 
    126122      ! 
    127123   END SUBROUTINE bdy_ice 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DIU/diu_coolskin.F90

    r13295 r13877  
    9595      !!---------------------------------------------------------------------- 
    9696      ! 
    97       IF( .NOT. ln_blk )   CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing") 
     97      IF( .NOT. (ln_blk .OR. ln_abl) )   CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing") 
    9898      ! 
    9999      DO_2D( 1, 1, 1, 1 ) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/IOM/iom.F90

    r13807 r13877  
    19351935      IF( iom_use(cdname) ) THEN 
    19361936#if defined key_iomput 
    1937          IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 
    1938             CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d 
    1939          ELSE 
    1940             CALL xios_send_field( cdname, pfield2d ) 
    1941          ENDIF 
     1937         CALL xios_send_field( cdname, pfield2d ) 
    19421938#else 
    19431939         WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19511947      IF( iom_use(cdname) ) THEN 
    19521948#if defined key_iomput 
    1953          IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 
    1954             CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d 
    1955          ELSE 
    1956             CALL xios_send_field( cdname, pfield2d ) 
    1957          ENDIF 
     1949         CALL xios_send_field( cdname, pfield2d ) 
    19581950#else 
    19591951         WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19671959      IF( iom_use(cdname) ) THEN 
    19681960#if defined key_iomput 
    1969          IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 
    1970             CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d 
    1971          ELSE 
    1972             CALL xios_send_field( cdname, pfield3d ) 
    1973          ENDIF 
     1961         CALL xios_send_field( cdname, pfield3d ) 
    19741962#else 
    19751963         WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19831971      IF( iom_use(cdname) ) THEN 
    19841972#if defined key_iomput 
    1985          IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 
    1986             CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d 
    1987          ELSE 
    1988             CALL xios_send_field( cdname, pfield3d ) 
    1989          ENDIF 
     1973         CALL xios_send_field( cdname, pfield3d ) 
    19901974#else 
    19911975         WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19991983      IF( iom_use(cdname) ) THEN 
    20001984#if defined key_iomput 
    2001          IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 
    2002             CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d 
    2003          ELSE 
    2004             CALL xios_send_field (cdname, pfield4d ) 
    2005          ENDIF 
     1985         CALL xios_send_field (cdname, pfield4d ) 
    20061986#else 
    20071987         WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    20151995      IF( iom_use(cdname) ) THEN 
    20161996#if defined key_iomput 
    2017          IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 
    2018             CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d 
    2019          ELSE 
    2020             CALL xios_send_field (cdname, pfield4d ) 
    2021          ENDIF 
     1997         CALL xios_send_field (cdname, pfield4d ) 
    20221998#else 
    20231999         WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    22252201      ! 
    22262202      CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0) 
    2227       CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0) 
     2203      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 
    22282204!don't define lon and lat for restart reading context.  
    22292205      IF ( .NOT.ldrxios ) & 
    … …  
    23282304      CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots) 
    23292305      CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0) 
    2330       CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0) 
     2306      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 
    23312307      CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp),   & 
    23322308         &                             latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))   
    2333       CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj_0) 
     2309      CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo) 
    23342310      ! 
    23352311      CALL iom_update_file_name('ptr') 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LBC/lib_mpp.F90

    r13571 r13877  
    516516            ALLOCATE(todelay(idvar)%y1d(isz)) 
    517517            todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp)   ! create %y1d, complex variable needed by mpi_sumdd 
     518            ndelayid(idvar) = MPI_REQUEST_NULL                             ! initialised request to a valid value 
    518519         END IF 
    519520      ENDIF 
    … …  
    523524         ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz))   ! allocate also %z1d as used for the restart 
    524525         CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr )   ! get %y1d 
    525          todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp)      ! define %z1d from %y1d 
    526       ENDIF 
    527  
    528       IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
     526         ndelayid(idvar) = MPI_REQUEST_NULL 
     527      ENDIF 
     528 
     529      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
    529530 
    530531      ! send back pout from todelay(idvar)%z1d defined at previous call 
    … …  
    535536      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    536537      CALL  mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) 
    537       ndelayid(idvar) = 1 
     538      ndelayid(idvar) = MPI_REQUEST_NULL 
    538539      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    539540# else 
    … …  
    596597            DEALLOCATE(todelay(idvar)%z1d) 
    597598            ndelayid(idvar) = -1                                      ! do as if we had no restart 
     599         ELSE 
     600            ndelayid(idvar) = MPI_REQUEST_NULL 
    598601         END IF 
    599602      ENDIF 
    … …  
    603606         ALLOCATE(todelay(idvar)%z1d(isz)) 
    604607         CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr )   ! get %z1d 
    605       ENDIF 
    606  
    607       IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
     608         ndelayid(idvar) = MPI_REQUEST_NULL 
     609      ENDIF 
     610 
     611      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
    608612 
    609613      ! send back pout from todelay(idvar)%z1d defined at previous call 
    … …  
    611615 
    612616      ! send p_in into todelay(idvar)%z1d with a non-blocking communication 
     617      ! (PM) Should we get rid of MPI2 option ? MPI3 was release in 2013. Who is still using MPI2 ? 
    613618# if defined key_mpi2 
    614619      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    615       CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ndelayid(idvar), ierr ) 
     620      CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ierr ) 
    616621      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    617622# else 
    … …  
    636641      !!---------------------------------------------------------------------- 
    637642#if defined key_mpp_mpi 
    638       IF( ndelayid(kid) /= -2 ) THEN   
    639 #if ! defined key_mpi2 
    640          IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    641          CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr )                        ! make sure todelay(kid) is received 
    642          IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    643 #endif 
    644          IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d 
    645          ndelayid(kid) = -2   ! add flag to know that mpi_wait was already called on kid 
    646       ENDIF 
     643      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
     644      ! test on ndelayid(kid) useless as mpi_wait return immediatly if the request handle is MPI_REQUEST_NULL 
     645      CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! after this ndelayid(kid) = MPI_REQUEST_NULL 
     646      IF( ln_timing ) CALL tic_tac( .FALSE., ld_global = .TRUE.) 
     647      IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d 
    647648#endif 
    648649   END SUBROUTINE mpp_delay_rcv 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcfwb.F90

    r13630 r13877  
    9494         snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass 
    9595         snwice_mass  (:,:) = 0.e0 
     96         snwice_fmass (:,:) = 0.e0 
    9697#endif 
    9798         ! 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90

    r13630 r13877  
    229229      CASE DEFAULT                     !- not supported 
    230230      END SELECT 
    231       IF( ln_diurnal .AND. .NOT. ln_blk )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" ) 
     231      IF( ln_diurnal .AND. .NOT. (ln_blk.OR.ln_abl) )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" ) 
    232232      ! 
    233233      !                       !**  allocate and set required variables 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/stpctl.F90

    r13571 r13877  
    6767      REAL(wp)                        ::   zzz                                   ! local real  
    6868      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
    69       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     69      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 
    7070      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
    7171      CHARACTER(len=20)               ::   clname 
    … …  
    125125      ! 
    126126      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain 
     127      ! 
     128      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain? 
     129      ! 
    127130      IF( ll_wd ) THEN 
    128131         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
    … …  
    149152      ENDIF 
    150153      zmax(9) = REAL( nstop, wp )                                                 ! stop indicator 
     154      ! 
    151155      !                                   !==               get global extrema             ==! 
    152156      !                                   !==  done by all processes if writting run.stat  ==! 
    153157      IF( ll_colruns ) THEN 
    154158         zmaxlocal(:) = zmax(:) 
    155          CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
     159         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true.  
    156160         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
    157       ENDIF 
     161      ELSE 
     162         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 
     163         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values... 
     164      ENDIF 
     165      ! 
     166      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage!  
     167      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage! 
     168      IF( ll_colruns ) THEN 
     169         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage!  
     170         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage! 
     171      ENDIF 
     172      ! 
    158173      !                                   !==              write "run.stat" files              ==! 
    159174      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    160175      IF( ll_wrtruns ) THEN 
    161          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
    162          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    163          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    164          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
    165          istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 
    166          istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 
    167          istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 
    168          IF( ln_zad_Aimp ) THEN 
    169             istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 
    170             istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 
    171          ENDIF 
     176         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 
     177         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 
     178            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 
     179         END DO 
    172180         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
    173181      END IF 
    … …  
    175183      !                                   !==  done by all processes at every time step  ==! 
    176184      ! 
    177       IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
    178          &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
    179          &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
    180          &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
    181          &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
    182          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
    183          &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     185      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
     186         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
     187         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
     188         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
     189         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
     190         & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
     191         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
    184192         ! 
    185193         iloc(:,:) = 0 
    … …  
    221229         ! 
    222230         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    223          CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
    224          CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
    225          CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
    226          CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
     231         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
     232         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
     233         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
     234         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
    227235         IF( Agrif_Root() ) THEN 
    228236            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SAS/stpctl.F90

    r13571 r13877  
    6868      REAL(wp)                        ::   zzz                                   ! local real  
    6969      REAL(wp), DIMENSION(4)          ::   zmax, zmaxlocal 
    70       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     70      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 
    7171      LOGICAL, DIMENSION(jpi,jpj)     ::   llmsk 
    7272      CHARACTER(len=20)               ::   clname 
    … …  
    124124      llmsk(:,Nje1: jpj) = .FALSE. 
    125125      ! 
    126       llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp   ! test only the inner domain 
    127       IF( COUNT( llmsk(:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    128          zmax(1) = MAXVAL(      vt_i (:,:)      , mask = llmsk )   ! max ice thickness 
    129          zmax(2) = MAXVAL( ABS( u_ice(:,:) )    , mask = llmsk )   ! max ice velocity (zonal only) 
    130          zmax(3) = MAXVAL(     -tm_i (:,:) + rt0, mask = llmsk )   ! min ice temperature (in degC) 
    131       ELSE 
    132          IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible 
    133             zmax(1:3) = -HUGE(1._wp) 
    134          ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 
    135             zmax(1:3) = 0._wp 
    136          ENDIF 
    137       ENDIF 
    138       zmax(4) = REAL( nstop, wp )                                     ! stop indicator 
     126      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp        ! test only the inner domain 
     127      ! 
     128      ll_0oce = .NOT. ANY( llmsk(:,:) )                                         ! no ocean point in the inner domain? 
     129      ! 
     130      zmax(1) = MAXVAL(      vt_i (:,:)      , mask = llmsk )                   ! max ice thickness 
     131      zmax(2) = MAXVAL( ABS( u_ice(:,:) )    , mask = llmsk )                   ! max ice velocity (zonal only) 
     132      zmax(3) = MAXVAL(     -tm_i (:,:) + rt0, mask = llmsk )                   ! min ice temperature (in degC) 
     133      zmax(4) = REAL( nstop, wp )                                               ! stop indicator 
     134      ! 
    139135      !                                   !==               get global extrema             ==! 
    140136      !                                   !==  done by all processes if writting run.stat  ==! 
    141137      IF( ll_colruns ) THEN 
    142138         zmaxlocal(:) = zmax(:) 
    143          CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
     139         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
    144140         nstop = NINT( zmax(4) )                 ! update nstop indicator (now sheared among all local domains) 
    145       ENDIF 
     141      ELSE 
     142         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 
     143         IF( ll_0oce )   zmax(1:3) = 0._wp       ! default "valid" values... 
     144      ENDIF 
     145      ! 
     146      zmax(3) = -zmax(3)                              ! move back from max(-zz) to min(zz) : easier to manage! 
     147      IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3)   ! move back from max(-zz) to min(zz) : easier to manage! 
     148      ! 
    146149      !                                   !==              write "run.stat" files              ==! 
    147150      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    148151      IF( ll_wrtruns ) THEN 
    149          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3) 
    150          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    151          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    152          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
     152         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3) 
     153         DO ji = 1, 3 
     154            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 
     155         END DO 
    153156         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
    154157      END IF 
    … …  
    158161      IF(   zmax(1) >  100._wp .OR.   &                   ! too large ice thickness maximum ( > 100 m) 
    159162         &  zmax(2) >   10._wp .OR.   &                   ! too large ice velocity ( > 10 m/s) 
    160          &  zmax(3) >  101._wp .OR.   &                   ! too cold ice temperature ( < -100 degC) 
     163         &  zmax(3) < -101._wp .OR.   &                   ! too cold ice temperature ( < -100 degC) 
    161164         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
    162165         &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
    … …  
    192195         ! 
    193196         WRITE(ctmp1,*) ' stp_ctl: ice_thick > 100 m or |ice_vel| > 10 m/s or ice_temp < -100 degC or NaN encounter in the tests' 
    194          CALL wrt_line( ctmp2, kt, 'ice_thick max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
    195          CALL wrt_line( ctmp3, kt, '|ice_vel| max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
    196          CALL wrt_line( ctmp4, kt, 'ice_temp  min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
     197         CALL wrt_line( ctmp2, kt, 'ice_thick max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
     198         CALL wrt_line( ctmp3, kt, '|ice_vel| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
     199         CALL wrt_line( ctmp4, kt, 'ice_temp  min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
    197200         IF( Agrif_Root() ) THEN 
    198201            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CANAL/MY_SRC/stpctl.F90

    r13632 r13877  
    6767      REAL(wp)                        ::   zzz                                   ! local real  
    6868      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
    69       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     69      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 
    7070      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
    7171      CHARACTER(len=20)               ::   clname 
    … …  
    125125      ! 
    126126      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain 
     127      ! 
     128      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain? 
     129      ! 
    127130      IF( ll_wd ) THEN 
    128131         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
    … …  
    149152      ENDIF 
    150153      zmax(9) = REAL( nstop, wp )                                                 ! stop indicator 
     154      ! 
    151155      !                                   !==               get global extrema             ==! 
    152156      !                                   !==  done by all processes if writting run.stat  ==! 
    153157      IF( ll_colruns ) THEN 
    154158         zmaxlocal(:) = zmax(:) 
    155          CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
     159         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true.  
    156160         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
    157       ENDIF 
     161      ELSE 
     162         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 
     163         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values... 
     164      ENDIF 
     165      ! 
     166      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage!  
     167      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage! 
     168      IF( ll_colruns ) THEN 
     169         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage!  
     170         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage! 
     171      ENDIF 
     172      ! 
    158173      !                                   !==              write "run.stat" files              ==! 
    159174      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    160175      IF( ll_wrtruns ) THEN 
    161          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
    162          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    163          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    164          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
    165          istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 
    166          istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 
    167          istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 
    168          IF( ln_zad_Aimp ) THEN 
    169             istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 
    170             istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 
    171          ENDIF 
     176         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 
     177         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 
     178            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 
     179         END DO 
    172180         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
    173181      END IF 
    … …  
    175183      !                                   !==  done by all processes at every time step  ==! 
    176184      ! 
    177       IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
    178          &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
    179 !!$         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
    180 !!$         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
    181 !!$         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
    182          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
    183          &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     185      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
     186         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
     187!!$         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
     188!!$         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
     189!!$         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
     190         & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
     191         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
    184192         ! 
    185193         iloc(:,:) = 0 
    … …  
    207215         ELSE                    ! find local min and max locations: 
    208216            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 
    209             llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain 
     217            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain 
    210218            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) ) 
    211219            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
    … …  
    221229         ! 
    222230         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    223          CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
    224          CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
    225          CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
    226          CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
     231         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
     232         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
     233         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
     234         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
    227235         IF( Agrif_Root() ) THEN 
    228236            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CPL_OASIS/EXPREF/file_def_nemo-ice.xml

    r12663 r13877  
    5353       <field field_ref="normstr"          name="normstr" /> 
    5454       <field field_ref="sheastr"          name="sheastr" /> 
    55        <field field_ref="isig1"            name="isig1"   /> 
    56        <field field_ref="isig2"            name="isig2"   /> 
    57        <field field_ref="isig3"            name="isig3"   /> 
     55       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     56       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5857        
    5958       <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_ADV2D/EXPREF/file_def_nemo-ice.xml

    r10516 r13877  
    5555        <field field_ref="normstr"          name="normstr" /> 
    5656        <field field_ref="sheastr"          name="sheastr" /> 
    57         <field field_ref="isig1"            name="isig1"   /> 
    58         <field field_ref="isig2"            name="isig2"   /> 
    59         <field field_ref="isig3"            name="isig3"   /> 
     57       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     58       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    6059 
    6160        <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_AGRIF/EXPREF/file_def_nemo-ice.xml

    r11159 r13877  
    5353        <field field_ref="normstr"          name="normstr" /> 
    5454        <field field_ref="sheastr"          name="sheastr" /> 
    55         <field field_ref="isig1"            name="isig1"   /> 
    56         <field field_ref="isig2"            name="isig2"   /> 
    57         <field field_ref="isig3"            name="isig3"   /> 
     55       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/> 
     56       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/> 
    5857 
    5958        <!-- heat fluxes --> 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/EXPREF/namelist_cfg

    r13632 r13877  
    477477/ 
    478478!----------------------------------------------------------------------- 
    479 !----------------------------------------------------------------------- 
    480 / 
    481 !----------------------------------------------------------------------- 
    482479&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    483480!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/dtatsd.F90

    r13295 r13877  
    191191         ENDIF 
    192192         ! 
    193          DO_2D( 1, 1, 1, 1 ) 
     193         DO_2D( 1, 1, 1, 1 )                  ! vertical interpolation of T & S 
    194194            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
    195195               zl = gdept_0(ji,jj,jk) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/eosbn2.F90

    r13295 r13877  
    182182   !! * Substitutions 
    183183#  include "do_loop_substitute.h90" 
     184#  include "domzgr_substitute.h90" 
    184185   !!---------------------------------------------------------------------- 
    185186   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    965966      IF( ln_timing )   CALL timing_start('bn2') 
    966967      ! 
    967       DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
     968      DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90 
    968969         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
    969970            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
    … …  
    17231724         ! 
    17241725      CASE( np_leos )                        !==  Linear ISOMIP EOS     ==! 
     1726 
     1727         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct) 
     1728          
    17251729         IF(lwp) THEN 
    17261730            WRITE(numout,*) 
    … …  
    17311735            WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0 
    17321736         ENDIF 
     1737         l_useCT = .TRUE.          ! Use conservative temperature 
    17331738         ! 
    17341739      CASE DEFAULT                     !==  ERROR in neos  ==! 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isf_oce.F90

    r12077 r13877  
    7575   ! 
    7676   ! 2.1 -------- ice shelf cavity parameter -------------- 
    77    LOGICAL , PUBLIC            :: l_isfoasis 
     77   LOGICAL , PUBLIC            :: l_isfoasis = .FALSE. 
    7878   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load 
    7979   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfcavgam.F90

    r12905 r13877  
    3030   PUBLIC   isfcav_gammats 
    3131 
     32#  include "domzgr_substitute.h90" 
    3233   !!---------------------------------------------------------------------- 
    3334   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfstp.F90

    r12905 r13877  
    1313   !!   isfstp       : compute iceshelf melt and heat flux 
    1414   !!---------------------------------------------------------------------- 
    15    ! 
    1615   USE isf_oce                                      ! isf variables 
    1716   USE isfload, ONLY: isf_load                      ! ice shelf load 
    … …  
    2120   USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables 
    2221 
    23    USE dom_oce, ONLY: ht, e3t, ln_isfcav, ln_linssh     ! ocean space and time domain 
     22   USE dom_oce        ! ocean space and time domain 
     23   USE oce      , ONLY: ssh                           ! sea surface height 
    2424   USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical 
    2525   USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef. 
    … …  
    3131 
    3232   IMPLICIT NONE 
    33  
    3433   PRIVATE 
    3534 
    3635   PUBLIC   isf_stp, isf_init, isf_nam  ! routine called in sbcmod and divhor 
    3736 
     37   !! * Substitutions 
     38#  include "domzgr_substitute.h90" 
    3839   !!---------------------------------------------------------------------- 
    3940   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    4142   !! Software governed by the CeCILL license (see ./LICENSE) 
    4243   !!---------------------------------------------------------------------- 
     44 
    4345CONTAINS 
    4446  
    … …  
    6062      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    6163      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
     64      !!---------------------------------------------------------------------- 
     65      INTEGER :: jk                               ! loop index 
     66      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t    ! e3t  
    6267      !!--------------------------------------------------------------------- 
    6368      ! 
    … …  
    7883         ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
    7984         rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:) 
    80          CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 
     85         DO jk = 1, jpk 
     86            ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 
     87         END DO  
     88         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 
    8189         ! 
    8290         ! 1.3: compute ice shelf melt 
    … …  
    100108         ! by simplicity, we assume the top level where param applied do not change with time (done in init part) 
    101109         rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:) 
    102          CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 
     110         DO jk = 1, jpk 
     111            ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 
     112         END DO 
     113         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 
    103114         ! 
    104115         ! 2.3: compute ice shelf melt 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/istate.F90

    r13295 r13877  
    2424   USE dom_oce        ! ocean space and time domain  
    2525   USE daymod         ! calendar 
    26    USE divhor         ! horizontal divergence            (div_hor routine) 
    2726   USE dtatsd         ! data temperature and salinity   (dta_tsd routine) 
    2827   USE dtauvd         ! data: U & V current             (dta_uvd routine) 
    … …  
    3534   USE lib_mpp         ! MPP library 
    3635   USE restart         ! restart 
     36#if defined key_agrif 
     37   USE agrif_oce_interp 
     38   USE agrif_oce 
     39#endif    
    3740 
    3841   IMPLICIT NONE 
    … …  
    4346   !! * Substitutions 
    4447#  include "do_loop_substitute.h90" 
     48#  include "domzgr_substitute.h90" 
    4549   !!---------------------------------------------------------------------- 
    4650   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5963      ! 
    6064      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     65      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute 
    6166!!gm see comment further down 
    6267      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace 
    … …  
    7075!!gm  Why not include in the first call of dta_tsd ?   
    7176!!gm  probably associated with the use of internal damping... 
    72                      CALL dta_tsd_init        ! Initialisation of T & S input data 
     77       CALL dta_tsd_init        ! Initialisation of T & S input data 
    7378!!gm to be moved in usrdef of C1D case 
    7479!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data 
    … …  
    8489#endif 
    8590 
     91#if defined key_agrif 
     92      IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN 
     93         numror = 0                           ! define numror = 0 -> no restart file to read 
     94         ln_1st_euler = .true.                ! Set time-step indicator at nit000 (euler forward) 
     95         CALL day_init  
     96         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent 
     97         ! 
     98         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)  
     99         ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
     100         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     101         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     102      ELSE 
     103#endif 
    86104      IF( ln_rstart ) THEN                    ! Restart from a file 
    87105         !                                    ! ------------------- 
    … …  
    100118            ! 
    101119            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest 
     120            uu  (:,:,:,Kbb) = 0._wp 
     121            vv  (:,:,:,Kbb) = 0._wp   
     122            ! 
    102123            IF( ll_wd ) THEN 
    103124               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
    … …  
    111132               END_2D 
    112133            ENDIF  
    113             uu  (:,:,:,Kbb) = 0._wp 
    114             vv  (:,:,:,Kbb) = 0._wp   
    115             ! 
     134             ! 
    116135         ELSE                                 ! user defined initial T and S 
    117             CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
     136            DO jk = 1, jpk 
     137               zgdept(:,:,jk) = gdept(:,:,jk,Kbb) 
     138            END DO 
     139            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    118140         ENDIF 
    119141         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    … …  
    121143         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    122144         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    123          hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
    124          CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value   
    125 !!gm                                    hdiv(:,:,:) = 0._wp 
    126145 
    127146!!gm POTENTIAL BUG : 
    128147!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    129 !!             as well as gdept and gdepw....   !!!!!  
     148!!             as well as gdept_ and gdepw_....   !!!!!  
    130149!!      ===>>>>   probably a call to domvvl initialisation here.... 
    131150 
    … …  
    151170         !  
    152171      ENDIF  
     172#if defined key_agrif 
     173      ENDIF 
     174#endif 
    153175      !  
    154176      ! Initialize "now" and "before" barotropic velocities: 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/sbcfwb.F90

    r13632 r13877  
    1717   USE dom_oce        ! ocean space and time domain 
    1818   USE sbc_oce        ! surface ocean boundary condition 
    19    USE isf_oce       ! ice shelf melting contribution 
     19   USE isf_oce , ONLY : fwfisf_cav, fwfisf_par, ln_isfcpl, ln_isfcpl_cons, risfcpl_cons_ssh ! ice shelf melting contribution 
    2020   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass 
    2121   USE phycst         ! physical constants 
    … …  
    7171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      - 
    7272      REAL(wp)   ,DIMENSION(1) ::   z_fwfprv   
    73       COMPLEX(wp),DIMENSION(1) ::   y_fwfnow   
     73      COMPLEX(dp),DIMENSION(1) ::   y_fwfnow   
    7474      !!---------------------------------------------------------------------- 
    7575      ! 
    … …  
    207207!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain ! 
    208208#if defined key_mpi3 
    209             CALL lbc_lnk_nc_multi( 'sbcfwb', zerp_cor, 'T', 1. ) 
     209            CALL lbc_lnk_nc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp ) 
    210210#else 
    211             CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1. ) 
     211            CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp ) 
    212212#endif 
    213213            ! 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP/EXPREF/namelist_cfg

    r13632 r13877  
    437437/ 
    438438!----------------------------------------------------------------------- 
    439 !----------------------------------------------------------------------- 
    440 / 
    441 !----------------------------------------------------------------------- 
    442439&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    443440!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/STATION_ASF/MY_SRC/stpctl.F90

    r13632 r13877  
    6565      REAL(wp)                        ::   zzz                                   ! local real  
    6666      REAL(wp), DIMENSION(4)          ::   zmax, zmaxlocal 
    67       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     67      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 
    6868      LOGICAL, DIMENSION(jpi,jpj)     ::   llmsk 
    6969      CHARACTER(len=20)               ::   clname 
    … …  
    115115      llmsk(:,Nje1: jpj) = .FALSE. 
    116116      ! 
    117       llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp   ! test only the inner domain 
    118       IF( COUNT( llmsk(:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    119          zmax(1) = MAXVAL(     taum(:,:)   , mask = llmsk )   ! max wind stress module 
    120          zmax(2) = MAXVAL( ABS( qns(:,:) ) , mask = llmsk )   ! max non-solar heat flux 
    121          zmax(3) = MAXVAL( ABS( emp(:,:) ) , mask = llmsk )   ! max E-P 
    122       ELSE 
    123          IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible 
    124             zmax(1:3) = -HUGE(1._wp) 
    125          ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 
    126             zmax(1:3) = 0._wp 
    127          ENDIF 
    128       ENDIF 
    129       zmax(4) = REAL( nstop, wp )                                     ! stop indicator 
     117      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp        ! test only the inner domain 
     118      ! 
     119      ll_0oce = .NOT. ANY( llmsk(:,:) )                                         ! no ocean point in the inner domain? 
     120      ! 
     121      zmax(1) = MAXVAL(     taum(:,:)  , mask = llmsk )                         ! max wind stress module 
     122      zmax(2) = MAXVAL( ABS( qns(:,:) ), mask = llmsk )                         ! max non-solar heat flux 
     123      zmax(3) = MAXVAL( ABS( emp(:,:) ), mask = llmsk )                         ! max E-P 
     124      zmax(4) = REAL( nstop, wp )                                               ! stop indicator 
     125      ! 
    130126      !                                   !==               get global extrema             ==! 
    131127      !                                   !==  done by all processes if writting run.stat  ==! 
    … …  
    134130         CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    135131         nstop = NINT( zmax(4) )                 ! update nstop indicator (now sheared among all local domains) 
    136       ENDIF 
     132      ELSE 
     133         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 
     134         IF( ll_0oce )   zmax(1:3) = 0._wp       ! default "valid" values... 
     135      ENDIF 
     136      !                                   !==               error handling               ==! 
    137137      !                                   !==              write "run.stat" files              ==! 
    138138      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    139139      IF( ll_wrtruns ) THEN 
    140140         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3) 
    141          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    142          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    143          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/ zmax(3)/), (/kt/), (/1/) ) 
     141         DO ji = 1, 3 
     142            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 
     143         END DO 
    144144         IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 
    145145      END IF 
Note: See TracChangeset for help on using the changeset viewer.