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

Changeset 13886


Ignore:
Timestamp:
2020-11-26T15:24:38+01:00 (4 years ago)
Author:
rlod
Message:

phasing with trunk at revision r13787

Location:
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
Files:
117 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source

    • Property svn:externals

      •  

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

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/INSTALL.rst

    r11734 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r13472 r13886  
    358358!!                                                                    !! 
    359359!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    360 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    361360!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    362361!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r13216 r13886  
    307307!!                                                                    !! 
    308308!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    309 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    310309!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    311310!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

    r13216 r13886  
    307307!!                                                                    !! 
    308308!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    309 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    310309!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    311310!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml

    r9896 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r13472 r13886  
    359359!!                                                                    !! 
    360360!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    361 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    362361!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    363362!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/AMM12/EXPREF/namelist_cfg

    r13461 r13886  
    345345!!                                                                    !! 
    346346!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    347 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    348347!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    349348!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/C1D_PAPA/EXPREF/namelist_cfg

    r13461 r13886  
    414414!!                                                                    !! 
    415415!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    416 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    417416!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    418417!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    … …  
    429428/ 
    430429!----------------------------------------------------------------------- 
    431 &namptr        !   Poleward Transport Diagnostic                        (default: OFF) 
    432 !----------------------------------------------------------------------- 
    433 / 
    434 !----------------------------------------------------------------------- 
    435430&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    436431!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/GYRE_BFM/EXPREF/namelist_cfg

    r13461 r13886  
    223223!!                                                                    !! 
    224224!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    225 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    226225!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    227226!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/GYRE_PISCES/EXPREF/namelist_cfg

    r13461 r13886  
    217217!!                                                                    !! 
    218218!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    219 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    220219!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    221220!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml

    r11945 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg

    r13461 r13886  
    401401!!                                                                    !! 
    402402!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    403 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    404403!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    405404!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml

    r12377 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r13472 r13886  
    400400!!                                                                    !! 
    401401!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    402 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    403402!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    404403!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg

    r13461 r13886  
    199199&namdrg        !   top/bottom drag coefficient                          (default: NO selection) 
    200200!----------------------------------------------------------------------- 
     201   ln_drg_OFF  = .true.   !  must select one drag coefficient to true even if we don't use it 
    201202/ 
    202203!----------------------------------------------------------------------- 
    … …  
    347348&namzdf        !   vertical physics manager                             (default: NO selection) 
    348349!----------------------------------------------------------------------- 
     350   ln_zdfcst   = .true.      !  must select one vertical physics to true even if we don't use it 
    349351/ 
    350352!----------------------------------------------------------------------- 
    … …  
    372374!!                                                                    !! 
    373375!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    374 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    375376!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    376377!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    … …  
    387388/ 
    388389!----------------------------------------------------------------------- 
    389 &namptr        !   Poleward Transport Diagnostic                        (default: OFF) 
    390 !----------------------------------------------------------------------- 
    391 / 
    392 !----------------------------------------------------------------------- 
    393390&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    394391!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg

    r13461 r13886  
    370370!!                                                                    !! 
    371371!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    372 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    373372!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    374373!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    … …  
    385384/ 
    386385!----------------------------------------------------------------------- 
    387 &namptr        !   Poleward Transport Diagnostic                        (default: OFF) 
    388 !----------------------------------------------------------------------- 
    389 / 
    390 !----------------------------------------------------------------------- 
    391386&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    392387!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml

    r9896 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg

    r13461 r13886  
    171171!!                                                                    !! 
    172172!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    173 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    174173!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    175174!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/field_def_nemo-ice.xml

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/namelist_ref

    r13812 r13886  
    459459!----------------------------------------------------------------------- 
    460460   ln_rnf_mouth = .false.   !  specific treatment at rivers mouths 
    461       rn_hrnf     =  15.e0    !  depth over which enhanced vertical mixing is used    (ln_rnf_mouth=T) 
    462       rn_avt_rnf  =   1.e-3   !  value of the additional vertical mixing coef. [m2/s] (ln_rnf_mouth=T) 
    463    rn_rfact    =   1.e0    !  multiplicative factor for runoff 
     461      rn_hrnf     =  15.e0     !  depth over which enhanced vertical mixing is used    (ln_rnf_mouth=T) 
     462      rn_avt_rnf  =   1.e-3    !  value of the additional vertical mixing coef. [m2/s] (ln_rnf_mouth=T) 
     463   rn_rfact     =   1.e0    !  multiplicative factor for runoff 
    464464   ln_rnf_depth = .false.   !  read in depth information for runoff 
    465    ln_rnf_tem  = .false.   !  read in temperature information for runoff 
    466    ln_rnf_sal  = .false.   !  read in salinity information for runoff 
    467    ln_rnf_depth_ini = .false. ! compute depth at initialisation from runoff file 
    468       rn_rnf_max  = 5.735e-4  !  max value of the runoff climatologie over global domain ( ln_rnf_depth_ini = .true ) 
    469       rn_dep_max  = 150.      !  depth over which runoffs is spread ( ln_rnf_depth_ini = .true ) 
    470       nn_rnf_depth_file = 0   !  create (=1) a runoff depth file or not (=0) 
    471    ln_rnf_icb  = .false.    ! freshwater flux from icebergs in runoff file 
    472  
    473    cn_dir      = './'      !  root directory for the runoff data location 
     465   ln_rnf_tem   = .false.   !  read in temperature information for runoff 
     466   ln_rnf_sal   = .false.   !  read in salinity information for runoff 
     467   ln_rnf_icb   = .false.   !  read iceberg flux 
     468   ln_rnf_depth_ini = .false.  ! compute depth at initialisation from runoff file 
     469      rn_rnf_max  = 5.735e-4   !  max value of the runoff climatologie over global domain ( ln_rnf_depth_ini = .true ) 
     470      rn_dep_max  = 150.       !  depth over which runoffs is spread ( ln_rnf_depth_ini = .true ) 
     471      nn_rnf_depth_file = 0    !  create (=1) a runoff depth file or not (=0) 
     472 
     473   cn_dir       = './'      !  root directory for the runoff data location 
    474474   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
    475475   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    476476   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    477477   sn_rnf      = 'runoff_core_monthly'   ,        -1.        , 'sorunoff',   .true.    , .true. , 'yearly'  , ''               , ''       , '' 
    478    sn_cnf      = 'runoff_core_monthly'   ,         0.        , 'socoefr0',   .false.   , .true. , 'yearly'  , ''               , ''       , '' 
     478   sn_cnf      = 'runoff_core_monthly'   ,       -12.        , 'socoefr0',   .false.   , .true. , 'yearly'  , ''               , ''       , '' 
    479479   sn_s_rnf    = 'runoffs'               ,        24.        , 'rosaline',   .true.    , .true. , 'yearly'  , ''               , ''       , '' 
    480480   sn_t_rnf    = 'runoffs'               ,        24.        , 'rotemper',   .true.    , .true. , 'yearly'  , ''               , ''       , '' 
    481    sn_dep_rnf  = 'runoffs'               ,         0.        , 'rodepth' ,   .false.   , .true. , 'yearly'  , ''               , ''       , '' 
     481   sn_i_rnf    = 'NOT USED'              ,        24.        , 'xxxxxxxx',   .true.    , .true. , 'yearly'  , ''               , ''       , '' 
     482   sn_dep_rnf  = 'runoffs'               ,       -12.        , 'rodepth' ,   .false.   , .true. , 'yearly'  , ''               , ''       , '' 
    482483/ 
    483484!----------------------------------------------------------------------- 
    … …  
    10501051   ln_dynrnf       =  .false.    !  runoffs option enabled (T) or not (F) 
    10511052   ln_dynrnf_depth =  .false.    !  runoffs is spread in vertical (T) or not (F) 
    1052 !   fwbcorr        = 3.786e-06   !  annual global mean of empmr for ssh correction 
     1053   fwbcorr         =    0.0      !  annual global mean of empmr for ssh correction 
    10531054 
    10541055   cn_dir      = './'      !  root directory for the ocean data location 
    … …  
    12341235!!                                                                    !! 
    12351236!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    1236 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    12371237!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    12381238!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    … …  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml

    r11536 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SPITZ12/EXPREF/namelist_cfg

    r13472 r13886  
    359359!!                                                                    !! 
    360360!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    361 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    362361!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    363362!!   namdiu       Cool skin and warm layer models                       (default: OFF) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/WED025/EXPREF/file_def_nemo-ice.xml

    r12905 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/WED025/EXPREF/namelist_cfg

    r13461 r13886  
    567567!!                                                                    !! 
    568568!!   namtrd       dynamics and/or tracer trends                         (default: OFF) 
    569 !!   namptr       Poleward Transport Diagnostics                        (default: OFF) 
    570569!!   namhsb       Heat and salt budgets                                 (default: OFF) 
    571570!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    … …  
    584583/ 
    585584!----------------------------------------------------------------------- 
    586 &namptr        !   Poleward Transport Diagnostic                        (default: OFF) 
    587 !----------------------------------------------------------------------- 
    588 / 
    589 !----------------------------------------------------------------------- 
    590585&namhsb        !  Heat and salt budgets                                 (default: OFF) 
    591586!----------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/ice.F90

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/ice1d.F90

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icecor.F90

    r13497 r13886  
    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 ! 
    … …  
    115118         CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
    116119      ENDIF 
    117  
    118       !                             !----------------------------------------------------- 
    119       SELECT CASE( kn )             !  Diagnostics                                       ! 
    120       !                             !----------------------------------------------------- 
    121       CASE( 1 )                        !--- dyn trend diagnostics 
    122          ! 
    123          IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
    124             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 
    125                &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
    126             diag_sice(:,:) =   SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    127             diag_vice(:,:) =   SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    128             diag_vsnw(:,:) =   SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
    129          ENDIF 
    130          !                       ! concentration tendency (dynamics) 
    131          IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
    132             zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice  
    133             CALL iom_put( 'afxdyn' , zafx ) 
    134          ENDIF 
    135          ! 
    136       CASE( 2 )                        !--- thermo trend diagnostics & ice aging 
    137          ! 
    138          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice   ! ice natural aging incrementation 
    139          ! 
    140          IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 
    141             diag_heat(:,:) = diag_heat(:,:) & 
    142                &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 
    143                &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
    144             diag_sice(:,:) = diag_sice(:,:) & 
    145                &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    146             diag_vice(:,:) = diag_vice(:,:) & 
    147                &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi 
    148             diag_vsnw(:,:) = diag_vsnw(:,:) & 
    149                &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
    150             CALL iom_put ( 'hfxdhc' , diag_heat )  
    151          ENDIF 
    152          !                       ! concentration tendency (total + thermo) 
    153          IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN  
    154             zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
    155             CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) 
    156             CALL iom_put( 'afxtot' , zafx ) 
    157          ENDIF 
    158          ! 
    159       END SELECT 
    160120      ! 
    161121      ! controls 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icectl.F90

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_adv_pra.F90

    r13497 r13886  
    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      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    132118      ! 
    … …  
    141127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
    142128         END WHERE 
    143       END DO 
    144       DO jl = 1, jpl 
    145          DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
    146             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), & 
    147                &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
    148                &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
    149                &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
    150          END_3D 
    151       END DO 
    152       DO jl = 1, jpl 
    153          DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
    154             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), & 
    155                &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
    156                &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
    157                &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
    158          END_3D 
    159       END DO 
    160       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    161       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
     132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
    162134      ! 
    163135      ! 
    … …  
    175147      ENDIF 
    176148      zdt = rDt_ice / REAL(icycle) 
     149      z1_dt = 1._wp / zdt 
    177150       
    178151      ! --- transport --- ! 
    … …  
    181154 
    182155      DO jt = 1, icycle 
     156 
     157         ! diagnostics 
     158         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     159         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     160         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     161            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    183162 
    184163         ! record at_i before advection (for open water) 
    … …  
    279258                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
    280259               ENDIF 
    281            ENDIF 
    282             ! 
     260            ENDIF 
     261            ! 
     262         ENDIF 
     263          
     264         ! --- Lateral boundary conditions --- ! 
     265         !     caution: for gradients (sx and sy) the sign changes 
     266         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume 
     267            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  & 
     268            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume 
     269            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  ) 
     270         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity 
     271            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  & 
     272            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration 
     273            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  ) 
     274         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age 
     275            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  ) 
     276         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy 
     277            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  )  
     278         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
     279            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
     280         IF ( ln_pnd_LEV ) THEN 
     281            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
     282               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
     283               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume 
     284               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )  
     285            IF ( ln_pnd_lids ) THEN 
     286               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume 
     287                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )  
     288            ENDIF 
    283289         ENDIF 
    284290 
    … …  
    312318         END_2D 
    313319         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp ) 
     320         ! 
     321         ! --- diagnostics --- ! 
     322         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     323            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     324         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     325            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     326         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     327            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     328            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    314329         ! 
    315330         ! --- Ensure non-negative fields --- ! 
    … …  
    350365      !!  
    351366      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     367      INTEGER  ::   jj0                                  ! dummy loop indices 
    352368      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    353369      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    354370      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     371      REAL(wp) ::   zpsm, zps0 
     372      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    355373      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    356374      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    357375      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    358376      !----------------------------------------------------------------------- 
     377      ! in order to avoid lbc_lnk (communications): 
     378      !    jj loop must be 1:jpj   if adv_x is called first 
     379      !                and 2:jpj-1 if adv_x is called second 
     380      jj0 = NINT(pcrh) 
    359381      ! 
    360382      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    363385         ! 
    364386         ! Limitation of moments.                                            
    365          DO_2D( 0, 0, 1, 1 ) 
    366             !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    367             psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    368             ! 
    369             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    370             zs1max  = 1.5 * zslpmax 
    371             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    372             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    373                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    374             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    375  
    376             ps0 (ji,jj,jl) = zslpmax   
    377             psx (ji,jj,jl) = zs1new         * rswitch 
    378             psxx(ji,jj,jl) = zs2new         * rswitch 
    379             psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    380             psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    381             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    382          END_2D 
    383  
    384          !  Calculate fluxes and moments between boxes i<-->i+1               
    385          DO_2D( 0, 0, 1, 1 )                   !  Flux from i to i+1 WHEN u GT 0 
    386             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    387             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    388             zalfq        =  zalf * zalf 
    389             zalf1        =  1.0 - zalf 
    390             zalf1q       =  zalf1 * zalf1 
    391             ! 
    392             zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    393             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    394             zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    395             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    396             zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    397             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    398             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    399  
    400             !  Readjust moments remaining in the box. 
    401             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    402             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    403             psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    404             psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    405             psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    406             psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    407             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    408          END_2D 
    409  
    410          DO_2D( 0, 0, 1, 0 )                   !  Flux from i+1 to i when u LT 0. 
    411             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    412             zalg  (ji,jj) = zalf 
    413             zalfq         = zalf * zalf 
    414             zalf1         = 1.0 - zalf 
    415             zalg1 (ji,jj) = zalf1 
    416             zalf1q        = zalf1 * zalf1 
    417             zalg1q(ji,jj) = zalf1q 
    418             ! 
    419             zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
    420             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
    421                &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
    422             zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
    423             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
    424             zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
    425             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
    426             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    427          END_2D 
    428  
    429          DO_2D( 0, 0, 0, 0 )                   !  Readjust moments remaining in the box. 
    430             zbt  =       zbet(ji-1,jj) 
    431             zbt1 = 1.0 - zbet(ji-1,jj) 
    432             ! 
    433             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    434             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    435             psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    436             psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    437             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    438             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    439             psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    440          END_2D 
    441  
    442          !   Put the temporary moments into appropriate neighboring boxes.     
    443          DO_2D( 0, 0, 0, 0 )                   !   Flux from i to i+1 IF u GT 0. 
    444             zbt  =       zbet(ji-1,jj) 
    445             zbt1 = 1.0 - zbet(ji-1,jj) 
    446             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    447             zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    448             zalf1         = 1.0 - zalf 
    449             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    450             ! 
    451             ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    452             psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    453             psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    454                &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
    455                &            + zbt1 * psxx(ji,jj,jl) 
    456             psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    457                &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    458                &            + zbt1 * psxy(ji,jj,jl) 
    459             psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    460             psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    461          END_2D 
    462  
    463          DO_2D( 0, 0, 0, 0 )                   !  Flux from i+1 to i IF u LT 0. 
    464             zbt  =       zbet(ji,jj) 
    465             zbt1 = 1.0 - zbet(ji,jj) 
    466             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    467             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    468             zalf1         = 1.0 - zalf 
    469             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    470             ! 
    471             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    472             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    473             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    474                &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    475                &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    476             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    477                &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    478             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    479             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
    480          END_2D 
    481  
     387         DO jj = Njs0 - jj0, Nje0 + jj0 
     388             
     389            DO ji = Nis0 - 1, Nie0 + 1 
     390 
     391               zpsm  = psm (ji,jj,jl) ! optimization 
     392               zps0  = ps0 (ji,jj,jl) 
     393               zpsx  = psx (ji,jj,jl) 
     394               zpsxx = psxx(ji,jj,jl) 
     395               zpsy  = psy (ji,jj,jl) 
     396               zpsyy = psyy(ji,jj,jl) 
     397               zpsxy = psxy(ji,jj,jl) 
     398 
     399               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     400               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     401               ! 
     402               zslpmax = MAX( 0._wp, zps0 ) 
     403               zs1max  = 1.5 * zslpmax 
     404               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     405               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
     406               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     407 
     408               zps0  = zslpmax   
     409               zpsx  = zs1new  * rswitch 
     410               zpsxx = zs2new  * rswitch 
     411               zpsy  = zpsy    * rswitch 
     412               zpsyy = zpsyy   * rswitch 
     413               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     414 
     415               !  Calculate fluxes and moments between boxes i<-->i+1               
     416               !                                !  Flux from i to i+1 WHEN u GT 0  
     417               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     418               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
     419               zalfq        =  zalf * zalf 
     420               zalf1        =  1.0 - zalf 
     421               zalf1q       =  zalf1 * zalf1 
     422               ! 
     423               zfm (ji,jj)  =  zalf  *   zpsm  
     424               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     425               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     426               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     427               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     428               zfxy(ji,jj)  =  zalfq *   zpsxy 
     429               zfyy(ji,jj)  =  zalf  *   zpsyy 
     430 
     431               !                                !  Readjust moments remaining in the box. 
     432               zpsm  =  zpsm  - zfm(ji,jj) 
     433               zps0  =  zps0  - zf0(ji,jj) 
     434               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     435               zpsxx =  zalf1  * zalf1q * zpsxx 
     436               zpsy  =  zpsy  - zfy (ji,jj) 
     437               zpsyy =  zpsyy - zfyy(ji,jj) 
     438               zpsxy =  zalf1q * zpsxy 
     439               ! 
     440               psm (ji,jj,jl) = zpsm ! optimization 
     441               ps0 (ji,jj,jl) = zps0  
     442               psx (ji,jj,jl) = zpsx  
     443               psxx(ji,jj,jl) = zpsxx 
     444               psy (ji,jj,jl) = zpsy  
     445               psyy(ji,jj,jl) = zpsyy 
     446               psxy(ji,jj,jl) = zpsxy 
     447               ! 
     448            END DO 
     449             
     450            DO ji = Nis0 - 1, Nie0 
     451               !                                !  Flux from i+1 to i when u LT 0. 
     452               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     453               zalg  (ji,jj) = zalf 
     454               zalfq         = zalf * zalf 
     455               zalf1         = 1.0 - zalf 
     456               zalg1 (ji,jj) = zalf1 
     457               zalf1q        = zalf1 * zalf1 
     458               zalg1q(ji,jj) = zalf1q 
     459               ! 
     460               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     461               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     462                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     463               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     464               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     465               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     466               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     467               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     468            END DO 
     469 
     470            DO ji = Nis0, Nie0 
     471               ! 
     472               zpsm  = psm (ji,jj,jl) ! optimization 
     473               zps0  = ps0 (ji,jj,jl) 
     474               zpsx  = psx (ji,jj,jl) 
     475               zpsxx = psxx(ji,jj,jl) 
     476               zpsy  = psy (ji,jj,jl) 
     477               zpsyy = psyy(ji,jj,jl) 
     478               zpsxy = psxy(ji,jj,jl) 
     479               !                                !  Readjust moments remaining in the box. 
     480               zbt  =       zbet(ji-1,jj) 
     481               zbt1 = 1.0 - zbet(ji-1,jj) 
     482               ! 
     483               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     484               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     485               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     486               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     487               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) ) 
     488               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     489               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     490 
     491               !   Put the temporary moments into appropriate neighboring boxes.     
     492               !                                !   Flux from i to i+1 IF u GT 0. 
     493               zbt   =       zbet(ji-1,jj) 
     494               zbt1  = 1.0 - zbet(ji-1,jj) 
     495               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     496               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     497               zalf1 = 1.0 - zalf 
     498               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     499               ! 
     500               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     501               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     502               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     503                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     504                  &            + zbt1 * zpsxx 
     505               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     506                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     507                  &            + zbt1 * zpsxy 
     508               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     509               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     510 
     511               !                                !  Flux from i+1 to i IF u LT 0. 
     512               zbt   =       zbet(ji,jj) 
     513               zbt1  = 1.0 - zbet(ji,jj) 
     514               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     515               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     516               zalf1 = 1.0 - zalf 
     517               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     518               ! 
     519               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     520               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     521               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     522                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     523                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     524               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     525                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     526               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     527               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     528               ! 
     529               psm (ji,jj,jl) = zpsm  ! optimization 
     530               ps0 (ji,jj,jl) = zps0  
     531               psx (ji,jj,jl) = zpsx  
     532               psxx(ji,jj,jl) = zpsxx 
     533               psy (ji,jj,jl) = zpsy  
     534               psyy(ji,jj,jl) = zpsyy 
     535               psxy(ji,jj,jl) = zpsxy 
     536            END DO 
     537            ! 
     538         END DO 
     539         ! 
    482540      END DO 
    483  
    484       !-- Lateral boundary conditions 
    485       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    486          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    487          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    488       ! 
     541      !       
    489542   END SUBROUTINE adv_x 
    490543 
    … …  
    507560      !! 
    508561      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     562      INTEGER  ::   ji0                                  ! dummy loop indices 
    509563      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    510564      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    511565      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     566      REAL(wp) ::   zpsm, zps0 
     567      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    512568      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    513569      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    514570      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    515571      !--------------------------------------------------------------------- 
     572      ! in order to avoid lbc_lnk (communications): 
     573      !    ji loop must be 1:jpi   if adv_y is called first 
     574      !                and 2:jpi-1 if adv_y is called second 
     575      ji0 = NINT(pcrh) 
    516576      ! 
    517577      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    520580         ! 
    521581         ! Limitation of moments. 
    522          DO_2D( 1, 1, 0, 0 ) 
    523             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    524             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    525             ! 
    526             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     582         DO_2D( 1, 1, ji0, ji0 ) 
     583            ! 
     584            zpsm  = psm (ji,jj,jl) ! optimization 
     585            zps0  = ps0 (ji,jj,jl) 
     586            zpsx  = psx (ji,jj,jl) 
     587            zpsxx = psxx(ji,jj,jl) 
     588            zpsy  = psy (ji,jj,jl) 
     589            zpsyy = psyy(ji,jj,jl) 
     590            zpsxy = psxy(ji,jj,jl) 
     591            ! 
     592            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     593            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     594            ! 
     595            zslpmax = MAX( 0._wp, zps0 ) 
    527596            zs1max  = 1.5 * zslpmax 
    528             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    529             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    530                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     597            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     598            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    531599            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    532600            ! 
    533             ps0 (ji,jj,jl) = zslpmax   
    534             psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    535             psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    536             psy (ji,jj,jl) = zs1new         * rswitch 
    537             psyy(ji,jj,jl) = zs2new         * rswitch 
    538             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    539          END_2D 
    540   
    541          !  Calculate fluxes and moments between boxes j<-->j+1               
    542          DO_2D( 1, 1, 0, 0 )                   !  Flux from j to j+1 WHEN v GT 0 
     601            zps0  = zslpmax   
     602            zpsx  = zpsx  * rswitch 
     603            zpsxx = zpsxx * rswitch 
     604            zpsy  = zs1new         * rswitch 
     605            zpsyy = zs2new         * rswitch 
     606            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     607 
     608            !  Calculate fluxes and moments between boxes j<-->j+1               
     609            !                                !  Flux from j to j+1 WHEN v GT 0    
    543610            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    544             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     611            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    545612            zalfq        =  zalf * zalf 
    546613            zalf1        =  1.0 - zalf 
    547614            zalf1q       =  zalf1 * zalf1 
    548615            ! 
    549             zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    550             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    551             zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    552             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    553             zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    554             zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    555             zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
    556             ! 
    557             !  Readjust moments remaining in the box. 
    558             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    559             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    560             psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    561             psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    562             psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    563             psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    564             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     616            zfm (ji,jj)  =  zalf  * zpsm 
     617            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     618            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     619            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     620            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     621            zfxy(ji,jj)  =  zalfq * zpsxy 
     622            zfxx(ji,jj)  =  zalf  * zpsxx 
     623            ! 
     624            !                                !  Readjust moments remaining in the box. 
     625            zpsm   =  zpsm  - zfm(ji,jj) 
     626            zps0   =  zps0  - zf0(ji,jj) 
     627            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     628            zpsyy  =  zalf1 * zalf1q * zpsyy 
     629            zpsx   =  zpsx  - zfx(ji,jj) 
     630            zpsxx  =  zpsxx - zfxx(ji,jj) 
     631            zpsxy  =  zalf1q * zpsxy 
     632            ! 
     633            psm (ji,jj,jl) = zpsm ! optimization 
     634            ps0 (ji,jj,jl) = zps0  
     635            psx (ji,jj,jl) = zpsx  
     636            psxx(ji,jj,jl) = zpsxx 
     637            psy (ji,jj,jl) = zpsy  
     638            psyy(ji,jj,jl) = zpsyy 
     639            psxy(ji,jj,jl) = zpsxy 
    565640         END_2D 
    566641         ! 
    567          DO_2D( 1, 0, 0, 0 )                   !  Flux from j+1 to j when v LT 0. 
     642         DO_2D( 1, 0, ji0, ji0 ) 
     643            !                                !  Flux from j+1 to j when v LT 0. 
    568644            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    569645            zalg  (ji,jj) = zalf 
    … …  
    584660         END_2D 
    585661 
    586          !  Readjust moments remaining in the box.  
    587          DO_2D( 0, 0, 0, 0 ) 
     662         DO_2D( 0, 0, ji0, ji0 ) 
     663            !                                !  Readjust moments remaining in the box. 
    588664            zbt  =         zbet(ji,jj-1) 
    589665            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    590666            ! 
    591             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    592             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    593             psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    594             psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    595             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    596             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    597             psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     667            zpsm  = psm (ji,jj,jl) ! optimization 
     668            zps0  = ps0 (ji,jj,jl) 
     669            zpsx  = psx (ji,jj,jl) 
     670            zpsxx = psxx(ji,jj,jl) 
     671            zpsy  = psy (ji,jj,jl) 
     672            zpsyy = psyy(ji,jj,jl) 
     673            zpsxy = psxy(ji,jj,jl) 
     674            ! 
     675            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     676            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     677            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     678            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     679            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     680            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     681            zpsxy = zalg1q(ji,jj-1) * zpsxy 
     682 
     683            !   Put the temporary moments into appropriate neighboring boxes.     
     684            !                                !   Flux from j to j+1 IF v GT 0. 
     685            zbt   =       zbet(ji,jj-1) 
     686            zbt1  = 1.0 - zbet(ji,jj-1) 
     687            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     688            zalf  = zbt * zfm(ji,jj-1) / zpsm  
     689            zalf1 = 1.0 - zalf 
     690            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     691            ! 
     692            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     693            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     694               &             + zbt1 * zpsy   
     695            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     696               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     697               &             + zbt1 * zpsyy 
     698            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     699               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     700               &             + zbt1 * zpsxy 
     701            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     702            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
     703 
     704            !                                !  Flux from j+1 to j IF v LT 0. 
     705            zbt   =       zbet(ji,jj) 
     706            zbt1  = 1.0 - zbet(ji,jj) 
     707            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     708            zalf  = zbt1 * zfm(ji,jj) / zpsm 
     709            zalf1 = 1.0 - zalf 
     710            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     711            ! 
     712            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     713            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     714            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     715               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     716               &                         + ( zalf1 - zalf ) * ztemp ) ) 
     717            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     718               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     719            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     720            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     721            ! 
     722            psm (ji,jj,jl) = zpsm ! optimization 
     723            ps0 (ji,jj,jl) = zps0  
     724            psx (ji,jj,jl) = zpsx  
     725            psxx(ji,jj,jl) = zpsxx 
     726            psy (ji,jj,jl) = zpsy  
     727            psyy(ji,jj,jl) = zpsyy 
     728            psxy(ji,jj,jl) = zpsxy 
    598729         END_2D 
    599  
    600          !   Put the temporary moments into appropriate neighboring boxes.     
    601          DO_2D( 0, 0, 0, 0 )                   !  Flux from j to j+1 IF v GT 0. 
    602             zbt  =       zbet(ji,jj-1) 
    603             zbt1 = 1.0 - zbet(ji,jj-1) 
    604             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    605             zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    606             zalf1         = 1.0 - zalf 
    607             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    608             ! 
    609             ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    610             psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    611                &             + zbt1 * psy(ji,jj,jl)   
    612             psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    613                &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    614                &             + zbt1 * psyy(ji,jj,jl) 
    615             psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    616                &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    617                &             + zbt1 * psxy(ji,jj,jl) 
    618             psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    619             psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
    620          END_2D 
    621  
    622          DO_2D( 0, 0, 0, 0 )                   !  Flux from j+1 to j IF v LT 0. 
    623             zbt  =       zbet(ji,jj) 
    624             zbt1 = 1.0 - zbet(ji,jj) 
    625             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    626             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    627             zalf1         = 1.0 - zalf 
    628             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    629             ! 
    630             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    631             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    632             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    633                &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    634                &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    635             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    636                &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    637             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    638             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
    639          END_2D 
    640  
     730         ! 
    641731      END DO 
    642  
    643       !-- Lateral boundary conditions 
    644       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    645          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    646          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    647732      ! 
    648733   END SUBROUTINE adv_y 
    … …  
    872957            ! 
    873958            !                                                        ! ice thickness 
    874             CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice ) 
    875             CALL iom_get( numrir, jpdom_auto, 'syice' , syice ) 
     959            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 
     960            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 
    876961            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 
    877962            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 
    878963            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 
    879964            !                                                        ! snow thickness 
    880             CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  ) 
    881             CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  ) 
     965            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp ) 
     966            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp ) 
    882967            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  ) 
    883968            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  ) 
    884969            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  ) 
    885970            !                                                        ! ice concentration 
    886             CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    ) 
    887             CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    ) 
     971            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp ) 
     972            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp ) 
    888973            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   ) 
    889974            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   ) 
    890975            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   ) 
    891976            !                                                        ! ice salinity 
    892             CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal ) 
    893             CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal ) 
     977            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 
     978            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 
    894979            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 
    895980            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 
    896981            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 
    897982            !                                                        ! ice age 
    898             CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage ) 
    899             CALL iom_get( numrir, jpdom_auto, 'syage' , syage ) 
     983            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 
     984            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 
    900985            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 
    901986            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 
    … …  
    904989            DO jk = 1, nlay_s 
    905990               WRITE(zchar1,'(I2.2)') jk 
    906                znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
    907                znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
     991               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
     992               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
    908993               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
    909994               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
    … …  
    913998            DO jk = 1, nlay_i 
    914999               WRITE(zchar1,'(I2.2)') jk 
    915                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
    916                znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
     1000               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
     1001               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
    9171002               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
    9181003               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
    … …  
    9211006            ! 
    9221007            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
    923                IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
    924                   CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap ) 
    925                   CALL iom_get( numrir, jpdom_auto, 'syap' , syap ) 
     1008               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     1009                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 
     1010                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 
    9261011                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
    9271012                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
    9281013                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
    9291014                  !                                                     ! melt pond volume 
    930                   CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp ) 
    931                   CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp ) 
     1015                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 
     1016                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 
    9321017                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
    9331018                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
    … …  
    9391024                  ! 
    9401025               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
    941                   IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
    942                      CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl ) 
    943                      CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl ) 
     1026                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     1027                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 
     1028                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 
    9441029                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 
    9451030                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 
    … …  
    10561141   END SUBROUTINE adv_pra_rst 
    10571142 
     1143   SUBROUTINE icemax3D( pice , pmax ) 
     1144      !!--------------------------------------------------------------------- 
     1145      !!                   ***  ROUTINE icemax3D ***                      
     1146      !! ** Purpose :  compute the max of the 9 points around 
     1147      !!---------------------------------------------------------------------- 
     1148      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1149      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1150      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1151      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1152      !!---------------------------------------------------------------------- 
     1153      DO jl = 1, jpl 
     1154         DO jj = Njs0-1, Nje0+1     
     1155            DO ji = Nis0, Nie0 
     1156               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1157            END DO 
     1158         END DO 
     1159         DO jj = Njs0, Nje0     
     1160            DO ji = Nis0, Nie0 
     1161               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1162            END DO 
     1163         END DO 
     1164      END DO 
     1165   END SUBROUTINE icemax3D 
     1166 
     1167   SUBROUTINE icemax4D( pice , pmax ) 
     1168      !!--------------------------------------------------------------------- 
     1169      !!                   ***  ROUTINE icemax4D ***                      
     1170      !! ** Purpose :  compute the max of the 9 points around 
     1171      !!---------------------------------------------------------------------- 
     1172      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1173      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1174      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1175      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1176      !!---------------------------------------------------------------------- 
     1177      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1178      DO jl = 1, jpl 
     1179         DO jk = 1, jlay 
     1180            DO jj = Njs0-1, Nje0+1     
     1181               DO ji = Nis0, Nie0 
     1182                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1183               END DO 
     1184            END DO 
     1185            DO jj = Njs0, Nje0     
     1186               DO ji = Nis0, Nie0 
     1187                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1188               END DO 
     1189            END DO 
     1190         END DO 
     1191      END DO 
     1192   END SUBROUTINE icemax4D 
     1193 
    10581194#else 
    10591195   !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_adv_umx.F90

    r13497 r13886  
    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      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
    136122      ! 
    … …  
    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 
    164       CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
    165       CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     133      END DO    
     134      CALL icemax4D( ze_i , zei_max ) 
     135      CALL icemax4D( ze_s , zes_max ) 
     136      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 
     137      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 
    166138      ! 
    167139      ! 
    … …  
    179151      ENDIF 
    180152      zdt = rDt_ice / REAL(icycle) 
     153      z1_dt = 1._wp / zdt 
    181154 
    182155      ! --- transport --- ! 
    … …  
    207180      !---------------! 
    208181      DO jt = 1, icycle 
     182 
     183         ! diagnostics 
     184         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     185         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     186         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     187            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    209188 
    210189         ! record at_i before advection (for open water) 
    … …  
    377356            ENDIF 
    378357         ENDIF 
     358 
     359         ! --- Lateral boundary conditions --- ! 
     360         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     361            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 & 
     362               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
     363         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     364            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 & 
     365               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
     366         ELSE 
     367            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 ) 
     368         ENDIF 
     369         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 
     370         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 
    379371         ! 
    380372         !== Open water area ==! 
    … …  
    384376               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    385377         END_2D 
    386          CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1.0_wp ) 
    387          ! 
     378         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1._wp ) 
     379         ! 
     380         ! --- diagnostics --- ! 
     381         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     382            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     383         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     384            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     385         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     386            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     387            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
    388388         ! 
    389389         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
    … …  
    445445      !!             work on H (and not V). It is partly related to the multi-category approach 
    446446      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    447       !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
    448       !!             since sv_i and e_i are still good. 
     447      !!             concentration is small). We also limit S and T. 
    449448      !!---------------------------------------------------------------------- 
    450449      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
    … …  
    490489      IF( pamsk == 0._wp ) THEN 
    491490         DO jl = 1, jpl 
    492             DO_2D( 1, 0, 1, 0 ) 
     491            DO_2D( 0, 0, 1, 0 ) 
    493492               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
    494493                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
    … …  
    499498               ENDIF 
    500499               ! 
     500            END_2D 
     501            DO_2D( 1, 0, 0, 0 ) 
    501502               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
    502503                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
    … …  
    533534      IF( PRESENT( pua_ho ) ) THEN 
    534535         DO jl = 1, jpl 
    535             DO_2D( 1, 0, 1, 0 ) 
    536                pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
    537                pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     536            DO_2D( 0, 0, 1, 0 ) 
     537               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 
     538               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 
     539            END_2D 
     540            DO_2D( 1, 0, 0, 0 ) 
     541               pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     542               pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
    538543            END_2D 
    539544         END DO 
    … …  
    549554         END_2D 
    550555      END DO 
    551       CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp ) 
    552556      ! 
    553557   END SUBROUTINE adv_umx 
    … …  
    588592            ! 
    589593            DO jl = 1, jpl              !-- flux in x-direction 
    590                DO_2D( 1, 0, 1, 0 ) 
     594               DO_2D( 1, 1, 1, 0 ) 
    591595                  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) 
    592596               END_2D 
    … …  
    594598            ! 
    595599            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    596                DO_2D( 0, 0, 0, 0 ) 
     600               DO_2D( 1, 1, 0, 0 ) 
    597601                  ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
    598602                     &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    601605               END_2D 
    602606            END DO 
    603             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    604607            ! 
    605608            DO jl = 1, jpl              !-- flux in y-direction 
    606                DO_2D( 1, 0, 1, 0 ) 
     609               DO_2D( 1, 0, 0, 0 ) 
    607610                  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) 
    608611               END_2D 
    … …  
    612615            ! 
    613616            DO jl = 1, jpl              !-- flux in y-direction 
    614                DO_2D( 1, 0, 1, 0 ) 
     617               DO_2D( 1, 0, 1, 1 ) 
    615618                  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) 
    616619               END_2D 
    … …  
    618621            ! 
    619622            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    620                DO_2D( 0, 0, 0, 0 ) 
     623               DO_2D( 0, 0, 1, 1 ) 
    621624                  ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
    622625                     &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    625628               END_2D 
    626629            END DO 
    627             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    628630            ! 
    629631            DO jl = 1, jpl              !-- flux in x-direction 
    630                DO_2D( 1, 0, 1, 0 ) 
     632               DO_2D( 0, 0, 1, 0 ) 
    631633                  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) 
    632634               END_2D 
    … …  
    677679         ! 
    678680         DO jl = 1, jpl 
    679             DO_2D( 1, 0, 1, 0 ) 
     681            DO_2D( 1, 1, 1, 0 ) 
    680682               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) ) 
     683            END_2D 
     684            DO_2D( 1, 0, 1, 1 ) 
    681685               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
    682686            END_2D 
    … …  
    695699            ! 
    696700            DO jl = 1, jpl              !-- flux in x-direction 
    697                DO_2D( 1, 0, 1, 0 ) 
     701               DO_2D( 1, 1, 1, 0 ) 
    698702                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    699703               END_2D 
    … …  
    702706 
    703707            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    704                DO_2D( 0, 0, 0, 0 ) 
     708               DO_2D( 1, 1, 0, 0 ) 
    705709                  ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
    706710                     &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    709713               END_2D 
    710714            END DO 
    711             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    712715 
    713716            DO jl = 1, jpl              !-- flux in y-direction 
    714                DO_2D( 1, 0, 1, 0 ) 
     717               DO_2D( 1, 0, 0, 0 ) 
    715718                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    716719               END_2D 
    … …  
    721724            ! 
    722725            DO jl = 1, jpl              !-- flux in y-direction 
    723                DO_2D( 1, 0, 1, 0 ) 
     726               DO_2D( 1, 0, 1, 1 ) 
    724727                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    725728               END_2D 
    … …  
    728731            ! 
    729732            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    730                DO_2D( 0, 0, 0, 0 ) 
     733               DO_2D( 0, 0, 1, 1 ) 
    731734                  ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
    732735                     &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    … …  
    735738               END_2D 
    736739            END DO 
    737             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp ) 
    738740            ! 
    739741            DO jl = 1, jpl              !-- flux in x-direction 
    740                DO_2D( 1, 0, 1, 0 ) 
     742               DO_2D( 0, 0, 1, 0 ) 
    741743                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    742744               END_2D 
    … …  
    895897         !         
    896898         DO jl = 1, jpl 
    897             DO_2D( 1, 0, 1, 0 ) 
     899            DO_2D( 0, 0, 1, 0 ) 
    898900               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    899901                  &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    … …  
    904906         ! 
    905907         DO jl = 1, jpl 
    906             DO_2D( 1, 0, 1, 0 ) 
     908            DO_2D( 0, 0, 1, 0 ) 
    907909               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    908910               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    … …  
    914916         ! 
    915917         DO jl = 1, jpl 
    916             DO_2D( 1, 0, 1, 0 ) 
     918            DO_2D( 0, 0, 1, 0 ) 
    917919               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    918920               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    928930         ! 
    929931         DO jl = 1, jpl 
    930             DO_2D( 1, 0, 1, 0 ) 
     932            DO_2D( 0, 0, 1, 0 ) 
    931933               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    932934               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    942944         ! 
    943945         DO jl = 1, jpl 
    944             DO_2D( 1, 0, 1, 0 ) 
     946            DO_2D( 0, 0, 1, 0 ) 
    945947               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    946948               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    … …  
    963965      IF( ll_neg ) THEN 
    964966         DO jl = 1, jpl 
    965             DO_2D( 1, 0, 1, 0 ) 
     967            DO_2D( 0, 0, 1, 0 ) 
    966968               IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    967969                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    … …  
    973975      !                                                     !-- High order flux in i-direction  --! 
    974976      DO jl = 1, jpl 
    975          DO_2D( 1, 0, 1, 0 ) 
     977         DO_2D( 0, 0, 1, 0 ) 
    976978            pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
    977979         END_2D 
    … …  
    10311033      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    10321034         DO jl = 1, jpl 
    1033             DO_2D( 1, 0, 1, 0 ) 
     1035            DO_2D( 1, 0, 0, 0 ) 
    10341036               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    10351037                  &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    10391041      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    10401042         DO jl = 1, jpl 
    1041             DO_2D( 1, 0, 1, 0 ) 
     1043            DO_2D( 1, 0, 0, 0 ) 
    10421044               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10431045               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    … …  
    10481050      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    10491051         DO jl = 1, jpl 
    1050             DO_2D( 1, 0, 1, 0 ) 
     1052            DO_2D( 1, 0, 0, 0 ) 
    10511053               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10521054               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    10611063      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    10621064         DO jl = 1, jpl 
    1063             DO_2D( 1, 0, 1, 0 ) 
     1065            DO_2D( 1, 0, 0, 0 ) 
    10641066               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10651067               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    10741076      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    10751077         DO jl = 1, jpl 
    1076             DO_2D( 1, 0, 1, 0 ) 
     1078            DO_2D( 1, 0, 0, 0 ) 
    10771079               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    10781080               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    … …  
    10951097      IF( ll_neg ) THEN 
    10961098         DO jl = 1, jpl 
    1097             DO_2D( 1, 0, 1, 0 ) 
     1099            DO_2D( 1, 0, 0, 0 ) 
    10981100               IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    10991101                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    … …  
    11051107      !                                                     !-- High order flux in j-direction  --! 
    11061108      DO jl = 1, jpl 
    1107          DO_2D( 1, 0, 1, 0 ) 
     1109         DO_2D( 1, 0, 0, 0 ) 
    11081110            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    11091111         END_2D 
    … …  
    11411143      ! -------------------------------------------------- 
    11421144      DO jl = 1, jpl 
    1143          DO_2D( 1, 0, 1, 0 ) 
     1145         DO_2D( 0, 0, 1, 0 ) 
    11441146            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
     1147         END_2D 
     1148         DO_2D( 1, 0, 0, 0 ) 
    11451149            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    11461150         END_2D 
    … …  
    12481252      ! --------------------------------- 
    12491253      DO jl = 1, jpl 
    1250          DO_2D( 1, 0, 1, 0 ) 
     1254         DO_2D( 0, 0, 1, 0 ) 
    12511255            zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    12521256            zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
    … …  
    12591263         END_2D 
    12601264 
    1261          DO_2D( 1, 0, 1, 0 ) 
     1265         DO_2D( 1, 0, 0, 0 ) 
    12621266            zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    12631267            zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
    … …  
    16161620   END SUBROUTINE Hsnow 
    16171621 
     1622   SUBROUTINE icemax3D( pice , pmax ) 
     1623      !!--------------------------------------------------------------------- 
     1624      !!                   ***  ROUTINE icemax3D ***                      
     1625      !! ** Purpose :  compute the max of the 9 points around 
     1626      !!---------------------------------------------------------------------- 
     1627      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1628      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1629      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1630      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1631      !!---------------------------------------------------------------------- 
     1632      DO jl = 1, jpl 
     1633         DO jj = Njs0-1, Nje0+1     
     1634            DO ji = Nis0, Nie0 
     1635               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1636            END DO 
     1637         END DO 
     1638         DO jj = Njs0, Nje0     
     1639            DO ji = Nis0, Nie0 
     1640               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1641            END DO 
     1642         END DO 
     1643      END DO 
     1644   END SUBROUTINE icemax3D 
     1645 
     1646   SUBROUTINE icemax4D( pice , pmax ) 
     1647      !!--------------------------------------------------------------------- 
     1648      !!                   ***  ROUTINE icemax4D ***                      
     1649      !! ** Purpose :  compute the max of the 9 points around 
     1650      !!---------------------------------------------------------------------- 
     1651      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1652      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1653      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1654      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1655      !!---------------------------------------------------------------------- 
     1656      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1657      DO jl = 1, jpl 
     1658         DO jk = 1, jlay 
     1659            DO jj = Njs0-1, Nje0+1     
     1660               DO ji = Nis0, Nie0 
     1661                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1662               END DO 
     1663            END DO 
     1664            DO jj = Njs0, Nje0     
     1665               DO ji = Nis0, Nie0 
     1666                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1667               END DO 
     1668            END DO 
     1669         END DO 
     1670      END DO 
     1671   END SUBROUTINE icemax4D 
    16181672 
    16191673#else 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_rdgrft.F90

    r13495 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_rhg_evp.F90

    r13497 r13886  
    129129      REAl(wp) ::   zbetau, zbetav 
    130130      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    131       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
     131      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    132132      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    133133      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    … …  
    138138      REAL(wp) ::   zshear, zdum1, zdum2 
    139139      ! 
    140       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
     140      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points 
    141141      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    142142      ! 
    … …  
    145145      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    146146      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    147       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     147      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
    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) 
    … …  
    368370 
    369371         END_2D 
    370          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 
    371  
    372          DO_2D( 0, 1, 0, 1 )   ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 ! no vector loop 
     372 
     373         DO_2D( 0, 0, 0, 0 ) 
    373374 
    374375            ! shear**2 at T points (doc eq. A16) 
    … …  
    390391             
    391392            ! delta at T points 
    392             zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    393  
    394             ! P/delta at T points 
    395             zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    396  
     393            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     394 
     395         END_2D 
     396         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 
     397 
     398         ! P/delta at T points 
     399         DO_2D( 1, 1, 1, 1 ) 
     400            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 
     401         END_2D 
     402 
     403         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 
     404 
     405            ! divergence at T points (duplication to avoid communications) 
     406            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     407               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     408               &    ) * r1_e1e2t(ji,jj) 
     409             
     410            ! tension at T points (duplication to avoid communications) 
     411            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     412               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     413               &   ) * r1_e1e2t(ji,jj) 
     414             
    397415            ! alpha for aEVP 
    398416            !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    … …  
    411429             
    412430            ! stress at T points (zkt/=0 if landfast) 
    413             zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    414             zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     431            zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 
     432            zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    415433           
    416434         END_2D 
    417          CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 
    418435 
    419436         ! Save beta at T-points for further computations 
    … …  
    711728            &   ) * r1_e1e2t(ji,jj) 
    712729         zdt2 = zdt * zdt 
     730 
     731         zten_i(ji,jj) = zdt 
    713732          
    714733         ! shear**2 at T points (doc eq. A16) 
    … …  
    726745          
    727746         ! delta at T points 
    728          zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    729          rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    730          pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     747         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta   
     748         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 
     749         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 
    731750 
    732751      END_2D 
    733       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 ) 
     752      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, & 
     753         &                                  zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp ) 
    734754       
    735755      ! --- Store the stress tensor for the next time step --- ! 
    736       CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp ) 
    737756      pstress1_i (:,:) = zs1 (:,:) 
    738757      pstress2_i (:,:) = zs2 (:,:) 
    … …  
    763782      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    764783 
    765       ! --- stress tensor --- ! 
    766       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    767          ! 
    768          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     784      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     785      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     786         ! 
     787         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    769788         !          
    770          DO_2D( 0, 0, 0, 0 ) 
    771             zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    772                &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    773                &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    774  
    775             zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    776  
    777             zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    778  
    779 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    780 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    781 !!               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 
    782 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    783             zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    784             zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    785             zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    786          END_2D 
    787          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 
    788          ! 
    789          CALL iom_put( 'isig1' , zsig1 ) 
    790          CALL iom_put( 'isig2' , zsig2 ) 
    791          CALL iom_put( 'isig3' , zsig3 ) 
    792          ! 
    793          ! Stress tensor invariants (normal and shear stress N/m) 
    794          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    795          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    796  
    797          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     789         DO_2D( 1, 1, 1, 1 ) 
     790             
     791            ! Ice stresses 
     792            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     793            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
     794            ! I know, this can be confusing... 
     795            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     796            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     797            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     798            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     799             
     800            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
     801            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     802            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     803                
     804         END_2D          
     805         ! 
     806         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
     807         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
     808         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     809          
     810         DEALLOCATE ( zsig_I, zsig_II ) 
     811          
     812      ENDIF 
     813 
     814      ! --- Normalized stress tensor principal components --- ! 
     815      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 
     816      ! Recommendation 1 : we use ice strength, not replacement pressure 
     817      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
     818      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     819         ! 
     820         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
     821         !          
     822         DO_2D( 1, 1, 1, 1 ) 
     823             
     824            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     825            !                        and **deformations** at current iterates 
     826            !                        following Lemieux & Dupont (2020) 
     827            zfac             =   zp_delt(ji,jj) 
     828            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
     829            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     830            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     831             
     832            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     833            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                            ! 1st stress invariant, aka average normal stress, aka negative pressure 
     834            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )   ! 2nd  ''       '', aka maximum shear stress 
     835             
     836            ! Normalized  principal stresses (used to display the ellipse) 
     837            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     838            zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     839            zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     840         END_2D               
     841         ! 
     842         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     843         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     844 
     845         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     846          
    798847      ENDIF 
    799848 
    … …  
    931980         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
    932981         ! close file 
    933          IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     982         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
    934983      ENDIF 
    935984       

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/iceitd.F90

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icestp.F90

    r13472 r13886  
    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      ! 
    … …  
    257261      END WHERE 
    258262      ! 
     263      CALL diag_set0                   ! set diag of mass, heat and salt fluxes to 0: needed for Agrif child grids 
     264      ! 
    259265      CALL ice_itd_init                ! ice thickness distribution initialization 
    260266      ! 
    261267      CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 
    262268      ! 
    263       !                                ! Initial sea-ice state 
    264       CALL ice_istate_init 
     269      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
     270      ! 
     271      CALL ice_istate_init             ! Initial sea-ice state 
    265272      IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 
    266273         CALL ice_rst_read( Kbb, Kmm, Kaa )         ! start from a restart file 
    … …  
    271278      CALL ice_var_agg(1) 
    272279      ! 
    273       CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters 
    274       ! 
    275280      CALL ice_dyn_init                ! set ice dynamics parameters 
    276281      ! 
    … …  
    280285      ! 
    281286      CALL ice_dia_init                ! initialization for diags 
     287      ! 
     288      CALL ice_drift_init              ! initialization for diags of conservation 
    282289      ! 
    283290      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
    … …  
    340347      ENDIF 
    341348      ! 
    342       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    343       ! 
    344349      rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse 
    345350      r1_Dt_ice = 1._wp / rDt_ice 
    … …  
    391396      !!               of the time step 
    392397      !!---------------------------------------------------------------------- 
    393       INTEGER  ::   ji, jj      ! dummy loop index 
    394       !!---------------------------------------------------------------------- 
    395       sfx    (:,:) = 0._wp   ; 
    396       sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
    397       sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    398       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    399       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    400       sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
    401       ! 
    402       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    403       wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    404       wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    405       wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    406       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    407       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
    408       wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 
    409       wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 
    410       wfx_snw_sni(:,:) = 0._wp  
    411       wfx_pnd(:,:) = 0._wp 
    412  
    413       hfx_thd(:,:) = 0._wp   ; 
    414       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    415       hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
    416       hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
    417       hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    418       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
    419       hfx_err_dif(:,:) = 0._wp 
    420       wfx_err_sub(:,:) = 0._wp 
    421       ! 
    422       diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp 
    423       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp 
    424  
    425       ! SIMIP diagnostics 
    426       qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 
    427       t_si       (:,:,:) = rt0   ! temp at the ice-snow interface 
    428  
    429       tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    430       cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
    431       qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 
    432       qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
    433       qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F) 
    434       ! 
    435       ! for control checks (ln_icediachk) 
    436       diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp 
    437       diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp 
    438       diag_trp_sv(:,:) = 0._wp 
    439        
     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 
    440459   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 
    441501 
    442502#else 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icetab.F90

    r10069 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd.F90

    r13472 r13886  
    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 
    134       CALL lbc_lnk( 'icethd', zfric, 'T', 1.0_wp ) 
     138      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
    135139      ! 
    136140      !--------------------------------------------------------------------! 
    … …  
    140144         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    141145         ! 
    142          !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    143          !           !  practically no "direct lateral ablation" 
    144          !            
    145          !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    146          !           !  temperature and turbulent mixing (McPhee, 1992) 
    147          ! 
    148146         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    149147         zqld =  tmask(ji,jj,1) * rDt_ice *  & 
    … …  
    151149            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    152150 
    153          ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 
     151         ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
     152         !     (mostly<0 but >0 if supercooling) 
    154153         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) 
    155154         zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
    156  
    157          ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 
     155         zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
     156 
     157         ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
     158         !     (mostly>0 but <0 if supercooling) 
    158159         zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin )  
    159          qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
    160  
    161          qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     160         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
     161          
    162162         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    163163         !                              the freezing point, so that we do not have SST < T_freeze 
    164          !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 
    165  
    166          !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 
    167          qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    168  
    169          ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
    170          ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
    171          IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    172             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 
    173             ELSE                    ;   fhld(ji,jj) = 0._wp 
    174             ENDIF 
     164         !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
     165         !                              The following formulation is ok for both normal conditions and supercooling 
     166         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     167 
     168         ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
     169         !     qlead is the energy received from the atm. in the leads. 
     170         !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
     171         !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
     172         IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     173            ! upper bound for fhld: fhld should be equal to zqld 
     174            !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
     175            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
     176            !                        The following formulation is ok for both normal conditions and supercooling 
     177            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 
     178               &                                 - qsb_ice_bot(ji,jj) ) 
    175179            qlead(ji,jj) = 0._wp 
    176180         ELSE 
    177181            fhld (ji,jj) = 0._wp 
     182            ! upper bound for qlead: qlead should be equal to zqld 
     183            !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
     184            !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
     185            !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
     186            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
     187            !                        The following formulation is ok for both normal conditions and supercooling 
     188            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    178189         ENDIF 
    179190         ! 
    180          ! Net heat flux on top of the ice-ocean [W.m-2] 
    181          ! --------------------------------------------- 
    182          qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
     191         ! If ice is landfast and ice concentration reaches its max 
     192         ! => stop ice formation in open water 
     193         IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
     194         ! 
     195         ! If the grid cell is almost fully covered by ice (no leads) 
     196         ! => stop ice formation in open water 
     197         IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
     198         ! 
     199         ! If ln_leadhfx is false 
     200         ! => do not use energy of the leads to melt sea-ice 
     201         IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
     202         ! 
    183203      END_2D 
    184204       
    … …  
    191211      ENDIF 
    192212 
    193       ! --------------------------------------------------------------------- 
    194       ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    195       ! --------------------------------------------------------------------- 
    196       !     First  step here              :  non solar + precip - qlead - qsensible 
    197       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)  
    198       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    199       qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    200          &             - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    201          &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux 
    202          &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt  
    203       !                                                                           !    (fhld should be 0 while bott growth) 
    204213      !-------------------------------------------------------------------------------------------! 
    205214      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
    … …  
    231240                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    232241                              CALL ice_thd_pnd                          ! Melt ponds formation 
    233                               CALL ice_thd_ent( e_i_1d(1:npti,:), .true. )      ! Ice enthalpy remapping 
     242                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    234243            ENDIF 
    235244                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !     
    … …  
    254263      ! 
    255264      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
     265      ! 
     266                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
     267      ! 
     268      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation      
    256269      ! 
    257270      ! convergence tests 
    … …  
    418431         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       ) 
    419432         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   ) 
    420          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     ) 
    421433         ! 
    422434         ! ocean surface fields 
    423435         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 
    424436         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 
     437         CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 
    425438         ! 
    426439         ! to update ice age 
    … …  
    510523         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     ) 
    511524         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 
    512          CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   ) 
    513525         ! 
    514526         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) ) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd_dh.F90

    r13472 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd_do.F90

    r13295 r13886  
    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)   & 
    … …  
    198200      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 
    199201      !------------------------------------------------------------------------------! 
    200       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 
     202      ! it occurs if cooling 
    201203 
    202204      ! Identify grid points where new ice forms 
    203205      npti = 0   ;   nptidx(:) = 0 
    204206      DO_2D( 1, 1, 1, 1 ) 
    205          IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 
     207         IF ( qlead(ji,jj)  <  0._wp ) THEN 
    206208            npti = npti + 1 
    207209            nptidx( npti ) = (jj - 1) * jpi + ji 
    … …  
    385387            END DO 
    386388            ! --- Ice enthalpy remapping --- ! 
    387             CALL ice_thd_ent( ze_i_2d(1:npti,:,jl), .false. )  
     389            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) )  
    388390         END DO 
    389391 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd_ent.F90

    r13472 r13886  
    3838CONTAINS 
    3939  
    40    SUBROUTINE ice_thd_ent( qnew, compute_hfx_err ) 
     40   SUBROUTINE ice_thd_ent( qnew ) 
    4141      !!------------------------------------------------------------------- 
    4242      !!               ***   ROUTINE ice_thd_ent  *** 
    … …  
    6464      !!------------------------------------------------------------------- 
    6565      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   qnew             ! new enthlapies (J.m-3, remapped) 
    66       LOGICAL, INTENT(in)                     ::   compute_hfx_err  ! determines whether to compute diag. 
    67                                                                     ! error or not 
    6866      ! 
    6967      INTEGER  :: ji         !  dummy loop indices 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/iceupdate.F90

    r13497 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/NST/agrif_user.F90

    r13472 r13886  
    405405         use_sign_north = .TRUE. 
    406406         sign_north = -1. 
     407         CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b)   ! must be called before unb_id to define ubdy 
     408         CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b)   ! must be called before vnb_id to define vbdy 
    407409         CALL Agrif_Bc_variable(        unb_id,calledweight=1.,procname=interpunb ) 
    408410         CALL Agrif_Bc_variable(        vnb_id,calledweight=1.,procname=interpvnb ) 
    409          CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 
    410          CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 
    411411         use_sign_north = .FALSE. 
    412412         ubdy(:,:) = 0._wp 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/BDY/bdyice.F90

    r13472 r13886  
    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 
    … …  
    110108      ! 
    111109      ! controls 
    112       IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
    113       IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    114       IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    115       IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
     110      IF( ln_icectl )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )   ! prints 
     111      IF( ln_timing )   CALL timing_stop ('bdy_ice_thd')                                       ! timing 
    116112      ! 
    117113   END SUBROUTINE bdy_ice 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DIA/diaptr.F90

    r13295 r13886  
    3636   END INTERFACE 
    3737 
    38    PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines 
    39    PUBLIC   ptr_sjk        !  
    40    PUBLIC   dia_ptr_init   ! call in memogcm 
    4138   PUBLIC   dia_ptr        ! call in step module 
    4239   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines 
    4340 
    44    !                                  !!** namelist  namptr  ** 
    4541   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.) 
    4642   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional) 
    4743 
    48    LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini) 
    49    INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc) 
     44   LOGICAL, PUBLIC ::   l_diaptr       !: tracers  trend flag 
    5045 
    5146   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
    … …  
    5954   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d 
    6055 
    61    LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini) 
     56   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag 
    6257    
    6358   !! * Substitutions 
    … …  
    8883      ! 
    8984      !overturning calculation 
    90       REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse 
    91       REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function 
    92  
    93       REAL(wp), DIMENSION(jpi,jpj,jpk,nptr)  :: z4d1, z4d2 
    94       REAL(wp), DIMENSION(jpi,jpj,nptr)      :: z3dtr ! i-mean T and S, j-Stream-Function 
     85      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::  sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse 
     86      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function 
     87 
     88      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  z4d1, z4d2 
     89      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr 
    9590      !!---------------------------------------------------------------------- 
    9691      ! 
    9792      IF( ln_timing )   CALL timing_start('dia_ptr') 
    9893 
    99       IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init 
    100       ! 
    101       IF( .NOT. l_diaptr )   RETURN 
    102  
     94      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init   ! -> will define l_diaptr and nbasin 
     95      ! 
     96      IF( .NOT. l_diaptr ) THEN 
     97         IF( ln_timing ) CALL timing_stop('dia_ptr') 
     98         RETURN 
     99      ENDIF 
     100      ! 
     101      ALLOCATE( z3dtr(jpi,jpj,nbasin) ) 
     102      ! 
    103103      IF( PRESENT( pvtr ) ) THEN 
    104104         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF 
    105             DO jn = 1, nptr                                    ! by sub-basins 
     105            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) ) 
     106            DO jn = 1, nbasin                                    ! by sub-basins 
    106107               z4d1(1,:,:,jn) =  ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )  ! zonal cumulative effective transport excluding closed seas 
    107108               DO jk = jpkm1, 1, -1  
    … …  
    113114            END DO 
    114115            CALL iom_put( 'zomsf', z4d1 * rc_sv ) 
     116            DEALLOCATE( z4d1 ) 
    115117         ENDIF 
    116118         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   & 
    … …  
    127129         ENDIF 
    128130         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN 
    129             DO jn = 1, nptr 
     131            DO jn = 1, nbasin 
     132               ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin),   & 
     133                  &                          zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) ) 
    130134               sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
    131135               r1_sjk(:,:,jn) = 0._wp 
    … …  
    137141               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 ) 
    138142               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 ) 
     143               DEALLOCATE( sjk, r1_sjk, v_msf, zt_jk, zs_jk ) 
    139144               ! 
    140145            ENDDO 
    141             DO jn = 1, nptr 
     146            DO jn = 1, nbasin 
    142147               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    143148               DO ji = 1, jpi 
    … …  
    146151            ENDDO 
    147152            CALL iom_put( 'sophtove', z3dtr ) 
    148             DO jn = 1, nptr 
     153            DO jn = 1, nbasin 
    149154               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    150155               DO ji = 1, jpi 
    … …  
    157162         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN 
    158163            ! Calculate barotropic heat and salt transport here  
    159             DO jn = 1, nptr 
     164            DO jn = 1, nbasin 
     165               ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) ) 
    160166               sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) ) 
    161167               r1_sjk(:,1,jn) = 0._wp 
    … …  
    167173               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn) 
    168174               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn) 
     175               DEALLOCATE( sjk, r1_sjk ) 
    169176               ! 
    170177            ENDDO 
    171             DO jn = 1, nptr 
     178            DO jn = 1, nbasin 
    172179               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    173180               DO ji = 1, jpi 
    … …  
    176183            ENDDO 
    177184            CALL iom_put( 'sophtbtr', z3dtr ) 
    178             DO jn = 1, nptr 
     185            DO jn = 1, nbasin 
    179186               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    180187               DO ji = 1, jpi 
    … …  
    190197         zts(:,:,:,:) = 0._wp 
    191198         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface  
     199            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) ) 
    192200            DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
    193201               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm) 
    … …  
    197205            END_3D 
    198206            ! 
    199             DO jn = 1, nptr 
     207            DO jn = 1, nbasin 
    200208               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
     209               DO ji = 1, jpi 
     210                  zmask(ji,:,:) = zmask(1,:,:) 
     211               ENDDO 
    201212               z4d1(:,:,:,jn) = zmask(:,:,:) 
    202213            ENDDO 
    203214            CALL iom_put( 'zosrf', z4d1 ) 
    204215            ! 
    205             DO jn = 1, nptr 
     216            DO jn = 1, nbasin 
    206217               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) & 
    207218                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 ) 
    … …  
    212223            CALL iom_put( 'zotem', z4d2 ) 
    213224            ! 
    214             DO jn = 1, nptr 
     225            DO jn = 1, nbasin 
    215226               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) & 
    216227                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 ) 
    … …  
    220231            ENDDO 
    221232            CALL iom_put( 'zosal', z4d2 ) 
     233            DEALLOCATE( z4d1, z4d2 ) 
    222234            ! 
    223235         ENDIF 
    … …  
    226238         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN   
    227239            !  
    228             DO jn = 1, nptr 
     240            DO jn = 1, nbasin 
    229241               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    230242               DO ji = 1, jpi 
    … …  
    233245            ENDDO 
    234246            CALL iom_put( 'sophtadv', z3dtr ) 
    235             DO jn = 1, nptr 
     247            DO jn = 1, nbasin 
    236248               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    237249               DO ji = 1, jpi 
    … …  
    244256         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN   
    245257            !  
    246             DO jn = 1, nptr 
     258            DO jn = 1, nbasin 
    247259               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    248260               DO ji = 1, jpi 
    … …  
    251263            ENDDO 
    252264            CALL iom_put( 'sophtldf', z3dtr ) 
    253             DO jn = 1, nptr 
     265            DO jn = 1, nbasin 
    254266               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    255267               DO ji = 1, jpi 
    … …  
    262274         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN   
    263275            !  
    264             DO jn = 1, nptr 
     276            DO jn = 1, nbasin 
    265277               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    266278               DO ji = 1, jpi 
    … …  
    269281            ENDDO 
    270282            CALL iom_put( 'sophteiv', z3dtr ) 
    271             DO jn = 1, nptr 
     283            DO jn = 1, nbasin 
    272284               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    273285               DO ji = 1, jpi 
    … …  
    287299             CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) ) 
    288300             CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) ) 
    289              DO jn = 1, nptr 
     301             DO jn = 1, nbasin 
    290302                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    291303                DO ji = 1, jpi 
    … …  
    294306             ENDDO 
    295307             CALL iom_put( 'sophtvtr', z3dtr ) 
    296              DO jn = 1, nptr 
     308             DO jn = 1, nbasin 
    297309               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    298310               DO ji = 1, jpi 
    … …  
    311323      ENDIF 
    312324      ! 
     325      DEALLOCATE( z3dtr ) 
     326      ! 
    313327      IF( ln_timing )   CALL timing_stop('dia_ptr') 
    314328      ! 
    … …  
    320334      !!                  ***  ROUTINE dia_ptr_init  *** 
    321335      !!                    
    322       !! ** Purpose :   Initialization, namelist read 
     336      !! ** Purpose :   Initialization 
    323337      !!---------------------------------------------------------------------- 
    324338      INTEGER ::  inum, jn           ! local integers 
    … …  
    326340      REAL(wp), DIMENSION(jpi,jpj) :: zmsk 
    327341      !!---------------------------------------------------------------------- 
    328  
    329       l_diaptr = .FALSE. 
    330       IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  & 
    331          &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  & 
    332          &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  & 
    333          &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &  
    334          &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  & 
    335          &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE. 
    336  
     342       
     343      ! l_diaptr is defined with iom_use 
     344      !   --> dia_ptr_init must be done after the call to iom_init 
     345      !   --> cannot be .TRUE. without cpp key: key_iom -->  nbasin define by iom_init is initialized 
     346      l_diaptr = iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  & 
     347         &       iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  & 
     348         &       iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  & 
     349         &       iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &  
     350         &       iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  & 
     351         &       iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' )  
    337352  
    338353      IF(lwp) THEN                     ! Control print 
    … …  
    340355         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 
    341356         WRITE(numout,*) '~~~~~~~~~~~~' 
    342          WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
    343357         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr 
    344358      ENDIF 
    … …  
    347361         ! 
    348362         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
    349  
     363         ! 
    350364         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt 
    351365         rc_ggram = rc_ggram * rho0              ! conversion from m3/s to Gg/s 
    … …  
    354368 
    355369         btmsk(:,:,1) = tmask_i(:,:)                  
    356          CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  ) 
    357          CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
    358          CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
    359          CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
    360          CALL iom_close( inum ) 
    361          btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin 
    362          DO jn = 2, nptr 
    363             btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only 
     370         IF( nbasin == 5 ) THEN   ! nbasin has been initialized in iom_init to define the axis "basin" 
     371            CALL iom_open( 'subbasins', inum ) 
     372            CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
     373            CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
     374            CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
     375            CALL iom_close( inum ) 
     376            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )            ! Indo-Pacific basin 
     377         ENDIF 
     378         DO jn = 2, nbasin 
     379            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)                 ! interior domain only 
    364380         END DO 
    365381         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations 
    … …  
    370386         END WHERE 
    371387         btmsk34(:,:,1) = btmsk(:,:,1)                  
    372          DO jn = 2, nptr 
    373             btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)               ! interior domain only 
     388         DO jn = 2, nbasin 
     389            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)                  ! interior domain only 
    374390         ENDDO 
    375391 
    … …  
    405421      IF( cptr == 'adv' ) THEN 
    406422         IF( ktra == jp_tem )  THEN 
    407              DO jn = 1, nptr 
     423             DO jn = 1, nbasin 
    408424                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    409425             ENDDO 
    410426         ENDIF 
    411427         IF( ktra == jp_sal )  THEN 
    412              DO jn = 1, nptr 
     428             DO jn = 1, nbasin 
    413429                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    414430             ENDDO 
    … …  
    418434      IF( cptr == 'ldf' ) THEN 
    419435         IF( ktra == jp_tem )  THEN 
    420              DO jn = 1, nptr 
     436             DO jn = 1, nbasin 
    421437                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    422438             ENDDO 
    423439         ENDIF 
    424440         IF( ktra == jp_sal )  THEN 
    425              DO jn = 1, nptr 
     441             DO jn = 1, nbasin 
    426442                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    427443             ENDDO 
    … …  
    431447      IF( cptr == 'eiv' ) THEN 
    432448         IF( ktra == jp_tem )  THEN 
    433              DO jn = 1, nptr 
     449             DO jn = 1, nbasin 
    434450                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    435451             ENDDO 
    436452         ENDIF 
    437453         IF( ktra == jp_sal )  THEN 
    438              DO jn = 1, nptr 
     454             DO jn = 1, nbasin 
    439455                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    440456             ENDDO 
    … …  
    444460      IF( cptr == 'vtr' ) THEN 
    445461         IF( ktra == jp_tem )  THEN 
    446              DO jn = 1, nptr 
     462             DO jn = 1, nbasin 
    447463                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    448464             ENDDO 
    449465         ENDIF 
    450466         IF( ktra == jp_sal )  THEN 
    451              DO jn = 1, nptr 
     467             DO jn = 1, nbasin 
    452468                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    453469             ENDDO 
    … …  
    467483      ierr(:) = 0 
    468484      ! 
     485      ! nbasin has been initialized in iom_init to define the axis "basin" 
     486      ! 
    469487      IF( .NOT. ALLOCATED( btmsk ) ) THEN 
    470          ALLOCATE( btmsk(jpi,jpj,nptr)    , btmsk34(jpi,jpj,nptr),   & 
    471             &      hstr_adv(jpj,jpts,nptr), hstr_eiv(jpj,jpts,nptr), & 
    472             &      hstr_ove(jpj,jpts,nptr), hstr_btr(jpj,jpts,nptr), & 
    473             &      hstr_ldf(jpj,jpts,nptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1)  ) 
     488         ALLOCATE( btmsk(jpi,jpj,nbasin)    , btmsk34(jpi,jpj,nbasin),   & 
     489            &      hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), & 
     490            &      hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), & 
     491            &      hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1)  ) 
    474492            ! 
    475493         ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2)) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DIU/diu_bulk.F90

    r13295 r13886  
    2222    
    2323   ! Namelist parameters 
    24    LOGICAL, PUBLIC :: ln_diurnal 
    25    LOGICAL, PUBLIC :: ln_diurnal_only 
     24   LOGICAL, PUBLIC :: ln_diurnal      = .false.   ! force definition if diurnal_sst_bulk_init is not called 
     25   LOGICAL, PUBLIC :: ln_diurnal_only = .false.   ! force definition if diurnal_sst_bulk_init is not called 
    2626 
    2727   ! Parameters 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DIU/diu_coolskin.F90

    r13295 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DOM/closea.F90

    r13286 r13886  
    3838   LOGICAL, PUBLIC :: ln_clo_rnf       !: closed sea treated as runoff (update rnf mask) 
    3939 
    40    LOGICAL, PUBLIC :: l_sbc_clo  !: T => net evap/precip over closed seas spread outover the globe/river mouth 
    41    LOGICAL, PUBLIC :: l_clo_rnf  !: T => Some closed seas output freshwater (RNF) to specified runoff points. 
    42  
    43    INTEGER, PUBLIC :: ncsg      !: number of closed seas global mappings (inferred from closea_mask_glo field) 
    44    INTEGER, PUBLIC :: ncsr      !: number of closed seas rnf    mappings (inferred from closea_mask_rnf field) 
    45    INTEGER, PUBLIC :: ncse      !: number of closed seas empmr  mappings (inferred from closea_mask_emp field) 
     40   ! WARNING: keep default definitions in the following lines as dom_clo is called only if ln_closea = .true. 
     41   LOGICAL, PUBLIC :: l_sbc_clo = .FALSE.   !: T => net evap/precip over closed seas spread outover the globe/river mouth 
     42   LOGICAL, PUBLIC :: l_clo_rnf = .FALSE.   !: T => Some closed seas output freshwater (RNF) to specified runoff points. 
     43 
     44   INTEGER, PUBLIC :: ncsg = 0   !: number of closed seas global mappings (inferred from closea_mask_glo field) 
     45   INTEGER, PUBLIC :: ncsr = 0   !: number of closed seas rnf    mappings (inferred from closea_mask_rnf field) 
     46   INTEGER, PUBLIC :: ncse = 0   !: number of closed seas empmr  mappings (inferred from closea_mask_emp field) 
    4647 
    4748   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_opnsea, mask_csundef  !: mask defining the open sea and the undefined closed sea 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DOM/daymod.F90

    r13286 r13886  
    8282      ndt05   = NINT( 0.5 * rn_Dt  ) 
    8383 
    84       IF( .NOT. l_offline )   CALL day_rst( nit000, 'READ' ) 
    85  
     84      lrst_oce = .NOT. l_offline   ! force definition of offline 
     85      IF( lrst_oce )   CALL day_rst( nit000, 'READ' ) 
     86       
    8687      ! set the calandar from ndastp (read in restart file and namelist) 
    8788      nyear   =   ndastp / 10000 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DOM/dom_oce.F90

    r13286 r13886  
    220220 
    221221   !!---------------------------------------------------------------------- 
     222   !! variable defined here to avoid circular dependencies... 
     223   !! --------------------------------------------------------------------- 
     224   INTEGER, PUBLIC ::   nbasin         ! number of basin to be considered in diaprt (glo, atl, pac, ind, ipc) 
     225 
     226   !!---------------------------------------------------------------------- 
    222227   !! agrif domain 
    223228   !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DOM/domain.F90

    r13472 r13886  
    120120         WRITE(numout,*)     '         cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg 
    121121      ENDIF 
    122       nn_wxios = 0 
    123       ln_xios_read = .FALSE. 
    124122      ! 
    125123      !           !==  Reference coordinate system  ==! 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN/divhor.F90

    r13497 r13886  
    7878      ! 
    7979      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                    !==  Horizontal divergence  ==! 
    80          hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      & 
     80         hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      & 
    8181            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
    8282            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN/dynspg_ts.F90

    r13497 r13886  
    917917               CALL iom_get( numror, jpdom_auto, 'ub2_i_b'  , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )    
    918918               CALL iom_get( numror, jpdom_auto, 'vb2_i_b'  , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 
     919            ELSE 
     920               ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif 
    919921            ENDIF 
    920922#endif 
    … …  
    922924            IF(lwp) WRITE(numout,*) 
    923925            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0' 
    924             ub2_b (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif 
    925             un_adv(:,:) = 0._wp   ;   vn_adv(:,:) = 0._wp   ! used in the 1st interpol of agrif 
    926             un_bf (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif 
     926            ub2_b  (:,:) = 0._wp   ;   vb2_b (:,:) = 0._wp   ! used in the 1st interpol of agrif 
     927            un_adv (:,:) = 0._wp   ;   vn_adv (:,:) = 0._wp   ! used in the 1st interpol of agrif 
     928            un_bf  (:,:) = 0._wp   ;   vn_bf (:,:) = 0._wp   ! used in the 1st update   of agrif 
    927929#if defined key_agrif 
    928             IF ( .NOT.Agrif_Root() ) THEN 
    929                ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif 
    930             ENDIF 
     930            ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif 
    931931#endif 
    932932         ENDIF 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN/dynvor.F90

    r13497 r13886  
    217217      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    218218      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
    219       REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace 
    220       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace 
     219      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx, zwy, zwt   ! 2D workspace 
     220      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz             ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 
    221221      !!---------------------------------------------------------------------- 
    222222      ! 
    … …  
    533533      REAL(wp) ::   zua, zva     ! local scalars 
    534534      REAL(wp) ::   zmsk, ze3f   ! local scalars 
    535       REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy , z1_e3f 
    536       REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse 
    537       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz 
     535      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy , z1_e3f 
     536      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse 
     537      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 
    538538      !!---------------------------------------------------------------------- 
    539539      ! 
    … …  
    677677      REAL(wp) ::   zua, zva       ! local scalars 
    678678      REAL(wp) ::   zmsk, z1_e3t   ! local scalars 
    679       REAL(wp), DIMENSION(jpi,jpj)     ::   zwx , zwy  
    680       REAL(wp), DIMENSION(jpi,jpj)     ::   ztnw, ztne, ztsw, ztse 
    681       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz 
     679      REAL(wp), DIMENSION(jpi,jpj)       ::   zwx , zwy  
     680      REAL(wp), DIMENSION(jpi,jpj)       ::   ztnw, ztne, ztsw, ztse 
     681      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwz   ! 3D workspace, avoid lbc_lnk on jpk that is not defined 
    682682      !!---------------------------------------------------------------------- 
    683683      ! 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN/wet_dry.F90

    r13497 r13886  
    5757   REAL(wp), PUBLIC  ::   ssh_ref     !: height of z=0 with respect to the geoid;  
    5858 
    59    LOGICAL,  PUBLIC  ::   ll_wd       !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl 
     59   LOGICAL,  PUBLIC  ::   ll_wd = .FALSE. !: Wetting/drying activation switch (ln_wd_il or ln_wd_dl) <- default def if wad_init not called 
    6060 
    6161   PUBLIC   wad_init                  ! initialisation routine called by step.F90 
    … …  
    111111 
    112112      r_rn_wdmin1 = 1 / rn_wdmin1 
    113       ll_wd = .FALSE. 
    114113      IF( ln_wd_il .OR. ln_wd_dl ) THEN 
    115114         ll_wd = .TRUE. 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/FLO/flo_oce.F90

    r11536 r13886  
    1919   !! ---------------- 
    2020   LOGICAL, PUBLIC ::   ln_floats   !: Activate floats or not 
    21    INTEGER, PUBLIC ::   jpnfl       !: total number of floats during the run 
     21   INTEGER, PUBLIC ::   jpnfl = 0   !: total number of floats during the run 
    2222   INTEGER, PUBLIC ::   jpnnewflo   !: number of floats added in a new run 
    2323   INTEGER, PUBLIC ::   jpnrstflo   !: number of floats for the restart 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbtrj.F90

    r13062 r13886  
    3535   PUBLIC   icb_trj_end     ! routine called in icbstp.F90 module 
    3636 
    37    INTEGER ::   num_traj 
     37   INTEGER ::   num_traj = 0 
    3838   INTEGER ::   n_dim, m_dim 
    3939   INTEGER ::   ntrajid 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/IOM/iom.F90

    r13472 r13886  
    123123      REAL(wp), DIMENSION(2,jpkam1)         :: za_bnds   ! ABL vertical boundaries 
    124124      LOGICAL ::   ll_closedef = .TRUE. 
     125      LOGICAL ::   ll_exist 
    125126      !!---------------------------------------------------------------------- 
    126127      ! 
    … …  
    230231          CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds ) 
    231232 
    232           CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) ) 
     233          CALL iom_set_axis_attr(  "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) ) 
    233234# if defined key_si3 
    234235          CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) ) 
    … …  
    243244          CALL iom_set_axis_attr( "iax_26C", (/ REAL(26,wp) /) )   ! strange syntaxe and idea... 
    244245          CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )   ! strange syntaxe and idea... 
    245           CALL iom_set_axis_attr( "basin"  , (/ (REAL(ji,wp), ji=1,5) /) ) 
     246          ! for diaprt, we need to define an axis which size can be 1 (default) or 5 (if the file subbasins.nc exists) 
     247          INQUIRE( FILE = 'subbasins.nc', EXIST = ll_exist ) 
     248          nbasin = 1 + 4 * COUNT( (/ll_exist/) ) 
     249          CALL iom_set_axis_attr( "basin"  , (/ (REAL(ji,wp), ji=1,nbasin) /) ) 
    246250      ENDIF 
    247251      ! 
    … …  
    19101914      IF( iom_use(cdname) ) THEN 
    19111915#if defined key_iomput 
    1912          IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 
    1913             CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d 
    1914          ELSE 
    1915             CALL xios_send_field( cdname, pfield2d ) 
    1916          ENDIF 
     1916         CALL xios_send_field( cdname, pfield2d ) 
    19171917#else 
    19181918         WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19261926      IF( iom_use(cdname) ) THEN 
    19271927#if defined key_iomput 
    1928          IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 
    1929             CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d 
    1930          ELSE 
    1931             CALL xios_send_field( cdname, pfield2d ) 
    1932          ENDIF 
     1928         CALL xios_send_field( cdname, pfield2d ) 
    19331929#else 
    19341930         WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19421938      IF( iom_use(cdname) ) THEN 
    19431939#if defined key_iomput 
    1944          IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 
    1945             CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d 
    1946          ELSE 
    1947             CALL xios_send_field( cdname, pfield3d ) 
    1948          ENDIF 
     1940         CALL xios_send_field( cdname, pfield3d ) 
    19491941#else 
    19501942         WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19581950      IF( iom_use(cdname) ) THEN 
    19591951#if defined key_iomput 
    1960          IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 
    1961             CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d 
    1962          ELSE 
    1963             CALL xios_send_field( cdname, pfield3d ) 
    1964          ENDIF 
     1952         CALL xios_send_field( cdname, pfield3d ) 
    19651953#else 
    19661954         WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19741962      IF( iom_use(cdname) ) THEN 
    19751963#if defined key_iomput 
    1976          IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 
    1977             CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d 
    1978          ELSE 
    1979             CALL xios_send_field (cdname, pfield4d ) 
    1980          ENDIF 
     1964         CALL xios_send_field (cdname, pfield4d ) 
    19811965#else 
    19821966         WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    19901974      IF( iom_use(cdname) ) THEN 
    19911975#if defined key_iomput 
    1992          IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 
    1993             CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d 
    1994          ELSE 
    1995             CALL xios_send_field (cdname, pfield4d ) 
    1996          ENDIF 
     1976         CALL xios_send_field (cdname, pfield4d ) 
    19971977#else 
    19981978         WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings 
    … …  
    22002180      ! 
    22012181      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) 
    2202       CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0) 
     2182      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 
    22032183!don't define lon and lat for restart reading context.  
    22042184      IF ( .NOT.ldrxios ) & 
    … …  
    22992279      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) 
    23002280      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) 
    2301       CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0) 
     2281      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 
    23022282      CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp),   & 
    23032283         &                             latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))   
    2304       CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj_0) 
     2284      CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo) 
    23052285      ! 
    23062286      CALL iom_update_file_name('ptr') 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/IOM/iom_def.F90

    r13286 r13886  
    3333   INTEGER, PUBLIC            ::   iom_open_init = 0   !: used to initialize iom_file(:)%nfid to 0 
    3434!XIOS write restart    
    35    LOGICAL, PUBLIC            ::   lwxios          !: write single file restart using XIOS 
    36    INTEGER, PUBLIC            ::   nxioso          !: type of restart file when writing using XIOS 1 - single, 2 - multiple 
     35   LOGICAL, PUBLIC            ::   lwxios = .FALSE.    !: write single file restart using XIOS 
     36   INTEGER, PUBLIC            ::   nxioso = 0          !: type of restart file when writing using XIOS 1 - single, 2 - multiple 
    3737!XIOS read restart    
    38    LOGICAL, PUBLIC            ::   lrxios          !: read single file restart using XIOS 
     38   LOGICAL, PUBLIC            ::   lrxios = .FALSE.     !: read single file restart using XIOS 
    3939   LOGICAL, PUBLIC            ::   lxios_sini = .FALSE. ! is restart in a single file 
    4040   LOGICAL, PUBLIC            ::   lxios_set  = .FALSE.  

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ISF/isf_oce.F90

    r12077 r13886  
    7474   ! 
    7575   ! 2.1 -------- ice shelf cavity parameter -------------- 
    76    LOGICAL , PUBLIC            :: l_isfoasis 
     76   LOGICAL , PUBLIC            :: l_isfoasis = .FALSE. 
    7777   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load 
    7878   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/LBC/lib_mpp.F90

    r13438 r13886  
    511511            ALLOCATE(todelay(idvar)%y1d(isz)) 
    512512            todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp)   ! create %y1d, complex variable needed by mpi_sumdd 
     513            ndelayid(idvar) = MPI_REQUEST_NULL                             ! initialised request to a valid value 
    513514         END IF 
    514515      ENDIF 
    … …  
    518519         ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz))   ! allocate also %z1d as used for the restart 
    519520         CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr )   ! get %y1d 
    520          todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp)      ! define %z1d from %y1d 
    521       ENDIF 
    522  
    523       IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
     521         ndelayid(idvar) = MPI_REQUEST_NULL 
     522      ENDIF 
     523 
     524      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
    524525 
    525526      ! send back pout from todelay(idvar)%z1d defined at previous call 
    … …  
    530531      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    531532      CALL  mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) 
    532       ndelayid(idvar) = 1 
     533      ndelayid(idvar) = MPI_REQUEST_NULL 
    533534      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    534535# else 
    … …  
    591592            DEALLOCATE(todelay(idvar)%z1d) 
    592593            ndelayid(idvar) = -1                                      ! do as if we had no restart 
     594         ELSE 
     595            ndelayid(idvar) = MPI_REQUEST_NULL 
    593596         END IF 
    594597      ENDIF 
    … …  
    598601         ALLOCATE(todelay(idvar)%z1d(isz)) 
    599602         CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr )   ! get %z1d 
    600       ENDIF 
    601  
    602       IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
     603         ndelayid(idvar) = MPI_REQUEST_NULL 
     604      ENDIF 
     605 
     606      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received 
    603607 
    604608      ! send back pout from todelay(idvar)%z1d defined at previous call 
    … …  
    606610 
    607611      ! send p_in into todelay(idvar)%z1d with a non-blocking communication 
     612      ! (PM) Should we get rid of MPI2 option ? MPI3 was release in 2013. Who is still using MPI2 ? 
    608613# if defined key_mpi2 
    609614      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    610       CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ndelayid(idvar), ierr ) 
     615      CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ierr ) 
    611616      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    612617# else 
    … …  
    631636      !!---------------------------------------------------------------------- 
    632637#if defined key_mpp_mpi 
    633       IF( ndelayid(kid) /= -2 ) THEN   
    634 #if ! defined key_mpi2 
    635          IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
    636          CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr )                        ! make sure todelay(kid) is received 
    637          IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 
    638 #endif 
    639          IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d 
    640          ndelayid(kid) = -2   ! add flag to know that mpi_wait was already called on kid 
    641       ENDIF 
     638      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 
     639      ! test on ndelayid(kid) useless as mpi_wait return immediatly if the request handle is MPI_REQUEST_NULL 
     640      CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! after this ndelayid(kid) = MPI_REQUEST_NULL 
     641      IF( ln_timing ) CALL tic_tac( .FALSE., ld_global = .TRUE.) 
     642      IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d 
    642643#endif 
    643644   END SUBROUTINE mpp_delay_rcv 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/LDF/ldftra.F90

    r13497 r13886  
    246246      ENDIF 
    247247      ! 
    248       IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                & 
    249            &            CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 
    250       IF( ln_isfcav .AND. ln_traldf_triad ) & 
    251            &            CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 
     248      IF( ln_isfcav .AND. ln_traldf_triad )   CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 
    252249           ! 
    253250      IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 
    … …  
    541538         IF( ln_traldf_blp )   CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 
    542539         ! 
     540         IF( .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )   & 
     541           &                  CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 
    543542         !                                != allocate the aei arrays 
    544543         ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/fldread.F90

    r13295 r13886  
    216216                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &             
    217217                     & sd(jf)%nrec(1,ibb), sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday 
    218                   WRITE(numout, *) '      zt_offset is : ',zt_offset 
     218                  IF( zt_offset /= 0._wp )   WRITE(numout, *) '      zt_offset is : ', zt_offset 
    219219               ENDIF 
    220220               ! temporal interpolation weights 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcfwb.F90

    r13286 r13886  
    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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcmod.F90

    r13477 r13886  
    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 
    … …  
    252252      sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting 
    253253      fmmflx(:,:) = 0._wp           !* freezing minus melting flux 
     254      cloud_fra(:,:) = pp_cldf      !* cloud fraction over sea ice (used in si3) 
    254255 
    255256      taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart) 
    … …  
    336337      IF( l_sbc_clo   )   CALL sbc_clo_init              ! closed sea surface initialisation 
    337338      ! 
    338       IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization 
    339  
    340       IF( ln_abl      )   CALL sbc_abl_init            ! Atmospheric Boundary Layer (ABL) 
    341  
    342       IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization 
     339      IF( ln_blk      )   CALL sbc_blk_init              ! bulk formulae initialization 
     340 
     341      IF( ln_abl      )   CALL sbc_abl_init              ! Atmospheric Boundary Layer (ABL) 
     342 
     343      IF( ln_ssr      )   CALL sbc_ssr_init              ! Sea-Surface Restoring initialization 
    343344      ! 
    344345      ! 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcwave.F90

    r13497 r13886  
    106106      !!--------------------------------------------------------------------- 
    107107      ! 
    108       ALLOCATE( ze3divh(jpi,jpj,jpk) ) 
     108      ALLOCATE( ze3divh(jpi,jpj,jpkm1) )   ! jpkm1 -> avoid lbc_lnk on jpk that is not defined 
    109109      ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) ) 
    110110      ! 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfdrg.F90

    r13497 r13886  
    383383      IF(ll_bot)   zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:)                         ! x seafloor mask 
    384384      ! 
     385      l_log_not_linssh = .FALSE.    ! default definition 
    385386      ! 
    386387      SELECT CASE( ndrg ) 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfgls.F90

    r13506 r13886  
    327327      ! at k=2, set de/dz=Fw 
    328328      !cbr 
    329       zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag 
    330       zd_lw(:,:,2) = 0._wp 
     329      DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo 
     330         zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 
     331         zd_lw(ji,jj,2) = 0._wp 
     332      END_2D 
    331333      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    332334      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    … …  
    419421         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    420422      END_3D 
    421       DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     423      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    422424         zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
    423425      END_3D 
    424       DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     426      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    425427         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    426428      END_3D 
    … …  
    537539         ! 
    538540         ! Neumann condition at k=2 
    539          zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag 
    540          zd_lw(:,:,2) = 0._wp 
     541         DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo 
     542            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 
     543            zd_lw(ji,jj,2) = 0._wp 
     544         END_2D 
    541545         ! 
    542546         ! Set psi vertical flux at the surface: 
    … …  
    613617         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    614618      END_3D 
    615       DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     619      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    616620         zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 
    617621      END_3D 
    618       DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     622      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    619623         psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    620624      END_3D 
    … …  
    811815         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri 
    812816         WRITE(numout,*) 
    813          WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
    814          WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top 
    815          WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot 
    816          WRITE(numout,*) 
    817817      ENDIF 
    818818 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfphy.F90

    r13472 r13886  
    337337      ! 
    338338   END SUBROUTINE zdf_phy 
     339 
     340 
    339341   INTEGER FUNCTION zdf_phy_alloc() 
    340342      !!---------------------------------------------------------------------- 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdftke.F90

    r13530 r13886  
    678678            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
    679679         END SELECT       
    680          IF( .NOT.ln_drg_OFF ) THEN 
    681             WRITE(numout,*) 
    682             WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
    683             WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top 
    684             WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot 
    685          ENDIF 
    686680         WRITE(numout,*) 
    687681         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/nemogcm.F90