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 12930 for NEMO/branches/2020 – NEMO

Changeset 12930 for NEMO/branches/2020


Ignore:
Timestamp:
2020-05-15T09:51:24+02:00 (4 years ago)
Author:
smasson
Message:

r12581_ticket2418: update with trunk @12929, see #2418

Location:
NEMO/branches/2020/r12581_ticket2418
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/r12581_ticket2418

    • Property svn:externals

      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@12798        sette 

  • NEMO/branches/2020/r12581_ticket2418/cfgs/WED025/EXPREF/file_def_nemo-ice.xml

    r11844 r12930  
    7878     </file> 
    7979      
    80      <file id="file22" name_suffix="_SBC_scalar" description="scalar variables" enabled=".true." > 
    81        <!-- global contents --> 
    82        <field field_ref="ibgvol_tot"     grid_ref="grid_1point"   name="ibgvol_tot"   /> 
    83        <field field_ref="sbgvol_tot"     grid_ref="grid_1point"   name="sbgvol_tot"   /> 
    84        <field field_ref="ibgarea_tot"    grid_ref="grid_1point"   name="ibgarea_tot"  /> 
    85        <field field_ref="ibgsalt_tot"    grid_ref="grid_1point"   name="ibgsalt_tot"  /> 
    86        <field field_ref="ibgheat_tot"    grid_ref="grid_1point"   name="ibgheat_tot"  /> 
    87        <field field_ref="sbgheat_tot"    grid_ref="grid_1point"   name="sbgheat_tot"  /> 
    88         
    89        <!-- global drifts (conservation checks) --> 
    90        <field field_ref="ibgvolume"      grid_ref="grid_1point"   name="ibgvolume"    /> 
    91        <field field_ref="ibgsaltco"      grid_ref="grid_1point"   name="ibgsaltco"    /> 
    92        <field field_ref="ibgheatco"      grid_ref="grid_1point"   name="ibgheatco"    /> 
    93        <field field_ref="ibgheatfx"      grid_ref="grid_1point"   name="ibgheatfx"    /> 
    94         
    95        <!-- global forcings  --> 
    96        <field field_ref="ibgfrcvoltop"   grid_ref="grid_1point"   name="ibgfrcvoltop" /> 
    97        <field field_ref="ibgfrcvolbot"   grid_ref="grid_1point"   name="ibgfrcvolbot" /> 
    98        <field field_ref="ibgfrctemtop"   grid_ref="grid_1point"   name="ibgfrctemtop" /> 
    99        <field field_ref="ibgfrctembot"   grid_ref="grid_1point"   name="ibgfrctembot" /> 
    100        <field field_ref="ibgfrcsal"      grid_ref="grid_1point"   name="ibgfrcsal"    /> 
    101        <field field_ref="ibgfrchfxtop"   grid_ref="grid_1point"   name="ibgfrchfxtop" /> 
    102        <field field_ref="ibgfrchfxbot"   grid_ref="grid_1point"   name="ibgfrchfxbot" /> 
    103      </file> 
    104       
    10580   </file_group> 
    10681    

  • NEMO/branches/2020/r12581_ticket2418/cfgs/WED025/EXPREF/namelist_cfg

    r12835 r12930  
    55!! namelists    2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 
    66!!                                    namsbc_sas, namtra_qsr, namsbc_rnf, 
    7 !!                                    namsbc_isf, namsbc_iscpl, namsbc_apr,  
     7!!                                    namisf, namsbc_apr,  
    88!!                                    namsbc_ssr, namsbc_wave, namberg) 
    99!!              3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) 
    … …  
    3838   nn_it000    =   1       !  first time step 
    3939   nn_itend    =  26280    !  last  time step (std 5475) 
    40    nn_date0    = 19760301  !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
     40   nn_date0    = 20000101  !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    4141   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    4242      nn_rstctl   =    2      !  restart control ==> activated only if ln_rstart=T 
    … …  
    6161   ln_tsd_init = .true.          !  ocean initialisation 
    6262   ln_tsd_dmp  = .false.         !  T-S restoring   (see namtra_dmp) 
    63     
     63 
    6464   cn_dir      = './'      !  root directory for the T-S data location 
    65    !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
    66    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    67    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    68    sn_tem = 'dta_temp_WED025'            ,         -12       , 'votemper',   .true.    , .true. , 'yearly'  ,    ''            ,    ''    ,    '' 
    69    sn_sal = 'dta_sal_WED025'             ,         -12       , 'vosaline',   .true.    , .true. , 'yearly'  ,    ''            ,    ''    ,    '' 
     65   !___________!_____________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 
     66   !           !  file name          ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
     67   !           !                     !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
     68   sn_tem = 'WED025_init_JRA_200001.nc',       -12       , 'votemper',   .false.   , .true. , 'yearly'  ,    ''            ,    ''    ,    '' 
     69   sn_sal = 'WED025_init_JRA_200001.nc',       -12       , 'vosaline',   .false.   , .true. , 'yearly'  ,    ''            ,    ''    ,    '' 
    7070/ 
    7171!----------------------------------------------------------------------- 
    … …  
    124124                     ! Misc. options of sbc :  
    125125   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr) 
    126    ln_dm2dc    = .true.    !  daily mean to diurnal cycle on short wave 
     126   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    127127   ln_ssr      = .false.   !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    128128   nn_fwb      = 0         !  FreshWater Budget: =0 unchecked 
    … …  
    141141   ln_NCAR     = .true.   ! "NCAR"      algorithm   (Large and Yeager 2008) 
    142142   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    143    ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    144    ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
     143   ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
     144   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 45r1) 
    145145 
    146146   cn_dir      = './'      !  root directory for the bulk data location 
    … …  
    148148   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !       weights filename               ! rotation ! land/sea mask ! 
    149149   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   ! 
    150    sn_wndi     = 'u10_core'              ,         6         , 'U_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bicubic_core.nc'           , 'Uwnd'   , '' 
    151    sn_wndj     = 'v10_core'              ,         6         , 'V_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bicubic_core.nc'           , 'Vwnd'   , '' 
    152    sn_qsr      = 'qsw_core'              ,        24         , 'SWDN_MOD',   .false.   , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    153    sn_qlw      = 'qlw_core'              ,        24         , 'LWDN_MOD',   .false.   , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    154    sn_tair     = 't10_core'              ,         6         , 'T_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    155    sn_humi     = 'q10_core'              ,         6         , 'Q_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    156    sn_prec     = 'precip_core'           ,        -1         , 'TPRECIP',    .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    157    sn_snow     = 'snow_core'             ,        -1         , 'SNOW'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
    158    sn_slp      = 'slp_core'              ,         6         , 'SLP'     ,   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , '' 
     150   sn_wndi     = 'u10_JRA'              ,         3         , 'uas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bicubic_JRA.nc'           , 'Uwnd'   , '' 
     151   sn_wndj     = 'v10_JRA'              ,         3         , 'vas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bicubic_JRA.nc'           , 'Vwnd'   , '' 
     152   sn_qsr      = 'rsds_JRA'             ,         3         , 'rsds'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     153   sn_qlw      = 'rlds_JRA'             ,         3         , 'rlds'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     154   sn_tair     = 't10_JRA'              ,         3         , 'tas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     155   sn_humi     = 'q10_JRA'              ,         3         , 'huss_10m',   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     156   sn_prec     = 'precip_JRA'           ,         3         , 'prto'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     157   sn_snow     = 'snow_JRA'             ,         3         , 'prsn'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
     158   sn_slp      = 'slp_JRA'              ,         3         , 'psl'     ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , '' 
    159159/ 
    160160!----------------------------------------------------------------------- 
    … …  
    201201   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    202202   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    203    sn_rnf      = 'runoff_WED025'         ,  -1               , 'runoff'  ,   .true.    , .false., 'yearly'  , ''               , ''       , '' 
     203   sn_rnf      = 'WED025_icb'         ,  -1               , 'runoff'  ,   .true.    , .false., 'yearly'  , ''               , ''       , '' 
    204204/ 
    205205!----------------------------------------------------------------------- 
    … …  
    221221         cn_isfcav_mlt = '3eq'   ! ice shelf melting formulation (spe/2eq/3eq/oasis) 
    222222         !                       ! spe = fwfisf is read from a forcing field 
    223          !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006) 
    224          !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 
     223         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006 for a short description) 
     224         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2016 for a short description) 
    225225         !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf 
    226226         !              !  cn_isfcav_mlt = 2eq or 3eq cases: 
    227227         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99) 
    228          !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    229          !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999) 
    230          rn_gammat0  = 1.4e-2    ! gammat coefficient used in blk formula 
    231          rn_gammas0  = 4.e-4    ! gammas coefficient used in blk formula 
     228         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0) 
     229         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description) 
     230         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description) 
     231         rn_gammat0  = 1.4e-2    ! gammat coefficient used in spe, vel and vel_stab gamma computation method 
     232         rn_gammas0  = 4.0e-4    ! gammas coefficient used in spe, vel and vel_stab gamma computation method 
    232233         ! 
    233234         rn_htbl     =  30.      ! thickness of the top boundary layer    (Losh et al. 2008) 
    … …  
    255256         sn_isfpar_zmin = 'isfmlt_par',      -12.      , 'sozisfmin' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    256257         !* 'spe' and 'oasis' case 
    257          sn_isfpar_fwf = 'isfmlt_par' ,      -12.      , 'sofwfisf' ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
     258         sn_isfpar_fwf = 'isfmlt_par' ,      -12.      ,'sofwfisf' ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
    258259         !* 'bg03' case 
    259          sn_isfpar_Leff = 'isfmlt_par',       0.       , 'Leff'     ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
     260         sn_isfpar_Leff = 'isfmlt_par',       0.       ,'Leff'     ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
    260261      ! 
    261262      ! ---------------- ice sheet coupling ------------------------------- 
    … …  
    300301   ln_tide     = .true.       ! Activate tides 
    301302      ln_tide_pot   = .false.               !  use tidal potential forcing 
    302       clname(1) = 'M2'  !  name of constituent - all tidal components must be set in namelist_cfg 
    303       clname(2) = 'S2' 
    304       clname(3) = 'K1' 
    305       clname(4) = 'O1' 
     303      sn_tide_cnames(1) = 'M2'  !  name of constituent - all tidal components must be set in namelist_cfg 
     304      sn_tide_cnames(2) = 'S2' 
     305      sn_tide_cnames(3) = 'K1' 
     306      sn_tide_cnames(4) = 'O1' 
    306307/ 
    307308!----------------------------------------------------------------------- 
    … …  
    340341   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    341342   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    342    bn_ssh      =    'bdyT_ssh_WED025'    ,         -1        , 'sossheig' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    343    bn_u2d      =    'bdyU_u2d_WED025'    ,         -1        , 'vobtcrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    344    bn_v2d      =    'bdyV_u2d_WED025'    ,         -1        , 'vobtcrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    345    bn_u3d      =    'bdyU_u3d_WED025'    ,         -1        , 'vozocrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    346    bn_v3d      =    'bdyV_u3d_WED025'    ,         -1        , 'vomecrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    347    bn_tem      =    'bdyT_tra_WED025'    ,         -1        , 'votemper' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    348    bn_sal      =    'bdyT_tra_WED025'    ,         -1        , 'vosaline' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     343   bn_ssh      =    'WED025_bdyT_ssh'    ,         -1        , 'sossheig' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     344   bn_u2d      =    'WED025_bdyU_u2d'    ,         -1        , 'vobtcrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     345   bn_v2d      =    'WED025_bdyV_u2d'    ,         -1        , 'vobtcrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     346   bn_u3d      =    'WED025_bdyU_u3d'    ,         -1        , 'vozocrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     347   bn_v3d      =    'WED025_bdyV_u3d'    ,         -1        , 'vomecrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     348   bn_tem      =    'WED025_bdyT_tra'    ,         -1        , 'votemper' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     349   bn_sal      =    'WED025_bdyT_tra'    ,         -1        , 'vosaline' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    349350!* for si3 
    350    bn_a_i      =    'bdyT_ice_WED025'    ,         -1        , 'ileadfra' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    351    bn_h_i      =    'bdyT_ice_WED025'    ,         -1        , 'iicethic' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    352    bn_h_s      =    'bdyT_ice_WED025'    ,         -1        , 'isnowthi' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     351   bn_a_i      =    'WED025_bdyT_ice'    ,         -1        , 'ileadfra' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     352   bn_h_i      =    'WED025_bdyT_ice'    ,         -1        , 'iicethic' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
     353   bn_h_s      =    'WED025_bdyT_ice'    ,         -1        , 'isnowthi' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     '' 
    353354/ 
    354355!----------------------------------------------------------------------- 
    355356&nambdy_tide   !  tidal forcing at open boundaries                      (default: OFF) 
    356357!----------------------------------------------------------------------- 
    357    filtide          = 'bdytide_WED025_'         !  file name root of tidal forcing files 
     358   filtide          = 'WED025_bdytide_'         !  file name root of tidal forcing files 
    358359/ 
    359360 
    … …  
    658659&namctl        !   Control prints                                       (default: OFF) 
    659660!----------------------------------------------------------------------- 
    660    sn_cfctl%l_runstat = .FALSE.   ! switches and which areas produce reports with the proc integer settings. 
     661   sn_cfctl%l_runstat = .true.    ! switches and which areas produce reports with the proc integer settings. 
    661662   ln_timing   = .true.           !  timing by routine write out in timing.output file 
    662663/ 

  • NEMO/branches/2020/r12581_ticket2418/cfgs/WED025/EXPREF/namelist_ice_cfg

    r11487 r12930  
    2626&namitd         !   Ice discretization 
    2727!------------------------------------------------------------------------------ 
     28   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05) 
     29      rn_himean     =   2.0           !  expected domain-average ice thickness (m) 
     30   rn_himin         =   0.01          !  minimum ice thickness (m) used in remapping 
    2831/ 
    2932!------------------------------------------------------------------------------ 
    3033&namdyn         !   Ice dynamics 
    3134!------------------------------------------------------------------------------ 
     35   ln_landfast_L16  = .true.          !  landfast: parameterization from Lemieux 2016 
    3236/ 
    3337!------------------------------------------------------------------------------ 
    … …  
    4246&namdyn_adv     !   Ice advection 
    4347!------------------------------------------------------------------------------ 
     48   ln_adv_Pra       = .false.         !  Advection scheme (Prather) 
     49   ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho) 
     50      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order) 
    4451/ 
    4552!------------------------------------------------------------------------------ 
    … …  
    6269&namthd_do      !   Ice growth in open water 
    6370!------------------------------------------------------------------------------ 
     71   rn_hinew         =   0.02          !  thickness for new ice formation in open water (m), must be larger than rn_himin 
     72   ln_frazil        = .true.          !  Frazil ice parameterization (ice collection as a function of wind) 
    6473/ 
    6574!------------------------------------------------------------------------------ 
    … …  
    7079&namthd_pnd     !   Melt ponds 
    7180!------------------------------------------------------------------------------ 
     81   ln_pnd           = .true.          !  activate melt ponds or not 
     82     ln_pnd_H12     = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
     83     ln_pnd_alb     = .true.          !  melt ponds affect albedo or not 
    7284/ 
     85 
    7386!------------------------------------------------------------------------------ 
    7487&namini         !   Ice initialization 
    7588!------------------------------------------------------------------------------ 
     89   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
     90   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F) 
     91   ! -- for ln_iceini_file = T 
     92   sn_hti = 'WED025_init_JRA_200001.nc', -12 ,'icethic_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
     93   sn_hts = 'WED025_init_JRA_200001.nc', -12 ,'icesnow_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
     94   sn_ati = 'WED025_init_JRA_200001.nc', -12 ,'ice_cover'  ,  .false.  , .true., 'yearly'  , '' , '', '' 
     95   sn_smi = 'NOT USED'              , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     96   sn_tmi = 'NOT USED'              , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     97   sn_tsu = 'NOT USED'              , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     98   sn_tms = 'NOT USED'              , -12 ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     99   !      melt ponds (be careful, sn_apd is the pond concentration (not fraction), so it differs from rn_apd) 
     100   sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     101   sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     102   cn_dir='./' 
    76103/ 
    77104!------------------------------------------------------------------------------ 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/BDY/bdydta.F90

    r12655 r12930  
    9191      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices 
    9292      INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers 
    93       INTEGER,   DIMENSION(jpbgrd)     ::   ilen1  
    9493      TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut 
    9594      TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias 
    … …  
    116115                  END DO 
    117116               ENDIF 
    118                IF( dta_bdy(jbdy)%lneed_dyn2d .AND. ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer 
     117               IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 
    119118                  igrd = 2 
    120                   DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)   ! u2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init 
     119                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)      ! u2d is used either over the whole bdy or only on the rim 
    121120                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    122121                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    123122                     dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)          
    124123                  END DO 
     124               ENDIF 
     125               IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 
    125126                  igrd = 3 
    126                   DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)   ! v2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init 
     127                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)      ! v2d is used either over the whole bdy or only on the rim 
    127128                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    128129                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    … …  
    210211         ! 
    211212         ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s) 
    212          IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d 
     213         IF( cn_tra(jbdy) == 'runoff' ) THEN   ! runoff 
    213214            ! 
    214             igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 
    215             DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    216                ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    217                ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    218                dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    219             END DO 
    220             igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 
    221             DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    222                ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    223                ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    224                dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
    225             END DO 
     215            IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 
     216               igrd = 2                         ! zonal flow (m3/s) to barotropic zonal velocity (m/s) 
     217               DO ib = 1, SIZE(dta_alias%u2d)   ! u2d is used either over the whole bdy or only on the rim 
     218                  ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     219                  ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     220                  dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
     221               END DO 
     222            ENDIF 
     223            IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain 
     224               igrd = 3                         ! meridional flow (m3/s) to barotropic meridional velocity (m/s) 
     225               DO ib = 1, SIZE(dta_alias%v2d)   ! v2d is used either over the whole bdy or only on the rim 
     226                  ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     227                  ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     228                  dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
     229               END DO 
     230            ENDIF 
    226231         ENDIF 
    227232 
    228233         ! tidal harmonic forcing ONLY: initialise arrays 
    229234         IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d  
    230             IF( dta_alias%lneed_ssh   .AND. ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp 
    231             IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp 
    232             IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp 
     235            IF( ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp 
     236            IF( ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp 
     237            IF( ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp 
    233238         ENDIF 
    234239 
    … …  
    330335            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop 
    331336               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 
    332                   IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=idx_bdy(jbdy)%nblen(:) 
    333                   ELSE                                 ;   ilen1(:)=idx_bdy(jbdy)%nblenrim(:) 
    334                   ENDIF 
    335                   IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
    336                   IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
    337                   IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
     337                  IF( ASSOCIATED(dta_bdy(jbdy)%ssh) ) dta_bdy_s(jbdy)%ssh(:) = dta_bdy(jbdy)%ssh(:) 
     338                  IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) dta_bdy_s(jbdy)%u2d(:) = dta_bdy(jbdy)%u2d(:) 
     339                  IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) dta_bdy_s(jbdy)%v2d(:) = dta_bdy(jbdy)%v2d(:) 
    338340               ENDIF 
    339341            END DO 
    340342         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
    341343            ! 
    342             ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step 
    343344            CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp ) 
    344345         ENDIF 
    … …  
    348349      ! 
    349350   END SUBROUTINE bdy_dta 
    350  
     351    
    351352 
    352353   SUBROUTINE bdy_dta_init 
    … …  
    380381      LOGICAL                                ::   llneed        ! 
    381382      LOGICAL                                ::   llread        ! 
     383      LOGICAL                                ::   llfullbdy     ! 
    382384      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
    383385      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    … …  
    494496               igrd = 2                                                    ! U point 
    495497               ipk = 1                                                     ! surface data 
    496                llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed 
     498               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%u2d will be needed 
    497499               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file 
    498500               bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy 
    499501               bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta 
    500                IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy 
    501                ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim 
     502               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need u2d over the whole bdy or only over the rim? 
     503               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd) 
     504               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd) 
    502505               ENDIF 
    503506            ENDIF 
    … …  
    506509               igrd = 3                                                    ! V point 
    507510               ipk = 1                                                     ! surface data 
    508                llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed 
     511               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%v2d will be needed 
    509512               llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file 
    510513               bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy 
    511514               bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta  
    512                IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy 
    513                ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim 
     515               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need v2d over the whole bdy or only over the rim? 
     516               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd) 
     517               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd) 
    514518               ENDIF 
    515519            ENDIF 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/BDY/bdyini.F90

    r12377 r12930  
    1919   USE oce            ! ocean dynamics and tracers variables 
    2020   USE dom_oce        ! ocean space and time domain 
     21   USE sbc_oce , ONLY: nn_ice 
    2122   USE bdy_oce        ! unstructured open boundary conditions 
    2223   USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine) 
    2324   USE bdytides       ! open boundary cond. setting   (bdytide_init routine) 
    2425   USE tide_mod, ONLY: ln_tide ! tidal forcing 
    25    USE phycst   , ONLY: rday 
     26   USE phycst  , ONLY: rday 
    2627   ! 
    2728   USE in_out_manager ! I/O units 
    … …  
    315316 
    316317         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 
     318 
     319         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN 
     320            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice 
     321            CALL ctl_stop( ctmp1 ) 
     322         ENDIF 
    317323 
    318324         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN  

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/BDY/bdytides.F90

    r12489 r12930  
    6565      !! namelist variables 
    6666      !!------------------- 
    67       CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
    68       LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
     67      CHARACTER(len=80)                         ::   filtide             ! Filename root for tidal input files 
     68      LOGICAL                                   ::   ln_bdytide_2ddta    ! If true, read 2d harmonic data 
    6969      !! 
    70       INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
    71       INTEGER                                   ::   ii, ij              !: dummy loop indices 
     70      INTEGER                                   ::   ib_bdy, itide, ib   ! dummy loop indices 
     71      INTEGER                                   ::   ii, ij              ! dummy loop indices 
    7272      INTEGER                                   ::   inum, igrd 
    73       INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
     73      INTEGER                                   ::   isz                 ! bdy data size 
    7474      INTEGER                                   ::   ios                 ! Local integer output status for namelist read 
    7575      INTEGER                                   ::   nbdy_rdstart, nbdy_loc 
    76       CHARACTER(LEN=50)                         ::   cerrmsg             !: error string 
    77       CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
    78       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data 
    79       REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      "  
     76      CHARACTER(LEN=50)                         ::   cerrmsg             ! error string 
     77      CHARACTER(len=80)                         ::   clfile              ! full file name for tidal input file  
     78      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            ! work space to read in tidal harmonics data 
     79      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !  "     "    "   "   "   "        "      "  
    8080      !! 
    81       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
     81      TYPE(TIDES_DATA), POINTER                 ::   td                  ! local short cut    
     82      TYPE(  OBC_DATA), POINTER                 ::   dta                 ! local short cut 
    8283      !! 
    8384      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta 
    … …  
    9394         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN 
    9495            ! 
    95             td => tides(ib_bdy) 
    96  
     96            td  => tides(ib_bdy) 
     97            dta => dta_bdy(ib_bdy) 
     98          
    9799            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
    98100            filtide(:) = '' 
    … …  
    130132            IF(lwp) WRITE(numout,*) ' ' 
    131133 
    132             ! Allocate space for tidal harmonics data - get size from OBC data arrays 
     134            ! Allocate space for tidal harmonics data - get size from BDY data arrays 
     135            ! Allocate also slow varying data in the case of time splitting: 
     136            ! Do it anyway because at this stage knowledge of free surface scheme is unknown 
    133137            ! ----------------------------------------------------------------------- 
    134  
    135             ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
    136             ! relaxation area       
    137             IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:) 
    138             ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:) 
    139             ENDIF 
    140  
    141             ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 
    142             ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 
    143  
    144             ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 
    145             ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 
    146  
    147             ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 
    148             ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
    149  
    150             td%ssh0(:,:,:) = 0._wp 
    151             td%ssh (:,:,:) = 0._wp 
    152             td%u0  (:,:,:) = 0._wp 
    153             td%u   (:,:,:) = 0._wp 
    154             td%v0  (:,:,:) = 0._wp 
    155             td%v   (:,:,:) = 0._wp 
    156  
     138            IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain 
     139               isz = SIZE(dta%ssh) 
     140               ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) ) 
     141               dta_bdy_s(ib_bdy)%ssh(:) = 0._wp   ! needed? 
     142            ENDIF 
     143            IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain 
     144               isz = SIZE(dta%u2d) 
     145               ALLOCATE( td%u0  ( isz, nb_harmo, 2 ), td%u  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) ) 
     146               dta_bdy_s(ib_bdy)%u2d(:) = 0._wp   ! needed? 
     147            ENDIF 
     148            IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain 
     149               isz = SIZE(dta%v2d) 
     150               ALLOCATE( td%v0  ( isz, nb_harmo, 2 ), td%v  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) ) 
     151               dta_bdy_s(ib_bdy)%v2d(:) = 0._wp   ! needed? 
     152            ENDIF 
     153 
     154            ! fill td%ssh0, td%u0, td%v0 
     155            ! ----------------------------------------------------------------------- 
    157156            IF( ln_bdytide_2ddta ) THEN 
     157               ! 
    158158               ! It is assumed that each data file contains all complex harmonic amplitudes 
    159159               ! given on the global domain (ie global, jpiglo x jpjglo) 
    … …  
    162162               ! 
    163163               ! SSH fields 
    164                clfile = TRIM(filtide)//'_grid_T.nc' 
    165                CALL iom_open( clfile , inum )  
    166                igrd = 1                       ! Everything is at T-points here 
    167                DO itide = 1, nb_harmo 
    168                   CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 
    169                   CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) )  
    170                   DO ib = 1, ilen0(igrd) 
    171                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    172                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    173                      IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? 
    174                      td%ssh0(ib,itide,1) = ztr(ii,ij) 
    175                      td%ssh0(ib,itide,2) = zti(ii,ij) 
    176                   END DO 
    177                END DO  
    178                CALL iom_close( inum ) 
     164               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain 
     165                  clfile = TRIM(filtide)//'_grid_T.nc' 
     166                  CALL iom_open( clfile , inum )  
     167                  igrd = 1                       ! Everything is at T-points here 
     168                  DO itide = 1, nb_harmo 
     169                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 
     170                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) )  
     171                     DO ib = 1, SIZE(dta%ssh) 
     172                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     173                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     174                        td%ssh0(ib,itide,1) = ztr(ii,ij) 
     175                        td%ssh0(ib,itide,2) = zti(ii,ij) 
     176                     END DO 
     177                  END DO 
     178                  CALL iom_close( inum ) 
     179               ENDIF 
    179180               ! 
    180181               ! U fields 
    181                clfile = TRIM(filtide)//'_grid_U.nc' 
    182                CALL iom_open( clfile , inum )  
    183                igrd = 2                       ! Everything is at U-points here 
    184                DO itide = 1, nb_harmo 
    185                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 
    186                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 
    187                   DO ib = 1, ilen0(igrd) 
    188                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    189                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    190                      IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? 
    191                      td%u0(ib,itide,1) = ztr(ii,ij) 
    192                      td%u0(ib,itide,2) = zti(ii,ij) 
    193                   END DO 
    194                END DO 
    195                CALL iom_close( inum ) 
     182               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain 
     183                  clfile = TRIM(filtide)//'_grid_U.nc' 
     184                  CALL iom_open( clfile , inum )  
     185                  igrd = 2                       ! Everything is at U-points here 
     186                  DO itide = 1, nb_harmo 
     187                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 
     188                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 
     189                     DO ib = 1, SIZE(dta%u2d) 
     190                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     191                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     192                        td%u0(ib,itide,1) = ztr(ii,ij) 
     193                        td%u0(ib,itide,2) = zti(ii,ij) 
     194                     END DO 
     195                  END DO 
     196                  CALL iom_close( inum ) 
     197               ENDIF 
    196198               ! 
    197199               ! V fields 
    198                clfile = TRIM(filtide)//'_grid_V.nc' 
    199                CALL iom_open( clfile , inum )  
    200                igrd = 3                       ! Everything is at V-points here 
    201                DO itide = 1, nb_harmo 
    202                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 
    203                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 
    204                   DO ib = 1, ilen0(igrd) 
    205                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    206                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    207                      IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove? 
    208                      td%v0(ib,itide,1) = ztr(ii,ij) 
    209                      td%v0(ib,itide,2) = zti(ii,ij) 
    210                   END DO 
    211                END DO   
    212                CALL iom_close( inum ) 
     200               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain 
     201                  clfile = TRIM(filtide)//'_grid_V.nc' 
     202                  CALL iom_open( clfile , inum )  
     203                  igrd = 3                       ! Everything is at V-points here 
     204                  DO itide = 1, nb_harmo 
     205                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 
     206                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 
     207                     DO ib = 1, SIZE(dta%v2d) 
     208                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     209                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     210                        td%v0(ib,itide,1) = ztr(ii,ij) 
     211                        td%v0(ib,itide,2) = zti(ii,ij) 
     212                     END DO 
     213                  END DO 
     214                  CALL iom_close( inum ) 
     215               ENDIF 
    213216               ! 
    214217               DEALLOCATE( ztr, zti )  
    … …  
    218221               ! Read tidal data only on bdy segments 
    219222               !  
    220                ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 
     223               ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) ) 
    221224               ! 
    222225               ! Open files and read in tidal forcing data 
    … …  
    225228               DO itide = 1, nb_harmo 
    226229                  !                                                              ! SSH fields 
    227                   clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 
    228                   CALL iom_open( clfile, inum ) 
    229                   CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    230                   td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
    231                   CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    232                   td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
    233                   CALL iom_close( inum ) 
     230                  IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain 
     231                     isz = SIZE(dta%ssh) 
     232                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 
     233                     CALL iom_open( clfile, inum ) 
     234                     CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     235                     td%ssh0(:,itide,1) = dta_read(1:isz,1,1) 
     236                     CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     237                     td%ssh0(:,itide,2) = dta_read(1:isz,1,1) 
     238                     CALL iom_close( inum ) 
     239                  ENDIF 
    234240                  !                                                              ! U fields 
    235                   clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 
    236                   CALL iom_open( clfile, inum ) 
    237                   CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    238                   td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
    239                   CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    240                   td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
    241                   CALL iom_close( inum ) 
     241                  IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain 
     242                     isz = SIZE(dta%u2d) 
     243                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 
     244                     CALL iom_open( clfile, inum ) 
     245                     CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     246                     td%u0(:,itide,1) = dta_read(1:isz,1,1) 
     247                     CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     248                     td%u0(:,itide,2) = dta_read(1:isz,1,1) 
     249                     CALL iom_close( inum ) 
     250                  ENDIF 
    242251                  !                                                              ! V fields 
    243                   clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 
    244                   CALL iom_open( clfile, inum ) 
    245                   CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    246                   td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
    247                   CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    248                   td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
    249                   CALL iom_close( inum ) 
     252                  IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain 
     253                     isz = SIZE(dta%v2d) 
     254                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 
     255                     CALL iom_open( clfile, inum ) 
     256                     CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     257                     td%v0(:,itide,1) = dta_read(1:isz,1,1) 
     258                     CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     259                     td%v0(:,itide,2) = dta_read(1:isz,1,1) 
     260                     CALL iom_close( inum ) 
     261                  ENDIF 
    250262                  ! 
    251263               END DO ! end loop on tidal components 
    … …  
    254266               ! 
    255267            ENDIF ! ln_bdytide_2ddta=.true. 
    256             ! 
    257             ! Allocate slow varying data in the case of time splitting: 
    258             ! Do it anyway because at this stage knowledge of free surface scheme is unknown 
    259             ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
    260             ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
    261             ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
    262             dta_bdy_s(ib_bdy)%ssh(:) = 0._wp 
    263             dta_bdy_s(ib_bdy)%u2d(:) = 0._wp 
    264             dta_bdy_s(ib_bdy)%v2d(:) = 0._wp 
    265268            ! 
    266269         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2 
    … …  
    283286      ! 
    284287      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step 
    285       INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices 
    286       INTEGER, DIMENSION(jpbgrd)   ::   ilen0  
    287       INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts 
     288      INTEGER  ::   itide, ib_bdy, ib         ! loop indices 
    288289      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset    
    289290      !!---------------------------------------------------------------------- 
    … …  
    310311         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN 
    311312            ! 
    312             nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd) 
    313             nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd) 
    314             ! 
    315             IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:) 
    316             ELSE                                   ;   ilen0(:) = nblenrim(:) 
    317             ENDIF      
    318             ! 
    319313            ! We refresh nodal factors every day below 
    320314            ! This should be done somewhere else 
    … …  
    337331            ! If time splitting, initialize arrays from slow varying open boundary data: 
    338332            IF ( PRESENT(kit) ) THEN            
    339                IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
    340                IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
    341                IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
     333               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:) 
     334               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:) 
     335               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:) 
    342336            ENDIF 
    343337            ! 
    … …  
    349343               z_sist = zramp * SIN( z_sarg ) 
    350344               ! 
    351                IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN 
    352                   igrd=1                              ! SSH on tracer grid 
    353                   DO ib = 1, ilen0(igrd) 
     345               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN   ! SSH on tracer grid 
     346                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh) 
    354347                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 
    355348                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & 
    … …  
    358351               ENDIF 
    359352               ! 
    360                IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 
    361                   igrd=2                              ! U grid 
    362                   DO ib = 1, ilen0(igrd) 
     353               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN  ! U grid 
     354                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d) 
    363355                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 
    364356                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 
    365357                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist ) 
    366358                  END DO 
    367                   igrd=3                              ! V grid 
    368                   DO ib = 1, ilen0(igrd)  
     359               ENDIF 
     360               ! 
     361               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN   ! V grid 
     362                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d) 
    369363                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 
    370364                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & 
    … …  
    372366                  END DO 
    373367               ENDIF 
     368               ! 
    374369            END DO              
    375          END IF 
     370         ENDIF 
    376371      END DO 
    377372      ! 
    … …  
    386381      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data 
    387382      ! 
    388       INTEGER ::   itide, igrd, ib       ! dummy loop indices 
    389       INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays) 
     383      INTEGER ::   itide, isz, ib       ! dummy loop indices 
    390384      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
    391385      !!---------------------------------------------------------------------- 
    392386      ! 
    393       igrd=1    
    394                               ! SSH on tracer grid. 
    395       ilen0(1) =  SIZE(td%ssh0(:,1,1)) 
    396       ! 
    397       ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) ) 
    398       ! 
    399       DO itide = 1, nb_harmo 
    400          DO ib = 1, ilen0(igrd) 
    401             mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 
    402             phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 
     387      IF( ASSOCIATED(td%ssh0) ) THEN   ! SSH on tracer grid. 
     388         ! 
     389         isz = SIZE( td%ssh0, dim = 1 ) 
     390         ALLOCATE( mod_tide(isz), phi_tide(isz) ) 
     391         ! 
     392         DO itide = 1, nb_harmo 
     393            DO ib = 1, isz 
     394               mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) ) 
     395               phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 
     396            END DO 
     397            DO ib = 1, isz 
     398               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     399               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 
     400            END DO 
     401            DO ib = 1, isz 
     402               td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     403               td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     404            END DO 
    403405         END DO 
    404          DO ib = 1 , ilen0(igrd) 
    405             mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
    406             phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 
    407          ENDDO 
    408          DO ib = 1 , ilen0(igrd) 
    409             td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    410             td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
    411          ENDDO 
    412       END DO 
    413       ! 
    414       DEALLOCATE( mod_tide, phi_tide ) 
     406         ! 
     407         DEALLOCATE( mod_tide, phi_tide ) 
     408         ! 
     409      ENDIF 
    415410      ! 
    416411   END SUBROUTINE tide_init_elevation 
    … …  
    424419      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data 
    425420      ! 
    426       INTEGER ::   itide, igrd, ib       ! dummy loop indices 
    427       INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays) 
     421      INTEGER ::   itide, isz, ib        ! dummy loop indices 
    428422      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
    429423      !!---------------------------------------------------------------------- 
    430424      ! 
    431       ilen0(2) =  SIZE(td%u0(:,1,1)) 
    432       ilen0(3) =  SIZE(td%v0(:,1,1)) 
    433       ! 
    434       igrd=2                                 ! U grid. 
    435       ! 
    436       ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) ) 
    437       ! 
    438       DO itide = 1, nb_harmo 
    439          DO ib = 1, ilen0(igrd) 
    440             mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 
    441             phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 
     425      IF( ASSOCIATED(td%u0) ) THEN   ! U grid. we use bdy u2d on this mpi subdomain 
     426         ! 
     427         isz = SIZE( td%u0, dim = 1 ) 
     428         ALLOCATE( mod_tide(isz), phi_tide(isz) ) 
     429         ! 
     430         DO itide = 1, nb_harmo 
     431            DO ib = 1, isz 
     432               mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) ) 
     433               phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 
     434            END DO 
     435            DO ib = 1, isz 
     436               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     437               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
     438            END DO 
     439            DO ib = 1, isz 
     440               td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     441               td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     442            END DO 
    442443         END DO 
    443          DO ib = 1, ilen0(igrd) 
    444             mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
    445             phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    446          ENDDO 
    447          DO ib = 1, ilen0(igrd) 
    448             td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    449             td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
    450          ENDDO 
    451       END DO 
    452       ! 
    453       DEALLOCATE( mod_tide , phi_tide ) 
    454       ! 
    455       igrd=3                                 ! V grid. 
    456       ! 
    457       ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) ) 
    458  
    459       DO itide = 1, nb_harmo 
    460          DO ib = 1, ilen0(igrd) 
    461             mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 
    462             phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 
     444         ! 
     445         DEALLOCATE( mod_tide, phi_tide ) 
     446         ! 
     447      ENDIF 
     448      ! 
     449      IF( ASSOCIATED(td%v0) ) THEN   ! V grid. we use bdy u2d on this mpi subdomain 
     450         ! 
     451         isz = SIZE( td%v0, dim = 1 ) 
     452         ALLOCATE( mod_tide(isz), phi_tide(isz) ) 
     453         ! 
     454         DO itide = 1, nb_harmo 
     455            DO ib = 1, isz 
     456               mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) ) 
     457               phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 
     458            END DO 
     459            DO ib = 1, isz 
     460               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     461               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
     462            END DO 
     463            DO ib = 1, isz 
     464               td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
     465               td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     466            END DO 
    463467         END DO 
    464          DO ib = 1, ilen0(igrd) 
    465             mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
    466             phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    467          ENDDO 
    468          DO ib = 1, ilen0(igrd) 
    469             td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    470             td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
    471          ENDDO 
    472       END DO 
    473       ! 
    474       DEALLOCATE( mod_tide, phi_tide ) 
    475       ! 
    476   END SUBROUTINE tide_init_velocities 
     468         ! 
     469         DEALLOCATE( mod_tide, phi_tide ) 
     470         ! 
     471      ENDIF 
     472      ! 
     473   END SUBROUTINE tide_init_velocities 
    477474 
    478475   !!====================================================================== 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/ISF/isfdiags.F90

    r12340 r12930  
    8888      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac  ! thickness of the tbl and fraction of last cell affected by the tbl 
    8989      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d        ! 2d var to map in 3d 
    90       CHARACTER(LEN=256), INTENT(in) :: cdvar 
     90      CHARACTER(LEN=*), INTENT(in) :: cdvar 
    9191      !!--------------------------------------------------------------------- 
    9292      INTEGER  :: ji, jj, jk                       ! loop indices 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk.F90

    r12844 r12930  
    673673         ! ... utau, vtau at U- and V_points, resp. 
    674674         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    675          !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    676          DO_2D_10_10 
     675         !     Note that coastal wind stress is not used in the code... so this extra care has no effect 
     676         DO_2D_00_00 
    677677            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    678678               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    … …  
    886886 
    887887      IF( ln_blk ) THEN 
    888          ! ------------------------------------------------------------ ! 
    889          !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    890          ! ------------------------------------------------------------ ! 
    891          ! C-grid ice dynamics :   U & V-points (same as ocean) 
    892          DO_2D_00_00 
    893             putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
    894                &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
    895                &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
    896             pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
    897                &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
    898                &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     888         ! ------------------------------------------------------------- ! 
     889         !    Wind stress relative to the moving ice ( U10m - U_ice )    ! 
     890         ! ------------------------------------------------------------- ! 
     891         zztmp1 = rn_vfac * 0.5_wp 
     892         DO_2D_01_01    ! at T point  
     893            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndi(ji,jj) - zztmp1 * ( puice(ji-1,jj  ) + puice(ji,jj) ) ) 
     894            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndj(ji,jj) - zztmp1 * ( pvice(ji  ,jj-1) + pvice(ji,jj) ) ) 
     895         END_2D 
     896         ! 
     897         DO_2D_00_00    ! U & V-points (same as ocean). 
     898            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology  
     899            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) ) 
     900            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) ) 
     901            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) ) 
     902            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) ) 
    899903         END_2D 
    900904         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_phy.F90

    r12623 r12930  
    3131   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg] 
    3232   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 
    33    REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
     33   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp  !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    3434   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...) 
    3535   REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice 

  • NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbccpl.F90

    r12623 r12930  
    14811481      INTEGER ::   ji, jj   ! dummy loop indices 
    14821482      INTEGER ::   itx      ! index of taux over ice 
     1483      REAL(wp)                     ::   zztmp1, zztmp2 
    14831484      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty  
    14841485      !!---------------------------------------------------------------------- 
    … …  
    15441545            p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V) 
    15451546            p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1) 
    1546          CASE( 'F' ) 
    1547             DO_2D_00_00 
    1548                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) ) 
    1549                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) ) 
    1550             END_2D 
    15511547         CASE( 'T' ) 
    15521548            DO_2D_00_00 
    1553                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 
    1554                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 
     1549               ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology  
     1550               zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) ) 
     1551               zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) ) 
     1552               p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 
     1553               p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 
    15551554            END_2D 
    1556          CASE( 'I' ) 
    1557             DO_2D_00_00 
    1558                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) ) 
    1559                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) ) 
    1560             END_2D 
     1555            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. ) 
    15611556         END SELECT 
    1562          IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN  
    1563             CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. ) 
    1564          ENDIF 
    15651557          
    15661558      ENDIF 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/EXPREF/file_def_nemo-oce.xml

    r11889 r12930  
    2121      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
    2222 
    23       <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
    24    <file id="file1" output_freq="1mo" name_suffix="_grid_T" description="ocean T grid variables" > 
    25      <field field_ref="toce"         name="votemper"  /> 
    26      <field field_ref="soce"         name="vosaline"  /> 
    27      <field field_ref="ssh"          name="sossheig"  /> 
     23   <file id="file1" output_freq="5d" name_suffix="_grid_T" description="ocean T grid variables" > 
     24     <field field_ref="toce"         name="votemper"  operation="average" freq_op="5d" > @toce_e3t / @e3t </field> 
     25     <field field_ref="soce"         name="vosaline"  operation="average" freq_op="5d" > @soce_e3t / @e3t </field> 
     26     <field field_ref="ssh"          name="sossheig" /> 
    2827          <!-- variable for ice shelf --> 
    29           <field field_ref="fwfisf_cav"       name="sowflisf"  /> 
    30           <field field_ref="isfgammat"    name="sogammat"  /> 
    31           <field field_ref="isfgammas"    name="sogammas"  /> 
     28          <field field_ref="fwfisf_cav"  name="sowflisf"  /> 
     29          <field field_ref="isfgammat"   name="sogammat"  /> 
     30          <field field_ref="isfgammas"   name="sogammas"  /> 
    3231          <field field_ref="ttbl_cav"    name="ttbl"  /> 
    33           <field field_ref="stbl"    name="stbl"  /> 
    34           <field field_ref="utbl"    name="utbl"  /> 
    35           <field field_ref="vtbl"    name="vtbl"  /> 
     32          <field field_ref="stbl"        name="stbl"  /> 
     33          <field field_ref="utbl"        name="utbl"  /> 
     34          <field field_ref="vtbl"        name="vtbl"  /> 
    3635        </file> 
    37    <file id="file2" output_freq="1mo" name_suffix="_grid_U" description="ocean U grid variables" > 
    38           <field field_ref="uoce"         name="vozocrtx" /> 
     36   <file id="file2" output_freq="5d" name_suffix="_grid_U" description="ocean U grid variables" > 
     37          <field field_ref="uoce"         name="vozocrtx" operation="average" freq_op="5d" > @uoce_e3u / @e3u </field> /> 
    3938        </file> 
    40    <file id="file3" output_freq="1mo" name_suffix="_grid_V" description="ocean V grid variables" > 
    41           <field field_ref="voce"         name="vomecrty" />  
     39   <file id="file3" output_freq="5d" name_suffix="_grid_V" description="ocean V grid variables" > 
     40          <field field_ref="voce"         name="vomecrty" operation="average" freq_op="5d" > @voce_e3v / @e3v </field> />  
    4241        </file> 
    4342      </file_group> 
     43 
     44      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
    4445      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files --> 
    4546      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files --> 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/EXPREF/namelist_cfg

    r12489 r12930  
    114114 
    115115   ln_usr      = .true.   !  user defined formulation                  (T => check usrdef_sbc) 
    116    nn_fwb      = 1 
     116   nn_fwb      = 4 
    117117/ 
    118118!----------------------------------------------------------------------- 
    … …  
    308308&nameos        !   ocean Equation Of Seawater                           (default: NO selection) 
    309309!----------------------------------------------------------------------- 
    310    ln_teos10   = .false.         !  = Use TEOS-10 
    311    ln_eos80    = .false.         !  = Use EOS80 
    312    ln_leos     = .true.          !  = Use S-EOS (simplified Eq.) 
     310   ln_leos     = .true.          !  = Use L-EOS (linear Eq.) 
    313311                                 ! 
    314312   !                     ! S-EOS coefficients (ln_seos=T): 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/dtatsd.F90

    r12077 r12930  
    3636   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsddmp ! structure of input SST (file informations, fields read) 
    3737 
     38   !! * Substitutions 
     39#  include "do_loop_substitute.h90" 
    3840   !!---------------------------------------------------------------------- 
    3941   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6769      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
    6870      ! 
    69       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
    7071      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    7172901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' ) 
    72       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    7373      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    7474902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' ) 
    … …  
    191191         ENDIF 
    192192         ! 
    193          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    194             DO ji = 1, jpi 
    195                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
    196                   zl = gdept_0(ji,jj,jk) 
    197                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    198                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    199                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    200                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    201                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    202                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    203                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    204                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    205                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    206                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    207                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    208                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
    209                         ENDIF 
    210                      END DO 
    211                   ENDIF 
    212                END DO 
    213                DO jk = 1, jpkm1 
    214                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    215                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    216                END DO 
    217                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    218                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     193         DO_2D_11_11 
     194            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     195               zl = gdept_0(ji,jj,jk) 
     196               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
     197                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
     198                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
     199               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
     200                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
     201                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
     202               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
     203                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
     204                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
     205                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
     206                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
     207                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     208                     ENDIF 
     209                  END DO 
     210               ENDIF 
    219211            END DO 
    220          END DO 
     212            DO jk = 1, jpkm1 
     213               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
     214               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
     215            END DO 
     216            ptsd(ji,jj,jpk,jp_tem) = 0._wp 
     217            ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     218         END_2D 
    221219         !  
    222220      ELSE                                !==   z- or zps- coordinate   ==! 
    … …  
    226224         ! 
    227225         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
    228             DO jj = 1, jpj 
    229                DO ji = 1, jpi 
    230                   ik = mbkt(ji,jj)  
    231                   IF( ik > 1 ) THEN 
    232                      zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    233                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
    234                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
    235                   ENDIF 
    236                   ik = mikt(ji,jj) 
    237                   IF( ik > 1 ) THEN 
    238                      zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
    239                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
    240                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
    241                   END IF 
    242                END DO 
    243             END DO 
     226            DO_2D_11_11 
     227               ik = mbkt(ji,jj)  
     228               IF( ik > 1 ) THEN 
     229                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     230                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
     231                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
     232               ENDIF 
     233               ik = mikt(ji,jj) 
     234               IF( ik > 1 ) THEN 
     235                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
     236                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
     237                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
     238               END IF 
     239            END_2D 
    244240         ENDIF 
    245241         ! 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/eosbn2.F90

    r12489 r12930  
    180180   REAL(wp) ::   BPE002 
    181181 
     182   !! * Substitutions 
     183#  include "do_loop_substitute.h90" 
    182184   !!---------------------------------------------------------------------- 
    183185   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    241243      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    242244         ! 
    243          DO jk = 1, jpkm1 
    244             DO jj = 1, jpj 
    245                DO ji = 1, jpi 
    246                   ! 
    247                   zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    248                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    249                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    250                   ztm = tmask(ji,jj,jk)                                         ! tmask 
     245         DO_3D_11_11( 1, jpkm1 ) 
     246            ! 
     247            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     248            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     249            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     250            ztm = tmask(ji,jj,jk)                                         ! tmask 
     251            ! 
     252            zn3 = EOS013*zt   & 
     253               &   + EOS103*zs+EOS003 
     254               ! 
     255            zn2 = (EOS022*zt   & 
     256               &   + EOS112*zs+EOS012)*zt   & 
     257               &   + (EOS202*zs+EOS102)*zs+EOS002 
     258               ! 
     259            zn1 = (((EOS041*zt   & 
     260               &   + EOS131*zs+EOS031)*zt   & 
     261               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     262               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     263               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     264               ! 
     265            zn0 = (((((EOS060*zt   & 
     266               &   + EOS150*zs+EOS050)*zt   & 
     267               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     268               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     269               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     270               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     271               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     272               ! 
     273            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     274            ! 
     275            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     276            ! 
     277         END_3D 
     278         ! 
     279      CASE( np_seos )                !==  simplified EOS  ==! 
     280         ! 
     281         DO_3D_11_11( 1, jpkm1 ) 
     282            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     283            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     284            zh  = pdep (ji,jj,jk) 
     285            ztm = tmask(ji,jj,jk) 
     286            ! 
     287            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     288               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     289               &  - rn_nu * zt * zs 
     290               !                                  
     291            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     292         END_3D 
     293         ! 
     294      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
     295         ! 
     296         DO_3D_11_11( 1, jpkm1 ) 
     297            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
     298            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
     299            zh  = pdep (ji,jj,jk) 
     300            ztm = tmask(ji,jj,jk) 
     301            ! 
     302            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     303            !                                  
     304            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     305         END_3D 
     306         ! 
     307      END SELECT 
     308      ! 
     309      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
     310      ! 
     311      IF( ln_timing )   CALL timing_stop('eos-insitu') 
     312      ! 
     313   END SUBROUTINE eos_insitu 
     314 
     315 
     316   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
     317      !!---------------------------------------------------------------------- 
     318      !!                  ***  ROUTINE eos_insitu_pot  *** 
     319      !! 
     320      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
     321      !!      potential volumic mass (Kg/m3) from potential temperature and 
     322      !!      salinity fields using an equation of state selected in the 
     323      !!     namelist. 
     324      !! 
     325      !! ** Action  : - prd  , the in situ density (no units) 
     326      !!              - prhop, the potential volumic mass (Kg/m3) 
     327      !! 
     328      !!---------------------------------------------------------------------- 
     329      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     330      !                                                                ! 2 : salinity               [psu] 
     331      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
     332      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     333      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
     334      ! 
     335      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     336      INTEGER  ::   jdof 
     337      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     338      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     339      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
     340      !!---------------------------------------------------------------------- 
     341      ! 
     342      IF( ln_timing )   CALL timing_start('eos-pot') 
     343      ! 
     344      SELECT CASE ( neos ) 
     345      ! 
     346      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     347         ! 
     348         ! Stochastic equation of state 
     349         IF ( ln_sto_eos ) THEN 
     350            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     351            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     352            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     353            DO jsmp = 1, 2*nn_sto_eos, 2 
     354              zsign(jsmp)   = 1._wp 
     355              zsign(jsmp+1) = -1._wp 
     356            END DO 
     357            ! 
     358            DO_3D_11_11( 1, jpkm1 ) 
     359               ! 
     360               ! compute density (2*nn_sto_eos) times: 
     361               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     362               ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     363               DO jsmp = 1, nn_sto_eos*2 
     364                  jdof   = (jsmp + 1) / 2 
     365                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     366                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     367                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     368                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     369                  ztm    = tmask(ji,jj,jk)                                         ! tmask 
    251370                  ! 
    252371                  zn3 = EOS013*zt   & 
    … …  
    263382                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    264383                     ! 
    265                   zn0 = (((((EOS060*zt   & 
     384                  zn0_sto(jsmp) = (((((EOS060*zt   & 
    266385                     &   + EOS150*zs+EOS050)*zt   & 
    267386                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    … …  
    271390                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    272391                     ! 
    273                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     392                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     393               END DO 
     394               ! 
     395               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     396               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     397               DO jsmp = 1, nn_sto_eos*2 
     398                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    274399                  ! 
    275                   prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
    276                   ! 
     400                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
    277401               END DO 
    278             END DO 
    279          END DO 
    280          ! 
    281       CASE( np_seos )                !==  simplified EOS  ==! 
    282          ! 
    283          DO jk = 1, jpkm1 
    284             DO jj = 1, jpj 
    285                DO ji = 1, jpi 
    286                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    287                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    288                   zh  = pdep (ji,jj,jk) 
    289                   ztm = tmask(ji,jj,jk) 
    290                   ! 
    291                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    292                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    293                      &  - rn_nu * zt * zs 
    294                      !                                  
    295                   prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
    296                END DO 
    297             END DO 
    298          END DO 
    299          ! 
    300       CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    301          ! 
    302          DO jk = 1, jpkm1 
    303             DO jj = 1, jpj 
    304                DO ji = 1, jpi 
    305                   zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
    306                   zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
    307                   zh  = pdep (ji,jj,jk) 
    308                   ztm = tmask(ji,jj,jk) 
    309                   ! 
    310                   zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    311                   !                                  
    312                   prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
    313                END DO 
    314             END DO 
    315          END DO 
    316          ! 
    317       END SELECT 
    318       ! 
    319       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
    320       ! 
    321       IF( ln_timing )   CALL timing_stop('eos-insitu') 
    322       ! 
    323    END SUBROUTINE eos_insitu 
    324  
    325  
    326    SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
    327       !!---------------------------------------------------------------------- 
    328       !!                  ***  ROUTINE eos_insitu_pot  *** 
    329       !! 
    330       !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
    331       !!      potential volumic mass (Kg/m3) from potential temperature and 
    332       !!      salinity fields using an equation of state selected in the 
    333       !!     namelist. 
    334       !! 
    335       !! ** Action  : - prd  , the in situ density (no units) 
    336       !!              - prhop, the potential volumic mass (Kg/m3) 
    337       !! 
    338       !!---------------------------------------------------------------------- 
    339       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
    340       !                                                                ! 2 : salinity               [psu] 
    341       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
    342       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
    343       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    344       ! 
    345       INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
    346       INTEGER  ::   jdof 
    347       REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
    348       REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
    349       REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    350       !!---------------------------------------------------------------------- 
    351       ! 
    352       IF( ln_timing )   CALL timing_start('eos-pot') 
    353       ! 
    354       SELECT CASE ( neos ) 
    355       ! 
    356       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    357          ! 
    358          ! Stochastic equation of state 
    359          IF ( ln_sto_eos ) THEN 
    360             ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
    361             ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
    362             ALLOCATE(zsign(1:2*nn_sto_eos)) 
    363             DO jsmp = 1, 2*nn_sto_eos, 2 
    364               zsign(jsmp)   = 1._wp 
    365               zsign(jsmp+1) = -1._wp 
    366             END DO 
    367             ! 
    368             DO jk = 1, jpkm1 
    369                DO jj = 1, jpj 
    370                   DO ji = 1, jpi 
    371                      ! 
    372                      ! compute density (2*nn_sto_eos) times: 
    373                      ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
    374                      ! (2) for t-dt, s-ds (with the opposite fluctuation) 
    375                      DO jsmp = 1, nn_sto_eos*2 
    376                         jdof   = (jsmp + 1) / 2 
    377                         zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    378                         zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
    379                         zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
    380                         zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
    381                         ztm    = tmask(ji,jj,jk)                                         ! tmask 
    382                         ! 
    383                         zn3 = EOS013*zt   & 
    384                            &   + EOS103*zs+EOS003 
    385                            ! 
    386                         zn2 = (EOS022*zt   & 
    387                            &   + EOS112*zs+EOS012)*zt   & 
    388                            &   + (EOS202*zs+EOS102)*zs+EOS002 
    389                            ! 
    390                         zn1 = (((EOS041*zt   & 
    391                            &   + EOS131*zs+EOS031)*zt   & 
    392                            &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    393                            &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    394                            &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    395                            ! 
    396                         zn0_sto(jsmp) = (((((EOS060*zt   & 
    397                            &   + EOS150*zs+EOS050)*zt   & 
    398                            &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    399                            &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    400                            &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    401                            &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    402                            &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    403                            ! 
    404                         zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
    405                      END DO 
    406                      ! 
    407                      ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
    408                      prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
    409                      DO jsmp = 1, nn_sto_eos*2 
    410                         prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    411                         ! 
    412                         prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
    413                      END DO 
    414                      prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
    415                      prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
    416                   END DO 
    417                END DO 
    418             END DO 
     402               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     403               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     404            END_3D 
    419405            DEALLOCATE(zn0_sto,zn_sto,zsign) 
    420406         ! Non-stochastic equation of state 
    421407         ELSE 
    422             DO jk = 1, jpkm1 
    423                DO jj = 1, jpj 
    424                   DO ji = 1, jpi 
    425                      ! 
    426                      zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    427                      zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    428                      zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    429                      ztm = tmask(ji,jj,jk)                                         ! tmask 
    430                      ! 
    431                      zn3 = EOS013*zt   & 
    432                         &   + EOS103*zs+EOS003 
    433                         ! 
    434                      zn2 = (EOS022*zt   & 
    435                         &   + EOS112*zs+EOS012)*zt   & 
    436                         &   + (EOS202*zs+EOS102)*zs+EOS002 
    437                         ! 
    438                      zn1 = (((EOS041*zt   & 
    439                         &   + EOS131*zs+EOS031)*zt   & 
    440                         &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    441                         &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    442                         &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    443                         ! 
    444                      zn0 = (((((EOS060*zt   & 
    445                         &   + EOS150*zs+EOS050)*zt   & 
    446                         &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    447                         &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    448                         &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    449                         &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    450                         &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    451                         ! 
    452                      zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    453                      ! 
    454                      prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    455                      ! 
    456                      prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
    457                   END DO 
    458                END DO 
    459             END DO 
    460          ENDIF 
    461           
    462       CASE( np_seos )                !==  simplified EOS  ==! 
    463          ! 
    464          DO jk = 1, jpkm1 
    465             DO jj = 1, jpj 
    466                DO ji = 1, jpi 
    467                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    468                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    469                   zh  = pdep (ji,jj,jk) 
    470                   ztm = tmask(ji,jj,jk) 
    471                   !                                                     ! potential density referenced at the surface 
    472                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
    473                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
    474                      &  - rn_nu * zt * zs 
    475                   prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
    476                   !                                                     ! density anomaly (masked) 
    477                   zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
    478                   prd(ji,jj,jk) = zn * r1_rho0 * ztm 
    479                   ! 
    480                END DO 
    481             END DO 
    482          END DO 
    483          ! 
    484       CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    485          ! 
    486          DO jk = 1, jpkm1 
    487             DO jj = 1, jpj 
    488                DO ji = 1, jpi 
    489                   zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
    490                   zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
    491                   zh  = pdep (ji,jj,jk) 
    492                   ztm = tmask(ji,jj,jk) 
    493                   !                                                     ! potential density referenced at the surface 
    494                   zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    495                   prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
    496                   !                                                     ! density anomaly (masked) 
    497                   prd(ji,jj,jk) = zn * r1_rho0 * ztm 
    498                   ! 
    499                END DO 
    500             END DO 
    501          END DO 
    502          ! 
    503       END SELECT 
    504       ! 
    505       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
    506       ! 
    507       IF( ln_timing )   CALL timing_stop('eos-pot') 
    508       ! 
    509    END SUBROUTINE eos_insitu_pot 
    510  
    511  
    512    SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
    513       !!---------------------------------------------------------------------- 
    514       !!                  ***  ROUTINE eos_insitu_2d  *** 
    515       !! 
    516       !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
    517       !!      potential temperature and salinity using an equation of state 
    518       !!      selected in the nameos namelist. * 2D field case 
    519       !! 
    520       !! ** Action  : - prd , the in situ density (no units) (unmasked) 
    521       !! 
    522       !!---------------------------------------------------------------------- 
    523       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
    524       !                                                           ! 2 : salinity               [psu] 
    525       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
    526       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    527       ! 
    528       INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    529       REAL(wp) ::   zt , zh , zs              ! local scalars 
    530       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
    531       !!---------------------------------------------------------------------- 
    532       ! 
    533       IF( ln_timing )   CALL timing_start('eos2d') 
    534       ! 
    535       prd(:,:) = 0._wp 
    536       ! 
    537       SELECT CASE( neos ) 
    538       ! 
    539       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    540          ! 
    541          DO jj = 1, jpjm1 
    542             DO ji = 1, fs_jpim1   ! vector opt. 
    543                ! 
    544                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    545                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    546                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     408            DO_3D_11_11( 1, jpkm1 ) 
     409               ! 
     410               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     411               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     412               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     413               ztm = tmask(ji,jj,jk)                                         ! tmask 
    547414               ! 
    548415               zn3 = EOS013*zt   & 
    … …  
    569436               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    570437               ! 
    571                prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly 
    572                ! 
    573             END DO 
    574          END DO 
    575          ! 
    576          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
    577          ! 
     438               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     439               ! 
     440               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     441            END_3D 
     442         ENDIF 
     443          
    578444      CASE( np_seos )                !==  simplified EOS  ==! 
    579445         ! 
    580          DO jj = 1, jpjm1 
    581             DO ji = 1, fs_jpim1   ! vector opt. 
    582                ! 
    583                zt    = pts  (ji,jj,jp_tem)  - 10._wp 
    584                zs    = pts  (ji,jj,jp_sal)  - 35._wp 
    585                zh    = pdep (ji,jj)                         ! depth at the partial step level 
    586                ! 
    587                zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    588                   &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    589                   &  - rn_nu * zt * zs 
    590                   ! 
    591                prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
    592                ! 
    593             END DO 
    594          END DO 
    595          ! 
    596          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
     446         DO_3D_11_11( 1, jpkm1 ) 
     447            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     448            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     449            zh  = pdep (ji,jj,jk) 
     450            ztm = tmask(ji,jj,jk) 
     451            !                                                     ! potential density referenced at the surface 
     452            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
     453               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
     454               &  - rn_nu * zt * zs 
     455            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
     456            !                                                     ! density anomaly (masked) 
     457            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
     458            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
     459            ! 
     460         END_3D 
     461         ! 
     462      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
     463         ! 
     464         DO_3D_11_11( 1, jpkm1 ) 
     465            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
     466            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
     467            zh  = pdep (ji,jj,jk) 
     468            ztm = tmask(ji,jj,jk) 
     469            !                                                     ! potential density referenced at the surface 
     470            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     471            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
     472            !                                                     ! density anomaly (masked) 
     473            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
     474            ! 
     475         END_3D 
     476         ! 
     477      END SELECT 
     478      ! 
     479      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
     480      ! 
     481      IF( ln_timing )   CALL timing_stop('eos-pot') 
     482      ! 
     483   END SUBROUTINE eos_insitu_pot 
     484 
     485 
     486   SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
     487      !!---------------------------------------------------------------------- 
     488      !!                  ***  ROUTINE eos_insitu_2d  *** 
     489      !! 
     490      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
     491      !!      potential temperature and salinity using an equation of state 
     492      !!      selected in the nameos namelist. * 2D field case 
     493      !! 
     494      !! ** Action  : - prd , the in situ density (no units) (unmasked) 
     495      !! 
     496      !!---------------------------------------------------------------------- 
     497      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     498      !                                                           ! 2 : salinity               [psu] 
     499      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
     500      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
     501      ! 
     502      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     503      REAL(wp) ::   zt , zh , zs              ! local scalars 
     504      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     505      !!---------------------------------------------------------------------- 
     506      ! 
     507      IF( ln_timing )   CALL timing_start('eos2d') 
     508      ! 
     509      prd(:,:) = 0._wp 
     510      ! 
     511      SELECT CASE( neos ) 
     512      ! 
     513      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     514         ! 
     515         DO_2D_11_11 
     516            ! 
     517            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     518            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     519            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     520            ! 
     521            zn3 = EOS013*zt   & 
     522               &   + EOS103*zs+EOS003 
     523               ! 
     524            zn2 = (EOS022*zt   & 
     525               &   + EOS112*zs+EOS012)*zt   & 
     526               &   + (EOS202*zs+EOS102)*zs+EOS002 
     527               ! 
     528            zn1 = (((EOS041*zt   & 
     529               &   + EOS131*zs+EOS031)*zt   & 
     530               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     531               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     532               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     533               ! 
     534            zn0 = (((((EOS060*zt   & 
     535               &   + EOS150*zs+EOS050)*zt   & 
     536               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     537               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     538               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     539               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     540               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     541               ! 
     542            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     543            ! 
     544            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly 
     545            ! 
     546         END_2D 
     547         ! 
     548      CASE( np_seos )                !==  simplified EOS  ==! 
     549         ! 
     550         DO_2D_11_11 
     551            ! 
     552            zt    = pts  (ji,jj,jp_tem)  - 10._wp 
     553            zs    = pts  (ji,jj,jp_sal)  - 35._wp 
     554            zh    = pdep (ji,jj)                         ! depth at the partial step level 
     555            ! 
     556            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     557               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     558               &  - rn_nu * zt * zs 
     559               ! 
     560            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
     561            ! 
     562         END_2D 
    597563         ! 
    598564      CASE( np_leos )                !==  ISOMIP EOS  ==! 
    599565         ! 
    600          DO jj = 1, jpjm1 
    601             DO ji = 1, fs_jpim1   ! vector opt. 
    602                ! 
    603                zt    = pts  (ji,jj,jp_tem)  - (-1._wp) 
    604                zs    = pts  (ji,jj,jp_sal)  - 34.2_wp 
    605                zh    = pdep (ji,jj)                         ! depth at the partial step level 
    606                ! 
    607                zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    608                   ! 
    609                prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
    610                ! 
    611             END DO 
    612          END DO 
    613          ! 
    614          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
     566         DO_2D_11_11 
     567            ! 
     568            zt    = pts  (ji,jj,jp_tem)  - (-1._wp) 
     569            zs    = pts  (ji,jj,jp_sal)  - 34.2_wp 
     570            zh    = pdep (ji,jj)                         ! depth at the partial step level 
     571            ! 
     572            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     573            ! 
     574            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
     575            ! 
     576         END_2D 
     577         ! 
    615578         ! 
    616579      END SELECT 
    617580      ! 
    618       IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
     581      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
    619582      ! 
    620583      IF( ln_timing )   CALL timing_stop('eos2d') 
    … …  
    648611      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    649612         ! 
    650          DO jk = 1, jpkm1 
    651             DO jj = 1, jpj 
    652                DO ji = 1, jpi 
    653                   ! 
    654                   zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    655                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    656                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    657                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    658                   ! 
    659                   ! alpha 
    660                   zn3 = ALP003 
    661                   ! 
    662                   zn2 = ALP012*zt + ALP102*zs+ALP002 
    663                   ! 
    664                   zn1 = ((ALP031*zt   & 
    665                      &   + ALP121*zs+ALP021)*zt   & 
    666                      &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    667                      &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    668                      ! 
    669                   zn0 = ((((ALP050*zt   & 
    670                      &   + ALP140*zs+ALP040)*zt   & 
    671                      &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    672                      &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    673                      &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    674                      &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    675                      ! 
    676                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    677                   ! 
    678                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
    679                   ! 
    680                   ! beta 
    681                   zn3 = BET003 
    682                   ! 
    683                   zn2 = BET012*zt + BET102*zs+BET002 
    684                   ! 
    685                   zn1 = ((BET031*zt   & 
    686                      &   + BET121*zs+BET021)*zt   & 
    687                      &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    688                      &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    689                      ! 
    690                   zn0 = ((((BET050*zt   & 
    691                      &   + BET140*zs+BET040)*zt   & 
    692                      &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    693                      &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    694                      &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    695                      &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    696                      ! 
    697                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    698                   ! 
    699                   pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
    700                   ! 
    701                END DO 
    702             END DO 
    703          END DO 
     613         DO_3D_11_11( 1, jpkm1 ) 
     614            ! 
     615            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     616            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     617            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     618            ztm = tmask(ji,jj,jk)                                         ! tmask 
     619            ! 
     620            ! alpha 
     621            zn3 = ALP003 
     622            ! 
     623            zn2 = ALP012*zt + ALP102*zs+ALP002 
     624            ! 
     625            zn1 = ((ALP031*zt   & 
     626               &   + ALP121*zs+ALP021)*zt   & 
     627               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     628               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     629               ! 
     630            zn0 = ((((ALP050*zt   & 
     631               &   + ALP140*zs+ALP040)*zt   & 
     632               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     633               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     634               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     635               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     636               ! 
     637            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     638            ! 
     639            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
     640            ! 
     641            ! beta 
     642            zn3 = BET003 
     643            ! 
     644            zn2 = BET012*zt + BET102*zs+BET002 
     645            ! 
     646            zn1 = ((BET031*zt   & 
     647               &   + BET121*zs+BET021)*zt   & 
     648               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     649               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     650               ! 
     651            zn0 = ((((BET050*zt   & 
     652               &   + BET140*zs+BET040)*zt   & 
     653               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     654               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     655               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     656               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     657               ! 
     658            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     659            ! 
     660            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
     661            ! 
     662         END_3D 
    704663         ! 
    705664      CASE( np_seos )                  !==  simplified EOS  ==! 
    706665         ! 
    707          DO jk = 1, jpkm1 
    708             DO jj = 1, jpj 
    709                DO ji = 1, jpi 
    710                   zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    711                   zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    712                   zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
    713                   ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
    714                   ! 
    715                   zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    716                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
    717                   ! 
    718                   zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    719                   pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
    720                   ! 
    721                END DO 
    722             END DO 
    723          END DO 
     666         DO_3D_11_11( 1, jpkm1 ) 
     667            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     668            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     669            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
     670            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
     671            ! 
     672            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     673            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
     674            ! 
     675            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     676            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
     677            ! 
     678         END_3D 
    724679         ! 
    725680      CASE( np_leos )                  !==  linear ISOMIP EOS  ==! 
    726681         ! 
    727          DO jk = 1, jpkm1 
    728             DO jj = 1, jpj 
    729                DO ji = 1, jpi 
    730                   zt  = pts (ji,jj,jk,jp_tem) - (-1._wp) 
    731                   zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    732                   zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point 
    733                   ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask 
    734                   ! 
    735                   zn  = rn_a0 * rho0 
    736                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
    737                   ! 
    738                   zn  = rn_b0 * rho0 
    739                   pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
    740                   ! 
    741                END DO 
    742             END DO 
    743          END DO 
     682         DO_3D_11_11( 1, jpkm1 ) 
     683            zt  = pts (ji,jj,jk,jp_tem) - (-1._wp) 
     684            zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     685            zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point 
     686            ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask 
     687            ! 
     688            zn  = rn_a0 * rho0 
     689            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
     690            ! 
     691            zn  = rn_b0 * rho0 
     692            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
     693            ! 
     694         END_3D 
    744695         ! 
    745696      CASE DEFAULT 
    … …  
    749700      END SELECT 
    750701      ! 
    751       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
    752          &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
     702      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
     703         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
    753704      ! 
    754705      IF( ln_timing )   CALL timing_stop('rab_3d') 
    … …  
    783734      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    784735         ! 
    785          DO jj = 1, jpjm1 
    786             DO ji = 1, fs_jpim1   ! vector opt. 
    787                ! 
    788                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    789                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    790                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    791                ! 
    792                ! alpha 
    793                zn3 = ALP003 
    794                ! 
    795                zn2 = ALP012*zt + ALP102*zs+ALP002 
    796                ! 
    797                zn1 = ((ALP031*zt   & 
    798                   &   + ALP121*zs+ALP021)*zt   & 
    799                   &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    800                   &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    801                   ! 
    802                zn0 = ((((ALP050*zt   & 
    803                   &   + ALP140*zs+ALP040)*zt   & 
    804                   &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    805                   &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    806                   &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    807                   &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    808                   ! 
    809                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    810                ! 
    811                pab(ji,jj,jp_tem) = zn * r1_rho0 
    812                ! 
    813                ! beta 
    814                zn3 = BET003 
    815                ! 
    816                zn2 = BET012*zt + BET102*zs+BET002 
    817                ! 
    818                zn1 = ((BET031*zt   & 
    819                   &   + BET121*zs+BET021)*zt   & 
    820                   &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    821                   &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    822                   ! 
    823                zn0 = ((((BET050*zt   & 
    824                   &   + BET140*zs+BET040)*zt   & 
    825                   &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    826                   &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    827                   &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    828                   &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    829                   ! 
    830                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    831                ! 
    832                pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
    833                ! 
    834                ! 
    835             END DO 
    836          END DO 
    837          !                            ! Lateral boundary conditions 
    838          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     736         DO_2D_11_11 
     737            ! 
     738            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     739            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     740            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     741            ! 
     742            ! alpha 
     743            zn3 = ALP003 
     744            ! 
     745            zn2 = ALP012*zt + ALP102*zs+ALP002 
     746            ! 
     747            zn1 = ((ALP031*zt   & 
     748               &   + ALP121*zs+ALP021)*zt   & 
     749               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     750               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     751               ! 
     752            zn0 = ((((ALP050*zt   & 
     753               &   + ALP140*zs+ALP040)*zt   & 
     754               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     755               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     756               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     757               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     758               ! 
     759            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     760            ! 
     761            pab(ji,jj,jp_tem) = zn * r1_rho0 
     762            ! 
     763            ! beta 
     764            zn3 = BET003 
     765            ! 
     766            zn2 = BET012*zt + BET102*zs+BET002 
     767            ! 
     768            zn1 = ((BET031*zt   & 
     769               &   + BET121*zs+BET021)*zt   & 
     770               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     771               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     772               ! 
     773            zn0 = ((((BET050*zt   & 
     774               &   + BET140*zs+BET040)*zt   & 
     775               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     776               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     777               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     778               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     779               ! 
     780            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     781            ! 
     782            pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
     783            ! 
     784            ! 
     785         END_2D 
    839786         ! 
    840787      CASE( np_seos )                  !==  simplified EOS  ==! 
    841788         ! 
    842          DO jj = 1, jpjm1 
    843             DO ji = 1, fs_jpim1   ! vector opt. 
    844                ! 
    845                zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    846                zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    847                zh    = pdep (ji,jj)                   ! depth at the partial step level 
    848                ! 
    849                zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    850                pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
    851                ! 
    852                zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    853                pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
    854                ! 
    855             END DO 
    856          END DO 
    857          !                            ! Lateral boundary conditions 
    858          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     789         DO_2D_11_11 
     790            ! 
     791            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     792            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     793            zh    = pdep (ji,jj)                   ! depth at the partial step level 
     794            ! 
     795            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     796            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
     797            ! 
     798            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     799            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
     800            ! 
     801         END_2D 
    859802         ! 
    860803      CASE( np_leos )                  !==  linear ISOMIP EOS  ==! 
    861804         ! 
    862          DO jj = 1, jpjm1 
    863             DO ji = 1, fs_jpim1   ! vector opt. 
    864                ! 
    865                zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0) 
    866                zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    867                zh    = pdep (ji,jj)                   ! depth at the partial step level 
    868                ! 
    869                zn  = rn_a0 * rho0 
    870                pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
    871                ! 
    872                zn  = rn_b0 * rho0 
    873                pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
    874                ! 
    875             END DO 
    876          END DO 
    877          ! 
    878          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    ! Lateral boundary conditions 
     805         DO_2D_11_11 
     806            ! 
     807            zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0) 
     808            zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     809            zh    = pdep (ji,jj)                   ! depth at the partial step level 
     810            ! 
     811            zn  = rn_a0 * rho0 
     812            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
     813            ! 
     814            zn  = rn_b0 * rho0 
     815            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
     816            ! 
     817         END_2D 
    879818         ! 
    880819      CASE DEFAULT 
    … …  
    884823      END SELECT 
    885824      ! 
    886       IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
    887          &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
     825      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
     826         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
    888827      ! 
    889828      IF( ln_timing )   CALL timing_stop('rab_2d') 
    … …  
    1026965      IF( ln_timing )   CALL timing_start('bn2') 
    1027966      ! 
    1028       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 ) 
    1029          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90 
    1030             DO ji = 1, jpi 
    1031                zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
    1032                   &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
    1033                   ! 
    1034                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
    1035                zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
    1036                ! 
    1037                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    1038                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    1039                   &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
    1040             END DO 
    1041          END DO 
    1042       END DO 
    1043       ! 
    1044       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
     967      DO_3D_11_11( 2, jpkm1 ) 
     968         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
     969            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
     970            ! 
     971         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     972         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
     973         ! 
     974         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
     975            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
     976            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     977      END_3D 
     978      ! 
     979      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
    1045980      ! 
    1046981      IF( ln_timing )   CALL timing_stop('bn2') 
    … …  
    10781013      z1_T0   = 1._wp/40._wp 
    10791014      ! 
    1080       DO jj = 1, jpj 
    1081          DO ji = 1, jpi 
    1082             ! 
    1083             zt  = ctmp   (ji,jj) * z1_T0 
    1084             zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
    1085             ztm = tmask(ji,jj,1) 
    1086             ! 
    1087             zn = ((((-2.1385727895e-01_wp*zt   & 
    1088                &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
    1089                &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
    1090                &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
    1091                &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
    1092                &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
    1093                &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
    1094                &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
    1095                ! 
    1096             zd = (2.0035003456_wp*zt   & 
    1097                &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
    1098                &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
    1099                ! 
    1100             ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
    1101                ! 
    1102          END DO 
    1103       END DO 
     1015      DO_2D_11_11 
     1016         ! 
     1017         zt  = ctmp   (ji,jj) * z1_T0 
     1018         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
     1019         ztm = tmask(ji,jj,1) 
     1020         ! 
     1021         zn = ((((-2.1385727895e-01_wp*zt   & 
     1022            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
     1023            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
     1024            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
     1025            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
     1026            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
     1027            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
     1028            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
     1029            ! 
     1030         zd = (2.0035003456_wp*zt   & 
     1031            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
     1032            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
     1033            ! 
     1034         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
     1035            ! 
     1036      END_2D 
    11041037      ! 
    11051038      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct') 
    … …  
    11331066         ! 
    11341067         z1_S0 = 1._wp / 35.16504_wp 
    1135          DO jj = 1, jpj 
    1136             DO ji = 1, jpi 
    1137                zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
    1138                ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
    1139                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
    1140             END DO 
    1141          END DO 
     1068         DO_2D_11_11 
     1069            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
     1070            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
     1071               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
     1072         END_2D 
    11421073         ptf(:,:) = ptf(:,:) * psal(:,:) 
    11431074         ! 
    11441075         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
    11451076         ! 
    1146       CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==! 
     1077      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    11471078         ! 
    11481079         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
    … …  
    11901121         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep 
    11911122         ! 
    1192       CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==! 
     1123      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    11931124         ! 
    11941125         ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   & 
    … …  
    12421173      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    12431174         ! 
    1244          DO jk = 1, jpkm1 
    1245             DO jj = 1, jpj 
    1246                DO ji = 1, jpi 
    1247                   ! 
    1248                   zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    1249                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    1250                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    1251                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    1252                   ! 
    1253                   ! potential energy non-linear anomaly 
    1254                   zn2 = (PEN012)*zt   & 
    1255                      &   + PEN102*zs+PEN002 
    1256                      ! 
    1257                   zn1 = ((PEN021)*zt   & 
    1258                      &   + PEN111*zs+PEN011)*zt   & 
    1259                      &   + (PEN201*zs+PEN101)*zs+PEN001 
    1260                      ! 
    1261                   zn0 = ((((PEN040)*zt   & 
    1262                      &   + PEN130*zs+PEN030)*zt   & 
    1263                      &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
    1264                      &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
    1265                      &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
    1266                      ! 
    1267                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1268                   ! 
    1269                   ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm 
    1270                   ! 
    1271                   ! alphaPE non-linear anomaly 
    1272                   zn2 = APE002 
    1273                   ! 
    1274                   zn1 = (APE011)*zt   & 
    1275                      &   + APE101*zs+APE001 
    1276                      ! 
    1277                   zn0 = (((APE030)*zt   & 
    1278                      &   + APE120*zs+APE020)*zt   & 
    1279                      &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
    1280                      &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
    1281                      ! 
    1282                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1283                   !                               
    1284                   pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 
    1285                   ! 
    1286                   ! betaPE non-linear anomaly 
    1287                   zn2 = BPE002 
    1288                   ! 
    1289                   zn1 = (BPE011)*zt   & 
    1290                      &   + BPE101*zs+BPE001 
    1291                      ! 
    1292                   zn0 = (((BPE030)*zt   & 
    1293                      &   + BPE120*zs+BPE020)*zt   & 
    1294                      &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
    1295                      &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
    1296                      ! 
    1297                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1298                   !                               
    1299                   pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 
    1300                   ! 
    1301                END DO 
    1302             END DO 
    1303          END DO 
     1175         DO_3D_11_11( 1, jpkm1 ) 
     1176            ! 
     1177            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     1178            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     1179            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1180            ztm = tmask(ji,jj,jk)                                         ! tmask 
     1181            ! 
     1182            ! potential energy non-linear anomaly 
     1183            zn2 = (PEN012)*zt   & 
     1184               &   + PEN102*zs+PEN002 
     1185               ! 
     1186            zn1 = ((PEN021)*zt   & 
     1187               &   + PEN111*zs+PEN011)*zt   & 
     1188               &   + (PEN201*zs+PEN101)*zs+PEN001 
     1189               ! 
     1190            zn0 = ((((PEN040)*zt   & 
     1191               &   + PEN130*zs+PEN030)*zt   & 
     1192               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
     1193               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
     1194               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
     1195               ! 
     1196            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1197            ! 
     1198            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm 
     1199            ! 
     1200            ! alphaPE non-linear anomaly 
     1201            zn2 = APE002 
     1202            ! 
     1203            zn1 = (APE011)*zt   & 
     1204               &   + APE101*zs+APE001 
     1205               ! 
     1206            zn0 = (((APE030)*zt   & 
     1207               &   + APE120*zs+APE020)*zt   & 
     1208               &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
     1209               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
     1210               ! 
     1211            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1212            !                               
     1213            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 
     1214            ! 
     1215            ! betaPE non-linear anomaly 
     1216            zn2 = BPE002 
     1217            ! 
     1218            zn1 = (BPE011)*zt   & 
     1219               &   + BPE101*zs+BPE001 
     1220               ! 
     1221            zn0 = (((BPE030)*zt   & 
     1222               &   + BPE120*zs+BPE020)*zt   & 
     1223               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
     1224               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
     1225               ! 
     1226            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1227            !                               
     1228            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 
     1229            ! 
     1230         END_3D 
    13041231         ! 
    13051232      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==! 
    13061233         ! 
    1307          DO jk = 1, jpkm1 
    1308             DO jj = 1, jpj 
    1309                DO ji = 1, jpi 
    1310                   zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
    1311                   zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
    1312                   zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
    1313                   ztm = tmask(ji,jj,jk)                ! tmask 
    1314                   zn  = 0.5_wp * zh * r1_rho0 * ztm 
    1315                   !                                    ! Potential Energy 
    1316                   ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
    1317                   !                                    ! alphaPE 
    1318                   pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
    1319                   pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
    1320                   ! 
    1321                END DO 
    1322             END DO 
    1323          END DO 
     1234         DO_3D_11_11( 1, jpkm1 ) 
     1235            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
     1236            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
     1237            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
     1238            ztm = tmask(ji,jj,jk)                ! tmask 
     1239            zn  = 0.5_wp * zh * r1_rho0 * ztm 
     1240            !                                    ! Potential Energy 
     1241            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
     1242            !                                    ! alphaPE 
     1243            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
     1244            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
     1245            ! 
     1246         END_3D 
    13241247         ! 
    13251248      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    13261249         ! 
    1327          DO jk = 1, jpkm1 
    1328             DO jj = 1, jpj 
    1329                DO ji = 1, jpi 
    1330                   zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0) 
    1331                   zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    1332                   zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point 
    1333                   ztm = tmask(ji,jj,jk)                  ! tmask 
    1334                   zn  = 0.5_wp * zh * r1_rho0 * ztm 
    1335                   !                                    ! Potential Energy 
    1336                   ppen(ji,jj,jk) = 0. 
    1337                   !                                    ! alphaPE 
    1338                   pab_pe(ji,jj,jk,jp_tem) = 0. 
    1339                   pab_pe(ji,jj,jk,jp_sal) = 0. 
    1340                   ! 
    1341                END DO 
    1342             END DO 
    1343          END DO 
     1250         DO_3D_11_11( 1, jpkm1 ) 
     1251            zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0) 
     1252            zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     1253            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point 
     1254            ztm = tmask(ji,jj,jk)                  ! tmask 
     1255            zn  = 0.5_wp * zh * r1_rho0 * ztm 
     1256            !                                    ! Potential Energy 
     1257            ppen(ji,jj,jk) = 0. 
     1258            !                                    ! alphaPE 
     1259            pab_pe(ji,jj,jk,jp_tem) = 0. 
     1260            pab_pe(ji,jj,jk,jp_sal) = 0. 
     1261            ! 
     1262         END_3D 
    13441263         ! 
    13451264      CASE DEFAULT 
    … …  
    13651284      INTEGER  ::   ioptio   ! local integer 
    13661285      !! 
    1367       NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS   , ln_LEOS, & 
    1368          &             rn_a0    , rn_b0   , rn_lambda1, rn_mu1 , & 
    1369          &                                  rn_lambda2, rn_mu2 , rn_nu 
    1370       !!---------------------------------------------------------------------- 
    1371       ! 
    1372       REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state 
     1286      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, ln_LEOS, rn_a0, rn_b0, & 
     1287         &             rn_lambda1, rn_mu1, rn_lambda2, rn_mu2, rn_nu 
     1288      !!---------------------------------------------------------------------- 
     1289      ! 
    13731290      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 
    13741291901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' ) 
    13751292      ! 
    1376       REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state 
    13771293      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 
    13781294902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' ) 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/isfcavgam.F90

    r12077 r12930  
    9191         pgs(:,:) = rn_gammas0 
    9292      CASE ( 'vel' ) ! gamma is proportional to u* 
    93          CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs ) 
     93         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs ) 
    9494      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 
    95          CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs ) 
     95         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pgt, pgs ) 
    9696      CASE DEFAULT 
    9797         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/isfstp.F90

    r12077 r12930  
    250250      IF ( l_isfoasis .AND. ln_isf ) THEN 
    251251         ! 
    252          CALL ctl_stop( ' ln_ctl and ice shelf not tested' ) 
     252         CALL ctl_stop( 'namelist combination ln_cpl and ln_isf not tested' ) 
    253253         ! 
    254254         ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation  
    … …  
    291291      !!---------------------------------------------------------------------- 
    292292      ! 
    293       REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    294293      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 
    295294901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' ) 
    296295      ! 
    297       REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    298296      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 
    299297902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' ) 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/istate.F90

    r12353 r12930  
    4141   PUBLIC   istate_init   ! routine called by step.F90 
    4242 
     43   !! * Substitutions 
     44#  include "do_loop_substitute.h90" 
    4345   !!---------------------------------------------------------------------- 
    4446   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7577      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    7678      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    77       ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk 
     79      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    7880      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    7981#if defined key_agrif 
    … …  
    9092         !                                    ! --------------- 
    9193         numror = 0                           ! define numror = 0 -> no restart file to read 
    92          neuler = 0                           ! Set time-step indicator at nit000 (euler forward) 
     94         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward) 
    9395         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    9496         !                                    ! Initialization of ocean to zero 
    … …  
    103105               ! Apply minimum wetdepth criterion 
    104106               ! 
    105                DO jj = 1,jpj 
    106                   DO ji = 1,jpi 
    107                      IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
    108                         ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    109                      ENDIF 
    110                   END DO 
    111                END DO  
     107               DO_2D_11_11 
     108                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
     109                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
     110                  ENDIF 
     111               END_2D 
    112112            ENDIF  
    113113            uu  (:,:,:,Kbb) = 0._wp 
    … …  
    159159      ! 
    160160!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 
    161       DO jk = 1, jpkm1 
    162          DO jj = 1, jpj 
    163             DO ji = 1, jpi 
    164                uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
    165                vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
    166                ! 
    167                uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
    168                vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
    169             END DO 
    170          END DO 
    171       END DO 
     161      DO_3D_11_11( 1, jpkm1 ) 
     162         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     163         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
     164         ! 
     165         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
     166         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
     167      END_3D 
    172168      ! 
    173169      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/sbcfwb.F90

    r12489 r12930  
    151151         ENDIF    
    152152         !                                         ! Update fwfold if new year start 
    153          ikty = 365 * 86400 / rn_Dt               !!bug  use of 365 days leap year or 360d year !!!!!!! 
     153         ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!! 
    154154         IF( MOD( kt, ikty ) == 0 ) THEN 
    155155            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow 

  • NEMO/branches/2020/r12581_ticket2418/tests/ISOMIP+/MY_SRC/tradmp.F90

    r12353 r12930  
    5151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1) 
    5252 
     53   !! * Substitutions 
     54#  include "do_loop_substitute.h90" 
    5355   !!---------------------------------------------------------------------- 
    5456   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    110112      CASE( 0 )                        !*  newtonian damping throughout the water column  *! 
    111113         DO jn = 1, jpts 
    112             DO jk = 1, jpkm1 
    113                DO jj = 2, jpjm1 
    114                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    115                      pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           & 
    116                         &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) ) 
    117                   END DO 
    118                END DO 
    119             END DO 
     114            DO_3D_00_00( 1, jpkm1 ) 
     115               pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           & 
     116                  &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) ) 
     117            END_3D 
    120118         END DO 
    121119         ! 
    122120      CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *! 
    123          DO jk = 1, jpkm1 
    124             DO jj = 2, jpjm1 
    125                DO ji = fs_2, fs_jpim1   ! vector opt. 
    126                   IF( avt(ji,jj,jk) <= avt_c ) THEN 
    127                      pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    128                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
    129                      pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
    130                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    131                   ENDIF 
    132                END DO 
    133             END DO 
    134          END DO 
     121         DO_3D_00_00( 1, jpkm1 ) 
     122            IF( avt(ji,jj,jk) <= avt_c ) THEN 
     123               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     124                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     125               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     126                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
     127            ENDIF 
     128         END_3D 
    135129         ! 
    136130      CASE ( 2 )                       !*  no damping in the mixed layer   *! 
    137          DO jk = 1, jpkm1 
    138             DO jj = 2, jpjm1 
    139                DO ji = fs_2, fs_jpim1   ! vector opt. 
    140                   IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
    141                      pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    142                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
    143                      pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
    144                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    145                   ENDIF 
    146                END DO 
    147             END DO 
    148          END DO 
     131         DO_3D_00_00( 1, jpkm1 ) 
     132            IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
     133               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     134                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     135               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     136                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
     137            ENDIF 
     138         END_3D 
    149139         ! 
    150140      END SELECT 
    … …  
    157147      ENDIF 
    158148      !                           ! Control print 
    159       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
    160          &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     149      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
     150         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    161151      ! 
    162152      IF( ln_timing )   CALL timing_stop('tra_dmp') 
    … …  
    178168      !!---------------------------------------------------------------------- 
    179169      ! 
    180       REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation 
    181170      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    182171901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' ) 
    183172      ! 
    184       REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation 
    185173      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    186174902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.